Kalman Filter and Variants for Estimation in 2DOF Serial Flexible Link and Joint Using Fractional Order PID Controller

Abstract

Robotic manipulators have been widely used in industries, mainly to move tools into different specific positions. Thus, it has become necessary to have accurate knowledge about the tool position using forward kinematics after accessing the angular locations of limbs. This paper presents a simulation study in which an encoder attached to the limbs gathers information about the angular positions. The measured angles are applied to the Kalman Filter (KF) and its variants for state estimation. This work focuses on the use of fractional order controllers with a Two Degree of Freedom Serial Flexible Links (2DSFL) and Two Degree of Freedom Serial Flexible Joint (2DSFJ) and undertakes simulations with noise and a square wave as input. The fractional order controllers fit better with the system properties than integer order controllers. The KF and its variants use an unknown and assumed process and measurement noise matrices to predict the actual data. An optimisation problem is proposed to achieve reasonable estimations with the updated covariance matrices.

1. Introduction

Until the mid 1900s, time-invariant as well as SISO (single-input–single-output) systems focused on classical approaches such as frequency domain analysis such that a transfer function model could represent the system. However, the necessities imposed by control engineering in the 1950s and 1960s significantly drove the research into the state space approach towards MIMO (multi-input–multi-output) techniques describing the dynamic evolution of the system using its states [1]. One can gain information about system properties and achieve control objectives using a suitably designed controller [2]. States are not always accessible, as these are internal system parameters. State estimation algorithms can overcome this drawback as it involves the fusing of the mathematical model and input–output data measurements. These algorithms have been a fundamental part of many control problems. For a naive engineer, it is easier to place the sensor on everything and everywhere, but this may not be practical and makes the system expensive and difficult to maintain. Many times, it is not practical to determine the internal state of a system from only a few measurements. State estimation methods are helpful in these scenarios. The presence of plant uncertainties and noises is another primary reason for state estimation [3]. It is not easy to model a physical system accurately as a physical system operates in the real environment. The initial condition of such a model is usually unknown. Even if a system model is initially accurate, the model parameters will vary due to the timely wear and tear of the physical model. Thus, it is quite impractical to determine system states, making state estimation methods necessary [4].

Figure 1 presents the general analogy for the state estimation techniques used in MIMO systems. Usually, linear state estimations are performed through the implementation of Kalman Filters, but this may not be a good option in the case of dynamic models, measurements that are highly nonlinear and non-Gaussian noise. Filters such as the particle filter (PF) [5], Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) [6,7] give a generalised solution for such systems, where low performance is caused due to intractable linearisation and Gaussian approximations. Traditionally, these Bayesian methods are significantly used in moving robot scenarios [8]. As an example, Jassemi et al. presented the application of the EKF to a rigid 2DOF robotic arm using measurements of tool acceleration and arm angle to improve trajectory tracking [9]. A similar extension of this work has been presented by Quigley et al. [10], where the EKF was used for a 7DOF robotic manipulator. In [11], the authors used the Kalman filter for single flexible links, where they used joint angle and tool acceleration measurements for estimation. Similarly, ref. [12] used an EKF for a manipulator with two links, with tool positions and joint action as the measurements. There are few robotic applications in which the applications have been embedded with sensor fusion techniques using particle filters [13,14] and the simulated data. The particle filter is additionally promoted as one can fully utilize the probability density function (PDF), which makes it easier to design control laws and therefore obtain more significant performance improvements. This also allows hard limits on system parameters and gives a benchmark for an easier solution, as in the case of the EKF.

Modern manufacturing requires robots that are capable of functioning at accuracies beyond those of present-day industrial robots. Accurately moving a tool on a predefined path is the aim. However, many modern industrial robot control approaches involve joint angle measurements of the manipulator [15]. Gunnarsson et al. [16] described a method to achieve improved accuracy through the identification of uncertain parameters in robot system by the use of the iterative learning control (ILC) method. Precise modelling has been a significant part of gaining information about systems [17]. The more accurate the modelling of the physical plant, the better the possible control of the system [18]. As the mathematical models have always been important, the use of fractional calculus gives users an upper hand when determining the model. Fractional calculus has been used in nearly all domains of science and has been shown to be advantageous in the recent literature [19]. It has been used in areas such as diode modelling [20] and in the modelling of chemical systems, such as to determine ethanol concentration [21]. Different fractional calculus methods and applications have been discussed [22], and the application of fractional calculus to robotic manipulators is one key area to prove its viability. The fractional order models, when compared to integer order models of any real-time system, are much more accurate [23]. As compared to integer order controllers, these provide better responses, as shown in [24]. Wang et al. [25] demonstrate a novel fractional sliding mode control structure for integer and fractional order systems. These studies suggest that such a controller also shows significantly improved accuracy and performance. In [26,27], a design of a fractional order controller is proposed for the improved response of the transient state of pendulum cart and gantry crane systems. Other examples in which a fractional control strategy has been used to improve system performance include [28]. By using fractional and integer order sliding mode control [29], a comparison study has been conducted with a planer robotic manipulator and has shown a significant improvement in stability. Similarly, in [30], the same control action has been applied to a nonlinear system based on Lyapunov stability, resulting in the rapid convergence of the system to an equilibrium point. Raouf et al. [31] presented a Lyapunov stability-based serial robotic manipulator whose controller was based on the fractional Proportional–Derivative–Integral (PID) controller. An effort for the control of a 2DOF flexible manipulator using a Linear Quadratic Regulator (LQR)-based controller was shown in [32]. Similarly, various other research works that involve the use of fractional order calculus are available, such as the FO sliding mode control of a Single Flexible Link Robotic Manipulator (SFLRM) [33], the use of a super twisting second-order algorithm ensuring robustness and high positional accuracy [34], and the use of discrete sliding Fourier transform to estimate tip deflection [35]. Fractional order controllers have the potential to improve SFLRM responses tremendously. In [36], the authors presented a fractional sliding mode observer based on the state estimation approach. An algorithm was proposed by Abhaya et al. [37] to choose an appropriate fractional order dynamic model for robotic manipulators. However, to the author’s knowledge, the control of a SFLRM and 2DSFJ using a fractional order PID controller with state estimation techniques is not present in the given literature. Therefore, this paper discusses this method of controller design in order to improve the system performance of robotic manipulators.

This work focuses on the estimation of the states of a 2DSFL and 2DSFJ using output measurements from encoders present at joints of the robotic manipulator. The research pertains to the fractional control and estimation of the orientation of robotic manipulators for better accuracy and low-cost implementation. This improved accuracy has been a major demand in applications such as laser cutting and CNC, where the first choices of sensors are low-cost accelerometers and encoders. The Bayesian state estimation methods such as the Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF), Kalman Filter (KF) and particle filter (PF) have been applied to a 2DOF Serial Flexible Link Robot and 2DOF Serial Flexible Joint Robot by Quanser. The accuracy of the tracking performance and the orientation of the robot are the parameters that have been used for evaluation. The key contributions and highlights in this paper are as follows:

• Utilization of joint angle measurements in the filters for state estimations;

• Implementation of a fractional order PID controller;

• An extensive comparison of filters.

A quick overview of the sections is presented as follows: Section 2 presents discussions and preliminaries of state estimation methods using Kalman Filters. The four basic filters of KF, EKF, UKF and PF are also discussed in this section. Section 3 presents a discussion about fractional calculus and its uses. It also presents the implementation of a fractional order system by approximating with an integer order system that can be implemented very easily. Section 4 shows the dynamic model of the system and a comparison between fractional order PID (FOPID) and integer order PID (IOPID) controllers applied in the simulation. The comparison is performed on the basis of phase and gain margins. The state space analysis is performed for the system model and as an example, the analysis of a link is presented. Similar approaches and functions have been applied for both 2DSFJ and 2DSFL. In Section 4, two distinct robotic manipulators are evaluated to verify the findings of the study: a 2DOF serial robot with flexible links and a 2DOF serial robot with flexible joints. The dynamics of the first flexible link of the two links available to the 2DSFL robotic manipulator are calculated in this article. The full derivation of the dynamics of 2DSFL and 2DSFJ can be found in our book [38]. The simulations results are shown in Section 5, and Section 6 presents concluding remarks.

2. Preliminaries of State Estimation Methods Using Kalman Filter

The Kalman Filter (KF) uses noisy observed data and data with other inconsistencies to estimate unknown states by the use of a mathematical model. It is a standard optimal estimation algorithm based on Bayesian filter theory [39] and was proposed early in the 1960s by R. E. Kalman [40]. As the algorithm offers quick, efficient, real-time and powerful anti-interference estimation, it has been widely used and accepted in the course of target tracking, orbit calculation and navigation. It has also gained importance within the fields of dynamic positioning, microeconomics, sensor fusion and digital image processing. The Extended and Unscented Kalman Filters are other major variants of Kalman Filters. Along with these, we present a comparison study with the particle filter (PF).

2.1. Kalman Filter (KF)

In linear systems, the KF primarily estimates the variables or performs state estimation. It works on dynamic models which have the state (1) and observation (2) equations, where the state equation is a linear combination of the process noise ($w_k$), inputs ($u_{k-1}$) and prior states ($x_{k-1}$), and the latter is combination of states ($x_k$) and measurement noise ($v_k$) [41]. Therefore, the state equation of the dynamic model is derived as follows:

(1)$x_k=A{x}_{k-1}+B{u}_{k-1}+w_k.$

Similarly, the dynamic model’s observation equation is derived as

(2)$y_k=H{x}_k+v_k,$

where the state and observation vectors are represented as ($x_k$, $y_k$), and the state transition matrix is represented by A, the observation matrix is represented by H and the input matrix is represented as B. The process and measurement noise vectors (3), which are represented by $w_k$ and $v_k$, respectively, are assumed to be normally distributed and positive definite:

(3)$q\left(w\right)\sim N(0,Q), q\left(v\right)\sim N(0,R).$

The states obtained during estimation at moments $k-1$ and k are called the prior state estimation (${\hat{x}}_k^{-}$) and posterior state estimation ($x_k$) [42]. Based on these prior and posterior states, we derived prior and posterior estimation errors (4) and the covariance equations of the priori and posterior states (5).

(4)$e_k^{-}=x_k-{\hat{x}}_k^{-}, e_k=x_k-{\hat{x}}_k.$

(5)$P_k^{-}=E\left[e_k^{-}{e}_k^{-T}\right], P_k=E\left[e_k{e}_k^T\right].$

The Kalman Filter equations require $x_k$ such that the posterior state estimation can be calculated based on a priori estimates and weighted differences between real and estimated measurement values, leading to a need for update and prediction equations. Therefore, the formulated prediction equations are given by (6) and (7):

(6)${\hat{x}}_k^{-}=A{\hat{x}}_{k-1}+B{u}_{k-1},$

(7)$P_k^{-}=A{P}_{k-1}{A}^T+Q,$

and, following are the update equations:

(8)$K_k=\frac{P_k^{-}{H}^T}{H{P}_k^{-}{H}^T+R},$

(9)${\hat{x}}_k={\hat{x}}_k^{-}+K_k(z_k-H{\hat{x}}_k^{-}),$

(10)$P_k=(I-K_{k}H)P_k^{-}.$

The Kalman gain matrix of (8), (9) and (10) is represented by $K_k$, and ($x_k$, $P_k,I$) are the optimum state filter vector $x_k$, process covariance matrix $P_k$ and identity matrix I, respectively. There are certain limitations of KF: it is helpful in linear systems only, and the noise should be Gaussian. Since nearly all the real-time systems behave in a nonlinear way, the KF has limited application practically [42].

2.2. Extended Kalman Filter (EKF)

As a well-known nonlinear filter, the EKF is quite simple to implement. The EKF uses Linear Taylor expansion, and when the observation and state equations become continuous and linear, the estimated value converges to the actual value [43]. The performance relies on process and measurement noise matrices, which may cause divergence in filter estimation if these matrices are not estimated accurately enough [44].

Let us begin with a nonlinear system as

(11)$x_k=f(x_{k-1},w_{k-1}),$

(12)$y_k=h(x_k,v_k).$

where $x_k$ is a state vector of the nth dimension and $y_k$ is an observation vector of m-dimensions. The state transition and observation functions f and h are found to be nonlinear. The formulated prediction equations are given by (13) and (14):

(13)$A=\frac{df}{dx}{|}_x={\hat{x}}_{k-1}, x_k^{-}=f\left({\hat{x}}_{k-1}\right),$

(14)$P_k=(I-K_{k}H)P_k^{-}.$

where A is the state matrix (13). Similarly, the update equations are

(15)$H=\frac{dh}{dx}{|}_x={\hat{x}}_{k-1}^{-},$

(16)$K_k=\frac{P_k^{-}{H}^T}{H{P}_k^{-}{H}^T+R},$

(17)${\hat{x}}_k={\hat{x}}_k^{-}+K_k(y_k-h({\hat{x}}_k^{-})),$

where H is the measurement Equation (15), K is the Kalman gain matrix (16), (17), and $P_k$ is the process covariance matrix (14). A case of conflict may arise because the EKF considers the system process and measurement noise to be non-Gaussian in nature, but this may not be the case practically. As a result, this could lead to a change in the process noise and measurement noise experiencing first-order Taylor series expansion; furthermore, the EKF may make assumptions regarding inconsistent noise. Another drawback of the EKF is that it has to recalculate the Jacobian matrices at each operational point, which may amount to a complicated and time-consuming EKF process.

2.3. Unscented Kalman Filter (UKF)

With methods similar to the EKF, the UKF generally has better performance and reduced complexity due to the sampling strategy and avoids the divergence of errors found in the case of EKF. It approximates the weighted density distribution of the nonlinear function and handles non-additive noise without requiring the calculation of Jacobian matrices [42], thus providing greater accuracy and faster convergence for a nonlinear problem.

Julier and Uhlmann devised the UKF and applied a strategy close to nonlinear distribution [45]. The sampling points in the UKF are called sigma points and are generally smaller in number. The 2n + 1 symmetric-sigma sampling strategy is the most widely accepted approach [46], where n is the state vector dimension. These sigma sampling points then have a non-linear transform applied on them to obtain a non-linearly transformed set of points with a mean and covariance of x and $P\left(x\right)$, respectively [47]. The UKF consists of the following processing steps: (i) Initialisation of state error covariance matrix as well as the state vector; (ii) selecting sigma sampling points at moment $k-1$ based on the state vector and error covariance matrix and calculation of weighted values; (iii) updating time through the chosen sampling points and propagation of the covariance and mean through the state equation; (iv) updating measurements by the selected sampling points through an observation equation that is non-linear; and (v) KF coefficients being updated. The state equations for this process are shown in Equations (11) and (12). The selection of sigma points is defined by (18)–(21) as

(18)$x_0={\hat{x}}_k, x_i={\hat{x}}_k+{\left(\sqrt{(n+\Delta )P}\right)}_k, x_{i+n}={\hat{x}}_k-{\left(\sqrt{(n+\Delta )P}\right)}_k, i=1,\dots ,n,$

(19)$W_0^m=\Delta /(\Delta +n), W_i^m=1/2(\Delta +n), i=1,\dots ,2n,$

(20)$W_0^m=W_0^c+(1-{\alpha }^2+\beta ), W_i^c=1/2(\Delta +n), i=1,\dots ,2n,$

(21)$\Delta ={\alpha }^2(k+n)-n,$

where the scaling constant is represented by $\Delta$, the spread of the selected sigma points is given by $\alpha$, k is the secondary scaling constant—mostly considered as zero (21)—$\beta$ incorporates previous information about the state variables distribution, and for the posterior sigma points, $W_i^m$ is the the ith weighted sample mean and $W_i^c$ is the ith weighted covariance (20). The steps for updating the time of the UKF are

(22)${ϵ}_i=f\left(x_i\right),$

(23)${\hat{x}}_{k+1/k}=\sum W_i^m{ϵ}_i,$

(24)$P_{k+1/k}=\sum W_i^c({ϵ}_i-{\hat{x}}_{k+1/k}){({ϵ}_i-{\hat{x}}_{k+1/k})}^T.$

where $ϵ$ is the function representing the non-linear system (22) and ${\hat{x}}_{k+1/k}$ is the prior state estimation (23), (24). Similarly, measurement updates of the UKF can be found according to (25)–(28):

(25)$Z_i=h\left({ϵ}_i\right),$

(26)${\hat{z}}_{k+1/k}=\sum W_i^m{Z}_i,$

(27)$P_{zz}=\sum W_i^c(Z_i-{\hat{Z}}_{k+1/k}){(Z_i-{\hat{Z}}_{k+1/k})}^T,$

(28)$P_{xz}=\sum W_i^c({ϵ}_i-{\hat{x}}_{k+1/k}){({ϵ}_i-{\hat{z}}_{k+1/k})}^T.$

Finally, the filter update of the UKF proceeds as follows:

(29)$K_{k+1}=P_{xz}{P}_{zz}^{-1},$

(30)${\hat{x}}_{k+1}={\hat{x}}_{k+1/k}+K_k(y_{k+1}-{\hat{z}}_{k+1/k}),$

(31)$P_{k+1/k+1}=P_{k+1/k}-K_{k+1}{P}_{zz}{K}_{k+1}^T.$

where K is the Kalman gain (29) and ${\hat{x}}_{k+1}$ is posterior state estimation (30), (31). Since the UKF deals with $2n+1$ sigma points, as well as the sampling and re-sampling of data, it is slower than the EKF.

2.4. Particle Filter (PF)

When processes and systems are non-Gaussian, the particle filter proves to be a very efficient predictive tool to represent and manage uncertainty and take account of the stochasticity of the process. The PF allows information to be processed from multiple measurement sources that are fused logically. However, there are a few drawbacks; e.g., the prediction results in PF greatly depend upon the number of particles or samples. One is more likely to obtain a close prediction of the estimated value to the actual value if there is a greater number of particles.

Let the discrete state space model design be

(32)$\begin{array}{ccc}\hfill P_{t+1}& =& g(x_t,u_t,w_t),\hfill \\ \hfill y_t& =& f\left(x_t\right)+e_t.\hfill \end{array}$

The state vector is given by $x_t\in {\mathbb{R}}^n$, and $u_t$ and $\mathbb{Y}={\left\{y_i\right\}}_{i=1}^t$ are the input vector and measurement vector, respectively. With probability density functions (PDFs) for process $p_w\left(w\right)$ and measurement noises $e_t$, and for the Bayesian inference [48], the filtering density given by $p\left(x_t\right|{\mathbb{Y}}_t)$, as well as the non-linear posterior prediction density represented as $p\left(x_{t+1}\right|{\mathbb{Y}}_t)$, are given as

(33)$x_{t+1}|{\mathbb{Y}}_t={\int }_{{\mathbb{R}}^n}p\left(x_{t+1}\right|x_t)p\left(x_t\right|{\mathbb{Y}}_t)d{x}_t,$

(34)$p\left(x_t\right|{\mathbb{Y}}_t)=\frac{p\left(y_t\right|x_t)p\left(x_t\right|{\mathbb{Y}}_{t-1})}{p\left(y_t\right|{\mathbb{Y}}_{t-1})}.$

An approximated solution for Equations (33) and (34) (the discrete-time Bayesian estimation problem) has been provided by Doucet et al. [5] and Gordon et al. [49] through a constant update in the posterior filtering density. The weight and location assigned to each particle represents the density in the state-space region because, unlike UKF, the particle filter estimates the filtering density $p\left(x_t\right|{\mathbb{Y}}_t)$ with a large set of weighted (${\psi }_t^{\left(i\right)}$) samples (particles) ${\{x_i^{\left(i\right)}\}}_{i=1}^N$ such that the sum of all weights is unity. The locations and corresponding weights of each particle are updated with each new measurement, and the corresponding weights and particles then can be used to solve Bayesian equations approximately; the re-sampling of samples is performed to avoid divergence [49]. The steps involved in the calculation of PF are (i) the generation of N samples ${\{x_i^{\left(i\right)}\}}_{i=1}^N$ from $p\left(x_0\right)$; (ii) computation of weights $\psi {\left(i\right)}_t=\psi {\left(i\right)}_{t-1}.\frac{p\left(y_t\right|x_t^{\left(i\right)})p\left(x_t^{\left(i\right)}\right|x_{t-1}^{\left(i\right)})}{p^{prop}\left(x_t^{\left(i\right)}\right|x_{t-1}^{\left(i\right)},y_t)}$, and normalisation of weights ${\overline{\psi }}_t^{\left(i\right)}={\psi }_t^{\left(i\right)}/{\sum }_{j=1}^N{\psi }_t^{\left(j\right)},i=1,\dots ,N$; (iii) generation of predictions from he proposed distribution $X_{t+1}^{\left(i\right)}\sim p^{prop}(x_{t+1}^{\left(i\right)}|x_t^{(i\ast )},y_{t+1}),i=1,\dots ,N$; and (iv) increment of t and repeated performance of step (ii).

The minimum mean square estimate (MMSE) (35) is often chosen to estimate each time t; that is,

(35)${\hat{x}}_{t|t}=\mathbb{E}\left(x_t\right)={\int }_{{\mathbb{R}}^n}{x}_{t}p\left(x_t\right|{\mathbb{Y}}_t)d{x}_t\approx \sum _{i=1}^N{\psi }_t^{\left(i\right)}{X}_t^{\left(i\right)},$

but there are also other available choices of estimation, such as the ML-estimate. From [5], we conclude that $N\to \infty$, meaning that the estimated PDF converges into the true value.

3. Preliminaries of Definitions of Fractional Calculus

From the literature review, it can be concluded that fractional order calculus represents a general version of ordinary integer order calculus. With fractional calculus, an operator D, of order h, needs to be established to generalise the ordinary concept of derivation (positive h) and integration (negative h), where h is not an integer number [50]. Riemann and Liouville [51] presented the most accepted definition of fractional derivatives which generalise the given relation with respect to integer orders as represented in (36),

(36)${}_0{D}_x^{-p}={\int }_x^c\frac{{(x-t)}^{p-1}}{(p-1)!}f\left(t\right)dt, p\in \mathbb{N}.$

For fractional order calculus, the general form of D becomes ${}_c{D}_x^{h}f\left(x\right)$. The operator D becomes $L\left[{}_0{D}_x^h\right]=s^{h}F\left(s\right)$ when applying a Laplace transform ($\mathcal{L}$) with zero initial conditions. The approximation to the integer order functions is the most usual way to implement such transfer functions, simulations and fractional order equations. An FO transfer function can be perfectly approximated to IO transfer function, if the IO transfer function has zeros (z) and poles (p) such that $(z,p)\to \infty$. However, the approximations are logically possible only with limited poles and zeros. To overcome this, Oustaloup discovered one of the best-known approximations, which uses a recursive distribution of poles and zeros:

(37)$s^v\approx a\prod _{n=1}^N\frac{1+(s/{\Omega }_{z,n})}{1+(s/{\Omega }_{p,n})}.$

The legitimacy of this approximation lies within the frequencies $[{\Omega }_l,{\Omega }_h]$. The gain (a) is chosen such that both sides of Equation (37) have a unit gain of 1 rad/s.

The desired performance of the system is highly dependent upon the constant (N), which is the order of transfer functions. Lower values of N are responsible for ripples in both phase and gain behaviours but are simpler to approximate. Conversely, the ripples can be removed by choosing higher values of N, but this may cause computationally heavier approximations. The poles and zeros obtained from Equation (37) have frequencies which are given by (38)–(40):

(38)${\Omega }_{z,1}={\Omega }_l\sqrt{ϵ},$

(39)${\Omega }_{p,n}={\Omega }_{(z,n)}\times \delta , n=1,\dots ,N,$

(40)${\Omega }_{z,n+1}={\Omega }_{(p,n)}\times \delta , n=1,\dots ,N-1,$

where $\delta ={({\Omega }_h/{\Omega }_l)}^{h/N}$ and $ϵ={({\Omega }_h/{\Omega }_l)}^{(1-h)/N}$. These equations are satisfactory for $h<1$ but are inconsistent for $h>1$. Thus, it is usual to isolate the fractional and integer orders of s (41):

(41)$s^h=s^i{s}^f, h=i+f, i\in \mathbb{Z},f\in [0,1].$

Therefore, only the term $s^f$ is to be approximated. The approximation in Equation (37) can be discretised if the case of a discrete fractional order is sought [52].

4. Dynamic Model

Let us consider a single link robotic manipulator with parameters $\theta$ and $\varphi$ that represent the angular position and angular shift of the link. The manipulator model has been taken from the course-ware manual by Quanser Inc. It is a flexible link which is able to rotate “to” and “fro” on its pivot. The mass of the link is m = 0.07 kg, with a length $l=0.42$ m, damping coefficient $B_2=0.004$ Nm/(rad/s) and link inertia $J_2=2.08\times {10}^{-3}$ kgm${}^2$ [19]. However, before considering the equations of the dynamic model, let us discuss the applied FO controller and its approximated version.

4.1. Fractional Order PID Controller

A fractional order controller ($P{I}^{\lambda }{D}^{\delta }$) is derived by the differential equation

(42)$u\left(t\right)=k_{p}e\left(t\right)+k_i{D}_t^{-\lambda }e\left(t\right)+k_d{D}_t^{\delta }e\left(t\right),$

where $e\left(t\right)$ is the error term. The continuous state transfer function is obtained using the Laplace transform ($\mathcal{L}$) on Equation (42):

(43)$G_c\left(s\right)=k_p+k_I{s}^{-\lambda }+k_D{s}^{\delta }.$

As shown in Equation (43), the design of the FOPID controller requires three gain parameters ($k_p,k_I,k_D$) and two order values ($\lambda ,\delta$). As mentioned earlier, the FO controller provides better flexibility and adaptability to achieving control objectives as it is a generalised form of traditional Integer order calculus. Clearly, $(\lambda =1),(\delta =1)$ is a case of an integer order controller.

The simulations on the 2DSFJ and 2DSFL controlled by the IOPID controller and FOPID controller are carried out on basis of the parameters provided in Section 5 with the limitation of maximum torque and input current (i). As stated earlier, N is the order of approximation on which the desired performance of the approximated transfer function strongly relies. As mentioned earlier, lower values of N cause ripples in both phase and gain behaviours but are simpler to approximate, and vice versa; this can be removed by increasing the order, leading to computationally heavier approximation. The fractional order transfer function of the controller is then approximated using an Outsaloup approximation of the integer order transfer function. It turns out that the order of the integer order transfer function is 15. On converting the transfer function into state space notation, the state matrix is of the order $15\times 15$. In order to reduce complexity, the order of the obtained transfer function is reduced using the inbuilt MATLAB model reducer toolbox. The “truncate” state elimination method is used in a frequency range of $({\omega }_l=0.01,{\omega }_h=10)$ with an absolute tolerance of $1\times {10}^{-6}$. The order of the reduced transfer function is chosen to be $4$.

The transfer function obtained after the reduction is given as

(44)$G_c\left(s\right)=\frac{2.131s^4+2.08s^3+0.2952s^2+0.0046s-6.29\times {10}^{-05}}{s^4+0.8902s^3+0.1049s^2+0.001004s-1.797\times {10}^{-04}}.$

Figure 2 and Figure 3 show the frequency response for the fractional order controller, the reduced version of the FO controller and the integer order controller. As our working limit is in frequency range of [0.01, 10], form Figure 2, it can be stated that the reduced transfer function of the FOPID and the original transfer function of the FOPID have nearly the same responses. Both provide significant gain and phase margins to the plant model, whereas in the case of the integer order PID controller, from Figure 3, it can be stated that there are insufficient margins in the desired frequency range of [0.01, 10]. Further comparisons and discussions are shown in Section 5.

4.2. State Space Model of 2DSFL and 2DSFJ

The dynamic model and equations of the system have been derived with the help of Euler–Lagrange formulation. The equations were then linearised using the Taylor expansion series to achieve representation in the state-space:

(45)$\dot{x}=Ax+Bu,$

$y=Cx+Du.$

The schematic diagram and the mechanical model equivalent are shown in Figure 4 and Figure 5, respectively.

The Lagrangian equation is applied to derive the linear state-space model. As we know, the Lagrangian is the difference between the total kinetic (KE) and potential energies (PE). To derive the equations in the form of (45), the following procedure was carried out. Using mechanics, we define the relations

(46)$KE=\frac{1}{2}{J}_1{\dot{\theta }}^2+\frac{1}{2}{J}_2{(\dot{\theta }+\dot{\varphi })}^2, PE=\frac{1}{2}{K}_s{\varphi }^2,$

where the constants $J_{1}and{K}_s$ have been explained in Table 1. From Equation (46), the Lagrangian L is derived as

(47)$L=Kinetic energy-Potential energy=\frac{1}{2}{J}_1{\dot{\theta }}^2+\frac{1}{2}{J}_2{(\dot{\theta }+\dot{\varphi })}^2-\frac{1}{2}{K}_s{\varphi }^2.$

Using Equation (47) and according to the Euler–Lagrangian formulation, we have (48):

(48)$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta }}-\frac{\partial L}{\partial \theta }=\tau -B_2\dot{\theta }, \frac{d}{dt}\frac{\partial L}{\partial \dot{\varphi }}-\frac{\partial L}{\partial \varphi }=0.$

From (47), we obtain the governing set of EOMs (equations of motion) as

(49)$(J_1+J_2)\ddot{\theta }+J_2\ddot{\varphi }+B_2\dot{\theta }=\tau ,$

(50)$J_1\ddot{\theta }+J_2\ddot{\varphi }+K_s\varphi =0.$

Differentiating the values of $\ddot{\theta }$ and $\ddot{\varphi }$ in (49) and (50), we obtain (51):

(51)$\ddot{\theta }=-\frac{B_1}{J_1}\dot{\theta }+\frac{K_s}{J_2}\dot{\varphi }+\frac{1}{J_1}\tau , \ddot{\varphi }=\frac{B_1}{J_1}\dot{\theta }+\frac{K_s(J_1+J_2)}{J_1{J}_2}\dot{\varphi }-\frac{1}{J_1}\tau .$

The linear state space model achieved from the above equations for a single flexible link is represented in Equation (52). Similarly, the state matrices and input matrices are found for Link 2. Using the same analogy, we have defined the state-space matrices for joints in Equation (53) as

(52)$A_l=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& \frac{K_s}{J_1}& -\frac{B_1}{J_1}& \frac{B_2}{J_1}\\ 0& -\frac{(J_1+J_2)K_s}{J_1{J}_2}& \frac{B_1}{J_1}& -\frac{(J_1+J_2)B_2}{J_1{J}_2}\end{array}\right], B_l=\left[\begin{array}{c}0\\ 0\\ \frac{K_t}{J_1}\\ -\frac{K_t}{J_1}\end{array}\right].$

(53)$A_j=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -\frac{K_s}{J_1}& \frac{K_s}{J_1}& -\frac{B_1}{J_1}& 0\\ \frac{K_s}{J_2}& -\frac{K_s}{J_2}& 0& -\frac{B_2}{J_2}\end{array}\right], B_j=\left[\begin{array}{c}0\\ 0\\ \frac{K_t}{J_1}\\ 0\end{array}\right].$

After the state-space model, various simulink (MATLAB) models have been designed for the following purposes: (i) first, the comparison of integer order and fractional order responses; (ii) second, the application of the Kalman Filter to remove noise and the prediction of state for the feedback purpose. Similar state space models have been formed for the remaining link and flexible joints. The physical parameters are the only difference between the state space model presented here and that of the remaining link and flexible joints [53,54].

5. Simulation and Results

Figure 6 shows the basic Simulink model that has been designed for FOPID and IOPID controllers.

The values of the constant chosen for the fractional order PID and integer order PID controllers are shown in Table 2, where $F_{lower}$ and $F_{upper}$ are the Oustaloup filter frequency range and the order is the order of approximation. The input was provided by a square wave generator and then fed to the state-space model with their respective controllers. From Table 2, it can be stated that a low value of actuator effort is required and better fits the system. However, from Table 3, it can also be concluded that lower actuator effort leads to a higher rise time. The parameters of the controllers are tuned by the tool available in Simulink to get the best response. The performance in terms of the rise time of the IO controller is better compared to the FO controller, but the FO controller outperforms the other for other parameters.

The response comparison for IO and FO systems for both the link and joint is shown in Figure 7, Figure 8, Figure 9 and Figure 10, indicating the improvement of performance parameters such as rise time, settling time and peak overshoot with the FOPID controller.

After implementing the FOPID controller, a noise signal arbitrarily introduced in the range of ${10}^{-4}$ helped to mimic real-world noise. Various types of Kalman fFlters have been applied and tested to ensure feedback with minimum noise. The filter parameters were chosen using the trial and error method to get the best responses. The figures below show the performance of the filters when applied. The responses are shown for both SFL and SJF, but only for Link 1 and Joint 1. The responses for Link 2 and Joint 2 are similar to their parent links and joint as the system considered is a cascaded system. Figure 11 and Figure 12 present the main response and its enlarged image for better visibility. A square wave input is given to the model.

Figure 11 and Figure 12 show that the state estimation methods as applied to integer order and fractional order controller systems provide better performance than the PF. The fractional order controller systems for both the UKF and PF provide better adaptability and manoeuvering. These contribute to lower actuator effort with better performances.

Figure 13 and Figure 14 present the comparison of all the filters studied in this manuscript for Link 1 and Joint 1, and the performance of filters can be analyzed from these. The number of particles used in the PF is 100. The PF provides a significantly improved performance than the other approaches.

As we know, the distribution considered in the EKF and UKF is Gaussian, and so the same is considered in the PF for a fair comparison. The filter matrices were chosen to be the same as the system state matrices without considering the noise. It can be concluded that the particle filter provides better estimations. Some recent and advanced filters can have a better response than the particle filter, but among all the filters compared, the PF provides a better response.

Quantitatively, the performance may vary for different control techniques in the case of IOPID and FOPID controllers (Figure 11 and Figure 12), but qualitatively, the PF has a better response. This is because the PF does not generalize all of the data; rather, it operates on the areas of data that are known as particles, which are both time-consuming and computationally intensive processes. A higher number of particles does not necessarily lead to a better performance, because two consecutive data points may have similar characteristics and they need to be operated together. Therefore, the number of particles is an important factor. In comparison to the PF, the UKF has a smaller number of particles, also known as sigma points. Another factor that leads to the better performance of the PF is the distribution of data. The KF, EKF and UKF consider the distribution of data as Gaussian, but there may be a case in which data is not distributed in a Gaussian form. The PF provides the flexibility to choose the type of distribution that better fits the data points. We have provided a review and results of different techniques for state estimation with the application of FO calculus. Overall, the combined algorithm offers improved performance for the 2DSFL and 2DSFJ concerning conventional control in terms of parameters such as tracking accuracy and the orientation of the manipulator.

6. Conclusions

The performance and merits–demerits of the KF, EKF, UKF and PF along with the implementation of a fractional order controller have been presented in this work. The methods have been validated by simulations in MATLAB. This research has focused on the improved performance of SFL and SFJ with respect to the conventional PID controller. The performance is mainly focused on the smooth and quick response of the model along with the noise removal, which can be validated form the responses shown in the graphs. In order to do this, fractional order controller and state estimation methods have been used for the following reasons: (i) a fractional order depicts the system more accurately, and (ii) state estimation methods provide better tracking of signals from noisy measurements. Therefore, it can be concluded that the combination of both fractional calculus and state estimation methods provides more robustness, performance and stability to the system than the conventional PID controller for our 2DSFL and 2DSFJ system.

Author Contributions

Conceptualization, A.P.S. and D.D.; methodology, S.G. and D.D.; software, S.G.; validation, S.G.; formal analysis, A.P.S. and S.G.; writing—original draft preparation, S.G. and D.D.; writing—review and editing, S.G., A.P.S., D.D. and S.O.; supervision, D.D.; funding acquisition, S.O. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the European Regional Development Fund in “A Research Platform focused on Industry 4.0 and Robotics in Ostrava Agglomeration” project, Reg. No. CZ.02.1.01/0.0/0.0/17_049/0008425 within the Operational Programme of Research, Development and Education.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Classic control strategy based on state estimation methods.

Figure 2. Bode plot for FO and reduced FO transfer function of controller for Link 1.

Figure 3. Bode plot for IO transfer function of controller for Link 1.

Figure 4. Single Flexible link schematic.

Figure 5. Single Flexible link mechanical model.

Figure 6. Basic process flow for 2DSFL and 2DSFJ, where θ: angular position, ϕ: angular shift, θ˙: angular velocity, ϕ˙: rate of angular shift.

Figure 7. Controller responses for angular position (θ) for Link 1.

Figure 8. Controller responses for angular shift (ϕ) for Link 1.

Figure 9. Controller responses for angular position (θ) for Joint 1.

Figure 10. Controller responses for angular shift (ϕ) for Joint 1.

Figure 11. Integer and fractional order controller comparison for different state estimation methods for Link 1.

Figure 12. Integer and fractional order controller comparison for different state estimation methods for Joint 1.

Figure 13. Comparison of filter responses with fractional order control for Link 1.

Figure 14. Comparison of filter responses with fractional order control for Joint 1.

Table 1

Parameters of robotic manipulator for link 1.

Symbol	Description	Units
$I_m$	Motor current for first link	A
$K_t$	Torque constant for first motor drive	Nm/A
T	Torque produced by Drive 1	Nm
$\theta$	Angular position of Link 1 drive shaft	rad
$\frac{d\theta }{dt}$	Angular velocity of Link 1 drive shaft	rad/s
$\varphi$	Link 1 Relative end-effector angular position	rad
$\frac{d\varphi }{dt}$	Link 1 relative end-effector angular velocity	rad/s
$J_1$	First joint equivalent moment of inertia	kg m${}^2$
$B_1$	First joint equivalent viscous damping coefficient	Nm s/rad
$J_2$	Link 1 end-effector equivalent moment of inertia (combined with second stage)	kg m${}^2$
$B_2$	Link 1 end-effector equivalent damping coefficient (combined with second stage)	Nm s/rad
$K_s$	First link torsional stiffness constant	Nm/rad

Table 2

PID parameters.

Controller Type	Parameter	Value	Controller Type	Parameter	Value
Integer order	P	6.94	Fractional order$F_{lower}$ = 0.01$F_{upper}$ = 10Order = 2	P	1.9
	I	0.9		I	0.3
	D	0.11		$\lambda$	0.4
	N	27.25		D	0.05
				$\mu$	0.3

Table 3

Performance parameters of Link 1 and Joint 1.

Controller →	Link 1			Joint 1
Parameter ↓	IO	FO	Improvement	IO	FO	Improvement
Settling time (in s) for 4% error	3.555	1.952	40.53%	3.194	2.419	24.26%
Peak overshoot (in %)	12.928	3.870	70.45%	53.8	8.124	84.90%
Rise time (5–95%) (in ms)	201.658	282.379	−28.57%	242.853	500.818	−51.51%