Modeling the dynamics and control of rehabilitative exoskeleton with robotic crutches

Abstract

Rehabilitative robotics is an area in the medical field, where one can see a variety of different robotic applications, one of which is the use of robotic exoskeleton in rehabilitation of paraplegics. Current developments are only able to support paraplegics at most and require manual operation of crutches by the patient. In order to overcome this limitation, a theoretical model of a robotic device with actuated robotic crutches is proposed, which can be used to support people with high-level disability, such as quadriplegics who cannot use the existing solutions to perform walking gaits. This work presents kinematic trajectory planning of the proposed model and dynamic analysis of main movement stages. Finally, we present an open-loop control scheme that uses time scaling in order to track the desired joint trajectory of the under-actuated motion stage of crutch swinging. A simulated robotic model has also been developed using Simulink [version R2012b to R2014a]/SimMechanics ([version first generation] developed by Mathworks inc.) environment and has been used for verifying dynamic computations and simulating the robotic device movement under the open-loop control commands of joint torques.

Introduction

Rehabilitative robotics is a continuously developing area of research. In this area, most robotic assisting devices are classified into three major applications:

Using smart prostheses on amputees or instead of severely damaged limbs in order to overcome mobility limitations, as shown by Burck et al. and Au et al. ^1,2

Reestablish limb’s functionality on people who experienced trauma that damaged the limb’s functionality, but the motor function still exists and can be recovered (such as stroke). In this case, it appears that a lot of attempts to recover the limb’s functionality are done using robotic assisting devices, such as the ones shown in the studies by Kahn et al., ³ Sale et al., ⁴ Bernhardt et al., ⁵ Ollinger et al., ⁶ Ferris et al., ⁷ Fan and Yin, ⁸ and Ferris and Lewis. ⁹

Replacing limb’s motor function in cases where it is damaged permanently, but the limb itself is intact (such as paraplegia). A solution to this problem is often the use of wheelchair, ¹⁰ but this solution leads to a variety of problems ranging from decrease in bone and muscle density to problems in cardiac and urinary systems. Another problem is the formation of bedsores, in addition to deterioration in the patient’s mental state. Another possible solution, which helps avoiding the mentioned problems, is the use of external exoskeleton.

Robotic exoskeletons are used in many fields ranging from military and industrial to entertainment and medical uses. In the military field, it is used for enhancing combat soldier’s ability to traverse rough terrains while lifting heavy objects, whereas in the industrial field, exoskeletons are used to support workers while carrying a heavy load. ¹¹ One field, which is continuously developing, is the use of exoskeleton for rehabilitation, where it is used to support people suffering from disabilities while performing gait. HAL, Cyberdyne’s (Tsukuba, Japan) device is an example for using exoskeletons to support and expand the physical capabilities of users, capable of moving by themselves. ¹¹ Rewalk’s, Ekso’s, and also HAL’s exoskeletons are able to support patients suffering from paraplegia while performing gait. ^{10 –16} These solutions, however, require upper-limb and lower-limb coordination, achieved by the patient’s ability to support himself using crutches during the gait, thus causing the solution to be effective only for several disabilities. For higher levels of disability, where the hand function is severely limited, we propose to add active robotic crutches. These crutches are essential for adding support and stabilization in order to maintain upright position of the patient, which is a major difficulty in upright walking of disabled patients. Thus, the primary goal of this article is to conduct a feasibility research on a theoretical planar model of an exoskeleton with an active robotic crutch, which is shown in Figure 1. Feasibility check of the robotic device includes design of kinematic trajectories, built based on field experiments with existing exoskeletons combined with the research developed on bipedal walking robots. ^{10 –13,17,18} These trajectories are fed back into the dynamic equations of motion and used for computing joint actuation torques and ground reaction forces, which are used to assure non-slipping contact while performing gait. The computed control torques are applied at the robot’s joints in open-loop control. A secondary objective of our research is the development of a simulation tool, based on Simulink/SimMechanics environment, ¹⁹ which can be used as a verification tool for various numeric computations. As shown in Figure 1, the robotic device includes an actuated crutch, which is used together with the robotic legs for support and maintaining gait, thus enabling usability for people suffering from wide range of disabilities, including quadriplegics.

Figure 1.

Robotic device carrying a patient.

While the dynamics of walking with manually operated crutches has been analyzed in previous works, ^{20 –22} the proposed concept of robotic exoskeleton with active crutches still requires theoretical feasibility investigation. A work, which is closer in spirit, is the augmented robotic limbs, ²³ which are mainly used for upright stabilization and load reduction for healthy persons carrying heavy loads.

A key factor in our analysis is the fact that the motion stage of crutch swinging is under-actuated, that is, the number of actuation input is less than the number of degrees of freedom (DoFs), which poses limitation on feasible dynamic motions. While set-point stabilization of equilibrium postures of under-actuated systems is well-studied, ²⁴ the problem of tracking a desired trajectory is mainly solved using control methods such as backstepping and feedback linearization, which are applicable for systems with very specific structure. ^25,26 In this work, we use the method of time scaling ^27,28 for obtaining open-loop inputs of actuation torques for tracking a desired kinematic path. This method is suitable for under-actuated systems and has been applied in previous works for multi-contact robotic manipulation, ²⁹ as well as biped robot control. ³⁰

The article is organized as follows. First, formulation of the general robotic model and decomposition of robots gait into three movement stages are presented. Kinematic trajectory planning of the robot’s gait is also shown next, thus followed by dynamic analysis of the main movement stages. A proposed open-loop control method for under-actuated sub-mechanism tracking is then shown, along with robot’s simulation, developed in SimMechanics environment for verifying direct computation results. Finally, the concluding section discusses the limitations of the results and lists some possible future extensions of the research.

Theoretical model

In this section, we present a simplified planar model of the robotic structure and classify the robot’s different stages of movement. Links’ lengths and mass notations are shown in Figure 2(a) and values are shown in Table 1. Each link is considered as a rod with uniform mass, having moment of inertia $I_c=\frac{m{l}^2}{12}$ . All joint motors are considered as point masses with mass value $m_i=0.2\text{\hspace{0.17em}kg}$ and are shown in black circles at the joints’ locations. Each mass index m_i has the same index as the link mass index M_i nearby, that is, m₀ can be found next to M₀, m₁ can be found next to M₁, and so on.

Figure 2.

(a). Gait parameters and (b). Degree of freedom (DoF) notation.

Table 1.

Links’ length and mass values.

Link	Link length (cm)	Mass	Mass value (kg)
L 0	10	M 0	0.2
L 1	40	M 1	3
L 2	40.8	M 2	4
L 3	10	M 3	6
L 4	40	M 4	5
L 5	50	M 5	4
L 6	100	M 6	4

The robot’s DoFs are shown in Figure 2(b). Each joint represents 1-DoF. All joints are revolute except d₆, which is prismatic. Joints description is shown in Table 2. The robotic device, shown in Figure 2(b), has 9-DoF represented by the vector of generalized coordinates q, where

$\mathbf{q}={\left[{\theta }_{1l},{\theta }_{2l},{\theta }_{3l},{\theta }_4,{\theta }_5,d_6,{\theta }_{3r},{\theta }_{2r},{\theta }_{1r}\right]}^T$

1

Table 2.

Joint’s type.

Joint	Type
θ 1l	Left leg ankle joint—passive in steps 2 and 3 and active in 1 (left leg uprises)
θ 2l	Left leg knee joint— active
θ 3l	Left leg hip joint—active
θ 4	Lower back joint—passive
θ 5	Crutch revolute joint—active
d 6	Crutch prismatic joint—active
θ 3r	Right leg hip joint—active
θ 2r	Right leg knee joint—active
θ 1r	Right leg ankle joint—passive in steps 2 and 3 and active in 1 (right leg uprises)

The roles of all joints and their actuation type (active/passive) in our proposed exoskeleton model are shown in Table 2. As shown in Table 2, some joints are actuated, while others are passive. One can see a particular interest in the type of ankle joints θ_1l and θ_1r, which varies depending on the state of foot contact, as detailed in the following assumption

Assumption 1

If the foot is not in contact with the ground during the movement stage—the ankle joint is considered as active. Otherwise, the joint is considered as passive.This assumption is based on the fact the foot’s mass is negligible compared to other links. Thus, a relatively weak torque is required at the ankle joint during foot swinging. (Practically, in exoskeleton devices used by paraplegics, the ankle joint contains a semi-passive spring-based mechanism with a limit stopper. ^31,32 This complication is not considered here, see further discussion in the concluding section.) On the other hand, when the foot is in contact with the ground, it forms a closed kinematic chain whose effective inertia is much larger and requires larger torques. Therefore, the torque of a weak elastic mechanism at the ankle joint becomes negligible and the joint is considered as being passively affected by the closed-chain kinematic constraints. This assumption results in simplifying the kinematic analysis and enabling the determination of the number of DoFs and actuation type in each stage of movement, as detailed next.

Stages of movement

The robot’s motion can be divided into three movement stages:

Movement stage 1: Moving the rear foot (swing foot) forward from toe off to heel strike, as shown in Figure 3.

Movement stage 2: Moving the crutch forward while maintaining the same legs’ position as was in the previous movement stage, as shown in Figure 4.

Movement stage 3: Front foot rotation from heel strike to full contact with the ground, followed by rear foot rotation from full contact with the ground to toe off, as shown in Figure 5.

Figure 3.

Movement stage 1.

Figure 4.

Movement stage 2.

Figure 5.

Movement stage 3.

Partition into these movement stages is determined by transitions in the feet and crutch’s contact states. It is also based on previous works on motion analysis for crutch-assisted rehabilitative exoskeleton, ^12,20 as well as healthy human gaits, ^17,18 and bipedal robots. ³³ Simulation movies of the robot’s different motion stages are included in the online supplementary material. In order to plan robot’s gait, one must determine the number of mechanism’s DoF at each motion stage, which depend on the contact state and kinematic constraints. This can be determined using Grübler’s equation ³⁴

$m=d(n-p-1)+{\displaystyle \sum _{i=1}^p}{f}_i$

2

where m is the mechanism’s number of DoFs also known as mobility, d classifies the problem to spatial or planar, and f_i represents the number of DoFs in each joint. Our model is planar so d = 3 and f_i = 1. n is the number of links and p is the number of joints in the mechanism. We use equation (2) to obtain the DoF in each movement stage and compare to the number of active joints (e.g. actuators). Based on the comparison, each movement stage is then classified as either under-actuated or fully actuated.

Movement stage 1

Using Grübler’s equation (2) on the mechanism shown in Figure 3, moving the robotic structure requires seven actuated DoFs ( $n=\text{1}0,$ p = 10, and f_i = 1). Following assumption 1, in this stage, when left leg uprises, the ankle joint of the swing foot θ_1l is considered as an actuated joint and that of the stance foot θ_1r is considered as a passive joint. Based on that, the robotic structure has seven actuated DoFs (as shown in Table 2). Therefore, the mechanism is fully actuated.

Movement stage 2

Using Grübler’s equation (2), seven DoFs are required to move the mechanism shown in Figure 4 ( $n=\text{\hspace{0.17em}1}0$ , $p=\text{1}0,\text{ and }{f}_i=\text{1}$ ). Following assumption 1, both ankle joints (θ_1l and θ_1r) are passive and the theoretical model has six actuated DoFs. Thus, movement stage 2 is considered as under-actuated.

Movement stage 3

Using Grübler’s equation (2) on the mechanism shown in Figure 5, moving the robotic structure requires six actuated DoFs (n = 11, p = 12, and f_i = 1), so at least 6-DoF must have servos in order to move the mechanism to the required position. As one can see from Table 2, and following assumption 1, both ankle joints (θ_1l and θ_1r) are passive and the theoretical model has six actuated DoFs. Thus, movement stage 3 is considered fully actuated.

As shown in Figure 5, robot’s movement in stage 3 is simpler than stages 1 and 2. Movement in stage 1 is based dynamically on the existence of no-slip and contact constraints as mentioned earlier, and stage 2 is an under-actuated movement stage. Based on these reasons, this article is focused on the analysis of movement stages 1 and 2 only.

Kinematic trajectory planning

This section describes kinematic motion planning used to move the robot along desired trajectories, dictated by the movement stages. For convenience of the analysis, the chosen kinematic trajectories are composed of basic motion primitives that are based on observations from human walking gaits.

Movement stage 1

This stage focuses on moving the swing foot forward. Kinematic trajectories, used to determine robot’s path by inverse kinematics, are defined by motion of the swing foot’s endpoint ( $x_a,z_a,\text{ and }{\theta }_a$ ), pelvis ( $x_{\text{hip}},z_{\text{hip}}$ ), and lower back orientation θ₄. An additional constraint is that the lower back is constantly in upright orientation ${\theta }_{1l}+{\theta }_{2l}+{\theta }_{3l}=0$ , as shown in Figure 2(b). Thus, 7-DoF movements are dictated. Joint space trajectory is the result of using inverse kinematics on these end-effector trajectories (we should mention that the dictated DoFs are not necessarily actuated).

Foot trajectory planning

As shown in Figure 2(a), foot trajectory is represented by $X_a={\left[\begin{array}{c}{x}_a(t),z_a(t)\end{array}\right]}^T.x_a\text{\hspace{0.17em}and\hspace{0.17em}}{z}_a$ are the ankle coordinates, and θ_a is the ankle angle. Planning a trajectory for movement stage 1 is done under the following constraints: (1) continuous first derivative, (2) differentiable second derivative, and (3) zero values of second derivative in the beginning and at the end of the movement stage. Planning a trajectory implementing the mentioned constraints using polynomial interpolation would produce a high-order polynomial, making the computation difficult. Planning the trajectory using third-order spline interpolation would help in avoiding these problems. ¹⁸ Demonstration of the design process using that technique is shown for the ankle angle. The ankle angle θ_a has three intermediate orientations, so the trajectory must pass through these three points during the movement phase. The trajectory is divided into two time intervals, which are then connected to represent the full path. Formulation of the mathematical equations is shown for one of the time intervals

${\theta }_{a_i}\left(t\right)=\left(1-t\right)y_{i-1}+t{y}_i+t\left(1-t\right)(a_i\left(1-t\right)+b_{i}t)$

3

where i = 1,2 represents the interval number and y_i is the value of ${\theta }_a\left(t=t_i\right)$ . Constraints for the ankle angle produce the following

$\begin{array}{l}{\theta }_a(t):\begin{array}{lll}{y}_0=q_b\hfill & ,\hfill & t_0=\left(k-1\right)T_c+T_d\hfill \\ y_1=12q_b\hfill & ,\hfill & t_1=\left(k-1\right)T_c+T_m\hfill \\ y_2=-q_f\hfill & ,\hfill & t_2=k{T}_c\hfill \end{array}\\ \begin{array}{lll}{q}_b=0.38\text{rad}\hfill & ,\hfill & q_f=0.38\text{rad}\hfill \end{array}\end{array}$

4

where $q_b\text{\hspace{0.17em}and\hspace{0.17em}}{q}_f$ are the values of the ankle angle θ_a at the beginning and final time of the motion stage, respectively. The angle values were chosen from observations of existing devices, used by patients with different types of disabilities. ^{10 –13} The choice of time duration of each motion stage is based on previous works on crutch or rail-assisted exoskeleton walking. ^10,35

In order to compute $a_i\text{\hspace{0.17em}and\hspace{0.17em}}{b}_i$ , needed to build the third-order spline trajectory in equation (3), ¹⁸ first one must compute the constants $c_0,c_1,\text{\hspace{0.17em}and\hspace{0.17em}}{c}_2$ . These constants can be computed using the following relation

$\left[\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ a_{21}& a_{22}& a_{23}\\ 0& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ c_2\end{array}\right]=\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\end{array}\right]$

5

where

$\begin{array}{l}{a}_{11}=\frac{2}{t_1-t_0}=\frac{2}{T_m-T_d}\hfill \\ a_{12}=a_{21}=\frac{1}{t_1-t_0}=\frac{1}{T_m-T_d}\hfill \\ a_{22}=2\cdot \left(\frac{1}{t_1-t_0}+\frac{1}{t_2-t_1}\right)=2\cdot \left(\frac{1}{T_m-T_d}+\frac{1}{T_c-T_m}\right)\hfill \\ a_{23}=a_{32}=\frac{1}{t_2-t_1}=\frac{1}{T_c-T_m}\hfill \\ a_{33}=\frac{2}{t_2-t_1}=\frac{2}{T_c-T_m}\hfill \\ \hfill \\ e_1=3\cdot \frac{y_1-y_0}{{\left(t_1-t_0\right)}^2}=3\cdot \frac{-0.5q_b}{{\left(T_m-T_d\right)}^2}\hfill \\ e_2=3\cdot \left(\frac{y_1-y_0}{{\left(t_1-t_0\right)}^2}+\frac{y_2-y_1}{{\left(t_2-t_1\right)}^2}\right)\hfill \\ =3\cdot \left(\frac{-0.5q_b}{{\left(T_m-T_d\right)}^2}-\frac{\left(q_f+0.5q_b\right)}{{\left(T_c-T_m\right)}^2}\right)\hfill \\ e_3=3\cdot \frac{y_2-y_1}{{\left(t_2-t_1\right)}^2}=-3\cdot \frac{q_f+0.5q_d}{{\left(T_c-T_m\right)}^2}\hfill \end{array}$

These constants can then be used to compute $a_i\text{\hspace{0.17em}and\hspace{0.17em}}{b}_i$ using the following

$\begin{array}{l}{a}_i=c_{i-1}\left(t_i-t_{i-1}\right)-\left(y_i-y_{i-1}\right)\hfill \\ b_i=-c_i\left(t_i-t_{i-1}\right)+\left(y_i-y_{i-1}\right)\hfill \end{array}$

6

All the parameters used for building the trajectory are mentioned and explained in Table 3. The same computation technique is applied for obtaining the ankle’s x and z-coordinates. Constraints for ankle’s coordinates, using the same t values as before, are shown in the following expressions

$\begin{array}{l}{x}_a(t):\begin{array}{l}{y}_0=(k-1)d+L_0\text{sin}{q}_b+l_{af}\left(1-\text{cos}{q}_b\right)\hfill \\ y_1=(k-1)d+L_{ao}\hfill \\ y_2=(k+1)d-L_0\text{sin}{q}_f-l_{ab}\left(1-\text{cos}{q}_f\right)\hfill \end{array}\\ z_a(t):\begin{array}{l}{y}_0=l_{af}\text{sin}{q}_b+L_0\text{cos}{q}_b\hfill \\ y_1=H_{ao}\hfill \\ y_2=L_0\text{cos}{q}_f+l_{ab}\text{sin}{q}_f\hfill \end{array}\end{array}$

7

Table 3.

Gait planning parameters.

Parameter	Value	Description
Tc	$\frac{v}{2d}=1.8\frac{s}{\text{step}}$	Period necessary for one walking step
Td	$2T_c=3.6s$	Interval of the double-support phase
T crutch	$0.2T_c=0.36s$	Period necessary for moving the crutch forward
d	0.27 m	Length of one step
v	0.3 m/s	walking velocity
k	$0,1,2,\dots [\text{step}]$	kth step
Tm	$\frac{t_d+t_c}{2}$ = 2.7s	Time when the swing foot is at its highest point
Lao	0.25 m	Position of the highest point along the x-axis
Hao	0.16 m	Position of the highest point along the z-axis

The reader can refer to Table 3 for explanation about the equation parameters.

Pelvis and lower back trajectory planning

During the walking cycle, the pelvis moves to balance patient’s body. Robot’s pelvis motion was based on an arched trajectory, followed by stance leg with locked knee acting as an inverted pendulum. ^{13,36 –38} θ₄ trajectory was planned using the third-order polynomial, based on the following initial and final values

$\begin{array}{lll}{\theta }_4(t_0)=\frac{\pi }{2}-\frac{\pi }{2.3}\hfill & ,\hfill & {\dot{\theta }}_4(t_0)=0\hfill \\ {\theta }_4(t_f)=0\hfill & ,\hfill & {\dot{\theta }}_4(t_f)=0\hfill \end{array}$

8

where t₀ and t_f are the initial and final times of the motion stage, respectively. Stance leg angle θ_sa acts as opposite to lower back angle θ₄, in the same angle’s range, according to relation

${\theta }_{\text{sa}}={\theta }_4-C$

9

while

$C=\frac{\pi }{2}-\frac{\pi }{2.3}$

The value in equation (9) is also taken from the works of Esquenazi et al., ^10,12 Mikolajewska and Mikolajewski, ¹¹ and Cass. ¹³ Pelvis and lower back movement is determined, using the fact that L₂ and L₁ are constant, on the following relations

$X_{hip}=\left[\begin{array}{c}{x}_{\text{hip}}\\ z_{\text{hip}}\end{array}\right]=\left[\begin{array}{c}kd+(L_2+L_1)\text{sin}{\theta }_{\text{sa}}\\ (L_2+L_1)\text{cos}{\theta }_{\text{sa}}\end{array}\right]$

10

Crutch’s endpoint is fixed during the movement.

Movement stage 2

This stage focuses on moving the crutch forward. Robot’s legs are fixed and form a closed kinematic chain (rigid link acting as parallel sub-mechanism), thus joint values θ_2l, θ_3l, θ_2r, and θ_3r are fixed to their final orientation given at the end of the previous stage. Crutch’s movement is determined using ${\theta }_4(t)$ , planned using the third-order polynomial as was explained previously and crutch’s end-effector trajectory ( $x_c,z_c$ ), thus dictating 7-DoF movements. Third-order polynomial for ${\theta }_4(t)$ is computed using the constraints shown in the following

$\begin{array}{lll}{\theta }_4(t_0)=0\hfill & ,\hfill & {\dot{\theta }}_4(t_0)=0\hfill \\ {\theta }_4(t_f)=\frac{\pi }{2}-\frac{\pi }{2.3}\hfill & ,\hfill & {\dot{\theta }}_4(t_f)=0\hfill \end{array}$

11

Crutch’s end effector was planned using the third-order spline according to the constraints, which are given by

$\begin{array}{l}{x}_c(t):\begin{array}{lll}{y}_0=\left(k+1\right)d\hfill & ,\hfill & t_0=k{T}_c\hfill \\ y_1=\left(k+1\right)d+0.5d\hfill & ,\hfill & t_1=k{T}_c+0.5T_{\text{crutch}}\hfill \\ y_2=\left(k+2\right)d\hfill & ,\hfill & t_2=k{T}_c+T_{\text{crutch}}\hfill \end{array}\\ z_c(t):\begin{array}{lll}{y}_0=0\hfill & ,\hfill & x_c(t_0)=\left(k+1\right)d{}_c\hfill \\ y_1=L_0\hfill & ,\hfill & x_c(t_1)=\left(k+1\right)d+0.5d\hfill \\ y_2=0\hfill & ,\hfill & x_c(t_2)=\left(k+2\right)d\hfill \end{array}\end{array}$

12

This concludes the kinematic trajectory planning for motion stages 1 and 2.

We now study the dynamics at movement stages 1 and 2 and compute the required actuation torques. While the motion of each movement stage is dynamic, we impose that the transitions between stages are reached at complete stop. This assumption is common in crutch-assisted motion and appears also in previous work. ¹² Thus, we do not consider here the effect of impacts at ground touchdown of feet or crutch, which has been studied by Carpentier et al. ³⁹ and Font-Llagunes et al. ⁴⁰

Dynamic analysis of movement stage 1

In this stage, the robot swings the rear leg forward while being supported by a crutch, placed on the ground with full sticking contact. The robot’s equations of motion are formulated using constrained Lagrange equations ³⁴ as

$M(q)\ddot{q}+B\text{(}q,\dot{q}\text{)}+G(q)\text{=}{F}_q+W{(q)}^{\text{T}}\Lambda$

13

where $F_q={\left[\begin{array}{c}{\tau }_{1l},{\tau }_{2l},{\tau }_{3l},{\tau }_4,{\tau }_5,f_6,{\tau }_{3r},{\tau }_{2r},{\tau }_{1r}\end{array}\right]}^T$ (generalized forces vector), q is the robot’s coordinates defined in equation (1) and Figure 2(b), $M(\mathbf{q})$ is the Robot’s inertia matrix, $B(q,\dot{q})$ is the vector-containing velocity-dependent forces, and Λ is the vector of Lagrange multipliers, enforcing contact and no-slip constraints (holonomic constraints).

The holonomic constraints represent the fact that the endpoint of the crutch is fixed. We express the vector of crutch’s endpoint using generalized coordinates $q:r_c={[x_c-x_0,z_c-z_0]}^T$ , x and z expressions are shown in the following

$z_c-z_0=L_1\text{cos}{\theta }_{1l}+L_2\text{cos}\left({\theta }_{1l}-{\theta }_{2l}\right)+L_3\text{cos}({\theta }_{1l}-\left({\theta }_{2l}+{\theta }_{3l}\right))+L_4\text{cos}({\theta }_{1l}+{\theta }_4-\left({\theta }_{2l}+{\theta }_{3l}\right))+\left(L_5+d_6\right)\text{cos}({\theta }_{1l}+{\theta }_4+{\theta }_5-\left({\theta }_{2l}+{\theta }_{3l}\right))+L_0=0x_c-x_0=L_1\text{sin}{\theta }_{1l}+L_2\text{sin}\left({\theta }_{1l}-{\theta }_{2l}\right)+L_3\text{sin}({\theta }_{1l}-\left({\theta }_{2l}+{\theta }_{3l}\right))+L_4\text{sin}({\theta }_{1l}+{\theta }_4-\left({\theta }_{2l}+{\theta }_{3l}\right))+\left(L_5+d_6\right)\text{sin}({\theta }_{1l}+{\theta }_4+{\theta }_5-\left({\theta }_{2l}+{\theta }_{3l}\right))=0$

14

Then, the constraint can be expressed as $W(q)\dot{q}=0,\text{ where }W(q)=d{r}_c/dq$ is a 2 × 9 Jacobian matrix of r_c with respect to q. The vector of Lagrange multipliers is $\Lambda =[{\lambda }_t,{\lambda }_n{]}^T$ which represent tangential and normal components of the ground reaction force at the crutch’s endpoint. The inertia matrix in equation (13), as well as the vector of velocity-dependent terms (such as Coriolis and centrifugal forces), and vector of gravitational forces acting on the robot, are evaluated using the following relations ³⁴

$M_{ij}=\frac{{\partial }^{2}T}{\partial {\dot{q}}_i\partial {\dot{q}}_j},B_i=\frac{{\partial }^{2}T}{\partial {\dot{q}}_i\partial q_j}{\dot{q}}_j-\frac{\partial T}{\partial q_i},G_i=\frac{\partial V_g}{\partial q_i}$

15

where T is the robot’s kinetic energy, computed using the following

$T={\displaystyle \sum _{j=1}^N}\frac{1}{2}{}_{\text{Mj}}\left(v_j\cdot v_j\right)+12{\omega }_j{}^{\text{T}}{I}_{cMj}{\omega }_j$

16

where $v_j/{\omega }_j$ is the linear/angular center of mass (C.O.M) velocities of joint/link j, and I_cMj is the inertia moment of joint/link j, which equals one of the following values

$I_{cMj}=\{\begin{array}{cc}\frac{M_j{L}_j^2}{12}& \text{robot’s\hspace{0.05em}\hspace{0.17em}links}\\ 0& \text{robot’s\hspace{0.17em}joints}\end{array}$

17

where M_j is joint/link j mass, and V_g is the robot’s potential energy, computed using the following

$V_g\text{=}{\displaystyle \sum _{j\text{\hspace{0.17em}=\hspace{0.17em}1}}^N}-{}_{\text{Mj}}{r}_j\cdot g$

18

where r _j is joint/link j C.O.M position vector. q_i is the ith generalized coordinate from q (equation (13)). Joint’s velocity $\dot{q}$ and acceleration $\ddot{q}$ are obtained using differentiation of the planned kinematic trajectories. Having the desired joint’s position/velocity/acceleration and required constraints, one can compute the required actuation torques and forces using inverse dynamics method, ³⁴ explained as follows.

The two holonomic constraints (no-slip and contact) remove 2-DoF from the 9-DoF, shown in the motion equations (equation (13)). The robotic structure has two unactuated and passive joints ${\theta }_4\text{\hspace{0.17em}and\hspace{0.17em}}{\theta }_{1l}$ . A close examination of the constrained motion equations shows that seven actuation torques and two constraint forces are needed in order to maintain the desired movement. These forces and torques are computed from the motion equations, evaluation and computation process is shown in the following

$\begin{array}{l}M(q)\ddot{q}+B(q,\dot{q})+G(q)=F_q+W{(q)}^{\text{T}}\Lambda =E\tau +W^{\text{T}}\Lambda \end{array}$

where

$\tau ={\left[\begin{array}{l}{\tau }_{2l}\hfill \\ {\tau }_{3l}\hfill \\ {\tau }_5\hfill \\ f_6\hfill \\ {\tau }_{3r}\hfill \\ {\tau }_{2r}\hfill \\ {\tau }_{1r}\hfill \end{array}\right]}_{7\times 1},E={\left[\begin{array}{lllllll}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]}_{9\times 7}$

19

Equation (19) can be rearranged as

$M(q)\ddot{q}+B(q,\dot{q})+G(q)=\underset{{\tilde{E}}_{\text{9\hspace{0.17em}}\times \text{\hspace{0.17em}9}}}{\underbrace{\left[E{W}^{\text{T}}\right]}}{\left[\begin{array}{l}\tau \hfill \\ \Lambda \hfill \end{array}\right]}_{9\times 1}$

20

The matrix $\tilde{E}$ in equation (20) is invertible; hence, the actuation torques and constraint forces can be extracted as

$\left[\begin{array}{l}\tau \hfill \\ \Lambda \hfill \end{array}\right]={\tilde{E}}^{-\text{1}}(\text{M}\left(q\right)\ddot{q}+B\left(q,\dot{q}\right)+G\left(q\right))$

21

By substituting the desired position, velocity, and acceleration vectors ( $q(t),\dot{q}(t),\ddot{q}(t)$ ) into equation (21), one can obtain the actuation torques/forces and constraint forces which are required for maintaining the trajectory defined in movement stage 1 (inverse dynamics). The resulting actuation torques and forces are shown in Figure 7 after comparison with those obtained from the robotic model, developed in SimMechanics environment from a collection of links and joints (see Figure 6).

Figure 6.

SimMechanics model of movement stages 1 (a) and 2(b). (a) is the left side and (b) the right side.

Source files and animation movies of the SimMechanics model’s motion are included in the online supplementary material. The SimMechanics environment provides an independent simulation of the robot’s dynamics, which is not based on our equations of motion. From the plots in Figure 7, the perfect matching between the results of the SimMechanics model and numerical integration of our equations of motion provides a strong indication for the correctness of our theoretical analysis.

Figure 7.

Actuation torques/forces during movement stage 1.

As shown in Figure 3, during movement stage 1, the robot is being supported by two contact points—crutch and stance foot. In order to ensure no-slip at both contact points, one must compute reaction forces occurring from crutch and foot interaction with the ground. While reaction forces in crutch’s tip, obtained from the simulation tool, were compared to those computed using direct computation of Λ(t), in stance foot interaction, reaction forces were obtained only from the simulation tool. Resulting normal contact forces in crutch’s tip and stance foot are shown in Figure 8, while SimMechanics model of stage 1 is shown in Figure 6(a).

Figure 8.

Normal contact force at (a) crutch’s tip and (b) stance foot.

As seen in Figure 8(a), contact force in normal direction is positive ( $\left({\lambda }_n>0\right)$ during all the movement stage, therefore contact is maintained. Figure 8(a) also shows that contact force magnitude is small in the beginning and increasing towards the end. This is logical because in the beginning of the movement stage, the patient supports himself partially using the crutch, while in the end, the crutch is being used for full patient’s support. Same logic holds for Figure 8(b), in the beginning of the movement stage, all the patient’s weight is on the stance foot, while in end of the movement stage, patient’s weight on the stance foot is partial. While tangential component of the contact force (λ_t) does not provide important information, a lot of information can be obtained from constraint forces ratio ${\lambda }_t/{\lambda }_n$ . No-slip contact is maintained by enforcing the inequality of Coulomb’s friction law, ³⁴ given by

$|{\lambda }_t|\le \mu {\lambda }_n\Rightarrow \mu \ge \frac{|{\lambda }_t|}{{\lambda }_n}$

22

By computing the ratio between two constraint forces (tangential and normal), one can obtain critical friction coefficients, both for crutch’s tip and stance foot, which are shown in Figure 9. As shown in Figure 9, applying a friction coefficient greater than μ = 0.2878 for crutch’s tip and a friction coefficient greater than μ = 0.1197 for stance foot will result in no-slip for both contacts during movement stage 1. These two values are physically achievable in practice.

Figure 9.

Ratio of tangential-to-normal contact force at (a) crutch’s tip and (b) stance foot.

Under-actuated dynamics and control for movement stage 2

In this stage, the robot swings the crutch forward in order to advance and prepare for the next step in the walking cycle, while maintaining the same legs position as was in the end of the previous movement stage. Given the rigid chain created by the robot’s legs (while all joint angles are kept constant), the full dynamic model could be reduced to 3-DoF, shown in Figure 10. Using equation (13) with

$q=\left[\begin{array}{l}{\theta }_4\hfill \\ {\theta }_5\hfill \\ d_6\hfill \end{array}\right]$

Figure 10.

Simplified model.

One obtains three motion equations with only two actuated DoFs $F_q={\left[\begin{array}{ll}{\tau }_5\hfill & f_6\hfill \end{array}\right]}^T$ , resulting in an under-actuated system. Lower back joint $\left({\theta }_4\right)$ is passive, so zero torque must be applied at that joint.

First, one should check if the crutch can track the desired motion trajectory described earlier, while using only two actuated joints $\left({\theta }_5,d_6\right)$ and a passive joint (both at the robot and patient) which must have zero torque, that is, τ₄ = 0. After accomplishing that, one must check if the crutch’s tip is able to detach from the ground at the beginning of this motion stage, using the same actuation torques computed in this section. The challenge here is to prescribe a trajectory for 3-DoF, while only two are actuated and the lower back joint $\left({\theta }_4\right)$ is passive. Additionally, one has to satisfy bounds on torques at the actuated joints. Practical bounds were chosen as 4100 Nm on τ₅ and 2700 N on f₆. These bounds were obtained by computing the maximal static loads and multiplying by a safety factor of 2.5. The key idea is to choose the appropriate time scaling for motion along this trajectory in a way that is achievable with the passive joint and also satisfies the bounds on actuation torques. The details of this method are as follows.

Under-actuated reference trajectory tracking during crutch swing phase

We define geometric parameterization of joints’ trajectory $q(s)$ . Time variation of s(t) then dictates the rate of motion along the trajectory. That is, we follow a specified trajectory of joint angles q(s), but the speed of motion along it can vary over time

$q(t)=\left[\begin{array}{l}{\theta }_4\hfill \\ {\theta }_5\hfill \\ d_6\hfill \end{array}\right]=q(s(t)),\dot{s}(t)>0$

23

The system must satisfy the following motion equations

$\text{M}\left(q\right)\ddot{q}+B\left(q,\dot{q}\right)+G\left(q\right)\text{=}{F}_q$

24

Defining equation (24) to be s-dependent through the following relations

$q=q(s\left(t\right))\dot{q}\left(t\right)=\frac{\partial q}{\partial t}=\frac{\partial q}{\partial s}\frac{\partial s}{\partial t}=q^{\prime }\left(s\right)\dot{s}\ddot{q}\left(t\right)=\frac{\partial }{\partial t}(\dot{q}\left(t\right))=\frac{\partial }{\partial t}(q^{\prime }\left(s\right)\dot{s})=q^{\prime }\left(s\right)\ddot{s}+\frac{\partial }{\partial s}(q^{\prime }\left(s\right)){\dot{s}}^2=q^{\prime }\left(s\right)\ddot{s}+q″\left(s\right){\dot{s}}^2$

we obtain

$\begin{array}{l}M(q\left(s\right))\cdot \left(q^{\prime }\left(s\right)\ddot{s}+q″\left(s\right){\dot{s}}^2\right)+B(q\left(s\right),q^{\prime }\left(s\right)\dot{s})\hfill \\ +G(q\left(s\right))=E\cdot u\left(s\right)=\left[\begin{array}{c}0\\ {\tau }_5\\ f_6\end{array}\right]\hfill \end{array}$

25

while

$E=\left[\begin{array}{lll}0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill \end{array}\right],u\left(s\right)=\left[\begin{array}{c}{\tau }_4(s\left(t\right))\\ {\tau }_5(s\left(t\right))\\ f_6(s\left(t\right))\end{array}\right]$

26

The first row of equation (25) contains the requirement that the torque at the passive hip joint is identically zero. Left-multiplication of equation (25) by $b=(1,0,0{)}^T$ and rearranging terms, one obtains

$A\left(s\right)\ddot{s}\left(t\right)+B\left(s\right){\dot{s}}^2\left(t\right)+C\left(s\right)=0$

27

where A, B, and C are s-dependent terms (which will not be shown here for brevity). Equation (27) is a second-order nonlinear scalar differential equation whose coefficients are s-dependent. A solution for the time scaling s(t) is obtained by numerical integration with events (ode45+event function) following the state-space representation

$\begin{array}{l}{x}_1=s\hfill \\ x_2=\dot{s}={\dot{x}}_1\hfill \\ \Rightarrow \dot{x}=\left[\begin{array}{l}{\dot{x}}_1\hfill \\ {\dot{x}}_2\hfill \end{array}\right]=\left[\begin{array}{l}{x}_2\hfill \\ -\frac{\left(B\left(x_1\right)x_2{}^2+C\left(x_1\right)\right)}{A\left(x_1\right)}\hfill \end{array}\right]\hfill \end{array}$

28

while:

Initial condition $s(t=t_0)=s_0$ and end condition $s(t=t_f)=s_f$ are known. The final time t_f is not known a priori, and is obtained by stopping the integration when reaching s_f, using event function.

$\dot{s}(t=t_0)={\dot{s}}_0$ is not unique; therefore, it is obtained using numerical search on interval limited by actuation torques practical bounds (4100 Nm on τ₅ and 2700 N on f₆) and requirement that $\dot{s}(t)\ge 0$ for all t.

Actuation torques/forces vector ${\tilde{F}}_q$ is computed by substituting the solution s(t) into equation (25).

Using the s(t) time scaling obtained from the solution of equation (27), the following results are obtained. Using open-loop control, based on the state-space representation shown in the following

$\begin{array}{l}{X}_1=q=\left[\begin{array}{l}{\theta }_4\hfill \\ {\theta }_5\hfill \\ d_6\hfill \end{array}\right],X_2={\dot{X}}_1=\left[\begin{array}{l}{\dot{\theta }}_4\hfill \\ {\dot{\theta }}_5\hfill \\ {\dot{d}}_6\hfill \end{array}\right]\hfill \\ \Rightarrow \{\begin{array}{l}{\dot{X}}_1=X_2\hfill \\ {\dot{X}}_2=\ddot{q}=M^{-\text{1}}\left[{\tilde{F}}_q-B-G\right]\hfill \end{array}\hfill \end{array}$

29

Joint’s actuation torques and forces were fed as an input to model’s actuators (τ₅ and f₆), resulting in throwing the crutch forward. Robot’s motion was simulated and verified using SimMechanics simulation, under the same open-loop control (forward dynamics method). Crutch’s tip trajectory, obtained from SimMechanics, was compared to desired trajectory and the trajectory obtained from direct computation. The results are shown in Figure 11. The successful tracking of the desired trajectory with under-actuated control proves the correctness of our time-scaling method. Agreement between numerical integration results and the SimMechanics model again verifies the correctness of our equations of motion.

Figure 11.

(a) s-method simulation results (ds/dt vs. time), (b) crutch tip tracking, (c) actuation torques/forces, and (d) lower back angle versus time.

Figure 11(c) shows that crutch’s actuation torques, obtained using the time-scaling method, are within actuation bound limits which were defined. One can also see that joint actuation torques for movement stage 1, shown in Figure 7, are also within these bounds. Because of the actuation torques magnitude in movement stage 1 (less than 1% of the smallest actuation torque bound), practical bounds are not shown in Figure 7. The difference in actuation torques magnitude is due to movement nature: In stage 1, movement nature is quasi-static and reaction forces carry most of the robot’s weight, and in stage 2, movement nature is dynamic, swinging the crutch creates a movement of trying to swing an inverted pendulum and thus the required torques are smaller.

Switch from full contact to flight phase in t = t₀

After succeeding in tracking the crutch’s desired trajectory under the shown constraints, mentioned previously, evaluation of the equations used for contact separation in t = t₀ is shown. These equations were used to create the permissible region of initial actuation force and torque. This region guarantees contact separation while applying zero torque to lower back joint (θ₄), which is passive. Two possible cases should be considered while switching from contact to flight phase:

1. Crutch is in sticking contact—contact + no-slip constraints. Using the motion equation (13) while maintaining no-slip contact results in constraint forces shown in the following

$\begin{array}{l}\Lambda \left(q,\dot{\mathbf{q}}\right)={\left(W{M}^{-\text{1}}{W}^T\right)}^{-\text{1}}\left(W{M}^{-\text{1}}\left(B+G-F_q\right)-\dot{W}\dot{\mathbf{q}}\right)\hfill \end{array}$

30

For further explanation, one can refer to the study by Murray et al. ³⁴ Note that the matrix $W{M}^{-\text{1}}{W}^{\text{T}}$ in equation (30) is always invertible.

2. Crutch is in slipping contact—contact constraints while slipping occurs. Using Coulomb’s friction law ${\lambda }_t=-\sigma \mu {\lambda }_n,\text{ where }\sigma \text{\hspace{0.17em}\hspace{0.05em}=\hspace{0.17em}sign}\left(V_t\right)$ , on equation (13) and time differentiation of $W_n\cdot \dot{q}=0$ , we get the constraint forces, shown in the following

$\begin{array}{l}A\left(q,\dot{\mathbf{q}}\right)\text{=}{W}_n{M}^{-\text{1}}\left(B+G-F_q\right)-{\dot{W}}_n\dot{q}\hfill \\ B\left(q\text{,}\mu \text{,}\sigma \right)\text{=}{W}_n{M}^{-\text{1}}{\left(W_n-\mu \sigma W_t\right)}^{\text{T}}\hfill \\ {\lambda }_n\text{=}\frac{A}{B}\text{,\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{\lambda }_{\text{t}}\text{=}-\sigma \mu {\lambda }_n\hfill \end{array}$

31

In order to guarantee detachment at t = t₀, one has to require that both contact modes are infeasible, that is, λ_n ≤ 0. By evaluating equations (30) and (31) in t = t₀ ( $q(t=t_0)=q_0\text{ and }\dot{q}(t=t_0)=0$ ) and satisfying the requirement λ_n ≤ 0 (contact separation), we obtain two equations of the form

$a{\tau }_5(t_0)+b{f}_6(t_0)+c\le 0$

32

In Figure 12, we show the region created by intersection of the two inequality constraints (equation (32)) obtained for the two contact modes. One can also see the point that represents actuation torques/forces using the time-scaling method in t = t₀ is marked by “*”. Figure 12 shows that the point is within the permissible region, thus contact separation at the beginning of this motion stage is guaranteed.

Figure 12.

Permissible actuation torques/forces region.

Summary

In this article, we performed a feasibility check on a simple planar model of a newly proposed robotic exoskeleton with active crutches. As part of the feasibility check, we numerically investigated the robot’s dynamics while trying to track our planned kinematic trajectories and formulated open-loop control inputs actuation forces and torques at the joints. We have found that for one of the major movements stages, the model’s movement can be achieved using sticking contact constraint forces. In the second major movement stage, we proposed an open-loop control scheme, designed for tracking of an under-actuated system. Our formulated open-loop control inputs of actuation forces and torques at the joints, based on time variation of joint’s trajectory geometric parameterization, achieved contact separation in the first moment and tracking of the planned kinematic trajectory. Thus, feasibility of implementing the gait’s major movement stages has been theoretically proven.

Discussion

We now briefly discuss the main limitations of our model and present several directions for future extensions of the research. First, our kinematic model of the robot is simplified. In particular, it does not account for details of the foot contact and ankle joint. Typical designs of ankle joint in lower-limb exoskeletons consist of a semi-passive dorsiflexion mechanism ^31,32 with a torsion spring and mechanical stoppers. Our model does not account for details of such a mechanism and uses instead the simplification of assumption 1, which states that the ankle joint is active at some motion stages and passive at others. In addition, it is assumed that the front foot is in heel contact during motion stage 2 of crutch swing, which is not fully realistic in healthy human gaits. These assumptions have been made in order to simplify the kinematic analysis of actuated DoF, without requiring a combination of both kinematic and dynamic considerations. Development of more detailed kinematic and dynamic models of the ankle and foot mechanics is a key challenge, which is essential for future realization of our proposed robotic exoskeleton. In addition, future generalizations should consider theoretical analysis of the proposed assistive device operated in transient motion of sitting and standing, which are highly essential tasks for disabled patients in everyday life.

Second, our kinematic model is limited to planar motion. The proposed kinematic trajectories are also simplified and are not fully based on human motion. We are currently working on the extension of the analysis in order to account for kinematic trajectories, which are adapted from measurement data of human motion collected during walking gaits. ⁴¹ Nevertheless, reliable projection of the measured data onto low-dimensional models is itself a challenging task. ³³ Additionally, our model considered a prismatic actuator for the crutch for simplicity of the analysis, and more general configurations of active crutches should be investigated. Finally, our analysis of the robot’s dynamics does not account for interaction between the human and the exoskeleton due to possible kinematic discrepancies and antagonistic efforts. Moreover, we provide only open-loop control of the actuation torques that are required in order to generate the desired kinematic trajectories. Augmenting our work with feedback control for stabilization of trajectory tracking is essential for any practical implementation and for improving the robustness under external perturbations and unmodeled effects. In lower-limb exoskeletons for paraplegics, ^14,16 an essential aspect is reactive control, which adapts to forces applied by the human and possibly regulates force-related quantities such as foot-ground center of pressure. ⁴² In the case of severely limited patients such as quadriplegics, the role of reactive feedback is still an open unknown issue. In summary, all differences abovementioned between the current simple theoretical model and real behavior of a rehabilitative device in practice should be resolved in future works by conducting clinical trials with disabled subjects in physical experiments.

Supplemental material

supplamental_material - Modeling the dynamics and control of rehabilitative exoskeleton with robotic crutches

supplamental_material for Modeling the dynamics and control of rehabilitative exoskeleton with robotic crutches by Adi Cohen and Yizhar Or in International Journal of Advanced Robotic Systems