The target platform is configured as a planar 2-DoF lower limb rehabilitation exoskeleton operating in the sagittal plane. Hip flexion extension is denoted by θ 1 and knee flexion extension by θ 2 ; the joint coordinate vector is θ = θ 1   θ 2 T with velocities θ ˙ and accelerations θ ¨ . Rigid thigh and shank frames are attached to the patient via adjustable cuffs and quick-release straps, and link lengths are adjusted to the user’s anthropometry. Typical range of motion envelopes are considered to guide control and safety limits ( h i p ≈ − 20 °   t o + 120 ° , k n e e ≈ 0 °   t o + 130 ° ), and mechanical end-stops are provided near the extremes to prevent over-travel. Each joint is actuated by an electric drive with high-ratio transmission (e.g., harmonic or planetary gearing), by which motor torque is amplified while reflected inertia is kept within clinically acceptable limits. Optional series elasticity or software torque limiting is employed to improve comfort during therapy. Joint-level torque and speed limits are enforced in firmware to maintain operation within safe bounds consistent with the limits used in the results section. Joint angles are measured using high-resolution absolute encoders, and joint velocities are obtained by numerical differentiation with appropriate filtering. Drive currents are monitored for torque estimation; inertial measurement units may be mounted on the thigh and shank for segment level orientation, and foot-contact sensing (e.g., insole force sensors) can be used for gait-phase or state detection. All safety critical signals (emergency stop, over-current, over-temperature) are handled by hardware interlocks in parallel with software supervision. Control is executed on a real-time embedded controller. A cascaded structure is adopted: an inner current/torque loop runs at high frequency to regulate actuator torque; a joint-position loop implements the MFOPID at a lower, yet real-time, rate; and a high-level trajectory generator with a safety supervisor coordinates task execution and enforces limits. The MFOPID parameters (gains and fractional orders) are tuned offline using a hybrid metaheuristic Improved Elk Herd Optimization hybridized with Grey Wolf and Multi-Verse Optimization algorithms (IElk-GM) so that exploration and exploitation are balanced during tuning while runtime complexity remains minimal.

The LLRER considered in this study is a planar 2-DoF structure consisting of two rigid links and two revolute joints, which correspond to the hip and knee joints of the human body. This configuration is designed to facilitate flexion and extension movement in the sagittal plane, thereby enabling gait rehabilitation for stroke and mobility-impaired patients (Al Rezage and Tokhi, 2016). The mechanical structure of the robot is shown in Figure 1. A dynamic model based on the anthropometric features of a human lower limb is used to describe the mobility of the robot. The model assumes a subject with a body mass of 74 kg and height of 1.69 m, with segment properties obtained from winter’s anthropometric data (Alshatti, 2019; Winter 2009). The robot dynamics are derived using the Euler-Lagrange method, capturing the effects of joint inertia, Coriolis and centrifugal forces, gravitational torque, control inputs, and external disturbances.

The general 2-DoF dynamics are given in Equation 1: M θ θ ¨ + C θ , θ ˙ θ ˙ + G θ + d t = u t (1)

The angle, angular velocity, and acceleration of a robot joint vector are denoted by the variables θ , θ ˙ , and θ ¨ , respectively. For every inertia M θ , Coriolis, and centrifugal torque C (θ, θ ˙ ) ∈ R 2 × 2 in human limb matrices. The one-dimensional vector of the torque of gravity G θ is ∈ R 2 × 1 , the control signal is represented by u t , and the vector of external disturbances is d t ∈ R n .

The dynamics of the robot are indicated by Equation 2: M 11 M 12 M 21 M 22 θ ¨ 1 θ ¨ 2 + C 11 C 12 C 21 C 22 θ ˙ 1 θ ˙ 2 + G 1 G 2 = u t 1 u t 2 (2)

The components of the inertia matrix M(θ) are depicted in Equation 3: M 11 = I 1 + I 2 + m 1 L C 1 2 + m 2 L 1 2 + m 2 L C 2 2 + 2 m 2 L 1 L C 2   cos θ 2 M 12 = M 21 = I 2 + m 2 L C 2 2 + m 2 L 1 L C 2   ⁡ cos θ 2 M 22 = I 2 + m 2 L C 2 2 (3)

The elements of C θ , θ ˙ are determined by Equation 4: C 11 = − m 2 L 1 L C 2   sin θ 2 θ ˙ 2 C 12 = − m 2 L 1 L C 2   sin θ 2 θ ˙ 1 + θ ˙ 2 C 21 = m 2 L 1 L C 2   sin θ 2 θ ˙ 1 C 22 = 0 (4)

The parameters of the gravitational vector G θ are specified by Equation 5: G 1 = m 1 L C 1 g   sin θ 1 + m 2 g L 1   sin θ 1 + m 2 g L C 2 ⁡ sin θ 1 + θ 2 G 2 = m 2 g L C 2 ⁡ cos θ 1 + θ 2 (5)

The variables of these equations are delineated by specific parameters presented in Table 1.

To improve performance and efficiency, a modified fractional order PID controller is suggested to reduce overshoot and steady-state error. Equation 13 shows the modified alpha function from Equation 6: γ L t = H P L e L t + H I L ∫ e L t d t + H D L e ˙ L t (13)

Accordingly, Equation 8 is to be modified in Equation 14 for the control law: u m o d L t = M θ , θ ˙ u F O P I D L t + u a u x t (14) where Equation 15 defines u a u x t : u a u x t = tanh   γ L t * K m o d L (15)

The IElk-GM algorithm is utilized to get system parameters for the MFOPID controller. The control goals are accomplished by the way this algorithm’s fitness function is set up. The candidate Lyapunov function is given by Equations 16, 17 as follows: V L t = 1 2   γ L 2 t (16) V L ˙ t = γ L t γ L ˙ , w h e r e   t < 0 (17) when V L < 0 it is guaranteed that γ L → 0 as t → ∞ . From Equation 13, it can be written: γ L ˙ t = H P L e ˙ L t + H I L ∫ e ˙ L t d t + H D L e L ¨ t (18)

Substituting Equation 18 into Equation 17: V L ˙ t = γ L t H P L e ˙ L t + H I L ∫ e ˙ L t d t + H D L e L ¨ t < 0 (19) The negativity condition used in the stability proof is expressed in Equation 19.

Since V L t ≥ 0 and V L ˙ t ≤ 0 , hence, according to the Lyapunov direct method, the system is Lyapunov globally stable. Moreover, due to the structure of γ L t , the convergence of e L t → 0 is smooth and without overshoot, as confirmed in simulation.

The IElk-GM algorithm is also used to calculate the optimal parameters of the controller u mod t of link1 ( H P 1 , H I 1 , H D 1 , β 11 , β 21 , β 31 , λ 1 , μ 1 , a n d   K mod ⁡ 1 , and link2 ( H P 2 , H I 2 , H D 2 , β 12 , β 22 , β 32 , λ 2 , μ 2 , a n d   K mod ⁡ 2 ) .

The nominal Lyapunov analysis in Section 4 establishes stability under ideal conditions. In rehabilitation, however, three non-ideal effects are unavoidable: (i) actuator saturation (torque limits and anti-windup), (ii) small I/O delays from sensing/actuation, and (iii) patient-induced disturbances (matched torques at the joints). Let u c denote the MFOPID control before limits and u = s a t u c the applied torque after saturation. Let Δ u : = u − u c be the saturation mismatch; τ d t bounded patient torque disturbances; and τ ∈ 0 , τ max a constant small delay in the loop (sensing/actuation). We use the same Lyapunov candidate V · as in Section 4 and the same error vector e . Closed-loop maps are locally Lipschitz.

The proposed IElk-GM algorithm is a hybrid metaheuristic that integrates three nature-inspired optimization strategies to balance global exploration and local exploitation.

• Improved Elk Herd Optimization (IElk-GM): The algorithm forms the backbone of the search process, offering enhanced population diversity and a structured Herd-based exploration mechanism. Unlike the standard EHO, the improved version incorporates elitism preservation, adaptive parameter control, and balanced harem assignment, which significantly enhance convergence speed and solution quality.

This hybridization of IEHO, GWO, and MVO capitalizes on the strengths of each algorithm: diversity, leadership-based exploitation, and randomness, resulting in a more robust and efficient optimization framework for tuning complex control parameters.

For more information about the Elk Herd Optimization Algorithm, the Grey Wolf Optimization Algorithm, and the Multi-Verse Optimization Algorithm, see (Al-Betar et al., 2024; Mirjalili et al., 2014; Mirjalili et al., 2016).

The overall working process of the proposed IElk-GM algorithm is presented in Figure 3.

The step response performance of the controlled LLRER subjected to a positive unit step input at link 1 (hip joint) and a negative unit step input at link 2 (knee joint) is illustrated in Figures 5, 6 for both FOPID and MFOPID controllers.

The results demonstrate that the MFOPID controller significantly enhances system performance, enabling the robot to accurately follow the desired trajectory with fast transient response, zero overshoot, and negligible steady-state tracking error. Specifically, the settling times are reduced (from 6.998 s to 0.430 s for the hip joint and from 7.150 s to 0.829 s for the knee joint). Moreover, the control signals remain smooth and within acceptable torque limits of less than 345 Nm for link 1 and less than 95 Nm for link 2. In contrast, the FOPID controller exhibits slower convergence, minor overshoot, and less stable control signals, indicating inferior performance under the same conditions. Table 5 lists the simulation results’ evaluation parameters for the FOPID and MFOPID.

The simulation results of the LLRER using the FOPID and MFOPID, tested with the non-linear cosine input signal ( x d 1 = π 4 + 1 − cos ⁡ 3 t for link1 and ( x d 2 = π 6 + 1 − cos ⁡ 5 t for link2, are shown in Figures 7, 8; these results demonstrate the performance to be reliable, despite the non-linearity of the input signal. The results demonstrate excellent performance parameters, negligible error, and a smooth control signal (less than 350 Nm for link 1 and less than 100 Nm for link 2) with MFOPID. In contrast, the standard FOPID controller exhibits inefficient tracking, with visible deviations from the reference trajectory and more oscillatory control actions, indicating its limited robustness under non-linear operating conditions.

To further validate the performance of the proposed MFOPID controller optimized via the IElk-GM algorithm, we conducted a direct comparison with the Adaptive Optimal Fractional-order Super-Twisting Sliding Mode (AOFSTSM) controller optimized using GWO, as proposed by Faraj et al. (2023). We benchmarked MFOPID + IElk-GM against AOFSTSM + GWO along four dimensions: robustness, control-signal smoothness, constraint handling, and computational burden.

i. Robustness: AOFSTSM attains high disturbance rejection via super-twisting on a fractional sliding surface and adaptive bound estimation, as reported in Faraj et al. (Faraj et al., 2023), whereas MFOPID achieves comparable tracking envelopes under parametric variations using smoothly shaped fractional actions without signum-type injections.

In summary, the AOFSTSM remains preferable when maximal invariance to severe uncertainties is paramount, whereas MFOPID + IElk-GM is advantageous when smooth control, embedded simplicity, and energy/comfort metrics are prioritized while maintaining competitive tracking accuracy. Table 6 summarizes the standard time-domain indices and control signals for hip and knee joints. MFOPID values come from our step-response results (Table 5), and AOFSTSM values are obtained from Faraj et al.‘s published plots by careful digitization (their paper does not tabulate these time-indices explicitly). In the last column of Table 6 we additionally report the percent improvement of MFOPID relative to the AOFSTSMC baseline, computed for lower-is-better indices as I m p r o v e m e n t   % = M F O P I D − A O F S T S M A O F S T S M * 100   %

Positive values indicate that MFOPID is lower/better than the AOFSTSMC baseline; negative values indicate the opposite. The same improvements are shown in Figure 9 for clarity.

The MFOPID was tuned with an ITAE-centric objective plus mild penalties on overshoot and torque bounds to favor comfort-oriented transients and smooth actuation. On Link 2 (knee), the combination of nonlinearities and anti-windup under actuator limits attenuates the effective integral action near steady state, yielding a small residual offset ( e s . s ≈ 0.011 rad ≈ 0.63 ° ). This trade-off is intentional, it yields zero overshoot, short/competitive settling times, and non-switching, low-ripple control torques (see Tables 6–8), which are clinically more relevant for repeated therapy than pushing e s . s to machine precision. If needed, the offset can be further reduced without altering the main conclusions by (i) slightly increasing the integral weight for Link 2 within the same torque limits, (ii) adding light feedforward compensation (gravity/friction), or (iii) adopting a 2-DoF set-point weighting that preserves the overshoot-free transient while tightening steady-state accuracy.

The data in Table 6 quantitatively substantiates the qualitative contrast above. For smoothness and comfort, MFOPID + IElk-GM yields overshoot-free transients in both joints with short settling times ( t s = 0.430   sec ⁡ . for hip; t s = 0.829   sec ⁡ . for knee), despite the absence of any switching terms. Tracking accuracy remains competitive: steady-state errors are small for both controllers; in the hip joint, AOFSTSM achieves a slightly lower e s . s (as expected for a sliding-mode design), while MFOPID maintains comparable accuracy without chattering mechanisms. Regarding control effort, MFOPID’s step-test peaks are bounded ( ≤ 345 N.m for hip; ≤90 N m for knee) and decay rapidly, whereas the peaks read for AOFSTSM in the uncertainty-case plots are ≈ 80 N.m for hip, and ≈ 60 N.m for knee; because the experimental contexts differ, these magnitudes are reported for completeness rather than as a like-for-like torque comparison, and the time-domain indices ( M P , t r , t s , e s . s ) should be taken as the primary evidence. Finally, the implementation burden of MFOPID remains lower due to its fixed, continuous fractional-PID structure with offline IElk-GM tuning, which simplifies embedded deployment while preserving the favorable transients summarized above.

Beyond the indices in Table 6, the benefits of MFOPID + IElk-GM can be assessed without relying on unavailable numeric data from (Faraj et al., 2023). First, MFOPID is continuous and non-switching, which is reflected in smoother torque traces (no high-frequency flicker or sharp corners) and is desirable for patient comfort and actuator longevity. Second, the overshoot-free transients observed in our plots indicate comfort-oriented behavior while maintaining competitive steady-state accuracy. Third, implementation is lighter: MFOPID uses a fixed-structure fractional PID tuned offline, whereas AOFSTSM requires super-twisting and online adaptive bounds. These points are summarized in Tables 7, 8.

Beyond the FOPID, we include two hybrid-optimized classical FOPID baselines: FOPID + PSO-GWO (PSO phase for global exploration followed by GWO refinement) and FOPID + GWO-MVO (GWO coarse search followed by MVO fine search). All controllers are evaluated under the same plant, references, actuator limits, and disturbance/uncertainty scenarios. The composite objective is ITAE with mild penalties on overshoot, settling time, and torque bounds; the optimization budget (population size × iterations) is matched across methods. We report settling time, overshoot, rising time, steady-state error, and control torque. For readability, we also provide percent improvements of MFOPID relative to each baseline using: I m p r o v e m e n t   % = M F O P I D − B a s e l i n e B a s e l i n e * 100   %

For all lower-is-better indices. Positive values indicate MFOPID is lower/better. Additionally, Tables 9, 10 summarize the robustness results, while Figures 10–17 depict the joint-position and control-torque responses for linear and nonlinear trajectories.