Robust adaptive PD-like control of lower limb rehabilitation robot based on human movement data

Analysis of the mechanism of human lower limb movement

Human lower limb movement data collection

The existing motion capture technology is mainly divided into four categories: mechanical motion capture system (Ivanov & AZhilenkova, 2019), acoustic motion capture system, electromagnetic motion capture system, and optical motion capture system (Ren et al., 2018; Pfister et al., 2014). Through the comparison of motion capture systems, this paper adopts the better NOKOV 3D infrared passive optical motion capture system. This system work principle is to capture the movement trajectory of the reflective Marker point attached to the human body in space by using passive infrared light principle. The 3D force measurement platform system mainly realizes the mechanical force measurement and analysis functions of the human body during the movement process. During the measurement, the force in the x, y and z axes are integrated into multiple channel data through 3D force sensors, and then collected by NI Card. Before collecting experimental data, site layout is essential. The cameras of the motion capture system are usually arranged according to the experimental requirements (Kurihara et al., 2002; Zeng et al., 2016). Because the experiment done in this paper is to collect human lower limb movement data, cameras are placed in a spac with 5 × 4 × 3 m. The cameras are placed in three equal intervals in each column on the left and right, and the 3D force measurement platforms are placed in the center, as shown in Fig. 2.

The NK_Cortex software can observe the movement state of the Marker points attached to the human body in real time. When the human body moves in the space of the motion capture system, as long as the Marker points at the joints can be captured by two or more cameras at the same time, the location information of the Marker points is recorded. Store the collected human movement data in real time and save it in required files format through the software. The quality of data collection plays an important role in the research. Therefore, after collecting the data, it is necessary to check whether the position, velocity, acceleration and other information of each Marker point are missing, and observe whether the curve is disconnected. If there is an obvious disconnection in the middle of the curve means that the data of the frame position is lost. If the lost data fluctuates greatly or the number of consecutive lost frames is large, data collection needs to be repeated, otherwise the data can be repaired.

Gait trajectory data processing

The collected human gait accuracy will affect the stability, coordination and driving torque of LLRR system. The ideal gait is the continuous and smooth motion trajectory of each joint, such as position, angle, angular velocity and angular acceleration. The 3D motion gait capture system mainly captures the spatial position coordinate information of the hip joint, knee joint and ankle joint. During the measurement process, NK_Cortex software can repair the curve disconnection problem, such as the missing coordinate data of some Marker points due to occlusion or other external factors. However, there are external factors such as marker jitter and environmental capture during human movement, spliced data needs to be filtered at this time. Filtering methods mainly include Least Squares filtering, Wavelet filtering, Kalman filtering, and others. In this paper, least square filtering and wavelet filtering are used to process the data of human left leg hip joint, knee joint and angle joint.

The kinematics of the robot is used to study the geometric relationship between the motion angle of each joint of the LLRR and its position in space. Similarly, the relationship between the spatial position coordinates of each joint of the human body's lower limb and the angle of each joint is found through the collected data, which plays a vital role in the control of the gait trajectory data and the joint angle of the LLRR in the later stage. Kinematics is divided into positive kinematics analysis method and inverse kinematics analysis method. Based on the needs of experimental research, the LLRR-like mode as shown in Fig. 3. this paper adopts the inverse kinematics analysis method, that is, using the position data of the end of the human lower limb hip joint $a\left(a_x,a_y\right)$, knee joint $b\left(b_x,b_y\right)$ and ankle joint $c\left(c_x,c_y\right)$, the hip joint angle $q_{d1}$ and knee joint angle $q_{d2}$ can be obtained through the tangent theorem,

(1) $$q_{\mathrm{d}1}=\arctan {\displaystyle \frac{b_{\mathrm{x}}-a_{\mathrm{x}}}{b_{\mathrm{y}}-a_{\mathrm{y}}}}$$

(2) $$\alpha =\arctan {\displaystyle \frac{c_{\mathrm{x}}-b_{\mathrm{x}}}{c_{\mathrm{y}}-b_{\mathrm{y}}}}$$

(3) $$q_{\mathrm{d}2}=q_{\mathrm{d}1}-\alpha =q_{\mathrm{d}1}-\arctan {\displaystyle \frac{c_{\mathrm{x}}-b_{\mathrm{x}}}{c_{\mathrm{y}}-b_{\mathrm{y}}}}$$From the results obtained from the above two formulas and analyzing the normal gait of the human body during walking, it can be obtained $q_{\mathrm{d}1}\in \left[\begin{array}{cc}0& {\displaystyle \frac{\pi }{2}}\end{array}\right]$, $q_{\mathrm{d}2}\in \left[\begin{array}{cc}0& {\displaystyle \frac{\pi }{2}}\end{array}\right]$. Therefore, the inverse kinematics equations of the hip and knee joints of the human body can be obtained as, (4) $$\{\begin{array}{c}{q}_{\mathrm{d}1}=\arctan {\displaystyle \frac{b_{\mathrm{x}}-a_{\mathrm{x}}}{b_{\mathrm{y}}-a_{\mathrm{y}}}}\\ q_{\mathrm{d}2}=q_{\mathrm{d}1}-\arctan {\displaystyle \frac{c_{\mathrm{x}}-b_{\mathrm{x}}}{c_{\mathrm{y}}-b_{\mathrm{y}}}}\end{array}$$Its movement trajectory is shown in Fig. 4.

Force data processing

Since the rhythm of the human body during walking is not fixed, under-sampling or over-sampling will occur when the movement is too fast or too slow in the process of force data acquisition (Xia et al., 2017; Zhou et al., 2020). Over-sampling will result in excessive information collection, which increases the burden of data analysis, and low utilization. Under-sampling will cause some essential information lost (Zhang & Cao, 2015). Therefore, this paper uses machine learning methods to process the collected 100 sets of force data. In order to deal with data over-sampling, the cubic spline interpolation method is used for processing, and the final processed data is 128 frames per second.

The 3D force measurement platform system collects the 3D force data of the end position of the experimenter’s lower limbs. This paper research is based on a two-dimensional plane, so the force in the x and y directions are selected, as shown in Fig. 4. And the force component of the position c can be calculated is shown in Fig. 5.

The force in the x direction is, (5) $$M_x=f_x(c_y-a_y)$$The force in the y direction is, (6) $$M_y=f_y\left(c_x-a_x\right)$$where, $f_x$, $f_y$ respectively present the component on the x-axis and y-axis during real-time change of the collected experimentation three-dimensional plantar force.

The Jacobian matrix $J$ of LLRR system establishes a mapping from joint speed to operating, which can better analyze the speed of system. Similarly, $J^T$ represents the force Jacobian matrix, which means a mapping linear relationship between operating force and joint force of LLRR system in the static equilibrium state.

The torque of the hip and knee joint positions can be obtained by the following equations, and its changes are shown in Fig. 6,

(7) $$\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\end{array}\right]=J^{\mathrm{T}}\left[\begin{array}{c}{M}_{\mathrm{x}}\\ M_{\mathrm{y}}\end{array}\right]$$Where

(8) $$J^{\mathrm{T}}=\left[\begin{array}{cc}{l}_{11}\cos q_{\mathrm{d}1}+l_{22}\cos \left(q_{\mathrm{d}1}-q_{\mathrm{d}2}\right)& -l_{11}\sin q_{\mathrm{d}1}-l_{22}\sin \left(q_{\mathrm{d}1}-q_{\mathrm{d}2}\right)\\ -l_{22}\cos \left(q_{\mathrm{d}1}-q_{\mathrm{d}2}\right)& \sin \left(q_{\mathrm{d}1}-q_{\mathrm{d}2}\right)\end{array}\right]$$Therefore, the experimentalist’s hip joint torque is, (9) $$\begin{array}{l}{\tau }_h=\left(l_{11}\cos q_{d1}+l_{22}\cos \left(q_{d1}-q_{d2}\right)\right)f_x\left(c_y-a_y\right)\\ +\left(-l_{11}\sin q_{d2}-l_{22}\sin \left(q_{d1}-q_{d2}\right)\right)f_y\left(c_x-a_x\right)\end{array}$$the experimentalist’s knee joint torque is, (10) $${\tau }_k=-l_{22}\cos \left(q_{d1}-q_{d2}\right)f_x\left(c_y-a_y\right)+\sin \left(q_{d1}-q_{d2}\right)f_y\left(c_x-a_x\right)$$where, $l_{11}$ is the length of experimenter’s thigh and $l_{11}=0.4m$. $l_{22}$ is the length of experimenter’s calf and $l_{22}=0.38m$.

Analysis of human movement mechanism

Human walking movement plays a vital role in daily life. The study of human gait parameters has a significant meaning for understanding human movement laws and limbs during walking coordination (Jackson et al., 2004). During walking, the time interval between two consecutive occurrences of the same action in a repetitive motion event is called the gait cycle. Generally speaking, the lower limbs of human body walking cycle is composed of support phase and swing phase (Barbareschi et al., 2015). As with most research results, human trajectory data demonstrate during a normal walking gait cycle, the support phase accounts for about 60% and the swing phase accounts for about 40%. And the proportion of support phase and swing phase in the gait cycle also changes with time.

Human biomechanics is committed to study bone strength, muscle strength, joint range, and various parts speed of the human body. The purpose of research on human biomechanics is to make human-computer interaction play a better role for patients, and avoid doing useless work. Ensure that people wearing LLRR can receive effective rehabilitation treatment and personal safety, and considering the safely and comfort of patients during rehabilitation training. The controller design must be based on the human biomechanics theories, such as the size of the human body's output, the trajectory of the limbs, the smoothness of the movement, and the movement direction of human body parts. In walking process, the clockwise moment required by human joints are found though analysis inverse dynamic. In collected dynamics data, the ground reaction force is the main force that acts on the body during experimenter movement. According to the research needs of this paper, the vertical ground reaction force and the shear force in one direction are studied.

The typical vertical ground reaction force profile of the collected single human walking steps, as shown in Fig. 5. When the heel is in contact with the ground, the vertical ground reaction force is zero and rises sharply within a few tenths of a second. At the moment experimenter’s foot is flat, her weight will dropp and move onto legs. In the middle of standing, the experimenter’s movement center of mass is actually upward. Finally, the experimenter’s feet lost contact with the ground and the force returns to zero when the toes left the ground. It can be seen than the experimenter’s foot force change is a typical M-shaped diagram. Numerous studies have shown that the M-shaped diagram or bimodal diagram are a typical display method of normal gait, which is used to display force fluctuations. The above analysis of the collected experimenter movement data is effective and feasible in this paper.

Dynamic model of LLRR

In this paper, the unilateral leg of LLRR is simplified as a two-bar rigid structure model, as shown in Fig. 7. The Lagrange equation of the LLRR system can be defined as the following form,

(11) $${\tau }_1={\displaystyle \frac{d}{dt}{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{\dot{q}}_1}-{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{q}_1},i=1,2}}}$$where, $L$ is the Lagrange function of the LLRR system in the generalized coordinate system, which expressed by the difference between the total kinetic energy and the total potential energy of the LLRR system. $q_i$ represent the joint angle, ${\dot{q}}_i$ represent the joint angle velocity, ${\tau }_i$represent the joint torque.

Suppose the thigh length of the LLRR is $l_1$, its mass is $m_1$, the calf length is $l_2$, its mass is $m_2$, the center of mass of the thigh rod is $c_1$, the center of mass of the calf rod is $c_2$, the distance from the hip joint to its center of mass is $d_1$, and the distance from the knee joint to its center of mass is $d_2$. The angle of the hip joint and the angle of the knee joint of the LLRR are $q_1$, $q_2$, the joints angular velocity of the LLRR are ${\dot{q}}_1$, ${\dot{q}}_2$, and the joint angular acceleration are ${\ddot{q}}_1$, ${\ddot{q}}_2$, respectively.

The joint torque of the LLRR system can be expressed as, $${\tau }_1={\displaystyle \frac{d}{dt}{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{\dot{q}}_1}-{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{q}_1}}}}$$

$$=\left(m_1{d}_1^2+m_2{l}_1^2+m_2{d}_2^2+2m_2{l}_1{d}_2\cos q_2\right){\ddot{q}}_1-\left(m_2{l}_1{d}_2\cos q_2+m_2{d}_2^2\right){\ddot{q}}_2$$

$$-2m_2{l}_1{d}_2\sin {\dot{q}}_2{q}_2{\dot{q}}_1+m_2{l}_1{d}_2{\dot{q}}_2\sin {\dot{q}}_2{q}_2-m_{2}g{d}_2\sin \left(q_1-q_2\right)$$

(12) $$-m_{1}g{d}_{1}sin{q}_1-m_{2}g{l}_{1}sin{q}_1$$

$${\tau }_2={\displaystyle \frac{d}{dt}{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{\dot{q}}_2}-{\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{q}_2}}}}$$

(13) $$=-m_2{d}_2\left(l_1\cos q_2+d_2\right){\ddot{q}}_1+m_2{l}_1{d}_2\sin {\dot{q}}_2{q}_1+m_2{d}_2^2{\ddot{q}}_2+m_{2}g{d}_2\sin \left(q_1-q_2\right)$$In order to simplify the expression of its equation, make

(14) $$D\left(q\right)=\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right]=\left[\begin{array}{cc}{m}_1{d}_1^2+m_2{l}_1^2+m_2{d}_2^2+2m_2{l}_1{d}_2\cos q_2& -\left(m_2{l}_1{d}_2\cos q_2+m_2{d}_2^2\right)\\ -m_2{d}_2\left(l_1\cos q_2+d_2\right)& m_2{d}_2^2\end{array}\right]$$

(15) $$C\left(q\right)=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]=\left[\begin{array}{cc}-2m_2{l}_1{d}_2{\dot{q}}_2\sin q_2& m_2{l}_1{d}_2{\dot{q}}_2\sin q_2\\ m_2{l}_1{d}_2{\dot{q}}_1\sin q_2& 0\end{array}\right]$$

(16) $$G\left(q\right)=\left[\begin{array}{c}{G}_1\\ G_2\end{array}\right]=\left[\begin{array}{c}-m_{1}g{d}_1\sin q_1-m_{2}g{l}_1\sin q_1-m_{2}g{d}_2\sin (q_1-q_2)\\ m_{2}g{d}_2\sin (q_1-q_2)\end{array}\right]$$Equations (12) and (13) can be changed to the following equation, (17) $$D\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=\tau$$where, $\tau ={\left[\begin{array}{cc}{\tau }_1& {\tau }_2\end{array}\right]}^{\mathrm{T}}$, ${\tau }_1\left(t\right)$ and ${\tau }_2\left(t\right)$ represent the hip joint torque and the knee joint torque of the controller. $D\left(q\right)\in R^{2\times 2}$ represents the inertia matrix, $C\left(q,\dot{q}\right)\in R^{2\times 2}$ represents the Coriolis force and the centripetal force matrix, $G\left(q\right)\in R^{2\times 1}$ represents gravity matrix.

The dynamic characteristics of the LLRR shown in formula Eq. (17) are,

1. $\dot{D}\left(\theta \right)-2C\left(\theta ,\dot{\theta }\right)$ is an oblique symmetric matrix.

2. The inertial matrix $D\left(q\right)$ is a $2\times 2$ symmetric positive definite matrix of. For any $2\times 1$ vector $\zeta$, there are positive numbers $m_1$ and $m_2$ satisfying the following inequalities, $$m_1\Vert \zeta {\Vert }^2\le {\zeta }^{T}D\left(q\right)\zeta \le m_2\Vert \zeta {\Vert }^2$$

3. There is a parameter vector $P$ related to the lower limb exoskeleton robot, so that the inertia matrix $D\left(q\right)$, Coriolis force matrix $C\left(q\right)$, and gravity terms $G\left(q\right)$ satisfy the following linear relationships

(19) $$D\left(q\right)\vartheta +C\left(q,\dot{q}\right)\rho +G\left(q\right)=\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P$$

In the formula, $\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)\in R^{2\times 5}$ is the regression matrix of the joint variables of the known LLRR. $P\in R^{5\times 1}$ is an unknown fixed-length parameter vector that expresses the characteristics of the robot's mass and inertia.

The RAPLC research on of the LLRR system in this paper, since its controller is designed based on an unknown model, it is necessary to perform online identification of the unknown model. The dynamic characteristics of the robot system are used to separate the known and unknown parameters of the robot model in the identification process. Make

(20) $$p_1=m_1{d}_1^2+m_2{l}_1^2+m_1{d}_2^2$$

(21) $$p_2=m_1{l}_1{d}_2$$

(22) $$p_3=m_2{d}_2^2$$

(23) $$p_4=-\left(m_1{d}_1+m_2{l}_1\right)g$$

(24) $$p_5=m_2{d}_{2}g$$The Eqs. (14)–(16) can be expressed as, (25) $$D\left(q\right)=\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right]=\left[\begin{array}{cc}{p}_1+2p_2\cos q_2& -\left(p_3+p_2\cos q_2\right)\\ -(p_3+p_2\cos q_2)& p_3\end{array}\right]$$

(26) $$C\left(q,\dot{q}\right)=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]=\left[\begin{array}{cc}-2p_2{\dot{q}}_2\sin q_2& p_2{\dot{q}}_2\sin q_2\\ p_2{\dot{q}}_1\sin q_2& 0\end{array}\right]$$

(27) $$G\left(q\right)=\left[\begin{array}{c}{G}_1\\ G_2\end{array}\right]=\left[\begin{array}{c}{p}_4\sin q_1-p_5\sin \left(q_1-q_2\right)\\ p_5\sin \left(q_1-q_2\right)\end{array}\right]$$Substituting Eqs. (25)–(27) into Eq. (19) can be obtained, (28) $$\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right]\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\end{array}\right]+\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\end{array}\right]+\left[\begin{array}{c}{G}_1\\ G_2\end{array}\right]=\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P$$Then the dynamic equation of the LLRR thigh shaft can be expressed as, $$\left[\begin{array}{cc}{D}_{11}& D_{12}\end{array}\right]\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\end{array}\right]+\left[\begin{array}{cc}{C}_{11}& C_{12}\end{array}\right]\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\end{array}\right]+G_1$$

$$=\left(p_1+2p_2\cos q_2\right){\vartheta }_1-\left(p_3+p_2\cos q_2\right){\vartheta }_2-2p_2{\dot{q}}_2\sin q_2{\rho }_1+p_2{\dot{q}}_2\sin q_2{\rho }_2$$

$$+p_4\sin q_1-p_5\sin (q_1-q_2)$$

$$={\vartheta }_1{p}_1+\left(2{\vartheta }_1\cos q_2-{\vartheta }_2\cos q_2{\vartheta }_2\cos q_2-2{\dot{q}}_2\sin q_2{\rho }_1+{\dot{q}}_2\sin q_2{\rho }_2\right)p_2-{\vartheta }_2{p}_3$$

$$+\sin q_1{p}_4-\sin \left(q_1-q_2\right)p_5$$

(29) $$=\left[\begin{array}{ccccc}{\mathrm{\Phi }}_{11}& {\mathrm{\Phi }}_{12}& {\mathrm{\Phi }}_{13}& {\mathrm{\Phi }}_{14}& {\mathrm{\Phi }}_{15}\end{array}\right]\left[\begin{array}{c}{p}_1\\ p_2\\ p_3\\ p_4\\ p_5\end{array}\right]$$where, $${\mathrm{\Phi }}_{11}={\vartheta }_1$$

$${\mathrm{\Phi }}_{12}=2{\vartheta }_1\cos q_2-{\vartheta }_2\cos q_2{\vartheta }_2\cos q_2-2{\dot{q}}_2\sin q_2{\rho }_1+{\dot{q}}_2\sin q_2{\rho }_2$$

(30) $${\mathrm{\Phi }}_{13}=-{\vartheta }_2$$

$${\mathrm{\Phi }}_{14}=\sin q_1$$

$${\mathrm{\Phi }}_{15}=-\sin \left(q_1-q_2\right)$$Similarly, the dynamic equation of the LLRR calf rod is, $$\left[\begin{array}{cc}{D}_{21}& D_{22}\end{array}\right]\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\end{array}\right]+\left[\begin{array}{cc}{C}_{21}& C_{22}\end{array}\right]\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\end{array}\right]+G_2$$

$$=-\left(p_3+p_2\cos q_2){\vartheta }_1+p_3{\vartheta }_2+p_2{\dot{q}}_1\sin q_2{\rho }_1+p_5\sin (q_1-q_2\right)$$

$$=\left({\dot{q}}_1\sin q_2{\rho }_1-{\vartheta }_1\cos q_2)p_2+\left({\vartheta }_2-{\vartheta }_1\right)p_3+\sin (q_1-q_2\right)p_5$$

(31) $$=\left[\begin{array}{ccccc}{\mathrm{\Phi }}_{21}& {\mathrm{\Phi }}_{22}& {\mathrm{\Phi }}_{23}& {\mathrm{\Phi }}_{24}& {\mathrm{\Phi }}_{25}\end{array}\right]\left[\begin{array}{c}{p}_1\\ p_2\\ p_3\\ p_4\\ p_5\end{array}\right]$$where, $${\mathrm{\Phi }}_{21}=0$$

$${\mathrm{\Phi }}_{22}={\dot{q}}_1\sin q_2{\rho }_1-{\vartheta }_1\cos q_2$$

(32) $${\mathrm{\Phi }}_{23}={\vartheta }_2-{\vartheta }_1$$

$${\mathrm{\Phi }}_{24}=0$$

$${\mathrm{\Phi }}_{25}=\sin (q_1-q_2)$$Derived from the above formula, the parameters in $\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)$ are only related to the joint angle and angular acceleration of the LLRR. These data are all known, that is, $\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)$ is the known parameter of the LLRR system model. $\vartheta$ represents the equation related to acceleration of the end position of the LLRR system. The parameters in $P$ are related to the joint quality, joint length, center of gravity and other data of the LLRR system. However, these parameters in the robot system are often inaccurate or completely unknown, especially the LLRR system. The joint quality is constantly changing real-time when the patient wear it for rehabilitation, so the parameters in $P$ are the unknown parameters.

Assumption 1: $q_d\in R^n$ is the desired position of joint, $q_d$ has first derivative and second derivative.

Assumption 2: The norms of error and disturbance satisfy

$$\Vert \omega \Vert \le z_1+z_2\Vert e\Vert +z_3\Vert \dot{e}\Vert$$where, $z_1$, $z_2$, $z_3$ are normal numbers, respectively. $e=q-q_d$ and $\dot{e}=\dot{q}-{\dot{q}}_d$ are tracking error and tracking error derivative, respectively.

Controller design

In the rehabilitation process, the comfort of human body wearing LLRR plays a significant role, so this paper proposes a comfort function, (33) $${\tau }_h=m\dot{e}+ne=m\left({\dot{q}}_d-\dot{q}\right)+n\left(q_d-q\right)$$where, $m$ and $n$ are positive number. $q$ and $\dot{q}$ are the angular displacement and angular velocity of the controller, respectively. $q_d$ and ${\dot{q}}_d$are the angular displacement and angular velocity of the experimenter, respectively.

The unilateral leg dynamics model of the two-link LLRR can be described by Eq. (17), namely, (34) $$D\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=\tau$$Introduce variables $y$ and $q_r$, and make, (35) $$y=\dot{e}+\gamma e$$

(36) $$q_r={\dot{q}}_d-\gamma e$$In the formula $\gamma >0$, and

(37) $$y=\dot{q}-{\dot{q}}_r$$In formula Eq. (19) make $\vartheta ={\ddot{q}}_r$, $\rho ={\dot{q}}_r$, (38) $$D\left(q\right){\ddot{q}}_r+C\left(q,\dot{q}\right){\dot{q}}_r+G\left(q\right)=\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P$$where, (39) $${\ddot{q}}_r=\ddot{q}-\dot{y}=\left[\begin{array}{c}{\ddot{q}}_1-\dot{y}\\ {\ddot{q}}_2-\dot{y}\end{array}\right]=\left[\begin{array}{c}{\ddot{q}}_r\left(1\right)\\ {\ddot{q}}_r\left(2\right)\end{array}\right]$$

(40) $${\dot{q}}_r=q-y=\left[\begin{array}{c}{\dot{q}}_1-y\\ {\dot{q}}_2-y\end{array}\right]=\left[\begin{array}{c}{\dot{q}}_r\left(1\right)\\ {\dot{q}}_r\left(2\right)\end{array}\right]$$Substitute formula Eq. (35) into formula Eq. (38), (41) $$D\left(q\right)(\ddot{q}-\dot{y})+C\left(q,\dot{q}\right)\left(\dot{q}-y\right)+G\left(q\right)=\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P$$Combining formula Eq. (41) and formula Eq. (17) can be obtained, (42) $$D\left(q\right)\dot{y}+C\left(q,\dot{q}\right)y=\tau -\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P-\omega$$For a complex LLRR system, due to the influence of external factors, some errors and deviations will inevitably occur. In order to ensure the LLRR system stability when the upper bound of the disturbance signal is known, the following RAPLC is designed, (43) $$\tau =\mathrm{\Phi }\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)\hat{P}-S_p{K}_{p}e-S_v{K}_v\dot{e}+u$$

(44) $$D\left(q\right)\dot{y}+C\left(q,\dot{q}\right)y=\tau -\mathrm{\Phi }\left(q,\dot{q},\rho ,\vartheta \right)P-\omega$$In the formula Eq. (43), $$K_p=K_{p1}+K_{p2}{B}_p\left(e\right),K_v=K_{v1}+K_{v2}{B}_v\left(\dot{e}\right)$$

$$K_{\mathrm{p}1}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(k_{\mathrm{p}11},k_{\mathrm{p}12}\right),K_{\mathrm{p}2}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(k_{\mathrm{p}21},k_{\mathrm{p}22}\right)$$

$$K_{\mathrm{v}1}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(k_{\mathrm{v}11},k_{\mathrm{v}12}\right),K_{\mathrm{v}2}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(k_{\mathrm{v}21},k_{\mathrm{v}22}\right)$$

$$B_{\mathrm{p}}\left(e\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\displaystyle \frac{1}{{\alpha }_1+\left|e_1\right|},{\displaystyle \frac{1}{{\alpha }_2+\left|e_2\right|}}}\right),B_{\mathrm{v}}\left(e\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\displaystyle \frac{1}{{\beta }_1+\left|{\dot{e}}_1\right|},{\displaystyle \frac{1}{{\beta }_2+\left|{\dot{e}}_2\right|}}}\right)$$where, ${\gamma }_1,{\gamma }_2$ are arbitrarily normal, $k_{\mathrm{p}1\mathrm{i}},k_{\mathrm{p}2\mathrm{i}},k_{\mathrm{v}1\mathrm{i}},k_{\mathrm{v}2\mathrm{i}},{\alpha }_{\mathrm{i}},{\beta }_{\mathrm{i}}$ are positive numbers in $i\in \left[1,2\right]$, and satisfy $k_{\mathrm{p}1\mathrm{i}}-{\displaystyle \frac{k_{\mathrm{v}2\mathrm{i}}}{2{\beta }_{\mathrm{i}}}>0}$, $k_{\mathrm{v}1\mathrm{i}}-{\displaystyle \frac{k_{\mathrm{p}2\mathrm{i}}}{2{\alpha }_{\mathrm{i}}}>0}$. $S_{\mathrm{p}}$ and $S_{\mathrm{v}}$ are S-shaped function matrix, and,

$S_{\mathrm{p}}=\left[\begin{array}{cc}{\displaystyle \frac{c_1}{1+e^{-{\mathrm{a}}_1\mathrm{t}}}}& 0\\ 0& {\displaystyle \frac{c_2}{1+e^{-{\mathrm{a}}_2\mathrm{t}}}}\end{array}\right]$, $S_{\mathrm{v}}=\left[\begin{array}{cc}{\displaystyle \frac{1}{1+e^{-{\mathrm{b}}_1\mathrm{t}}}}& 0\\ 0& {\displaystyle \frac{1}{1+e^{-{\mathrm{b}}_2\mathrm{t}}}}\end{array}\right]$

where, $a_1,a_2,b_1,b_2$ are the decay factor of the exponential function, $c_1$and $c_2$ are the coefficients of the exponential function.

The parameter estimation law of $\hat{P}$, (45) $$\stackrel{.}{\stackrel{⌢}{P}}=-{\text{ΓΦ}}^T\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)y$$where, $\mathrm{\Gamma }$ is a symmetric positive definite matrix. For the LLRR system shown in Eq. (17), and when the disturbance error signal is bounded, the designed control method can ensure the global gradual stability of the system.

Stability analysis

Make $\tilde{P}=\dot{P}-P$, and take the Lyapunov function of LLRR as, (46) $$V={\displaystyle \frac{1}{2}\left(y^{\mathrm{T}}D\left(q\right)y+e^{\mathrm{T}}\left(S_{\mathrm{p}}{K}_{\mathrm{p}1}+\gamma S_{\mathrm{v}}{K}_{\mathrm{v}1}\right)e+{\tilde{P}}^{\mathrm{T}}{\mathrm{\Gamma }}^{-1}\tilde{P}\right)}$$

According to the dynamic characteristics of the LLRR, $D\left(q\right)$ is a symmetric positive definite matrix, $S_{\mathrm{p}}$ and $S_{\mathrm{v}}$ are also symmetric positive definite matrix. So, (47) $${\left(y^{\mathrm{T}}D\left(q\right)y\right)}^{{}^{\prime }}=2y^{\mathrm{T}}D\left(q\right)\dot{y}+y^{\mathrm{T}}\dot{D}\left(q\right)y$$

(48) $${\left(e^{\mathrm{T}}\left(S_{\mathrm{p}}{K}_{\mathrm{p}1}+\gamma S_{\mathrm{v}}{K}_{\mathrm{v}1}\right)e\right)}^{{}^{\prime }}=2e^{\mathrm{T}}\left(S_{\mathrm{p}}{K}_{\mathrm{p}1}+\gamma S_{\mathrm{v}}{K}_{\mathrm{v}1}\right)\dot{e}$$

(49) $${\left({\tilde{P}}^T{\mathrm{\Gamma }}^{-1}\tilde{P}\right)}^{{}^{\prime }}=2{\tilde{P}}^T{\mathrm{\Gamma }}^{-1}\dot{\tilde{P}}$$Take the derivative of Eq. (46) to get, $$\dot{V}=y^{\mathrm{T}}D\left(q\right)\dot{y}+{\displaystyle \frac{1}{2}{y}^{\mathrm{T}}\stackrel{\dot{y}}{D\left(q\right)}{e}^{\mathrm{T}}+e^{\mathrm{T}}\left(S_{\mathrm{p}}{K}_{\mathrm{p}1}+\gamma S_{\mathrm{v}}{K}_{\mathrm{v}1}\right)\dot{e}+{\tilde{P}}^{\mathrm{T}}{\mathrm{\Gamma }}^{-1}\dot{\tilde{P}}}$$

$$=-{\dot{e}}^T{S}_v{K}_v\dot{e}-\gamma e^T{S}_p{K}_{p}e-e^T\left(S_p{K}_{p2}{B}_p\left(e\right)+\gamma S_v{K}_{v2}{B}_p\left(\dot{e}\right)\right)\dot{e}$$

(50) $$+{\tilde{P}}^T{\mathrm{\Gamma }}^{-1}\stackrel{\dot{+}}{\tilde{P}}{y}^T\mathrm{\Phi }\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)\tilde{P}+y^{\mathrm{T}}\left(u-\omega \right)$$Combine $y^T\mathrm{\Phi }\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)\tilde{P}={\tilde{P}}^T{\mathrm{\Phi }}^T\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)y$, $\dot{\tilde{P}}=\dot{\hat{P}}$ and Eq. (45), (51) $$y^T\mathrm{\Phi }\left(q,\dot{q},{\dot{q}}_r,{\ddot{q}}_r\right)\tilde{P}+{\tilde{P}}^T{\mathrm{\Gamma }}^{-1}\dot{\tilde{P}}=0$$Then, $$\dot{V}=-{\dot{e}}^{\mathrm{T}}{S}_{\mathrm{v}}{K}_{\mathrm{v}}\dot{e}-\gamma e^{\mathrm{T}}{S}_{\mathrm{p}}{K}_{\mathrm{p}}e-e^{\mathrm{T}}\left(S_{\mathrm{p}}{K}_{\mathrm{p}2}{B}_{\mathrm{p}}\left(e\right)+\gamma S_{\mathrm{v}}{K}_{\mathrm{v}2}{B}_{\mathrm{p}}\left(\dot{e}\right)\right)\dot{e}+y^{\mathrm{T}}\left(u-\omega \right)$$

$$=-\sum _{i=1}^2\left(\gamma \left(S_p{k}_{p1i}+{\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}}\right){e_i}^2\right)-\sum _{i=1}^2\left(\left(S_v{k}_{v1i}+{\displaystyle \frac{S_v{k}_{v2i}}{{\beta }_i+\left|\dot{e_i}\right|}}\right){\dot{e_i}}^2\right)$$

(52) $$-\sum _{i=1}^2\left(\left({\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}+{\displaystyle \frac{\gamma S_v{k}_{v2i}}{{\beta }_i+\left|\dot{e_i}\right|}}}\right)e_i\dot{e_i}\right)+y^{\mathrm{T}}\left(u-\omega \right)$$

The following formula can be obtained due to $e_i{\dot{e}}_i\le {\displaystyle \frac{1}{2}\left({e_i}^2+{{\dot{e}}_i}^2\right)}$, $$\dot{V}\le -\sum _{i=1}^2\left(\gamma \left(S_p{k}_{p1i}+{\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}}\right)-{\displaystyle \frac{1}{2}\left({\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}+{\displaystyle \frac{\gamma S_v{k}_{v2i}}{{\beta }_i+\left|{\dot{e}}_i\right|}}}\right)}\right){e_i}^2$$

(53) $$-\sum _{i=1}^2\left(\left(S_v{k}_{v1i}+{\displaystyle \frac{S_v{k}_{v2i}}{{\beta }_i+\left|\dot{e_i}\right|}}\right)-{\displaystyle \frac{1}{2}\left({\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}+{\displaystyle \frac{\gamma S_v{k}_{v2i}}{{\beta }_i+\left|\dot{e_i}\right|}}}\right){\dot{e_i}}^2}\right)+y^T\left(u-\omega \right)$$When ${\displaystyle \frac{1}{2}\le \gamma \le 2}$, $$\gamma {\displaystyle \frac{S_{\mathrm{p}}{k}_{\mathrm{p}2\mathrm{i}}}{{\alpha }_{\mathrm{i}}+\left|e_{\mathrm{i}}\right|}-{\displaystyle \frac{1}{2}{\displaystyle \frac{S_{\mathrm{p}}{k}_{\mathrm{p}2\mathrm{i}}}{{\alpha }_{\mathrm{i}}+\left|e_{\mathrm{i}}\right|}\ge 0}}}$$

(54) $${\displaystyle \frac{S_{\mathrm{v}}{k}_{\mathrm{v}2\mathrm{i}}}{{\beta }_{\mathrm{i}}+\left|{\dot{e}}_{\mathrm{i}}\right|}-{\displaystyle \frac{1}{2}{\displaystyle \frac{\gamma S_{\mathrm{v}}{k}_{\mathrm{v}2\mathrm{i}}}{{\beta }_{\mathrm{i}}+\left|{\dot{e}}_{\mathrm{i}}\right|}\ge 0}}}$$Formula Eq. (55) can be expressed, $$\dot{V}\le -\sum _{i=1}^2\left(\gamma \left(S_p{k}_{p1i}-{\displaystyle \frac{1}{2}{\displaystyle \frac{S_v{k}_{v2i}}{{\beta }_i+\left|{\dot{e}}_i\right|}}}\right){e_i}^2+\left(S_v{k}_{v1i}-{\displaystyle \frac{1}{2}{\displaystyle \frac{S_p{k}_{p2i}}{{\alpha }_i+|e_i|}}}\right)\right)+y^T\left(u-\omega \right)$$

$$\le -\sum _{\mathrm{i}=1}^2\left(\gamma \left(S_{\mathrm{p}}{k}_{\mathrm{p}1\mathrm{i}}-{\displaystyle \frac{1}{2}{\displaystyle \frac{S_{\mathrm{v}}{k}_{\mathrm{v}2\mathrm{i}}}{{\beta }_i}}}\right){e_{\mathrm{i}}}^2+\left(S_{\mathrm{v}}{k}_{\mathrm{v}1\mathrm{i}}-{\displaystyle \frac{1}{2}{\displaystyle \frac{S_{\mathrm{p}}{k}_{\mathrm{p}2\mathrm{i}}}{{\alpha }_{\mathrm{i}}}}}\right)\right)+y^{\mathrm{T}}\left(u-\omega \right)$$

(55) $$\le y^{\mathrm{T}}\left(u-\omega \right)$$Since, $$y^{\mathrm{T}}u=\sum _{\mathrm{i}=1}^2{y}_{\mathrm{i}}\left[-\left(z_1+z_{2}e+z_3\dot{e}\right)\mathrm{s}\mathrm{g}\mathrm{n}\left(y_{\mathrm{i}}\right)\right]$$

$$\le \sum _{\mathrm{i}=1}^2-\left(||\omega ||.||y_{\mathrm{i}}||\right)\mathrm{s}\mathrm{g}\mathrm{n}\left(y_{\mathrm{i}}\right)-y^{\mathrm{T}}\omega$$

(56) $$\le ||y^{\mathrm{T}}||.||\omega ||$$Then, (57) $$y^{\mathrm{T}}\left(u-\omega \right)\le \sum _{\mathrm{i}=1}^2\left(-||\omega ||\cdot ||y_{\mathrm{i}}||\right)+||y^{\mathrm{T}}||\cdot ||\omega ||=0$$Namely, (58) $$\dot{V}\le 0$$It can be seen from the above proof results, $\underset{t\to \mathrm{\infty }}{lim}e=0$, $\underset{t\to \mathrm{\infty }}{lim}\dot{e}=0$, the LLRR control system is globally asymptotically stable. Therefore, the LLRR system can realize the desired trajectory tracking no matter where it starts from the initial position.

Comfort evaluation function

When the human body wears the LLRR, too much fluctuation of the angle of the system will cause discomfort to the wearer. The comfort of patients is usually known through subjective feelings, but most patients with lower limb movement disorders are not sensitive. Most studies usually evaluate changes in angle, joint torque, changes in patient muscle strength, etc. This paper uses the similarity formula Eq. (59) between the trajectory of the lower limb exoskeleton robot and the trajectory of the human body to objectively evaluate the comfort of the human body wearing LLRR.

(59) $$C={\displaystyle \frac{1}{1+{({\displaystyle \frac{q_1-q_{\mathrm{d}1}}{q_{1\mathrm{m}\mathrm{a}\mathrm{x}}-q_{\mathrm{d}1\mathrm{m}\mathrm{i}\mathrm{n}}})}}^2+{({\displaystyle \frac{q_2-q_{\mathrm{d}2}}{q_{2\mathrm{m}\mathrm{a}\mathrm{x}}-q_{\mathrm{d}2\mathrm{m}\mathrm{i}\mathrm{n}}})}}^2}}$$where, $q_1$and $q_{\mathrm{d}1}$ are LLRR hip joint and experimenter hip joint. $q_2$ and $q_{\mathrm{d}2}$ are LLRR knee joint and experimenter knee joint. $q_{1\mathrm{m}\mathrm{a}\mathrm{x}}$ and $q_{2\mathrm{m}\mathrm{a}\mathrm{x}}$ are the maximum values of LLRR hip and knee joints. $q_{\mathrm{d}1\mathrm{m}\mathrm{i}\mathrm{n}}$ and $q_{\mathrm{d}2\mathrm{m}\mathrm{i}\mathrm{n}}$ are the minimum values of the hip and knee joints of experimenter. When C is close to 1, it means that LLRR trajectory is very similar to the trajectory of the human body, which indirectly reflects that the comfort of the human body wearing LLRR will be higher.