Optimizing FOPID controller utilizing combination of hedge algebra and wolf optimizer: A new solution to trajectory tracking for autonomous mobile robots

Abstract

Improving the path tracking performance of two-wheel differential mobile robots is very important, especially in problems that consider dynamic nonlinearities and motor torque constraints. This paper proposes a hybrid controller that combines backstepping control (BSC) with fractional order PID controller (FOPID). The parameters of BSC after being proven globally stable by Lyapunov will be determined by the hedge algebra method. The parameters of the FOPID controller are optimized using an improved metaheuristic algorithm, which is a combination of the wolf optimizer (GWO) algorithm with the slime molding algorithm. This is intended to improve trajectory tracking accuracy. In this work, the optimization process for FOPID is based on a cost function consisting of Integral absolute error and integral squared error, which helps to reduce the position and velocity errors. The controller performance is validated through MATLAB-Simulink with various trajectory scenarios and compared with conventional optimization methods such as standard PSO and GWO. Simulation results show superior trajectory tracking, error reduction and improved control performance.

Introduction

In robotic fields, differential wheeled mobile robots (DWMRs) are widely deployed due to their simple structures and practical performance with different types of complex navigation tasks. Different from omnidirectional robots that require sophisticated mechanical systems, DWMRs operate based on a minimalist configuration with two independently driven wheels and one passive caster wheel. ¹

Improving trajectory tracking capabilities is important for the DWMR in various applications. For example, autonomous cleaning robots may follow predefined paths precisely to ensure that the area in dynamic environments is covered. ² Similarly, in leader-follower control scenarios, the robot follows behind a human or another robot. And precise path tracking becomes necessary.^3,4 These practical applications usually involve nonlinear and dynamic trajectories. Hence, which necessitate advanced control strategies to ensure stable and accurate tracking performance.

The performance of DWMRs performance is constrained by actuator limitations, especially in terms of velocity and torque that result in nonlinear system dynamics.^5,6 As a result, many different control techniques have been proposed to deal with these challenges. In some related work, Li et al. ⁷ provided an adaptive neural network (NN)-based controller tailored for trajectory tracking in DWMRs with full-state constraints. Another study demonstrated the effectiveness of the backstepping control (BSC) technique in considering nonlinearities in DWMR dynamic control, improving the overall system performance.^8–11

Fuzzy set theory has proven effective in constructing flexible, highly adaptable control systems, especially in uncertain environments. However, some limitations remain, such as the heavy reliance on expert experience in establishing control laws, a lack of accuracy in handling large and complex systems, and limited scalability. Classical control techniques, despite their solid mathematical foundation, often require precise models, making them difficult to apply to real-world systems with significant disturbances or nonlinearities. Hedge algebra (HA) is a useful mathematical tool for modeling the linguistic values ​​of a linguistic variable. These values ​​are quantified by real numbers from 0 to 1. Therefore, HA has many advantages over traditional fuzzy set theory (FC) based controllers in terms of setup steps, control performance, computation time and optimization. ¹² In Mac et al., ¹³ the authors used HA as a position controller for the trajectory tracking problem of non-uniform mobile robots. Where HA is combined with improved PSO (APSO) resulting in an optimized version due to better tuned HA. However, the application of HA as a solution to optimize parameters for other controllers has not been mentioned by many studies.

Beyond kinematic control, velocity control plays a critical role in overall trajectory tracking performance. Several studies have focused on improving conventional PID controllers for DWMRs. For instance, Khan et al. ¹⁴ optimized a traditional PID controller to enhance velocity tracking, while Ahmad ¹⁵ applied the H∞ control technique to improve PID robustness against parameter variations caused by power fluctuations. Although there are many improvements proposed, conventional PID controllers are facing inherent limitations in handling nonlinear system dynamics and transient behavior, primarily due to their reliance on integer-order derivatives. To solve further problems, many researchers introduced controllers based on fractional-order calculus. Abed et al. ¹⁶ demonstrated that a neural network-based FOPID controller achieves better velocity tracking than its integer-order counterpart. In addition, Gheisarnejad and Khooban ¹⁷ provided the practical implementation of FOPID controllers in mobile robots.

A key challenge in both PID and FOPID controllers is the choice of optimal control parameters. Both manual random tuning methods often result in suboptimal system responses. To solve this problem, there are some methods that use bio-inspired swarm intelligence algorithms that have gained popularity. The methods require only the evaluation of a cost function rather than analytical gradients or integrals. ¹⁸ Several papers have focused on the problem of intelligent control and optimization for systems with nonlinear, uncertain, and time-varying characteristics. To reduce dependence on model inaccuracies, studies often use non-model-based or semi-model-based intelligent controllers, such as BELBIC, fuzzy logic, and FOPID.^19,20 These are then combined with modern metaheuristic optimization algorithms such as NSGA-III, GA, COA, MOQSOA, etc., to tune controller parameters, aiming to balance global exploration and local exploitation.^21–23

In the context, some work used Particle Swarm Optimization (PSO) to enhance path planning for DWMR under velocity constraints. ²⁴ This work proposed by Aner et al. ²⁵ presents a PSO-optimized PID controller that overcomes a PSO-optimized fuzzy logic controller (FLC) in trajectory tracking. Moreover, grey wolf optimizer (GWO) has been used to tune adaptive fractional-order fuzzy PID controllers, ²⁶ and hybrid schemes combining GWO with LQR and FOPID have also been reported to enhance DWMR performance. ²⁷

Although numerous studies have investigated either kinematic controllers (such as backstepping controller (BSC)) or dynamic controllers (such as PID/FOPID) separately, integrating a kinematic controller with PID or FOPID controllers significantly enlarges the parameter space to be optimized. This results in greater computational demands. As noted by the authors, the optimization of such hybrid control architecture remains an underexplored area in existing literature. Based on the existing limitations as presented, our research focuses on developing a more effective control strategy, with highlights including.

• A hybrid optimization strategy combining the algorithm with the Slime Mushroom (SMA) algorithm is developed for 10 parameters of the FOPID controller of the dynamic element.

• A cost function integrating integral absolute error (IAE) and integral squared error (ISE) is proposed. The weights of the cost function consider the effects of positional and angular errors.

• Hedge algebra is designed and used to determine the parameters of the BSC.

• A closed-loop control architecture integrating the BSC with a fractional order PID controller (FOPID) is proposed, where the FOPID parameters are optimized with velocity constraints.

• Simulations were performed with complex trajectories, including Lemniscate, Cloverleaf and Sinusoidal trajectories. Data showed that the new optimization techniques improved over other conventional methods such as PSO ²⁵ and GWO ²⁸ in terms of accuracy and cost function minimum target.

The rest of this paper is organized as follows. The background and research methods including Monitoring Control System, Backstepping based Trajectory Controller and FOPID based speed controller are provided in Section 2. In section 3, the combination of GWO and SMA is addressed. Section 4 provides the Improved GWO optimization algorithm. Simulation results are discussed in Section 5. Finally, Section 6 provides conclusions and future work.

Background and research methods

The DWMR is equipped with two independently driven rear wheels and a passive front wheel that provides balance, as depicted in Figure 1. In this figure, Q represents the midpoint of the wheel axis, G denotes the center of mass, and b is the distance between points G and B. Point Q = [xQ, yQ]T, lies on the XY plane. ${\theta }_l$ and ${\theta }_r$ are the angular velocities of the left and right wheels.

Figure 1.

TWMR mobile robot.

The standard kinematic model describing the relationship between the velocity (u, ω) and state (x, y, Ψ) of a differential robot is expressed as equation (1):

$\stackrel{·}{x}=u\hspace{0.17em}\cos \hspace{0.17em}\mathrm{\Psi }; \stackrel{·}{y}=u\hspace{0.17em}\sin \hspace{0.17em}\mathrm{\Psi }; \stackrel{·}{\mathrm{\Psi }}=\omega$

(1)

where x, y are robot coordinates in the global coordinate system, Ψ is orientation angle, u is linear velocity, and $\omega$ is angular velocity.

Additionally, a denotes the wheelbase of the robot, Ψ represents its posture (orientation), and γ corresponds to the angular velocity. The mathematical model of the DWMR is defined as follows ²⁹ :

$\left[\begin{array}{l}\begin{array}{l}\stackrel{·}{x}\\ \stackrel{·}{y}\end{array}\\ \stackrel{·}{\mathrm{\Psi }}\\ \stackrel{·}{u}\\ \stackrel{·}{\omega }\end{array}\right]=\left[\begin{array}{l}\begin{array}{l}u\hspace{0.17em}\cos \hspace{0.17em}\mathrm{\Psi }-a\omega \hspace{0.17em}\sin \hspace{0.17em}\mathrm{\Psi }\\ u\hspace{0.17em}\sin \hspace{0.17em}\mathrm{\Psi }+a\omega \hspace{0.17em}\cos \hspace{0.17em}\mathrm{\Psi }\end{array}\\ \omega \\ \frac{{\theta }_3}{{\theta }_1}{\omega }^2-\frac{{\theta }_4}{{\theta }_1}u\\ -\frac{{\theta }_5}{{\theta }_2}u\omega -\frac{{\theta }_6}{{\theta }_2}\omega \end{array}\right]+\left[\begin{array}{ll}\begin{array}{l}0\\ 0\\ 0\\ \frac{1}{{\theta }_1}\end{array}& \begin{array}{l}0\\ 0\\ 0\\ 0\end{array}\\ 0& \frac{1}{{\theta }_2}\end{array}\right]\left[\begin{array}{l}{u}_{ref}\\ u_{ref}\end{array}\right]$

(2)

where u and $\omega$ are the linear velocity and angular velocity of the robot, respectively. Similarly, and are the velocity inputs of the DWMR system. It should be noted that the dynamic model of the DWMR includes a passive rudder. In (2), are the internal parameters of the robot model that form a defined parameter vector. The entire robot can be described by the following kinematic and dynamic models respectively as follows:

$\stackrel{·}{\eta }=J\left(\eta \right)\varsigma$

(3)

$\stackrel{·}{\varsigma }=D\left(\varsigma \right)+n\left(\varsigma \right)\tau$

(4)

where $\eta ={\left[x,y,\mathrm{\Psi }\right]}^T,\varsigma ={\left[u,\omega \right]}^T$ and $\tau ={\left[u_{ref},{\omega }_{ref}\right]}^T$ denote the output and input of the dynamic model. The matrix J has dimensions of 3x2, D, and n are calculated as follows:

$\begin{array}{l}J\left(\eta \right)=\left[\begin{array}{ll}\begin{array}{l}\cos \hspace{0.17em}\mathrm{\Psi }\\ \sin \hspace{0.17em}\mathrm{\Psi }\end{array}& \begin{array}{l}-a\hspace{0.17em}\sin \hspace{0.17em}\mathrm{\Psi }\\ a\hspace{0.17em}\cos \hspace{0.17em}\mathrm{\Psi }\end{array}\\ 0& 1\end{array}\right]\\ D\left(\varsigma \right)=\left[\begin{array}{l}\frac{{\theta }_3}{{\theta }_1}{\omega }^2-\frac{{\theta }_4}{{\theta }_1}u\\ \frac{{\theta }_5}{{\theta }_2}u\omega -\frac{{\theta }_6}{{\theta }_2}\omega \end{array}\right],n\left(\varsigma \right)=\left[\begin{array}{ll}\frac{1}{{\theta }_1}& 0\\ 0& \frac{1}{{\theta }_2}\end{array}\right]\end{array}$

(5)

The above DWMR kinematic and the dynamic model only show the motion and speed characteristics of the mobile robot under ideal conditions. It does not consider the practical limitations due to mechanical constraints. The inputs of the linear velocity and angular velocity of the vehicle for the dynamic model are limited to the following ranges:

$\left|u\left(t\right)\right|\le u_{\max },\left|\omega \left(t\right)\right|\le {\omega }_{\max }$

(6)

where $u_{\max }$ and ${\omega }_{\max }$ are the maximum linear velocity and angular velocity of the mobile robot output, respectively.

Monitoring Control System

In this study, a closed-loop control structure is proposed to control the motion of the DWMR. The inner loop uses a FOPID controller for angular velocity and linear velocity. Given the current velocity of the robot $\varsigma ={\left[u,\omega \right]}^T$ , the inner loop controller aims to determine the control signal $\tau ={\left[u_{ref},{\omega }_{ref}\right]}^T$ so that ${\lim }_{t\to \infty }\hspace{0.17em}\varsigma \left(t\right)=c_{\varsigma }\left(t\right)$ .

The outer loop is the BSC controller for the kinematic model of the mobile robot. The BSC ensures stability of the position and attitude tracking of the mobile robot. When the reference trajectories for the robot are set to $r_d={\left[x_d,y_d,{\mathrm{\Psi }}_d\right]}^T$ and the current trajectory is $\eta ={\left[x,y,\mathrm{\Psi }\right]}^T$ the BSC searches for the appropriate output $c_{\varsigma }={\left[u_c,{\omega }_c\right]}^T$ so that ${\lim }_{t\to \infty }\hspace{0.17em}\eta \left(t\right)=r_d\left(t\right)$ .

The Saturation blocks shown in Figure 2 are used to constrain the velocity outputs from the FOPID controller. This reflects practical limitations, as motors and actuators are subject to torque constraints, which in turn restrict both angular and linear velocities.

Figure 2.

Overview of the cascade closed-loop control structure.

Backstepping based Trajectory Controller

Signals from the reference orbit, in addition to providing information about the desired position and attitude coordinates, also provide the desired linear and angular velocities:

$r_d=\left[x_d,y_d,{\mathrm{\Psi }}_d,{\stackrel{·}{x}}_d,{\stackrel{·}{y}}_d,u_d,{\omega }_d\right]$

(7)

BSC is a nonlinear control technique based on Lyapunov theory. It allows to construct control laws from Lyapunov functions. For the position and attitude of the DWMR, the position error transformation is expressed as follows:

$\left[\begin{array}{l}{e}_x\\ e_y\\ e_{\mathrm{\Psi }}\end{array}\right]=\left[\begin{array}{lll}\cos \left(\mathrm{\Psi }\right)& \sin \left(\mathrm{\Psi }\right)& 0\\ -\sin \left(\mathrm{\Psi }\right)& \cos \left(\mathrm{\Psi }\right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{l}{x}_d-x\\ y_d-y\\ {\mathrm{\Psi }}_d-\mathrm{\Psi }\end{array}\right]=\left[\begin{array}{l}{e}_x\\ e_y\\ e_{\mathrm{\Psi }}\end{array}\right]=\left[\begin{array}{l}\omega e_y-u+u_d\hspace{0.17em}\cos \left(e_{\mathrm{\Psi }}\right)\\ -u_d\hspace{0.17em}\sin \left(e_{\mathrm{\Psi }}\right)-\omega \left(e_x\right)\\ {\omega }_d-\omega \end{array}\right]$

(8)

To ensure the stability of the learning controller, the Lyapunov function is chosen as follows:

$L=\frac{e_x^2+e_y^2}{2}+\frac{1-\cos \left(e_{\mathrm{\Psi }}\right)}{k_2},$

where L is positive for all $k_2$ > 0. The time derivative of L is

$\begin{array}{l}\stackrel{·}{L}=e_x{\stackrel{·}{e}}_x+e_y{\stackrel{·}{e}}_y+\frac{\stackrel{·}{e_{\mathrm{\Psi }}}}{k_2}\sin \left(e_{\mathrm{\Psi }}\right)\\ =e_x\left(-u+u_d\hspace{0.17em}\cos \left(e_{\mathrm{\Psi }}\right)\right)+u_d\hspace{0.17em}\sin \left(e_{\mathrm{\Psi }}\right)\left(e_y+\frac{{\omega }_d-\omega }{k_2{u}_d}\right)\end{array}$

(9)

Determine two variables $u_1=-u+u_d\hspace{0.17em}\cos \left(e_{\mathrm{\Psi }}\right)$ and $u_2=e_y+\frac{{\omega }_d-\omega }{k_2{u}_d}$ , choose the error response such that and of system (8) are as follows:

$\begin{array}{l}{u}_1=-k_1{e}_x\\ u_2=-\frac{k_3}{k_2}\sin \left(e_{\mathrm{\Psi }}\right)\end{array}$

(10)

Substitute equation (10) into equation (9) and then simplify as follows:

$\stackrel{·}{L}=-k_1{e}_x^2-\frac{k_3}{k_2}{u}_d\hspace{0.17em}{\sin }^2\left(e_{\mathrm{\Psi }}\right)$

(11)

where $\stackrel{·}{L}\le 0$ and are continuous at $k_1,k_2,k_3$ are both >0 with $u_d$ > 0. And the positive Lyapunov function L, ${\left[e_x,e_y,e_{\mathrm{\Psi }}\right]}^T$ is simultaneously bounded. Hence the position errors ${\mathrm{e}}_{\mathrm{x}}$ and the angular errors ${\mathrm{e}}_{\mathrm{\Psi }}$ converge to 0. Hence the trajectory errors converge to 0 globally.

To determine the three parameters $k_1,k_2,k_3$ in the Backstepping controller of DWMR adaptively and nonlinearly according to the current state, the Hedge Algebra (HA) method is applied to convert the quantitative errors into semantic values, then infer back to obtain the corresponding control parameters.

The HA calculation process is carried out in six main steps including:

• Determine three input variables which are position and direction errors ( $e_x,e_y,e_{\varphi }$ ) and three corresponding output variables ( $k_1,k_2,k_3$ ).

• Normalize the input variables in the domain [-1,1] and determine the semantic domain for the output variables in [0,1].

• Divide the semantic domain into 9 linguistic value levels corresponding to the fuzzy levels “very low”, “low”, “slightly low”, “medium”, “slightly high”, “high”, “very high”, etc. Each level is quantified by the SQM value φ_i∈{0.0,”“ 0.125,”“ 0.250,”“ 0.375,”“ 0.500,”“ 0.625,”“ 0.750,”“ 0.875,”“ 1.000}.

• Construct functions belonging to triangle ${\mu }_i\left(x\right)$ such that the intersecting domains are continuous and ensure the normalized sum ${\sum }_i{\mu }_i=1$ .

• Establish quantitative semantic mapping for each control parameter through the relationship in formula (12):

$k_j={\displaystyle \sum }_{i=1}^9{\mu }_i\left(e_x,e_y,e_{\varphi }\right)·k_{j.level}\left(i\right),j=1,2,3$

(12)

where $k_{j.level}\left(i\right)$ is the base value corresponding to the ith language level.

• • Solve the approximate reasoning problem to determine the output quantitative semantic value, from which interpolate the quantitative semantic curve and determine the actual control value $k_j^{*}$ .

The selection of SQM values ​​ensures compatibility with linear scaling in the control parameter domain, simplifying implementation and improving stability during real-time computation. The semantic quantitative metrics (SQM) reference table for the nine levels is shown below Table 1:

Table 1.

Reference semantic quantitative metrics (SQM) for 9 levels.

Level	1	2	3	4	5	6	7	8	9
ϕᵢ	0.000	0.125	0.250	0.375	0.500	0.625	0.750	0.875	1.000
Language label	vvL	vL	lvL	L	W	H	lvH	vH	vvH

From three input errors ( $e_x,e_y,e_{\varphi }$ ), the HA rule base is established according to the semantic matrix (table below) – in which the row represents $e_y$ , the column represents $e_x$ , and each cell represents the semantic level of the corresponding control function (when considering the same influence of $e_{\varphi }$ ) Table 2:The SQM values ​​from this table are mapped linearly into the control domain of each parameter:

$k_1\in \left[0,10\right],k_2\in \left[0,30\right],k_3\in \left[0,50\right]$

(13)

In which, the parameter $k_3$ is adjusted according to the deflection angle ${\theta }_e$ by the enhancement factor:

$scale=1+0.5\left|{\theta }_e\right|,$

(14)

Table 2.

HA rule basis for 2 variables e_x and e_y.

e_y e_x	vvL	vL	lvL	L	W	H	lvH	vH	vvH
vvL	vvL	vvL	vL	lvL	L	H	lvH	vH	vvH
vL	vvL	vL	lvL	L	W	H	lvH	vH	vvH
lvL	vL	lvL	L	W	H	lvH	vH	vvH	vvH
L	lvL	L	W	H	lvH	vH	vH	vvH	vvH
W	L	W	H	lvH	vH	vvH	vvH	vvH	vvH
H	W	H	lvH	vH	vvH	vvH	vvH	vvH	vvH
lvH	H	lvH	vH	vvH	vvH	vvH	vvH	vvH	vvH
vH	lvH	vH	vvH	vvH	vvH	vvH	vvH	vvH	vvH
vvH	vH	vvH	vvH	vvH	vvH	vvH	vvH	vvH	vvH

FOPID based speed controller

Two FOPID controllers are designed to control the dynamics of the DWMR. The linear velocity and angular velocity are controlled separately.

The transfer function of the FOPID controller is (15):

$C\left(s\right)=K_p+K_i{\left(\frac{1}{s}\right)}^{\lambda }+K_d{s}^{\mu }$

(15)

The FOPID controller consists of five parameters $\left[K_p,K_i,K_d\right]$ and $\left[\lambda ,\mu \right]$ . If $\lambda =\mu =1$ , the FOPID controller is simplified into a classical PID controller. In the FOPID fractional controller, the fractional operators cannot be implemented directly in the simulation, so they need to be approximated by standard fractional transfer functions, in the form of products of first order (zero-pole) filters in the frequency domain. The FOPID controller and the outer-loop filter are related through fractional order operators. The outer-loop filter is used to implement fractional order operators in the frequency domain. It aims to approximate fractional order functions in the frequency range $\left({\omega }_b,{\omega }_h\right)$ and order N. Improved control of nonlinear systems with large inertia, or systems that are not accurately modeled.

The schematic diagram of the FOPID controller is shown in Figure 3.

Figure 3.

Internal configuration of the loop speed control.

The outer-loop recursive filter is given by :

$s^{\alpha }=K{\prod }_k^N=-N\frac{s+{\omega }_k^{\prime }}{s+{\mathrm{w}}_k}$

(16)

with $K={\omega }_h^{\alpha }$ ,

${\omega }_k={\omega }_b{\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{k+N+\frac{1}{2}\left(1+\alpha \right)}{2N+1}};{\omega }_k^{\prime }={\omega }_b{\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{k+N+\frac{1}{2}\left(1-\alpha \right)}{2N+1}}$

where N is the order of approximation representing the number of zero-pole pairs (2N + 1). In this problem N is chosen to be 3 and the filter range is $\left[{\omega }_b,{\omega }_h\right]$ = [0.01, 100].

The transfer functions of the linear velocity controller and the angular velocity controller based on FOPID for mobile robots can be expressed by the following equations:

$f\left(s\right)=\frac{u_x\left(s\right)}{e_u\left(s\right)}=K_{p1}+K_{i1}{\left(\frac{1}{s}\right)}^{\lambda 1}+K_{d1}{s}^{\mu 1}$

(17)

$g\left(s\right)=\frac{{\omega }_x\left(s\right)}{e_{\omega }\left(s\right)}=K_{p2}+K_{i2}{\left(\frac{1}{s}\right)}^{\lambda 2}+K_{d2}{s}^{\mu 2}\hspace{0.17em}$

(18)

The velocity saturation components are limited by the following piecewise function:

$u_y={\phi }_u\left(u_x\right)=\{\begin{array}{l}-u_{\max }\\ u_x\\ u_{\max }\end{array}\begin{array}{l}{u}_x<-u_{\max }\\ -u_{\max }\le u_x\le u_{\max }\\ u_x>u_{\max }\end{array}$

(19)

${\omega }_y={\phi }_{\omega }\left({\omega }_x\right)=\{\begin{array}{l}-{\omega }_{\max }\\ {\omega }_x\\ {\omega }_{\max }\end{array}\begin{array}{l}{\omega }_x<-{\omega }_{\max }\\ -{\omega }_{\max }\le {\omega }_x\le {\omega }_{\max }\\ {\omega }_x>{\omega }_{\max }\end{array}$

(20)

FOPID velocity controller under the velocity constraints as above. The DWMR closed-loop system is represented by state-space equations:

$\left[\begin{array}{l}{\stackrel{·}{e}}_u\\ {\stackrel{·}{e}}_{\omega }\end{array}\right]=\left[\begin{array}{l}{\stackrel{·}{u}}_c-\frac{{\theta }_3}{{\theta }_1}{\left({\omega }_c-e_{\omega }\right)}^2+\frac{{\theta }_4}{{\theta }_1}\left(u_c-e_u\right)-\frac{u_{ref}}{{\theta }_1}\\ {\stackrel{·}{\omega }}_c-\frac{{\theta }_5}{{\theta }_2}\left({\omega }_c-e_{\omega }\right)\left(u_c-e_u\right)+\frac{{\theta }_6}{{\theta }_2}\left({\omega }_c-e_{\omega }\right)-\frac{{\omega }_{ref}}{{\theta }_2}\end{array}\right]$

(21)

where $e_u=u_c-u,e_{\omega }={\omega }_c-\omega$ are the state variables of the velocity closed-loop system. The vector $c_{\varsigma }={\left[u_c,{\omega }_c\right]}^T$ represents the desired velocity vector in the inner loop. The maximum velocity limit used in this paper is $u_{\max }$ = 0.8 m/s and ${\omega }_{\max }$ = 1.7 rad/s. The main objective in designing the FOPID controller is to determine the parameter $\left[K_{p1},K_{i1},{\lambda }_1,K_{d1},{\mu }_1\right]$ sets for the linear velocity and $\left[K_{p2},K_{i2},{\lambda }_2,K_{d2},{\mu }_2\right]$ for the angular velocity.

Therefore, optimizing the parameters of the velocity controller is essential to ensure accurate trajectory tracking by the robot. To this end, we propose a hybrid optimization strategy that integrates the Grey Wolf Optimizer (GWO) with the Slime Mold Algorithm (SMA).

Optimal combination of GWO and SMA

GWO and SMA are two commonly used hyper-visual optimization algorithms in solving engineering problems. However, they operate on different search mechanisms. SMA is inspired by the behavior of slime molds in nature, while GWO is inspired by the social structure and hunting patterns of gray wolves. Although both are effective, GWO usually has good sensitivity and early convergence parameters. SMA strikes a balance between exploration and exploitation. In the context of trajectory tracking for DWMR, the optimization task becomes more complicated. This study combines GWO and SMA into a hybrid optimization method to improve the trajectory tracking performance of DWMR.

The Grey Wolf Optimizer (GWO) was first introduced by Mirjalili et al. ³⁰ and has since been widely adopted in various scientific and engineering optimization problems. In the algorithm, each grey wolf represents a possible solution in a D-dimensional searching space. There are top three fittest wolves, denoted as α (alpha), β (beta), and δ (delta), identified during each iteration and responsible for guiding the rest of the pack towards the optimal solution. In the following equations, the position of the i-th wolf ${\overrightarrow{X}}_i$ is updated based on the influence of these three leaders as:

${\overrightarrow{D}}_{\alpha }=\left|{\overrightarrow{C}}_1·{\overrightarrow{X}}_{\alpha }-{\overrightarrow{X}}_i\right|,{\overrightarrow{X}}_1={\overrightarrow{X}}_{\alpha }-{\overrightarrow{A}}_1·{\overrightarrow{D}}_{\alpha }$

(22)

${\overrightarrow{D}}_{\beta }=\left|{\overrightarrow{C}}_2·{\overrightarrow{X}}_{\beta }-{\overrightarrow{X}}_i\right|,{\overrightarrow{X}}_2={\overrightarrow{X}}_{\beta }-{\overrightarrow{A}}_2·{\overrightarrow{D}}_{\beta }$

(23)

${\overrightarrow{D}}_{\delta }=\left|{\overrightarrow{C}}_3·{\overrightarrow{X}}_{\delta }-{\overrightarrow{X}}_i\right|,{\overrightarrow{X}}_3={\overrightarrow{X}}_{\delta }-{\overrightarrow{A}}_3·{\overrightarrow{D}}_{\delta }$

(24)

${\overrightarrow{X}}_i\left(t+1\right)=\frac{{\overrightarrow{X}}_1+{\overrightarrow{X}}_2+{\overrightarrow{X}}_3}{3}$

(25)

where t represents the current iteration; ${\overrightarrow{D}}_{\alpha ,\beta ,\delta }$ and ${\overrightarrow{X}}_{1,2,3}$ are the intermediate vectors of the wolf. $\overrightarrow{A}$ is the wolf’s exploitative signature while following the prey determined by a random value in the interval [-a, a]; This value of “a” reduces linearly from 2 to 0 during the iteration; |A| > 1 presents the wolf’s search process while |A| < 1 denotes the exploitation process; $\overrightarrow{C}$ dynamically reconciles the distance between the prey and the wolf according to the weights. Typically, $\overrightarrow{A}=2\overrightarrow{a}\overrightarrow{r_1}\overrightarrow{a}$ and $\overrightarrow{C}=2\overrightarrow{r_2}$ with $r_1$ , $r_2$ are two uniform distributed random number within [0,1].

The SMA provides the foraging behavior of slime molds, balancing between exploration and exploitation. ³¹ The updating position mechanism of the SMA is managed by adaptive weights. The adaptive weights are calculated based on the weight vector W, which affects the search direction and the intensity:

$\overrightarrow{W}(\text{SmellIndex}\left(i\right)=\{\begin{array}{l}1+r\hspace{0.17em}\log \left(\frac{bF-S\left(i\right)}{bF-wF}+1\right),\mathrm{c}\text{onditions}\\ 1-r\hspace{0.17em}\log \left(\frac{bF-S\left(i\right)}{bF-wF}+1\right),\text{others}\end{array}$

(26)

$\text{SmellIndex}=\text{sort}\left(\mathrm{S}\right),$

(27)

where S(i) represents the fitness of the i-th agent; wF and bF are the worst and best fitness obtained in the current iteration, respectively; the SmellIndex is a sorted array that contains the agents’ fitness values in ascending or descending order.

Initially, the agents adjust their positions based on the adaptive weights to explore the search space effectively and move toward promising regions. They then proceed to search for optimal paths leading to the food source. Once potential food sources are located, agents surround the best one and update their positions accordingly. This process ensures a balance between exploration and exploitation, improving convergence towards the optimal solution. The position update equation is as follows:

$\overrightarrow{X}\left(t+1\right)=\{\begin{array}{l}\in ·\left(UB-LB\right)+LB,\\ {\overrightarrow{X}}_b\left(t\right)+\overrightarrow{vb}·(\overrightarrow{W}·{\overrightarrow{X}}_A\left(t\right)-{\overrightarrow{X}}_B\left(t\right)\\ \overrightarrow{vc}·\overrightarrow{X}\left(t\right),\end{array}\begin{array}{l}\in <z\\ \hspace{1em}r<p\\ \hspace{0.5em}r\ge p\end{array}$

(28)

where $X_b$ denotes the location of highest concentration, ${\overrightarrow{X}}_A$ and ${\overrightarrow{X}}_B$ denote two randomly chosen agents, LB and UB are the search range boundaries, z and p are the selection conditions that determine the transitions of X(t), ε and random numbers in the interval [0, 1]. The value of $\overrightarrow{vb}$ randomly changes between [-a, a] and decreases to zero as the number of iterations increases, the value of $\overrightarrow{vc}$ has the same property as $\overrightarrow{vb}$ , but the boundary is limited by [-1, 1]. The parameters p and a are defined as follows:

$p=\tanh \left|S\left(i\right)-DF\right|$

(29)

$a=\mathrm{arctanh}\left(-\left(\frac{t}{T}\right)+1\right)$

(30)

where DF is the best fit computed over all iterations, and T denotes the maximum number of iterations.

Improved GWO optimization algorithm

To effectively solve the optimization problem presented in Part 3. The operating mechanism of the hybrid algorithm as well as the objective function will be presented in this section. To improve the optimization process, the search space $\mathrm{\Omega }$ (parameters of the two FOPID controllers) is divided into ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ . GWO tries to determine the three best solutions in $\alpha ,\beta$ and in $\delta$ . SMA tries to achieve the best fit $X_b$ in ${\mathrm{\Omega }}_2$ successive iterations. The agent $X_i^{*}$ in ${\mathrm{\Omega }}_2$ can be expressed by the following equation:

$X_i^{*}=UB+LB-X_i,X_i\in {\mathrm{\Omega }}_1,X_i^{*}\in \mathrm{\Omega }$

(31)

The pseudocode for the hybrid optimization technique is shown in Algorithm 1:

The SMA is carried out in ${\mathrm{\Omega }}_2$ at the same time. $X_b$ is the best fit that will be recorded at each iteration In order to evaluate the performance of GWO-SMA for trajectory tracking of DWMR, a cost function F is suggested by combining both IAE and ISE for position and angle. Denoted as $IA{E}_{XY},IA{E}_{\mathrm{\Psi }},IS{E}_{XY}$ and $IS{E}_{\mathrm{\Psi }}$ respectively. By combining both the IAE and the ISE into the cost function, we earn some of their complementary benefits. The IAE supports reducing persistent steady-state errors over time. And the ISE supports minimizing larger deviations in system response.^22,32 This dual-metric method allows the DWMR to earn improved accuracy in both position and orientation tracking, and to promote a balanced performance between precision and system stability. The definitions are provided mathematically as follows:

$\begin{array}{l}IA{E}_{XY}={\int }_0^{tf}\left|e_x\left(t\right)\right|+\left|e_y\left(t\right)\right|dt\\ IA{E}_{\mathrm{\Psi }}={\int }_0^{tf}\left|e_{\mathrm{\Psi }}\left(t\right)\right| dt\\ IS{E}_{XY}={\int }_0^{tf}\left(e_x^2\left(t\right)\right)+\left(e_x^2\left(t\right)\right)dt\\ IS{E}_{\mathrm{\Psi }}={\int }_0^{tf}\left(e_{\mathrm{\Psi }}^2\left(t\right)\right) dt\end{array}$

(32)

where tf denotes the total simulation time. Hence, the cost function is as follows:

$F=w_{1}IA{E}_{XY}+w_{2}IA{E}_{\mathrm{\Psi }}+w_{3}IS{E}_{XY}+w_{4}IS{E}_{\mathrm{\Psi }}$

(33)

where, $w_1$ = $w_3$ = 0,4, $w_2$ = $w_4$ = 0,1 are the weights used in the optimization process. The weights are set to support the robot’s position on the XY plane over its angle.

In order to validate the effectiveness of the proposed method, the hybrid method for BSC-FOPID optimization, there is a comparison with PSO, ²⁵ and GWO. ²⁸ The maximum number of iterations and the population size are set to be uniform to 25 and 60, respectively. The parameters of the DWMR system are:

$\theta =\left[0.2604,0.2509,-0.000499,0.9965,0.00263,1.0768\right],$

where a = 0.2 m. ²⁹

The actual position of the DWMR in the global coordinate system is $\eta$ (0) = [0.2, 0, 0]. A lemniscate trajectory is chosen to optimize the initial control parameters, which can be mathematically described by (34):

$x_d=\sin \left(2{\omega }_{d}t\right),y_d=\sin \left({\omega }_{d}t\right),$

(34)

where ${\omega }_d$ = 0.25 rad/s and t ∈ [0, 40].

In order to demonstrate the generalization ability of the optimized controller and ensure objective evaluation, the Cloverleaf and the Sin trajectory is added to validate the controller’s performance. This trajectory is mathematically described as follows ³³ :

${\stackrel{·}{x}}_d=0.052\hspace{0.17em}\cos \left(0.13t\right)+0.056\hspace{0.17em}\cos \left(0.39t\right)$

(35)

${\stackrel{·}{y}}_d=-0.052\hspace{0.17em}\sin \left(0.13t\right)+0.156\hspace{0.17em}\sin \left(0.39t\right),t\in \left[0,50\right]$

(36)

${\stackrel{·}{x}}_d=0.3t,$

(37)

${\stackrel{·}{y}}_d=0.05\hspace{0.17em}\sin \left(0.5t\right),t\in \left[0,50\right]$

(38)

Results and discussion

Before proceeding with optimization, it is necessary to set initial values ​​for 10 parameters of 02 FOPID sets. Details of the initial setup values before optimization and the 10-parameter limit of the two FOPIDs for controlling angular velocity and linear velocity are given in Table 3.

Table 3.

Setting initial values and limits of FOPID parameters.

	$k_{p1}$	$k_{i1}$	$lambd{a}_1$	$k{d}_1$	$m{u}_1$	$k_{p2}$	$k_{i2}$	$lambd{a}_2$	$k{d}_2$	$m{u}_2$
Initial value	1	1	0.7	1	0.7	1	1	0.7	1	0.7
Optimization limits	[0,30]	[0,30]	[0.01,1]	[0,30]	[0.01,1]	[0,30]	[0,30]	[0.01,1]	[0,30]	[0.01,1]

The controller simulation diagram is built on Simulink as shown in Figures 4 and 5.

Figure 4.

Simulation diagram on Simulink.

Figure 5.

Building cost function in Simulink.

The surface has a “symmetrical fold” shape, showing a clear change of $k_1$ when $e_x$ changes sign from negative to positive. When the robot is behind the desired trajectory (e_x < 0), $k_1$ is small (blue area), the control velocity is reduced to avoid overshooting the trajectory. When the robot is ahead or deviates positively (e x > 0), $k_1$ increases rapidly (red area), helping to increase the correction velocity and bring the robot back to the standard trajectory. The change in $e_y$ is smaller, this is because $k_1$ mainly adjusts the linear velocity component according to the longitudinal error, while $e_y$ mainly affects the angular control ( $k_2,k_3$ ). The range of values $k_1$ ∈[0,10] shows that the HA set has interpolated the coefficient within the desired limit, ensuring smooth and stable control, avoiding “mutation” or oscillation when the deviation is large. The surface changes continuously without break, demonstrating that the hedge algebra has interpolated semantically properly while maintaining smoothness in the control domain (Figure 6).

Figure 6.

Semantic aspect of k1 according to e_x and e_y when determined by HA.

The surface has a symmetrical “V” shaped fold, similar to $k_1$ but with a larger slope because the value of $k_2$ varies between 0 and 30. The colors range from black to red, yellow and the largest is white, representing a sharp increase in $k_2$ in both the variables $e_x$ and $e_y$ . When $e_x$ and $e_y$ are both positive or both negative, $k_2$ reaches its maximum value (representing the strongest adjustment strength in situations where DWMR deviates from the orbit in both directions). As the robot deviates laterally from the trajectory (|e_y| is large), the coefficient $k_2$ increases to increase the angular velocity ω and the DWMR rotates faster towards the standard direction. When the longitudinal error (e_x) is also large, $k_2$ is further amplified by HA showing the coordination between the two errors, helping to simultaneously adjust both position and direction (Figure 7).

Figure 7.

Semantic aspect of k2 according to e_x and e_y when determined by HA.

The 3D surface is smoothly curved and symmetrical, showing a soft nonlinear relationship between angular velocity and deviations. When the deviation $e_y$ >0 (DWMR deviates to the left from the desired trajectory), the angular velocity ω increases positively, corresponding to the orange-red region in the upper part of Figure 8). The robot will turn right to correct the direction. Conversely, when $e_y$ <0 (right bias), the angular velocity is negative (blue region) for the robot to rotate backwards, which demonstrates the controller’s symmetrical and stable response. The influence of $e_x$ (vertical deviation) is relatively small, as shown by the parallel contour lines near the $y_e$ axis. This is consistent with the theoretical model, since the directional control depends mainly on the horizontal deviation $e_y$ and the yaw angle ${\theta }_e$ . The value domain ω ∈[-20,20] is within the physical allowable range of the robot, showing that the HA set has adjusted $k_2,k_3$ reasonably, avoiding over-adjustment or strong oscillation.

Figure 8.

Backstepping output angular velocity when determining parameters using HA.

The smooth, continuous surface form demonstrates the smooth interpolation capability of the HA system during the determination of the $k_1$ parameter (Figure 9). There were no sudden breaks or fluctuations, which demonstrated that the HA system was stable and had good language-to-quantity conversion capabilities. When the deviation $e_x{e}_x$ >0 (the robot deviates ahead of the desired point), the linear velocity ν increases linearly (light yellow region) corresponding to large positive values. This makes the robot move forward faster to reduce the deviation along the x-axis, which is consistent with the feedback rule of Backstepping, and vice versa when $e_x{e}_x$ <0 the long velocity is negative, showing the ability to decelerate or back slightly to adjust the position—represented by the blue area on the opposite side.

Figure 9.

Backstepping output velocity when parameterizing with HA.

From the trajectory responses of DWMR and reference trajectory, the position and angle errors are calculated and used to build the cost function for the optimization process. Details are shown in Figures 10,11, and 12. In all three cases, our solution for robot trajectory closely followed the reference trajectory better than GWO and PSO, especially in areas with large curvature and rapid turns. Compared to the other two methods, GWO-SMA reduced positional error, limited phase lag and amplitude oscillation, and improved the system’s convergence speed.

Figure 10.

Lemniscate orbit response of optimization techniques.

Figure 11.

Cloverleaf orbit response of optimization techniques.

Figure 12.

Sin orbit response of optimization techniques.

The cost function, in addition to combining the two indexes ISE and IAE to evaluate the optimization process as presented in Part II, the author also builds the ITAE index to be able to expand the evaluation and build new cost functions for subsequent studies.

The trajectory tracking simulation results show that the optimal solution with improved GWO (red line) tracks the set trajectory better than GWO (blue line) and PSO (yellow line). Especially in the Sin trajectory, the obtained trajectory has high agreement with the set trajectory.

The x- and y-axis deviation results are shown and compared as in Figure 13(a),(b). It can be seen that the GWO-SMA method shows superior performance compared to PSO and GWO in all aspects. Specifically, the positional error along the x and y axes of GWO-SMA is mainly in the range of ±0.03–0.05 m, while GWO and PSO achieve approximately ±0.10–0.12 m and ±0.18–0.22 m respectively, representing a reduction in error of up to 70–80% compared to PSO. Although there is a fluctuation of the deviation, the convergence speed to zero of the y-axis deviation is under 10 s.

Figure 13.

Error results along x-axis, y-axis and theta angle deviation comparison. a) Error results along x-axis b) Error results along y-axis c) Theta angle deviation comparison.

While the pure GWO optimization and PSO have much slower convergence speeds.

Comparing the Theta angle error (Figure 13(c)) shows that our method, although having a larger initial error than GWO (approximately 0.75 rad), has the fastest convergence speed (7.02 s). PSO has the largest error and convergence speed among the three methods.

The angular orientation error is also significantly reduced, from approximately ±1.0 rad (PSO) and ±0.5 rad (GWO) to ±0.15–0.2 rad with GWO-SMA. Regarding velocity error, the maximum amplitudes of linear and angular velocities decreased from 1.1 to 1.2 m/s and ±0.8–1.0 rad/s (PSO), respectively, to approximately 0.1–0.15 m/s and ±0.1–0.15 rad/s when using GWO-SMA (Figure 14). The inertial velocity and angular velocity errors show that the improved GWO gives better results than GWO and PSO; especially the angular velocity error shows superiority. Although the inertial velocity error does not show a clear improvement, the oscillation amplitude and the state transition in the trajectory sections with sharp turns show smoother results than the pure PSO and GWO.

Figure 14.

Comparison of linear velocity and angular velocity between optimization methods. a) Linear velocity error b) Angular velocity error.

With three relatively complex trajectories selected, with many sharp curves, the controller is required to operate stably and efficiently. In both trajectories, the improved GWO optimization shows better trajectory tracking ability than the basic PSO and GWO optimization. The parameters of the selected optimization techniques are detailed as shown in Table 4.

Table 4.

Setting parameters for each optimization algorithm.

Optimization algorithm	Parameters
PSO	$wϵ,\left[0.4,2\right],c_1=c_2=2$
GWO	$aϵ,\left[0,2\right]$
GWO-SMA	$sϵ,\left[0,1\right],PSRs=0.34$

After the first 4 iterations, the objective function has changed significantly, and the final Best Fitness is 23.1754 (Figure 15). The parameters of the 2 FOPID sets controlling the linear velocity and the obtained angular velocity are: [kp1, ki1, kd1, λ1, μ1, kp2, ki2, kd2, λ2, μ2] = [15.8028, 29.9719, 1.2855, 0.9937, 0.4726, 16.2090, 7.7351, 9.0467, 0.3469, 0.9689].

Figure 15.

Change of objective function over 25 iterations with improved GWO - SMA optimization technique.

GWO optimization shows that the objective function starts to decrease sharply from iteration 15 onwards. However, the final Best Cost value obtained after 25 iterations is still 73.391023 (Figure 16). The parameters of the final FOPID obtained are: [kp ₁ , ki ₁ , kd ₁ , λ ₁ , μ ₁ , kp ₂ , ki ₂ , kd ₂ , λ ₂ , μ ₂ ] = [7.2634, 4.5317, 19.7268, 0.3842, 0.9968, 0.0491, 0.0752, 0.7005, 0.0621, 0.9125].

Figure 16.

Change of objective function over 25 iterations with GWO optimization technique.

The PSO technique shows that the objective function has a much smaller change than GWO and improved GWO. This may be because the number of iterations 25 is small compared to the PSO technique while the optimization space and the number of parameters of the optimization process are large (Figure 17). The final Global Best Cost of PSO is 41.0424 and the FOPID parameters obtained after optimization are. The PSRs (Probability of Random Strategy Switching) parameter we chose is 0.34.

Figure 17.

Change of objective function over 25 iterations with PSO optimization technique.

[kp ₁ , ki ₁ , kd ₁ , λ ₁ , μ ₁ , kp ₂ , ki ₂ , kd ₂ , λ ₂ , μ ₂ ] = [23.4820, 18.2709, 14.3248 0.6651, 0.7969, 14.9474, 13.5385, 8.0807, 0.3314, 0.7090].

A summary of the comparative indices of the optimization solutions is presented in detail in Table 5.

Table 5.

Summary comparison of performance indicators of optimization methods.

Performance index	PSO	GWO	GWO - SMA
Maximum position error (x, y) (m)	±0.18-0.22	±0.10-0.12	±0.03-0.05
Maximum orientation error (rad)	≈ ± 1.0	≈ ± 0.5	±0.15-0.20
Maximum linear velocity error (m/s)	1.1-1.2	≈0.3-0.4	0.10-0.15
Maximum angular velocity error (rad/s)	±0.8-1.0	±0.4-0.5	±0.10-0.15
Final objective function value	41.0424	≈73-75	23.1754

Conclusions and future work

This study proposes an innovative solution for trajectory tracking of DWMRs by integrating a Backstepping Controller (BSC) with a Fractional Order PID (FOPID) controller within a closed-loop control system. We have applied hedge algebra to determine the parameters for the Backstepping controller. The improvement of an enhanced Grey Wolf Optimization (GWO) strategy, in integration with the Slime Mold Algorithm (SMA), proven that it leads to improved trajectory tracking performance. By integrating both ISE and IAE into the cost function with carefully selected weights, the proposed method achieves better balance between exploration and exploitation while compared to conventional optimization techniques such as PSO and standard GWO. The results are verified with the efficiency of the proposed BSC-FOPID control approaches. Despite its contributions, our BSC–FOPID method, combined with HA-GWO–SMA optimization techniques, still has some limitations. The selection of weights in the ISE–IAE objective function is still based on experience, and we have only compared and evaluated it with other techniques based on the objective function we developed.

For future work, there are potential research directions that can expand upon these findings. The improved GWO optimization framework could be deployed to other types of mobile robots or robotic systems to evaluate the adaptability and the effectiveness across different applications. In addition, further improvement of the weighting parameters or integration of alternative metrics could reduce late-stage errors in the optimization process and enhance system performance in complex environments. Additionally, the improved GWO against other nature-inspired algorithms such as the Whale Optimization Algorithm (WOA) and Ant Colony Optimization (ACO) could be further exploited. Finally, we could integrate machine learning techniques to detect and mitigate wheel slip in real-time to increase the robustness and adaptability of DWMRs in unpredictable and dynamic settings.