A Realistic Model Reference Computed Torque Control Strategy for Human Lower Limb Exoskeletons

1. Introduction

Physical disability significantly impairs a person’s mobility and physical capabilities, affecting both individuals and society. Nearly a billion people worldwide face challenges in accessing adequate treatment due to limited resources. Exoskeleton robots offer a promising solution for rehabilitation, providing personalized physiotherapy and continuous performance assessment throughout recovery. This technology holds great potential for enhancing physical therapy and improving neurorehabilitation outcomes for individuals with disabilities affecting both upper and lower extremities.

Effective rehabilitation with exoskeleton robots requires precise maneuverability, which depends on an integrated control system combining sensing, actuation, and computation units that communicate via standard protocols and are driven by control algorithms. The high-performance control execution server functions as the brain of the rehabilitation robot, enabling it to provide passive, active, active-assist, and active-resist therapies. Sensors monitor the robot and facilitate interaction with the user, while actuators, functioning like muscles, move the mechanical components. However, even advanced sensors and actuators are insufficient without efficient control algorithms to ensure proper robot maneuverability. Researchers have developed various linear and nonlinear control algorithms to enhance maneuverability and optimize the delivery of diverse physical therapies.

Table 1 summarizes the features of recent human lower limb exoskeleton robots, detailing their actuated joints, control algorithms, and key characteristics.

Table 1 illustrates that both linear and nonlinear control algorithms are utilized to address nonlinear robot dynamics. Most nonlinear control algorithms rely on a mathematical model of the robot’s dynamics, making it essential to create an accurate model as part of the control architecture. Although PD and PID techniques do not strictly require a model, having one aids in gain selection and performance evaluation. Developing an accurate robot model is challenging, but feedback control mechanisms offer some protection against modeling inaccuracies. The acceptable modeling error range varies by control technique, with some algorithms being sensitive to errors while others are more robust against uncertainties [33].

This paper presents a 7-DOF dynamic model of the human lower extremity developed using the Lagrange energy method. To control this model, a novel Realistic Model Reference Computed Torque Controller (RMRCTC) is introduced, providing enhanced protection against discrepancies between plant and model parameters. Additionally, a realistic friction model is incorporated to simulate dynamic frictional effects. The paper is organized into seven sections. Section two covers human lower extremity anatomy, which is crucial for anthropomorphic modeling of the exoskeleton robot. Section three describes the development of kinematic and dynamic models and friction modeling. Section four details the Realistic Model Reference Computed Torque control architecture and controller stability analysis. Section five presents the simulation results along with a discussion, while section six verifies the robustness of the developed controller concerning modeling discrepancies. Finally, section eight concludes the research with summary remarks.

2. Human Lower Extremity Anatomy

The human lower extremity is a complex structure comprising bones, joints, and major muscles that work together to facilitate movement and support the body. The primary bones include the femur (thigh bone), patella (kneecap), tibia (shin bone), and fibula. The hip joint connects the femur to the pelvis and allows for a wide range of motion in multiple planes, providing approximately three degrees of freedom (flexion/extension, abduction/adduction, and internal/external rotation). The knee joint, formed by the femur, tibia, and patella, enables flexion and extension and internal-external rotation, contributing two degrees of freedom. The ankle joint connects the tibia and fibula to the foot, allowing for dorsiflexion and plantarflexion, with additional slight movements of inversion and eversion.

Key muscles include the quadriceps, which extend the knee; the hamstrings, responsible for knee flexion; and the gastrocnemius and soleus muscles in the calf, which facilitate plantarflexion of the foot. The gluteal muscles, including the gluteus maximus, medius, and minimus, play vital roles in hip stabilization and movement. Together, this intricate anatomy supports weight-bearing, locomotion, and balance, allowing for efficient and coordinated movement during various activities. Figure 1 shows the human lower extremity degrees of freedoms.

2.1. Human Lower Extremity Ranges of Motion

The range of motion in the human lower extremity is essential for determining the joint ranges of motion in exoskeleton robots. Most healthy individuals exhibit similar ranges of motion. However, for physical therapy patients, the acceptable range values differ based on each individual’s health condition. Average joint ranges of motion have been documented from various sources [34,35,36].

The human lower extremity exhibits a wide range of motion across various joints, essential for mobility and functionality. The hip joint allows approximately 120 degrees of flexion, 30 degrees of extension, 45 degrees of abduction, and 30 degrees of adduction, along with internal and external rotation (both are 45 degrees). The knee joint permits about 135 degrees of flexion and limited extension. The ankle joint facilitates dorsiflexion of around 20 degrees and plantarflexion of approximately 50 degrees, along with inversion/pronation and eversion/supination movements (35 degrees and 15 degrees). This extensive range of motion enables activities such as walking, running, jumping, and climbing, contributing to overall mobility and stability.

2.2. Human Lower Extremity Anthropometric Parameters

Anthropometric properties are essential for developing kinematic and dynamic models of the human lower extremity. Nikolova and Toshev conducted a study analyzing these properties using data from 5,290 individuals, including 2,435 males and 2,855 females, with their findings published in 2008 [37]. Additionally, Contini created empirical equations to determine anthropometric parameters of the human upper and lower extremity based on the subject’s weight and height. For this article, we used the empirical formulas proposed by Contini. Equations (1) to (21) outline the calculations used for determining the anthropometric parameters necessary for dynamic modeling and simulation.

$$B_d=0.6905+0.0297C{\displaystyle \frac{lb}{f{t}^3}},\mathrm{}where\mathrm{}C=H{W}^{-{\displaystyle \frac{1}{3}}},$$

(1)

In Equation (1) $B_d$is the body density ($lbs/f{t}^3)$, $H$ is the height of the subject ($inch)$, $W$ is the body weight of the subject ($lbs)$. The thigh, shank and foot densities are determined by the Equations (2) to (4),

$$T_d=1.035+0.814\ast B_d{\displaystyle \frac{lb}{f{t}^3}}$$

(2)

$$S_d=1.065+B_d{\displaystyle \frac{lb}{f{t}^3}}$$

(3)

$$F_d=1.071+B_d{\displaystyle \frac{lb}{f{t}^3}}$$

(4)

In Equations (2)–(4), $T_d$is the thigh density, $S_d$ is the shank density and $F_d$ is the foot density. The whole-body volume ${(B}_v)$ was calculated using body weight and body density.

$$B_v={\displaystyle \frac{W}{B_d}}f{t}^3$$

(5)

The volume of thigh, shank, and foot were calculated as follows:

$$T_v=0.0922\ast B_{v}f{t}^3$$

(6)

$$S_V=0.0464\ast B_{v}f{t}^3$$

(7)

$$F_v=0.0124\ast B_{v}f{t}^3$$

(8)

Approximate weights of the thigh, shank, and foot were calculated as given using Equations (9)–(11)

$$T_m=T_v\ast T_{d}lbs$$

(9)

$$S_m=S_v\ast S_{d}lbs$$

(10)

$$F_m=F_v\ast F_{d}lbs$$

(11)

The length of the thigh $\left(T_l\right)$, shank $\left(S_l\right)$, foot $\left(F_l\right)$ and ankle to lower face of the foot $\left(A_g\right)$ were calculated as shown in Equations (12)–(15),

$$T_l=0.245\ast Hinch,$$

(12)

$$S_l=0.285\ast Hinch,$$

(13)

$$F_l=0.152\ast Hinch,$$

(14)

$$A_g=0.043\ast Hinch$$

(15)

The locations of the center of the mass from the proximal joint (for thigh $\left(T_{cm}\right)$, shank $\left(S_{cm}\right),$ foot $\left(F_{cm}\right)$) are given by Equations (16)–(18),

$$T_{cm}=0.41\ast T_{l}inch,$$

(16)

$$S_{cm}=0.393\ast S_{l}inch,$$

(17)

$$F_{cm}=0.445\ast F_{l}inch$$

(18)

The empirical equations for the inertial properties of the thigh $T_i,$ shank $S_i,$ and foot $F_i$ are given in Equation (19) to Equation (21),

$$T_i=\left[\begin{array}{ccc}{T}_m{\left(0.124\ast T_l\right)}^2& 0& 0\\ 0& T_m{\left(0.267\ast T_l\right)}^2& 0\\ 0& 0& T_m{\left(0.267\ast T_l\right)}^2\end{array}\right]$$

(19)

$$S_i=\left[\begin{array}{ccc}{S}_m{\left(0.281\ast S_l\right)}^2& 0& 0\\ 0& S_m{\left(0.114\ast S_l\right)}^2& 0\\ 0& 0& S_m{\left(0.275\ast S_l\right)}^2\end{array}\right]$$

(20)

$$F_i=\left[\begin{array}{ccc}{F}_m{\left(0.124\ast F_l\right)}^2& 0& 0\\ 0& F_m{\left(0.245\ast F_l\right)}^2& 0\\ 0& 0& F_m{\left(0.257\ast F_l\right)}^2\end{array}\right]$$

(21)

3. Human Lower Extremity Kinematic and Dynamic Modeling

The human lower extremity exoskeleton robot is designed to be closely attached to the human body, aiming to replicate human joints and movements. The distribution of moment of inertia in the exoskeleton mirrors that of the human lower limb. To maintain a general approach rather than focusing on a specific configuration, the anatomical structure and anthropometric parameters of the human lower limb were used into this simulation. We are actively developing an exoskeleton robot for lower extremity rehabilitation. Figure 5 illustrates the human lower extremity exoskeleton robot currently under development.

Figure 2. Human lower extremity rehabilitation exoskeleton robot.

Figure 2. Human lower extremity rehabilitation exoskeleton robot.

3.1. Kinematic and Dynamic Modeling

The kinematic model was developed using the modified Denavit-Hartenberg (DH) parameters method. A separate link frame was assigned to each degree of freedom. Figure 3 illustrates the link frame assignments according to the modified DH parameters. The modified D-H parameters for each joint are given in Table 2.

The general form of the transformation matrix is presented in Equation (22):

$${}_{\mathrm{I}}{}^{\mathrm{i}-1}\mathrm{T}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\theta }}_{\mathrm{i}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\theta }}_{\mathrm{i}}& 0& {\mathsf{\alpha }}_{\mathrm{i}-1}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\theta }}_{\mathrm{i}}\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\alpha }}_{\mathrm{i}-1}& \mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\theta }}_{\mathrm{i}}\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\alpha }}_{\mathrm{i}-1}& -\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{\mathrm{i}-1}& -\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{\mathrm{i}-1}{\mathrm{d}}_{\mathrm{i}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\theta }}_{\mathrm{i}}\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{\mathrm{i}-1}& \mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\theta }}_{\mathrm{i}}\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\alpha }}_{\mathrm{i}-1}& \mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\alpha }}_{\mathrm{i}-1}& \mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\alpha }}_{\mathrm{i}-1}{\mathrm{d}}_{\mathrm{i}}\\ 0& 0& 0& 1\end{array}\right]$$

(22)

By substituting the modified DH parameter values from Table 3 into Equation. (22), the resulting transformation matrices, Equation (23) to Equation (29), were obtained.

$${}^{0}T_1=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}(\begin{array}{c}{\theta }_1)\end{array}& -sin\left({\theta }_1\right)& 0& 0\\ sin\left({\theta }_1\right)& cos\left({\theta }_1\right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(23)

$${}^{1}T_2=\left[\begin{array}{cccc}sin\left({\theta }_2\right)& cos\left({\theta }_2\right)& 0& 0\\ 0& 0& 0& 0\\ cos\left({\theta }_2\right)& -sin\left({\theta }_2\right)& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(24)

$${}^{2}T_3=\left[\begin{array}{cccc}cos\left({\theta }_3\right)& -sin\left({\theta }_3\right)& 0& 0\\ 0& 0& 0& -l_1\\ -sin\left({\theta }_3\right)& -cos\left({\theta }_3\right)& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(25)

$${}^{3}T_4=\left[\begin{array}{cccc}cos\left({\theta }_4\right)& -sin\left({\theta }_4\right)& 0& 0\\ 0& 0& -1& 0\\ sin\left({\theta }_4\right)& cos\left({\theta }_4\right)& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(26)

$${}^{4}T_5=\left[\begin{array}{cccc}cos\left({\theta }_5\right)& -sin\left({\theta }_5\right)& 0& 0\\ 0& 0& 1& -l_2\\ -sin\left({\theta }_5\right)& -cos\left({\theta }_5\right)& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(27)

$${}^{5}T_6=\left[\begin{array}{cccc}\mathit{sin}\left({\theta }_6\right)& \mathit{cos}\left({\theta }_6\right)& 0& 0\\ 0& 0& -1& 0\\ -\mathit{cos}\left({\theta }_6\right)& \mathit{sin}\left({\theta }_6\right)& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(28)

$${}^{6}T_7=\left[\begin{array}{cccc}cos\left({\theta }_7\right)& -sin\left({\theta }_7\right)& 0& a_1\\ 0& 0& 1& 0\\ -sin\left({\theta }_7\right)& -cos\left({\theta }_7\right)& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(29)

The homogeneous transformation matrix determines the position and orientation of the end frame relative to the base frame. The homogeneous transformation matrix was constructed by multiplying the individual transformation matrices from Equation (23) to Equation (29). For the lower extremity, the homogeneous transformation matrix represents the position and orientation of the foot relative to the hip joint.

$${\displaystyle \genfrac{}{}{0pt}{}{0}{7}}T=\left[{\displaystyle \genfrac{}{}{0pt}{}{0}{1}}T{\displaystyle \genfrac{}{}{0pt}{}{1}{2}}T{\displaystyle \genfrac{}{}{0pt}{}{2}{3}}T{\displaystyle \genfrac{}{}{0pt}{}{3}{4}}T{\displaystyle \genfrac{}{}{0pt}{}{4}{5}}T{\displaystyle \genfrac{}{}{0pt}{}{5}{6}}T{\displaystyle \genfrac{}{}{0pt}{}{6}{7}}T\right]$$

(30)

Lagrange’s method is a widely used mathematical approach for deriving dynamic equations of motion. It involves differentiating energy terms with respect to the system variables and time. The dynamic model of the human lower extremity was developed using Lagrange’s method [38]. The kinetic energy of a link at any given moment can be calculated using Equation (31).

$$k_i=\left[{\displaystyle \frac{1}{2}}{m}_i{.v}_{ci}^T.v_{ci}+{{\displaystyle \frac{1}{2}}}^i{\omega }_i^T{.}^{ci}{I}_i.{\omega }_i\right]$$

(31)

In Equation (31), $m_i$ represents the mass of the link, while $v_{ci}$ and${\omega }_{ci}$ denote the linear and angular velocities of the link at its center of mass. $I_i$ is the moment of inertia of the link about its center of mass. The total kinetic energy of the model can be calculated using Equation (32).

$$k=\sum _{i=1}^n{k}_i$$

(32)

In Equation, (32), $n$ represents the total number of degrees of freedom, $n=7$ for the developed human lower extremity dynamic model.

The potential energy of any link can be calculated using Equation (33),

$$u_i=-m_i{.}^0{g}^T{.}^0{P}_{ci}+u_{ref}$$

(33)

$$u=\sum _{i=1}^n{u}_i$$

(34)

In Equation (33)

, $m_i$ represents the mass of the link, $g$ is the gravitational acceleration, and${}_{0}P_{ci}$ denotes the position of the link’s center of mass relative to the ground reference. Finally, by using the Lagrange energy method, the required joint torque can be determined as follows:

$${\tau }_i=\left[{\displaystyle \frac{d}{dt}}{\displaystyle \frac{\delta k}{\delta {\dot{\theta }}_i}}-{\displaystyle \frac{\delta k}{\delta {\theta }_i}}+{\displaystyle \frac{\delta u}{\delta {\theta }_i}}\right]$$

(35)

The dynamic equations of motion for the robot are expressed in Equation (36),

$${\tau }_{Joint}=\left[M\left(\theta \right)\ddot{\theta }+V\left(\theta ,\dot{\theta }\right)+G\left(\theta \right)\right]$$

(36)

In Equation (36), $M\left(\theta \right)$ is the mass matrix, a $\left(7x7\right)$ symmetric positive definite matrix, $V(\theta ,\dot{\theta })$ is $\left(7x1\right)$ matrix represents the Coriolis and the centripetal terms, and $G\left(\theta \right)$ is $\left(7x1\right)$ matrix representing the gravitational term. ${\tau }_{Joint}$, $\left(7x1\right)$ matrix, represents the joint torque requirements. The robot’s dynamic equation of motion can be rewrite as:

$$\ddot{\theta }=\left[M{\left(\theta \right)}^{-1}\left({\tau }_{Joint}-V\left(\theta ,\dot{\theta }\right)-G\left(\theta \right)\right)\right]$$

(37)

Since $M\left(\theta \right)$ is a positive definite matrix, ${M\left(\theta \right)}^{-1}$ always exists. Figure 4 shows a schematic diagram of the ideal/model robot dynamics without counting the joint frictions.

Friction between two mating parts with relative motion is unavoidable. After considering joint friction torques, the robot dynamical equations become,

$${\tau }_{Joint}=\left[M\left(\theta \right)\ddot{\theta }+V\left(\theta ,\dot{\theta }\right)+G\left(\theta \right)+{\tau }_{friction}\right]$$

(38)

where,

$${\tau }_{friction}=\left[F\left(\dot{\theta }\right)\right]$$

(39)

Equation (38) can also be written as,

$$\ddot{\theta }=\left[M{\left(\theta \right)}^{-1}\left({\tau }_{Joint}-V\left(\theta ,\dot{\theta }\right)-G\left(\theta \right)-F\left(\dot{\theta }\right)\right)\right]$$

(40)

Figure 5 presents the schematic diagram of the robot dynamics with friction effects.

Figure 5. Internal architecture of the physical Robot model.

Figure 5. Internal architecture of the physical Robot model.

3.2. Friction Modeling

In robotic manipulators, links are connected by bearings, transmissions, or seals, where relative motion at the joints leads to the generation of friction forces. These friction forces can account for as much as 50% of the transmitted force or torque [39]. A robust control system is essential for compensating for the effects of friction. The magnitude of friction force or torque depends on various factors, such as surface roughness, lubricant viscosity, transmitted load, temperature, and relative velocity between contact surfaces. Since these parameters are variable, it is challenging to develop a comprehensive theoretical friction model.

Most friction models are empirical and have been successfully used for years. Researchers have described friction through models like the Coulomb friction model, Viscous friction model, Stribeck effect, pre-sliding behavior, stiction phase, and hysteresis effects. Advanced models like the Dahl, LuGre, and Karnopp models [40] also exist. In this paper, a friction model combining the Coulomb, Viscous, and Stribeck effects, equivalent to the LuGre model [41], is used. The next section will explain the formulation of the friction model in detail.

Coulomb friction $\left({\mathit{T}}_{\mathit{C}}\right)$ is characterized by a constant friction torque regardless of the velocity or movement of the system, meaning that it remains uniform at any given time. In contrast, viscous friction $\left({\mathit{T}}_{\mathit{V}}\right)$ produces a resistive torque that is directly proportional to the relative velocity between the surfaces in contact. Stribeck friction $\left({\mathit{T}}_{\mathit{S}}\right)$ describes a more complex behavior observed at low velocities, where there is a noticeable decrease in friction force, resulting in a negative slope in the friction-velocity relationship. Together, these friction models describe the various types of resistive forces that can occur in mechanical systems.

Equation (41) presents the friction model, with Equations (41) to (43) were used to calculate the joint friction torque:

$$T=\sqrt{\left(2e\right)}\left(T_{brk}-T_C\right).\exp \left(-{\left({\displaystyle \frac{\omega }{{\omega }_{St}}}\right)}^2\right).{\displaystyle \frac{\omega }{{\omega }_{St}}}+T_C.\tanh \left({\displaystyle \frac{\omega }{{\omega }_{\mathit{Coul}}}}\right)+f\omega$$

(41)

$${\omega }_{St}={\omega }_{brk}\sqrt{2}$$

(42)

$${\omega }_{Coul}={\displaystyle \frac{{\omega }_{brk}}{10}}$$

(43)

where,

$\mathit{T}$ is the total friction torque

${\mathit{T}}_{\mathit{C}}$ is the Coulomb friction torque

${\mathit{T}}_{\mathit{b}\mathit{r}\mathit{k}}$ is the breakaway friction torque: This is the sum of the Coulomb and Stribeck friction torques near zero velocity, often referred to as breakaway friction.

${\mathit{\omega }}_{\mathit{b}\mathit{r}\mathit{k}}$ is the breakaway friction velocity: The velocity at which Stribeck friction reaches its peak, and the sum of the Stribeck and Coulomb friction equals the breakaway friction torque.

${\mathit{\omega }}_{\mathit{S}\mathit{t}}$ is the Stribeck velocity threshold

${\mathit{\omega }}_{\mathit{C}\mathit{o}\mathit{u}\mathit{l}}$ is the Coulomb velocity threshold

$\mathit{\omega }$ is the input angular velocity

$\mathit{f}$ is the viscous friction coefficient: A proportionality constant between the friction torque and angular velocity, which must have a positive value.

Figure 6 and Figure 7 presents the relation between the angular velocity and the generated friction torque.

Figure 7 shows the simulated friction torque based on Equations (41) to (43)

In the above simulation, the following assumptions were made:

$T_{Peak}=100Nm,$ $\omega =100rad/sec,$ $f=5Nm/(rad/sec)$, $T_C=0.1\ast T_{Peak}Nm$, ${\omega }_{brk}=0.01rad/sec$, $T_{brk}=0.15\ast T_{Peak}Nm$

Friction between two interacting parts with relative motion is unavoidable. By considering the joint friction torques, the robot’s dynamics are modified as follows:

$${\tau }_{Joint}=\left[M\left(\theta \right)\ddot{\theta }+V\left(\theta ,\dot{\theta }\right)+G\left(\theta \right)+{\tau }_{friction}\right]$$

(44)

where,

$${\tau }_{friction}=F\left(\dot{\theta }\right)$$

(45)

Equation (38) can be written in the form of Equation (46)

$$\ddot{\theta }=M{\left(\theta \right)}^{-1}\left[\tau -V\left(\theta ,\dot{\theta }\right)-G\left(\theta \right)-F\left(\dot{\theta }\right)\right]$$

(46)

The next section (section 4) will explain the construction of the Realistic Model Reference Computed Torque Controller (RMRCTC).

4. Realistic Model Reference Computed Torque Control

The realistic model reference computed torque controller is an enhanced version of the standard model reference computed torque controller. One of the key limitations of the traditional model reference computed torque controller is that it requires the dynamic model of the plant to be executed twice during each sampling period, demanding significant computational power [42]. The computed torque control scheme, theoretically based on the plant’s inverse dynamics, is utilized within the model reference computed torque controller for real-time torque estimation. The following section will discuss the computed torque controller, the model reference computed torque controller, and the realistic model reference computed torque controller in detail.

4.1. Computed Torque Controller:

The computed torque controller is widely recognized as one of the most effective model-based control methods, leveraging inverse dynamics. It consists of two key components: a linearization mechanism and a control action. However, its primary limitation is the need for an accurate dynamic model of the system, making it more suitable for applications such as industrial robotics, where model parameters are well-defined. In contrast, rehabilitation exoskeletons must adapt to users with varying body sizes and shapes, making it impractical to accurately measure the mass and inertia of individual limbs. Since the controller’s stability and trajectory tracking performance rely heavily on precise modeling, this presents a significant challenge. Furthermore, the entire dynamic model must be recalculated during each sampling period, placing a heavy computational burden on the control system..

It performs two key functions: first, it linearizes the nonlinear system using feedback linearization, and second, it generates a control input that meets performance criteria [38]. Figure 8 illustrates the architecture of the CTC scheme.

In Figure 8, $M\left(\theta \right)\in {\mathbb{R}}^{n\times n}$ represents the robot mass matrix, $V\left(\theta ,\dot{\theta }\right)\in {\mathbb{R}}^{n\times 1}$ indicates the Coriolis and centrifugal forces, and $G\left(\theta \right)\in {\mathbb{R}}^{n\times 1}$ denotes the gravitational force. Figure 1 clearly shows that the CTC relies on the robot’s dynamic model for both linearizing the robot dynamics and generating the control input. Consequently, any discrepancies between the robot and its model can result in instability and uncontrolled behavior, making it essential to protect the system from such discrepancies.

4.2. Model Reference Computed Torque Controller

A promising alternative to the computed torque controller is the model reference computed torque controller. In this approach, the plant’s model is executed at every sampling time with the help of a computed torque controller for estimating the torque required to track the provided trajectories. The predicted torque is then applied to the physical system. The trajectory tracking error between the model and the physical plant is measured, and any discrepancies are corrected using an additional PID controller. Figure 9 shows the schematic diagram of a model reference computed torque controller. In essence, the model reference computed torque control scheme requires the robot’s dynamics to be executed twice in each sampling period (Robot’s dynamic model and M, G, V matrices). This significantly increases the computational load on the digital controlle. High speed multiple core digital computer is required to run the control algorithms. The control system consumes a significant amount of energy to run the control algorithms.

4.3. Realistic Model Reference Computed Torque Controller

In order to make the Model Reference Computed Torque Controller more efficient and realistic multiple changes were made in the Realistic Model Reference Computed Torque Controller.

Equation (36) shows that the total torque required to operate the robot is used for accelerating the robot’s links and payloads (M matrix), compensating for gravitational effects (G matrix), and counteracting Coriolis and centrifugal effects (V matrix). Often joints friction are considered as the disturbances. It has been observed that the magnitude of Coriolis and centrifugal forces is significantly smaller compared to the torque needed for acceleration and gravity compensation.

Figure 10. Comparison between the total torque, torque required for link acceleration (M), to overcome the gravity (G), Coriolis and centrifugal forces (V).

Figure 10. Comparison between the total torque, torque required for link acceleration (M), to overcome the gravity (G), Coriolis and centrifugal forces (V).

Calculating the Coriolis and centrifugal terms at every sampling period for the robot’s specific position and velocity demands more computational power than executing the mass and gravity matrices (M and G matrices).

To optimize the realistic model reference computed torque controller, we excluded the Coriolis and centrifugal terms from the dynamic model, treating them as disturbances. Additionally, we divided the control algorithm into two loops: slower loop for joint torque estimation and faster loop contains the PID controller for minimizing the tracking error between the reference trajectory and the plant output. The slower loop run at (100Hz), while the faster loop run at 1kHz. Running 2 loops at two different rates reduces the computational load and energy consumption for real time computing. Predicting torque at slower speed cause some inaccuracies. The PID controller was used to remove the errors and improve the robustness.The developed realistic model reference computed torque controller showed excellent trajectory tracking performance, maintaining small tracking errors even in the presence of disturbances. By excluding the Coriolis and centrifugal force terms $V\left(\theta ,\dot{\theta }\right)$, it achieves significantly higher energy efficiency compared to traditional computed torque control and model reference computed torque control methods. Figure 11 shows the schematic diagram of the realistic model reference computed toque controller.

Stability Analysis of the Proposed Realistic Model Reference Computed Torque Controller:

Ensuring controller stability is a very fundamental requirement of a control system. The following section will prove the combined system stability by computing the stability of the individual loops.

Loop1:

The dynamic model of the robot can be represented by Equation (47),

$$\left[M\left(\theta \right)\ddot{\theta }+V\left(\theta ,\dot{\theta }\right)+G\left(\theta \right)\right]={\tau }_{Joint}$$

(47)

The control torque for the model is calculated by using Equation (48),

$${\tau }_M=M\left(\theta \right)[{\ddot{\theta }}_d+K_v({\dot{\theta }}_d-\dot{\theta })+K_P({\theta }_d-\theta \left)\right]$$

(48)

By comparing Equations (47) and (48),

$$(\ddot{{\theta }_d}-\ddot{\theta })+K_V(\dot{{\theta }_d}-\dot{\theta })+K_P({\theta }_d-\theta )=0$$

(49)

The robot tracking error is defined as $E=({\theta }_d-\theta )$. Using Equation (49), the error dynamics can be expressed as:

$$\ddot{E}+K_V\dot{E}+K_{P}E=0$$

(50)

By applying the Laplace transform on Equation (50), the closed-loop characteristic polynomial can be written as,

$${\mathsf{\Delta }}_{\mathrm{C}}\left(S\right)=|S^{2}I+K_{V}S+K_P|$$

(51)

$$\mathrm{where},K_P=diag\left\{K_{Pi}\right\}$$

(52)

$$K_V=diag\left\{K_{Vi}\right\}$$

(53)

According to the Routh-Hurwitz stability criteria in Table 4, Equation (51) is asymptotically stable as long as the gain matrices $(K_P,K_V)$ are positive definite.

Loop 2:

Calculates the deviation of the plant output from the model $(e,\dot{e})$ and forces the plant to compensate for the tracking errors in a controlled way. Based on the developed RMRCTC scheme, the plant input torque $\left({\tau }_P\right)$ can be calculated as,

$${\tau }_P={\tau }_M+K_I\int \left({\theta }_M-{\theta }_P\right)+K_P\left({\theta }_M-{\theta }_P\right)+K_V\left(\dot{{\theta }_M}-\dot{{\theta }_P}\right)+{\tau }_f$$

(54)

In developing this controller, it is assumed that the model and the plant share the same dynamical structure, differing only in their parameters.

The tracking error $(E,\dot{E})$ between the model and the plant is defined as,

$$E=({\theta }_M-{\theta }_P)$$

(55)

$$\dot{E}=({\dot{\theta }}_M-{\dot{\theta }}_P)$$

(56)

$${\tau }_P={\tau }_M+K_V\dot{E}+K_I\int E+K_{P}E+{\tau }_f$$

(57)

$${\tau }_P\left(S\right)={\tau }_M\left(S\right)+{\displaystyle \frac{K_{I}E\left(S\right)}{S}}+K_{P}E\left(S\right)+K_{V}SE\left(S\right)+{\tau }_f\left(S\right)$$

(58)

$${\tau }_P\left(S\right)-{\tau }_M\left(S\right)-{\tau }_f\left(S\right)={\displaystyle \frac{\left({S^{2}K}_V+S{K}_P+K_I\right)E\left(S\right)}{S}}$$

(59)

Using Equation (59), the error dynamics can be expressed as:

$${\displaystyle \frac{E\left(S\right)}{{\tau }_P\left(S\right)-{\tau }_M\left(S\right)-{\tau }_f\left(S\right)}}={\displaystyle \frac{S}{S^2{K}_V+S{K}_P+K_I}}$$

(60)

Equation (60) represents the transfer function of the tracking error

Table 5 shows that the control system will be asymptotically stable if the gain matrices $K_P,K_V,K_I$ are positive definite. Here, $K_P,K_I,K_V$ are the $\left(7x7\right)$diagonal matrices.

$$K_P=diag\left\{K_P\right\}$$

(61)

$$K_I=diag\left\{K_I\right\}$$

(62)

$$K_V=diag\left\{K_V\right\}$$

(63)

An analysis of the stability of Loop1 and Loop 2, quickly leads to the conclusion that the proposed control system is asymptotically stable as long as the gain matrices are positive definite.

5. Simulation Results and Discussion

The trajectory tracking simulation was conducted in a MATLAB-Simulink^® environment. The mass and inertial properties of the human lower extremities used in the simulation are listed in Table 6. Controller gains for both the computed torque control (Loop 1) and the PID controller (Loop 2) are also provided in Table 6. The simulation was performed for both sequential and simultaneous joint movements. Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 display the results of the developed controller for the human lower extremity’s sequential and simultaneous joint movements.

The parameters used for modeling dynamic friction are listed below,

$$f=0.1Nm/(rad/sec),{\tau }_{coulomb}=0.1\ast {\tau }_{peak}Nm,{\tau }_{brake}=0.15\ast {\tau }_{peak}Nm,$$

$${\omega }_{stribeck}=\sqrt{2}{\omega }_{brake}rad/sec,{\omega }_{coulomb}=0.10\ast {\omega }_{brake}rad/sec$$

The peak torques were calculated based on the torque required to track the same trajectory using the dynamic model of the human lower extremity.

Figure 12, Figure 13, Figure 14 and Figure 15 illustrate the tracking performance, including trajectory tracking, tracking error, joint torque, and friction torque generated during sequential joint movement.

Figure 12 shows the input trajectory (Trajec. Joint #) supplied to the robot and the resultant output (Output Joint #) trajectories. The simulation result indicates the output closely follows the input trajectories with minimal error. Figure 13 provides a more detailed view of the tracking errors.

Figure 13 shows that the maximum tracking errors for sequential joint movements hip abduction-adduction, hip flexion-extension, hip internal-external rotation, knee flexion-extension, shank internal rotation, and ankle flexion-extension and pronation-supination were [0.38°, 0.68°, 0.12°, 0.29°, -0.19°, 0.17°, 0.16°] respectively.

Figure 14 shows the joint torque needed for the model (Orange color) and the physical robot joints (blue) during the tracking of sequential joint movements. It is observed that the plant requires more torque than the model. The torque required by the plant consists of two components: firstly, the torque necessary for the model for tracking the trajectory and secondly the additional torque required for compensating the differences between the plant input trajectory and the plant’s output.

Figure 15 presents the frictional torque generated in the joints during sequential joint movements. Equation (41) was used for estimating the dynamic friction torque generated in the joints due to relative motion between the mating surfaces. The simulation results indicate that joint 2 experienced the highest friction, followed by joints 1 and 4. In comparison, joints 3, 5, 6, and 7 produced significantly lower friction torque than joints 1, 2, and 4.

Figure 16, Figure 17, Figure 18 and Figure 19 present the tracking performance of the RMRCTC during simultaneous joint movements. Figure 16 shows the input (Trajec. Joint #) and output (Output Joint #) trajectories, demonstrating that the output trajectories closely followed the input trajectories.

Figure 17 illustrates the tracking errors that occurred during simultaneous joint movements, with all joints starting simultaneously and tracking the input trajectories. The maximum tracking errors recorded were [0.38°, 0.65°, 0.22°, 0.33°, 0.13°, 0.16°, 0.17°].

Figure 18 illustrates the joint torque required by both the model and the plant during simultaneous joint movements. Similar to sequential joint movements, the plant required more torque than the model for simultaneous movements.

The torque required by the plant consists of the torque needed by the model, along with additional torque to compensate for the relative errors between the input trajectory and the plant’s output.

Figure 19 shows the friction torque generated during simultaneous joint movements. Similar to the sequential joint movements joint 2 experienced the highest friction torque, followed by joints 1, 4, and the others.

All the PID controller gains used to make the plant behave like the model were determined through trial and error. Fine-tuning the controller gains may further reduce the tracking error.

Figure 19 shows the joint friction torque developed during the simultaneous joint movements.

6. Controller Performance Verification

Statistical analysis was conducted to assess the controller’s robustness and tracking performance under extreme conditions. Of the dynamic model’s various parameters, only two were variable: the subject’s weight and height. These two parameters were adjusted (as shown in Table 6), and different values were applied while collecting the resulting trajectory tracking errors.

Figure 20 illustrates the joint trajectory tracking errors for various weights and heights of the subject. The error histogram reveals that changes in the subject’s weight have negligible impact on the trajectory tracking errors, while variations in height do affect them. The median and standard deviation of the joint trajectory tracking errors were also calculated and maximum possible trajectory tracking error were tabulated in Table 7.

Based on the obtained mean/median and the standard deviation the maximum possible trajectory tracking errors were calculated. Syms like joint 2 and joint 5 have higher standard deviation that may lead to tracking errors up to 5°. It can also be concluded that based on the subject’s height and weight the controller’s gain can be tuned to get better accuracy.