Lyapunov-based control for enhanced trajectory tracking in robotic arms using multi-objective particle swarm optimization

Fundamental equations

In the classical PSO algorithm, each particle i has a position vector $x_i$ and a velocity $v_i$. At iteration $t$, these are updated as follows (Zeedan, Attiya & El-Fishawy, 2022; Li, Yan & Qu, 2020; Yang et al., 2023).

(1) $$V_i^{t+1}=w{v}_i^t+c_1{r}_1\left(pbes{t}_i-x_i^{\left(t\right)}\right)+c_2{r}_2\left(gbest-x_i^{\left(t\right)}\right),$$where $w$ is the inertial weight, $c_1$ is the cognitive coefficient, $c_2$ is the social coefficient, $r_1$ and $r_2$ are random scalars in [0, 1], $pbes{t}_i$ is the best position found by the particle $i$, and $gbest$ is the global best among the swarm (or a neighborhood best in local PSO variants).

Position update step

$$x_i^{\left(t+1\right)}=x_i^{\left(t\right)}+v_i^{\left(t+1\right)}.$$Proper selection of $w$, $c_1$, and $c_2$ strongly affects global exploration vs local exploration (Lobato-Larios, Starostenko & Alarcon-Aquino, 2024; Li, Yan & Qu, 2020; Yang et al., 2023).

PSO for multi-objective problems

In many real-world control applications, it is common to have multiple objectives that often conflict with each other. For instance, in robotics, minimizing tracking error may conflict with minimizing energy consumption: pushing the system to follow a path with very high precision typically requires additional torque, thus consuming more energy. When there are many aims to be optimized simultaneously, the Pareto front idea (or a set of non-dominating solutions) is gaining significance (Wang & Zhou, 2024; Sun, Deng & Li, 2020; Li, Zhang & Hu, 2023; Ma et al., 2025).

Non dominated solutions

A solution $x^{\left(1\right)}$ is said to dominate another solution $x^{\left(2\right)}$ if and only if:

• 1.

  $x^{\left(1\right)}$ is not worse than $x^{\left(2\right)}$ in all objectives, and

• 2.

  $x^{\left(2\right)}$ is strictly better in at least one objective.

In contrast, a solution is non-dominated if no other solution in the population (or archive) dominates it. The collection of all such non-dominated points constitutes the Pareto front. Each point on the front represents a distinct trade-off between objectives. Key enhancements in MOPSO include diversity preservation (via crowding distance) to maintain a uniform approximation of the Pareto front. Adaptive coefficients (w, $c_1,c_2$) to favor exploration at early stages and exploitation at later ones (Barba-González et al., 2020; Xu, Zhang & Chen, 2020). Mutation or perturbation operators to avoid stagnation in suboptimal regions (Wang et al., 2024c). Rather than finding a single “best” set of improvements, it is possible to choose between a range of equivalent bests. That is, a marginally increased tracking error can be tolerated when it reduces power consumption (or vice versa). By viewing the entire Pareto front, it is possible to determine which performance criteria hold more weight (precision vs energy, accuracy vs cost, etc.).

MOPSO algorithm

MOPSO is an extension of the original PSO algorithm developed to deal with several objectives and produce a Pareto front instead of a single global optimum. Classical PSO aims for a single optimal (global best) solution, but MOPSO retains and updates an external archive (or repository) of non-dominated solutions. The archive then leads to the swarm during movement in the search space. At each iteration, the particles are evaluated considering various objectives (e.g., tracking error and energy consumption). The algorithm adds all newly found non-dominated particles to the archive and eliminates existing solutions found to be dominated by the newcomers. The MOPSO tends to choose one or several of the “leaders” from the archive to influence the velocity update of the particles. Instead of a single $gbest$ leader, there might be multiple leaders, each of which is equally valid for Pareto dominance. Each particle still maintains its personal best ( $pbes{t}_i$), but when searching for the global best, MOPSO randomly chooses one solution (a “leader”) from the external archive. This leader selection can use strategies, such as roulette wheel based on crowding distance, or simply random selection among all non-dominated solutions. Over iterations, individual experience ( $pbes{t}_i$) and diversity of leaders (the non-dominated solutions in the archive) are utilized by the swarm.

Ideally, particles will settle down to a well-diversified Pareto front, offering several trade-off solutions. Throughout the optimization process, MOPSO may adapt the inertial weight w and the cognitive/social factors $c_1$, $c_2$. A high inertial weight induces overall search in early iterations, and a low inertia weight aids in refinement locally (exploitation) in later iterations. This mechanism can become a critical feature in evading early convergence to a local Pareto front and in having an ideal balance between searching new regions and enhancing the current best collection of solutions (Barba-González et al., 2020; Xu, Zhang & Chen, 2020; Ma et al., 2025; Wang & Zhou, 2024; Yang et al., 2023), MOPSO naturally handles distinct objectives in parallel. Maintaining an external archive gives a global vision of search progress, pointing out which solutions are currently non dominated and ensuring that valuable trade-offs are preserved. Once the algorithm terminates, the final Pareto front can be inspected and the solution that best fits the practical constraints can be selected (e.g., a midpoint solution balancing moderate energy usage with acceptable accuracy).

Kinematics and dynamics of a 2-DOF robotic manipulator

Robot manipulation is commonly in the form of nonlinear, coupled dynamics, presenting difficulties for typical control schemes (Alattas et al., 2021; Yan et al., 2023; Craig, 2022). A 2-DOF planar manipulator is composed of two rigid links that are connected through revolute joints, as depicted in Fig. 1. The following variables should be defined:

The angles of the joints, ${\theta }_1$ and ${\theta }_2$, in terms of angles relative to a reference axis; the length of the links as the quantities $L_1$ and $L_2$; the mass of the first and second links as the quantities $m_1$ and $m_2$, respectively; the vector $\overrightarrow{g}$ is the gravity, acting upon the masses in the negative Z; the location of the centers of mass for the first and second links as the quantities $c_1$ and $c_2$, respectively; the control torques acting on the first and second joints, so that ${\tau }_1$ and ${\tau }_2$ respectively, are the quantities that drive the motion of the first and second links, so that the motion of the end of the second link follows the desired motion in the plane, and $F_e$ is the external force acting upon the end-effector, decomposed along the $F_{ex}$ and $F_{ez}$.

Kinematic and dynamic modeling of the 2-DOF robot manipulator

Gain optimization with MOPSO

To overcome heuristic or trial-and-error tuning, MOPSO can be employed as follows (Wang et al., 2024c):

• 1.

  Define objectives:

  • $J_T\left(K_p,K_d,K_i\right)$: tracking error.

  • $J_C\left(K_p,K_d,K_i\right)$: control error.

• 2.

  Set feasible ranges for the gains $\left(K_p,K_{d,}{K}_i\right)$.

• 3.

  Optimization loop:

• Each particle represents a candidate $\left(K_p,K_d,K_i\right)$.

The robot’s dynamic simulation computes $(J_T,J_C)$ for that candidate.

Non-dominated solutions are stored in Pareto repository, guiding the swarm’s evolution.

By combining Lyapunov stability with MOPSO’s capability to search multi-objective spaces, one obtains robust and efficient gain values (Barba-González et al., 2020).

Trajectory generation

We consider various end effector paths. Besides simple shapes, we include more intricate patterns to challenge the controller across wide operating ranges. Given (x_d(t), z_d(t)), we obtain ${\theta }_d\left(t\right),{\dot{\theta }}_d\left(t\right),{\ddot{\theta }}_d\left(t\right)$ through forward differentiation and inverse kinematics (Hernandez-Barragan et al., 2022).

Multi-objective fitness formulation

Tracking Error $(E_T)$:

(24) $$J_T={\displaystyle \frac{1}{T}\underset{0}{\overset{T}{\int }}\Vert e_{\theta }\left(t\right){\Vert }^{2}dt,}$$quantifying joint space deviation (Dorf & Bishop, 2022).

Control Effort $(E_C)$:

(25) $$J_C={\displaystyle \frac{1}{T}\underset{0}{\overset{T}{\int }}\Vert \tau \left(t\right)\Vert dt,}$$where T represents the total time interval for which the integral is calculated representing torque or energy usage (Dorf & Bishop, 2022).

$J_T$ represents the mean squared joint-space deviation and directly quantifies tracking accuracy; $J_C$ measures the time-averaged absolute torque, capturing the total mechanical impulse delivered by the actuators.

Hence, the MOPSO-based optimization problem is:

Let

$$K={\left(K_p,K_d,K_i\right)}^Tϵ{\mathbb{R}}^3,$$denote the vectors of proportional, derivate, and integral gains.

The multi-objective optimization problem is stated as

$$\underset{Kϵ\mathrm{\Omega }}{min}\left(J_T\left(K\right),J_{C\left(K\right)}\right),$$

(26) $$\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\Omega }=\left\{\begin{array}{c}Kϵ{\mathbb{R}}^3:K_p^{min}\le K_p\le K_p^{max},\\ K_d^{min}\le K_d\le K_d^{max},\\ K_i^{min}\le K_i\le K_i^{max}\end{array}\right\},$$The feasible set a MOPSO algorithm is employed to construct a Pareto front of non-dominated solutions, providing a spectrum of gain combinations that balance tracking accuracy against actuator effort.

Optimization objectives

Minimize Total Error ( $E_T$): quantified using the mean squared error (MSE) between the desired and actual trajectories.

Minimize Control Effort ( $E_C$): Measured by the integrated absolute torque applied over the simulation period.

A PID based on a Lyapunov function is employed in this work (see Fig. S1 in Supplemental Files). The input parameters are set first. The system designs a desired path (gear-shape and rose curve), followed by a robotic arm. Inverse kinematics is utilized by the system to identify the needed angles and positions to track a desired path. For the flowing and smooth motion of the system, resulting paths in angles are differentiated and smoothed, computing the accelerations along with the velocities of the angles. The system is then controlled using a Lyapunov PID based, which is created on the error signals, their derivatives, and their integrals that have been computed and utilized in the dynamic adjustment of the control. The system is then tracked using the computation of parameters such as overshoot. Finally, the system plots the computed results.

One of the proposed strategies in the process of research is the optimization of parameters and Lyapunov based PID control using MOPSO (see Fig. 2). The process begins with the loading of parameters and trajectory. Inverse kinematic computations are conducted in the determination of joint angles to go along a given path. The process begins with the initialization of an initial set of candidate solutions using random generation based on given limits. Each of these is given an initial cost using the given performance metrics. The iterative MOPSO optimization is conducted to determine the best controlling gains. In each of these iterations, it checks for the satisfaction of the stop criterion, i.e., the number of iterations. If not, all these particles are simulated, and each of these is cost-computed, compared against past best, as well as compared against the Pareto front, and these are updated when the current solution dominates the past solution. One of these is chosen as a leader based on random coefficients of ( $r_1$, $r_2$, Eq. (1)), and each of these velocities is updated for $K_p$, $K_d$, and $K_i$.

The new parameters are constrained using given limits of $\left[K_{pmin},K_{pmax}\right]$, $\left[K_{dmin},K_{dmax}\right]$ and $\left[K_{imin},K_{imax}\right]$ for controlling parameters. The counter is incremented. At the termination of optimization, a Pareto optimal point is selected. The optimized gains are input in each structure of a Lyapunov PID based. Performance metrics are recorded for evaluating the optimization of the controller. Finally, it terminates, which concludes MOPSO optimization.

Simulation setup

All the simulations for the present work were executed on MATLAB R2023a (The MathWorks, Inc., Natick, MA, USA). The simulation environment was designed using MATLAB script files, exploiting maximum use of its computational and visualization facilities. The simulations were executed in a Windows 11 environment installed on an Intel Core i9 CPU and 32 GB RAM. The simulation protocol followed an ordered sequence of the following steps.

Trajectory generation: A gear-like and rose curve reference trajectories of the end effector of the robot in the ZX plane were designed such that all the target points are within the reachable workspace of the manipulator.

Inverse kinematics: The joint angles for every point in the trajectory are computed using a special analytical inverse kinematic function, specially designed for the two-link planar manipulator.

Modeling and implementation of dynamic control: The nonlinear dynamic model of the two-link robot was designed, considering physical parameters, the effect of friction, and external agencies. A Lyapunov-based PID controller has been employed in the simulation loop to track of the trajectories.

Tuning of the controller using MOPSO: Multi-objective particle swarm optimization algorithm is used to fine-tune of the PID gains ( $K_p,K_d$, $K_i$) in the Lyapunov–PID controller. The optimization is for the simultaneous minimization of the tracking error and the control effort by the composite cost function.

Simulation and performance measure: The closed loop system is simulated for the total trajectory. The performance parameters, such as mean squared error (MSE), overshoot, steady state error, and actuator torque, are employed to quantify the performance. Also, they are plotted for pictorial representation. A schematic of simulation architecture is presented in Fig. S2 in the Supplemental Files section. It shows the key modules and the flow of information.

The simulation parameters are shown in Table 2 (Craig, 2022; Spong, Hutchinson & Vidyasagar, 2020; Sciavicco & Siciliano, 2000).

Three distinct paths have been taken to evaluate the performance of the controller (the spiral and rose curves are included in the Supplemental Files). Gear path, with sharp intersection and rotational trends. Let $n$ be the number of gear teeth, $\alpha$ the tooth fraction, $N_{arc}$ the number of points for the base arc, $R_1$ the root (base) radius, $R_2$ the tip radius, and $\left(c_x,c_z\right)$ the center coordinates of the gear.

Let $j=0,1,\dots ,n-1$ index the gear teeth. For each $j$, the following parameters are defined:

(29) $$\mathrm{\Delta }\theta =\frac{2\pi }{n},{\theta }_j=j\mathrm{\Delta }\theta ,{\theta }_{j,tip}={\theta }_j+\alpha \mathrm{\Delta }\theta ,{\theta }_{j,end}={\theta }_j+\mathrm{\Delta }\theta .$$The gear shaped trajectory is constructed by sequentially connecting the following segments for each tooth. Radial ascent: From base radius $R_1$ to tip radius $R_2$ at angle ${\theta }_j$. Tip arc: an arc of radius $R_2$ spanning the angular range $\left[{\theta }_j,{\theta }_{j,tip}\right]$. Radial descent: From tip radius $R_2$ back to the base radius $R_1$ at angle ${\theta }_{j,tip}$. Base arch. An arc of radius $R_1$ spanning the angular range $\left[{\theta }_{j,tip},{\theta }_{j,end}\right]$ discretized into $N_{arc}$ points for smoothness. Each segment is converted to Cartesian coordinates via (Stewart, Clegg & Watson, 2021).

(30) $$x=c_x+rsin\theta ,$$

(31) $$z=c_{z}rcos\theta ,$$where $\left(c_x,c_z\right)$ denote the trajectory center, and $r$ and $\theta$ are the segment specific radius and angle.

Performance metrics

Appropriate procedures were implemented to maintain fairness, such as testing all the controllers using the same system dynamics, initial conditions, and perturbation scenarios. To maintain an open assessment of all the control strategies implemented within the present research work, the next methodology is followed. First, the mathematical representations and implementation aspects of all the controllers have been provided in fine details, along with the structure of the controllers, the parameter settings, and algorithmic routines of the controllers. The same conditions of the system dynamics, initial states, and sequences of perturbations, illustrated in the experimental setup, have been used for comparing all the controllers. The same set of performance metrics are employed for comparison to facilitate direct and objective evaluation. They comprise classical control metrics e.g., overshoot, steady state error, and mean squared error as well as some additional metrics of special interest in robotics, such as the total control effort, maximum per joint torque, and recovery time after perturbation. Each of the metrics is defined rigorously in mathematics as follows (Song, Li & Zhang, 2024; Siciliano, Sciavicco & Oriolo, 2021; Rivera et al., 2020; Elsisi et al., 2021).

Overshoot (OS) measures the maximum amount by which the system output exceeds the desired value following a reference change, expressed as a percentage of the total reference change. It indicates the tendency of the system to “overreact” before setting.

(32) $$\mathrm{O}\mathrm{S}={\displaystyle \frac{y_{max}-r_{final}}{\left|r_{final}-r\left(0\right)\right|}\cdot 100\mathrm{\%},}$$where $y_{max}$ is the output value after the step change or reference shift, $r_{final}$ is the value attained by the system output after the reference shift, $r\left(0\right)$ is the initial reference value (prior to the reference change).

Steady-State Error $\left(e_{ss}\right)$ is the net error between the output and the reference after the system is in steady state, i.e., after all the transient effects have faded. It indicates the ultimate precision of the system.

(33) $$e_{ss}=\underset{t\to \mathrm{\infty }}{lim}\left|e\left(t\right)\right|\approx {\displaystyle \frac{1}{T_{ss}}{\int }_{t_{f-T_{ss}}}^{t_f}\left|e\left(t\right)\right|dt,}$$where $e\left(t\right)$ is the error in time $t$, $T_{ss}$ is the steady-state interval, $t_f$ is the final time of the simulation/experiment.

MSE quantifies the average of the squared tracking error over the entire operation or simulation. This metric globally assesses tracking accuracy, with larger errors weighted more heavily.

(34) $$MSE=\frac{1}{N}\sum _{i=1}^N{\left|r\left(t_i\right)-y\left(t_i\right)\right|}^2,$$where $N$ is the number of samples, $r\left(t_i\right)$ is the required value or path of the $i-th$ sample, $y\left(t_i\right)$ is the true system output or the sampled value of the system output of the $i-th$ sample, $t_i$is the value of the $i-th$ sample time.

Total Control Effort $\left(U_{tot}\right)$ is the total of the actuators’ effort over time and is the integral (or sum) of the absolute value of the control inputs. It is the “energy cost” of control.

(35) $$U_{tot}=\sum _{i=1}^N\sum _{j=1}^m\lceil u_j\left(t_i\right)\rceil \cdot \mathrm{\Delta }t,$$where $m$ is the number of joints, $N$ is the number of time steps, $u_j\left(t_i\right)$ is the controller output (torque command to joint $j$ at time $t_i$, $\mathrm{\Delta }t$ is the size of the time step.

Maximum Torque per Joint $\left({\tau }_{max,j}\right)$ encompasses the experiment or simulation specified maximum value of torque being put on the robot’s individual joints. It is used to measure actuator requirements and to ensure the operational constraints are never violated.

${\tau }_{max,j}={max}_{1\le i\le N}\left|{\tau }_j\left(t_i\right)\right|$, where ${\tau }_j\left(t_i\right)$ is the control torque at time $t_i$ corresponding to the joint $j$.

Recovery Time after Disturbance $\left(t_{rec,X},t_{rec,Z}\right)$ is the time needed for the system to reach again a small acceptable error (say, two times the static error) after the occurrence of an external disturbance, considered independently for the Z and X axes. It is an indication of the robustness and the disturbance rejection of the system. For the given disturbance applied between $t_{dist}^{start}$ to $t_{dist}^{end}$.

For axis X:

(36) $$t_{rec,X}=min\left\{t>t_{dist}^{end}:\left|e_X\left(t\right)\right|<2\cdot e_{ss,X}\right\}-t_{dist}^{end},$$for axis Z:

(37) $$t_{rec,Z}=min\left\{t>t_{dist}^{end}:\left|e_Z\left(t\right)\right|<2\cdot e_{ss,Z}\right\}-t_{dist}^{end},$$where $t_{rec,X}$, $t_{rec,Z}$ are the recovery time after perturbation for the axes Z and X, respectively; $t_{dist}^{end}$ and $t_{dist}^{end}$ are the initial and the end time, respectively of the applied perturbation; $\left|e_X\left(t\right)\right|$, $\left|e_Z\left(t\right)\right|$ are the absolute value of the tracking error for the axes Z and X, respectively at moment $t$; $e_{ss,X}$ and $e_{ss,Z}$ are the steady-state error for the axes Z and X, respectively.

To verify the performance for the optimization of the MOPSO, the nonoptimized Lyapunov-based PID and the PID controller are also contrasted. The PID controller and the Lyapunov PID applied the initial tuning with preliminary experiments resulting in constant $K_p$, $K_i$, and $K_d$ gains. Comparing the simulation results reveals significant improvements in controller performance from MOPSO, underscoring the value of using MOPSO to design effective robot control strategies. The necessity of a systematic and scientific method of selecting the gains of the PID and of stability analysis according to Lyapunov based on the theory of robot manipulators has also been emphasized in literature. For example, Hernández-Guzmán & Orrante-Sakanassi (2019) presented a globally asymptotically stable scheme of PID for multi-joint manipulators using Lyapunov theory for which they demonstrated stability and, through simulations, verified the performance of the control with properly selected gains of the PID. Their findings characterize support for the necessity of stability-guaranteed PID, as well as for introducing methodical selection techniques of gains, like in the current work, using MOPSO.

In addition, Amiri & Ramli (2022) developed a Lyapunov-based adaptive PID controller for multi jointed robotic arms, comparing several gain initialization strategies and showing that correct gain selection enhances simultaneously tracking performance and robustness of the system. Similarly, Qiao et al. (2021) developed two adaptive PID controllers for robotic manipulators that ensure local asymptotic convergence of joint position and velocity tracking errors without requiring any equality or inequality constraints on the control gains. Unlike classical and other adaptive PID approaches, their proposed controllers offer enhanced robustness against dynamic uncertainties, even when the upper bound of such disturbances is unknown. Both controllers were analyzed using Lyapunov stability theory, which formally proved closed-loop stability and boundedness of the parameter estimation errors. Simulations on a two-degree-of-freedom planar manipulator demonstrated superior performance in tracking accuracy and disturbance rejection compared to classical PID, adaptive PD, LADRC, and NDOB-PID controllers. Their findings also highlight the potential importance in including advanced optimization and stability analysis, as discussed in the current contribution, for obtaining robust and high-performance controllers for robotic manipulators. By combining Lyapunov-based analysis with multi-objective optimization, the current contribution moves towards a systematic and dependable controller design procedure for nonlinear robotic systems, surpassing the limitations of heuristic or solely trial and error adjustment procedures.

Gear shape trajectory

The comparative results demonstrate performance improvement when utilizing the MOPSO-Lyapunov-PID controller over results utilizing Lyapunov-PID and conventional PID control schemes. Specifically, steady-state errors of the X and Z axes are minimized by at least 43–50%, resulting in improved accuracy for trajectory tracing. MOPSO-Lyapunov-PID achieves notable accuracy gains, the MSE decreases by up to 93% in X and 77% in Z, indicating closer alignment with the desired trajectory. Transient behavior is also improved, overshoot is reduced by more than 80% in X and by 57–64% in Z relative to non-optimized controllers, yielding a smoother and more stable response. These improvements are achieved at a modest cost, the total control effort increases by only ~1.6%, a minor trade-off given the pronounced improvements in tracking accuracy and robustness. The recovery rates following the occurrence of a disturbance are significantly enhanced by the MOPSO-Lyapunov-PID, after the introduction of a disturbance in X, they are roughly 85–90% faster, and in Z roughly 90% faster than for the other schemes.

The simulation plots of the track of gear shape are given in Fig. 3. This figure shows the desired trajectory and the subsequent trajectory for three controllers: (a) the optimized Lyapunov-PID, (b) the non-optimized Lyapunov-PID, and (c) the PID-only controller. The desired trajectory is represented by dashed lines, and the subsequent trajectories of the non-optimized and optimized controllers can be observed from this figure.

The Pareto front for the gear-shaped path is illustrated in Fig. 4. It presents the Pareto optimal trade-off between tracking error (ET) and control effort (EC) achieved through the MOPSO method. The solutions that lie within this frontier represent a range of trade-offs, allowing for flexibility based on the specific requirements of the application.

In the following descriptions the principal findings are discussed.

Figure 5A presents the tracking performance metrics of a gear-shaped trajectory by a Lyapunov-based PID controller optimized by MOPSO. The cumulative mean squared error in the X and Z directions, shown in the upper plots, quickly decreases and stabilizes around zero, revealing successful convergence of the trajectories. The overshoot in the X and Z directions, as shown in the lower plot, remains continuously at its lowest value, indicating high accuracy control with negligible deviation from the intended track. Figure 5B shows the Joint torque patterns for gear shaped end-effector trajectory tracking with the optimized Lyapunov PID controller. Joint 1 has larger magnitudes of the torque, in conformity to its functionality in producing the dominant motion at the end-effector, where the other Joint 2 has lower magnitudes of the torque.

Figure 6A depicts the force response tracking the gear shaped path by the PID controller optimized by MOPSO. Top plot: the norm on the resultant end-effector force F; middle plot: the X axis components; bottom plot: the z axis components. The force magnitudes remain consistently below 50 N during the task. The force component in the X direction varies with the curved portions in the path, and the component for Z indicates vertical motion and support. The result is the stable dynamic interaction force during operation by the controller. Figure 6B depicts the steady-state tracking errors in the X and Z directions by the optimized controller for the gear shaped trajectory. Both the components of the error drop drastically when the trajectory settles down and remains on the lower side for the simulation, with excursions less than 3–4 mm. The result confirms the controller in maintaining systematic tracking errors at minimum, the evidence for its robustness and precision in sustained dynamic motion.

Figure 7A presents the metrics of performance during the tracking of gear shaped trajectory by the application of the Lyapunov-based PID controller in the absence of optimization. The first row depicts the cumulative MSE in the Z and X directions. The error reaches its peak and reduces gradually, indicating convergence but at a slower pace than the optimized one. The second row presents the overshooting in both directions, where it is more at the start of the trajectory and stabilizes in the later stages.

Figure 7B shows the joint torque profiles for the gear-shaped pattern with the unoptimized Lyapunov PID controller. Joint 1 has torques in the range of approximately 20 to 40 Nm, and Joint 2 is in ±10 Nm. These observations highlight the trade-off between the simplicity of the controller and its performance when tuning is not accomplished through optimization.

Figures 8A and 8B show the end effector force response and steady state tracking error, respectively, for the implementation of the gear-shaped trajectory by using the non-optimized version of the Lyapunov-based PID controller.

The following figures show the performance of the robotic manipulator during the tracking of a gear shaped trajectory using a classical PID controller. Figure S3A in the Supplemental Files shows the performance of the robotic manipulator during the tracking of a gear-shaped trajectory using a classical PID controller. The cumulative MSE and overshoot in the X and Z directions are presented. As seen, the tracking error starts with high peaks and gradually decreases over time. Overshoot is also more pronounced, especially during the initial moments of motion. Figure S3B in the Supplemental Files depicts the joint torques of the two actuators.

Figures S4A and S4B in the Supplemental Files present the performance of the robotic manipulator during the execution of a gear-shaped trajectory using a classical PID controller. Trajectory (a) depicts the end effector force response as a function of time. The magnitude of the total force $(\left|F_{ee}\right|$) oscillates between about 10 and 35 N during the task. The X component of the force $\left(\left|F_{eex}\right|\right)$ oscillates between 20 and 10 N, exhibiting a distinct modulation structure which adheres to the geometry of the trajectory. Along the Z direction ( $\left|F_{eez}\right|)$, the force oscillates between 0 N and about 30 N, exhibiting relatively consistent behavior with slight oscillations as a function of time. Figure S4B shows the steady-state tracking error along the two cartesian axes. The X error peaks at around 0.022 m during the first few seconds and has a decline, settling mainly below 0.003 m after 50 s. Analogously, the Z error peaks at around 0.007 m at the beginning and stays mainly below 0.002–0.003 m throughout the rest of the path duration.

Table 3 gives a quantitative comparison between the three considered control schemes: MOPSO optimized Lyapunov PID, non-optimized Lyapunov PID, and conventional PID for gear shape. Typical performance metrics such as steady-state error, MSE, overshoot, total control input, max torque, as well as recovery time from disturbances are provided to critically compare and evaluate each scheme’s weaknesses and strengths.