1. Introduction
One of the main causes of lower limb disability is a stroke, which affects a person’s coordination, balance, and mobility [Reference Gorelick1]. Rehabilitation of lower limbs is currently done with manual techniques that involve repetitive lower limb movements along a natural gait trajectory under the constant supervision of physiotherapists. However, due to high demand and an insufficient number of specialists available, many patients do not receive the necessary training on time, leaving them untreated [Reference Chen, Qi, Guo and Yu2]. Over the past two decades, there has been an increase in the development of robotic devices designed for the lower limb training of stroke patients. This approach has gained significant attention as it addresses the limitations of manual therapy and assists patients in regaining their ability to walk [Reference Dollar and Herr3–Reference Zou, Zeng, Huang, Peng, Zhang and Cheng5]. Popular treadmill-based gait rehabilitation robots, such as LOKOMAT [Reference Colombo, Joerg, Fau - Schreier, Schreier, Fau - Dietz and Dietz6], Lower Extremity Powered Exoskeleton (LOPES) [Reference Veneman, Kruidhof, Hekman, Ekkelenkamp, Asseldonk and v. d. Kooij7], Active Leg Exoskeleton (ALEX) [Reference Banala, Kim, Agrawal and Scholz8], ReoAmbulator [Reference Fisher, Lucas and Thrasher9] and the others [Reference Hussain, Jamwal, Ghayesh and Xie10–Reference Copilusi, Ceccarelli and Carbone14] are equipped with several actuators that actively control the motions of lower limb joints to facilitate walking. Most of the assistive motion supplied by these multi-DOF (Degree-of-Freedom) lower limb rehabilitation robots occurs in the sagittal plane. Thus, essential training of the natural gait can still be achieved, avoiding redundant DOFs. This balance allows users to walk without feeling overly constrained by an excessively intricate robotic system. By employing designs that use fewer components, simplified control algorithms, and less complex assembly methods, it is possible to significantly lower the cost of gait rehabilitation robots without compromising the quality of gait training.
Recently, there has been increased attention on developing less complex designs for lower limb rehabilitation robots based on single DOF linkage mechanisms. These mechanisms usually based on four-bar [Reference Ji and Manna15, Reference Singh, Chaudhary and Singh16], five-bar [Reference Sun, Hu, Gao, Gao, Ma and Wang17], six-bar [Reference Tsuge, Plecnik and Michael McCarthy18–Reference Li, Fau - Yan, Yan, Fau - Zhao, Zhao, Fau - Ma, Ma, Fau - Li and Li20], eight-bar [Reference Shin, Deshpande and Sulzer21–Reference Wei, Zhang and Zhang24], and ten-bar [Reference Tsuge and McCarthy25] linkage systems, have demonstrated their capability to generate human-like gait motion trajectories. The majority of these linkage-based gait rehabilitation robots are in the development stage, with primarily theoretical concepts presented. All manufactured proof-of-concept prototypes are reported to be controlled using constant speed. The actuator of the linkage is rotated at a constant speed, moving the legs in a looped trajectory generated by the mechanism without considering the stiffness of the lower limbs. Such a control method may not be suitable for stroke survivors, as the level of disability varies among them. The movement intention, in terms of torque supplied by the patient, needs to be estimated to adjust the assistive robotic torque accordingly, allowing the patient to engage more fully in a rehabilitation process [Reference Cai, Fong, Otoshi, Liang, Burdick, Roy and Edgerton26]. Active human participation during the gait training can significantly enhance the effectiveness and speed of recovery [Reference Lotze, Braun, Fau - Birbaumer, Birbaumer, Fau - Anders, Anders, Fau - Cohen and Cohen27]. For this purpose, terms such as Assist-as-Needed (AAN), human-cooperative, compliant, and interactive control schemes have been extensively used in the literature [Reference Riener, Lünenburger, Jezernik, Anderschitz, Colombo and Dietz28–Reference Viteckova, Kutilek, de Boisboissel, Krupicka, Galajdova, Kauler, Lhotska and Szabo32]. A recent review covering the patient-cooperative control strategies implemented for multi-DOF lower limb exoskeletons is presented in ref. [Reference Masengo, Zhang, Dong, Alhassan, Hamza and Mudaheranwa33].
Human-cooperative control strategy, which is sometimes referred to as AAN can be realized by controlling the stiffness of the robot depending on the external forces/torques during training. By adjusting the robot’s stiffness responding to varying external torques due to the patient’s disability conditions, customized assistance can be achieved. Stiffness control is part of a general well-known impedance control. The concept of impedance control, pioneered by Hogan [Reference Hogan34], allows for the precise adjustment of robot impedance, enabling the robot to respond to external forces/torques during training. So far, impedance control has been implemented among various multi-DOF gait rehabilitation robots such as Lokomat [Reference Riener, Lünenburger, Jezernik, Anderschitz, Colombo and Dietz28]. However, the heavy weight and high friction, limit the robot to operate properly in impedance-controlled mode. This control scheme has been presented in pneumatically driven gait rehabilitation robots [Reference Hussain, Xie and Jamwal35, Reference Jamwal, Hussain, Ghayesh and Rogozina36], however, such types of actuators are not precise and the level of assistive force may be incorrectly supplied. For patients with mild impairments, a stiffness control between the robot and the human accommodates natural movement patterns, while for individuals unable to move their legs, a stiffer robot offers the necessary stability and support. Establishing an active interaction between the human and the robot is an important aspect of the gait rehabilitation process as it can expedite the recovery process.
This research focuses on the development of stiffness control for the Stephenson III six-bar linkage-based gait rehabilitation robot for providing gait training based on the principles of AAN. We simplify the impedance control concept to stiffness control between the human and the robot as the developed robot doesn’t require quick velocity changes and the friction is low. Readers are referred to [Reference Kapsalyamov, Hussain, Brown, Goecke, Hayat and Jamwal37, Reference Kapsalyamov, Hussain, Brown, Goecke and Jamwal38] for more information on the developed robot, its dimensional synthesis, and its dynamic characteristics. The constructed gait rehabilitation robot has only one actuator to drive the system. The motor is located at the hip region and is able to provide naturalistic lower limb trajectories for knee and ankle joints relative to the hip joint in sagittal plane as shown in Figure 1(c). It is important to assure that the trajectories generated by the mechanism provide natural walking angles to avoid joints following unintended motion patterns against body’s natural style [Reference Mihcin39]. The optimization of the linkage mechanism to meet the naturalistic walking requirements is presented in our previous work [Reference Kapsalyamov, Hussain, Brown, Goecke, Hayat and Jamwal37]. The device has a telescopic feature which allows to extend the link to fit people with different leg lengths. The device can assist humans weighing up to 130 kg, which is an integral part of the robotic gait rehabilitation system. The torque corresponding to one that is required to move the dead weight of the lower limbs of the person who is completely paralyzed is weight normalized.
Figure 1. (a) The gait rehabilitation system and its major components with passive DOFs labeled, (b) Loadcells’ placement behind the braces for capturing the human/robot interaction force, (c) Trajectories generated by the mechanism for knee and ankle joints relative to the hip joint in the sagittal plane.
This paper is organized as follows: Section II presents the dynamical formulation for the mechanism, the stiffness control scheme, and human intent torque estimation using computational neural model. Section III presents the experimental protocol and the results obtained from eight healthy participants walking in the rehabilitation robot controlled by the implemented control scheme. Finally, the discussion and conclusions from the present research are provided in Section IV.
The contributions of this work include the development of a customized stiffness control strategy for the Stephenson III six-bar linkage robot, which is actuated by a single motor to provide customized training. To the authors’ best knowledge, no prior work has studied or implemented a human-cooperative control strategy on single DOF mechanisms designed for gait rehabilitation purposes. Thus, the novelty of this work lies in demonstrating the applicability of stiffness control techniques on such types of mechanisms. This work will also serve as a foundation for building human-cooperative control strategies for single DOF linkage-based rehabilitation robots.
2. Methodology
2.1. Stephenson III six-bar mechanism
The Stephenson III six-bar linkage mechanism developed is a bilateral system for both right and left leg as shown in Figure 1(a). It can provide lower limb motions required for walking in a sagittal plane and be operating previously using constant speed control [Reference Kapsalyamov, Hussain, Brown, Goecke and Jamwal38]. To further implement human-cooperative control, it is necessary to monitor the interaction forces between the robot and the human. Therefore, the prototype was appended with additional tension/compression load cells added at the thigh and shank regions in the prototype as indicated in Figure 1(b).
2.2. Safety
To ensure the safety of individuals engaged in experimental activities, the robotic system is equipped with essential safety provisions. Emergency stop buttons, integrated into an independent circuit, are incorporated within the system, accessible to both the researcher and the participant during the experiment. These buttons serve the purpose of promptly deactivating the power source in the event of an emergency or when the subject experiences discomfort. While the developed mechanism is designed to follow predetermined trajectories, supplementary mechanical stops are installed as an additional safety measure at the rocker links to avoid going beyond the range of motion. These mechanical stops serve to restrict the robot from surpassing physiological ranges of motion, particularly in unforeseen circumstances. They can withstand the peak torque delivered by the motor, thereby safeguarding against potential hazards. To ensure the safety of paralyzed individuals, especially those with weakened bone structures due to conditions such as osteoporosis, our rehabilitation system is meticulously designed to prevent any risk of bone fractures. The motor does not directly apply force to the bones; instead, it uses intermediary links to move the legs. This method significantly reduces the risk of injury as the motor’s torque is translated through these links rather than being applied directly to the bones. For instance, the motor rotates a crank, which then moves other links by pulling or pushing them, thereby providing a naturalistic movement trajectory that aligns with the human body’s natural angles. This design consideration minimizes the risk of unnatural forces that could lead to fractures. Furthermore, the system calculates the overall torque at the hip joint by considering the leg’s deadweight and applies the torque through the joints, ensuring it remains within safe limits. The peak torque supplied by our motor is 120 Nm, which is well below the torque required to fracture even weakened bones. It is reported that the compressive strength of a healthy adult femur translates to a fracture torque of approximately 30,750 Nm [Reference Morgan, Unnikrisnan and Hussein40]. This significant difference provides a substantial safety margin, ensuring that our system operates safely even for individuals with severely weakened bones. By maintaining this safety margin, we effectively prevent any risk of fractures, providing a secure and effective rehabilitation solution. Precautions are also taken into consideration for electricity and wiring of the machine to assure the device does not cause any electrocution of the patients during training. Given the importance of meeting regulatory standards for medical devices, this research acknowledges the necessity of adhering to ISO 13,485 regulations intended for medical devices and adheres to the safety requirements set for such devices [Reference Mihcin, Strehlow, Demedts, Schwenke, Levy and Melzer41, Reference Mihcin, Karakitsios, Le, Strehlow, Demedts, Schwenke, Haase, Preusser and Melzer42].
2.3. Dynamic formulation of the mechanism
The overall motion of the mechanism can be described using the Lagrange modeling methodology [Reference Goldstein, Poole and Safko43]. The notations and configurations defined for solving the dynamics of the mechanism are shown in Figure 2. The lengths between the joints are defined to be
$L_{i}$
. The masses of links
$\overline{GH}, \overline{HEK}$
(triangular platform
$p_{1}$
),
$\overline{BC}, \overline{CDE}$
(triangular platform
$p_{2}$
) and
$\overline{AD}$
are defined as
$m_{1}, mp_{1}, {m_{5}}, mp_{2}$
and
$m_{9}$
, respectively, and
$I_{i}$
is the mass moment of inertia for the corresponding link.
$V_{i}$
is the velocity at the mass center of the corresponding link. The distance from the joint to the center of mass of the ternary link is defined by
$R_{i}$
. The derivations of the assumed angles that result to the input angle motion
$\unicode{x1D6F3}$
at crank
$\overline{AD}$
are provided at the end of the manuscript in Appendix section. The whole system’s Lagrangian is defined as the difference between the total kinetic energy and the total potential energy.
Figure 2. Scheme of stephenson III six-bar linkage for motion evaluation.
\begin{equation} \mathcal{L}=K-P \end{equation}
The definition of kinetic energy
$K$
of the mechanism can be expressed as follows:
\begin{align} K=&\frac{1}{2}\left[\left(m_{1}V_{1}^{2}+I_{1}\dot{\unicode{x1D6F3}^{2}}_{1}\right)+\left(mp_{1}Vp_{1}^{2}+Ip_{1}\dot{\unicode{x1D6F3}^{2}}_{3}\right)+\left(m_{5}V_{5}^{2}+I_{5}\dot{\unicode{x1D6F3}^{2}}_{5}\right)+\left(mp_{2}Vp_{2}^{2}+Ip_{2}\dot{\unicode{x1D6F3}^{2}}_{6}\right)\right.\\&\qquad\qquad +\left.\left(m_{9}V_{9}^{2}+I_{9}\dot{\unicode{x1D6F3}^{2}}\right)\right] \nonumber\end{align}
where the substitution variables are as follows:
\begin{gather*} \unicode{x1D6F3}_{1}=\alpha _{2}+\beta _{2}; \\\unicode{x1D6F3}_{5}=\beta -\gamma;\\\unicode{x1D6F3}_{3}=\alpha _{3}+\beta _{3};\\\unicode{x1D6F3}_{6}=\alpha _{1}+\beta _{1};\\V_{1}^{2}=\frac{L_{1}^{2}}{4}\dot{\unicode{x1D6F3}^{2}}_{1};\\V_{5}^{2}=\frac{L_{5}^{2}}{4}\dot{\unicode{x1D6F3}^{2}}_{5};\\V_{9}^{2}=\frac{L_{9}^{2}}{4}\dot{\unicode{x1D6F3}^{2}};\\Vp_{1}^{2}=L_{1}^{2}\dot{\unicode{x1D6F3}^{2}}_{1}+R_{2}^{2}\dot{\gamma _{2}^{2}}+2L_{1}R_{2}\dot{\unicode{x1D6F3}_{1}}\dot{\gamma }_{2}\cos (\unicode{x1D6F3}_{1}-\gamma _{2});\\Vp_{2}^{2}=L_{9}^{2}\dot{\unicode{x1D6F3}^{2}}+R_{6}^{2}\dot{\gamma ^{2}}_{6}+2L_{9}R_{6}\dot{\unicode{x1D6F3}}\dot{\gamma }_{6}\cos (\unicode{x1D6F3}-\gamma _{6});\end{gather*}
Then, the potential energy
$P$
is determined as:
\begin{align} P=&\; m_{1}g\frac{L_{1}}{2}\sin \unicode{x1D6F3}_{1}+mp_{1}g\left(L_{1}\sin \unicode{x1D6F3}_{1}+R_{2}\sin \gamma _{2}\right)+mp_{2}g\left(L_{9}\sin \unicode{x1D6F3}+R_{6}\sin \unicode{x1D6F3}_{6}\right)\\&+m_{5}g\frac{L_{5}}{2}\sin \unicode{x1D6F3}_{5}+m_{9}g\frac{L_{9}}{2}\sin \unicode{x1D6F3} \nonumber\end{align}
where
$g$
is the gravitational acceleration.
After reconfiguring Equation (3), the final expression for the Lagrangian of the mechanism has the following form:
\begin{align}\mathcal{L}= &\; J_{1}\dot{\unicode{x1D6F3}^{2}}_{1}+J_{3}\dot{\unicode{x1D6F3}^{2}}_{3}+J_{5}\dot{\unicode{x1D6F3}^{2}}_{5}+J_{6}\dot{\unicode{x1D6F3}^{2}}_{6}+J_{0}\dot{\unicode{x1D6F3}^{2}}+M_{1}\dot{\gamma _{2}^{2}}+M_{2}\dot{\gamma ^{2}}_{6}+P_{1}C_{1}\left(\theta _{1},\gamma _{2}\right)\dot{\unicode{x1D6F3}_{1}}\dot{\gamma }_{2}\\& +P_{6}C_{6}\left(\theta, \gamma _{6}\right)\dot{\unicode{x1D6F3}}\dot{\gamma }_{6}+G\left(\unicode{x1D6F3}_{1},\unicode{x1D6F3}_{5},\unicode{x1D6F3}_{6},\unicode{x1D6F3},\gamma _{2}\right)\nonumber\end{align}
The definitions of variables used in the above equation
$J_{0}$
,
$J_{1}$
,
$J_{3}$
,
$J_{5}$
,
$J_{6}$
,
$M_{1}$
,
$M_{2}$
,
$P_{1}$
,
$C_{1}(\theta _{1},\gamma _{2})$
,
$P_{6}$
,
$C_{6}(\theta, \gamma _{6})$
and
$G(\unicode{x1D6F3}_{1},\unicode{x1D6F3}_{5},\unicode{x1D6F3}_{6},\unicode{x1D6F3},\gamma _{2})$
are presented in the Appendix. By introducing function
$S$
which is the derivative of the assumed angles
$\dot{\gamma }_{2},\dot{\gamma }_{6},\dot{\unicode{x1D6F3}_{1}},\dot{\unicode{x1D6F3}_{3}},\dot{\unicode{x1D6F3}_{5}},\dot{\unicode{x1D6F3}}_{6}$
with respect to
$\dot{\unicode{x1D6F3}}$
we could replace these velocity terms and present Lagrangian in the terms of
$\theta$
:
\begin{align} \mathcal{L}=&\left(J_{0}+J_{1}S_{1}^{2}\left(\theta, \alpha _{2}, \beta _{2}\right)+J_{3}S_{3}^{2}\left(\theta, \alpha _{3}, \beta _{3}\right)+J_{5}S_{5}^{2}\left(\theta, \beta, \gamma \right)+J_{6}S_{6}^{2}\left(\theta, \alpha _{1}, \beta _{1}\right)+M_{1}S_{2}^{2}\left(\theta, \gamma \right)\right.\\&+\left. M_{2}S_{7}^{2}\left(\theta, \gamma \right)+P_{1}C_{1}\left(\theta _{1},\gamma _{2}\right)S_{8}\left(\theta, \alpha _{2}, \beta _{2},\alpha _{3}\right)+P_{6}C_{6}\left(\theta, \gamma _{6}\right)S_{9}\left(\theta, \alpha _{1}\right)\right)\dot{\unicode{x1D6F3}^{2}}\nonumber\\&+G\left(\unicode{x1D6F3}_{1},\unicode{x1D6F3}_{5},\unicode{x1D6F3}_{6},\unicode{x1D6F3},\gamma _{2}\right) \nonumber\end{align}
The expression for the Lagrangian contains parameters that are also functions of
$\theta$
. The motor rotates the crank (input link which the motor drives continuously) by input angle
$\theta$
. Next, we need to identify the equation of motion by:
\begin{equation} \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\unicode{x1D6F3}}}\right)-\frac{\partial \mathcal{L}}{\partial \theta }=T_{rob} \end{equation}
$T_{rob}$
is the torque supplied by the robot,
$\dot{\unicode{x1D6F3}}$
is the angular velocity of the crank
$\overline{AD}$
. Since the damping in the mechanism is relatively small, and bearings are used between the links to minimize the friction between them, dissipative energy can be ignored and is not considered in this work. To finalize the equation of motion, first, let’s determine the following:
\begin{align} \frac{\partial \mathcal{L}}{\partial \dot{\unicode{x1D6F3}}}=&\; 2\left(J_{0}+J_{1}S_{1}^{2}\left(\theta, \alpha _{2}, \beta _{2}\right)+J_{3}S_{3}^{2}\left(\theta, \alpha _{3}, \beta _{3}\right)+J_{5}S_{5}^{2}\left(\theta, \beta, \gamma \right)+J_{6}S_{6}^{2}\left(\theta, \alpha _{1}, \beta _{1}\right)+M_{1}S_{2}^{2}\left(\theta, \gamma \right)\right.\\&+\left. M_{2}S_{7}^{2}\left(\theta, \gamma \right)+P_{1}C_{1}\left(\theta _{1},\gamma _{2}\right)S_{8}\left(\theta, \alpha _{2}, \beta _{2},\alpha _{3}\right)+P_{6}C_{6}\left(\theta, \gamma _{6}\right)S_{9}\left(\theta, \alpha _{1}\right)\right) \nonumber\end{align}
Next, the derivations for
$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\unicode{x1D6F3}}}\right)$
and
$\frac{\partial \mathcal{L}}{\partial \theta }$
together with the definitions of assumed angles can be found in the Appendix. The overall dynamic equation of motion for the system can be described using the following expression:
\begin{align}&2\left(J_{0}+J_{1}S_{1}^{2}+J_{3}S_{3}^{2}+J_{5}S_{5}^{2}+J_{6}S_{6}^{2}+M_{1}S_{2}^{2}+M_{2}S_{7}^{2}+P_{1}C_{1}S_{8}+P_{6}C_{6}S_{9}\right)\ddot{\theta }\nonumber\\&+\left[2J_{1}S_{1}\left(\frac{\partial S_{1}}{\partial \theta }+\frac{\partial S_{1}}{\partial \alpha _{2}}\frac{\partial \alpha _{2}}{\partial \theta }+\frac{\partial S_{1}}{\partial \beta _{2}}\frac{\partial \beta _{2}}{\partial \theta }\right)+2J_{3}S_{3}\left(\frac{\partial S_{3}}{\partial \theta }+\frac{\partial S_{3}}{\partial \alpha _{3}}\frac{\partial \alpha _{3}}{\partial \theta }+\frac{\partial S_{3}}{\partial \beta _{3}}\frac{\partial \beta _{3}}{\partial \theta }\right)\right.\nonumber\\& +2J_{5}S_{5}\!\left(\frac{\partial S_{5}}{\partial \theta }+\frac{\partial S_{5}}{\partial \beta }\frac{\partial \beta }{\partial \theta }+\frac{\partial S_{5}}{\partial \gamma }\frac{\partial \gamma }{\partial \theta }\right)\!+2J_{6}S_{6}\left(\frac{\partial S_{6}}{\partial \theta }+\frac{\partial S_{6}}{\partial \alpha _{1}}\frac{\partial \alpha _{1}}{\partial \theta }+\frac{\partial S_{6}}{\partial \beta _{1}}\frac{\partial \beta _{1}}{\partial \theta }\!\right)\!+2M_{1}S_{2}\left(\frac{\partial S_{2}}{\partial \theta }+\frac{\partial S_{2}}{\partial \gamma }\frac{\partial \gamma }{\partial \theta }\!\right)\nonumber\\& +2M_{2}S_{7}\left(\frac{\partial S_{7}}{\partial \theta }+\frac{\partial S_{7}}{\partial \gamma }\frac{\partial \gamma }{\partial \theta }\right)+P_{1}\left(C_{1}\left(\frac{\partial S_{8}}{\partial \theta }+\frac{\partial S_{8}}{\partial \alpha _{2}}\frac{\partial \alpha _{2}}{\partial \theta }+\frac{\partial S_{8}}{\partial \beta _{2}}\frac{\partial \beta _{2}}{\partial \theta }+\frac{\partial S_{8}}{\partial \alpha _{3}}\frac{\partial \alpha _{3}}{\partial \theta }\right)\right.\nonumber\\& + \left. S_{8}\left(\frac{\partial C_{1}}{\partial \theta _{1}}\frac{\partial \theta _{1}}{\partial \theta }+\frac{\partial C_{1}}{\partial \gamma _{2}}\frac{\partial \gamma _{2}}{\partial \theta }\right)\right) \left.+P_{6}\left(C_{6}\left(\frac{\partial S_{9}}{\partial \theta }+\frac{\partial S_{9}}{\partial \alpha _{1}}\frac{\partial \alpha _{1}}{\partial \theta }\right)+S_{9}\left(\frac{\partial C_{6}}{\partial \theta }+\frac{\partial C_{6}}{\partial \gamma _{6}}\frac{\partial \gamma _{6}}{\partial \theta }\right)\right)\right]\dot{\theta }^{2}\nonumber\\&-\frac{\partial G}{\partial \theta }+\frac{\partial G}{\partial \unicode{x1D6F3}_{1}}\frac{\partial \theta _{1}}{\partial \theta }+\frac{\partial G}{\partial \unicode{x1D6F3}_{5}}\frac{\partial \theta _{5}}{\partial \theta }+\frac{\partial G}{\partial \unicode{x1D6F3}_{6}}\frac{\partial \theta _{6}}{\partial \theta }+\frac{\partial G}{\partial \gamma _{2}}\frac{\partial \gamma _{2}}{\partial \theta }=T_{rob}\end{align}
2.4. Stiffness control architecture
Depending on the level of lesion, the disability in the lower limbs among stroke patients may vary substantially. Stiffness control can adjust and provide the necessary assistive motion to the subject based on his/her needs. The assistance from the robot is modified by considering the human joint forces supplied by the subject. When the subject does not move, the robot increases the assistance to guide the subject along the physiological gait trajectories. On the other hand, if the subject makes an effort to walk, the robot reduces stiffness. This approach is effective in patients requiring gait training, as it encourages active participation from them, which facilitates the recovery process of impaired gait capabilities.
The architecture of the implemented stiffness controller is illustrated in Figure 3. The angular positions of the robot’s joints are measured using the rotary encoders. The load cells are installed on the thigh and shank regions to track the interaction force between the human and the robot during walking. The human-robot interaction forces occurring at the shank
$F_{sh}$
and thigh segments
$F_{th}$
, measured from loadcells, crank position
$\unicode{x1D6F3}$
and its angular velocity
$\dot{\unicode{x1D6F3}}$
serve as essential inputs for the overall control scheme. The dynamic model from (6) determines the robot’s torque based on the inputs such as crank angle
$\unicode{x1D6F3}$
and its angular velocity
$\dot{\unicode{x1D6F3}}$
. The computational neural model, denoted as
$h(\cdot ),$
accepts measured four inputs from the human-robot system and produces a single output of predicted human torque exerted on the system as follows:
Figure 3. Overall stiffness control scheme for the gait rehabilitation. The robot outputs the crank position, crank angular velocity, and interaction forces at shank and thigh regions. The robot dynamics block calculates the robot’s torque, while a computational neural model block estimates the human’s torque. This scheme enables dynamic adjustment of the robot’s assistance based on the human’s input.
\begin{equation} h\left(F_{sh}, F_{th}, \unicode{x1D6F3}, \dot{\unicode{x1D6F3}}\right)\approx \hat{\tau }_{h} \end{equation}
Where
$\hat{\tau }_{h}$
is the predicted human torque (explained in the following section) obtained from the computational neural model. Further, the robot’s torque
$T_{rob}$
and the predicted human torque
$\hat{\tau }_{h}$
are added together producing the total human-robot torque.
2.5. Computational neural model for human torque estimation
This section describes a computational neural model-based system identification approach to estimate the combined torque supplied by the human during the training process on the robot. Solving analytically for the human torque, based on the position and interaction forces for this mechanism, is a challenging and time-consuming process due to the non-linearity and incomplete knowledge of the biomechanical model. Instead, existing machine learning techniques can be utilized to solve this task. These techniques do not require strict knowledge of the system’s model, can manage uncertainties, and can establish additional criteria to evaluate the system’s performance [Reference Nguyen and Cheah44]. Therefore, a computational neural model was used in this work, since it is a powerful tool for nonlinear function approximation. The input layer is represented by the vector
$x=[F_{sh}, F_{th}, \unicode{x1D6F3}, \dot{\unicode{x1D6F3}}]$
. The network outputs a single value, the predicted human torque. The model is expected to converge to the biomechanically applied torque values, despite the absence of ground truth for the human torque values. This expectation is based on the use of measured interaction forces occurring at the shank and thigh regions as inputs for the model. Various hyperparameters, including the depth of the network, number of nodes/neurons per layer, learning rate, momentum, etc., were selected heuristically through a grid search approach [Reference Jiang and Xu45]. ReLU is used as the activation function to avoid the vanishing gradients problem and save computational time. The network is trained using the Adam (Adaptive Moment Estimation) gradient-based optimizer [Reference Kingma and Ba46] with a learning rate of 0.001. The architecture of the model used can be presented as follows:
\begin{equation} h\left(x\right)=W_{2}\cdot ReLU\left(W_{1}\cdot x+b_{1}\right)+b_{2} \end{equation}
where ReLU is the Rectified Linear Unit activation function [Reference Nair and Hinton47], represented by
\begin{align} ReLU\left(x\right)=\begin{cases} x,\; if\ x\geq 0\\ \textit{negative}\_ \textit{slope}\times x, \textit{otherwise} \end{cases} \end{align}
$W_{1}\in \mathbb{R}^{32\times 4}, W_{2}\in \mathbb{R}^{1\times 32}$
are weight matrices, and
$b_{1}\in \mathbb{R}^{32\times 1}$
,
$b_{2}\in \mathbb{R}$
are bias vectors.
$\boldsymbol{\psi }=(W_{1},W_{2}, b_{1},b_{2})$
are learnable parameters of the neural model and are optimized using an objective function. The objective function is to minimize the MSE between the predicted and approximated human torque values. We compute approximated human intention torque values by subtracting the measured torque from the motor driving the human/robot system and the torque required to drive empty robot without a human subject. We have also added
$L_{2}$
regularized to balance variance-bias tradeoff. The objective loss function is defined as follows:
\begin{equation} \mathcal{L}_{MSE}=\frac{1}{N}\sum _{j=1}^{N}\sum _{i=1}^{T}\left({\hat{\tau }_{{h_{\left(i,j\right)}}}}-\left(\tau _{moto{r_{\left(i,j\right)}}}-\tau _{{r_{\left(i,j\right)}}}\right)\right)^{2}+\lambda \boldsymbol{\psi } \end{equation}
Here
$N$
is the number of samples,
$T$
is the number of instances recorded in a single gait cycle,
$\hat{\tau }_{h}$
is the predicted value of the human torque,
$\tau _{\textit{motor}}$
is the actual torque provided by the motor to drive the human/robot system and
$\tau _{r}$
is the torque supplied to drive the robot without a human subject. In (12), the subscripts
$(i,j)$
correspond to two different aspects of the data. The ‘
$i$
’ index refers to the specific sample within the batch and can take on any value within the total number of samples, N. On the other hand, the
$'j'$
index refers to the time-specific instance within a single gait cycle torque trajectory, with ‘
$j$
’ being any value within the total instances,
$T$
.
$\lambda$
is a regularization parameter used to reduce the risk of overfitting.
The value of
$\lambda$
was determined empirically through a series of preliminary experiments, evaluating the model’s performance on a validation set for different values of
$\lambda$
. The selection was based on achieving an optimal balance between bias and variance, thereby minimizing the validation error. In an overfit model, the coefficients are generally inflated. Regularization adds penalties to the parameters and avoids them weighing heavily to avoid overfitting. Each crank angle value in the mechanism corresponds to a specific point in the gait cycle, allowing for a direct comparison of torque values with and without a human subject at the same gait cycle instance.
Before feeding the data into the model, it was preprocessed to remove outliers and scaled accordingly, ensuring that all features were on the same scale. The dataset consists of 330,397 points, recorded at a sampling rate of 100 Hz from the loadcells and encoders. According to the work done by Goyal & Bengio, the inductive bias introduced with prior knowledge helps to learn better in many tasks [Reference Goyal and Bengio48]. The processed data was then split into a training set (60%), validation set (20%), and testing set (20%) to assess the model’s performance. It was split on shuffled samples to strengthen the model’s generalizability on seen subjects’ different trials. The dataset is obtained from eight participants (7m, 1f) who were asked to stay passive while wearing the robot and allow it to move their legs in a predefined motion and walk concurrently while wearing the robot. This was to ensure capturing both end extremities in the interaction forces’ values.
3. Experimental evaluation
3.1. Experimental protocol
Before the start of the experiments, corresponding ethics approval to experiment on eight neurologically intact subjects (7 m, 1f) (age 22–45 years) was granted by the Human Research Ethics Committee at the University of Canberra. The subjects reported no history of having neurologic disorders. During the experiments, participants’ legs were fixed to the robot legs with the braces located in the thigh and shank segments using Velcro straps. Subjects’ feet were attached to the footplate. The body weight support system was used to balance the weight of human subjects against gravity, thus allowing them to keep their legs relaxed during the first stage of the experiment. During this stage, the robot moved the subjects’ legs on a predefined trajectory at a constant speed. In passive mode, participants were asked to stay relaxed without moving their legs and allowing the robot to direct the lower limbs along the predefined trajectory. During the active mode, participants were asked to perform the usual walking while wearing the robot. The walking speed for the participants was set by the speed of the treadmill. The study aimed to determine if significant differences existed between the interaction forces (at the robot and human interface) observed during passive and active modes. We utilized the Wilcoxon Rank-Sum test [Reference Wilcoxon49], a non-parametric method, to perform this analysis as it offers robustness against data not adhering to a normal distribution. Wilcoxon Rank-Sum test was used for both thigh and shank interaction forces and the torque supplied by the robot during active and passive
3.2. Results
Initial results from the stiffness control scheme developed for the proposed mechanism are presented in this section. The model was trained using the Pytorch Lightning framework, which was further converted to the Open Neural Network Exchange (ONNX) format and then imported to Matlab R2022B. It took the model less than 1 hour to be trained on 2nd Gen Intel(R) Core(TM) i7-1265U. The training and testing loss results are presented in Figure 4. The training Mean Squared Error (MSE) is a value of approximately 0.0835 Nm2. The small value for MSE indicates a high level of accuracy in predicting torque values based on the given input features, which include joint angular positions, angular velocity, and interaction forces. The validation set was introduced to avoid common problems in machine learning such as overfitting. The validation loss was monitored throughout the training process to ensure that the model was learning generalizable patterns instead of memorizing the training data.
Figure 4. Performance of the model: training loss, validation loss, and testing results.
The smooth convergence to a steady state, as depicted in Figure 4, validates the model selection, demonstrating that the optimized parameters are suitable for the proposed architecture. Regularization, introduced along with the cross-validation technique using validation sets, helped to avoid overfitting. Parameters of the network are fine-tuned using the validation set, the validation set is evaluated at every epoch. Once the learning curve stabilizes, training is halted, and the model is then evaluated on the test dataset. The test was done once at the last stage to verify it is the same as validation loss. These findings demonstrate the model’s ability to generalize well to unseen data and make accurate human torque predictions for robotic applications. The model’s ability to generate accurate torque predictions with a low error suggests that it has the potential to contribute to stable robot control.
The change of assistive torque supplied by the robot during passive and active modes over the gait cycles together with the torque to drive the robot without human subject are indicated in Figure 5. The torque supplied during active mode is the torque required to keep the robot moving while being aligned to the walking of the person. The values for the torques are the average across eight healthy participants. The crank of the robot was rotating at around 3 rad/s, driving a stride frequency of 0.47 Hz, and the walking speed of the treadmill was in the range of 0.2 – 0.4 m/s. The transition of the change in torque supplied by the robot in response to switching from active to passive and from passive to active modes is shown in Figure 6.
Figure 5. The adjusted torque from robot during the passive and active phases averaged across eight subjects with±1 st.dev. (shaded).
Further, the interaction forces that occurred during active and passive phases are plotted in Figure 7. The interaction forces occurring at the shank and thigh regions of the leg during passive phase are much larger than the interaction forces happening during the active phase. Each spike in the graph representing the crank angular position corresponds to a single gait cycle. As the torque adjusts the supplied torque based on the torque of human subject the human-cooperative control strategy is achieved.
Table 1. Maximum absolute values parameters. Standard deviations ± are presented for within-subject variability.
Figure 6. The overall torque supplied by the robot during transition from active to passive and passive to active modes.
Figure 7. The human/robot interaction forces occurring during the passive and active modes.
4. Discussion and conclusion
In this work, a customized stiffness control was designed for a Stephenson III six-bar linkage-based gait rehabilitation robot. The experimental results demonstrated the feasibility of the proposed control scheme and its applicability to such types of mechanisms. From the experimental results, significant differences were observed in the data from both the human-robot interaction forces and the torques supplied by the robot during the active and passive stages. The medians of the interaction forces and torques differed between active and passive modes, with both measures yielding a p-value of less than 0.01. This finding signifies a statistically significant difference in both interaction forces and robot-supplied torques between the two modes. The maximum absolute values recorded for the interaction forces at shank and thigh regions and
$T_{rob}$
provided during active and passive stages with standard deviation are presented in Table I.
The peak value for the robot’s torques is close between the active and passive phases, which indicates the instantaneous torque recorded at the starting point when the robot needs to start the motion from the resting stage. Further, the standard deviation in the active phase is much smaller than in the passive phase.
The boxplots that provide a graphical representation of the shank and thigh parts’ interaction forces observed during the passive and active modes are drawn in Figure 8. For the shank interaction forces, the boxplot reveals a distinct difference between the active and passive modes. It can be noticed that the median force in active mode is lower when compared to the passive mode. The interquartile range in active mode appears to be smaller, which means there is a smaller spread of the data. Since there are few outliers, most of the data are close to the median. A similar trend is observed in the boxplot for the thigh interaction forces. However, there is a large difference in median forces between the active and passive modes. The active mode has a smaller spread of data, as evidenced by the height of the box. The outliers in this case are few, suggesting that the data is more evenly distributed. This visual representation supports our statistical findings and highlights the significant differences in forces during the two different modes.
Figure 8. Distribution and comparison of human/robot interaction forces occurring at thigh and shank regions during active and passive phases.
A body weight support system was used to unload the subjects’ weight to help them relax their legs. This was done to simulate as closely as possible the impaired lower limbs of stroke survivors. We quantified the differences in the mechanical torque generated in these active/passive modes taken during multiple gait cycles. The results demonstrated that the robot could vary the torque supplied at the crank depending on the torque supplied by the human joints during. Providing assistive torque during specific phases of gait can have several clinical implications. For example, providing assistive torque during the stance phase can help stabilize the patient and reduce the risk of falls. On the other hand, providing assistive torque during the swing phase can help facilitate a more natural and efficient gait pattern. The key difference between the stiffness control implemented and a standard trajectory controller lies in the ability of the system to adapt to the user’s effort dynamically. While any trajectory controller will indeed reduce the robot’s torque when the human effort is helping, a stiffness controller will also adapt to the user’s effort in real time across the entire gait cycle. This means that if the user is unable to exert sufficient effort, the robot will supply more torque to keep the legs along the naturalistic lower limb trajectories. Conversely, if the user exerts more effort, the robot will reduce its assistance accordingly. This dynamic adaptation to the user’s effort helps to promote active participation from the user, which is crucial for effective rehabilitation.
One of the limitations of this work is that implemented computational neural model estimates the combined torque exerted by the lower limb on the robot, which causes challenges in identifying torque of individual joints. For that purpose, we will upgrade the design to be able to quickly disengage the robot, to tackle only either hip, knee, or ankle joint separately. As for the future work, different controllers will be also tested and compared to each other to obtain the most optimum approach for people who had stroke.
Author contributions
Akim Kapsalyamov: Writing – original draft, Editing, Coding, Dynamic analysis, Experimentation set-up and management, Prototype Design, Data Collection and Interpretation. Shahid Hussain: Experimentation set-up, Conceptualization, Methodology, Reviewing, Supervision, Resources, Project administration. Nicholas A.T. Brown: Conceptualization, Writing – review & editing. Roland Goecke: Writing – review & editing. Prashant K. Jamwal: Conceptualization, Methodology, Writing – review & editing.
Financial support
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Competing interests
The authors declare no conflicts of interest exist.
Ethical approval
Human Research Ethics Committee (HREC) at the University of Canberra has approved this study.
Appendix
The angle definitions from Figure 2 are derived as follows:
\begin{equation}{\unicode[Arial]{x03B1}} =\sin ^{-1} \left(\frac{\mathrm{L}9\sin \left(\Pi -\unicode{x1D6F3}-{\unicode[Arial]{x03B3}} \right)}{\sqrt{\rho _{1}^{2}+\rho _{2}^{2}}}\right)\end{equation}
where
\begin{align*} \rho _{1} & = A_{x}+L9cos(\unicode{x1D6F3})-B_{x} \\ \rho _{2} & = A_{y}+L9sin(\unicode{x1D6F3})-B_{y} \end{align*}
\begin{align}{\unicode[Arial]{x03B2}} =\cos ^{-1} \left(\frac{\left(\rho _{1}\right)^{2}+\left(\rho _{2}\right)^{2}+\mathrm{L}5^{2}-\mathrm{L}7^{2}}{2\mathrm{*}\sqrt{\rho _{1}^{2}+\rho _{2}^{2}}\;*\;\mathrm{L}5}\right)\end{align}
\begin{align}{\unicode[Arial]{x03B1}} _{1}=\tan ^{-1} \left(\frac{\rho _{3}-\mathrm{L}5\sin \left({\unicode[Arial]{x03B1}} +{\unicode[Arial]{x03B2}} -{\unicode[Arial]{x03B3}} \right)}{\rho _{4}+\mathrm{L}5\cos \left({\unicode[Arial]{x03B1}} +{\unicode[Arial]{x03B2}} -{\unicode[Arial]{x03B3}} \right)}\right)\end{align}
where
\begin{gather} \rho _{3}=\mathrm{A}_{\mathrm{y}}+\mathrm{L}9\sin (\unicode{x1D6F3})-\mathrm{B}_{\mathrm{y}} \nonumber\\ \rho _{4}=-\mathrm{A}_{\mathrm{x}}+\mathrm{L}9\cos (\unicode{x1D6F3})+\mathrm{B}_{\mathrm{x}} \\{\unicode[Arial]{x03B2}} _{1}=\cos ^{-1} \left(\frac{\mathrm{L}6^{2}+\mathrm{L}7^{2}-\mathrm{L}8^{2}}{\left(2\ \mathrm{L}6\ \mathrm{L}7\right)}\right)\nonumber\end{gather}
\begin{align}{\unicode[Arial]{x03B1}} _{2}=\tan ^{-1} \left(\frac{{\mathrm{B}_{\mathrm{y}}}+\rho _{6}-\mathrm{L}6\sin \left(\rho _{8}\right)-\mathrm{Gy}}{-\mathrm{B}_{\mathrm{x}}-\rho _{7}-\mathrm{L}6\cos \left(\rho _{8}\right)+\mathrm{Gx}}\right)\end{align}
where
\begin{gather} \rho _{6}=\mathrm{L}5\sin ({\unicode[Arial]{x03B1}} +{\unicode[Arial]{x03B2}} -{\unicode[Arial]{x03B3}} ) \nonumber\\ \rho _{7}=\mathrm{L}5\cos ({\unicode[Arial]{x03B1}} +{\unicode[Arial]{x03B2}} -{\unicode[Arial]{x03B3}} ) \nonumber\\ \rho _{8}={\unicode[Arial]{x03B1}} 1+{\unicode[Arial]{x03B2}} 1\nonumber\\ {\unicode[Arial]{x03B2}} _{2}=\cos ^{-1} \left(\frac{\mathrm{L}1^{2}+\overline{\mathrm{EG}}^{2}-\mathrm{L}4^{2}}{2\mathrm{L}1\ \overline{\mathrm{EG}}}\right)\end{gather}
\begin{align}\overline{\mathrm{EG}}=\sqrt{(\rho _{6}-\mathrm{L}6\sin (\rho _{8})^{2}+(\rho _{7}+L6\cos (\rho _{8}))^{2}}\end{align}
\begin{align}{\unicode[Arial]{x03B1}} _{3}=\tan ^{-1} \left(\frac{\rho _{9}+\rho _{6}-L6\sin \left(\rho _{8}\right)}{\rho _{10}+\rho _{7}+L6\cos \left(\rho _{8}\right)}\right)\end{align}
where
\begin{gather} {\rho _{9}}=G_{y}-L1\sin ({\unicode[Arial]{x03B1}} _{2}+{\unicode[Arial]{x03B2}} _{2})-B_{y} \nonumber\\ {\rho _{10}}=G_{x}+L1\cos ({\unicode[Arial]{x03B1}} _{2}+{\unicode[Arial]{x03B2}} _{2})-B_{x} \nonumber\\ \beta _{3}=\cos ^{-1} \left(\frac{L4^{2}+L2^{2}-L3^{2}}{2\ L4\ L2}\right) \end{gather}
\begin{align}\gamma _{2}=\cos ^{-1} \left(\frac{L4+R_{2}^{2}-R_{1}^{2}}{2R_{2}L4}+{\unicode[Arial]{x03B1}} _{3}\right)\end{align}
\begin{align}\gamma _{6}=\cos ^{-1} \left(\frac{R_{6}^{2}-R_{5}^{2}+L7^{2}}{2R_{6}L7}+\Pi -{\unicode[Arial]{x03B1}} _{1}\right)\end{align}
the substation variables used are as follows:
\begin{gather*}J_{0}=\frac{1}{2}\left[I_{9}\right] \\ J_{1}=\frac{1}{2}\left[m_{1}\frac{L_{1}^{2}}{4}+I_{1}+mp_{1}L_{1}^{2}\right] \\ J_{3}=\frac{1}{2}\left[Ip_{1}\right] \\ J_{5}=\frac{1}{2}\left[m_{5}\frac{L_{5}^{2}}{4}+I_{5}\right] \\ J_{6}=\frac{1}{2}\left[mp_{2}L_{9}^{2}+m_{9}\frac{L_{9}^{2}}{4}+Ip_{2}\right] \\ J_{0}=\frac{1}{2}\left[I_{9}\right] \\ M_{1}=\frac{1}{2}\left[mp_{1}R_{2}^{2}\right] \\ M_{2}=\frac{1}{2}\left[mp_{2}R_{6}^{2}\right] \\ P_{1}=\frac{1}{2}\left[mp_{1}2L_{1}R_{2}\right] \\\end{gather*}
\begin{gather*} C_{1}\left(\theta _{1},\gamma _{2}\right)=\cos \left(\unicode{x1D6F3}_{1}-\gamma _{2}\right) \\ P_{6}=\frac{1}{2}\left[mp_{2}2L_{9}R_{6}\right] \\ C_{6}\left(\theta, \gamma _{6}\right)=\cos \left(\unicode{x1D6F3}-\gamma _{6}\right)\end{gather*}
\begin{align*}G({\theta _1},{\theta _5},{\theta _6},\theta, {\gamma _2}) = &\; {m_1}g\frac{{{L_1}}}{2}\sin {\theta _1} + m{p_1}g({L_1}\sin {\theta _1} + {R_2}\sin ({\gamma _2})) + m{p_2}g({L_9}\sin \theta + {R_6}\sin {\theta _6})\\& + {m_5}g\frac{{{L_5}}}{2}\sin {\theta _5} + {m_9}g\frac{{{L_9}}}{2}\sin \theta\end{align*}
Open access
