2. Design of prototype system

Three-dimensional six-degree-of-freedom motion of human lumbar spine—three-dimensional refers to three motion axes (coronal axis, sagittal axis, vertical axis), six-degree-of-freedom motion refers to three angular displacements and three linear displacements. This paper is mainly to achieve the lumbar rehabilitation training of patients through the wire-driven moving platform. As shown in Fig. 1, it is a three-dimensional model of the wire-driven lumbar rehabilitation training parallel robot.

Figure 1. 3D model of WWRTPR.

In this paper, the wire-driven lumbar rehabilitation training parallel robot is mainly composed of brace, wire-driven moving platform, wire-driven system, control system, weight reduction system, lower limb fixed adjustment device, universal pulley, lifting platform, etc., to realize the rehabilitation training of patients with three-dimensional six-degree-of-freedom movement.

The brace is made of 1.5 and 2 m square aluminum profiles, respectively. The wire-driven system is mainly composed of four groups of servo motors, ball screw and moving slider, etc., which are used to drive four wires, so as to change the length and direction of the wire through the pulley. The weight reduction system is mainly composed of two groups of servo motor, ball screw, fixed pulley, weight reduction spring, load-bearing vest, etc. the wire is driven by the servo motor, so as to stretch the spring to generate tension on the human body and achieve the purpose of weight reduction during rehabilitation training. The lower limb fixed adjustment device is fixed on the lifting platform, which is mainly composed of thigh adjustment module and calf adjustment module. The belt of thigh and calf is adjusted according to the patient’s weight, and the rise and fall of the lifting platform are adjusted according to the patient’s height, so as to ensure the safety and comfort of patients in rehabilitation training. The designed wire-driven rehabilitation training parallel robot can provide different contact forces to meet the use requirements of different patients according to the needs of patients [Reference Zi, Yin and Zhang14].

3. Kinematics modeling analysis

In the rehabilitation training operation of WWRTPR, the motion trajectory planning of rehabilitation training is essentially an inverse kinematics problem. The inverse kinematics modeling is established according to the relationship between the wire length and the pose of the end effector (patient), that is, the variation of the wire length causes the adjustment of the patient’s pose. The kinematic model of the system is established to lay a theoretical foundation for the dynamic modeling. The structural schematic diagram of WWRTPR is shown in Fig. 2.

Figure 2. Structural schematic diagram of WWRTPR.

In Fig. 2, the static coordinate system (OXYZ) is fixed on the lifting platform and the moving coordinate system (PXYZ) is established at the waist of the patient, and the initial direction of the coordinate axes of the two coordinate systems is the same. Four traction wires (wire 1, wire 2, wire 3, wire 4) are used to drive patients to achieve spatial movement of six-degree-of-freedom, and the movement of the four traction wires is independent of each other. Two traction wires (O ₁ I ₁, O ₂ I ₂) and springs are used for weight reduction.

Figure 3. Kinematic diagram of WWRTPR.

The kinematic diagram of the system is shown in Fig. 3. Where, set ${\boldsymbol{P}}_{i}=\overline{\overline{{\mathrm{O}P_{i}}}}, {\boldsymbol{B}}_{i}=\overline{\overline{{\mathrm{O}B_{i}}}}$ , in the static coordinate system, the wire length vector ${\boldsymbol{L}}_{i}$ satisfies

(1) \begin{equation} {\boldsymbol{L}}_{i}={\boldsymbol{X}}_{P}+{\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}-{\boldsymbol{B}}_{i} \end{equation}

where ${\boldsymbol{X}}_{P} (X_{P},Y_{P},Z_{P})^{\mathrm{T}}$ is the coordinate of the origin P of the moving coordinate system in the static coordinate system; ${\boldsymbol{C}}_{0} (x_{{P_{i}}},y_{{P_{i}}},z_{{P_{i}}})^{\mathrm{T}}$ is the coordinate of P _i in the moving coordinate system Pxyz (i = 1,…, 4, i is the value of the wire); ${\boldsymbol{R}}$ is the rotation transformation matrix from the moving coordinate system to the static coordinate system, which can be expressed as

(2) \begin{equation} {\boldsymbol{R}}=\left[\begin{array}{ccc} \mathrm{c}\theta \mathrm{c}\psi &\quad \mathrm{s}\mathit{\phi }\mathrm{s}\theta \mathrm{c}\psi -\mathrm{c}\phi \mathrm{s}\psi &\quad \mathrm{c}\mathit{\phi }\mathrm{s}\theta \mathrm{c}\psi +\mathrm{s}\phi \mathrm{s}\psi \\ \mathrm{c}\theta \mathrm{s}\psi &\quad \mathrm{s}\mathit{\phi }\mathrm{s}\theta \mathrm{s}\psi +\mathrm{c}\phi \mathrm{c}\psi &\quad \mathrm{c}\mathit{\phi }\mathrm{s}\theta \mathrm{s}\psi -\mathrm{s}\phi \mathrm{c}\psi \\ -\mathrm{s}\theta &\quad \mathrm{s}\phi \mathrm{c}\theta &\quad \mathrm{c}\phi \mathrm{c}\theta \end{array}\right] \end{equation}

where $\phi, \theta, \psi$ , respectively, represent the forward and backward swing angle, left and right swing angle and roll angle of the patient’s waist rotating around the OX axis, OY axis and OZ axis in the static coordinate system (OXYZ).

Set the length of the ith wire as L _i (i = 1, …, 4), and its expression is

(3) \begin{equation} L_{i}=\sqrt{\left({\boldsymbol{X}}_{P}+{\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}-{\boldsymbol{B}}_{i}\right)^{\mathrm{T }}\!\left({\boldsymbol{X}}_{P}+{\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}-{\boldsymbol{B}}_{i}\right)} \end{equation}

where ${\boldsymbol{B}}_{i}(X_{{B_{i}}},Y_{{B_{i}}},Z_{{B_{i}}})^{\mathrm{T}}$ is the coordinate of each hinge point B _i in the static coordinate system.

The time derivative of the rotation transformation matrix R can be obtained as follows:

(4) \begin{equation} \dot{{\boldsymbol{R}}}={\boldsymbol{R}}^{\mathrm{o}}{\boldsymbol{R}} \end{equation}

where ${\boldsymbol{R}}^{\mathrm{o}}=\left[\begin{array}{ccc} 0 &\quad -\omega _{Z} &\quad \omega _{Y}\\ \omega _{Z} &\quad 0 &\quad -\omega _{X}\\ -\omega _{Y} &\quad \omega _{X} &\quad 0 \end{array}\right], \boldsymbol{\omega }$ is the angular velocity, ${\boldsymbol{R}}^{\mathrm{o}}=\boldsymbol{\omega }^{\times }$ , ( $\boldsymbol{\omega }^{\times }$ is the antisymmetric matrix of $\boldsymbol{\omega }$ ).

By substituting Eq. (4) into Eq. (3) and simplifying, the following results can be obtained:

(5) \begin{equation} L_{i}\dot{L}_{i}=\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\!\left(\dot{{\boldsymbol{X}}}_{P}+{\boldsymbol{R}}^{\mathrm{o}}{\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}\right) \end{equation}

Let ${\boldsymbol{r}}_{i}={\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}$ , then,

(6) \begin{equation} L_{i}\dot{L}_{i}=\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\dot{{\boldsymbol{X}}}_{P}+\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\!\left(\boldsymbol{\omega }^{\times }{\boldsymbol{r}}_{i}\right) \end{equation}

Since $\dot{{\boldsymbol{r}}}_{i}={\boldsymbol{R}}^{\mathrm{o}}{\boldsymbol{R}}{\boldsymbol{x}}_{{P_{i}}}, {\boldsymbol{R}}^{\mathrm{o}}=\boldsymbol{\omega }^{\times }$ , then $\dot{{\boldsymbol{r}}}_{i}=\boldsymbol{\omega }^{\times }{\boldsymbol{r}}_{i}$ . Through the example of MATLAB, it is proved that:

(7) \begin{equation} \boldsymbol{\omega }^{\times }{\boldsymbol{r}}_{i}=-{{\boldsymbol{r}}_{i}}^{\times }\boldsymbol{\omega }\mathit{,}\,\,\boldsymbol{\omega }^{\times }{\boldsymbol{r}}_{i}=\boldsymbol{\omega }\times {\boldsymbol{r}}_{i} \end{equation}

Therefore,

(8) \begin{equation} L_{i}\dot{L}_{i}=\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\dot{{\boldsymbol{X}}}_{P}+\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\!\left(\boldsymbol{\omega }\times {\boldsymbol{r}}_{i}\right) \end{equation}

Through sorting out Eq. (8), it can be concluded that:

(9) \begin{equation} L_{i}\dot{L}_{i}=\left({\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\dot{{\boldsymbol{X}}}_{P}+\left({\boldsymbol{r}}_{i}\times {\boldsymbol{L}}_{i}\right)^{\mathrm{T}}\boldsymbol{\omega } \end{equation}

When Eq. (9) is written in the form of matrix, it can be obtained:

(10) \begin{equation} \tilde{{\boldsymbol{L}}}\dot{{\boldsymbol{L}}}=\tilde{{\boldsymbol{J}}}_{A}\dot{{\boldsymbol{X}}}_{\omega } \end{equation}

The matrices in Eq. (10) are successively defined as

\begin{align*} &\dot{{\boldsymbol{X}}}_{\omega }=\left[\begin{array}{c} \dot{{\boldsymbol{X}}}_{{\boldsymbol{P}}}\\ \boldsymbol{\omega } \end{array}\right]=\left[\begin{array}{l} {\boldsymbol{v}}\\ \boldsymbol{\omega } \end{array}\right]={\boldsymbol{G}}\dot{{\boldsymbol{X}}},\ \dot{{\boldsymbol{L}}}=\left[\dot{L}_{1}\quad \dot{L}_{2}\quad \dot{L}_{3}\quad \dot{L}_{4} \right]^{\mathrm{T}}={\boldsymbol{JG}}\dot{{\boldsymbol{X}}},\ \tilde{{\boldsymbol{L}}}=\mathrm{diag}\left[\begin{array}{cccc} L_{1} &\quad L_{2} &\quad L_{3} &\quad L_{4} \end{array}\right],\\& {\tilde{{\boldsymbol{J}}}}_{A}^{\mathrm{T}}=\left[\begin{array}{cccc} {\boldsymbol{L}}_{1} &\quad {\boldsymbol{L}}_{2} &\quad {\boldsymbol{L}}_{3} &\quad {\boldsymbol{L}}_{4}\\ {\boldsymbol{r}}_{1}\times {\boldsymbol{L}}_{1} &\quad {\boldsymbol{r}}_{2}\times {\boldsymbol{L}}_{2} &\quad {\boldsymbol{r}}_{3}\times {\boldsymbol{L}}_{3} &\quad {\boldsymbol{r}}_{4}\times {\boldsymbol{L}}_{4} \end{array}\right].\end{align*}

where $\dot{{\boldsymbol{X}}}_{\omega }$ is the motion velocity vector of the patient model, which includes the linear velocity ${\boldsymbol{v}}$ and angular velocity $\boldsymbol{\omega }$ of the patient model relative to the brace, ${\boldsymbol{J}}$ is the Jacobian matrix of $4\times 6, {\boldsymbol{G}}$ is the transformation matrix based on attitude angle, and ${\boldsymbol{G}}=\left[\begin{array}{cc} {\boldsymbol{I}}_{3\times 3} &\quad \boldsymbol{0}_{3\times 3}\\ \boldsymbol{0}_{3\times 3} &\quad {\boldsymbol{H}} \end{array}\right], {\boldsymbol{H}}=\left[\begin{array}{ccc} \mathrm{c}\theta \mathrm{c}\psi &\quad -\mathrm{s}\psi &\quad 0\\ \mathrm{c}\theta \mathrm{s}\psi &\quad \mathrm{c}\psi &\quad 0\\ -\mathrm{s}\theta & \quad 0 & \quad 1 \end{array}\right]$ .

4. Dynamic modeling analysis

The purpose of building dynamic model is to study the relationship between the mechanism motion and the force of WWRTPR. The main purpose of the research is to design the motion control, so as to realize the specific motion of the end effector (patient) of the parallel robot and achieve the optimal control goal or better dynamic performance.

The dynamic model of WWRTPR includes three parts: dynamic model of moving platform, elastic dynamic model of wire and dynamic model of driving system, and the dynamic model of the moving platform (the end effector) is the basis of motion control design. The dynamic symbol schematic diagram of defining the system is shown in Fig. 4.

Figure 4. Dynamic diagram of WWRTPR.

Since the mass of the wire is much smaller than the sum of the mass and load of the end effector (the patient’s moving platform), the mass and inertia of the driving wire are ignored, and the elastic model of the wire is not considered, so as to simplify the dynamic model of the system. The dynamic model of the whole system is composed of the dynamic model of the patient’s moving platform and the dynamic model of the driving system. The two dynamic models are combined to obtain the dynamic model of the whole system.

The Newton–Euler method is used to establish the dynamic equation of the wire-driven lumbar rehabilitation training parallel robot relative to the origin P of the moving coordinate system, the established dynamic equation is as follows:

(11) \begin{align} \begin{cases} m\ddot{{\boldsymbol{P}}}={\boldsymbol{f}}_{e}+\mathop{\sum}\limits_{i=1}^{4}\left({-}T_{i}{\boldsymbol{u}}_{i}\right)+m{\boldsymbol{g}}\\ {\boldsymbol{A}}_{G}\dot{\boldsymbol{\omega }}=\boldsymbol{\tau }_{e}+\mathop{\sum}\limits_{i=1}^{4}{\boldsymbol{r}}_{i}\times \left({-}T_{i}{\boldsymbol{u}}_{i}\right)-\boldsymbol{\omega }\times \left({\boldsymbol{A}}_{G}\boldsymbol{\omega }\right) \end{cases} \end{align}

where m is the mass of the end effector (patient), $\ddot{{\boldsymbol{P}}}=\left(\ddot{X}_{P},\ddot{Y}_{P},\ddot{Z}_{P}\right)^{\mathrm{T}}$ is the linear acceleration vector of the reference point p of the end effector, ${\boldsymbol{f}}_{e}$ and $\boldsymbol{\tau }_{e}$ are, respectively, the external forces and torques acting on the patient, ${\boldsymbol{u}}_{i}=\frac{{\boldsymbol{L}}_{i}}{L_{i}}$ is the tension vector of the wire, $T_{i}$ is the tension of the ith wire, and ${\boldsymbol{g}}$ is the gravity acceleration vector, ${{\boldsymbol{A}}}_{G}^{\ast }$ is the inertia matrix of the end effector with respect to the center of gravity, ${\boldsymbol{A}}_{G}\boldsymbol{=}{\boldsymbol{R}}{{\boldsymbol{A}}}_{G}^{\ast }{\boldsymbol{R}}^{\mathrm{T}}$ and $\dot{\boldsymbol{\omega }}$ is the angular acceleration vector.

According to the dynamic torque balance equation, the dynamic equation of the driving system is as follows:

(12) \begin{equation} {\boldsymbol{M}}_{0}\ddot{\boldsymbol{\theta }}_{m}+{\boldsymbol{C}}_{0}\dot{\boldsymbol{\theta }}_{m}+\boldsymbol{\tau }_{l}=\boldsymbol{\tau } \end{equation}

where ${\boldsymbol{M}}_{0}$ is the inertia matrix equivalent to the actuator; ${\boldsymbol{C}}_{0}$ is the viscous friction coefficient matrix equivalent to the actuator; $\boldsymbol{\theta }_{m}$ is the motor rotation angle; $\boldsymbol{\tau }_{l}=\mu {\boldsymbol{u}}$ is the load moment generated by the wire tension; $\boldsymbol{\tau }$ is the output torque vector of the driver.

The generalized force acting on the system is ${\boldsymbol{F}}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l@{\quad}l} F_{x} & F_{y} & F_{z} & M_{x} & M_{y} & M_{z} \end{array}\right]$ , then ${\boldsymbol{F}}$ has the following relationship with the matrix and the tension of the wire:

(13) \begin{equation} {\boldsymbol{u}}=({\boldsymbol{J}}^{\mathrm{T}})^{+}{\boldsymbol{F}} \end{equation}

where ${\boldsymbol{u}}=\left[u_{1} \quad u_{2} \quad u_{3} \quad u_{4}\right]^{\mathrm{T}}$ is the tension vector of four wires, and $({\boldsymbol{J}}^{\mathrm{T}})^{+}$ is the pseudo inverse matrix of ${\boldsymbol{J}}^{\mathrm{T}}$ [Reference Vashista, Khan and Agrawal15, Reference Qian, Zi and Ding16].

Equations (11) and (12), as well as the relationship among wire length ${\boldsymbol{L}}$ , motor rotation angle $\boldsymbol{\theta }_{m}$ and patient posture ${\boldsymbol{X}}$ , are used to construct the whole dynamic equation of the system (as shown in Eq. (14)), which lays a theoretical foundation for the subsequent studies.

(14) \begin{equation} \begin{array}{l} \left({\boldsymbol{M}}({\boldsymbol{X}})-\dfrac{1}{\mu ^{2}}\cdot {\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{JG}}\right)\!\ddot{{\boldsymbol{X}}}-\dfrac{1}{\mu ^{2}}\left({\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}\dot{{\boldsymbol{J}}}{\boldsymbol{G}}+{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{J}}\dot{{\boldsymbol{G}}}+{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{C}}_{0}{\boldsymbol{JG}}\right)\!\dot{{\boldsymbol{X}}}\\ =-\dfrac{1}{\mu }\cdot {\boldsymbol{J}}^{\mathrm{T}}\boldsymbol{\tau }+{\boldsymbol{w}}_{g}+{\boldsymbol{w}}_{e}-{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right) \end{array} \end{equation}

where ${\boldsymbol{M}}({\boldsymbol{X}})$ is the inertia matrix of the end effector (patient), ${\boldsymbol{X}}$ is the actual pose of the end effector (patient), $\mu$ is the transmission coefficient, ${\boldsymbol{J}}$ is the Jacobian matrix of the system, ${\boldsymbol{w}}_{g}$ is the gravity vector of the end effector, ${\boldsymbol{w}}_{e}=[f_{e};\tau _{e}]$ , and ${\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)$ is the nonlinear Coriolis centrifugal force matrix.

Figure 5. Structure schematic diagram of force/position hybrid control of WWRTPR.

5. Control law design and stability analysis

The research object of this article, the WWRTPR, is a complex strong coupling, multi-input and multi-output, and nonlinear time-varying system, which has strict requirements for flexibility, safety and comfort. Therefore, the design of its motion control strategy (considering force and pose) is extremely important.

5.1. Design of control law

The ideal motion trajectory of waist rehabilitation training is set as ${\boldsymbol{X}}_{d}$ , the motion tracking error is defined as ${\boldsymbol{e}}={\boldsymbol{X}}-{\boldsymbol{X}}_{d}$ , the velocity error is defined as $\dot{{\boldsymbol{e}}}=\dot{{\boldsymbol{X}}}-\dot{{\boldsymbol{X}}}_{d}$ and the acceleration error is defined as $\ddot{{\boldsymbol{e}}}=\ddot{{\boldsymbol{X}}}-\ddot{{\boldsymbol{X}}}_{d}$ . In order to improve the motion control accuracy of patients and the control performance of the system, referring to the design advantages of literature [Reference Gao, Wang, Wu and Lin17–Reference Li, Yang and Song24] and combining with the rehabilitation training characteristics of WWRTPR, an intelligent control method of force/position hybrid control based on RBF neural network is proposed, as shown in Eq. (15), and the structural schematic diagram of force/position hybrid control is designed, as shown in Fig. 5.

(15) \begin{equation} \boldsymbol{\tau }=\boldsymbol{\tau }_{1}+\boldsymbol{\tau }_{2} \end{equation}

where $\boldsymbol{\tau }$ is the total control law, $\boldsymbol{\tau }_{1}$ is the position control law, and $\boldsymbol{\tau }_{2}$ is the force control law. Among them,

(16) \begin{align} \boldsymbol{\tau }_{1}= & -\mu \!\left({\boldsymbol{J}}^{\mathrm{T}}\right)^{+}\left\{\left({\boldsymbol{M}}({\boldsymbol{X}})-\frac{1}{\mu^{2}}\cdot {\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{JG}}\right)\!\ddot{{\boldsymbol{X}}}_{d}+\left[-\frac{1}{\mu ^{2}}\!\left({\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}\dot{{\boldsymbol{J}}}{\boldsymbol{G}}+{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{J}}\dot{{\boldsymbol{G}}} +{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{C}}_{0}{\boldsymbol{JG}}\right)\right]\!\dot{{\boldsymbol{X}}}_{d} \right.\nonumber\\& +{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)+\boldsymbol{\lambda }-{\boldsymbol{K}}_{p}{\boldsymbol{e}}-{\boldsymbol{K}}_{d}\dot{{\boldsymbol{e}}}\bigg\} \end{align}

Let ${\boldsymbol{K}}_{1}\boldsymbol{=}{\boldsymbol{M}}({\boldsymbol{X}})-\dfrac{1}{\mu ^{2}}\cdot {\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{JG}}, {\boldsymbol{K}}_{2}=-\dfrac{1}{\mu ^{2}}\!\left({\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}\dot{{\boldsymbol{J}}}{\boldsymbol{G}}+{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{M}}_{0}{\boldsymbol{J}}\dot{{\boldsymbol{G}}}+{\boldsymbol{J}}^{\mathrm{T}}{\boldsymbol{C}}_{0}{\boldsymbol{JG}}\right), {\boldsymbol{K}}_{3}={\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)-{\boldsymbol{w}}_{g}-{\boldsymbol{w}}_{e}$ , then,

(17) \begin{equation} \boldsymbol{\tau }_{1}=-\mu \!\left({\boldsymbol{J}}^{\mathrm{T}}\right)^{+}\left({\boldsymbol{K}}_{1}\ddot{{\boldsymbol{X}}}_{d}+{\boldsymbol{K}}_{2}\dot{{\boldsymbol{X}}}_{d}+{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)+\boldsymbol{\lambda }-{\boldsymbol{K}}_{p}{\boldsymbol{e}}-{\boldsymbol{K}}_{d}\dot{{\boldsymbol{e}}}\right) \end{equation}

(18) \begin{equation} \boldsymbol{\tau }_{2}=-\mu \!\left({\boldsymbol{J}}^{\mathrm{T}}\right)^{+}\left({\boldsymbol{F}}_{d}-{\boldsymbol{F}}_{s}-{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)\right) \end{equation}

In Eqs. (17) and (18), ${\boldsymbol{u}}_{d}=({\boldsymbol{J}}^{\mathrm{T}})^{+}{\boldsymbol{F}}_{d}$ is the expected wire tension vector, ${\boldsymbol{u}}_{s}=({\boldsymbol{J}}^{\mathrm{T}})^{+}{\boldsymbol{F}}_{s}$ is the actual wire tension vector, ${\boldsymbol{e}}_{F}={\boldsymbol{F}}_{d}-{\boldsymbol{F}}_{s}, \boldsymbol{\lambda }$ is external interference, dynamic coupling and patient diversity, etc., ${\boldsymbol{K}}_{p}$ is the proportional gain matrix, and ${\boldsymbol{K}}_{d}$ is the differential gain matrix.

By substituting $\boldsymbol{\tau }$ into the whole dynamic equation (14) of the system, the error closed-loop equation of the system can be obtained:

(19) \begin{equation} {\boldsymbol{K}}_{1}\ddot{{\boldsymbol{e}}}+\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\!\dot{{\boldsymbol{e}}}+{\boldsymbol{K}}_{p}{\boldsymbol{e}}=\boldsymbol{\lambda }+{\boldsymbol{w}}_{g}+{\boldsymbol{w}}_{e}-{\boldsymbol{e}}_{F}-{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right) \end{equation}

where $({\boldsymbol{J}}^{\mathrm{T}})({\boldsymbol{J}}^{\mathrm{T}})^{+}={\boldsymbol{I}}_{6\times 6}, ({\boldsymbol{J}})^{+}{\boldsymbol{J}}={\boldsymbol{I}}_{6\times 6}$ .

Let ${\boldsymbol{u}}_{1}=-{\boldsymbol{N}}\!\left({\boldsymbol{X}},\dot{{\boldsymbol{X}}}\right)-{\boldsymbol{e}}_{F}, {\boldsymbol{u}}_{1}$ is the feedback control law; let $\Delta \delta =\boldsymbol{\lambda }+{\boldsymbol{w}}_{g}+{\boldsymbol{w}}_{e}, \Delta \delta$ is the uncertainty and external disturbance of the system, and RBF neural network is used to approximate compensation $\Delta \delta$ , the expression is as follows:

(20) \begin{equation} \Delta \delta ={{\boldsymbol{W}}_{f}}^{\ast }{\boldsymbol{h}}({\boldsymbol{x}})+\xi _{f} \end{equation}

where $\xi _{f}$ is the approximation error, ${{\boldsymbol{W}}_{f}}^{\ast }$ is the ideal neural network weight, and ${\boldsymbol{h}}({\boldsymbol{x}})$ is the RBF Gaussian function.

By defining ${\boldsymbol{x}}_{1}={\boldsymbol{e}}$ and ${\boldsymbol{x}}_{2}=\dot{{\boldsymbol{e}}}+\varepsilon {\boldsymbol{e}}$ , Eq. (19) is transformed into:

(21) \begin{equation} \begin{cases} \dot{{\boldsymbol{x}}}_{1}={\boldsymbol{x}}_{2}-\varepsilon {\boldsymbol{x}}_{1}\\ \begin{array}{l} {\boldsymbol{K}}_{1}\dot{{\boldsymbol{x}}}_{2}=-\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right){\boldsymbol{x}}_{2}-{\boldsymbol{K}}_{p}{\boldsymbol{e}}+\Delta \delta \\ \qquad\quad +{\boldsymbol{u}}_{1}+\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\varepsilon {\boldsymbol{e}}+{\boldsymbol{K}}_{1}\varepsilon \dot{{\boldsymbol{e}}} \end{array} \end{cases} \end{equation}

where $\varepsilon \gt 0$ .

Based on HJI theorem, Eq. (21) is written in the form of Eq. (22):

(22) \begin{align} \begin{cases} {\boldsymbol{x}}={\boldsymbol{f}}({\boldsymbol{x}})+{\boldsymbol{g}}({\boldsymbol{x}})\boldsymbol{\lambda }\\ {\boldsymbol{z}}={\boldsymbol{l}}({\boldsymbol{x}}) \end{cases} \end{align}

where ${{\boldsymbol{f}}}\!\left({{\boldsymbol{x}}}\right)=\left[\begin{array}{c} {{\boldsymbol{x}}}_{2}-\varepsilon {{\boldsymbol{x}}}_{1}\\[5pt] {{{\boldsymbol{K}}}_{1}}^{-1}({-}({{\boldsymbol{K}}}_{2}+{{\boldsymbol{K}}}_{d}){{\boldsymbol{x}}}_{2}-{{\boldsymbol{K}}}_{p}{{\boldsymbol{e}}}+{\boldsymbol{W}_{f}}^{\ast }\boldsymbol{h}({{\boldsymbol{x}}})\\[5pt] +\left.\boldsymbol{u}_{1}+\left({{\boldsymbol{K}}}_{2}+{{\boldsymbol{K}}}_{d}\right)\varepsilon {{\boldsymbol{e}}}+{{\boldsymbol{K}}}_{1}\varepsilon \dot{{{\boldsymbol{e}}}}\right) \end{array}\right], {{\boldsymbol{g}}}\!\left({{\boldsymbol{x}}}\right)=\left[\begin{array}{l} 0\\[5pt] {{{\boldsymbol{K}}}_{1}}^{-1} \end{array}\right], \boldsymbol{\lambda }=\xi _{f}$ .

Define the evaluation index ${\boldsymbol{z}}$ : since $\boldsymbol{\lambda }=\xi _{f}$ , that is, the approximation error is considered as external interference $\xi _{f}$ , and the evaluation index is defined as ${\boldsymbol{z}}={\boldsymbol{x}}_{2}=\dot{{\boldsymbol{e}}}+\varepsilon {\boldsymbol{e}}$ , then L ₂ gain is $J=\underset{\left|\left|\xi _{f}\right|\right|\neq 0}{\sup }\frac{\left|\left|{\boldsymbol{z}}\right|\right|_{2}}{\left|\left|\xi _{f}\right|\right|_{2}}$ . HJI theorem is for a positive number $\eta$ , if there is a positive definite and differentiable function $V(x)\geqslant 0$ and $\dot{V}\!\left(x\right)\leqslant \dfrac{1}{2}\left\{\eta ^{2}\left|\left|\boldsymbol{\lambda }\right|\right|^{2}-\left|\left|{\boldsymbol{z}}\right|\right|^{2}\right\}\left(\forall \boldsymbol{\lambda }\right)$ , then $J\leqslant \eta$ .

Based on Eq. (21), the following adaptive control law is designed

(23) \begin{equation} \dot{\hat{{\boldsymbol{W}}}}_{f}=\alpha {\boldsymbol{x}}_{2}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)^{\mathrm{T}} \end{equation}

where $\alpha \gt 0$ .

The feedback control law is designed as follows:

(24) \begin{equation} {\boldsymbol{u}}_{1}=-\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\varepsilon {\boldsymbol{e}}-{\boldsymbol{K}}_{1}\varepsilon \dot{{\boldsymbol{e}}}-\frac{1}{2}{\boldsymbol{x}}_{2}+{\boldsymbol{K}}_{p}{\boldsymbol{e}}-\frac{1}{2\eta ^{2}}{\boldsymbol{x}}_{2}-\frac{1}{2}\dot{{\boldsymbol{K}}}_{1}{\boldsymbol{x}}_{2}-\hat{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right) \end{equation}

where $\hat{{\boldsymbol{W}}}_{f}$ is the estimated weight of the neural network. Then the closed-loop system, that is, Eq. (19) satisfies $J\leq \eta$ .

5.2. Stability analysis

In order to verify the stability of the closed-loop system, based on refs. [Reference Liu25–Reference Lou, Lin, Quan, Wei and Di30], the following Lyapunov function is constructed

(25) \begin{equation} V=\frac{1}{2}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{K}}_{1}{\boldsymbol{x}}_{2}+\frac{1}{2\alpha }\left|\left|{{\boldsymbol{W}}}_{f}\right|\right|^{2} \end{equation}

where $\tilde{{\boldsymbol{W}}}_{f}=\hat{{\boldsymbol{W}}}_{f}-{{\boldsymbol{W}}_{f}}^{\mathrm{*}}, \alpha \gt 0$ . Since ${\boldsymbol{K}}_{1}$ is a positive definite matrix, $V\gt 0$ .

Let $\left|\left|{\boldsymbol{P}}\right|\right|^{2}=\mathrm{tr}({\boldsymbol{P}}{\boldsymbol{P}}^{\mathrm{T}})=\mathrm{tr}({\boldsymbol{P}}^{\mathrm{T}}{\boldsymbol{P}})$ , where $\mathrm{tr}({\cdot})$ is the matrix trace, then $\left|\left|{{\boldsymbol{W}}}_{f}\right|\right|^{2}=\mathrm{tr}\!\left({{{\boldsymbol{W}}}_{f}}^{\mathrm{T}}{{\boldsymbol{W}}}_{f}\right)$ .

By deriving $V$ , the following results can be obtained:

(26) \begin{align} \dot{V} & = {{\boldsymbol{x}}_{2}} {}^{\mathrm{T}}{\boldsymbol{K}}_{1}\dot{{\boldsymbol{x}}}_{2}+\frac{1}{2}{{\boldsymbol{x}}_{2}} {}^{\mathrm{T}}\dot{{\boldsymbol{K}}}_{1}{\boldsymbol{x}}_{2}+\frac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}} {}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right)\nonumber\\ &={{\boldsymbol{x}}_{2}} {}^{\mathrm{T}}\!\left({-}\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\!{\boldsymbol{x}}_{2}-{\boldsymbol{K}}_{p}{\boldsymbol{e}}+\Delta \delta +{\boldsymbol{u}}_{1}+\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\varepsilon {\boldsymbol{e}}+{\boldsymbol{K}}_{1}\varepsilon \dot{{\boldsymbol{e}}}\right)+\frac{1}{2}{{\boldsymbol{x}}_{2}} {}^{\mathrm{T}}\dot{{\boldsymbol{K}}}_{1}{\boldsymbol{x}}_{2}+\frac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}} {}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right)\nonumber\\ &={{\boldsymbol{x}}_{2}}^{\mathrm{T}}\!\left({-}\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\!{\boldsymbol{x}}_{2}-{\boldsymbol{K}}_{p}{\boldsymbol{e}}+{{\boldsymbol{W}}_{f}}^{\ast }{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)+\xi _{f}-\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\varepsilon {\boldsymbol{e}}-{\boldsymbol{K}}_{1}\varepsilon \dot{{\boldsymbol{e}}}-\frac{1}{2\eta ^{2}}{\boldsymbol{x}}_{2}\right.\nonumber\\ &\quad -\left.\hat{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)+\left({\boldsymbol{K}}_{2}+{\boldsymbol{K}}_{d}\right)\varepsilon {\boldsymbol{e}}+{\boldsymbol{K}}_{1}\varepsilon \dot{{\boldsymbol{e}}}-\frac{1}{2}{\boldsymbol{x}}_{2}+{\boldsymbol{K}}_{p}{\boldsymbol{e}}\right)+\frac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}} {}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right) \end{align}

By simplifying Eq. (26), it can be concluded:

(27) \begin{equation} \dot{V}=-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\!\left({\boldsymbol{K}}_{d}+{\boldsymbol{K}}_{2}\right){\boldsymbol{x}}_{2}-\frac{1}{2\eta ^{2}}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2}-\frac{1}{2}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2}-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)+{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\xi _{f}+\frac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}} {}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right) \end{equation}

Definition: $L=\dot{V}-\frac{1}{2}\eta ^{2}\left|\left|\xi _{f}\right|\right|^{2}+\frac{1}{2}\left|\left|{\boldsymbol{z}}\right|\right|^{2}$ , then,

(28) \begin{align} L= & -{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{K}}_{d}{\boldsymbol{x}}_{2}-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{K}}_{2}{\boldsymbol{x}}_{2}-\frac{1}{2\eta ^{2}}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2}-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)+\frac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right)+{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\xi _{f}-\frac{1}{2}\eta ^{2}\left|\left|\xi _{f}\right|\right|^{2}\nonumber \\& +\frac{1}{2}\left|\left|{\boldsymbol{z}}\right|\right|^{2}-\frac{1}{2}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2} \end{align}

In Eq. (28),

(29) \begin{equation} -{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{K}}_{d}{\boldsymbol{x}}_{2}-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{K}}_{2}{\boldsymbol{x}}_{2}\leqslant 0 \end{equation}

(30) \begin{equation} -\frac{1}{2\eta ^{2}}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2}+{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\xi _{f}-\frac{1}{2}\eta ^{2}\left|\left|\xi _{f}\right|\right|^{2}=-\frac{1}{2}\left|\left|\frac{1}{\eta }{{\boldsymbol{x}}_{2}}^{\mathrm{T}}-\eta \xi _{f}\right|\right|^{2}\leqslant 0 \end{equation}

(31) \begin{equation} \frac{1}{2}\left|\left|{\boldsymbol{z}}\right|\right|^{2}-\frac{1}{2}{{\boldsymbol{x}}_{2}}^{\mathrm{T}}{\boldsymbol{x}}_{2}=0 \end{equation}

(32) \begin{equation} \begin{array}{l} \dfrac{1}{\alpha }\mathrm{tr}\!\left({\dot{\tilde{{\boldsymbol{W}}}}_{f}} {}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}\right)-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)=\dfrac{1}{\alpha }\mathrm{tr}\!\left({\tilde{{\boldsymbol{W}}}_{f}} {}^{\mathrm{T}}\alpha {\boldsymbol{x}}_{2}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)^{\mathrm{T}}\right)-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)\\ =\mathrm{tr}\!\left({\tilde{{\boldsymbol{W}}}_{f}}^{\mathrm{T}}{\boldsymbol{x}}_{2}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)^{\mathrm{T}}\right)-{{\boldsymbol{x}}_{2}}^{\mathrm{T}}\tilde{{\boldsymbol{W}}}_{f}{\boldsymbol{h}}\!\left({\boldsymbol{x}}\right)=0\\ \end{array} \end{equation}

Based on the above conditions, $L\leqslant 0$ can be obtained. According to the definition of $L, \dot{V}-\frac{1}{2}\eta ^{2}\left|\left|\xi _{f}\right|\right|^{2}+\frac{1}{2}\left|\left|{\boldsymbol{z}}\right|\right|^{2}\leqslant 0$ can be obtained, that is, $\dot{V}\leqslant 0$ .

Therefore, based on the HJI theorem, $J\leqslant \eta$ can be obtained, so that $\left|\left|{\boldsymbol{z}}\right|\right|$ satisfies the performance index, ${\boldsymbol{e}}$ and $\dot{{\boldsymbol{e}}}$ satisfy the convergence requirements.

6. Simulation analysis and results

6.1. Motion trajectory planning

When WWRTPR carries out lumbar rehabilitation training for the patients, it makes the patient’s waist perform rotational rehabilitation training around the three coordinate axis (X axis, Y axis, Z axis) in the effective range according to the different degree of lumbar disease, that is, back and forth swing, left and right swing and rotation. Based on the principle of motion equivalence, the human lumbar joint is simplified as a ball hinge, and the human body is simplified as two elliptical columns connected by the ball hinge, as shown in Fig. 6.

Table I. Coordinates of P _i point and B _i point.

Figure 6. Simplified schematic diagram of WWRTPR.

In rehabilitation training, the planning of lumbar rehabilitation motion trajectory plays an important role. Unreasonable motion trajectory will cause theoretical conflict between the length and force of the wire, which will cause secondary injury to the patients. Therefore, the rehabilitation motion trajectory is reasonably planned, and the rehabilitation motion trajectory of the wire-driven lumbar rehabilitation training robot is defined as follows:

(33) \begin{align} \begin{cases} \theta =15^{^{\circ}}\sin \left(\frac{\Pi }{10}t+\Pi \right)\\[5pt] \phi =15^{^{\circ}}\sin \left(\frac{\Pi }{10}t+\Pi \right)\\[5pt] \psi =7^{^{\circ}}\sin \left(\frac{\Pi }{10}t+\Pi \right) \end{cases} \end{align}

where $\theta, \phi, \psi$ are the rotation angles around the y-axis, x-axis and z-axis, respectively.

6.2. Simulation analysis

In order to verify the correctness and effectiveness of the designed control strategy, the force/position hybrid control strategy-based RBF neural network is simulated and analyzed for the four-wire three-degree-of-freedom wire-driven parallel robot for rehabilitation training. The simulation analysis in this section is based on the 2.5:1 scale model, and the simulation is carried out by MATLAB software. The coordinates of connection point P _i of the end effector (patient) and pulley hinge point B _i are shown in Table I, and the coordinates of P _i point and B _i point are basically symmetrically distributed.

The composite motion simulation results of lumbar rehabilitation training are shown in Figs. 7–13.

Figure 7. Angle trajectory curve of waist compound movement.

Figures 7 and 8 are, respectively, the angle trajectory curve and angle error curve of lumbar compound motion. From Fig. 7, it can be concluded that the curve of the angle trajectory of the lumbar compound movement is smooth and continuous, the angle trajectory is 0 in the median, and it reaches the expected limit rehabilitation value in the anterior and posterior position. From Fig. 8, it can be concluded that the angle error is within ± 0.13°, the error is small and changes within a certain range, meeting the conditions of waist compound motion.

Figure 8. Angle error curve of waist compound movement.

Figures 9 and 10 are, respectively, the expected motor angular speed curve and the motor speed error curve. It can be concluded from Fig. 9 that the expected angular velocity curve of the motor is smooth and continuous, and has certain regularity. From Fig. 10, it can be concluded that the motor speed error is within ± 0.15 rad/s, and the error changes within a certain range, which meets the basic conditions of waist compound motion simulation.

Figure 9. Motor desired angular speed of waist compound movement.

Figure 10. Error variation of motor speed.

Figures 11, 12 and 13 are, respectively, the wire length variation curve, wire tension variation curve and control input curve. It can be concluded from Fig. 11 that the wire length variation curve is smooth and continuous, which is similar to the trigonometric function variation, and the wire length variation is within a certain range and has certain regularity. It can be seen from Fig. 12 that the wire tension variation is also within a certain range, and the variation is gentle, which can reduce the risk of wire interruption during rehabilitation training. It can be seen from Fig. 13 that the control input varies within the rated torque range of the motor (1.27 N. M), which conforms to the physical meaning of the design.

Figure 11. Variation curve of wire length.

Figure 12. Variation curve of wire tension.

Figure 13. Control input curve.

In conclusion, the simulation results verify the correctness, feasibility and effectiveness of the designed control strategy, and the simulation results show that WWRTPR can achieve the rehabilitation training of patients’ waist, and has a certain effect on patients’ waist rehabilitation.