2.1. Mechanical design

The mechanical design of the proposed 1-DOF hip exoskeleton is represented in Fig. 1. The main components of the exoskeleton include the Bowden cable, the elastic spring, the thigh linkage, the thigh support part, the torso support part, the sliding block, the angle sensor, and the force sensor. The pulling force $F_P$ generated by the PAMs is transmitted through the Bowden cable to drive the thigh linkage to rotate counterclockwise, and the force $F_T$ generated by the spring actuates the linkage to rotate clockwise. The Bowden cable transmission allows remote placement of the PAM actuators so as to minimize gravitational and inertial loading on the user. The sliding block connected with the torso support component and the thigh support components can be adjusted to match the height of the subject before the exoskeleton is put on, which enables the rotation axis of the robot to be aligned with the axis of the human hip joint. Furthermore, there is a sliding block at the back of the torso support component so that the width of the exoskeleton can be adjusted to match the width of the subject’s trunk.

Figure 1. The mechanical design of the lower limb exoskeleton.

The side view of the exoskeleton is represented in Fig. 2. When the exoskeleton is at the flexion state shown in Fig. 2(a), the elastic spring is at the original length, and the thigh linkage is at the end of the swing phase. When the exoskeleton is at the extension state shown in Fig. 2(b), the thigh linkage is at the beginning of the swing phase and also the ending of the standing phase.

Figure 2. The side view of the exoskeleton.

When the pneumatic actuating part acts on the thigh linkage through the Bowden cable, the thigh linkage of the exoskeleton rotates counterclockwise, and the spring is stretched to store the elastic potential energy. When the pneumatic muscle is deflated, the force in the Bowden cable is released. Thus, the spring releases its elastic potential energy to actuate the thigh linkage rotates clockwise and returns to the initial position.

2.2. Control architecture

As represented in Fig. 3, an xPC control system is constructed in this paper to realize the real-time control of the proposed lower limb exoskeleton. The control architecture consists of the following three parts: (1) the compressed air part (yellow); (2) the driving force part (green); (3) the circuit part (blue). In the compressed air part, the air compressor (Fusheng, SA08+) provides the compressed air to the proportional valve (Festo, VPPM-6L-L-1-G18-0L10H-V1N), and the proportional valve gives the compressed air with the desired air pressure to the PAMs (Festo, DMSP-40-600N-RM-RM). The air pressure of the valve is proportional to the control signal. In the driving force part, the PAMs actuated by the compressed air generate a pulling force transmitted by the Bowden cable to actuate the exoskeleton hip joints’ movement. In the control circuit part, the control strategy is programmed in the host PC with the Simulink module of Matalb2013b, and the xPC target (Advantech, IPC-610L) downloads the control strategy in the form of machine codes from the host PC. Then, the xPC target equipped with the data acquisition card (NI, PCI-6229) receives the real-time data of the angle sensors (Rep-avago,18S-1024-2MD-2-15-00E) and calculates the control signal. The control signal is then sent to the proportional valve, and the inside air pressure of the PAMs can be controlled. Thus, the real-time control of the exoskeleton’s movement is realized.

Figure 3. The real-time control architecture of the lower limb exoskeleton system.

2.3. Dynamic model

To facilitate the dynamic analysis, we construct the schematic diagram of the exoskeleton represented in Fig. 4. The motion of the thigh linkage is equated as a rigid, uniformly distributed mass connecting linkage BC rotating around the fixed axis O.

Figure 4. The schematic diagram of the exoskeleton dynamics.

The black dotted line segments AB, BC, and DE represent the positions of the elastic spring, the thigh rod, and the Bowden cable at the initial flexion state, respectively. The point O is the axis of rotation of the thigh linkage. $l_T$ denotes the length of OB, and $m$ represents the mass of the whole thigh linkage. The dotted line OZ denotes the reference line perpendicular to the ground. $\theta _0$ denotes the angle between the thigh linkage and OZ at the initial state, and $\theta$ indicates the rotation angle of the thigh linkage relative to the initial state. According to Lagrange’s theorem, the dynamics of this single-link model can be expressed as:

(1) \begin{equation} \tau _{hip}=m\ddot{\theta }(l_P+l_T)^2+I\ddot{\theta }+\Delta G, \end{equation}

where the first term is related to the linear velocity, the second term is associated with the angular velocity, and the third term is related to the gravitational potential energy. The effect of friction on the joint motion would be regarded as modeling uncertainties later.

The mass of a unit length of the linkage can be expressed as:

(2) \begin{equation} \lambda =\frac{m}{l_P+l_T}. \end{equation}

The moment of inertia $I$ can be expressed as:

(3) \begin{equation} \begin{aligned} I&=\int _{-l_P}^{l_T}\lambda x^2dx\\[5pt] &=\frac{m}{3(l_P+l_T)}(l_P^3+l_T^3) \end{aligned} \end{equation}

The gravitational potential energy term $\Delta G$ can be expressed as:

(4) \begin{equation} \Delta G=[\cos{\theta _0}-\cos{(\theta -\theta _0)}]\frac{mgl_P(l_P-l_T)}{2(l_P+l_T)}, \end{equation}

where $m_P$ and $m_T$ denote the mass of the OB part and the mass of the OC part with the axis O as the dividing line. Thus, $m_P$ and $m_T$ can be expressed as:

(5) \begin{equation} m_P=\frac{ml_P}{l_P+l_T}, \end{equation}

(6) \begin{equation} m_T=\frac{ml_T}{l_P+l_T}. \end{equation}

The simplified single-link dynamic model of the exoskeleton can be expressed as:

(7) \begin{equation} \tau _{hip}=m\ddot{\theta }(l_P+l_T)^2+\frac{m\ddot{\theta }}{3(l_P+l_T)}(l_P^3+l_T^3)+[\cos{\theta _0}-\cos{(\theta -\theta _0)}]\frac{mgl_P(l_P-l_T)}{2(l_P+l_T)}. \end{equation}

The moments produced by the Bowden cable pulling force and the elastic spring force on the joint are analyzed below. AB’, B’C’, and DE’ represent the current positions of the elastic spring, the exoskeleton thigh rod, and the Bowden cable when the PAMs are inflated and the Bowden cable actuates the exoskeleton to rotate at an angle $\theta$ . $F_P$ and $F_T$ in the red part of the figure represent the force transmitted by the Bowden cable and the force generated by the spring stretching, respectively. The red dashed lines represent the force arm $d_T$ of the force $F_P$ and the force arm $d_P$ of the spring force $F_T$ , respectively. It can be found that the force arm $d_T$ of the spring hardly changes during the rotation, so it is assumed that the force arm $d_T$ of the spring is constant and $d_T=l_T$ . The approximate length of the spring stretch is expressed as:

(8) \begin{equation} \Delta x=\theta l_T. \end{equation}

Then, the moment generated by the force of the spring can be expressed as:

(9) \begin{equation} \begin{aligned} \tau _T&= \Delta x k_T d_T\\[5pt] &=\theta k_T l_T^2. \end{aligned} \end{equation}

For the force arm $d_P$ of the Bowden cable, approximated here as $\textsf{OD}=\textsf{OE}$ , the force arm of the Bowden cable can be calculated as:

(10) \begin{equation} d_P=\frac{l_P}{2}\cos{\left(\frac{\theta _1-\theta }{2}\right)}, \end{equation}

where $\theta _1$ is the angle between OC and OD at the initial state.

The tension in the Bowden cable $F_P$ can be obtained through the three-element PAM model which has been extensively used in the previous research [Reference Chen, Huang, Wu and Tu18, Reference Huang, Cao, Xiong and Zhang19]. Therefore, $F_P$ can be expressed as:

(11) \begin{equation} F_P=F(P)-B(P)\Delta \dot{x}-K(P)\Delta x, \end{equation}

where $P$ is the internal air pressure of the PAM, $B(P)$ denotes the damping coefficient, $K(P)$ denotes the elasticity coefficient, and the values of $B(P)$ and $K(P)$ are related to the states of the PAM. $\Delta x$ denotes the displacement of the PAM, which is also the contraction length of the Bowden cable. Therefore, the relationship between $\Delta x$ and the rotation angle $\theta$ can be expressed as:

(12) \begin{equation} \Delta x=l_P\left(\sin{\frac{\theta _1}{2}} -\sin{\frac{\theta _1-\theta }{2}}\right), \end{equation}

and:

(13) \begin{equation} \Delta \dot{x}= \frac{1}{2}\dot{\theta }\cos{\frac{(\theta _1-\theta )}{2}}. \end{equation}

Substituting (12) and (13) into (11), the pulling force in the Bowden cable can be expressed as:

(14) \begin{equation} F_P=F(P)-B(P)\frac{1}{2}\dot{\theta }\cos{\frac{(\theta _1-\theta )}{2}}-K(P)l_P\left(\sin{\frac{\theta _1}{2}} -\sin{\frac{\theta _1-\theta }{2}}\right). \end{equation}

Therefore, the torque generated by the Bowden cable can be expressed as:

(15) \begin{equation} \begin{aligned} \tau _P&=F_Pd_P\\[5pt] &=\frac{l_P}{2}\cos\left(\frac{\theta _1-\theta }{2}\right)\left[F(P)-B(P)\frac{1}{2}\dot{\theta }\cos{\frac{(\theta _1-\theta )}{2}}-K(P)l_P\left(\sin{\frac{\theta _1}{2}} -\sin{\frac{\theta _1-\theta }{2}}\right)\right]. \end{aligned} \end{equation}

The total torque generated by the Bowden cable and the antagonistic elastic spring can be expressed as:

(16) \begin{equation} \tau _{hip}=\tau _P-\tau _T. \end{equation}

Substituting (7), (9) and (15) into (16) yields:

(17) \begin{equation} M\dot{x_2}=H_1(x_1,x_2)+H_2(x_1,x_2)u, \end{equation}

where

(18) \begin{equation} M=\left(\frac{4}{3}l_P^2+\frac{4}{3}l_T^2+\frac{5}{3}l_Pl_T\right)m, \end{equation}

(19) \begin{equation} \begin{aligned} H_1(x_1,x_2)&=\frac{l_P}{2}\cos{\frac{\theta _1-x_1}{2}}\left[F_0-\frac{B_0x_2}{2}\cos{\frac{\theta _1-x_1}{2}}-K_0l_P\left(\sin{\frac{\theta _1}{2}}-\sin{\frac{\theta _1-x_1}{2}}\right)\right]\\[5pt] &-\frac{mgl_P(l_P-l_T)}{2(l_P+l_T)}\left[\cos{\theta _0}-\cos{(x_1-\theta _0)}\right]-x_1k_Tl_T^2, \end{aligned} \end{equation}

(20) \begin{equation} H_2(x_1,x_2)=\frac{l_P}{2}\cos{\frac{\theta _1-x_1}{2}}\left[F_1-B_1\frac{1}{2}x_2\cos{\frac{\theta _1-x_1}{2}}-K_1l_P\left(\sin{\frac{\theta _1}{2}}-\sin{\frac{\theta _1-x_1}{2}}\right)\right]. \end{equation}

In summary, the dynamic model of the exoskeleton can be expressed as follows:

(21) \begin{equation} \left \{ \begin{aligned} \dot{x}_1 &=x_2\\[5pt] \dot{x}_2 &=f(\bar{x})+g(\bar{x})u(t)+d(t)\\[5pt] y &=x_1, \end{aligned} \right. \end{equation}

(22) \begin{equation} g(\bar{x})=H_1(x_1,x_2)/M, \end{equation}

(23) \begin{equation} f(\bar{x})=H_2(x_1,x_2)/M, \end{equation}

where $y$ denotes the output of the system, $x_1$ and $x_2$ are system states and $\bar{x}=[x_1,x_2]^T$ . Due to the different sizes of the subjects, the human-robot interaction, and the errors of the system parameter identification, there exist errors in the above constructed dynamic model, and $d(t)$ denotes a lumped term including the dynamic modeling uncertainties of the exoskeleton system and the external disturbances.

3.1. Preliminaries

Lemma 3.1. [Reference Li, Tong and Li20] For a continuous unknown function $f(x)$ defined on a compact set $\Omega$ , there exists a fuzzy system $\Theta ^T\Phi (x)$ such that for any given constant $\varepsilon \gt 0$ , the following inequality holds:

(24) \begin{equation} \mathop{sup}\limits _{x \in \Omega }|f(x)-\Theta ^T\Phi (x)|\leq \varepsilon. \end{equation}

Lemma 3.2. [Reference Wang, Chen, Lin, Zhang and Meng21] For $z_k\in \mathbb{R}$ , $k=1,2,\ldots,n$ , $0\lt l\lt 1$ , the following inequality holds:

(25) \begin{equation} \left(\sum _{k=1}^n|z_k|\right)^l\leq \sum _{k=1}^n|z_k|^l\leq n^{1-l}\left(\sum _{k=1}^n|z_k|\right)^l. \end{equation}

Lemma 3.3. [Reference Wang, Chen, Liu and Lin22, Reference Zhu, Xia and Fu23] Considering the system $\dot{\varsigma }=f(\varsigma,u)$ , for smooth positive definite function $V(\varsigma )$ , if there exist scalar $c\gt 0$ , $0\lt l\lt 1$ and $\varrho \gt 0$ , such that:

(26) \begin{equation} \dot{V}(\varsigma )\leq -cV^l(\varsigma )+\varrho,\;\;\;t\geq 0, \end{equation}

the nonlinear system $\dot{\varsigma }=f(\varsigma,u)$ is semi-global practical finite time stable (SGPFS).

3.2. The single-parameter adaptive fuzzy controller (SAFC) design

According to Lemma 3.1, the lumped term $d(t)$ can be expressed in the form of a fuzzy system as follows:

(27) \begin{equation} d(t)=\Theta ^{*T}\Phi (\bar{x})+\epsilon, \end{equation}

where $*$ means the ideal fuzzy consequence parameter value, $\bar{x}$ denotes the vector containing the system states, and $\epsilon$ is the approximation error.

The fuzzy system approximating the actual $d(t)$ is expressed as:

(28) \begin{equation} \hat{d}=\hat{\Theta }^T\Phi (\bar{x}), \end{equation}

where $\hat{}$ denotes the approximate fuzzy consequence parameter value.

We define the single adaptation parameter as:

(29) \begin{equation} \hat{\omega }=\hat{\Theta }^T\hat{\Theta }. \end{equation}

The approximate error can be expressed as:

(30) \begin{equation} \begin{aligned} \tilde{d}(t)&=d(t)-\hat{d}\\[5pt] &=\Theta ^{*T}\Phi (\bar{x})-\hat{\Theta }^T\Phi (\bar{x})+\epsilon \\[5pt] &=\tilde{\Theta }^T\Phi (\bar{x})+\epsilon. \end{aligned} \end{equation}

The main results of the controller are summarized in the following theorem.

Theorem 3.4. For the nonlinear system (21) with full state constraints, the lumped term $d(t)$ can be approximated and compensated by a fuzzy system with one adaptation parameter $\hat{\omega }$ , all signals of the close-loop system are SGPFS in a finite time with the following virtual control signal $\alpha _1$ , controller $u(t)$ and adaptive law $\dot{\hat{\omega }}$ :

(31) \begin{equation} \alpha _1=\dot{y}_r-\frac{e_1^{2l-1}}{(k_{b1}^2-e_1^2)^{l-1}}-\frac{e_1}{2(k_{b1}^2-e_1^2)}, \end{equation}

(32) \begin{equation} u(t)=\frac{1}{g(\bar{x})}[-f(\bar{x})+\dot{\alpha }_1-\frac{e_2^{2l-1}}{(k_{b2}^2-e_2^2)^{l-1}}-\frac{e_2}{2(k_{b2}^2-e_2^2)}-\frac{\hat{\omega }\Phi (\bar{x})^T\Phi (\bar{x})e_2}{2l_1^2(k_{b2}^2-e_2^2)}-\frac{e_2}{2}(k_{b2}^2-e_2^2)], \end{equation}

(33) \begin{equation} \dot{\hat{\omega }}=\gamma [\frac{\Phi (\bar{x})^T\Phi (\bar{x})e_2^2}{2l_1^2(k_{b2}^2-e_2^2)^2}-l_2\hat{\omega }], \end{equation}

where $0\lt l\lt 1$ , $l_1\gt 0$ , $l_2\gt 0$ are parameters to be designed, and $k_{b1}\gt 0$ and $k_{b2}\gt 0$ are constraints of the errors $e_1$ and $e_2$ .

Proof. We define the error terms as follows:

(34) \begin{equation} e_1=x_1-y_r, \end{equation}

(35) \begin{equation} e_2=x_2-\alpha _1, \end{equation}

where $\alpha _1$ denotes the virtual control signal of $x_2$ .

Taking the derivative of $e_1$ yields:

(36) \begin{equation} \begin{aligned} \dot{e}_1&=x_2-\dot{y}_r. \end{aligned} \end{equation}

The first Lyapunov function with a barrier function disposing of the constraint of $e_1$ is designed as:

(37) \begin{equation} V_1=\frac{1}{2}\log{\frac{k^2_{b1}}{k^2_{b1}-e_1^2}}, \end{equation}

where $\log$ denotes the nature logarithm function, $k_{b1}\gt 0$ denotes the constraint of $e_1$ .

Taking the derivative of $V_1$ yields:

(38) \begin{equation} \begin{aligned} \dot{V}_1&=\frac{e_1(\dot{y}_r-\dot{x}_1)}{k_{b1}^2-e_1^2}\\[5pt] &=\frac{e_1(\dot{y}_r-x_2)}{k_{b1}^2-e_1^2}. \end{aligned} \end{equation}

Substituting (35) into (38) yields:

(39) \begin{equation} \dot{V}_1=\frac{e_1(\dot{y}_r-e_2-\alpha _1)}{k_{b1}^2-e_1^2}. \end{equation}

Substituting the virtual signal $\alpha _1$ into (39) yields:

(40) \begin{equation} \begin{aligned} \dot{V}_1=\frac{e_1e_2}{k_{b1}^2-e_1^2}-\left[\frac{e_1^2}{k_{b1}^2-e_1^2}\right]^l-\frac{e_1^2}{2(k_{b1}^2-e_1^2)^2}. \end{aligned} \end{equation}

According to Young’s inequality, we have:

(41) \begin{equation} \frac{e_1e_2}{k_{b1}^2-e_1^2}\leq \frac{e_1^2}{2(k_{b1}^2-e_1^2)^2}+\frac{1}{2}e_2^2. \end{equation}

Substituting (41) into (40), we get:

(42) \begin{equation} \dot{V}_1\leq -\left[\frac{e_1^2}{k_{b1}^2-e_1^2}\right]^l+ \frac{1}{2}e_2^2. \end{equation}

It is obvious that $\log \frac{k_{b1}^2}{k_{b1}^2-e_1^2}\leq \frac{e_1^2}{k_{b1}^2-e_1^2}$ when $|e_1|\lt k_{b1}$ , which proves that $[\log \frac{k_{b1}^2}{k_{b1}^2-e_1^2}]^l\leq [\frac{e_1^2}{k_{b1}^2-e_1^2}]^l$ holds for $0\lt l\lt 1$ . Then, we can obtain:

(43) \begin{equation} \dot{V}_1\leq -V_1^l+\frac{1}{2}e_2^2. \end{equation}

We design another Lyapunov function as:

(44) \begin{equation} V_2=\log{\frac{k_{b2}^2}{k_{b2}^2-e_2^2}}+\frac{1}{2\gamma }\tilde{\omega }^T\tilde{\omega }, \end{equation}

where $\omega ^*=\Theta ^{*T}\Theta ^*$ , $\tilde{\omega }=\omega ^*-\hat{\omega }$ , $\gamma \gt 0$ , and $k_{b2}\gt 0$ denotes the constraint of $e_2$ .

Taking the derivative of $V_2$ and substituting $\dot{x}_2$ into the result yield:

(45) \begin{equation} \begin{aligned} \dot{V}_2&=\frac{e_2\dot{e}_2}{k_{b2}^2-e_2^2}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}\\[5pt] &=\frac{e_2(\dot{x}_2-\dot{\alpha }_1)}{k_{b2}^2-e_2^2}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}\\[5pt] &=\frac{e_2(f(\bar{x})+g(\bar{x})u(t)+d(t)-\dot{\alpha }_1)}{k_{b2}^2-e_2^2}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}. \end{aligned} \end{equation}

Substituting the control signal (32) and (27) into (45) yields:

(46) \begin{equation} \begin{aligned} \dot{V}_2&=\frac{e_2}{k_{b2}^2-e_2^2}\Bigg\{f(\bar{x})+\Bigg[-f(\bar{x})+\dot{\alpha }_1+\frac{-e_2^{2l-1}}{(k_{b2}^2-e_2^2)^{l-1}}-\frac{e_2}{2(k_{b2}^2-e_2^2)}-\frac{\hat{\omega }\Phi (\bar{x})^T\Phi (\bar{x})e_2}{2l_1^2(k_{b2}^2-e_2^2)}\\[5pt] &-\frac{e_2}{2}(k_{b2}^2-e_2^2)\Bigg]+\Theta ^{*T}\Phi (\bar{x})+\epsilon -\dot{\alpha }_1)\Bigg\}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}\\[5pt] &=\frac{e_2}{k_{b2}^2-e_2^2}\Bigg\{\frac{-e_2^{2l-1}}{(k_{b2}^2-e_2^2)^{l-1}}-\frac{e_2}{2(k_{b2}^2-e_2^2)}-\frac{\hat{\omega }\Phi (\bar{x})^T\Phi (\bar{x})e_2}{2l_1^2(k_{b2}^2-e_2^2)}-\frac{e_2}{2}(k_{b2}^2-e_2^2)+\Theta ^{*T}\Phi (\bar{x})\\[5pt] &+\epsilon \Bigg\}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}\\[5pt] &=-\left[\frac{e_2^2}{k_{b2}^2-e_2^2}\right]^l-\frac{e_2^2}{2(k_{b2}^2-e_2^2)^2}-\frac{\hat{\omega }\Phi (\bar{x})^T\Phi (\bar{x})e_2^2}{2l_1^2(k_{b2}^2-e_2^2)^2}-\frac{e_2^2}{2}+\frac{e_2\Theta ^{*T}\Phi (\bar{x})+e_2\epsilon }{k_{b2}^2-e_2^2}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}. \end{aligned} \end{equation}

With Young’s inequality, we have:

(47) \begin{equation} \frac{e_2\Theta ^{*T}\Phi (\bar{x})}{k_{b2}^2-e_2^2}\leq \frac{\omega ^*\Phi (\bar{x})^T\Phi (\bar{x})e_2^2}{2l_1^2(k_{b2}^2-e_2^2)^2}+\frac{l_1^2}{2}, \end{equation}

(48) \begin{equation} \frac{e_2\epsilon }{k_{b2}^2-e_2^2}\leq \frac{e_2^2}{2(k_{b2}^2-e_2^2)^2}+\frac{\epsilon ^2}{2}. \end{equation}

Then, we can obtain:

(49) \begin{equation} \begin{aligned} \dot{V}_2\leq &-\left[\frac{e_2^2}{k_{b2}^2-e_2^2}\right]^l+\frac{\tilde{\omega }\Phi (\bar{x})^T\Phi (\bar{x})e_2^2}{2l_1^2(k_{b2}^2-e_2^2)^2}+\frac{l_1^2}{2}+\frac{\epsilon ^2}{2}-\frac{e_2^2}{2}-\frac{1}{\gamma }\tilde{\omega }\dot{\hat{\omega }}. \end{aligned} \end{equation}

Substituting the adaptive law (33) into (49) yields:

(50) \begin{equation} \dot{V}_2\leq -\left[\frac{e_2^2}{k_{b2}^2-e_2^2}\right]^l+\frac{l_1^2}{2}+\frac{\epsilon ^2}{2}-\frac{e_2^2}{2}-l_2\hat{\omega }\tilde{\omega }. \end{equation}

Using the relationship of $\hat{\omega }=\omega ^*-\tilde{\omega }$ and Young’s inequality, we have:

(51) \begin{equation} \tilde{\omega }\hat{\omega }=\tilde{\omega }(\omega ^*-\tilde{\omega })=-\tilde{\omega }^2+\tilde{\omega }\omega ^*\leq -\frac{1}{2}\tilde{\omega }^2+\frac{1}{2}\omega ^{*2}. \end{equation}

Then, we can obtain:

(52) \begin{equation} \begin{aligned} \dot{V}_2&\leq -\left[\frac{e_2^2}{k_{b2}^2-e_2^2}\right]^l+\frac{1}{2}\omega ^{*2}-\frac{1}{2}\tilde{\omega }^2+\frac{1}{2}l_1^2+\frac{1}{2}\epsilon ^2-\frac{1}{2}e_2^2\\[5pt] &\leq -V_2^l+\frac{1}{2}\omega ^{*2}+\frac{1}{2}l_1^2+\frac{1}{2}\epsilon ^2-\frac{1}{2}e_2^2. \end{aligned} \end{equation}

Define the final Lyapunov function as:

(53) \begin{equation} V=V_1+V_2. \end{equation}

Taking the derivative of $V$ and substituting (43) and (52) into it, we get:

(54) \begin{equation} \begin{aligned} \dot{V}&=\dot{V}_1+\dot{V}_2\\[5pt] &\leq -V_1^l-V_2^l+\frac{1}{2}\omega ^{*2}+\frac{1}{2}l_1^2+\frac{1}{2}\epsilon ^2. \end{aligned} \end{equation}

Applying Lemma 3.2, we obtain:

(55) \begin{equation} \begin{aligned} \dot{V}&\leq -(V_1+V_2)^l+\frac{1}{2}\omega ^{*2}+\frac{1}{2}l_1^2+\frac{1}{2}\epsilon ^2\\[5pt] &\leq -V^l+\rho, \end{aligned} \end{equation}

where $\rho =\frac{1}{2}\omega ^{*2}+\frac{1}{2}l_1^2+\frac{1}{2}\epsilon ^2$ is a positive constant.

Define $T^*$ as:

(56) \begin{equation} T^*=\frac{1}{\eta (1-l)}\left[V^{1-l}(0)-\left(\frac{\rho }{(1-\eta )^{\frac{1-l}{l}}}\right)\right], \end{equation}

where $0\lt \eta \lt 1$ , $V(0)$ is the initial value of $V(t)$ . According to Lemma 3.3, for $\forall t\gt T^*$ , $V^l\leq \frac{\rho }{(1-\eta )*1}$ , that is, all the signals of the closed-loop system are SGPFS. Furthermore, based on the definition of $V(t)$ in ref. [Reference Li, Zhao, He and Lu16], for $\forall t\gt T^*$ , we have:

(57) \begin{equation} |y-y_r|\leq k_{b1}\left[1-e^{-2\left(\frac{\rho }{(1-\eta )}\right)^{\frac{1}{l}}}\right]^{\frac{1}{2}}, \end{equation}

which implies that after the finite time $T^*$ , the tracking error remain in a small neighborhood of the origin. We assume that $|y_r|\lt B_0$ ,and thus, it is obvious that $|x_1|\leq |e_1|+|y_r|\lt k_{b1}+B_0$ . We define that $k_{b1}=k_{c1}-B_0$ , and then, it has $|x_1|\lt k_{c1}$ . Similarly, with the assumption that $\dot{y}_r\lt B_1$ , we can finally get $|x_2|\lt k_{c2}$ and $k_{c2}=k_{b2}+B_1$ . Consequently, the system states are confined to the predefined regions effectively. This completes the proof.

5.1. Experiment platform setup

The practical exoskeleton experiment platform is represented in Fig. 9. The main components of the system include an air compressor, an xPC target, a host PC, PAM actuators, Bowden cables, a treadmill, a healthy subject wearing the exoskeleton, and a body weight support system compensating for the subject’s weight. In addition, two pairs of airbags are placed on the front and back sides of the exoskeleton leg straps part to detect the subject’s intention to move forward or backward. The mechanism and control system details have been explained in Section 2.

Figure 9. The practical exoskeleton system.

5.2. Trajectory tracking experiments with no-load

First of all, we conduct a set of experiments with no-load to validate the feasibility of the proposed SAFC method in practical implementation. The no-load experiments can protect the subjects from terrible control results. Then, the proposed SAFC method is applied to the exoskeleton control in healthy subjects’ passive and active gait training tasks, and the supplementary video file shows the experimental results. The sampling time of the real-time system is set as $0.0005\;{\rm s}$ . The desired trajectory $y_r$ of the exoskeleton is obtained by fitting the human lower hip motion data with $1.0m/s$ walking speed on the treadmill.

It is worth noting that since the proposed SAFC in this paper and the standard AFC in the above numerical simulation show similar control performances, and the proposed SAFC achieves better computational efficiency, the standard AFC is not considered in the experimental comparison. A conventional PID controller commonly utilized in practical engineering applications is adopted to compare with the SAFC to further investigate its control performance.

The designed parameters of the proposed SAFC are set as: $k_{b1}=5$ , $b_{b2}=300.0$ , $l=0.5$ , $l_1=10.0$ , $l_2=1.0$ , $\gamma =10.0$ . The initial adaptation parameter $\omega _0=10.0$ . The parameters of the PID controller are set as: $P=1050.0$ , $I=305.0$ , $D=10.0$ . For a fair comparison, we have tried to adjust the control parameters so that the controllers can perform their best in experiments.

The tracking trajectories are represented in Fig. 10, and the corresponding tracking errors are shown in Fig. 11. Only the tracking results of the right hip joint are presented here, since the left and right limbs of the exoskeleton conduct symmetrical movements. We can find that the tracking trajectory of the proposed SAFC nearly overlapped with the desired trajectory in Fig. 10, while that of PID deviates a lot from the desired trajectory, especially in the areas of peaks and bottoms of waves. From Fig. 11, we can find that the maximum tracking error of the proposed SAFC is less than $5\;{\rm deg}$ , which is bounded by the predefined boundary $k_{b1}$ , while that of PID is more than $5\;{\rm deg}$ and the maximum error is about $8\;{\rm deg}$ . Obviously, the proposed SAFC performs better than the conventional PID controller in the practical experiments with no-load. Therefore, we can conclude that the proposed SAFC method can precisely track the desired gait trajectory of the pneumatic lower limb exoskeleton to conduct the following gait training tasks, which indicates that this theoretically safe method guaranteeing the full state constraints can ensure training security of the exoskeleton in practical implementation.

Figure 10. Trajectories of experiments with no-load.

Figure 11. Tracking errors of experiments with no-load.

5.3. Passive gait training

From Figs. 10 and 11, we can find that the tracking trajectory of PID variated a lot from the desired trajectory and the maximum tracking error of PID was larger than 5 deg, which indicated that the PID was unable to constraint the tracking error to a predefined range and achieve safety control compared to the proposed SAFC. Therefore, the PID was not adopted for the poor tracking performance in the following gait tracking experiments.

To further validate the effectiveness of the proposed SAFC method in the task of passive gait training, four healthy volunteer subjects are invited to participate in our experiments. The subjects are asked to remain relaxed and walk on the treadmill with as much assistance as possible from the exoskeleton, thus achieving the effect of simulating hemiplegic patients. It is worth noting that here we consider the influence of the subject’s participation on the exoskeleton mechanism and modify the constraint $k_{b1}$ to $8$ , thus allowing a slightly larger tracking error than the no-load experiments for the successful implementation of the gait training experiments. In addition, we use the maximum absolute error(MAE) and the root mean square error (RMSE) in $10\sim 30\;s$ to evaluate the control performance of the proposed method. The RMSE is defined as:

(59) \begin{equation} RMSE=\sqrt{\frac{1}{n}\sum _{t=1}^{n}e_t^2}, \end{equation}

where $e_t$ denotes the tracking error at the $t$ th sampling time point.

Table II. The information of the healthy subjects and the gait training performances.

Figure 12. Trajectories of passive gait training.

The information of these four subjects and the details of the tracking errors are represented in Table II. The tracking results of the exoskeleton in the passive gait training are represented in Figs. 12 and 13. It can be clearly found that the tracking error is large in the first $10\;{\rm s}$ , while the tracking error can be maintained within $5\;{\rm deg}$ in the following $20\;{\rm s}$ , which indicates that after two gait periods of adjustment, the subjects can adapt well to the exoskeleton assistance to complete the passive gait training. Therefore, with the control results represented above, we can also conclude that although the four subjects vary in their heights and weights, the proposed SAFC method can still achieve effective tracking control to conduct the passive gait training task.

Figure 13. Tracking errors of passive gait training.

5.4. Active gait training

For the rehabilitation of partial patients, increasing the patient’s participation in the rehabilitation will be more beneficial than purely passive gait training [Reference Brahmi, Saad, Luna, Archambault and Rahman25–Reference Pehlivan, Sergi, Erwin, Yozbatiran, Francisco and O’Malley27]. Therefore, in this section, we conduct an active gait training experiment on one healthy subject to verify the feasibility of the proposed pneumatic lower limb exoskeleton with the SAFC method. The active control diagram is based on a zero-force control strategy, and the corresponding control diagram is represented in Fig. 14. The air bags detect the subject’s movement intention, the admittance model generates the desired trajectory $\Delta x_{ref}$ based on the interaction force $F_{int}$ and the desired force $F_d$ , and the proposed SAFC method is used to control the joint angle of the exoskeleton to achieve active control. Here we assume that the interaction force $F_{int}$ between the air bags and the subject’s leg and the air pressure difference between the air bag pairs are proportional, and the proportionality gain is $\kappa$ . We set $F_d=0$ and $x_{ref}=0$ so that the exoskeleton is completely passive and the subject is completely active. The parameters of the admittance model are set as: $m=1/500$ , $b=0.05$ , and $k=0.01$ . The subject walks on the treadmill at a speed of $1.0\,{\rm m/s}$ .

Figure 14. The block diagram of the active gait training control strategy.

The pressure difference of the air bags is represented in Fig. 15, and the trajectories are represented in Fig. 16. From Fig. 15, we can find that the air pressure difference of the air bags varied periodically, indicating that the subject’s intention of movement can be successfully detected by the air bags positioned in both the front and back of the subject’s legs. From Fig. 16, we can find that the trajectory generated by the active control strategy shown in the blue curve is close in magnitude and frequency to the reference trajectory in the passive gait training mentioned previously, which indicates that the admittance model can successfully generate the desired reference trajectory used in active training. The red curve shows the movement trajectory of the subject wearing the exoskeleton under the active control strategy. It can be found that the exoskeleton with the SAFC method can effectively track the trajectory generated by the admittance model to achieve the active gait training task.

Figure 15. Pressure difference of the air bags.

Figure 16. Trajectories of active gait training.