1. Introduction
The exoskeleton robot can fully adapt to the physiological structure of the human body by being worn on the outside of the human body, moving in coordination with the human to achieve various purposes, such as exercise assistance, ability improvement, and medical rehabilitation [Reference Brahmi, Saad, Rahman and Ochoa-Luna1–Reference Lee, Liu, Perroud, Chavarriaga and Millana4]. With the development of lower limbs exoskeleton robots, the demand for high-precision dynamic models of the exoskeleton is increasing. Therefore, it is necessary to identify the dynamic parameters when the robot is working [Reference Anwar and Al Jumaily5, Reference Mallat, Bonnet, Huo, Karasinski, Amirat, Khalil and Mohammed6].
The main task of dynamic parameters identification is to obtain a parameter set that integrates various impact factors through experimental means so that the torque calculated by the parameter set is consistent with the actual torque required by the robot. The dynamic parameter identification mainly includes dynamic modeling, linearization, identification trajectory optimization, identification algorithm construction, parameter acquisition and processing, and experimental verification [Reference Jianxiong7]. The design of the identification algorithm is the most critical. There are many mature research identification algorithms, such as genetic algorithm identification [Reference Jina, Tsai and Kao8], maximum likelihood estimation identification [Reference Swevers, Verdonck, Naumer and Pieters9], least square method identification [Reference Kinsheel and Taha10, Reference Jonker, Zheng and Aravkin11]. Due to the dynamic coupling of different joints, the identification parameters are not independent. When performing parameter identification, the mutually independent parameters must be obtained firstly. The independence among the identification parameters is determined by the minimum parameter set. M. Gautier has proposed methods with deriving the minimum parameter set of tandem robots [Reference Gautier and Khalil12, Reference Gautier13]. To make identification parameter in a good ‘incentive’ state, it is necessary to select the reasonable identification trajectory. Chen proposed a dynamic parameter identification method to analyze the human-machine coupling relationship [Reference Jing, Yuping, Xibin, Shit, ong and Longfei14]. Kyung-Jo Park [Reference Park15] and Azeddien Kinsheel [Reference Kinsheel and Taha10] applied this method to optimize the trajectory coefficients for different identification objects and further improved the identification results. Swevers uses the maximum likelihood estimation method to calculate the value of the identification parameter [Reference Swevers, Verdonck, Naumer and Pieters9]. The least square method is a traditional parameter identification method. Zhang Tie used Newton Euler’s method to establish the SCARA robot dynamics model and used the least square method to identify high-precision robot dynamics parameters [Reference Zhang, Liang and Qin16]. The least square method has high identification accuracy, but it cannot meet the real-time requirements when the amount of calculation is large [Reference Zhiyu, Bin and Chaohui17]. Ding Li identified the dynamic parameters of industrial robots with an identification algorithm based on the combination of the weighted least squares method and artificial bee colony algorithm and verified the accuracy of the algorithm in torque identification through experiments [Reference Li, Hongtao and Yu18]. BLEEX exoskeleton used the least square method for parameter identification. This algorithm could obtain the inertia parameter, the joint friction parameter, damping parameter, and stiffness parameter [Reference Ghan, Steger and Kazerooni19]. However, the parameter identification method is very cumbersome, which makes the process time-consuming. The dynamic model of the lower extremity exoskeleton was established based on Lagrange method, and the unknown parameters in the dynamic model are identified by the least square method to obtain a more accurate system dynamics model [Reference Pei, Pei, Gu and Tang20]. Compared with the traditional form, the recursive form has higher computational efficiency. Deep learning has powerful capabilities for data analysis and mining and is now gradually being used to solve robot-related problems. Bargsten used the dynamic parameters as the training weights of the back-propagation neural network and obtained the optimal dynamic parameter values through the feedback adjustment of the network training [Reference Bargsten, De Gea Fernandez and Kassahun21]. Liu established the relationship between robot joint motion and torque through long short-term memory technology and achieved good results in predicting the torque of robot joint motion [Reference Liu22]. The long short-term memory technology was combined with an optimized recurrent neural network which is proposed to compensate for the dynamic model of the proposed 6-DoF collaborative robot based on the consideration of gravity, coriolis force, inertial force, and friction force [Reference Gong, Zheng, Zhicheng and Wenlin23]. The algorithm based on long short-term memory algorithm and recurrent neural network model has a good prediction on the actual torque. And the root-mean-square (RMS) error between predicted and actual torques is reduced by 61.8–78.9% compared to the traditional dynamic model. Chen deduced the dynamics equation of the robot servo motor and used a neural network for dynamic compensation [Reference Enwei24]. However, the deep learning network requires a large amount of data, and the recognition accuracy largely depends on the quality of the data [Reference Rueckert, Nakatenus and Tosatto25].
The optimization algorithm based on swarm intelligence is also the main research object currently used to identify parameters [Reference Verdel, Bastide, Vignais, Bruneau and Berret26, Reference Chang, An and Ma27]. Chang modeled legged robots by intelligent decoupling modeling method. This method groups the DoF according to the number of links and treats each link as a tandem manipulator structure, which greatly reduces the number of identified dynamic parameters at the same time [Reference Chang, An and Ma27]. The representative one is the PSO algorithm proposed by Kennedy and Eberhart [Reference Kennedy and Eberhart28]. The research results show that the standard PSO algorithm can effectively identify the parameters of the Bouc-Wen model under high-noise conditions [Reference Ye and Wang29]. The modified PSO algorithm was proposed by some researchers to implement some meaningful novelties. Wang Yueling used an improved particle swarm genetic algorithm to identify the dynamics of flexible joint robots. The algorithm adopts a dynamic adaptive adjustment strategy, which effectively improves the convergence speed of the particle swarm algorithm [Reference Yueling, Yue and Qi30]. Sheng Wentao proposed a parameter identification method that combines recursive least squares and PSO. Compared with the identification algorithm for BLEEX, the method has higher accuracy and lower complexity [Reference Wentao31]. Most dynamic parameter identification methods adopt the identification strategy for the combined parameter’s value [Reference Qin, Baron and Birglen32–Reference Goher Khaled and Fadlallah Sulaiman34]. Qin proposed a novel sequential identification procedure to identify the dynamic parameters of robotic manipulators. The method improved the identification precision compared to conventional methods [Reference Qin, Baron and Birglen32]. The methods are only suitable for robot control with the constant load. When the load changes, the robot system must re-identify the parameter set and adjust the robot control model, which is time-consuming and labor-intensive. Therefore, this study combines group optimization algorithm with beetle foraging mechanism, and a novel BSO algorithm was designed and applied to the dynamic parameter identification. The structure of the paper is shown below. Section 2 describes the proposed exoskeleton robot. Section 3 describes the specific content of dynamic parameter identification. Section 4 describes the test results of the proposed BSO algorithm for multi-modal function and the actual exoskeleton robot model. Section 5 is the discussion.
2. Exoskeleton robot system
2.1. The proposed 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 DoF for the hip joint and knee joint, respectively, and one passive DoF for the ankle joint. The thigh rod and calf rod are adjustable, from 402 to 505 mm and 397 to 506 mm, respectively. The motion range of the hip joint is 2.17°–50.18°. The motion range of the knee joint is 6.85°–95.53°.
Figure 1. Overall structure diagram (1) back support; (2) driving source of hip joint; (3) hip joint component; (4) brace of thigh; (5) driving source 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.
The identification object in this paper is a 6-DOF exoskeleton robot. The D-H method was used to establish a joint coordinate system for the robot. The coordinate system of each joint is shown in Fig. 2, and the D-H parameters of each joint are shown in Table I. a i represents the length of connecting rod. d i represents the distance between adjacent coordinate systems. α i represents the angle between the connecting rod and the x-axis.
Table I. D-H parameters.
Figure 2. Lower limb exoskeleton model.
2.2. Exoskeleton robot control system
The exoskeleton robot control system is composed of a host computer controller, joint drive control unit, and power supply. Each joint drive control unit is connected through the bus and communicates with the controller, as shown in Fig 3. Each joint drive control unit includes a servo motor, a harmonic reducer, an encoder, a drive board, and a control board. The motion and current information collected by the sensors were feedback to the controller. For safety reasons, an emergency stop button is added between the power supply and the controller to perform emergency stop operations when the exoskeleton robot behaves abnormally.
Figure 3. Control system.
2.3. Dynamic modeling
For a multi-joint tandem robot, the joints are numbered from 1 to n in turn. The joint torque of the i-th joint
${}^i\tau $
can be modeled by Newton’s Euler method.
\begin{align}{}^i\tau = M(\theta )\theta '' + C(\theta ,\theta ')\theta ' + G(\theta ) + {}^i{\tau _f}\end{align}
where
$\theta $
,
$\theta '$
, and
$\theta ''$
represent the joint angle, angular velocity, and angular acceleration, respectively. M(θ) is inertia matrix.
$C(\theta ,\theta ')$
is Coriolis force and Centrifugal force matrix, G(θ) is gravity vector.
${}^i{\tau _f}$
is the friction torque. Joint friction is a complex nonlinear factor, but for most robot systems, a simplified linear model can be used to describe it.
\begin{align}{}^i{\tau _f} = {}^i{f_c}{\mathop{\rm sgn}} (\theta ') + {}^i{f_v}\theta '\end{align}
where the friction parameter
${}^i{p_f}$
is included in the term
${}^i{\tau _f}$
, and
${}^i{p_f} = {({}^i{f_c},{}^i{f_v})^T}$
.
${}^i{f_c}$
and
${}^i{f_v}$
are coulomb and viscous friction coefficients, respectively. sgn () represents symbolic function.
3. Dynamic parameter identification
There are many factors that affect robot dynamics, but many related parameters are not provided by the manufacturer. And it must be considered in the control system. Therefore, it is necessary to identify the dynamic parameters of the robot and obtain a set of comprehensive parameter values that include various influencing factors.
Each joint of the exoskeleton robot has 10 inertial parameters, and the parameters of the i-th joint are represented by a vector as:
\begin{align}{}^iP = {({}^i{I_{XX}},{}^i{I_{XY}},{}^i{I_{XZ}},{}^i{I_{YY}},{}^i{I_{YZ}},{}^i{I_{ZZ}},{}^i{H_X},{}^i{H_Y},{}^i{H_Z},{}^im)^T}\end{align}
where
${}^iP$
is the parameter to be identified for dynamic parameter identification.
${}^i{I_{XX}}$
∼
${}^i{I_{ZZ}}$
are the six parameters of the inertia matrix.
${}^i{H_X}$
,
${}^i{H_Y}$
, and
${}^i{H_Z}$
are the three components of
${}^iH$
, and
${}^iH = {}^im \times {}^i{r_c} = {}^im({}^i{r_{cX}},{}^i{r_{cY}},{}^i{r_{cZ}})$
,
${}^i{r_c}$
is the centroid vector,
${}^im$
is weight parameter, and is included in
$G(\theta )$
.
3.1. Model linearization
The dynamic model obtained by the Newton-Euler method is not linear. The joint drive torque expression contains
${}^im{}^ir_c^2$
, that is
${}^iH{}^i{r_c}$
, which is caused by
${}^i{r_{ci}} \times {}^i{f_{ci}}$
. The parallel axis theorem is applied to convert the inertia matrix
${}^iI$
from the relative centroid description to the relative coordinate origin description, which can eliminate the nonlinear term and obtain a linear model.
The parallel axis theorem is:
\begin{align}{}^{oi}{I_i} = {}^{ci}{I_i} + {m_i}({}^ir_{ci}^T \cdot {}^i{r_{ci}} \cdot I - {}^ir_{ci}^T \cdot {}^i{r_{ci}})\end{align}
where
${}^{oi}{I_i}$
is the inertia matrix described by relative centroid.
${}^{ci}I$
is the inertia matrix described by relative ordinate origin.
Then, the dynamic equation is converted to
\begin{align}{}^i{\tau _i} = {}^{oi}{I_i}{}^i{\omega ^{\prime}_i} + {}^i{\omega _i} \times ({}^{oi}{I_i}{}^i{\omega _i}) + {m_i}{}^i{r_{ci}} \times {}^i{v^{\prime}_i} + {}_i^{i + 1}{R^{i + 1}}{\tau _{i + 1}} + {}^i{p_{i + 1}}{}_{\phantom{+1}i}^{i + 1}{R^{i + 1}}{f_{i + 1}}\end{align}
The
${}^iH{}^i{r_c}$
term is not included in the dynamic equation anymore. And the Eq. (1) could convert to the linear form for the identified parameter.
\begin{align}{}^i\tau & = {}^i{c_1}{}^i{I_{XX}} + {}^i{c_2}{}^i{I_{XY}} + \cdot \cdot \cdot + {}^i{c_{10}}{}^im + \cdot \cdot \cdot + {}^n{c_1}{}^n{I_{XX}} + {}^n{c_2}{}^n{I_{XY}} + \cdot \cdot \cdot + {}^n{c_{10}}{}^nm \nonumber\\[3pt] &= {}^i{c_{iner}}^T \cdot {}^i{p_{iner}}{\rm{ + }}{}^i{c_f}^T \cdot {}^n{p_f} + {}^{i + 1}{c_{iner}}^T \cdot {}^{i + 1}{p_{iner}} + \cdot \cdot \cdot + {}^n{c_{iner}}^T \cdot {}^n{p_{iner}} + {}^n{c_f}^T \cdot {}^n{p_f}\nonumber\\[3pt] &\quad {}^{i + 1}{p_{iner}} \sim {}^n{p_{iner}}\end{align}
where
${}^i{c_{iner}}^T$
is the coefficient vector of the inertial parameter of the i-th joint,
${}^i{c_{iner}}^T = {({}^i{c_1},{}^i{c_2}, \cdot \cdot \cdot ,{}^i{c_{10}})^T}$
. Each component in the vector
${}^i{c_{iner}}^T$
is the function of joint angle
$\theta $
, angular velocity
$\theta '$
, and angular acceleration
$\theta ''$
.
${}^i{c_f}$
is the coefficient vector of the friction parameter of the i-th joint,
${}^i{c_f} = {({\mathop{\rm sgn}} (\theta '),\theta ')^T}$
. In the process of dynamic parameter identification, these items could be calculated by collecting the joint motion parameters.
3.2. Parameter independence analysis
According to Eq. (6), due to the coupling of serial robot system dynamics, the identified parameters may not be independent. The torque
${}^i\tau $
is not only affected by
${}^i{p_{iner}}$
but also by the inertial parameters of
${}^{i + 1}{p_{iner}} \sim {}^n{p_{iner}}$
.The mutual independence of parameters is the prerequisite for obtaining higher identification accuracy. Therefore, a set of mutually independent parameters, namely the minimum parameter set, must be found. According to the minimum parameter set theory [Reference Gautier and Khalil12], the minimum parameter set corresponding to the 10 inertial parameters of any joint is
\begin{align}{}^i{P_m} = {({}^i{I_{XX}}{\rm{ - }}{}^i{I_{YY}},{}^i{I_{XY}},{}^i{I_{XZ}},{}^i{I_{YZ}},{}^i{I_{ZZ}},{}^i{H_X},{}^i{H_Y})^T}\end{align}
where the three parameters
${}^i{I_{YY}}$
,
${}^i{H_Z}$
, and
${}^im$
are unidentifiable torque parameters for
${}^i\tau $
. It must be identified by the torque
${}^{i - 1}\tau $
of the i-1th joint. However, the parameters
${({}^i{I_{YY}},{}^i{H_Z},{}^im)^T}$
and the parameters in the minimum parameter set
${}^{i - 1}{P_m}$
do not satisfy mutual independence. M. Gautier uses the method of symbol derivation to give a general recursive formula for
${({}^i{I_{YY}},{}^i{H_Z},{}^im)^T}$
and
${}^{i - 1}{P_m}$
:
\begin{align}{}^{i - 1}{I^m}_{XX} = {}^{i - 1}{I_{XX}} + {}^{i - 1}{I^m}_{YY} + 2{d_i}{}^i{H^m}_Z + 2{d_i}^2{}^i{m^m}\end{align}
\begin{align}{}^{i - 1}{I^m}_{XY} = {}^{i - 1}{I_{XY}} + {a_i}\sin {\alpha _i}({}^i{H^m}_Z + {d_i}{}^i{m^m})\end{align}
\begin{align}{}^{i - 1}{I^m}_{XZ} = {}^{i - 1}{I_{XZ}} + {a_i}\cos {\alpha _i}({}^i{H^m}_Z + {d_i}{}^i{m^m})\end{align}
\begin{align}{}^{i - 1}{I^m}_{YY} = {}^{i - 1}{I_{YY}} + {\cos ^2}{\alpha _i}({}^i{I^m}_{YY} + 2{d_i}{}^i{H^m}_Z + 2{d_i}^2{}^i{m^m}) + {a^2}_i{}^i{m^m}\end{align}
\begin{align}{}^{i - 1}{I^m}_{YZ} = {}^{i - 1}{I_{YZ}} + \cos {\alpha _i}\sin {\alpha _i}({}^i{I^m}_{YY} + 2{d_i}{}^i{H^m}_Z + 2{d_i}^2{}^i{m^m})\end{align}
\begin{align}{}^{i - 1}{I^m}_{ZZ} = {}^{i - 1}{I_{ZZ}} + {\sin ^2}{\alpha _i}({}^i{I^m}_{YY} + 2{d_i}{}^i{H^m}_Z + 2{d_i}^2{}^i{m^m}) + {a^2}_i{}^i{m^m}\end{align}
\begin{align}{}^{i - 1}{H^m}_Z = {}^{i - 1}{H_Z} + \cos {\alpha _i}({}^i{H^m}_Z + {d_i}{}^i{m^m})\end{align}
\begin{align}{}^{i - 1}{H^m}_X = {}^{i - 1}{H_X} + {a_i}{}^i{m^m}\end{align}
\begin{align}{}^{i - 1}{H^m}_Y = {}^{i - 1}{H_Y} - \sin {\alpha _i}({}^i{H^m}_Z + {d_i}{}^i{m^m})\end{align}
\begin{align}{}^{i - 1}{m^m} = {}^{i - 1}m + {}^i{m^m}\end{align}
After combination, the parameters
$({}^i{I^m}_{XX}{\rm{ - }}{}^i{I^m}_{YY},{}^i{I^m}_{XY},{}^i{I^m}_{XZ},{}^i{I^m}_{YZ},{}^i{I^m}_{ZZ},{}^i{H^m}_X,{}^i{H^m}_Y)$
were mutual independent. The superscript m on the right represents the parameter in the smallest parameter set.
${a_i}$
,
${d_i}$
, and
${\alpha _i}$
represent the D-H parameter of the i-th joint.
3.3. Dynamic parameter identification algorithm
3.3.1. Particle swarm optimization algorithm.
The PSO algorithm finds the optimal solution through the movement of the particle swarm. Each particle has two properties, velocity, and position, and the optimal solution can be searched separately. The individual particle will share information with the others to find the globally optimal solution.
In the D-dimensional space, the position and velocity of particles at the t-th iteration are
$x_k^t = (x_{k1}^t,x_{k2}^t, \cdot \cdot \cdot ,x_{kd}^t, \cdot \cdot \cdot ,x_{kD}^t)$
and
$v_k^t = (v_{k1}^t,v_{k2}^t, \cdot \cdot \cdot ,v_{kd}^t, \cdot \cdot \cdot ,v_{kD}^t)$
. The optimal position of individuals and groups are
$P_k^t = (P_{k1}^t,P_{k2}^t, \cdot \cdot \cdot ,P_{kd}^t, \cdot \cdot \cdot ,P_{kD}^t)$
and
$P_G^t = (P_{G1}^t,P_{G2}^t, \cdot \cdot \cdot ,P_{Gd}^t, \cdot \cdot \cdot ,P_{GD}^t)$
. Update of the d-th dimension position and velocity of the k-th particle are
\begin{align}x_{kd}^{t + 1} = x_{kd}^t + v_{kd}^{t + 1}\end{align}
\begin{align}v_{kd}^{t + 1} = \omega v_{kd}^t + {c_1}rand(P_{kd}^t - x_{kd}^t) + {c_2}rand(P_{gd}^t - x_{kd}^t)\\[-5pt] \nonumber \end{align}
where
$\omega$
is inertia weight, c
1
is individual learning factor, c
2
is global learning factor, rand () is random number between (0,1).
This paper adopts the strategy of decreasing inertia weight, and the formula is as follows:
\begin{align}\omega {\rm{ = }}{\omega _{\max }} - \frac{{{\omega _{\max }} - {\omega _{min}}}}{M}t\end{align}
where
$\omega$
max and
$\omega$
min represent the maximum and minimum value of
$\omega$
. M is the maximum number of iterations.
The PSO algorithm has the advantages of fast convergence speed and simple calculation, but it is easy to fall into the optimal local solution. BAS was proposed in 2017 [Reference Xiangyuan and Shuai35–Reference Wu, Ma, Xu, Li and Chen38]. The algorithm is designed to imitate the beetle relying on the strength of the food smell to forage. The simplified model of the beetle is shown in Fig. 4.
Figure 4. Simplified model of beetle.
In the D-dimensional space, the center position of the beetle at the t-th iteration is
${x^t} = ({x_1},{x_2}, \cdot \cdot \cdot ,{x_D})$
. The left whisker and right whisker of beetle are defined as
$x_l^t$
and
$x_r^t$
, and
\begin{align}\left\{ {\begin{array}{*{20}{c}}{x_l^t = x + {l^t}\mathop d\limits^ \to }\\ \\[-7pt] {x_r^t = x - {l^t}\mathop d\limits^ \to }\end{array}} \right.\end{align}
where
${l^t}$
is the distance from the center of the beetle to the whisker at the t-th iteration, and the expression is:
\begin{align}{l^t} = {\eta _l}{l^{t - 1}} + {l_0}\end{align}
$\mathop d\limits^ \to $
is a random unit vector, and the expression is:
\begin{align}\mathop d\limits^ \to = \frac{{rand(D,1)}}{{\left\| {rand(D,1)} \right\|}}\end{align}
where
${\eta _l}$
is the attenuation coefficient of the search distance.
${l_0}$
is a constant.
$rand(D,1)$
is D-dimensional vector composed of random numbers between -1 and 1.
The position of the beetle will be updated by comparing the scent of the left and right whisker.
\begin{align}{x^{t + 1}} = {x^t} - {\delta ^t}\mathop d\limits^ \to * {\mathop{\rm sgn}} (f(x_l^t) - f(x_r^t))\end{align}
where
$f(x_l^t)$
and
$f(x_r^t)$
are the fitness function value of the left whisker and right whisker, respectively. sgn() is symbolic function.
${\delta ^t}$
is step size of the t-th iteration.
\begin{align}{\delta ^t} = {e_1} * {\delta ^{t - 1}} + {e_0}\end{align}
where e 1 is attenuation coefficient of step size, e 0 is adjustment parameters of step size.
\begin{align}{l^t} = \frac{{{\delta ^t}}}{{{e_2}}}\end{align}
where e 2 is the relationship parameter between the step length and the distance from the centroid of longhorn beetle to the whisker.
The biggest advantage of the BAS algorithm is simplicity and high efficiency. However, the algorithm ignores the information sharing between individuals, and the identification results for high-dimensional complex optimization problems are not ideal.
3.3.2. Dynamic parameter identification based on BSO algorithm
There is a population of N beetles represented as X=(X
1
, X
2
, …, X
N
), the maximum number of iterations is M, the initial iteration step length of beetle is
${\delta ^0}$
, the initial distance between the centroid of beetle and the whisker is
${l^0}$
, the inertia weight is
$\omega$
, the learning factors is c
1
, c
2
, the initial position is
$x_k^0 = (x_{k1}^0,x_{k2}^0, \cdot \cdot \cdot ,x_{kd}^0, \cdot \cdot \cdot ,x_{kD}^0)$
, initial velocity is
$v_k^0 = (v_{k1}^0,v_{k2}^0, \cdot \cdot \cdot ,v_{kd}^0, \cdot \cdot \cdot ,v_{kD}^0)$
, D is the dimension of the vector.
By calculating the value of fitness function
$f(x_k^t)$
of each beetle, the current optimal position of the individual
$P_k^t = (P_{k1}^t,P_{k2}^t, \cdot \cdot \cdot ,P_{kd}^t, \cdot \cdot \cdot ,P_{kD}^t)$
and the optimal position of the group
$P_G^t = (P_{G1}^t,P_{G2}^t, \cdot \cdot \cdot ,P_{Gd}^t, \cdot \cdot \cdot ,P_{GD}^t)$
are obtained.
The position coordinates of the left and right whiskers of beetle are
$x_{kl}^t$
and
$x_{kr}^t$
. By calculating the values of fitness function and
$f(x_{kr}^t)$
of the two whiskers of each beetle, the speed of the BAS is updated. The update rule is Eq. (27).
\begin{align}\left\{ {\begin{array}{*{20}{c}}{x_{kl}^t = x_k^t + v_k^t{l^t}}\\ \\[-7pt] {x_{kr}^t = x_k^t - v_k^t{l^t}}\end{array}} \right.\end{align}
\begin{align}\psi _k^{t + 1} = - {\delta ^t}v_k^t * {\mathop{\rm sgn}} (f(x_{kl}^t) - f(x_{kr}^t)) \\[8pt] \nonumber \end{align}
After updating the position separately according to the methods of PSO and BAS, the new position is weighted to allow for the variations of the parameter set. The direction of whisker and step length of BAS and the speed of PSO are used to derive the updated position of the beetle. The update rule is:
\begin{align}x_i^{t + 1} = x_i^t + \lambda v_k^{t + 1} + (1 - \lambda )v_{kd}^{t + 1}\end{align}
where λ represents the weight of the PSO algorithm.
In this study, the flow chart of BSO algorithm is shown in Table II.
4. Results and discussion
4.1. Experimental parameter settings
The number of particle swarms N is 20, the maximum number of iterations M is 1000, inertia weight w is 0.8, individual learning factor c 1 is 1.4, global learning factor c 2 is 1.4. The initial step length is 1, the initial distance from the centroid of beetle to the whisker is 0.1, the attenuation coefficient e 1 of the step length is 0.95, and the step length adjustment parameter e 0 is 0.005. The relationship parameter e 2 between the distance from the center of mass of the beetle to the whisker is 10. In the iterative process, the condition for the iteration termination is that the error of parameter identification is less than 0.1.
4.2. Comparative analysis of different identification algorithms in four modal test functions
In order to verify the performance of the BSO algorithm, the test function f i (x) with an optimal solution was selected to perform function optimization and convergence test. To avoid a one-sided conclusion, the different functions were used to test the performance of the proposed method. The test functions are given in Table III.
Table II. The pseudocode of the BSO algorithm.
Table III. The different modal test functions.
Figure 5. Behavior of BSO 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.
This paper selects uni-modal function, multi-modal function, and fixed-dimensional multi-modal function to observe the optimization behavior of the BSO algorithm, respectively. The Fig. 5 shows the contour map of the optimal trajectory, search path, and convergence curve of each test function. The optimal trajectory gives the optimal route of the beetles. Since the initial position of the beetle is randomly generated, the optimal trajectory may be different when copying the result. The contour map of the search path can more intuitively display the trajectory of the beetle and connect the same z value on the x and y planes, making it easier to observe the movement of the beetle. The convergence curve shows the function value of the best solution obtained so far. The beetles gradually move to the optimal point and eventually gather near the global optimal point. The BSO algorithm firstly initializes a set of random solutions. The search unit updates its location based on its own search mechanism and the best solution currently available at each iteration. The BSO algorithm with the information sharing mechanism can avoid local optima more intelligently. The combination of these two algorithms can accelerate the iteration speed of the population and reduce the probability of the population falling into local optimum, making it more stable in dealing with high-dimensional problems.
It can be seen from Table IV that when processing the test function, the processing speed of BSO is equivalent to that of PSO, but it is significantly better than the BSA algorithm. In addition, compared with the other two algorithms, the performance of the BSO algorithm is more stable. In the optimization process, the information sharing search mechanism is added, so that the algorithm has better global optimization performance, speeds up the convergence speed of the algorithm, and effectively avoids the phenomenon of local optimization.
Table IV. Comparison of optimization results of three algorithms for different test functions.
Figure 6. Human wearing exoskeleton in ADAMS.
4.2. BSO for dynamic parameters identification
In order to verify the superiority of the BSO algorithm proposed in this research in the identification of dynamic parameters, the traditional one-time identification method is selected for comparative analysis. The one-time identification method is widely applied in robotic field [Reference Jing, Yuping, Xibin, Shit, ong and Longfei14, Reference Park15, Reference Wentao31]. In the process of parameter identification, all joints are moved together, and then, the currents and rotation angles of all joints are collected and brought into the overall identification model, and all the parameters to be identified are calculated at one time.
Simulation adopts the method by combining ADAMS with Mathematica. The dynamics simulation of the robot is carried out in ADAMS, the human-machine model is shown in Fig. 6. The 3D model of the human wearable exoskeleton is imported into Adams. The material properties, constraints, and angular displacement driving functions of the exoskeleton and human model are given, respectively. The lower extremity exoskeleton model is shown in Fig. 6. The human body model was established by Solidworks. According to China’s national adult human body size standard GB10000-88, the height of the established human body model is 1750 mm, and the quality is 75 kg. During a gait cycle, the maximum force on the thighs and calf is 950 and 560 N, respectively. Through the dynamic analysis of the lower extremity exoskeleton system, the driver of the hip joint and knee joint adopted the servo motor with a rated power of 400 W. The length of the ball screw is 4 mm, and the rated thrust is 2000 N. The maximum linear speedup is 100 mm/s. The torque and joint motion parameters obtained from the simulation are output to simulate the sampling values of the parameters in the actual identification experiment. The sample value is substituted into the identification algorithm constructed in Mathematica, and the value of the identification parameter can be solved.
The traditional common parameter identification method is one-time identification method. We use the one-time identification method as the reference object. The real value X rea , identification value X ide and error e of each joint inertia parameter and friction parameter measured by the proposed dynamic parameter identification method and traditional one-time identification method were shown in Tables V and VI, respectively.
Table V. The real value and the identified value of the identification parameters obtained by BSO algorithm.
Table VI. The real value and the identified value of the identification parameters obtained by the traditional one-time identification method.
The error e is defined as
\begin{align}e = \frac{{\left| {{X_{rea}} - {X_{ide}}} \right|}}{{{X_{rea}}}}\end{align}
Since the first joint (hip joint) of the object identified in this paper is parallel to the direction of gravity, there is only one identified parameter for the first joint. Through the comparison of Tables V and VI, the method proposed in this paper can identify the independent values of all the parameters of the second and the third joints. However, the traditional one-time identification method can only identify the minimum parameter set, that is, the combined value of the parameters. This is because the combined value is susceptible to load changes, and it needs to be re-identified when the load changes. The accuracy of the corresponding parameters identified by the BSO method is higher than that of the one-time identification, especially for the parameters of the first joint. The traditional method is to identify all joints together, and the parameter identification value of the i+1∼n th joint needs to be substituted into the parameter identification equation of the i-th joint, and the error accumulation is large, especially when the number of joints is large.
The goal of parameter identification is to accurately predict the torque required for the robot to perform target motion through the parameter identification value. Therefore, the verification of the identification results must be achieved through the accuracy of the torque prediction. According to the human body motion trajectory obtained from the gait experiment, this paper takes the verification trajectory as:
\begin{align}{\theta _{hip}}(t) = 118.7\sin(0.092t + 0.016) + 3.085\sin(0.025t + 1.632) + 111.8\sin(0.0938t + 3.02)\end{align}
\begin{align}{\theta _{knee}}(t) = - 5089{e^{ - {{\left(\frac{{p - 76.25}}{{13.44}}\right)}^2}}} + {\rm{5677}}{e^{ - {{\left(\frac{{p - 75.76}}{{12.88}}\right)}^2}}} + {\rm{5}}{\rm{.3288}}{e^{ - {{\left(\frac{{p - 19.74}}{{23.41}}\right)}^2}}}\end{align}
Under this trajectory, when the simulation time is 4.8 s, the angular acceleration of the robot reaches 5.2 rad/s2. When the robot moves along the verification trajectory Eqs. (31) and (32), the predicted torque value τ P can be calculated from the parameter identification value x ide (Table V) through Eq. (6). Given the verification trajectory for each joint in ADAMS, the real torque value τ R could be obtained. The verification criterion takes the RMS value of τ P and τ R.
\begin{align}{}^i\varepsilon RMS = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{(\tau _P^{(i)} - \tau _R^{(i)})}^2}} } \end{align}
where
$\tau _P^{(i)}$
and
$\tau _R^{(i)}$
represent the i-th sampling value of predicted torque value and real torque value, respectively. The values of RMS for the three joints are
${}^1\varepsilon RMS = 23.85$
,
${}^2\varepsilon RMS = 20.39$
and,
${}^3\varepsilon RMS = 9.26$
, respectively. They account for 7.42%, 6.5%, and 4.91% of the absolute mean of their joint torques, respectively.
It can be seen from Fig. 7 that the τ
P
and τ
R
of each joint are very close, and the RMS value of
${}^i\varepsilon RMS$
(i=1, 2, 3) is smaller than the amplitude of the noise. This indicates that although the identification accuracy of some dynamic parameters is not high, it has little effect on the prediction torque, and the prediction accuracy of torque mainly depends on those parameters with a great contribution. The error of torque is less than the added noise, which indicates that the parameter identification method proposed in this paper can suppress the noise interference to a certain range.
Figure 7. Experimental results of joint torque.
5. Discussion and conclusion
This paper proposed a dynamic parameter identification algorithm based on the optimization of the long-horned beetle herd, which can identify the independent values of all parameters and overcome the shortcomings of traditional identification methods that the identification results are affected by load changes. The proposed identification method greatly reduces the amount of identification calculation, and the identification accuracy of each time is independent of each other, and there is no error accumulation phenomenon. The experimental results show that the 10 dynamic identification parameters can be accurately predicted. The identification method can accurately predict the driving torque of the robot. The RMS values are all smaller than the added noise level, which shows that although the identification accuracy of some dynamic parameters is not high, it has little effect on the predicted torque. The accuracy of torque prediction mainly depends on the parameters that contribute a lot. The error of the torque is smaller than the added noise, which shows that the parameter identification method proposed in this paper has a certain inhibitory effect on noise interference. In future work, we will consider real-time issues to speed up the response speed of the control system.
Conflicts of Interest
The work described has not been submitted elsewhere for publication, in whole or in part, and all the authors listed have approved the manuscript that is enclosed.
Financial Support
This research is partly supported by Baicheng Bayuan Project – nutrition and health key technology and intelligent manufacturing (21ZYQCSY00050) and Tianjin key laboratory of integrated design and online monitoring of light industry and food engineering machinery equipment.
Authors Contributions
All authors have made substantial contribution. The first author Peng Zhang designed the experiment scheme, analyzed the experimental results, and wrote the paper. Corresponding author Junxia Zhang modified the paper.
