2. Structural prototype and its biomimetic functional mechanism

The human ankle is commonly described as a three-degree-of-freedom joint complex, comprising dorsiflexion–plantarflexion in the sagittal plane, adduction–abduction in the transverse plane, and pronation–supination as composite multiplanar motions, as illustrated in Figure 1 [Reference Conconi, Sancisi, Leardini and Belvedere27]. These motions do not occur independently during gait. Plantarflexion, for example, is often accompanied by eversion, reflecting the coupled nature of ankle motion. The resulting kinematic pattern also varies considerably among individuals because of differences in bone morphology and orientation, articular-surface geometry, and the constraints imposed by ligaments, joint capsules, tendons, and other soft tissues.

Figure 1.

Schematic diagram of human ankle-joint motion.

Drawing on these anatomical and biomechanical characteristics, a passive compliant ankle-exoskeleton concept was developed, together with the structural demonstration prototype shown in Figure 2(b) and (c). (The structural demonstration prototype was fabricated to examine the proposed structural arrangement, manufacturability, and assembly feasibility. It was not used to validate motion performance, wearing adaptability, stiffness, damping, energy storage, interaction force, output torque, or user comfort. Its presence should therefore not be interpreted as experimental evidence of functional performance.) Rather than relying on a single assistance element, the proposed design combines several biomimetic subsystems to accommodate passive compliant motion during the prescribed stance-related stepping condition. These subsystems comprise (I) a soleus muscle-inspired elastic-damping energy storage assistance mechanism; (II) an ankle soft tissue-mimetic elastic load-bearing structure; (III) a foot–ankle complex motion-inspired buffering-propulsion unit; and (IV) an ankle-joint motion-inspired passive compliant support.

Figure 2.

Structural concept of the bioinspired passive compliant ankle exoskeleton: (a) human lower-limb musculoskeletal model; (b) overall exoskeleton structure; (c) structural demonstration prototype; and (d) three-dimensional models of the biomimetic subsystems: I, soleus muscle-inspired elastic-damping energy storage assistance mechanism; II, ankle soft tissue-mimetic elastic load-bearing structure; III, foot–ankle complex motion-inspired buffering-propulsion unit; and IV, ankle-joint motion-inspired passive compliant support.

From the perspective of local motion mechanisms, the primary motion of the ankle joint is coordinated by a pair of antagonistic muscle groups: the plantarflexor group, mainly including the gastrocnemius and soleus, and the dorsiflexor group, mainly including the tibialis anterior. Through coordinated contraction and relaxation, these two muscle groups enable flexible joint actuation and stable control. The human musculoskeletal system exhibits good compliance, efficient energy conversion, and strong environmental adaptability [Reference Herr and Grabowski28]. To reduce muscle load and improve structural compliance, existing studies have commonly employed elastic elements, overrunning clutches, or low-power actuation structures to achieve push-off assistance in ankle–foot systems [Reference Collins, Wiggin and Sawicki24, Reference Mazzarini, Fantozzi, Papapicco, Fagioli, Lanotte, Baldoni, Dell’Agnello, Ferrara, Ciapetti, Molino Lova, Gruppioni, Trigili, Crea and Vitiello29]. However, such mechanical clutch mechanisms usually suffer from structural complexity, relatively large mass, and noticeable motion inertia, which limit their lightweight integration and application adaptability in wearable devices. To address these limitations, this study proposes a soleus muscle-inspired elastic-damping energy-storage assistance mechanism. As illustrated in Figure 2(d-I), the core of this subsystem is a rigid-flexible-soft composite mechanism engineered for energy-storage-enabled geometric accommodation. In the present study, the subsystem is represented only by its geometric and kinematic constraints, and the analysis is limited to its motion-transmission behavior under the prescribed stepping condition. Its mass, force–displacement response, energy-storage efficiency, and impact-attenuation performance are not evaluated. These aspects will be investigated in future work through dynamic modeling, bench testing, and experimental validation.

On the other hand, from the perspective of foot motion mechanisms, the ankle joint mainly contributes to mechanical work output, whereas the foot structure primarily dissipates impact energy [Reference Takahashi, Worster and Bruening30, Reference Isvilanonda, Li, Williams, Cavanagh, Haynor, Chu and Ledoux31]. It is noteworthy that ankle soft tissues can partially absorb and dissipate impact energy owing to their damping characteristics, and this function is independent of the contraction of specific antagonistic muscle groups [Reference Dewolf, Ivaniski-Mello, Peyre-Tartaruga and Mesquita32–Reference Zhuang, Gu, Li, Shen, Liu, Li, Wang, Wang, Liu and Ba34]. Meanwhile, considering the flexible load-bearing characteristics of ankle soft tissues, an ankle soft tissue-mimetic elastic load-bearing structure is proposed, as shown in Figure 2(d-II). In terms of structural configuration, this unit serves to accommodate posture-related geometric variations and mitigate the risk of constraint interference. Building on this energy-metabolism characteristic of the foot–ankle complex, this study proposes a foot–ankle complex motion-inspired buffering–propulsion unit, as shown in Figure 2(d-III). The mechanism follows a modular design concept and incorporates a dual-stage spring structure. Through sequential compression of the dual-stage spring structure, the system is configured to provide displacement buffering during initial ground contact and support propulsion-related passive motion during late push-off, thereby improving the kinematic continuity of the mechanism during the analyzed stance-related motion cycle. Specifically, the primary spring in the outer tube provides initial displacement buffering during early ground contact; when the inner tube reaches its travel limit, the secondary spring is engaged, forming a progressive stroke-limited geometric response. This design is intended to improve geometric accommodation during the analyzed stance-related stepping condition.

In view of the multi-degree-of-freedom coupled motion characteristics of the ankle joint during the gait cycle [Reference Conconi, Sancisi, Leardini and Belvedere27], this study designs an ankle-joint motion-inspired passive compliant support, as shown in Figure 2(d-IV), with the aim of providing a structural basis for compliant ankle motion, although the quantitative analysis in this study is limited to sagittal-plane dorsiflexion–plantarflexion-related kinematics. Through the frictional rolling and adaptive sliding interactions between rollers and guides in the sagittal and frontal planes, this module is designed to provide a potential geometric compliance path for ankle dorsiflexion–plantarflexion and pronation-supination motions and to reduce additional geometric constraints that may arise during movement. The structure is intended to provide necessary structural support while preserving, as far as possible, the degrees of freedom and flexibility required for physiological ankle motion, thereby offering a structural basis for improving human–exoskeleton motion coordination and motion adaptability under varying foot–ground contact conditions. Given that the primary research scope of this work is centered on sagittal-plane kinematic characteristics, all subsequent numerical analyses are conducted predominantly within this plane. The ability of the proposed structure to accommodate full three-dimensional coupled ankle motion will be systematically investigated as an independent research direction in future work.

In summary, the proposed bioinspired passive compliant ankle exoskeleton does not rely on a single assistive unit to realize its intended functions. Instead, the structural demonstration prototype is designed to support passive and compliant motion during the stance-related stepping condition through the structural coupling of multiple biomimetic subsystems. Specifically, the roller-guide-rail structure is designed to provide compliant motion allowance for ankle dorsiflexion–plantarflexion and pronation-supination. The two-stage spring structure is designed to provide displacement accommodation during landing-related and late-stance-related motions. The rigid-rope-elastic-rope-extension-spring assembly is used to accommodate large relative displacements between components during gait-related motion. The elastic load-bearing structure is used to emulate the flexible support and posture-adaptation characteristics of ankle soft tissues at the structural-design level. Without external energy input, these structures provide a geometric and structural basis for investigating passive support, motion transmission, and configuration adaptability of the exoskeleton during the stance-related stepping condition. The following model represents these subsystems according to their geometric roles and does not evaluate their mechanical properties.

A complete gait cycle consists of two phases: the swing phase and the stance phase. The ankle joint plays a key role in motion regulation during this process. Through coordinated postural changes during push-off, it contributes to the temporal matching between ground contact and propulsion. The design objective of the exoskeleton system is to match the primary motion characteristics of the human lower limb during the stance phase, thereby reducing geometric constraint conflicts in human–exoskeleton motion transmission.

Since the exoskeleton remains in a follower or unlocked state during the swing phase and does not impose additional motion constraints, this study focuses on a prescribed stance-related stepping condition, as shown in Figure 3. In this study, the analyzed condition is not intended to represent the entire gait cycle; rather, it is used to describe the main loading and unloading process associated with foot–ground contact, during which the ankle undergoes alternating dorsiflexion–plantarflexion. This condition is selected to investigate the kinematic transmission relationship between human input and exoskeleton response, as well as the smoothness of the passive compliant response.

To clearly distinguish the structural design objectives from the analytical scope of the adopted kinematic model, the intended physical functions, simplified kinematic formulations, and unvalidated performance attributes of the compliant and elastic subsystems are elaborated as follows. The kinematic model established in this work is primarily developed to investigate geometric constraint relations, configuration-transmission characteristics, and interval variations of key kinematic responses under a prescribed stance-related stepping condition. At the structural design level, components including springs, elastic ropes, compliant load-bearing units, buffering-propulsion modules, and roller-guide-rail assemblies are designed to deliver elastic energy storage, compliant connection, load transfer, buffering-related motion accommodation, rolling guidance, and passive motion support. Nevertheless, no detailed mechanical model is constructed for these components in terms of stiffness, damping, mass, inertia, contact force, friction, or force–displacement relationships. Instead, they are equivalently represented as stroke-limited prismatic joints, variable-length kinematic elements, ideal guiding constraints, or angular/displacement response variables within the current kinematic framework. Accordingly, the results of this study cannot be directly applied to evaluate elastic force output, damping performance, stiffness characteristics, energy storage efficiency, impact attenuation potential, nominal design-stage load-bearing function, human–exoskeleton interaction force, or practical assistive effectiveness. All the aforementioned mechanical and functional performances require further systematic investigation via dynamic modeling and experimental validation.

Figure 3.

Analysis of the prescribed stance-related stepping condition.

Figure 3. Long description

A diagram of the stepping sequence during walking. The diagram shows the phases of dorsiflexion and plantarflexion, with labeled steps including Complete Heel-off, Heel Strike, Foot Flat, Heel Off, and Complete Heel-off. The analyzed stepping sequence is divided into First Dorsiflexion, First Plantarflexion, Second Dorsiflexion, and Second Plantarflexion. The energy release and storage processes are also indicated, with Dorsiflexion Energy Release, Plantarflexion Energy Release, Plantarflexion Energy Release or Weakness, Dorsiflexion Energy Release, and Plantarflexion Energy Storage. The diagram visually represents the biomechanical roles of the ankle joint in generating and regulating propulsion, maintaining gait stability, and adjusting step frequency.

3. Sagittal-plane human–exoskeleton–environment equivalent model

3.1. Sagittal-plane kinematic model of the exoskeleton

Although the human ankle exhibits multi-degree-of-freedom coupled motion, the present study focuses only on the sagittal-plane human–exoskeleton motion-transmission relationship in the dorsiflexion–plantarflexion direction under a prescribed stance-related stepping condition. Accordingly, a human–exoskeleton–environment equivalent kinematic model with prescribed foot–ground boundary constraints is established. These expressions refer to the same human–exoskeleton–environment equivalent multibody kinematic framework but emphasize different aspects of the model, namely its sagittal-plane scope, human–exoskeleton motion mapping, or prescribed foot–ground boundary constraints. The model describes the geometric mapping relationships among the prescribed human input, exoskeleton configuration responses, and foot–ground support constraints. In this study, the term “environment” refers exclusively to the prescribed foot–ground support constraints and boundary conditions associated with the prescribed stance-related stepping condition. Terrain geometry, friction, normal contact force, impact response, ground reaction force and environmental dynamics are not incorporated into the present model. Since the exoskeleton remains in a follower or unlocked state during the swing phase and does not impose additional motion constraints, this study mainly analyses the configuration changes and key motion responses induced by human input during the representative stance-related stepping condition. The equivalent kinematic model established based on the bioinspired passive compliant ankle exoskeleton structural demonstration prototype is shown in Figure 4. The human–exoskeleton coordination equivalent model based on multibody kinematics [Reference Shabana35] is constructed, as shown in Figure 4(b), where $O^{\mathrm{0}}X^{\mathrm{0}}Y^{\mathrm{0}}$ represents the global coordinate system and $O^iX^iY^i$ denotes the body-fixed coordinate system of component $i$ = 1,2, … , 14. The position vector ${\boldsymbol{R}}^{i}$ of any point on component $i$ can be described as follows:

(1) \begin{equation} {\boldsymbol{r}}^{i}={\boldsymbol{R}}^{i}+{\boldsymbol{A}}^{i}\overline {{\boldsymbol{u}}}^{i} \end{equation}

where ${\boldsymbol{R}}^{i}$ is the position vector of the origin of the body-fixed coordinate system of component $i$ relative to the origin of the global coordinate system; $\overline {{\boldsymbol{u}}}^{i}$ is the position vector of any point on component $i$ measured in its body-fixed coordinate system $O^iX^iY^i$ ; ${\boldsymbol{A}}^{i}$ is the pose transformation matrix of $O^iX^iY^i$ relative to the global coordinate system:

(2) \begin{equation} {\boldsymbol{A}}^{i}= \begin{bmatrix} \cos \theta ^{i} & -\sin \theta ^{i}\\[5pt] \sin \theta ^{i} & \cos \theta ^{i} \end{bmatrix} \end{equation}

where $\theta ^{i}$ is the rotation angle of the body-fixed coordinate system of component $i$ relative to the global coordinate system. In summary, the overall system of constraint equations for the exoskeleton is

(3) \begin{equation} \boldsymbol{C} = \left [ \boldsymbol{C}_{\text{1}}^{\mathrm{T}} \quad \boldsymbol{C}_{\text{2}}^{\mathrm{T}} \quad \boldsymbol{C}_{\text{3}}^{\mathrm{T}} \right ]_{\text{1}\times \text{41}}^{\mathrm{T}} \end{equation}

in the equation, $\boldsymbol{C}_{\text{1}}$ represents the kinematic pair constraints of the exoskeleton system, specifically expressed as follows:

(4) \begin{equation} {\boldsymbol{C}}_{1}= \left [ \begin{array}{c} {\boldsymbol{R}}^{{\text{1e}}}+{\boldsymbol{A}}^{\text{ 1e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle B}^{\text{ 1e}}-{\boldsymbol{R}}^{\text{2e}}-{\boldsymbol{A}}^{\text{ 2e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle B}^{\text{2e}}\\[5pt] {\boldsymbol{R}}^{\text{2e}}+{\boldsymbol{A}}^{\text{ 2e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle C}^{\text{ 2e}}-{\boldsymbol{R}}^{\text{3e}}-{\boldsymbol{A}}^{\text{ 3e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle C}^{\text{3e}}\\[5pt] {\boldsymbol{R}}^{\text{5e}}+{\boldsymbol{A}}^{\text{ 5e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle F}^{\text{ 5e}}-{\boldsymbol{R}}^{\text{6e}}-{\boldsymbol{A}}^{\text{ 6e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle F}^{\text{6e}}\\[5pt] {\boldsymbol{R}}^{\text{6e}}-{\boldsymbol{A}}^{\text{ 6e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle G}^{\text{ 6e}}-{\boldsymbol{R}}^{\text{7e}}+{\boldsymbol{A}}^{\text{ 7e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle G}^{\text{7e}}\\[5pt] {\boldsymbol{R}}^{\text{8e}}+{\boldsymbol{A}}^{\text{ 8e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle H}^{\text{ 8e}}-{\boldsymbol{R}}^{\text{9e}}-{\boldsymbol{A}}^{\text{ 9e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle H}^{\text{9e}}\\[5pt] {\boldsymbol{R}}^{\text{1e}}+{\boldsymbol{A}}^{\text{ 1e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle B}^{\text{ 1e}}-{\boldsymbol{R}}^{1\text{2e}}-{\boldsymbol{A}}^{1\text{2e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle B}^{1\text{2e}}\\[5pt] {\boldsymbol{R}}^{\text{2e}}+{\boldsymbol{A}}^{\text{ 2e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle C}^{\text{ 2e}}-{\boldsymbol{R}}^{\text{10e}}-{\boldsymbol{A}}^{\text{ 10e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle C}^{\text{10e}}\\[5pt] {\boldsymbol{R}}^{\text{10e}}+{\boldsymbol{A}}^{\text{ 10e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle I}^{\text{ 10e}}-{\boldsymbol{R}}^{1\text{1e}}\\[5pt] {\boldsymbol{R}}^{\text{3e}}+{\boldsymbol{A}}^{\text{ 3e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle D}^{\text{ 3e}}-{\boldsymbol{R}}^{\text{9e}}-{\boldsymbol{A}}^{\text{ 9e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle D}^{\text{9e}}\\[5pt] {\boldsymbol{R}}^{\text{4e}}+{\boldsymbol{A}}^{\text{ 4e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle E}^{\text{ 4e}}-{\boldsymbol{R}}^{1\text{4e}}\\[5pt] {\boldsymbol{R}}^{\text{3e}}+{\boldsymbol{A}}^{\text{ 3e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle F}^{\text{313e}} \end{array} \right ] = \mathbf{0}_{22}\times \mbox{1} \end{equation}

in the equation, ${\boldsymbol{C}}_{2}$ represents the component pose constraints of the exoskeleton system, specifically expressed as:

(5) \begin{equation} {\boldsymbol{C}}_{2}= \left [ \begin{array}{c} \theta ^{\text{1e}}-\theta ^{\text{2e}}\\[5pt] \theta ^{\text{4e}}-\theta ^{\text{5e}}\\[5pt] \theta ^{\text{7e}}-\theta ^{\text{8e}}\\[5pt] \theta ^{\text{3e}}-\theta ^{1\text{3e}}\\[5pt] \theta ^{1\text{1e}}-\theta ^{1\text{2e}}\\[5pt] \theta ^{13re}-\theta ^{1\text{4e}}\\[5pt] \theta ^{\text{2e}}-\theta ^{\text{3e}}-\phi ^{\text{1e}}\\[5pt] \theta ^{\text{3e}}-\theta ^{\text{9e}}-\phi ^{\text{2e}}\\[5pt] \theta ^{\text{5e}}-\theta ^{\text{6e}}-\phi ^{\text{3e}}\\[5pt] {\boldsymbol{R}}^{1\text{3e}}+\overline {{\boldsymbol{u}}}_{\scriptscriptstyle E}^{1\text{3e}}-{\boldsymbol{R}}^{1\text{4e}}-\overline {{\boldsymbol{u}}}_{\scriptscriptstyle E}^{1\text{4e}}\\[5pt] {\boldsymbol{R}}^{\text{7e}}+{\boldsymbol{A}}^{\text{7e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle G}^{7\text{8e}}-{\boldsymbol{R}}^{\text{8e}}-{\boldsymbol{A}}^{\text{8e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle H}^{\text{8e}}\\[5pt] {\boldsymbol{R}}^{\text{4e}}+{\boldsymbol{A}}^{\text{4e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle E}^{4\text{5e}}-{\boldsymbol{R}}^{\text{5e}}-{\boldsymbol{A}}^{\text{5e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle F}^{\text{5e}}\\[5pt] r_{2}^{1\text{1e}}-r_{2}^{1\text{2e}}-\tan \theta ^{1\text{2e}}\left (r_{1}^{1\text{1e}}-r_{1}^{1\text{2e}}\right ) \end{array} \right ] = \mathbf{0}_{16}\times \mbox{1} \end{equation}

Figure 4.

Schematic diagram of multibody kinematics of the ankle exoskeleton.

The external constraints of the exoskeleton system are represented by ${\boldsymbol{C}}_{3}$ :

(6) \begin{equation} {\boldsymbol{C}}_{3}= \left [ \begin{array}{c} {\boldsymbol{R}}^{\text{1e}}+{\boldsymbol{A}}^{\text{ 1e}}\overline {{\boldsymbol{u}}}_{\scriptscriptstyle A}^{\text{1e}}\\[5pt] r_{2}^{1\text{2e}}-\left (L^{1\text{2e}}/2\right )\sin \theta ^{1\text{2e}} \end{array} \right ]=\mathbf{0}_{3\times 1} \end{equation}

The system configuration variables can be described as follows:

(7) \begin{equation} {\boldsymbol{q}}_{n}=\left [{\boldsymbol{q}}_{i}^{\mathrm{T}}\quad {\boldsymbol{q}}_{d}^{\mathrm{T}}\right ]_{1\times 45}^{\mathrm{T}} \end{equation}

The independent coordinates in the constraint equation set can be described as:

(8) \begin{equation} {\boldsymbol{q}}_{i}=\left [\theta ^{\text{1e}}\quad {\boldsymbol{q}}^{6\mathrm{T}}\right ]_{1\times 4}^{\mathrm{T}} \end{equation}

and the dependent coordinates can be described as:

(9) \begin{equation} {\boldsymbol{q}}_{d}=\left [{\boldsymbol{R}}^{\text{1e}\mathrm{T}}\quad {\boldsymbol{q}}^{i\mathrm{T}}\quad \theta ^{13re}\quad L^{4\text{5e}}\quad L^{7\text{8e}}\right ]_{1\times 41}^{\mathrm{T}} \end{equation}

where ${\boldsymbol{q}}_{d}$ denotes the responses of all moving components in the exoskeleton. By taking the total differential of Eq. (3) [Reference Yuwei, Bin, Haibo, Xinhua, Yang, Liang and Lei36], we obtain the following:

(10) \begin{equation} {\boldsymbol{C}}_{{\boldsymbol{q}}_{i}}\dot {{\boldsymbol{q}}}_{i}+{\boldsymbol{C}}_{{\boldsymbol{q}}_{d}}\dot {{\boldsymbol{q}}}_{d}=\mathbf{0}_{41\times 1} \end{equation}

In the equation, ${\boldsymbol{C}}_{{\boldsymbol{q}}_{i}}$ and ${\boldsymbol{C}}_{{\boldsymbol{q}}_{d}}$ represent the partial derivatives of the constraint equations with respect to ${\boldsymbol{q}}_{i}$ and ${\boldsymbol{q}}_{d}$ , respectively, where ${\boldsymbol{C}}_{{\boldsymbol{q}}_{d}}$ constitutes the Jacobian matrix of the system’s motion transmission. They are expressed as follows:

(11) \begin{equation} \left \{ \begin{aligned} {\boldsymbol{C}}_{{\boldsymbol{q}}_{i}}&=\frac {\partial {\boldsymbol{C}}}{\partial {\boldsymbol{q}}_{i}}= \left [\frac {\partial {\boldsymbol{C}}_{1}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{i}}\quad \frac {\partial {\boldsymbol{C}}_{2}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{i}}\quad \frac {\partial {\boldsymbol{C}}_{3}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{i}}\right ]_{3\times 41}^{\mathrm{T}}\\[5pt] {\boldsymbol{C}}_{{\boldsymbol{q}}_{d}}&=\frac {\partial {\boldsymbol{C}}}{\partial {\boldsymbol{q}}_{d}}= \left [\frac {\partial {\boldsymbol{C}}_{1}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{d}}\quad \frac {\partial {\boldsymbol{C}}_{2}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{d}}\quad \frac {\partial {\boldsymbol{C}}_{3}^{\mathrm{T}}}{\partial {\boldsymbol{q}}_{d}}\right ]_{41\times 41}^{\mathrm{T}} \end{aligned} \right . \end{equation}

Thus, the mapping relationship between $\dot {{\boldsymbol{q}}}_{i}$ and $\dot {{\boldsymbol{q}}}_{d}$ can be expressed as:

(12) \begin{equation} \boldsymbol{C}_{di} = -\boldsymbol{C}_{q_d}^{-1}\boldsymbol{C}_{q_i}, \qquad \dot {\boldsymbol{q}}_{d} = \boldsymbol{C}_{di}\dot {\boldsymbol{q}}_{i} \end{equation}

where ${\boldsymbol{C}}_{di}$ is the velocity Jacobian matrix of the system, obtained by combining $\dot {{\boldsymbol{q}}}_{i}$ and $\dot {{\boldsymbol{q}}}_{d}$ , and can be expressed as follows:

(13) \begin{equation} \dot {\boldsymbol{q}}_{n} = \frac {\mathbf{d}\boldsymbol{q}_{n}}{\mathbf{d}t} = \boldsymbol{B}_{di}\dot {\boldsymbol{q}}_{i} \end{equation}

where ${\boldsymbol{q}}_{n}$ is given in Eq. (7), and the expression for ${\boldsymbol{B}}_{di}$ is as follows:

(14) \begin{equation} {\boldsymbol{B}}_{di}=\left [{\boldsymbol{I}}\quad {\boldsymbol{C}}_{di}^{\mathrm{T}}\right ]_{4\times 45}^{\mathrm{T}} \end{equation}

where ${\boldsymbol{I}}$ denotes the identity matrix.

3.2. Forward kinematics

In the exoskeleton design, link 6 serves as the wearable component for human–exoskeleton interaction and is fixed to the human shank through straps. Therefore, the swing angle of the shank and that of the corresponding exoskeleton component are kinematically consistent. The lower end of this component is connected to the foot-support region, and the motion of the human foot is transmitted to the lower end of the exoskeleton through the support-contact relationship. As a result, the swing angle of the lower end of the exoskeleton remains synchronized with the foot motion, establishing a kinematic equivalent relationship between human input and exoskeleton response. In this model, three augmented configuration variables are introduced, namely, $L^{\text{45e}}$ , $L^{\text{78e}}$ , and $\theta ^{\text{13re}}$ . Here, $L^{\text{45e}}$ represents the equivalent displacement response of an ankle soft tissue-mimetic elastic load-bearing structure; $L^{\text{78e}}$ represents the equivalent displacement response of a soleus muscle-inspired elastic-damping energy-storage assistance mechanism; and $\theta ^{\text{13re}}$ represents the relative angular response between the roller and the contact guide rail in an ankle-joint motion-inspired passive compliant support.

The derivation begins with the length of $L^{\text{45e}}$ , as shown in Figure 4(c). By introducing angle $\alpha$ and angle $\beta$ , their relationship can be expressed as:

(15) \begin{equation} \tan \alpha =\frac {r_{2}^{1\text{3e}}-r_{2}^{Fe}}{r_{1}^{1\text{3e}}-r_{1}^{Fe}} \end{equation}

There is a geometric relationship between angles $\alpha$ and $\beta$ , which can be expressed as:

(16) \begin{equation} \beta =\pi -(\pi -\alpha )-\theta ^{\text{5e}}=\alpha -\theta ^{\text{5e}} \end{equation}

Here, the link $L^{EFe}$ is formed by three segments, $L^{\text{4e}}$ , $L^{\text{5e}}$ , and $L^{\text{45e}}$ ; their relationship is given as follows:

(17) \begin{equation} L^{EFe}=L^{\text{4e}}+L^{4\text{5e}}+L^{\text{5e}} \end{equation}

Based on the law of cosines, $L^{\text{45e}}$ is derived as follows:

(18) \begin{equation} \begin{aligned} L^{4\text{5e}}={}& \frac {\cos \theta ^{\text{5e}}+\tan \alpha \sin \theta ^{\text{5e}}}{\sqrt {1+\tan ^{2}\alpha }}L^{FBe}\\[5pt] &\pm \sqrt {(L^{1\text{3e}}-L^{1\text{4e}})^{2} -\left (L^{FBe}\frac {\tan \alpha \cos \theta ^{\text{5e}}-\sin \theta ^{\text{ 5e}}}{\sqrt {1+\tan ^{2}\alpha }}\right )^{2}} -L^{\text{4e}}-L^{\text{5e}} \end{aligned} \end{equation}

From the geometric relationship between ${\boldsymbol{R}}^{Ge}$ and ${\boldsymbol{R}}^{He}$ , as shown in Figure 4(d), the expression for $L^{\text{78e}}$ is

(19) \begin{equation} L^{7\text{8e}}=\sqrt {\left (r_{1}^{Ge}-r_{1}^{He}\right )^{2}+\left (r_{2}^{Ge}-r_{2}^{He}\right )^{2}} \end{equation}

From the derivation of $\theta ^{\text{13re}}$ shown in Figure 4(e), because $\overline {{\boldsymbol{u}}}^{13re}=\left [L^{1\text{3e}}\cos \theta ^{13re}\quad L^{1\text{ 3e}}\sin \theta ^{13re}\right ]^{\mathrm{T}}$ , a constraint can be established:

(20) \begin{equation} {\boldsymbol{R}}^{\text{4e}}+{\boldsymbol{A}}^{\text{ 4e}}\overline {{\boldsymbol{u}}}_{E}^{\text{ 4e}}-{\boldsymbol{R}}^{1\text{3e}}-\overline {{\boldsymbol{u}}}^{13re}=\mathbf{0}_{2\times 1} \end{equation}

By rearranging Eq. (20), we obtain $\tan \theta ^{13re}$ :

(21) \begin{equation} \tan \theta ^{13re}= \frac {\sin \theta ^{13re}}{\cos \theta ^{13re}}= \frac {r_{2}^{\text{4e}}-\frac {L^{\text{4e}}}{2}\sin \theta ^{\text{4e}}-r_{2}^{1\text{3e}}} {r_{1}^{\text{4e}}-\frac {L^{\text{4e}}}{2}\cos \theta ^{\text{4e}}-r_{1}^{1\text{3e}}}. \end{equation}

By simplifying Eq. (21), we obtain $\theta ^{\text{13re}}$ :

(22) \begin{equation} \theta ^{13re}= \arctan \left ( \frac {r_{2}^{\text{4e}}-\frac {L^{\text{4e}}}{2}\sin \theta ^{\text{4e}}-r_{2}^{1\text{3e}}} {r_{1}^{\text{4e}}-\frac {L^{\text{4e}}}{2}\cos \theta ^{\text{4e}}-r_{1}^{1\text{3e}}} \right ) \end{equation}

In this study, $r_{1}^{Ie}$ is also an unknown that needs to be solved by forward kinematics, as shown in Figure 4(f). Using the tangent theorem, we obtain the following:

(23) \begin{equation} \tan \theta ^{1\text{2e}} = \frac {r_{2}^{Be}-r_{2}^{\,Je}}{r_{1}^{Be}-r_{1}^{\,Je}} \end{equation}

Meanwhile, $r_{2}^{\,Je}=0$ ; thus, it can be inferred that:

(24) \begin{equation} r_{2}^{\,Je}+L^{1\text{2e}}\sin \theta ^{1\text{2e}}=r_{2}^{Be} \end{equation}

where $\sin \theta ^{1\text{2e}}$ is taken under $\theta ^{1\text{2e}}\lt \pi /2$ .

Given $\theta ^{1\text{1e}}=\theta ^{1\text{2e}}$ , the coordinate relationship between the points $I$ and $C$ can be established. Consequently, $L^{\text{10e}}$ can be expressed as follows:

(25) \begin{equation} L^{\text{10e}}=\sqrt {\left (r_{2}^{\,Ie}-r_{2}^{Ce}\right )^{2}+\left (r_{1}^{Ie}-r_{1}^{Ce}\right )^{2}} \end{equation}

Let $k=\tan \theta ^{1\text{2e}}$ ; then the relationship can be expressed as:

(26) \begin{equation} r_{2}^{\,Je}=kr_{1}^{\,Je}+b \end{equation}

The relationship is as follows:

(27) \begin{equation} r_{2}^{\,Ie}=kr_{1}^{Ie}+b \end{equation}

by substituting Eq. (26) into Eq. (27), we obtain the following:

(28) \begin{equation} r_{2}^{\,Ie}=k\left (r_{1}^{Ie}-r_{1}^{\,Je}\right )+r_{2}^{\,Je} \end{equation}

by substituting Eq. (28) into Eq. (25), we obtain the following:

(29) \begin{equation} \left (r_{1}^{Ie}-r_{1}^{Ce}\right )^{2} +\left (k\left (r_{1}^{Ie}-r_{1}^{\,Je}\right )+r_{2}^{\,Je}-r_{2}^{Ce}\right )^{2} =\left (L^{\text{10e}}\right )^{2} \end{equation}

by rearranging the above equation, we obtain:

(30) \begin{equation} r_{1}^{Ie}= \frac {r_{1}^{Ce}-\tan \theta ^{1\text{2e}}\left (b-r_{2}^{Ce}\right )}{\left (\tan \theta ^{1\text{2e}}\right )^{2}+1} \pm \frac {\sqrt {\left (\left (\tan \theta ^{1\text{2e}}\right )^{2}+1\right )\left (L^{\text{10e}}\right )^{2} -\left (\tan \theta ^{1\text{2e}}r_{1}^{Ce}+b-r_{2}^{Ce}\right )^{2}}} {\left (\tan \theta ^{1\text{2e}}\right )^{2}+1} \end{equation}

by rearranging the equation, $r_{1}^{Ie}$ can be obtained under $\theta ^{\text{10e}}\lt \pi /2$ .

4. Chebyshev polynomial response-interval method

Based on the kinematic mapping relationship between human input and exoskeleton motion responses established in the previous section, the Chebyshev polynomial response-interval method [Reference Wu, Zhang, Chen and Luo37] is introduced to quantify the influence of wearing-position uncertainty on key motion response variables. The considered uncertainty mainly arises from small variations in human input caused by individual differences or gait fluctuations, which are transmitted through the kinematic mapping to the displacement, velocity, and acceleration responses of the exoskeleton. By constructing the mapping relationship between uncertain inputs and response intervals, the sensitivity of the exoskeleton motion response under given structural parameters can be evaluated. A narrower response interval indicates lower kinematic sensitivity of the current configuration parameters to external disturbances. In the following section, uncertainty response-interval analysis is conducted for key response variables such as $L^{\text{45e}}$ , $L^{\text{78e}}$ , and $\theta ^{\text{13re}}$ , and the method is briefly introduced using a one-dimensional Chebyshev polynomial as an example.

For a one-dimensional uncertain variable, the $n$ -th-order Chebyshev polynomial can be expressed as follows:

(31) \begin{equation} \left \{ \begin{array}{@{}l@{}} C_{n}(\xi )=\cos (n\theta )\\[5pt] \theta =\arccos \left ( \dfrac { \mbox{2}\xi -\left (\chi _{\mathrm{up}}+\chi _{\mathrm{low}}\right ) }{ \chi _{\mathrm{up}}-\chi _{\mathrm{low}} } \right ) \end{array} \right . \end{equation}

In the above equation, $\xi \in [\chi _{low},\chi _{up}]$ , $\chi _{up}$ and $\chi _{low}$ represent the upper and lower bounds of the uncertain parameter, and $\theta \in [0,\pi ]$ . Then, the continuous function $f(\xi )$ can be approximated as follows:

(32) \begin{equation} f(\xi )\approx \frac {1}{2}f_{0}+\sum _{i=1}^{n}f_{i}C_{i}(\xi ) \end{equation}

where $f_{i}$ denotes the constant coefficients of the polynomial, and $n$ represents the highest order of the Chebyshev polynomial. Based on the orthogonality of Chebyshev polynomials, the coefficients can be derived using Mehler’s integral method:

(33) \begin{equation} f_{i}=\frac {2}{\pi }\int _{0}^{\pi }f(\xi )\cos \!(i\theta )\,d\theta \approx \frac {2}{m}\sum _{j=1}^{m}f(\xi ^{j})\cos (i\theta ^{j}) \end{equation}

where

(34) \begin{equation} \left \{ \begin{aligned} \theta ^{j}&=\frac {2j-1}{2m}\pi \\[5pt] \xi ^{j}&=\frac {\chi _{up}+\chi _{low}}{2}+\frac {\chi _{up}-\chi _{low}}{2}\cos \theta ^{j} \end{aligned} \right . \end{equation}

Here, $j$ = 1,2, … , $m$ , with $m$ representing the number of interpolation points. To ensure accuracy in approximation, it is generally required that $m\geq n+1$ . For multidimensional problems, Chebyshev polynomials can be seen as the tensor product of equations, yielding $k$ -dimensional Chebyshev polynomials. In this design, where the dimensionality $k$ is 2, the motion response can be expressed as:

(35) \begin{equation} f(\xi _{c},\theta ^{6})\approx \sum _{i_{1}=0}^{n}\sum _{i_{2}=0}^{n} \left (\frac {1}{2}\right )^{s} f_{i_{1},i_{2}}(\theta ^{6})C_{i_{1},i_{2}}(\xi _{c}) \end{equation}

Wearing-position uncertainty is equivalently represented as an initial positional disturbance at the human–exoskeleton connection point, including horizontal and vertical deviations of the wearing position in the sagittal plane. Accordingly, the uncertain input vector is defined as:

(36) \begin{equation} \xi _{c} = \Delta \boldsymbol{R}^{{\scriptscriptstyle 6}\mathrm{e}} = \left [ r_{1}^{{\scriptscriptstyle 6}\mathrm{e}}\pm d_x \quad r_{2}^{{\scriptscriptstyle 6}\mathrm{e}}\pm d_y \right ]_{1\times 2}^{\mathrm{T}} \end{equation}

where $d_x$ and $d_y$ denote the horizontal and vertical deviations of the wearing position, respectively. The motion response interval is defined as:

(37) \begin{equation} \left \{ \begin{aligned} {\boldsymbol{y}}_{low}&=\left (\frac {1}{2}\right )^{2}{\boldsymbol{y}}_{0,0}(\theta ^{6}) -\sum _{i_{1}=0}^{n}\sum _{i_{k}=0}^{n} \left |\left (\frac {1}{2}\right )^{s}{\boldsymbol{y}}_{i_{1},\ldots ,i_{k}}(\theta ^{6})\right |\\[5pt] {\boldsymbol{y}}_{up}&=\left (\frac {1}{2}\right )^{2}{\boldsymbol{y}}_{0,0}(\theta ^{6}) +\sum _{i_{1}=0}^{n}\sum _{i_{k}=0}^{n} \left |\left (\frac {1}{2}\right )^{s}{\boldsymbol{y}}_{i_{1},\ldots ,i_{k}}(\theta ^{6})\right | \end{aligned} \right . \end{equation}

In the equation, ${\boldsymbol{y}}_{up}$ and ${\boldsymbol{y}}_{low}$ denote the upper and lower bounds of the system response, respectively. By introducing Chebyshev polynomials, the mapping relationship between the uncertainty factors of the exoskeleton system and the system response is established.

5. Biomimetic kinematic response containment-smoothness

When human limb motion is used as input, the exoskeleton system generates the corresponding motion responses. If the motion response of the exoskeleton mechanism is not coordinated with the natural movement tendency of the human body, additional geometric constraints may be introduced at the wearing interfaces, thereby affecting human–exoskeleton motion coordination. Conventional evaluation methods usually assess motion smoothness using the range, mean value, or standard deviation of the response of a single component. However, because different response quantities have different physical dimensions, direct comprehensive comparison among them is difficult. Therefore, based on the multibody kinematic constraint model and the Chebyshev response-interval method, this study proposes a comprehensive kinematic evaluation index, termed the BKRCS index. In this study, response containment refers to the ability of the mechanism to bound and attenuate uncertainty-induced variations in key kinematic responses. The proposed BKRCS characterizes both uncertainty-induced response containment and motion-transmission smoothness of the exoskeleton mechanism.

According to Eq. (37), the upper and lower bounds of the key motion responses of the exoskeleton under wearing-position uncertainty are obtained. For a passive ankle exoskeleton, within a prescribed range of wearing-position uncertainty, a larger response interval associated with a structural parameter indicates higher sensitivity to wearing-position deviations. In contrast, a smaller response interval implies that the structure can attenuate the amplification of uncertain inputs in the motion response, thereby reflecting better kinematic containment. Accordingly, the interval width of the key motion responses is used as the basis for the containment evaluation, and the response range is defined as:

(38) \begin{equation} W_{s}(k) = \bigl (\,\boldsymbol{y}_{\mathrm{up}}(k)\bigr )_{s} - \left (\,\boldsymbol{y}_{\mathrm{low}}(k)\right )_{s} \end{equation}

where $s$ indexes the key motion response and $k$ represents the discrete point in the gait cycle. Due to differences in physical dimensions and numerical scales among different key response quantities, their original response-interval widths cannot be directly used for weighted averaging. Therefore, $S_{s}$ is adopted as the reference scale:

(39) \begin{equation} S_{s} = \max \, \bigl (\,\boldsymbol{y}_{\mathrm{up}}(k)\bigr )_{s} - \min \left (\,\boldsymbol{y}_{\mathrm{low}}(k)\right )_{s} \end{equation}

To ensure consistency among response quantities with different physical dimensions, their response-interval widths are uniformly non-dimensionalized as follows:

(40) \begin{equation} \varphi _{s}(k)=\frac {W_{s}(k)}{S_{s}} \end{equation}

An arctangent mapping function is further introduced to project the non-dimensionalized response elements, originally associated with different physical dimensions, into a unified bounded interval. The mapping function is expressed as:

(41) \begin{equation} y=\frac {2}{\pi }\arctan \!(x) \end{equation}

where $x$ is the non-dimensionalized value and $y$ is the scaled value. After the response quantities are non-dimensionalized, the arctangent mapping function is applied to project them into a unified bounded interval, thereby enabling dimensionally consistent comparison among different responses. The mean value over all discrete gait points is then computed as:

(42) \begin{equation} B_{Ck}=\frac {1}{n}\sum _{k=1}^{n}\left |\frac {2}{\pi }\arctan \left (\varphi _{s}(k)\right )\right | \end{equation}

Here, $n$ denotes the number of discrete points in the normalized prescribed stance-related stepping condition. Considering m key motion responses, the containment term of the BKRCS associated with wearing-position uncertainty is defined as:

(43) \begin{equation} B_{C}=\sum _{s=1}^{m}\omega _{s}B_{Ck} \end{equation}

in this expression, $B_{C}$ represents the response-containment term in the BKRCS, and $\omega _{s}$ denotes the weighting coefficient of the $s$ -th response. A smaller $B_{C}$ indicates that, under the same wearing-position uncertainty, the interval width of the key motion responses becomes narrower and the mechanism is less sensitive to wearing-position uncertainty, which corresponds to better kinematic compatibility.

The partial derivative of the system Jacobian matrix ${\boldsymbol{C}}_{di}$ with respect to the configuration variable ${\boldsymbol{q}}_{n}$ is given by:

(44) \begin{equation} \left ({\boldsymbol{C}}_{di}^{m}\right )_{{\boldsymbol{q}}_{n}}\delta {\boldsymbol{q}}_{n} =\frac {\partial {\boldsymbol{C}}_{di}^{m}}{\partial {\boldsymbol{q}}_{n}}\delta {\boldsymbol{q}}_{n} =\frac {\partial {\boldsymbol{C}}_{di}^{m}}{\partial {\boldsymbol{q}}_{n}}{\boldsymbol{B}}_{di}\delta {\boldsymbol{q}}_{i} ={\boldsymbol{C}}_{Cdi}^{m}\delta {\boldsymbol{q}}_{i} \end{equation}

where $m=1,\ldots ,r$ , ${\boldsymbol{C}}_{di}^{m}$ represents the $m$ -th row element of ${\boldsymbol{C}}_{di}$ , and $r$ denotes the number of elements in ${\boldsymbol{q}}_{di}$ . The sensitivity matrix of ${\boldsymbol{C}}_{di}^{m}$ with respect to the configuration variables is denoted as ${\boldsymbol{C}}_{Cdi}^{m}$ . The larger the absolute value of an element in this matrix, the more significant the transmission and scaling effect of the corresponding velocity component in $\dot {{\boldsymbol{q}}}_{i}$ on the elements of the Jacobian matrix during motion. At a given instantaneous configuration, the numerical values of the elements in ${\boldsymbol{C}}_{Cdi}^{m}$ and ${\boldsymbol{B}}_{di}$ remain constant. The values in ${\boldsymbol{C}}_{Cdi}^{m}$ reflect the intensity of changes in the elements of the Jacobian matrix under the current configuration. However, even for the same system configuration, the velocities of the independent coordinates may vary. To prevent the velocity itself from influencing the relationship between scaling parameters and motion smoothness, ${\boldsymbol{C}}_{Cdi}^{m}$ can be used to evaluate whether the motion response undergoes abrupt changes under the current configuration. This matrix also suffers from dimensional inconsistency among its elements. The evaluation criterion for this matrix is that its elements should be as close to zero as possible. Therefore, to eliminate the effect of dimensional inconsistency, data with different physical dimensions are first non-dimensionalized and then mapped into a unified bounded interval via the arctangent function. On this basis, the local configuration-response smoothness of each system can be expressed as:

(45) \begin{equation} B_{KCl}=\frac {1}{lp}\sum _{m=1}^{l}\sum _{j=1}^{p}\left |\frac {2}{\pi }\arctan \left (C_{i}^{mj}\right )\right | \end{equation}

in the equation, let $B_{KCl}\in [0,1]$ , let $C_{i}^{mj}$ represent the average value of the $j$ -th element in the dimensionless matrix, and let $p$ denote the total number of elements. Then, the BKRCS of the system under a global instantaneous configuration can be expressed as follows:

(46) \begin{equation} B_{KC}=\frac {1}{n}\sum _{k=1}^{n}\left (B_{KC}\right )_{k} \end{equation}

in the equation, $B_{KC}\in [0,1]$ , and $\left (B_{KC}\right )_{k}$ represents the local configuration-response smoothness of the system under the $k$ -th configuration.

Assuming the system has n configurations in space, the scaled average variation distance of all elements in the Jacobian matrix can be expressed as:

(47) \begin{align} B_{KJ} = \frac {1}{n} \sum _{j=1}^{r} \frac {2}{\pi } \arctan \Bigg( \Bigg | \left (C_{di}^{j}\right )_{k} - \frac { \displaystyle \sum _{k=1}^{n}\left (C_{di}^{j}\right )_{k} }{n} \Bigg | \Bigg ) \end{align}

Here, $B_{KJ}$ serves as an intuitive indicator of the system’s motion smoothness. From the above equation, it follows that $B_{KJ}\in [0,1]$ . A smaller value of $B_{KJ}$ indicates smaller variations of the elements in the system’s Jacobian matrix across different configurations. Here, $\left ({\boldsymbol{C}}_{di}^{j}\right )_{k}$ denotes the $j$ -th element in the Jacobian matrix under the $k$ -th configuration, and $l$ represents the total number of elements in the Jacobian matrix.

For concise notation, the BKRCS is denoted by $B_{K}$ . The BKRCS of the exoskeleton can therefore be formulated as:

(48) \begin{equation} B_{K}=\mu _{\scriptscriptstyle BKJ}\cdot B_{KJ}+\mu _{\scriptscriptstyle BKC}\cdot B_{KC}+\mu _{\scriptscriptstyle BC}\cdot B_{C} \end{equation}

where $\mu _{\scriptscriptstyle BKJ}$ , $\mu _{\scriptscriptstyle BKC}$ , and $\mu _{\scriptscriptstyle BC}$ denote the weighting coefficients of the BKRCS, which balance the relative contributions of the local configuration-response smoothness $B_{\scriptscriptstyle KC}$ , the global configuration variation smoothness $B_{KJ}$ , and the biomimetic response-containment term $B_{C}$ . Specifically, $B_{\scriptscriptstyle KC}$ describes the instantaneous variation of the motion-transmission Jacobian matrix around the current configuration and reflects local abrupt changes in the mechanism response. $B_{KJ}$ represents the average variation distance of the Jacobian matrix among different configurations over the stance phase, thereby evaluating the overall continuity of motion transmission during configuration evolution. $B_{C}$ quantifies the containment capability of the mechanism under wearing-position uncertainty by characterizing the response variation within the uncertainty interval.

The weighted combination of these three terms enables BKRCS to evaluate both local configuration-induced response variations and global motion-transmission smoothness. When $\mu _{\scriptscriptstyle BKJ}$ , $\mu _{\scriptscriptstyle BKC}$ , $\mu _{\scriptscriptstyle BC}\gt 0$ , and $\mu _{\scriptscriptstyle BKJ}+\mu _{BKC}+\mu _{\scriptscriptstyle BC}=1$ , a smaller BKRCS value indicates that the motion-transmission Jacobian matrix of the exoskeleton fluctuates less with configuration changes during the stance phase. A lower BKRCS indicates a lower weighted combination of response-interval width and Jacobian-variation terms within the prescribed model.

6. Multi-objective kinematic dimensional synthesis and simulation analysis of the exoskeleton

Based on the sagittal-plane human–exoskeleton–environment equivalent kinematic model, the Chebyshev response-interval model, and the BKRCS, multibody kinematic dimensional synthesis is further performed for the key structural parameters of the exoskeleton. By adjusting the structural dimensions, the synthesis aims to improve the smoothness of the mechanism response under a prescribed human input and attenuate the influence of wearing-position uncertainty on the key motion responses. A hybrid optimization strategy combining particle swarm optimization (PSO) [Reference Wang, Tan and Liu38] with a genetic algorithm (GA) [Reference Schoenauer and Xanthakis39] is adopted for multi-objective kinematic dimensional synthesis. In this strategy, GA enhances population diversity and reduces the risk of local convergence, whereas PSO improves the search efficiency in the parameter space. The hybrid strategy therefore provides both global exploration and local search capabilities, making it suitable for the multivariable and nonlinear dimensional synthesis problem addressed in this study. The optimized parameter combination is obtained by balancing motion-response smoothness and structural dimensional constraints. The kinematic performance of the exoskeleton is evaluated using $L^{\text{45e}}$ , $L^{\text{78e}}$ , and $\theta ^{\text{13re}}$ as the core responses. The multibody kinematic dimensional synthesis is described below.

The parameters for optimization in ankle exoskeleton design are as follows:

(49) \begin{equation} \left [L^{ie}\quad \phi ^{je}\right ]_{1\times 18}^{\mathrm{T}}, i=1\sim 10,12,13,14, j=1\sim 5 \end{equation}

In the equation, $L^{ie}$ represents the length of the $i$ -th rod member, and $\phi ^{je}$ denotes the angle between the $j$ -th fixed pair. The expression for the total structural parameter length $L_{s}$ is as follows:

(50) \begin{equation} L_{s}=\sum _{i=1}^{13}L^{ie},i=1\sim 10,12\sim 14. \end{equation}

To remove the dimensional influence of $L_{s}$ on the objective function, the total link length is normalized as:

(51) \begin{equation} \overline {L}_s=\frac {L_s-L_{s,\min }}{L_{s,\max }-L_{s,\min }},\qquad \overline {L}_s\in [0,1]. \end{equation}

In this study, the total structural length is used only as an indirect indicator of structural compactness and lightweight potential under the assumption of identical material properties.

This study performs the dimensional synthesis of the ankle exoskeleton system, and the objective function can be expressed as follows:

(52) \begin{equation} f=\mu _{\scriptscriptstyle \overline {L}_s}\overline {L}_s+\mu _{\scriptscriptstyle BK}B_{K}, \end{equation}

where $\mu _{\scriptscriptstyle BK}$ and $\mu _{\scriptscriptstyle \overline {L}_s}$ are the optimization weighting coefficients. The BKRCS characterizes the uncertainty-induced response containment and smoothness of the exoskeleton mechanism under a prescribed human input, whereas the total structural length is used as an indirect measure of structural compactness and lightweight potential. Since the main objective of the dimensional synthesis is to improve the kinematic smoothness of the passive ankle exoskeleton under wearing-position uncertainty, the BKRCS is treated as the primary optimization term. The total structural length is included as an auxiliary constraint term to limit excessive structural dimensions and promote lightweight design.

The constraint intervals of the structural parameters in Table I are determined by considering the anthropometric dimensions of the ankle–foot complex, the available installation space, the mechanism motion range, component-interference avoidance, stroke constraints, and manufacturing and assembly feasibility. The lower bounds of the length parameters are imposed to prevent joint interference or singular configurations caused by excessively short links, whereas the upper bounds restrict the overall structural envelope and improve wearable integration. The angular parameter ranges are specified to ensure reasonable installation postures of the fixed joints and to avoid infeasible configurations during stance-phase motion. The kinematic dimensional synthesis procedure is illustrated in Figure 5.

Table I.

Constraint ranges of the design parameters.

Table I. Long description

The table presents constraint intervals for various structural parameters of an ankle-foot complex. It has two main sections, each with two columns labeled Parameters and Interval. The first section lists parameters L 1e through L 1ce with intervals ranging from millimeters to degrees. The second section continues with parameters L 12e through dy with similar interval ranges. Each row specifies a parameter and its corresponding interval, detailing the constraints for length and angular parameters to ensure proper functioning and integration of the ankle-foot mechanism.

During the stepping process, relative slip occurs at the shank interface of the ankle exoskeleton. To characterize the relative motion between the human shank and the exoskeleton interface, the slip displacement is specified as 15 mm in the horizontal or anterior–posterior direction and 25 mm in the vertical direction, based on existing studies on human–exoskeleton interface relative motion in ankle–foot exoskeletons [Reference Akiyama, Okamoto, Yamada and Ishiguro40]. The results shown in Tables I–III indicate that, under the prescribed parameter bounds and multi-objective kinematic formulation, the objective-function value decreased from 0.2552 to 0.2276, corresponding to a reduction of approximately 10.8%. The BKRCS also decreased from 0.2055 to 0.1902, representing a reduction of approximately 7.45%. Among its component terms, $B_{\mathrm{\scriptscriptstyle \textit {KJ}}}$ decreased from 0.2509 to 0.2106, representing a reduction of approximately 16.06%, whereas $B_{\mathrm{\scriptscriptstyle \textit {KC}}}$ and $B_{\mathrm{\scriptscriptstyle \textit {C}}}$ decreased only slightly, from 0.1302 to 0.1293 and from 0.2356 to 0.2307, respectively. Therefore, the reduction in BKRCS is mainly attributable to the decreased global variation of the motion-transmission Jacobian, while the local configuration-response smoothness and response-interval containment show comparatively small changes after optimization. These results indicate a component-dependent improvement in the selected comprehensive kinematic metric rather than a uniform improvement in all three BKRCS components.

Figure 5.

Procedure for kinematic dimensional synthesis (CIM, Chebyshev polynomial response-interval method; BKRCS, biomimetic kinematic response containment–smoothness index; PSO, particle swarm optimization; GA, genetic algorithm).

To further examine the relationship between the reduction in BKRCS and conventional trajectory-smoothness measures, LDLJ [Reference Balasubramanian, Melendez-Calderon, Roby-Brami and Burdet41] and RMS jerk [Reference van Meulen, Klaassen, Held, Reenalda, Buurke, van Beijnum, Luft and Veltink42] were introduced as auxiliary comparison metrics. It should be noted that these metrics are used only to support the interpretation of the optimization results and do not constitute an independent validation of BKRCS. As shown in Table IV, the RMS jerk of $L^{\text{45e}}$ decreases from 2.142 to 0.946 after optimization, indicating a marked reduction in jerk magnitude. For $L^{\text{78e}}$ , LDLJ changes from −4.854 to −3.845, moving closer to zero, whereas RMS jerk increases only slightly from 2.047 to 2.054. Therefore, the smoothness-related change in $L^{\text{78e}}$ is relatively limited, and the two metrics exhibit slightly different trends. For $\theta ^{\text{13re}}$ , LDLJ changes from −7.388 to −3.819, while RMS jerk decreases from 255.679 to 17.861, indicating that both metrics reflect a marked improvement in the numerical kinematic smoothness of this response. Overall, the conventional smoothness metrics exhibit trends that are broadly consistent with the reduction in BKRCS, providing auxiliary support for the rationality and feasibility of using BKRCS as a numerical design metric for kinematic dimensional synthesis.

Figure 6 compares the displacement, velocity, and acceleration responses predicted by the Chebyshev interval method (CIM) with those obtained from direct kinematic calculations. The two sets of curves exhibit generally consistent temporal trends, indicating that the CIM captures the principal variation characteristics of the selected kinematic responses.

For $L^{\text{45e}}$ , kinematic fluctuations are consistently and substantially attenuated at all three levels. The displacement range narrows from 0.021 to 0.002 m, with its standard deviation falling from 0.007 to 0.001 m, corresponding to reductions of 86.12% and 85.40%, respectively. The velocity range decreases from 0.124 to 0.053 m/s, alongside a decline in standard deviation from 0.047 to 0.011 m/s. The acceleration range also contracts from 0.765 to 0.416 m/s $^{2}$ , confirming an overall suppression of numerical kinematic variations for this response. For $L^{\text{78e}}$ , a contrasting trend emerges. Its displacement range expands from 0.057 to 0.084 m, with the corresponding standard deviation increasing from 0.01888 to 0.02722 m, and the velocity range likewise rises from 0.361 to 0.601 m/s. Conversely, the acceleration range decreases from 0.819 to 0.713 m/s $^{2}$ , and its standard deviation falls from 0.265 to 0.222 m/s $^{2}$ . Thus, optimization does not simultaneously reduce all response metrics for $L^{\text{78e}}$ but instead yields a trade-off between amplified displacement/velocity variations and moderately mitigated acceleration fluctuations. For $\theta ^{\text{13re}}$ , the optimization effect is most pronounced at the velocity and acceleration levels. The displacement range narrows from 0.558 to 0.503 rad, with displacement standard deviation decreasing from 0.186 to 0.157 rad. The velocity range falls from 12.525 to 5.248 rad/s. More strikingly, the acceleration range plummets from 106.604 to 7.599 rad/s $^{2}$ , while the corresponding standard deviation drops from 27.923 to 2.169 rad/s $^{2}$ , reflecting a marked reduction in velocity and acceleration fluctuations for the roller-guide relative-motion response.

Table II.

Optimization results of the structural parameters.

Table II. Long description

The table presents a comparison of initial and optimized parameters for an ankle exoskeleton. It has two main sections, each with multiple rows and columns. The first section lists parameters with their initial and optimized values, while the second section lists parameters with only their optimized values. Each row represents a different parameter, and the columns show the initial and optimized values for each parameter. The parameters include various measurements in millimeters and degrees, as well as some dimensionless values. Notable parameters include L 1 to the power of 6 in millimeters, phi 2 in degrees, and B K R C, among others. The table shows the changes in these parameters after optimization, indicating improvements in the kinematic performance of the ankle exoskeleton.

Table III.

PSO-GA optimization settings and weighting coefficients.

Table III. Long description

A table titled PSO and GA parameters and weighting coefficients with 5 rows and 8 columns. The columns are labeled as follows: N subscript p, G, omega, c subscript 1, c subscript 2, T subscript GA, rho subscript GA, and omega subscript theta subscript 1 subscript M subscript C. The row values are as follows: Row 1: N subscript p, 50; G, 100; omega, 0.8; c subscript 1, 1.5; c subscript 2, 1.8; T subscript GA, 20; rho subscript GA, 0.3; omega subscript theta subscript 1 subscript M subscript C, 1/3. Row 2: mu subscript B subscript K subscript C, 1/3; mu subscript B subscript K subscript S, 1/3; mu subscript B subscript C, 1/3; mu subscript T subscript r, 0.2; mu subscript B subscript K, 0.8.

Table IV.

Comparison of LDLJ and RMS jerk for key response quantities before and after optimization.

Table IV. Long description

A table with four rows and five columns comparing initial and optimized LDLJ and RMS jerk values for three responses. The columns are labeled Response, Initial LDLJ, Optimized LDLJ, Initial RMS jerk, and Optimized RMS jerk. The rows are labeled with different responses and their corresponding values. For the first response, initial LDLJ is -7.045, optimized LDLJ is -7.143, initial RMS jerk is 2.142, and optimized RMS jerk is 0.946. For the second response, initial LDLJ is -4.854, optimized LDLJ is -3.845, initial RMS jerk is 2.047, and optimized RMS jerk is 2.054. For the third response, initial LDLJ is -7.388, optimized LDLJ is -3.819, initial RMS jerk is 255.679, and optimized RMS jerk is 17.861.

Overall, the optimization produces response-dependent effects. The numerical variations of $L^{\text{45e}}$ and $\theta ^{\text{13re}}$ are generally reduced, whereas $L^{\text{78e}}$ exhibits increased displacement and velocity variations together with reduced acceleration variation. The results should therefore be interpreted as an improvement in selected kinematic response characteristics under the prescribed sagittal-plane model, weighting scheme, and constraint conditions, rather than as a uniform improvement in all response quantities or as direct evidence of reduced dynamic loading, improved cushioning performance, or enhanced wearing comfort.

Figure 6.

Response intervals of selected displacement, velocity, and acceleration variables during the heel-off to foot-strike condition.

Figure 6. Long description

The image contains nine line graphs arranged in a 3x3 grid, each showing different biomechanical characteristics of human gait. Each graph plots data against the gait cycle percentage on the horizontal axis. Panel A: The first graph shows the displacement L 45 (cm) with values ranging from -0.08 to 0.06. The second graph shows the displacement L 78 (cm) with values ranging from 0.06 to 0.18. The third graph shows the angle θ 36 (rad) with values ranging from -1.3 to -0.7. Panel B: The first graph shows the velocity dL 45 (cm/s) with values ranging from -0.08 to 0.06. The second graph shows the velocity dL 78 (cm/s) with values ranging from -0.3 to 0.3. The third graph shows the angular velocity dθ 36 (rad/s) with values ranging from -8 to 8. Panel C: The first graph shows the acceleration ddL 45 (m/s²) with values ranging from -0.4 to 0.5. The second graph shows the acceleration ddL 78 (m/s²) with values ranging from -0.6 to 0.4. The third graph shows the angular acceleration ddθ 36 (rad/s²) with values ranging from -80 to 80. Each graph compares optimization parameters for CIM and KIM, with initial parameters shown in blue and optimization parameters shown in orange.

Figure 7.

Response intervals of selected displacement variables from heel-off to heel-strike.

Figure 7 compares the response intervals of the key variables before and after optimization under the prescribed wearing-position uncertainty. A narrower interval indicates lower kinematic sensitivity to uncertainty. For $L^{\text{45e}}$ , the mean and maximum interval widths increase from 0.05561 m and 0.05631 to 0.068 m and 0.080 m, respectively. For $L^{\text{78e}}$ , they increase only slightly from 0.053 m and 0.053 to 0.054 m and 0.055 m. In contrast, the mean and maximum interval widths of $\theta ^{\text{13re}}$ decrease from 4.302 and 4.619 rad to 1.417 and 1.876 rad, corresponding to reductions of 67.05% and 59.38%, respectively. Therefore, the optimization mainly improves the response containment of $\theta ^{\text{13re}}$ , while $L^{\text{45e}}$ and $L^{\text{78e}}$ exhibit different interval trends. The result should be interpreted as a reduction in the selected comprehensive BKRCS objective rather than a uniform narrowing of all response intervals.

Complementing the time-domain response intervals in Figures 6, and 8 further characterizes the configuration dependence of uncertainty sensitivity by plotting response intervals against input angle $\theta ^{\text{1e}}$ . In the initial configuration, the $\theta ^{\text{13re}}$ interval expands sharply at large input angles (corresponding to the late push-off phase), identifying this pose as the primary region of uncertainty amplification. Following dimensional synthesis, the most substantial interval contraction occurs precisely within this high-sensitivity regime, confirming that optimization selectively mitigates uncertainty amplification in critical poses. For $L^{\text{78e}}$ , interval expansion concentrates at intermediate input angles (mid-stance), reflecting the pose-level trade-off between increased effective stroke of the soleus-mimetic subsystem and reduced acceleration fluctuations. Across the full angular range, all direct kinematic samples lie within CIM-predicted bounds, verifying the method’s sufficient reliability for comparative analysis in dimensional synthesis.

Figure 8.

Comparison of response intervals with respect to the input angle $\theta ^{\text{1e}}$ before and after optimization.

To examine the influence of weighting coefficients on the optimization outcome, a sensitivity analysis was conducted for both the internal weights of BKRCS and the global objective-function weights. Different weighting schemes were tested, including the equal-weighting case, $B_C$ -dominant case, $B_{KC}$ -dominant case, $B_{KJ}$ -dominant case, kinematic-dominant case, and compactness-dominant case, while satisfying $\mu _{BKJ},\mu _{BKC},\mu _{BC}\gt$ 0 and $\mu _{BKJ}+\mu _{BKC}+\mu _{BC}=$ 1. For all cases, the same design-variable bounds, wearing-position uncertainty range, response variables, and PSO-GA optimization settings were adopted. As shown in Table V, the optimized configurations consistently achieved lower BKRCS and objective-function values than the corresponding initial configurations under different weighting combinations. For example, under the equal-weighting scheme, BKRCS decreased from 0.2055 to 0.1902, and the objective function decreased from 0.2552 to 0.2276. When $B_C$ , $B_{KC}$ , and $B_{KJ}$ were assigned dominant weights, the objective function decreased from 0.2648 to 0.2259, from 0.2311 to 0.1523, and from 0.2697 to 0.1255, respectively. These results indicate that the proposed dimensional-synthesis framework is not effective only for a single prescribed weighting scheme but can obtain improved solutions under different design preferences. Nevertheless, the final optimized dimensions and relative improvement magnitudes vary with the selected weights. A larger $\mu _{BC}$ makes the optimization more sensitive to response-interval containment under wearing-position uncertainty, whereas larger $\mu _{BKC}$ or $\mu _{BKJ}$ places greater emphasis on local or global motion-transmission smoothness. In addition, increasing $\mu _{L_s}$ promotes structural compactness, although it may reduce the improvement margin of BKRCS. Therefore, the selected baseline weighting coefficients should be interpreted as a design compromise among response containment, kinematic smoothness, and structural compactness under the prescribed constraints, rather than as universally optimal weights. Overall, the consistent reduction in the objective function under all tested weighting schemes demonstrates the robustness of the proposed optimization procedure with respect to weighting-coefficient variations.

Table V.

Sensitivity analysis results for different weighting coefficients.

Table V. Long description

A table comparing different weighting schemes and their effects on optimization outcomes. The table has 7 rows and 9 columns. The columns are labeled as Case, State, μ_BK, μ_BC, μ_BKC, μ_BKJ, μ_Ls, BKRCS, Ls, and Objective. The rows are labeled as Equal, B_C-dominant, B_KC-dominant, B_KJ-dominant, Kinematic-dominant, and Compactness-dominant. Each row provides values for the initial and optimized states under different weighting schemes. For example, under the equal-weighting scheme, the initial state has μ_BC of 0.80, μ_BKC of 0.20, μ_BKJ of 0.20, μ_Ls of 0.2055, BKRCS of 0.99500, Ls of 0.2552, and the optimized state has μ_BC of 1/3, μ_BKC of 1/3, μ_BKJ of 1/3, μ_Ls of 0.1902, BKRCS of 0.95368, Ls of 0.2276. Similar data is provided for other weighting schemes.

Figure 9.

(a) Comparison of $B_{\scriptscriptstyle KC}$ before and after optimization. (b) Link-length variation before and after optimization.

According to the updated $B_{KC}$ results under the selected representative stance-related stepping condition, as shown in Figure 9(a), the optimized $B_{KC}$ curve is generally located below the initial curve. The mean $B_{KC}$ value decreases from 0.1205 to 0.1134, corresponding to a reduction of approximately 5.92%. Since a smaller $B_{KC}$ value indicates weaker local configuration-response coupling and better local kinematic smoothness, this downward shift demonstrates that the optimized configuration improves the local motion-transmission characteristics of the exoskeleton under the analyzed condition. Therefore, the $B_{KC}$ result provides further evidence that the proposed optimization strategy can effectively reduce the local kinematic response sensitivity and enhance the configuration-response smoothness of the mechanism.

Finally, the rod-length optimization results shown in Figure 9(b) were comprehensively analyzed. The variation trends of the structural parameters indicate that incorporating BKRCS into the optimization objective effectively guided the design toward a specific set of parameter combinations. This parameter combination not only improved kinematic coordination but also promoted structural compactness by using the total structural length as an indirect constraint on structural size. These results further demonstrate, from a kinematic perspective, the effectiveness and applicability of the PSO-GA hybrid optimization strategy in solving multivariable and nonlinear structural dimensional synthesis problems of this type.

It should be noted that several optimized length and angular parameters reach or approach the upper or lower limits of their prescribed constraint intervals. This indicates that, under the selected weighting scheme and design-variable bounds, the optimizer tends to exploit the boundary of the feasible design space to reduce the objective function. Therefore, the optimized configuration should be interpreted as a boundary-feasible solution under the prescribed kinematic and geometric constraints, rather than as a globally optimal structure in an unconstrained design space. To assess its geometric feasibility, the optimized configuration was checked against the predefined link-length ranges, angular ranges, prismatic-joint stroke limits, assembly-space requirements, component-interference constraints, manufacturing-related dimensional bounds, and wearable-envelope constraints. The check showed that all optimized parameters remained within the prescribed feasible intervals and did not violate the predefined geometric or assembly bounds. Nevertheless, this assessment is limited to the kinematic and geometric levels. It does not validate load-bearing capacity, structural strength, stiffness, fatigue resistance, or wearing comfort. These aspects require further tolerance and interference analysis, mechanical strength assessment, dynamic simulation, bench testing, and human–subject experiments.