3.1. Equivalent multi-body mechanics model of the human spine

Due to the force characteristics of the human lumbar region and the pathogenesis of lumbar diseases, many researchers have conducted extensive studies on an abstract model capable of delineating the force parameters of the human spine (mechanical modeling). This approach is preferred due to safety concerns and the complexity associated with human experiments. Roberts and Chen developed a finite element model for the human spine from T1 to S1 [Reference Roberts and Chen27]. They used Euler-Bernoulli beam elements to depict the static elastic deformation problem of the spine under force. However, the model does not account for the influence of soft tissues present in the spinal region. Shirazi et al. considered the vertebral body a rigid structure and connected the cones to a beam unit to simulate intervertebral discs [Reference Shirazi-Adl and Pamianpour28]. They analyzed the load-bearing capacity of the lumbar spine under various motion postures by applying a bending moment to L1 to mimic the bending force. The analysis took into account the soft tissues of the spine, but did not consider the intricate muscular forces of the lumbar and dorsal regions. Matteo et al. developed a detailed finite element model of the musculoskeletal spine [Reference Panico, Bassani, Villa and Galbusera29]. They applied simplified loading conditions to evaluate bending deformation and determine the bending angle and stress distribution of spine segments under specific loads. However, complex spine geometry and muscle models present challenges for dynamic analysis. At the initiation of the stance phase, activate the stepper motor to drive the internal ratchet clockwise. With the ongoing rotation, the linear spring progressively extends, initiating the energy accumulation.

In summary, most of the current literature concentrates on finite element-based static force analysis of the spine [Reference Roberts and Chen27–Reference Panico, Bassani, Villa and Galbusera29]. However, there are few studies related to the dynamic characteristics of the spine and its human-exoskeleton wearable parallel system with exoskeleton. Hence, this study utilized a multi-body system dynamics approach [Reference Shabana30] to simplify the human lumbar region, taking into account spinal and back muscle tension. The lumbar region is represented as multiple sections of Euler–Bernoulli beams [Reference Shirazi-Adl and Pamianpour28, Reference Panico, Bassani, Villa and Galbusera29] with varying stiffness (refer to Table II), as shown in Fig. 2.

Figure 2.

Study sketch of the equivalent multi-body system.

3.2. Analysis of system kinematics

Due to the broad applicability of kinematics and dynamics modeling techniques and space constraints, this study employs Cartesian coordinates to construct kinematics and dynamics models of the system. It specifically examines the motion of the human-exoskeleton parallel waist exoskeleton in the sagittal plane of the human body. Fig. 2(a) shows a diagram illustrating the mechanism of a human–exoskeleton parallel waist exoskeleton in the sagittal plane. Here, $O^{0}X^{0}Y^{0}$ represents the global coordinate system of the system, and $O^{i}X^{i}Y^{i}$ represents the body coordinate system of component $i(=1,2,\cdots ,8)$ . Component 8 represents the component of the human-exoskeleton interaction, while component 7 simulates the rigidity of the human chest. Mobile vice constraints are placed between components 7 and 8 to simulate displacements at the point where human-exoskeleton interaction forces occur during human movement. Most studies on human motion using finite element analysis have typically employed simplified loading conditions, which are easier to implement in the human body than complex muscle forces. In this research, we specifically apply the simplified loading conditions of the pure moment type [Reference Panico, Bassani, Villa and Galbusera29] M to model the motion of the lumbar spine (refer to Fig. 2), with $F_{n}$ denoting the load supported by the human body. The design of flexible structures and actuators that mimic human muscle functions allows joint assistance without rigid constraints. Based on this foundation, a flexible body finite element model of the lumbar spine, identified as component 6 in Fig. 2(a), is put forth. Flexible actuation offers advantages by circumventing the constraints imposed by rigid mechanisms on joints, decreasing system weight, and enhancing the exoskeleton’s wearability. The structure of Component 5 resembles a flexible beam [Reference Shirazi-Adl and Pamianpour28], which is equivalent to a human waist, as depicted in Fig. 2(b). To depict the bending motion of the human lumbar region, a model is proposed in which five discrete flexible components are linked to a rigid component 5. The rigid motion of component $5_{j}\,(j=1,2,\cdots ,5)$ represents the movement of human vertebrae L1-L5, while the flexible deformation of component $6_{k}\,(k=1,2,\cdots ,5)$ describes the flexibility of the intervertebral disc. The figure shows $O_{j}^{5}X_{j}^{5}Y_{j}^{5}$ as the body coordinate system of the rigid equivalent vertebrae $j(\!=1,2,\cdots ,5)$ and $O_{k}^{6}X_{k}^{6}Y_{k}^{6}$ as the body coordinate system of the flexible equivalent intervertebral disc $k(\!=1,2,\cdots ,5)$ . The flexible body coordinates are defined using the floating coordinate method [Reference Shabana30]. The position vector $\mathbf{r}^{i}$ at any point on component $i$ can be described as:

(1) \begin{equation} \mathbf{r}^{i}=\mathbf{R}^{i}+\mathbf{A}^{i}\mathbf{S}^{ij}\left(\overline{\mathbf{q}}_{o}^{ij}+\overline{\mathbf{q}}_{f}^{ij}\right) \end{equation}

The equation uses an upper horizontal line to denote measurements in the body coordinate system and the absence of the line to denote measurements in the global coordinate system. $\mathbf{R}^{\mathbf{i}}$ represents the position vector of the origin of component $i$ ’s body coordinate system in the global coordinate system. The coordinates of the component unit nodes before and after deformation are represented by $\overline{\mathbf{q}}_{o}^{ij}$ and $\overline{\mathbf{q}}_{f}^{ij}$ , respectively. The variable j indicates the serial number of the discrete unit. For rigid components, $\overline{\mathbf{q}}_{f}^{ij}=\mathbf{0}_{3\times 3}$ . $\mathbf{A}^{i}$ is the transformation matrix from the body coordinate system to the global coordinate system and can be defined as:

(2) \begin{equation} \mathbf{A}^{i}=\left[\begin{array}{l@{\quad}l} \cos \theta ^{i} & -\!\sin \theta ^{i}\\ \sin \theta ^{i} & \cos \theta ^{i} \end{array}\right]_{2\times 2} \end{equation}

The rotation angle of component $i$ ’s body coordinate system relative to the global coordinate system is denoted as $\theta ^{i}$ . The flexible body shape function, denoted by $\mathbf{S}^{ij}$ , is defined as follows:

(3) \begin{equation} \mathbf{S}^{ij}=\left[\begin{array}{c@{\quad}c@{\quad}c} \xi ^{ij} & 0 & 0\\ 0 & 3\xi ^{{ij^{2}}}-2\xi ^{{ij^{3}}} & l^{ij}\!\left(\xi ^{{ij^{3}}}-\xi ^{{ij^{2}}}\right) \end{array}\right]_{2\times 3} \end{equation}

Within equation $\xi ^{ij}=\,\overline{x}^{ij}/l^{ij} , \overline{x}^{ij}$ denotes the x-coordinate of a point on discrete element j of flexible component $i$ , measured in its body coordinate system before deformation, while $l^{ij}$ represents the length of the component.

The elastic coordinates for each unit of flexible component 6 are denoted as $\overline{\mathbf{q}}_{f}^{\mathrm{6}j}=[\begin{array}{llll} q_{f\mathrm{1}}^{6j} & q_{f\mathrm{2}}^{6j} & \cdots & q_{f\mathrm{6}}^{6j} \end{array}]_{1\times 6}^{\mathrm{T}}$ following the Euler–Bernoulli beam theory. According to the cantilever beam boundary condition, one side of the node has zero degrees of freedom, whereas the elastic deformation of flexible component 6 is described using only three elastic coordinates on the other side: $\overline{\mathbf{q}}_{f}^{\mathrm{6}j}=[\begin{array}{lll} q_{f\mathrm{4}}^{6j} & q_{f\mathrm{5}}^{6j} & q_{f\mathrm{6}}^{6j} \end{array}]_{1\times 3}^{\mathrm{T}}$ . Therefore, the independent elastic coordinate vector for representing the equivalent flexible beam at the waist of the human body is:

(4) \begin{equation} \overline{\mathbf{q}}_{f}^{\mathrm{6}}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \overline{q}_{f}^{61} & \overline{q}_{f}^{62} & \cdots & \overline{q}_{f}^{65} \end{array}\right]^{\mathrm{T}}=\left[\begin{array}{lllllll} q_{f\mathrm{4}}^{61} & q_{f\mathrm{5}}^{61} & q_{f\mathrm{6}}^{61} & \cdots & q_{f\mathrm{4}}^{65} & q_{f\mathrm{5}}^{65} & q_{f\mathrm{6}}^{65} \end{array}\right]_{1\times 15}^{\mathrm{T}} \end{equation}

The configuration vector for component $i$ can be characterized as:

(5) \begin{equation} \left\{\begin{array}{l@{\quad}l} \,\,\mathbf{q}^{i} =\left[\begin{array}{l@{\quad}l@{\quad}l} r_{x}^{i} & r_{y}^{i} & \theta ^{i} \end{array}\right]_{1\times 3}^{\mathrm{T}} &i\neq 6\\ \mathbf{q}^{\mathrm{6}j} =\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} r_{x}^{\mathrm{6}j} & r_{y}^{\mathrm{6}j} & \theta ^{\mathrm{6}j} & \left(\overline{\mathbf{q}}_{f}^{\mathrm{6}j}\right)^{\mathrm{T}} \end{array}\right]_{1\times 6}^{\mathrm{T}} &i=6 \end{array}\right. \end{equation}

where $r_{x}^{i}$ and $r_{y}^{i}$ represent the $x$ and $y$ coordinates of component $i$ ’s body coordinate system origin with respect to the global coordinate system, respectively.

The kinematics constraints of the human-exoskeleton parallel waist exoskeleton consist of both kinematics constraints and rigid coupling constraints between the discrete units of building blocks 5 and 6. The kinematics constraint equations ( $\mathbf{C}_{1}=\mathbf{0}_{10\times 1}$ ) from the kinematics vice are included:

(6) \begin{equation} \left[\begin{array}{l} \mathbf{R}^{1}+\mathbf{A}^{1}\overline{\mathbf{u}}_{\mathbf{O}^{1}}^{1}-\mathbf{r}_{{\mathbf{O}^{1}}}\\[3pt] \mathbf{R}^{2}+\mathbf{A}^{2}\overline{\mathbf{u}}_{\mathbf{O}^{2}}^{2}-\mathbf{R}^{1}-\mathbf{A}^{1}\overline{\mathbf{u}}_{\mathbf{O}^{2}}^{1}\\[3pt] \mathbf{R}^{3}+\mathbf{A}^{3}\overline{\mathbf{u}}_{\mathbf{O}^{3}}^{3}-\mathbf{R}^{2}-\mathbf{A}^{2}\overline{\mathbf{u}}_{\mathbf{O}^{3}}^{2}\\[3pt] \mathbf{R}^{4}+\mathbf{A}^{4}\overline{\mathbf{u}}_{\mathbf{O}^{4}}^{4}-\mathbf{R}^{3}-\mathbf{A}^{3}\overline{\mathbf{u}}_{\mathbf{O}^{4}}^{3}\\[3pt] r_{x}^{7}-\left(r_{y}^{8}-r_{y}^{7}\right)\times \sin \theta ^{7}-r_{x}^{8}\\[3pt] \theta ^{7}-\theta ^{8} \end{array}\right]=\mathbf{0}_{10\times 1} \end{equation}

The vector $\overline{\mathbf{u}}_{O^{j}}^{i}$ describes the position of point $O^{j}$ relative to component $i$ in the body coordinate system, whereas $\mathbf{r}_{{O^{\mathrm{1}}}}$ represents the position vector of point $O^{1}$ in the global coordinate system. The rotational sub-hinge points between the components of the exoskeleton are denoted as $O^{1} , O^{2} , O^{3}$ , and $O^{4}$ . The x and y coordinates of components 7 and 8 in the global coordinate system are denoted as $r_{x}^{7} , r_{y}^{7} , r_{x}^{8}$ , and $r_{y}^{8}$ . Although the discrete units of components 5 and 6 are rigidly coupled, the flexible deformation of component 6 enables rigid body motion between them. Based on the provisions of the unit coordinate system in Fig. 2, the rigid coupling constraint equations ( $\mathbf{C}_{2}=\mathbf{0}_{33\times 1}$ ) for each discrete unit and between components 5 and 7 should be added as follows:

(7) \begin{equation} \left[\begin{array}{l} \mathbf{R}^{65}+\mathbf{A}^{65}\overline{\mathbf{u}}_{O_{5}^{6}}^{65}-\mathbf{r}_{O_{5}^{6}}\\[3pt] \mathbf{R}^{\mathrm{5}j}+\mathbf{A}^{5j}\overline{\mathbf{u}}_{O_{j}^{5}}^{5j}-\mathbf{R}^{6j}-\mathbf{A}^{6j}\overline{\mathbf{u}}_{O_{j}^{5}}^{6j}\\[3pt] \mathbf{R}^{6k}+\mathbf{A}^{6k}\overline{\mathbf{u}}_{O_{k}^{6}}^{6k}-\mathbf{R}^{5g}-\mathbf{A}^{5g}\overline{\mathbf{u}}_{O_{k}^{6}}^{5g}\\[3pt] \mathbf{R}^{7}+\mathbf{A}^{7}\overline{\mathbf{u}}_{O^{7}}^{7}-\mathbf{R}^{51}-\mathbf{A}^{51}\overline{\mathbf{u}}_{O^{7}}^{51}\\[3pt] \theta ^{65}-\frac{\pi }{2}\\[3pt] \theta ^{55}-q_{f\mathrm{6}}^{65}\\[3pt] \theta ^{6k}-\theta ^{5g}-\frac{\pi }{2}\\[3pt] \theta ^{5k}-\left(\theta ^{6k}+q_{f\mathrm{6}}^{6k}-\frac{\pi }{2}\right)\\[3pt] \theta ^{7}-\theta ^{51} \end{array}\right]=\mathbf{0}_{33\times 1} \end{equation}

In the given equation, $\mathbf{r}_{O_{5}^{6}}$ denotes the position vector of point $O_{5}^{6}$ within the global coordinate system, with $j(\!=1,2,\cdots ,5)$ and $k(\!=1,2,3,4)$ representing the discrete element numbers, and $g=k+1$ . Thus, the following system constraint equations can be obtained:

(8) \begin{equation} \mathbf{C}=\left[\mathbf{C}_{1}^{\mathrm{T}}\quad \mathbf{C}_{2}^{\mathrm{T}}\right]_{1\times 43}^{\mathrm{T}}=\mathbf{0}_{43\times 1} \end{equation}

Given the constraints imposed by system motion and rigid coupling, the independent coordinate method [Reference Korayem, Rahimi and Nikoobin31] offers a viable approach for isolating the independent variables of the system:

(9) \begin{equation} \mathbf{q}_{i}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l@{\quad}l} \theta ^{1} & \theta ^{2} & \theta ^{3} & \theta ^{4} & \left(\overline{\mathbf{q}}_{f}^{\mathrm{6}}\right)^{\mathrm{T}} & r_{y}^{8} \end{array}\right]_{1\times 20}^{\mathrm{T}} \end{equation}

Utilizing system-independent variables to represent system configuration variables can result in:

(10) \begin{equation} \mathbf{q}=f\!\left(\mathbf{q}_{i},t\right) \end{equation}

The expression for the system configuration variable, denoted as $\boldsymbol{q}$ , is as follows:

(11) \begin{equation} \mathbf{q}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \mathbf{q}^{{1^{{^{\mathrm{T}}}}}} & \mathbf{q}^{{2^{{^{\mathrm{T}}}}}} & \cdots & \mathbf{q}^{{8^{{^{\mathrm{T}}}}}} \end{array}\right]_{1\times 63}^{\mathrm{T}} \end{equation}

The variables associated with the system include the following:

(12) \begin{equation} \mathbf{q}_{d}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \mathbf{R}^{{1^{\mathrm{T}}}} & \cdots & \mathbf{R}^{{4^{\mathrm{T}}}} & \mathbf{q}^{{5^{\mathrm{T}}}} & \mathbf{q}_{c}^{6^{\mathrm{T}}} & \mathbf{q}^{{7^{\mathrm{T}}}} & r_{x}^{8} & \theta ^{8} \end{array}\right]_{1\times 43}^{\mathrm{T}} \end{equation}

Where

(13) \begin{align} \mathbf{q}^{5}&=\left[\begin{array}{l@{\quad}l@{\quad}l} \mathbf{q}^{{51^{\mathrm{T}}}} & \cdots & \mathbf{q}^{{55^{\mathrm{T}}}} \end{array}\right]_{1\times 15}^{\mathrm{T}} \end{align}

(14) \begin{align} \mathbf{q}_{c}^{6}&=\left[\begin{array}{l@{\quad}l@{\quad}l} \mathbf{q}_{c}^{61^{\mathrm{T}}} & \cdots & \mathbf{q}_{c}^{65^{\mathrm{T}}} \end{array}\right]_{1\times 15}^{\mathrm{T}} \end{align}

(15) \begin{align} \mathbf{q}_{c}^{6j}&=\left[\begin{array}{l@{\quad}l@{\quad}l} r_{x}^{6j} & r_{y}^{6j} & \theta ^{6j} \end{array}\right]_{1\times 3}^{\mathrm{T}}\left(j=1,\cdots ,5\right)\\[8pt]\nonumber \end{align}

3.3. Analysis of system dynamics

Derivation of both sides of Eq. (1) with respect to time yields the velocity vector $\dot{\mathbf{r}}^{i}$ at any point on component $i$ :

(16) \begin{equation} \dot{\mathbf{r}}^{i}=\mathbf{B}^{i}\dot{\mathbf{q}}^{i} \end{equation}

In Eq. (16), when $i=6 , \mathbf{B}^{i}=[\begin{array}{lll} \mathbf{I} & \mathbf{A}_{\theta }^{i}\overline{\mathbf{u}}^{i} & \mathbf{A}^{i}\mathbf{S}^{ij} \end{array}]_{2\times 3}$ , where $\mathbf{I}$ represents the identity matrix. For $i\neq 6 , \mathbf{B}^{i}=[\begin{array}{ll} \mathbf{I} & \mathbf{A}_{\theta }^{i}\overline{\mathbf{u}}^{i} \end{array}]_{2\times 3}$ . The derivative of the configuration vector of component $i$ with respect to time is indicated by the variable $\dot{\mathbf{q}}^{i}$ . The derivative of the attitude transformation matrix $\mathbf{A}^{i}$ with respect to the rotation angle $\theta ^{i}$ in the body coordinate system is denoted as $\mathbf{A}_{\theta }^{i}$ . The expression for the kinetic energy of component $i$ is as follows:

(17) \begin{equation} \mathbf{T}^{i}=\frac{1}{2}\int _{{V^{i}}}\dot{\mathbf{r}}^{{i^{\mathrm{T}}}}\dot{\mathbf{r}}^{i}\rho ^{i}dV^{i}=\frac{1}{2}\dot{\mathbf{q}}^{{i^{{^{\mathrm{T}}}}}}\mathbf{M}^{i}\dot{\mathbf{q}}^{i} \end{equation}

where $m^{i}$ represents the mass of component $i , V^{i}$ denotes the volume of component $i , \rho ^{i}$ signifies the density of component $i$ , and $\mathbf{M}^{i}$ stands for the inertia tensor of component $i$ , the expression for $\mathbf{M}^{i}$ is provided as follows:

(18) \begin{equation} \mathbf{M}^{i}=\int V^{i}\rho ^{i}\mathbf{B}^{i\,\mathrm{T}}\mathbf{B}^{i}dV^{i} \end{equation}

Then the inertia tensor of this human-exoskeleton wearable parallel exoskeleton system is:

(19) \begin{equation} \mathbf{M}_{t}=\left[\begin{array}{c@{\quad}c@{\quad}c} \mathbf{M}^{1} & \cdots & 0\\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathbf{M}^{8} \end{array}\right]_{63\times 63} \end{equation}

$\mathbf{M}^{i}$ can be found in Eq. (18) among the aforementioned elements.

In literature [Reference Lu, Wang and Feng32], a human-exoskeleton interaction force model was introduced. This model characterizes the human-exoskeleton interaction force between the exoskeleton and the human body by incorporating a spring-damper unit at the contact interface of the exoskeleton and the human body. In this study, we adopt a consistent methodology to develop the human-exoskeleton interaction force (see Fig. 2). Spring-damper units $K^{5}$ are introduced between the exoskeleton’s chest straps and the human chest equivalent beam to model the supplementary force applied. Additionally, spring-damper units $K^{6}-K^{9}$ are incorporated between the human waist equivalent beam and the exoskeleton to characterize the compressive force that the exoskeleton may apply to the human body during motion, denoted as the redundant human-exoskeleton interaction forces. For the exoskeleton to provide a booster effect to the human body during bending, spring damper units $K^{1}-K^{4}$ are added to the exoskeleton structure. The spring damper force vector $\mathbf{F}_{k}^{i}$ in the system is:

(20) \begin{align} \mathbf{F}_{k}^{i}&=\left(K^{i}\times \!\left(l^{{i^{o}}}-l^{i}\right)-C^{i}\times \dot{l}^{i}\right)\times \frac{\left({\unicode[Arial]{x0394}} \mathbf{r}\,^{abi}\right)}{\left| {\unicode[Arial]{x0394}} \mathbf{r}\,^{abi}\right| }\nonumber\\ {\unicode[Arial]{x0394}} \mathbf{r}^{abi}&=\mathbf{r}^{ai}-\mathbf{r}^{bi}\end{align}

In Eq. (20), $K^{i}$ and $C^{i}$ represent the stiffness and damping coefficients of the i-th spring damper unit, respectively. In addition, $C^{i}=k_{c}\times K^{i}$ , where $k_{c}$ denotes the damping ratio. Denoted as $l^{io} , l^{i}$ , and $\dot{l}^{i}$ , the original length, current length, and rate of change in length pertain to the i-th spring. $\mathbf{r}^{ai}$ and $\mathbf{r}^{bi}$ are the position vectors of the i-th spring damper element at both ends within the global coordinate system. Upon equating the spring damper element to a generalized force, the system’s total force model is illustrated in Fig. 3(a). The derivation of the generalized force for component 1, along with its force diagram depicted in Fig. 3(b), is needed. The external forces acting on component 1 include gravity ( $G_{1}$ ) at the center of mass point( $P^{1}$ ), the spring force ( $\mathbf{F}_{k}^{6}$ ) due to human-exoskeleton interactions at point $O^{1}$ , and the spring forces at point $A_{1}$ from the exoskeleton structure ( $\mathbf{F}_{k}^{1}$ and $\mathbf{F}_{k}^{2}$ ). In accordance with the principle of variation, the system exhibits the following characteristics:

Figure 3.

System of force models.

(21) \begin{equation} \left\{\begin{array}{l@{\quad}l} \delta \mathbf{r}_{A_{1}}^{1} & =\mathbf{B}_{A_{1}}^{1}\delta \mathbf{q}^{1}=\left[\begin{array}{l@{\quad}l} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{A_{1}}^{1} \end{array}\right]\delta \mathbf{q}^{1}\\[3pt] \delta \mathbf{r}_{O^{1}}^{1} & =\mathbf{B}_{O^{1}}^{1}\delta \mathbf{q}^{1}=\left[\begin{array}{l@{\quad}l} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{O^{1}}^{1} \end{array}\right]\delta \mathbf{q}^{1}\\[3pt] \delta \mathbf{r}_{P^{1}}^{1} & =\mathbf{B}_{P^{1}}^{1}\delta \mathbf{q}^{1}=\left[\begin{array}{l@{\quad}l} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{P^{1}}^{1} \end{array}\right]\delta \mathbf{q}^{1} \end{array}\right. \end{equation}

In Eq. (21), $\overline{\mathbf{u}}_{A_{1}}^{1} , \overline{\mathbf{u}}_{O^{1}}^{1}$ and $\overline{\mathbf{u}}_{P^{1}}^{1}$ represent the position vectors of points $A_{1} , O^{1}$ and $P^{1}$ in the coordinate system of component 1. By adhering to the principle of virtual work, the generalized force can be calculated as:

(22) \begin{equation} \begin{array}{l@{\quad}l} \mathbf{Q}_{1}^{1^{{^{{^{\mathrm{T}}}}}}} & =\left(\left(\mathbf{F}_{k}^{1^{{^{{^{\mathrm{T}}}}}}}+\mathbf{F}_{k}^{2^{{^{\mathrm{T}}}}}\right)\left[\begin{array}{ll} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{O^{1}}^{1} \end{array}\right]\right)^{\mathrm{T}}\\[3pt] \mathbf{Q}_{2}^{1^{{^{{^{\mathrm{T}}}}}}} & =\left(\mathbf{F}_{k}^{6^{{^{{^{\mathrm{T}}}}}}}\left[\begin{array}{l@{\quad}l} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{A_{1}}^{1} \end{array}\right]\right)^{\mathrm{T}}\\ \mathbf{Q}_{g}^{1^{{^{\mathrm{T}}}}} & =\left(\left[\begin{array}{l@{\quad}l} 0 & G^{1} \end{array}\right]\left[\begin{array}{l@{\quad}l} \mathbf{I} & \mathbf{A}_{\theta }^{1}\overline{\mathbf{u}}_{P^{1}}^{1} \end{array}\right]\right)^{\mathrm{T}}\\[3pt] \mathbf{Q}^{{1^{{^{\mathrm{T}}}}}} & =\mathbf{Q}_{1}^{1^{{^{\mathrm{T}}}}}+\mathbf{Q}_{2}^{1^{{^{\mathrm{T}}}}}+\mathbf{Q}_{g}^{1^{{^{\mathrm{T}}}}} \end{array} \end{equation}

Correspondingly, the generalized forces $\mathbf{Q}^{{2^{{^{\mathrm{T}}}}}} , \mathbf{Q}^{{3^{{^{\mathrm{T}}}}}} , \mathbf{Q}^{{4^{{^{\mathrm{T}}}}}} , \mathbf{Q}^{{5^{{^{\mathrm{T}}}}}} , \mathbf{Q}^{{6^{{^{\mathrm{T}}}}}} , \mathbf{Q}^{{7^{{^{\mathrm{T}}}}}}$ and $\mathbf{Q}^{{8^{{^{\mathrm{T}}}}}}$ for the other elements can be calculated, but their derivations are omitted here due to space constraints. In conclusion, the generalized force of the human–exoskeleton parallel waist exoskeleton system [Reference Korayem, Rahimi and Nikoobin31] can be summarized as follows:

(23) \begin{equation} \mathbf{Q}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \mathbf{Q}^{{1^{{^{\mathrm{T}}}}}} & \mathbf{Q}^{{2^{{^{\mathrm{T}}}}}} & \cdots & \mathbf{Q}^{{8^{{^{\mathrm{T}}}}}} \end{array}\right]_{1\times 60}^{\mathrm{T}} \end{equation}

In the simplified multibody dynamics model of the human lumbar spine, component 5 (the lumbar vertebrae) is represented as a rigid beam, while component 6 (the intervertebral disc) is represented as a flexible beam. These two elements are rigidly coupled in the model. By utilizing strain energy and the Euler–Bernoulli principle (without accounting for shear deformation effects), the stiffness matrix of discrete flexible elements in component 6 can be defined. The deformation energy of a discrete element $j$ on flexible component 6 is expressed as follows:

(24) \begin{align}U^{6j} & =\frac{1}{2}\int _{0}^{l^{6j}}\left[\begin{array}{l@{\quad}l} {\overline{\mathbf{u}}}^{\prime}_{f6j} & {\overline{\mathbf{u}}}^{\prime\prime}_{f6j} \end{array}\right]\left[\begin{array}{l@{\quad}l} Ea & 0\\ 0 & EI \end{array}\right]\left[\begin{array}{l} {\overline{\mathbf{u}}}^{\prime}_{f6j}\\ {\overline{\mathbf{u}}}^{\prime\prime}_{f6j} \end{array}\right]dx\nonumber\\ & =\frac{1}{2}\overline{\mathbf{q}}_{f}^{6j^{{^{\mathrm{T}}}}}\mathbf{K}_{ff}^{6j}\overline{\mathbf{q}}_{f}^{6j} \end{align}

In Eq. (24), ${\overline{\mathbf{u}}}^{\prime}_{f\mathrm{6}j}=\partial \overline{\mathbf{u}}_{f\mathrm{6}j}/\partial (\overline{x}^{\mathrm{6}j}) , {\overline{\mathbf{u}}}^{\prime\prime}_{f\mathrm{6}j}=\partial ^{2}\overline{\mathbf{u}}_{f\mathrm{6}j}/\partial (\overline{x}^{\mathrm{6}j})^{2}$ , and $\overline{\mathbf{u}}_{f\mathrm{6}j}=\mathbf{S}^{\mathrm{6}j}\overline{\mathbf{q}}_{f}^{f\mathrm{6}j}$ . The elastic modulus, moment of inertia, and cross-sectional area of the flexible component are denoted by $E , I$ and $a$ , respectively. This paper employs a circular cross-section equivalent to the axial cross-section of the human vertebrae; thus, $a=\pi D^{2}/4$ and $I=\pi D^{4}/64$ . The parameter $D$ represents the circular cross-sectional diameter, while $\mathbf{K}_{ff}^{6j}$ denotes the stiffness matrix of a discrete element $j$ located on flexible component 6. It can be described as:

(25) \begin{equation} \mathbf{K}_{ff}^{6j}=\int _{0}^{l^{6j}}\left(\mathbf{S}_{x}^{6j}\right)^{\mathrm{T}}\left[\begin{array}{l@{\quad}l} Ea & 0\\ 0 & EI \end{array}\right]\mathbf{S}_{x}^{6j}dk \end{equation}

In Eq. (25), $\mathbf{S}_{x}^{\mathrm{6}j}=\partial \mathbf{S}^{\mathrm{6}j}/\partial \overline{x}^{\mathrm{6}j}$ . The stiffness matrix pertaining to component 6 is as follows:

(26) \begin{equation} \mathbf{M}_{t}=\left[\begin{array}{l@{\quad}l@{\quad}l} \mathbf{K}_{ffc}^{61} & \cdots & 0\\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathbf{K}_{ffc}^{65} \end{array}\right]_{30\times 30} \end{equation}

In Eq. (26), $\mathbf{K}_{ffc}^{\mathrm{6}j}=\left[\begin{array}{l@{\quad}l} \mathbf{0}_{3\times 3} & \mathbf{0}_{3\times 3}\\[3pt] \mathbf{0}_{3\times 3} & \mathbf{K}_{ff}^{\mathrm{6}j} \end{array}\right]_{6\times 6}$ . Using the bending deformation data from Matteo’s established musculoskeletal spine model [Reference Panico, Bassani, Villa and Galbusera29], one can calculate the equivalent stiffness of the flexible beam in the simplified lumbar spine model during flexion motion. Reference [Reference Bazrgari, Shirazi-Adl and Parnianpour33] conducted measurements on the dimensions of human spinal vertebrae and intervertebral discs, as presented in Table I. Reference [Reference Pearsall, Reid and Livingston34] conducted measurements on the masses of vertebrae and intervertebral discs. The average of the above data was calculated in a simplified lumbar spine model in this paper, which also integrates the mass of the human lumbar, head, and upper limbs into the spine model. The lumbar vertebrae measure 22 mm in length, while the intervertebral discs measure 10 mm. A circular cross-section with a 40 mm diameter is utilized to substitute for the axial cross-section of both the vertebrae and intervertebral discs. The equivalent beam parameters of the spine in the simplified human body model, as outlined in Table II, are presented herein.

Table I.

Range of vertebral bone and disc size parameters in the literature [Reference Bazrgari, Shirazi-Adl and Parnianpour33].

Table II.

Simplification of spine equivalent beam parameters in the model.

Based on the dynamics of multibody systems [Reference Shabana30, Reference Yuwei, Bin, Haibo, Xinhua, Yang, Liang and Lei35], one can derive the dynamic equations for each component of the human-exoskeleton parallel waist exoskeleton system as follows:

(27) \begin{equation} \mathbf{M}_{t}\overset{..}{\mathbf{q}}+\mathbf{K}_{f}\mathbf{q}+\mathbf{C}_{f}\overset{.}{\mathbf{q}}+\mathbf{C}_{q}^{\mathrm{T}}{\unicode[Arial]{x03BB}} =\mathbf{Q} \end{equation}

In Eq. (27), $\overset{..}{\mathbf{q}}$ represents the second-order derivative of the system configuration variable with respect to time, $\mathbf{C}_{\mathbf{q}}^{\mathrm{T}}$ denotes the derivative of the system constraint equation with respect to the system configuration variable, and $\boldsymbol{\lambda }$ signifies the Lagrangian multiplier, and there is:

(28) \begin{equation} \left\{\begin{array}{l} \mathbf{K}_{f}=\left[\begin{array}{l@{\quad}l@{\quad}l} \mathbf{0}_{27\times 27} & \mathbf{0}_{27\times 30} & \mathbf{0}_{27\times 6}\\ \mathbf{0}_{30\times 27} & \mathbf{K}_{ff}^{6} & \mathbf{0}_{30\times 6}\\ \mathbf{0}_{6\times 27} & \mathbf{0}_{6\times 30} & \mathbf{0}_{6\times 6} \end{array}\right]_{63\times 63}\\[3pt] \mathbf{C}_{f}=\left[\begin{array}{l@{\quad}l@{\quad}l} \mathbf{0}_{27\times 27} & \mathbf{0}_{27\times 30} & \mathbf{0}_{27\times 6}\\ \mathbf{0}_{30\times 27} & \mathbf{C}_{ff}^{6} & \mathbf{0}_{30\times 6}\\ \mathbf{0}_{6\times 27} & \mathbf{0}_{6\times 30} & \mathbf{0}_{6\times 6} \end{array}\right]_{63\times 63} \end{array}\right. \end{equation}

The damping matrix of the flexible body is denoted as $\mathbf{C}_{ff}^{6}$ . The intricate nature of damping mechanisms presents challenges in accurately quantifying damping levels in practical analyses. Typically, the damping matrix of a flexible body is derived through a linear combination of its inertia tensor and stiffness matrix [Reference Shabana36, Reference Trautt and Bayo37] as follows:

(29) \begin{equation} \mathbf{C}_{ff}^{6}=\alpha \mathbf{M}_{t}^{6}+\beta \mathbf{K}_{ff}^{6} \end{equation}

The determination of the Lagrangian multiplier ${\unicode[Arial]{x03BB}}$ in Eq. (27) poses challenges due to the presence of second-order differential equations and algebraic equations within the equation, hindering its solvability. The system dynamics equations are simplified in this paper using the independent coordinate method [Reference Yuwei, Bin, Haibo, Xinhua, Yang, Liang and Lei35]. The independent variables and configuration variables of the system can be determined from Eqs. (9) and (10), respectively. This is followed by a reorganization of the configuration variables of the system, which can be obtained as follows:

(30) \begin{equation} \mathbf{q}_{n}=\left[\begin{array}{l@{\quad}l} \mathbf{q}_{i} & \mathbf{q}_{d} \end{array}\right] \end{equation}

By employing system-independent variables, Eq. (27) can be simplified to:

(31) \begin{equation} \mathbf{M}_{ni}\,\ddot{\mathbf{q}}_{i}=\mathbf{Q}_{i} \end{equation}

In Eq. (31), $\mathbf{M}_{ni}$ represents the inertia tensor associated with the system’s independent variables, $\mathbf{Q}_{i}$ denotes the system’s generalized force and the quadratic velocity inertia force linked to $\mathbf{q}_{i}$ and $\dot{\mathbf{q}}_{i}$ , and $\ddot{\mathbf{q}}_{i}$ signifies the second-order derivative of the system’s independent variables with respect to time. By defining the state vectors $\mathbf{Y}=[\begin{array}{l@{\quad}l}\mathbf{{q}_{i}}^{\mathrm{T}} & \dot{\mathbf{q}}_{_{i}}^{\mathrm{T}} \end{array}]^{\mathrm{T}}$ and $\dot{\mathbf{Y}}=[\begin{array}{l@{\quad}l} \dot{\mathbf{q}}_{_{i}}^{\mathrm{T}} & \ddot{\mathbf{q}}_{_{i}}^{\mathrm{T}} \end{array}]^{\mathrm{T}}$ as per Eq. (10) and subsequently substituting them into Eq. (31) with necessary adjustments, we are able to derive the following result:

(32) \begin{equation} \overset{.}{\mathbf{Y}}=f\!\left(\mathbf{Y},t\right) \end{equation}

The simplified system dynamics equations consist of a set of first-order differential equations. By providing the initial values of the independent variables ( $\mathbf{{q}}_{\boldsymbol{i}}$ ) and solving the equations within the simulation time interval, one can derive the evolving positions, velocities, and accelerations of the system components over time.

4.1. The biomimetic load-bearing comfort of the exoskeleton

The biomimetic load-bearing comfort of the exoskeleton depends on the kinematic compatibility (reflected in the matching of the exoskeleton to the natural motion posture of the human body without affecting it) and its dynamic load-bearing effectiveness. Among them, the dynamic load-bearing effectiveness of the exoskeleton is mainly reflected in the magnitude of the equivalent force on the waist during human movement. The smaller the equivalent force, the better the exoskeleton’s motion-carrying performance. Meanwhile, the kinematic compatibility of the exoskeleton is primarily discernible in the discrepancy between the human trajectories prior to and following the exoskeleton’s utilization during human movement. The lesser the discrepancy, the higher the calibration of the motion bionic matching of the exoskeleton. Based on the above multi-body system dynamics model of the human-exoskeleton parallel wearable equivalent research prototype, the equivalent flexible beam at the waist of the human body during bending motion experiences a flexible force of:

(33) \begin{equation} \mathbf{F}_{d}^{j}=\mathbf{K}_{ff}^{6j}\mathbf{q}^{6j}+\mathbf{C}_{ff}^{6j}\dot{\mathbf{q}}^{6j}=\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \!\mathbf{0}_{1\times 3} & F_{rx}^{6j} & F_{ry}^{6j} & F_{rM}^{6j} \end{array}\right]_{1\times 6}^{\mathrm{T}} \end{equation}

In Eq. (33), $j$ represents the discrete unit number that corresponds to the lumbar spine of the $Lj\sim L(j+1)$ segment within the human body. $F_{rx}^{6j} , F_{ry}^{6j}$ , and $F_{rM}^{6j}$ denote the axial force, tangential force, and torque exerted on the flexible beam at the waist, respectively. Transforming these forces to the global coordinate system yields:

(34) \begin{equation} \left[\begin{array}{l@{\quad}l@{\quad}l} F_{dx}^{6j} & & F_{dy}^{6j} \end{array}\right]_{1\times 2}^{\mathrm{T}}=\mathbf{A}^{6j}\left[\begin{array}{l@{\quad}l@{\quad}l} F_{rx}^{6j} & & F_{ry}^{6j} \end{array}\right]_{1\times 2}^{\mathrm{T}} \end{equation}

In Eq. (34), $F_{dx}^{6j}$ and $F_{dy}^{6j}$ represent the forces corresponding to the horizontal and vertical orientations of the waist-equivalent flexible beam within the global coordinate system. As the direction of the force from the body’s weight and handling loads is perpendicular to the ground, the average force in the vertical direction was used as the equivalent comprehensive force on the waist; that is, the evaluation index of the exoskeleton’s dynamic load capacity can be defined as:

(35) \begin{equation} F_{d}=\sum _{j=1}^{5}F_{dy}^{6j} \end{equation}

In Fig. 3(a), $\mathbf{F}_{k}^{62}-\mathbf{F}_{k}^{92}$ represent the redundant human-exoskeleton interaction forces the exoskeleton may exert on the waist during movement, while $\mathbf{F}_{k}^{12}-\mathbf{F}_{k}^{42}$ indicate the exoskeleton’s structural spring damping forces. Both types of forces can impede movement, increasing physical burden and affecting natural posture, thereby reducing comfort of wearing the exoskeleton. The exoskeleton’s drive relies on the energy generated by the body’s metabolism, which is expressed as the work done to overcome the structure spring damping force of the exoskeleton during movement; thus, lower energy consumption implies greater dexterity and comfort of the exoskeleton. This paper uses the redundant interaction force magnitude and the work needed to overcome structural damping forces to quantify hindrance effect on human motion. The effect on the natural movement posture of the human body is described by the difference in the human trajectory before and after wearing the exoskeleton. Consequently, the biomimetic cooperation index of the exoskeleton can be described as follows:

(36) \begin{equation} \left\{\begin{array}{l} F_{s}=\sum\limits _{i=6}^{9}\left| \mathbf{F}_{k}^{i}\right| \\[3pt] E_{c}=\sum\limits _{i=1}^{4}\left| \mathbf{F}_{k}^{i}\right| \times {\unicode[Arial]{x0394}} \left| \mathbf{r}^{ai}-\mathbf{r}^{bi}\right| \\[3pt] R_{z}=\sum\limits _{j=1}^{5}\left(\left| r_{x}^{5j}-r_{xo}^{5j}\right| +\left| r_{y}^{5j}-r_{yo}^{5j}\right| \right) \end{array}\right. \end{equation}

In Eq. (36), $r_{xo}^{5j}$ and $r_{yo}^{5j}$ denote the x and y coordinates of discrete element $j$ of component 5 in the global coordinate system when the exoskeleton is not worn. To achieve effective achievement of bionic design for the exoskeleton, and considering human-exoskeleton harmony, comprehensive bionic performance optimization is deeply carried out in this paper. This optimization is based on the dynamics loading indices, wearing comfort indices, and Eq. (32). Because it is difficult to reflect the dynamic characteristics of the system by directly calculating the index value at the end of the simulation time, uniform time elements with an amount of $n$ are used to discretize the simulation time, and the indices are calculated at the end of each time element. By accumulating these indices, global performance indices reflecting the system’s dynamic characteristics are obtained :

(37) \begin{equation} \min \left\{\begin{array}{c} F_{d}=\sum\limits _{i=1}^{n}F_{d}^{i},F_{s}=\sum\limits _{i=1}^{n}F_{s}^{i}\\[3pt] E_{c}=\sum\limits _{i=1}^{n}E_{c}^{i},R_{z}=\sum\limits _{i=1}^{n}R_{z}^{i} \end{array}\right. \end{equation}

Where Eq. (37) is the optimal target function of exoskeleton bionic comprehensive performance optimization, and $F_{d} , F_{s} , E_{c}$ , and $R_{z}$ are indices of the equivalent comprehensive force on the waist, the redundant human-exoskeleton interaction force, the energy consumption of the exoskeleton, and the difference in the human body’s trajectory, respectively. Furthermore, in order to more intuitively evaluate the kinematic compatibility and dynamic load-bearing effectiveness of the exoskeleton, this paper proposes a dynamics performance evaluation index for optimizing the dynamic characteristics of human-exoskeleton wearable parallel systems, termed biomimetic load-bearing comfort, defined as:

(38) \begin{equation} KF=\mu _{1}\times K_{{F_{d}}}+\mu _{2}\times K_{{F_{s}}}+\mu _{3}\times K_{{E_{c}}}+\mu _{4}\times K_{{R_{z}}} \end{equation}

In Eq. (38), $K_{{F_{d}}} , K_{{F_{s}}} , K_{{E_{c}}}$ , and $K_{{R_{z}}}$ represent the assessment criteria for the equivalent force on the human waist, the force of human-exoskeleton interaction, trajectory variance, and overall energy expenditure of the exoskeleton, respectively. The weights assigned to these indicators are denoted by $\mu _{i}(i=1,2,3,4)$ . And there is:

(39) \begin{equation} \sum _{i=1}^{4}\mu _{i}=1 \end{equation}

The reason for adopting this weighted summation method to calculate the evaluation indexes is that we seek to optimize the overall performance of the exoskeleton’s dynamic characteristics. This approach accounts for the mutual influences and constraints among performance indicators, reflecting the exoskeleton’s performance from various perspectives. By assigning appropriate weights to each index, we can assess their impact on overall performance and achieve global optimization.

Nevertheless, during the calculation process, if the order of magnitude of an indicator differs significantly from that of the other indicators, it will have a detrimental impact on the final optimization results. Hence, during the weight assignment process, we comprehensively consider the importance and influence of each basic item in the evaluation system. Subsequently, we conduct appropriate normalization procedures to guarantee the comparability of indicators during the calculation. By employing this thorough performance evaluation approach, we can accurately assess the overall performance of the exoskeleton in an objective manner, mitigating the undue impact of any single indicator on the final evaluation outcome. This method ensures the scientificity and reliability of the evaluation results. Not only does this method assist in accurately assessing the dynamic characteristics of the exoskeleton but it also establishes a theoretical basis for optimizing and improving exoskeleton design. Adaptive effects are difficult to achieve due to the difficulty of controlling order-of-magnitude differences using the size of the weights. In order to eliminate the effects of order of magnitude, this study conducts dimensionless processing on four performance indicators. The effect of optimizing individual metrics is embodied by calculating the percentage in the reduced value of the metric after optimization relative to the value before optimization or before the exoskeleton was worn. As shown in Eq. (40).

(40) \begin{equation} \left\{\begin{array}{l} K_{{F_{d}}}=\dfrac{F_{d}^{ex}-F_{d}}{F_{d}^{ex}}\\[3pt] K_{{F_{s}}}=\dfrac{F_{s}^{ex}-F_{s}}{F_{s}^{ex}}\\[3pt] K_{{E_{c}}}=\dfrac{E_{c}^{ex}-E_{c}}{E_{c}^{ex}}\\[3pt] K_{{R_{z}}}=\dfrac{\left(\left(R_{z}^{ex}-R_{z}^{o}\right)-\left(R_{z}-R_{z}^{o}\right)\right)}{R_{z}^{ex}-R_{z}^{o}} \end{array}\right. \end{equation}

In Eq. (40), each performance indicator is denoted by a $ex$ superscript indicating the pre-optimization value, a $o$ superscript indicating the value before wearing the exoskeleton, and no superscript indicating the post-optimization value. The analysis of Eqs. (37) and (38) reveals a direct relationship between the biomimetic load-bearing comfort of the exoskeleton and various biomechanical factors. Specifically, a higher comfort level correlates with increased axial force on the human waist, greater resistance needed to move with the exoskeleton, greater compression force on the waist, and minimal alteration in posture. These findings suggest that enhanced the kinematic compatibility and dynamic load-bearing effectiveness are associated with exoskeleton performance. Additionally, the results suggest that optimizing these biomechanical factors can lead to better overall exoskeleton performance.

4.2. The biomimetic load-bearing comfort of the exoskeleton

Given the challenge of customizing exoskeleton structural parameters to accommodate individual variations in wearers and the convenience of replacing structural springs, the stiffness coefficients $K^{1}-K^{4}$ of exoskeleton structural springs are chosen as system-dimension parameters. By combining Eqs. (32) and (38), it becomes evident that the exoskeleton’s biomimetic load-bearing comfort follows a nonlinear, strongly coupled model with multiple inputs and outputs, requiring the solution of differential equations. Traditional analytical methods struggle to derive the global optimal solution. This paper employs the method of multiple parameter combination optimization to enhance the dynamic performance of the exoskeleton [Reference Ping, Yihao, Peitao and Bin38], focusing on optimizing the biomimetic load-bearing comfort. By utilizing optimization algorithms, this study aims to identify the optimal dimensional parameters of the exoskeleton to improve its kinematic compatibility and dynamic load-bearing effectiveness. The dynamic dimensional synthesis problem in exoskeletons involves optimizing the dimension parameters $K^{1}-K^{4}$ to enhance the overall dynamic performance, as outlined in reference [Reference Li, Mei, Liu and Huang39].

The multidimensional space search method involves selecting parameters within specified optimization parameter ranges using a predetermined step size to create a new parameter value space. These parameters within each newly defined space are then systematically integrated to optimize the objective function and identify the global optimal value. This paper explores the range of exoskeleton dimension parameters as the parameter space, with the biomimetic load-bearing comfort ( $KF$ ) serving as the optimization objective function. To enhance the precision of the global optimal value determination, the algorithm should iterate multiple times, progressively decreasing the step size in each iteration to achieve a finer parameter value space. When the discrepancy between the optimal values derived from successive iterations falls below the predetermined threshold, the algorithm is deemed to have converged, and the resultant optimal value represents the global optimum solution.

Although multidimensional space search is straightforward and practical, it suffers from drawbacks such as low algorithm efficiency and limited global search capability, increasing the risk of becoming trapped in local optima. To achieve the best possible outcomes, this paper utilizes a multi-algorithm parallel optimization technique. Specifically, it employs the particle swarm optimization algorithm for optimization and assesses the optimal values of system dimension parameters. Employing the particle swarm optimization algorithm to explore the global optimal value of the biomimetic load-bearing comfort ( $KF$ ) evaluation index for enhancing exoskeleton performance as introduced earlier, the exoskeleton dimension parameter range serves as the particle search space. The algorithm initiates the process by generating $m$ random particles within the particle search space. These particles represent the initial population and possess a dimension of $d$ , which aligns with the number of system optimization parameters. The equation employed for updating the particle velocity and position in the iteration process is as follows:

(41) \begin{align} v_{id}^{\left(k+1\right)}&=\lambda v_{id}^{\left(k\right)}+c_{1}r_{1}\!\left(P_{Lid}^{\left(k\right)}-x_{id}^{\left(k\right)}\right)+c_{2}r_{2}\left(P_{Gid}^{\left(k\right)}-x_{id}^{\left(k\right)}\right)\nonumber\\ x_{id}^{\left(k+1\right)}&=x_{id}^{\left(k\right)}+v_{id}^{\left(k+1\right)} \end{align}

In Eq. (41), $P_{Lid}^{(k)}$ represents the individual optimal value, while $P_{Gid}^{(k)}$ denotes the group optimal value. The learning factors are represented by $c_{1}$ and $c_{2}$ , where $r_{1}$ and $r_{2}$ denote random numbers within the range of [0, 1]. The parameters $k , v_{id}^{(k)}$ , and $x_{id}^{(k)}$ correspond to the number of iterations, speed, and position of the $i$ -th particle’s flight in the $k$ -th iteration, respectively.

The compression factor method is applied to refine the particle swarm algorithm. ${\unicode[Arial]{x03BB}}$ is the compression factor with:

(42) \begin{equation} {\unicode[Arial]{x03BB}} =\frac{2}{\left| 2-\varphi -\sqrt{\varphi ^{2}-4\varphi }\right| } \end{equation}

In Eq. (42), $\varphi =c_{1}+c_{2}$ . A diagram illustrating the system optimization process is shown in Fig. 4.

Figure 4.

Boundary conditions of finite element model.

4.3. Research methodology

Considering the encompassing kinematic compatibility, and dynamic load-bearing effectiveness of biomimetic coordination, we present a systematic research approach that integrates qualitative and quantitative analyses, as illustrated in Fig. 5:

Figure 5.

Flowchart of a systematic research approach based on the biomimetic coordination.

For qualitative analyses, we propose a novel unpowered rigid-flexible coupling waist exoskeleton with multiple degrees of freedom and its human-exoskeleton parallel wearable equivalent research prototype, based on the biomechanical properties of the lumbar vertebrae in the human body.

For quantitative analyses, we conduct a dynamic dimensional synthesis that integrates the multi-body dynamics model of the human-exoskeleton wearable multi-body dynamics models, comprehensive quantifiable indicators, optimization algorithms, and numerical analyses. Initially, based on experimental biomechanical data concerning the human lumbar region in the sagittal plane, a multi-body system dynamics model is meticulously developed for the aforementioned research prototype. Subsequently, a novel comprehensive optimization index, termed biomimetic load-bearing comfort, is proposed to effectively quantify the biomimetic coordination of the exoskeleton. Thereafter, based on the numerical simulation of system dynamics and the biomimetic load-bearing comfort index, the optimization algorithm of multi-dimensional spatial search and compression factor particle swarm in parallel is comprehensively applied to finally determine the global optimal dimension parameter group of the exoskeleton system. Finally, we conduct a comparative analysis of simulation data obtained before and after optimization.

4.4. System optimization results and numerical simulation

Tables III and IV display the initial conditions, simulation outcomes, and optimization results concerning the structural parameters and independent variables of the exoskeleton.

Table III.

Exoskeleton parameters and independent variable initial conditions.

Table IV.

Simulation and optimization results (two decimal places are retained).

A comparison of the equivalent integrated force acting on the waist at different time intervals with three distinct sets of dimension parameters (as displayed in Fig. 6) reveals a notable abrupt change in waist force between 1.4 s and 1.5 s when the stiffness of the four spring sets is closely matched (initial values 1 and 2, as outlined in Table IV). This shift aligns with the point at which the angle of the exoskeleton connects to the rod $\theta ^{4}\lt \theta ^{3}$ (as shown in Fig. 7), signifying concave dynamic instability within the exoskeleton configuration. Given that achieving non-dimensionality for each indicator using Eq. (40) necessitates utilizing pre-optimization outcomes, enhancing the optimization precision of system dynamic dimensional synthesis requires careful selection of initial parameter values to prevent the aforementioned scenarios. However, with an initial value 3, all domains exhibit $\theta ^{4}\gt \theta ^{3}\gt \theta ^{2}\gt \theta ^{1}$ , demonstrating the stability of the human-exoskeleton parallel waist exoskeleton system. Consequently, this study proceeds to perform dynamic dimensional synthesis using the results obtained under the initial value 3, as illustrated in Table IV.

Figure 6.

Equivalent combined force on the human lumbar region.

Figure 7.

Angle of each linkage of the exoskeleton.

Figure 8.

Redundant human-exoskeleton interaction forces.

Figure 9.

Total exoskeleton energy consumption.

Figure 10.

Exoskeleton cumulative trajectory difference.

Figures 6–10 illustrate the changes in the equivalent comprehensive force exerted on the waist, the angles of different components of the exoskeleton, the total redundant human-exoskeleton interaction force, the overall energy consumption of the exoskeleton, and the cumulative difference in human trajectory over time. In Fig. 7, the optimal dimension parameters obtained through multidimensional space search and particle swarm optimization demonstrate that the exoskeleton angles follow the order $\theta ^{4}\gt \theta ^{3}\gt \theta ^{2}\gt \theta ^{1}$ consistently throughout the entire domain. This consistent pattern signifies the stability of the human–exoskeleton parallel waist exoskeleton system, thereby validating the rational and efficient choice of the initial value 3. From Table IV, Figs. 6, and 8, it is clear that before optimization, the cumulative equivalent comprehensive force on the waist decreased by 1.30% under the parameters of initial value 3. After optimization, reductions in two sets of metrics were 1.25% and 1.22%, indicating improved dynamic load-bearing capabilities of the exoskeleton throughout the operational range. The dynamic load-bearing effect is superior before optimization; however, post-optimization, the redundant human-exoskeleton interaction force $F_{s}$ of the exoskeleton decreased significantly, by an order of magnitude. Despite the exoskeleton’s initial strong load-bearing capacity, kinematic compatibility was inadequate. Analysis of Figs. 9 and 10 shows that the optimized exoskeleton outperformed the pre-optimized version in cumulative energy consumption ( $E_{c}$ ) and cumulative human trajectory difference ( $R_{z}$ ), though the optimization effect was limited due to low weights assigned to these metrics. Kinematic compatibility under particle swarm optimization surpassed that of pre-optimization and multidimensional space search. Despite a slight compromise in dynamic load-bearing performance, the overall optimal dynamic performance of the exoskeleton is achieved using the particle swarm optimization parameters. Evaluating this performance objectively requires calculating the instantaneous biomimetic load-bearing comfort over simulation time. Thus:

(43) \begin{equation} KF^{i}=\mu _{1}\times K_{F_{d}}^{i}+\mu _{2}\times K_{F_{s}}^{i}+\mu _{3}\times K_{E_{c}}^{i}+\mu _{4}\times K_{R_{z}}^{i} \end{equation}

Figure 11.

Instantaneous trajectory difference.

In Eq. (43), $i$ denotes the serial number of the discrete time unit, $K_{F_{d}}^{i}$ signifies the axial force evaluation of the waist at the conclusion of the time unit, and $K_{E_{c}}^{i}$ represents the cumulative energy consumption evaluation index of the exoskeleton as of the $i$ -th discrete time unit.

Figure 11 illustrates the curve depicting the instantaneous difference in human trajectory, denoted as $R_{z}^{i}$ , as a function of time. According to Figs. 8 and 11, they are evident that both $F_{s}$ and $R_{z}^{i}$ exist zero points before optimization. Nonetheless, upon applying dimensionless processing as per Eq. (40), the pre-optimization outcome is utilized as the denominator, resulting in an infinite value. To circumvent this issue, the sum of indicator values from the preceding time interval is computed to prevent infinite calculations. To wit:

(44) \begin{equation} F_{s}^{i}=\sum _{n=1}^{i-1}F_{s}^{n}, \quad R_{z}^{i}=\sum _{n=1}^{i-1}R_{z}^{n} \end{equation}

Due to the immediate application of the driving force and load, the transient response prior to 0.2 s is disregarded, with a focus solely on the steady-state response during the subsequent simulation period. Figure 12 illustrates the temporal fluctuation in the comfort level resulting from the exoskeleton’s instantaneous biomimetic load-bearing comfort. Figure 12 clearly shows that the instantaneous biomimetic load-bearing comfort of the exoskeleton following particle swarm optimization exceeds 0. Moreover, the instantaneous biomimetic load-bearing comfort of the exoskeleton under these specific parameters significantly surpasses that achieved under alternative conditions, suggesting superior dynamic performance within this parameter framework. An analysis of the optimization outcomes detailed in Table IV supports the assertion that the particle swarm optimization parameters represent a globally optimal solution. The simulation results corroborate the reasonableness and effectiveness of the exoskeleton dynamic dimensional synthesis and its research methodology.

Figure 12.

Instantaneous biomimetic load-bearing comfort.

Based on Figs. 6–8, the presence of spring preload in exoskeleton structures causes the system to be impacted by the redundant human-exoskeleton interaction force at 0.4 s. This generates elastic vibrations in both the exoskeleton and the human body, resulting in oscillations in the equivalent combined force on the human waist. To guarantee the stability, reliability, and safety of the exoskeleton during human movement, it is imperative to perform vibration suppression analysis on these devices. It is therefore necessary to investigate the dynamic properties of the human-exoskeleton parallel wearable equivalent research prototype under different exoskeleton damping conditions. Figures 13, 14, and 15 depict the variations in the equivalent comprehensive force on the human waist, the redundant human-exoskeleton interaction force, and the exoskeleton energy consumption concerning the exoskeleton damping ratio. Figure 14 reveals a direct relationship between the damping ratio of the spring in the exoskeleton structure and the speed of oscillation convergence of the equivalent comprehensive force on the human waist. Nevertheless, it is noteworthy that the amplitude at the onset of oscillation also increases as the damping ratio increases. The data presented in Fig. 15 illustrate that an increase in the spring damping ratio of the exoskeleton structure leads to a higher equivalent stiffness of the exoskeleton, thereby amplifying the redundant human-exoskeleton interaction force imposed by the exoskeleton on the human body. Consequently, at higher damping ratios, the amplitude of the equivalent combined force applied to the human waist exceeds that experienced at lower damping ratios.

Figure 13.

Equivalent combined human lumbar force at different damping ratios.

Figure 14.

Redundant human-exoskeleton interaction forces across varying damping ratios.

Figure 15.

Energy consumption of the exoskeleton under varying damping ratios.

When damping is low, an inadequate equivalent stiffness of the exoskeleton results in the redundant human-exoskeleton interaction force at 1.8 s. An increase in the damping ratio, as depicted in Fig. 15, correlates with heightened energy consumption by the exoskeleton, thereby augmenting motion resistance experienced by the human body. Furthermore, the force of redundant human-exoskeleton interaction impact is heightened by higher damping ratios, resulting in a decrease in kinematic compatibility. Hence, the influence of the damping ratio on the exoskeleton’s comprehensive performance is significant and warrants attention. To fully understand how the damping ratio affects the integrated dynamic characteristics of the exoskeleton, it is crucial to examine how biomimetic load-bearing comfort in the exoskeleton changes with varying damping ratios.

The biomimetic load-bearing comfort of the exoskeleton, taking into account the damping ratio, is depicted in Fig. 16. Observationally, it is evident that the biomimetic load-bearing comfort of the exoskeleton experiences a noticeable shift in the vicinity of a damping ratio of 0.02, indicating that this is the critical value for the spring damping ratio of the exoskeleton structure. Below this threshold, the exoskeleton structure becomes unstable, leading to decreased kinematic compatibility due to heightened elastic vibration and redundant human-exoskeleton interaction force. Conversely, when the damping ratio exceeds this value, the exoskeleton structure remains stable; however, with an increasing damping ratio, the redundant human-exoskeleton interaction force rises to a significant level that cannot be ignored, consequently diminishing the kinematic compatibility of the exoskeleton. In conclusion, the spring damping ratio of the exoskeleton structure should be taken as a value greater than 0.02, and considering that the vibration of the system converges slowly at low damping ratios and the redundant human-exoskeleton interaction force occurs at 1.8 s, a slight increase in the exoskeleton structure’s spring damping ratio is essential to enhance the kinematic compatibility.

Figure 16.

Biomimetic load-bearing comfort taking into account the damping ratio.