Optimal Backstepping-FOPID Controller Design for Wheeled Mobile Robot

2.1 Modeling the wheeled mobile robot

The WMR in Figure 1 contains two pairs of DC motors with driving wheels of radius R, where φL and φR are the left and right rotating angles of wheels, respectively. θ is the robot orientation and L is the width of the robot body. At the distance d from the mind-point A, where (Xa, Ya) is the coordinate of A in the inertial frame (X,Y), and the coordinates of any point in the robot frame are defined by (Xr, Yr).

Three steps are necessary to obtain the mathematical model of the robot mobile: kinematic modeling, dynamic modeling, and actuator modeling, which are defined as follows: 

2.1.1 Kinematic model

This section aims to define a relationship between the linear and angular speeds of the mechanical systems without taking into account the forces affecting the motion [10, 11]. The linear speed of each driving wheel is:

$\left\{\begin{array}{l}\mathcal{v}_{R}=R \dot{\varphi}_{R} \\ \mathcal{v}_{L}=R \dot{\varphi}_{R}\end{array}\right.$        (1)

1.png

Figure 1. Schematic design of the WMR

The linear and angular velocities of the WMR are given by Eqns. (2) and (3), respectively:

$v=\frac{v_{R}+v_{L}}{2}=R \frac{\left(\dot{\varphi}_{R}+\dot{\varphi}_{L}\right)}{2}$         (2)

$\omega=\frac{v_{R}-v_{L}}{2 L}=R \frac{\left(\dot{\varphi}_{R}-\dot{\varphi}_{L}\right)}{2 L}$       (3)

The kinematic constraint can express by the following equations [36]:

No slip constraint:

$-\dot{x}_{a} \sin \theta+\dot{y}_{a} \cos \theta=0$         (4)

Pur rolling constraint:

$\dot{x}_{a} \cos \theta+\dot{y}_{a} \sin \theta+L \dot{\theta}=R \dot{\varphi}_{R}$

$\dot{x}_{a} \cos \theta+\dot{y}_{a} \sin \theta-L \dot{\theta}=R \dot{\varphi}_{L}$        (5)

The three constraint equations are:

$\Lambda(q) \dot{q}=0$        (6)

where:

$(q)=\left[\begin{array}{ccccc}-\sin \theta & \cos \theta & 0 & 0 & 0 \\ \cos \theta & \sin \theta & L & -R & 0 \\ \cos \theta & \sin \theta & -L & 0 & -R\end{array}\right]$       （7）

and

$\dot{q}=\left[\dot{x}_{a} \dot{Y}_{a} \dot{\theta}_{R} \dot{\varphi}_{L}\right]^{T}$       （8）

So, the kinematic model obtained is:

$\left[\begin{array}{c}\dot{x}_{a} \\ \dot{y}_{a} \\ \dot{\theta} \\ \dot{\varphi}_{R} \\ \dot{\varphi}_{L}\end{array}\right]=\left[\begin{array}{cc}\cos \theta & 0 \\ \sin \theta & 0 \\ 0 & 1 \\ \frac{1}{R} & \frac{-L}{R} \\ \frac{1}{R} & \frac{-L}{R}\end{array}\right]\left[\begin{array}{c}v \\ \omega\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}\cos \theta & 0 \\ \sin \theta & 0 \\ 0 & 1 \\ \frac{1}{R} & \frac{-L}{R} \\ \frac{1}{R} & \frac{-L}{R}\end{array}\right]\left[\begin{array}{c}\dot{\varphi}_{R} \\ \dot{\varphi}_{L}\end{array}\right]$        (9)

which may be written as:

$\dot{q}=S(q)_{\eta}$       (10)

where, $n=\left[\begin{array}{c}\dot{\varphi}_{\mathrm{R}} \\ \dot{\varphi}_{\mathrm{L}}\end{array}\right]$ is the vector of the angular velocities of two wheels.

2.1.2 Dynamical model

The purpose of dynamic modeling is to take into account all the different forces and energies that affect the mechanical system motion. The motion equation of the WMRs is given by [36, 37]:

$\begin{aligned} M(q) \ddot{q}+V(q, \dot{q}) \dot{q}+F(\dot{q})+G(q)+\tau_{d} \\=& B(q) \tau-\Lambda^{T}(q)^{\lambda} \end{aligned}$       (11)

For simulation purposes and control, Eq. (11) should be transformed into an alternative form, using the kinematic model (10):

$S^{T}(q) \Lambda^{T}(q)=0$        (12)

The new matrices are express as follows:

$\left\{\begin{array}{c}\bar{M}(q)=S^{T}(q) M(q) S(q) \\ \bar{V}=S^{T}(q) M(q) \dot{S}(q)+S^{T}(q) V(q, \dot{q}) S(q) \\ \bar{B}=S^{T}(q) B(q)\end{array}\right.$         (13)

The reduced form of the dynamic equations is expressed as:

$\{\bar{M}(q) \dot{\eta}+\bar{V}(q, \dot{q}) \eta=\bar{B}(q) \tau$        (14)

where,

$\bar{M}(q)=\left[\begin{array}{cc}I_{w}+\frac{R^{2}}{4 L^{2}}\left(m L^{2}+I\right) & \frac{R^{2}}{4 L^{2}}\left(m L^{2}-I\right) \\ \frac{R^{2}}{4 L^{2}}\left(m L^{2}-I\right) & I_{w}+\frac{R^{2}}{4 L^{2}}\left(m L^{2}+I\right)\end{array}\right]$        (14a)

and

$\begin{aligned} \bar{V}(q, \dot{q})=& {\left[\begin{array}{cc}0 & \frac{R^{2}}{2 L} m_{c} d \dot{\theta} \\ -\frac{R^{2}}{2 L} m_{c} d \dot{\theta} & 0\end{array}\right] }, \\ \bar{B}(q)=\left[\begin{array}{ll}1 & 0 \\ 0 & 1\end{array}\right] \end{aligned}$        (14b)

Eq. (14) may be rewritten in a compact form:

$\left\{\begin{array}{c}\left(m+\frac{2 I_{w}}{R^{2}}\right) \dot{v}-m_{c} d \omega^{2}=\frac{1}{R}\left(\tau_{R}+\tau_{L}\right) \\ \left(I+\frac{2 L^{2}}{R^{2}} I_{w}\right) \omega+m_{c} d \omega v=\frac{L}{R}\left(\tau_{R}-\tau_{L}\right)\end{array}\right.$        (15)

where, $m=m_{c}+2 m_{w}$ is the total mass of the robot, and $I=I_{c}+m_{c} d^{2}+2 m_{w} L^{2}+2 I_{m}$ represents the total equivalent inertia.

2.1.3 Actuator modeling

Two pairs of dc motors are considered as actuators for producing the torque to control the inputs for driving the wheels of the WMR as presented by Figure 2. The dynamic model of the actuators can be represented as [38]:

$\left\{\begin{aligned} v_{a}=R_{a} i_{a}(t) &+L_{a} \frac{d i_{a}(t)}{d t}+e_{a}(t) \\ e_{a}(t) &=K_{b} \omega_{m}(t) \\ \tau_{m}=j \frac{d w_{m}(t)}{d t} &+f w_{m}(t)+K_{t} i_{a}(t) \\ \tau &=N \tau_{m} \end{aligned}\right.$        (16)

where, ia is the armature current, (Ra, La) is the resistance and inductance of the armature winding, respectively, ea is the back emf, ωm is the rotor angular speed, τm is the motor torque, (Kt, Kb) are the torque constant and back emf constant, respectively, N is the gear ratio, and τ is the output torque applied to the wheel Since the robot motors are mechanically coupled to wheels through the gears.

2.png

Figure 2. Equivalent electrical schema of the dc motor

Therefore, each dc motor will have:

$\left\{\begin{aligned} \omega_{m R}=N \dot{\varphi}_{w R} & & \text { and } \tau_{R}=N \tau_{m R} \\ \omega_{m L}=N \dot{\varphi}_{w L} & & \text { and } \tau_{L}=N \tau_{m L} \end{aligned}\right.$           (17)

For the two motors, the dynamic model is expressed as:

$\begin{cases}\frac{1}{\left(R+L_{a} P\right)}\left(e_{a r}-K_{b} N \dot{\varphi}_{w R}\right)=i_{R} & \text { with } \tau_{R}=N k_{t} i_{R} \\ \frac{1}{\left(R+L_{a} P\right)}\left(e_{a l}-K_{b} N \dot{\varphi}_{w L}\right)=i_{L} & \text { with } \tau_{L}=N k_{t} i_{L}\end{cases}$      (18)

The robot physical parameters are provided in Table 1 [38].

Table 1. Physical parameter of the robot

2.2 Controller design and strategy

Figure 3 presents the schematic design of the control strategy based on the dynamic and kinematic of the wheeled mobile robot. The trajectory generator produces the reference coordinate x, y and θ, a back-stepping kinematic controller is located in the external loop, here the controller experiences the difference between the reference data received from the bloc trajectory generator and actual value from the robot then generate his own angular and linear speeds that are sent into the internal loop where two blocs of FOPID controller are used as a dynamic controller. Here the two FOPID controllers yield their own signals of angular and linear speeds based on the reference values obtained from the kinematic controller.

3.png

One of the goals here is to realize path tracking. The back-stepping technique has been widely used and is well known as a stable tracking control rule which makes it adopted for this purpose [11, 14]. The controller structure is highlighted in Figure 3, where the input error and velocity vector (v_c) are:

$\left[\begin{array}{l}e_{x} \\ e_{y} \\ e_{\theta}\end{array}\right]=\left[\begin{array}{ccc}\cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1\end{array}\right]\left[\begin{array}{cc}X_{r}- & X \\ Y_{r}- & Y \\ \theta_{r}- & \theta\end{array}\right]=T_{e} e_{r}$        (19)

$v_{c}=\left[\begin{array}{c}v_{r} \cos \left(e_{\theta}\right)+k_{x} e_{x} \\ w_{r}+k_{y} v_{r} e_{y}+k_{\theta} v_{r} \sin \left(e_{\theta}\right)\end{array}\right\rfloor$         (20)

where, k_x, k_y and k_θ are the tuning parameters.

2.2.2 Dynamic controller description

The fractional proportional integral derivative FOPID was first proposed in 1999 by Podlubny [39]. Introducing the two new fractional components λ and μ to the classic PID controller under the name of Fractional integrator and differentiator respectively, the operator of non-integer-order is provided in Eq. (21) [40].

$D_{n}^{t}= \begin{cases}\frac{d^{n}}{d t^{n}} & n>0 \\ 1 & n=0 \\ \int_{t}^{a} d t^{n} & n<0\end{cases}$        (21)

where, a and t are the lower limit and upper limit of the process and $n \in \mathbb{R}$  the constant integral differential operator.

The fractional calculus has three main definitions Riemann Liouville (RL), Grünwald Letnikov (GL), and Caputo. The most common definition utilized in engineering applications is the Caputo method [30, 40] given in Eq. (22).

$D_{n}^{t} f(t)=\frac{1}{\Gamma(\alpha-n)} \int_{a}^{t} \frac{f^{n}(\tau)}{(t-\tau)^{n-\alpha+1}} d \tau$       (22)

for $n-1 \leq n \leq \alpha$ . The term s^α doesn’t have an analytical solution, but it has a fractional order, which is hard to implement. Therefore, numerical solutions such as Oustaloup approximation are adopted [40].

$s^{\alpha} \approx K \prod_{n=-N}^{N} \frac{1+\frac{s}{\omega_{z, n}}}{1+\frac{s}{\omega_{p, n}}} \quad \alpha>0$       (23)

where, N is the number of poles/zeros. Figure 4 presents the parallel structure of the FOPID controller block diagram. Where E(S) distinguishes for the input and errors (E(S) and U(S)).

4.png

Figure 4. Block diagram of fractional order

5.png

Figure 5. Generalized FOPID controller

The representation of the transfer function for the Fractional Order PID (FOPID) in the time domain (Eq. (24)) and the transfer function corresponds to the fractional order in the Laplace domain (Eq. (25)).

$u(t)=K_{p} e(t)+k_{i} D_{t}^{-\lambda} e(t)+K_{d} D_{t}^{\mu} e(t)$         (24)

$C(s)=\frac{U(s)}{E(s)}=K_{p}+\frac{K_{i}}{s^{\lambda}}+K_{d} S^{\mu}$         (25)

where, K_p, K_i and K_d are proportional, integral and derivative gains constants, respectively, $\lambda$ and $\mu$ are factional order of the integral and derivative term. All type of the classic PID may be obtained by taking $(\lambda, \mu)=\{(1,1),(1,0),(0,1),(0,0)\}$, as shown in Figure 5.

2.2.3 Tuning parameters for FOPID

The WMR block diagram for tuning parameters of the Two FOPID controllers is highlighted in Figure 6. The first controller receives the difference between the reference and actual velocity from the right wheel as an input where the desired velocity is constant: Ur =1 m/s.

6.png

Figure 6. Bloc diagram for parameters tuning

The second controller in the left wheel experiences the difference between the actual and desired orientations as input. Where the reference angle was represented by a constant value: θ_r = 0.785 rad. The following Eq. (26) represent the relation between the input voltage of the right and left DC actuators, respectively U_R, U_L with the outputs of the first and second controller, respectively U_V, U_θ:

$\left\{\begin{array}{l}U_{R}=\frac{U_{V}+U_{\theta}}{2} \\ U_{L}=\frac{U_{V}-U_{\theta}}{2}\end{array}\right.$       (26)

The cost function used in the study is the integral square error (ISE) which is demonstrated by Eq. (27):

$I S E=\left(\int_{0}^{\infty}\left[e_{v}(t)\right]^{2} d t+\int_{0}^{\infty}\left[e_{\theta}(t)\right]^{2} d t\right)$       (27)

where:

$\left\{\begin{array}{l}e_{v}=U_{r}-U_{m} \\ e_{\theta}=\theta_{r}-\theta_{m}\end{array}\right.$       (28)

U_r is the desired velocity, U_m the actual velocity, e_v the velocity error, θ_r the desired orientation, θ_m the measured orientation, and e_θ is the orientation error. The optimization algorithms used to tune the parameters (K_p, K_I, K_d, λ, µ) are described in the next section.

2.3 Meta-heuristic algorithms

2.3.1 Particle swarm optimization (PSO)

In this method and for each iteration, all of the particles in the swarm update their positions with the use of the following equations [34, 41]:

$\begin{aligned} V_{i}^{k+1}=w V_{i}^{k}+C_{1} & R_{1}\left(p_{i}^{k}+x_{i}^{k}\right) \\ &+C_{2} R_{2}\left(g b e s t_{i}+x_{i}^{k}\right) \end{aligned}$         (29)

$X_{i}^{k+1}=X_{i}^{k}+V_{i}^{k+1}$        (30)

Here, R₁ and R₂ are random numbers in the range [0, 1]. i refer to the particle in the swarm. k is the iteration step carried out, while w is the inertia weight parameter. The coefficients C₁ and C₂ are the optimization parameters, X: position vector, and $p_{i}^{k}$: best position information achieved by the i^th particle, gbest $_{i}$: best position information available in the swarm [33].

2.3.2 Grey wolf optimizer (GWO)

The mathematical models of GWO consist of social hierarchy, enriching prey, search for prey, attacking prey, and hunting [42-44]. Social hierarchy: The leader of wolves is called alpha $(\alpha)$, Beta $(\beta)$ and delta $(\delta)$ are the second and third levels in the group, respectively.

$d=\left|c \cdot x_{p}(n)-x(n)\right|$       (31)

$(n+1)=x_{p}(n)-a \cdot d$        (32)

where, x_p: position vector of the prey, n: current iteration, and x: position vector of a grey wolf. The vectors a and c are mathematically formulated as follows:

$\vec{a}_{(.)}=2 \vec{l} \cdot \vec{r}_{1}-\vec{l}$        (33)

$\vec{c}_{(.)}=2 \cdot \vec{r}_{2}$       (34)

where, r₁ and r₂ are random numbers in [0, 1]. The following equations for hunting are used:

$\left\{\begin{array}{l}\vec{d}_{\alpha}=\left|\vec{c}_{1} \cdot \vec{x}_{\alpha}-\vec{x}\right| \\ \vec{d}_{\beta}=\left|\vec{c}_{2} \cdot \vec{x}_{\beta}-\vec{x}\right| \\ \vec{d}_{\delta}=\left|\vec{c}_{3} \cdot \vec{x}_{\delta}-\vec{x}\right|\end{array}\right.$       (35a)

$\left\{\begin{array}{l}\vec{x}_{1}=\vec{x}_{\alpha}-\vec{a}_{1} \cdot\left(\vec{d}_{\alpha}\right) \\ \vec{x}_{2}=\vec{x}_{\beta}-\vec{a}_{2} \cdot\left(\vec{d}_{\beta}\right) \\ \vec{x}_{3}=\vec{x}_{\delta}-\vec{a}_{3} \cdot\left(\vec{d}_{\delta}\right)\end{array}\right.$        (35b)

$x(n+1)=\frac{\vec{x}_{1}+\vec{x}_{2}+\vec{x}_{3}}{3}$         (35c)

2.3.3 Hybrid PSO-GWO algorithm

HPSOGWO is the hybridize Particle Swarm Optimization with the Grey Wolf Optimizer algorithm that uses the functionalities of both variants to enhance exploitation in PSO and GWO and to generated both variants’ strength [34].

In this algorithm, the position of the first three agents is updated in the search space using Eq. (34). The control of the grey wolf exploration and exploitation is conducted with the inertia constant. Using the following modified equations:

$\begin{aligned} \vec{d}_{\alpha} &=\left|\vec{c}_{1} \cdot \vec{x}_{\alpha}-\omega * \vec{x}\right| \\ \vec{d}_{\beta} &=\left|\vec{c}_{1} \cdot \vec{x}_{\beta}-\omega * \vec{x}\right| \\ \vec{d}_{\delta} &=\left|\vec{c}_{1} \cdot \vec{x}_{\delta}-\omega * \vec{x}\right| \end{aligned}$        (36)

The following equation is for combining the variants of PSO and GWO:  

$v_{i}^{k+1}=\omega *\left(\begin{array}{c}v_{i}^{k}+c_{1} r_{1}\left(x_{1}-x_{1}^{k}\right)+c_{2} r_{2}\left(x_{2}-x_{i}^{k}\right) \\ +c_{3} r_{3}\left(x_{3}-x_{i}^{k}\right)\end{array}\right)$

$x_{i}^{k+1}=x_{i}^{k}+v_{i}^{k+1}$         (37)

2.3.4 Whale optimization algorithm

This Algorithm consists of three stages [45-47]:

Bubble-net attacking method:

Using the bubble-net strategy, the humpback whales attack the prey. This is shown in the following two methods [45]:

• Shrinking encircling mechanism: it is applicable by reducing the amount of l in Eq. (33).
• Spiral updating position: the mathematical description of the helix-shape motion of humpback whales is:

$\vec{X}(n+1)=\vec{D}^{\prime} \cdot e^{b v} \cdot \cos (2 \pi v)+\vec{X}^{*}(n)$        (38)

$\vec{D}^{\prime}=\left|\vec{X}^{*}(n)-\overrightarrow{X(n)}\right|$        (39)

A probability of 50% to choose either a shrinking mechanism or a helical path down below the mathematical formula:

$\vec{X}(n+1)$

$=\left\{\begin{array}{l}\vec{X}^{*}(n)-\vec{A} \cdot \vec{D} \text { if } p<0.5 \\ \vec{D}^{\prime} \cdot e^{b v} \cdot \cos (2 \pi v)+\vec{X}^{*}(n) \text { if } p>0.5\end{array}\right.$        (40)

where, p is a random number in [0, 1].

Search for prey:

The humpback whales search randomly for prey. Exploitation phase fellow those rules [45]:

$\vec{D}=\left|\vec{C} \cdot \vec{X}_{\text {rand }}-\vec{X}\right|$        (41)

$\vec{X}(n+1)=\vec{X}_{\text {rand }}-\vec{A} \cdot \vec{D}$        (42)

2.3.5 Hybrid WOA-GWO algorithm

The HWGO algorithm aims to enhance the efficiency of the WOA algorithm by applying the leadership hierarchy of GWO. The next step is to apply the outcome to the WOA attacking strategy. The HWGO selected the three best candidate solutions: the 1^st, 2^nd, and 3^rd levels in the group are alpha (a), beta (b), and delta (d).

The mathematical model for updating the position of whales using the leadership hierarchy of GWO during optimization is given as:

• The updating of the position of whales using the hierarchy leadership of GWO is expressed by Eq. (35).
• The updating of the position along of a spiral-shape path of humpback whales is as follows:

$\left\{\begin{array}{l}\vec{D}_{\alpha}^{\prime}=\left|\vec{X}_{\alpha}(n)-\vec{X}\right| \\ \vec{D}_{\beta}^{\prime}=\left|\vec{X}_{\beta}(n)-\vec{X}\right| \\ \vec{D}_{\delta}^{\prime}=\left|\vec{X}_{\delta}(n)-\vec{X}\right|\end{array}\right.$           (43)

$\left\{\left\{\begin{array}{l}\vec{X}_{1}(n)=\vec{X}_{\alpha}(n)+\vec{D}_{\alpha} \cdot e^{b v} \cdot \cos (2 \pi v) \\ \vec{X}_{2}(n)=\vec{X}_{\beta}(n)+\vec{D}_{\beta} \cdot e^{b v} \cdot \cos (2 \pi v) \\ \vec{X}_{3}(n)=\vec{X}_{\delta}(n)+\vec{D}_{\delta} \cdot e^{b v} \cdot \cos (2 \pi v)\end{array}\right.\right.$       (44)

$\left\{\vec{X}(n+1)=\frac{\vec{X}_{1}+\vec{X}_{2}+\vec{X}_{3}}{3}\right.$       (45)

This paper presents an advanced optimization method to obtain the turning parameters of fractional order proportional integral derivative controller (FOPID controller). A comparison study was performed for the first time in the robotics field, where several tuning methods have been employed to set the details of two FOPID controllers, which are designed to make sure that the linear and angular speeds of the WMR follow the input reference. Those five meta-heuristics approached namely particle swarm optimizer (PSO), grey wolf optimizer (GWO), hybrid particle swarm grey wolf optimizer (HPSOGWO), whale optimizer (WOA) and hybrid grey wolf whale optimizer (HWGO) were found to be all statistically and graphically approved, by considering integral square error ISE as a cost function. A back-stepping kinematic controller was adopted to guarantee that tracking errors convergence to zero to achieve the trajectory tracking. A transient response in the time domain has been made to examine the performance of selected algorithms. Simulation results have successfully confirmed the over performance of HWGO in comparison to PSO, GWO, HPSOGWO, and WOA in terms of lowest error value obtained, convergence rapidity, with the minimum overshoot, received rise time, peak, peak time, and settling. Finally, a star trajectory demonstrated the robustness of the mentioned controller.