2.1. Exoskeleton robot

The developed lower limb exoskeleton is shown in Fig. 1. The exoskeleton robot designed in this paper mainly considers the motion in the sagittal plane. For each leg, it has one active degree of freedom (DoF) for the hip joint and knee joint, respectively, and one passive DoF for the ankle joint. The ankle joint has a high structural complexity, and it has an important influence on human walking stability [Reference Maharaj30–Reference Sato32]. Most of existing exoskeleton robots are designed with passive DoF ankle joint, such as Rewalk robot [Reference Esquenazi, Talaty and Packel4] and Mina exoskeleton robot [Reference Mummolo, Peng, Agarwal, Griffin, Neuhaus and Kim5]. On the one hand, it is to reduce the design difficulty and cost of the exoskeleton robot; on the other hand, the study found that when the ankle joint motion is disturbed by external interference, the walking stability of human body will be greatly reduced [Reference Sato32]. Hence, this study makes full use of the free movable property of the ankle joint to improve the walking stability.

Figure 1.

(a) Overall structure diagram (1) back support; (2) actuators of hip joint; (3) hip joint component; (4) brace of thigh; (5) actuators of knee joint; (6) brace of calf; (7) ankle joint component; (8) flexible belt of waist; (9) waist component; (10) thigh component; (11) flexible belt of thigh; (12) knee joint component; (13) calf component; (14) flexible belt of calf; (15) pedal; (b). Exoskeleton mechanical diagram.

2.2 Dynamic modeling

Lower limb exoskeleton robots are nonlinear, strongly coupled and multiple DoFs. In order to facilitate the completion of the center of gravity calculation and spatial state description, the structure needs to be appropriately simplified. The human reference coordinate system was established based on the human anatomical posture. The exoskeleton robot was simplified to a five-link model, just as shown in Fig. 2. The white circle is the joint and the black solid circle is the CoM of the linkage.

Figure 2.

Five-link model.

The lower limb exoskeleton robot positional state is represented by a Cartesian coordinate system and generalized coordinates. The robot is supposed to walk in the sagittal XZ plane, and the ankle joint coordinates of the support leg and swing leg are (x _a , z _a ) and (x _b , z _b ), respectively. (x _ci , z _ci ) is the CoM coordinate of the linkage i, i = 1, 2,…, 5. l _i is the length of each linkage of the lower limb; d _i is the distance from the CoM of each linkage to the joint. m _i is the sum of the mass of linkage i and its corresponding human lower limb. θ ₁ , θ ₂ ,…, θ ₅ are the angles between each linkage and its vertical direction, and clockwise direction is positive. M ₁ , M ₂ ,…, M ₅ are the joint torques of each linkage.

According to the geometric relationship in Fig. 2, the CoM coordinates of each rod are

(1) \begin{equation} \left\{\begin{array}{l} x_{ci}=\sum _{j=1}^{i-1}\left(s_{j}l_{j}\sin \theta _{j}\right)+d_{i}\sin \theta _{i}+x_{a}\\[8pt] z_{ci}=\sum _{j=1}^{i-1}\left(s_{j}l_{j}\cos \theta _{j}\right)+d_{i}\cos \theta _{i}+z_{a} \end{array}\right. \end{equation}

where i = 3, s _j = 0; i = 1, 2, 4, 5, s _j = 1.

Differentiating Eq. (1) with respect to time t, the velocity of each linkage CoM is obtained:

(2) \begin{equation} \left\{\begin{array}{l} \dot{x}_{ci}=\sum _{j=1}^{i-1}\left(s_{j}l_{j}\cos \theta _{j}\right)+d_{i}\dot{\theta }_{i}\cos \theta _{i}\\[8pt] \dot{z}_{ci}=\sum _{j=1}^{i-1}\left(-s_{j}l_{j}\sin \theta _{j}\right)-d_{i}\dot{\theta }_{i}\sin \theta _{i} \end{array}\right. \end{equation}

The total potential energy E _p of the lower limb exoskeleton robot system is

(3) \begin{equation} E_{p}=\sum _{i=1}^{5}m_{i}gz_{ci}=\sum _{i=1}^{5}\left\{m_{i}g\left[\sum _{j=1}^{i-1}\left(s_{j}l_{j}\cos \theta _{j}\right)+d_{i}\cos \theta _{i}+z_{a}\right]\right\} \end{equation}

The total kinetic energy E _k of the system is

(4) \begin{equation} \begin{array}{l} E_{k}=\sum _{i=1}^{5}\left\{\dfrac{1}{2}\left[m_{i}\left({\dot{x}^{2}}_{ci}+{\dot{z}^{2}}_{ci}\right)+I_{i}{\dot{\theta }^{2}}_{i}\right]\right\}=\sum _{i=1}^{5}\left\{\dfrac{1}{2}m_{i}\left[\sum _{j=1}^{i-1}\left(s_{j}l_{j}\cos \theta _{j}\right)+d_{i}\dot{\theta }_{i}\cos \theta _{i}\right]^{2}\right\}+\\[14pt] \sum _{i=1}^{5}\left\{\dfrac{1}{2}m_{i}\left[\sum _{j=1}^{i-1}\left(-s_{j}l_{j}\sin \theta _{j}\right)-d_{i}\dot{\theta }_{i}\sin \theta _{i}\right]^{2}\right\}+\sum _{i=1}^{5}\left\{\dfrac{1}{2}m_{i}{d^{2}}_{i}{\dot{\theta }^{2}}_{i}\sum _{j=1}^{i-1}\left[s_{j}l_{j}\dot{\theta }_{j}\cos \left(\theta _{i}-\theta _{j}\right)\right]\right\} \end{array} \end{equation}

where I _i is the rotational inertia of the linkage i with respect to the CoM.

Substituting Eqs. (3) and (4) into the Lagrangian equation of the dynamical state of the lower limb exoskeleton robot system, the joint torque M _i acting on joint i during robot walking is obtained from Eq. (5).

(5) \begin{equation} M_{i}=\frac{d}{d_{t}}\frac{\partial E_{k}}{\partial \dot{\theta }_{i}}-\frac{\partial E_{k}}{\partial \theta _{i}}+\frac{\partial E_{p}}{\partial \theta _{i}},\quad i = 1,2,\ldots,5 \end{equation}

3.1. Stability criteria

In order to comprehensively measure the degree of stability during robot motion, the ZMP stability margin J _z was used to characterize the degree of stability in the horizontal direction, and the expression is shown in Eq. (8).

(8) \begin{equation} J_{z}=a\sum _{i=1}^{n}r_{x}\left(i\right)+b\sum _{i=1}^{n}r_{y}\left(i\right) \end{equation}

where $r_{x}\left(i\right), r_{y}\left(i\right)$ denote the distance from the ZMP point to the center of the foot in the x-axis and y-axis, respectively. n is the total number of the rod of exoskeleton robot. a and b are the weighting coefficients in the X-axis and Y-axis, respectively, and $a+b=1, \frac{a}{b}=\frac{L_{f}}{W_{f}}$ . L _f denotes foot length, and W _f denotes foot width, just as shown in Fig. 3.

Figure 3.

ZMP stability margin diagram.

In the absence of external disturbance, exoskeleton robots mainly rely on the frictional torque to offset the influence of the deflection torque and maintain the balance of the torque in the vertical direction. Considering that the static friction torque is complicated to model and difficult to obtain accurately, the square of the deflection torque is used to characterize the stability in the equilibrium direction.

(9) \begin{equation} J_{b}=\left\| M_{b}\right\| _{2}^{2} \end{equation}

where M _b is the deflection torque.

The ZMP is an important criterion used to check the stability of the robot. According to the definition of ZMP, the ZMP criterion can only ensure that the horizontal component of the torque around ZMP is 0, but its vertical component is generally not 0. As the torso and body parts will generate the swaying torque around the support leg during the forward motion, when the swaying torque exceeds the maximum frictional torque on the ground, the robot will rotate around the support leg, and the deflection caused by the excessive swaying torque is more obvious, so it is necessary to consider the torque balance between horizontal and vertical directions in the gait planning process. The M _b is defined as:

(10) \begin{equation} M_{b}=\sum _{i=1}^{n}m_{i}\ddot{y}_{i}\left(x_{i}-x_{zmp}\right)-\sum _{i=1}^{n}m_{i}\ddot{x}_{i}\left(y_{i}-y_{zmp}\right) \end{equation}

(11) \begin{equation} x_{zmp}=\frac{\sum _{i}m_{i}\left(\ddot{z}_{i}+G_{\mathrm{g}}\right)x_{i}-\sum _{i}m_{i}\ddot{x}_{i}z_{i}}{\sum _{i}m_{i}\left(\ddot{z}_{i}+G_{\mathrm{g}}\right)} \end{equation}

(12) \begin{equation} y_{zmp}=\frac{\sum _{i}m_{i}\left(\ddot{z}_{i}+G_{g}\right)y_{i}-\sum _{i}m_{i}\ddot{y}_{i}z_{i}}{\sum _{i}m_{i}\left(\ddot{z}_{i}+G_{g}\right)} \end{equation}

where (x _i , y _i , z _i ) is the coordinates of the CoM position of the i-th linkage, m _i is the mass of the i-th linkage, n is the total number of links in the robot model and (x _zmp , y _zmp ) denotes the coordinates position of the ZMP.

3.2. Power consumption analysis

In terms of energy consumption, the minimum input energy of the driving joint is taken as the goal for the optimal solution of gait parameters. The average power method is an important measure of energy consumption analysis. Considering that the power of the machine is the product of motor drive torque and joint angular velocity, the average power P _ave of the robot during a walking cycle can be characterized as:

(13) \begin{equation} P_{ave}=\frac{1}{T}\sum _{i=1}^{5}\int _{0}^{T}\left| M_{i}\left(t\right)\cdot \dot{\theta }_{i}\left(t\right)\right| dt \end{equation}

where M _i (t) is the torque, $\dot{\theta }_{i}\left(t\right)$ is the joint velocity, and T is the gait cycle.

In the process of robot motion, the exoskeleton robot may be affected by external disturbances and its own structural deformation, resulting in the instantaneous power approaching infinity, but the average power at this time may be very small. The instantaneous power peak will reduce the stability of the entire system. In response to this situation, this paper introduces the root mean square deviation of the instantaneous power, which could accurately reflect the discrete degree between instantaneous power and average power.

(14) \begin{equation} P_{i}=\sum _{i=1}^{5}\left| M_{i}\left(t\right)\cdot \dot{\theta }_{i}\left(t\right)\right| \end{equation}

(15) \begin{equation} \sigma =\sqrt{\frac{1}{\left(T\right)}\int _{0}^{T}\left(P_{i}-P_{ave}\right)^{2}dt} \end{equation}

where P _i is the instantaneous power. σ is root mean square deviation of instantaneous power.

3.3. The multi-objectives function design

We constructed the overall multi-objectives function $F\!\left(\mathbf{x}\right)$ for optimizing the exoskeleton robot system.

(16) \begin{equation} F\!\left(\mathbf{x}\right)=\alpha J_{z}+\beta J_{b}+\mu \left(P_{ave}+\sigma \right) F\!\left(\mathbf{x}\right)=\alpha J_{z}+\beta J_{b}+\mu \left(P_{ave}+\sigma \right) \end{equation}

where x is the vector, and $\mathbf{x}=\left[S,\rho _{1},\rho _{2},\rho _{3},\theta _{2},\theta _{3},\theta _{4},\theta _{5}\right]$ . α, β, μ∈ [0,1], the value is set according to the proportion of gait ZMP stability margin and driving energy consumption. The smaller the $F\!\left(\mathbf{x}\right)$ , the more stable and smooth the robot gait trajectory. In this study, the ideal gait could be obtained through optimizing the gait parameters. Therefore, multi-objectives function $F\!\left(\mathbf{x}\right)$ can be expressed as the optimized object function GoalF with gait parameters, and its expression is

(17) \begin{equation} \textit{GoalF}=\min \left[F\!\left(\mathbf{x}\right)\right] \end{equation}

4.1. Overview of BSO-EOLLFF algorithm

Based on the beetle antennae search algorithm and particle swarm optimization (PSO) algorithm, the BSO algorithm was designed [Reference Zhang and Zhang21, Reference Jiang and Li22]. The BSO algorithm has the characteristics of fast solving speed and high accuracy and has been successfully applied in the fields of signal positioning [Reference Liu, Qian and Jia23] and data classification [Reference Wu, Ma, Xu, Li and Chen24]. However, the BSO algorithm is easy to fall into the local optimum when searching in a high-dimensional space. In this regard, this paper introduced the EOL strategy, the LFF strategy and the dynamic mutation strategy to improve the optimization performance. In this study, we used EOL to initialize the population to promote a uniform distribution of initial population information and improve search efficiency. In addition, the LFF strategy enables the individuals to learn both their own and group experiences, allowing each individual to move purposefully and instructively and improving the convergence performance of the algorithm. Finally, a dynamic mutation strategy was introduced to increase the population diversity at the end of the iteration and prevent the algorithm from falling into the local optimum.

4.2. Population initialization based on EOL strategy

Tizhoosh introduced the concept of opposition-based learning in 2005 and showed that the inverse solution has a nearly 50% higher probability of being close to the global optimum than the current solution [Reference Tizhoosh25]. The main idea was that by generating reverse individuals in the area where the current individual was located and by competing for the reverse individuals with the current individuals, the best individuals will enter the next generation.

Definition 3. (Elite inverse solution). Let the inverse solution of an elite individual $X_{best}=\left(x_{o1},x_{o2},\cdots,x_{on}\right)$ in a population p be $X'_{\!\!best}=\left(x'_{\!\!o1},x'_{\!\!o2},\cdots,x'_{\!\!on}\right)$ in an n-dimensional search space, and the elite inverse solution could be defined as $x'_{\!\!oi}=\lambda *\left(de_{i}+du_{i}\right)-x_{oi}$ . $\left[de_{i},du_{i}\right]$ is the dynamic boundary of the population P in the i-th dimensional search space and could be calculated according to (18).

(18) \begin{equation} de_{i}=\min _{1\leq j\leq \left| P\right| }\left(x_{oi}\right), du_{i}=\max _{1\leq j\leq \left| P\right| }\left(x_{oi}\right) \end{equation}

where $\left| P\right|$ is the size of the current group.

The dynamic boundary of the search space was used instead of the fixed boundary to preserve the search experience so that the generated inverse solutions could be located in the gradually shrinking search space and promote the faster convergence of the algorithm. Since the elite inverse solution may also jump out of the boundary $\left[e_{i},u_{i}\right]$ and become non-feasible, the method in the literature [Reference Zhou, Wu and Wang26] was used here to reset the transgressed values, as shown in (19):

(19) \begin{equation} x'_{\!\!oi}=rand\left(e_{i},u_{i}\right), e_{i}\lt x'_{\!\!oi} \;\;\textrm{or}\;\; u_{i}\gt x'_{\!\!oi} \end{equation}

By forming the elite inverse solution, the detection of elite individuals’ neighborhoods can be enhanced. The BSO algorithm performs reverse learning for the elite individuals of each iteration to generate the reverse population of elite individuals and participate in the evolutionary competition.

4.3. LFF strategy

In the standard BSO algorithm, the search range of beetle individuals was limited, and it is difficult to transfer the search position from the global optimal to the local optimal. Although the group search could expand the search range, there is no information exchange and feedback between the individuals, which will affect the convergence accuracy of the algorithm. The individuals in the group need to continuously learn the experience of the historical information, which plays a decisive role in the improvement of the algorithm’s convergence speed. Therefore, the LFF strategy was introduced. LFF strategy was proposed by Levy in the 1930s, and Mandelbrotb described it in detail [Reference Jiadong, Xinqian, Qian, Haitao and Xiaolin27]. The strategy obeys the Levy distribution. LFF strategy could increase the diversity of the population, expand the search range and prevent the algorithm from falling into the local optimum. The Levy distribution satisfies:

The calculation formula of the step length s is

(20) \begin{equation} s=\frac{\psi }{\left| \nu \right| ^{\frac{1}{\varphi }}} \end{equation}

where $\psi$ and ν obey normal distribution, $\mu,v\sim N\left(0,\sigma ^{2}\right)\mu \sim N\left(0,\sigma _{\mu }^{2}\right), \nu \sim N\left(0,\sigma _{\nu }^{2}\right)$ . $\varphi$ is constant

(21) \begin{equation} \sigma _{\psi }=\left[\frac{\Gamma \left(1+\varphi \right)\ast \sin \dfrac{\pi \varphi }{2}}{\Gamma \left(\dfrac{1+\varphi }{2}\right)\ast \varphi \ast 2^{\frac{\varphi -1}{2}}}\right]^{\frac{1}{\varphi }} \end{equation}

where Γ() is the standard gamma distribution.

The update rule of beetle direction is

(22) \begin{equation} d\!\left(t+1\right)=\omega \ast d\!\left(t\right)+C_{1}\ast Levy\!\left(\theta \right)\ast \left(\textit{gbest}\left(t\right)-X\!\left(t\right)\right)+C_{2}\ast Levy\!\left(\theta \right)\ast \left(\textit{zbest}\!\left(t\right)-X\!\left(t\right)\right) \end{equation}

The update rule of beetle position is

(23) \begin{equation} X\!\left(t+1\right)=X\!\left(t\right)+k_{1}*s*d\!\left(t\right)*sign\!\left(f\!\left(X_{r}\!\left(t\right)\right)-f\!\left(X_{l}\!\left(t\right)\right)\right)+k_{2}d\!\left(t\right) \end{equation}

where d(t) represents the movement direction of the t-th generation. X(t) represents the beetle position of the t-th generation. gbest() represents the individual extreme value of the t-th generation beetle. zbest() represents the global extreme value. ω is the inertia weight, and C ₁ and C ₂ are the learning factors. Levy(θ) is a Levy random number. f(x _r (t)) and f(x _l (t)) are the fitness function value of the right whisker and left whisker of the t-th generation longhorn beetle. k ₁ and k ₂ are the scale factor. sign() is a sign function.

4.4. Dynamic mutation strategy

In the later stage of the iteration, the population diversity of the BSO algorithm will become lower, which will reduce the search ability of the algorithm. To avoid premature phenomena, a dynamic mutation strategy was introduced to increase the diversity of the population and improve convergence accuracy. At present, related scholars have proposed a variety of mutation algorithms, and the typical ones are Gaussian mutation [Reference Ni, Ling, He and Hai28] and Cauchy mutation [Reference Rudolph29]. Compared with the Gaussian mutation, the Cauchy mutation could generate a more extensive range of random numbers so that the algorithm has a greater chance to jump out of the local optimum. Meanwhile, when the peak value is small, the Cauchy mutation only takes less time to search nearby areas. The Cauchy mutation was selected for the second optimization of the population, and the mutation operation was performed on X(t):

(24) \begin{equation} X^{*}\!\left(t\right)=X\!\left(t\right)+\zeta *C\!\left(0,1\right) \end{equation}

(25) \begin{equation} \zeta =e^{-\gamma \frac{t}{H}} \end{equation}

where ζ is the weight of mutation, which will decrease with increasing iteration. H is the maximum number of iteration. γ is constant. C(0,1) is a random number generated by Cauchy operator.

4.5. The basic steps and verification of BSO-EOLLFF algorithm

Because of a large number of parameters to be optimized for the objective function in this paper, the process of gradually approximating the optimal value is complicated, the accuracy of the optimized solution is not high, and the convergence speed is relatively slow. During the search process of the traditional PSO algorithm, the position of each particle is approached toward the global optimal position in a decreasing sine wave manner. During this process, a particle converges to the currently found global optimum position, and the other particles will quickly converge to this local optimum solution. If this optimal position is only a local optimal point, then most of the particles are confined by this local optimal solution. To overcome the shortcomings of traditional algorithms, the BSO-EOLLFF optimization algorithm is proposed in this paper. According to the gait parameters to be optimized, the population individuals are selected as $x_{p}=\left[S,\rho _{1},\rho _{2},\rho _{3},\theta _{2},\theta _{3},\theta _{4},\theta _{5}\right]$ in the improved algorithm. Since the ankle joint has a passive DoF, θ ₁ is not used as an optimization parameter. According to Eq. (16), the fitness value corresponding to the multi-objectives function $F\!\left(\mathbf{x}\right)$ at each point x _p was calculated.

4.5.1. The steps of BSO-EOLLFF algorithm

The specific steps of the BSO-EOLLFF algorithm are:

1. 1. Initialize the algorithm parameters: the beetle scale, iterative step length and the maximum number of iterations. Use the EOL strategy to initialize the beetle population and initialize the direction.

2. 2. Calculate the corresponding fitness function value of the beetle and determine the individual extreme value and the global extreme value according to the fitness function value.

3. 3. Use (22) and (23) to update the direction and position of the beetle individuals and perform cross-border processing of beetle populations.

4. 4. Use (24) to perform mutation operation on population.

5. 5. Determine whether the algorithm meets the iteration termination condition. If it is satisfied, output the optimal global solution and its corresponding position, otherwise return to step 2.

The pseudo code of BSO-EOLLFF algorithm is shown in Table I.

Table I.

Psedo code of BSO-EOLLFF algorithm.

Table II.

Parameters setting of the two algorithms.

4.5.2. Performance analysis of BSO-EOLLFF algorithm

In order to verify the performance of the BSO-EOLLFF algorithm, the Rastrigin function $z=-20e^{-0.2\sqrt{0.5\left(x^{2}+y^{2}\right)}}-e^{0.5\cos \left(2\pi x\right)+\cos \left(2\pi y\right)}+20+2.72$ with lots of local solutions was selected to perform function optimization and convergence test. The parameter setting of the two algorithms is shown in Table II. The test results are shown in Fig. 4. The BSO-EOLLFF algorithm could accurately search for the optimal solution and would not fall into the local optimization. The algorithm’s search time is 1.35 s. To verify the superiority of the BSO-EOLLFF algorithm, the PSO algorithm was used to perform parameter search for the same test function, and the experimental results are shown in Fig. 5. The experimental results showed that although the PSO algorithm can eventually obtain the global optimal solution as well, there are four local optimal solutions in the process of finding the optimal solution, and the time to find the optimal solution is 5.73 s, which is significantly slower than the algorithm proposed in this study. The experimental results showed that the algorithm proposed in this study has a significant effect enhancement for the regular test function.

Figure 4.

Optimization results of BSO-EOLLFF algorithm for the test functions. (a) Test function diagram, (b) optimal trace of each test function, (c) the contour map of the search path and (d) the convergence curves.

Figure 5.

Optimization results of PSO algorithm on the test functions. From left to right: optimal trace of each test function, the contour map of the search path, and the convergence curves.

5.1. Parameter setting

The mechanism parameters of the lower limb exoskeleton robot are set as: m ₁ = 2 kg, m ₂ = 2.5 kg, L ₁ = 0.37 m, L ₂ = 0.4 m. The number of beetles n is 50, the maximum number of iterations H is 1000, inertia weight w is 0.8, individual learning factor C ₁ is 1.8, and global learning factor C ₂ is 1.4. The ratio k ₁ and k ₂ are 0.4 and 0.6, respectively. The initial weight of mutation is 0.7. The initial step length s ₀ is 1. The gait cycle T is 1 s. The constant γ and $\varphi$ are 10 and 1.5, respectively. In the iterative process, the condition for the iteration termination is that the error of parameter optimization is less than 0.1. In order to avoid the chance of experimental results, we conducted five simulation experiments and took the average of the experimental results for analysis.

5.2. Convergence analysis

In order to avoid the BSO algorithm from falling into the local optimum, the EOL algorithm with global search capability, the LFF algorithm with local search capability and the dynamic mutation strategy with high population diversity were introduced to improve the optimization performance. To verify the convergence of the BSO-EOLLFF algorithm, the convergence curves of each algorithm were compared with the PSO algorithm and BSO algorithm, respectively, as shown in Fig. 7.

Figure 7.

Convergence curve of each algorithm.

It can be seen from Fig. 7 that the traditional PSO algorithm stops searching after 40 iterations and falls into the local optimum due to its weak global search ability [Reference Zhang and Zhang21]. When the BSO algorithm iterates 23 times, the function fitness value tends to be stable. However, it has three local optimal solutions (gray square), which will affect the convergence. Although it is lower than the global optimal value found by the PSO algorithm, the convergence speed is slower. The BSO-EOLLFF algorithm searches for the optimal solution when iterated 8 times. The BSO-EOLLFF could find the only global optimal value, which is the lowest. And the averaged convergence time is shown in Table III. The time consumption of the algorithm is 9.76 s, which is significantly faster than the other two algorithms. After the comparative analysis of the search time, the efficiency of the BSO-EOLLFF algorithm was not weakened after the addition of the EOL algorithm and LFF algorithm, which indicates that the fit and conversion between the algorithms are good. The strategy proposed in this paper increases the diversity of the population, which speeds up the algorithm search speed and further improves the optimization accuracy.

5.3. Stability analysis

To further verify the stability of the BSO-EOLLFF algorithm proposed in this paper, the gait ZMP stability margin and joint driving energy consumption of the lower limb exoskeleton robot optimized by the three algorithms were compared and analyzed. The exoskeleton robot’s ZMP stability margin curves were optimized using the PSO algorithm, the BSO algorithm and the BSO-EOLLFF algorithm, as shown in Fig. 8.

It can be seen from Fig. 8, in one gait cycle, after the optimization of the PSO algorithm, the gait ZMP of the robot only passes through the stable center-line once. The trajectory of the center of gravity of the robot deviates from the bipedal support area, and the human body is prone to fall when walking. After optimizing by the BSO algorithm, the gait ZMP trajectory is close to the boundary. The minimum distance between the gait ZMP trajectory and the upper boundary is 0.0063 m. There are two times passing the stability center-line. The robot’s center of gravity trajectory is roughly distributed in the middle of the bipedal support domain, which has an ideal stability margin. The ZMP optimized by the BSO-EOLLFF algorithm has better performance. The closest gait ZMP point to the boundary line of the support domain is 0.1272m. At this time, the minimum distance between the two is 0.0228 m. The gait ZMP trajectory is uniformly in the middle of the support domain, which can realize stable walking.

5.4. Energy consumption analysis

As can be seen from Fig. 9, the drive energy consumption of the robot increases and then decreases with each walking step, which is because the magnitude of the torque required by the joints in the oscillation process is mainly determined by the work done by gravity. To thoroughly verify the effectiveness of the proposed optimization algorithm, the minimum, maximum and average values of the drive energy consumption of the robot in 100 walking cycles optimized by the three algorithms were selected for comparison, and the algorithm energy consumption is shown in Table IV.

Table III.

The averaged convergence time.

Figure 8.

ZMP stability margin curve after algorithm optimization.

Figure 9.

Gait drive energy consumption curve.