Adaptive Particle Swarm Optimization of PID Gain Tuning for Lower-Limb Human Exoskeleton in Virtual Environment

Abstract

Tuning of a proportional-integral-derivative (PID) controller for a complex multi-joint structure, such as an exoskeleton, using conventional methods is difficult and imprecise. In this paper, an optimal PID tuning method for a 3-dimensional model of a lower-limb human exoskeleton in gait training condition is presented. The dynamic equation of the human-exoskeleton is determined using a Lagrangian approach, and its transfer function is established in a closed-loop control system. PID controller gains, initialized by the Ziegler–Nichols (Z-N) method, are used as the input to an adaptive particle swarm optimization (APSO) algorithm for minimizing the multi-joint trajectory error. The optimized controller is tested in the Gazebo virtual environment and compared with the Z-N and conventional optimization methods. The numerical analysis shows that the PID controller tuned by a combination of Z-N and APSO improves the performance of a lower-limb human exoskeleton in gait training.

1. Introduction

Stroke and spinal cord injury are the main causes of disability in the elderly [1,2], creating problems in daily life. For many years, rehabilitation therapies have been concentrated on recovering their motion function using various manual methods. In recent years, the lower-limb exoskeleton, which is a wearable robot, has been used to aid patients with mobility problems [3,4,5].

For the exoskeleton wearer to walk properly, the desired trajectory of the joints requires a robust control system to minimize steady-state errors to reduce disturbances. Huo et al. [6] proposed an active impedance control strategy for controlling a lower-limb exoskeleton for the sit-to-stand motion. Many researchers focused on adaptive-based control approaches to design a controller with uncertainties and non-linear systems. For instance, Daachi et al. [7] controlled one degree of freedom (DoF) lower-limb knee joint orthosis using a neural network based on adaptive control using the Lyapunov approach for rehabilitation reasons. Chen et al. [8] studied an adaptive robust cascade force control for a one DoF exoskeleton with a hydraulic actuator to minimize the effect of uncertainties. Lu et al. [9] described an adaptive learning control system with unknown parameters using the Lyapunov approach for stability. In another study, Yang et al. [10] designed a robust learning controller with hybrid electro-hydraulic actuators, and verified their exoskeleton using a periodic reference.

For many years, proportional-integral-derivative (PID) controllers have been commonly used because of their reliable performance, and convenient and robust operation [11,12,13]. In our approach, we used particle swarm optimization (PSO) to tune the PID controller’s gains. PSO has a simple concept and is easy to implement in different optimization problems. Hence, it has been used to control the robustness of PID parameters [14]. PSO has few parameters that need to be adjusted, and it has high probability and efficiency in finding the global optima, but initializing its parameters can be difficult. Although PSO has fast error convergence, it can converge prematurely and be trapped in the local optima [15]. Therefore, researchers have combined different tuning methods to optimize its gains. For example, Belkadi et al. [16] developed an optimal method to obtain the parameters of the PID controller and verified the performance in a simulation using PSO with random initialization. Similarly, Durgadevi et al. [17] used PSO for tuning the parameters of proportional-integral (PI) and PID controllers of the temperature of a continuous stirred tank reactor.

As an alternative to PSO, other optimization methods have been used for tuning controller parameters [18]. Qianhao et al. [19] developed the whale optimization algorithm, which is a type of swarm intelligence optimization algorithm, to improve the performance of a linear–quadratic regulator controller. Pan et al. [20] presented the fuzzy-PID algorithm for the nonlinear dynamics and uncertainties of the human exoskeleton robot. Their method was verified in a co-simulation technique combining ADAMS and MATLAB/Simulink. Moharam et al. [21] used a hybrid differential evolution and PSO algorithm for tuning a PID controller and compared it with the genetic algorithm (GA) and the Ziegler–Nichols (Z-N) method to prove the hybrid method is robust to disturbances. Lal et al. [22] implemented the grasshopper optimization algorithm (GOA) to obtain the PID and fuzzy-PID controller parameters for a multi-area interconnected microgrid power system. They compared two different PID controllers and concluded that fuzzy-PID produced better dynamic responses compared to a PID controller tuned by GOA.

This paper presents the tuning of PID controller gains by combining Z-N and adaptive particle swarm optimization (APSO). A simulation and the dynamic equation of the human exoskeleton model walking with the aid of a parallel bar in a virtual environment are described. The transfer functions of each joint based on the dynamic equation of human exoskeleton model are established in a closed-loop control system. An objective function formulated by the error of the control system is defined. Our main objective was to optimize a PID controller using initialization with Z-N and optimized by APSO. The tuned PID was verified in the 3-dimensional virtual simulation environment and compared with the conventional method.

2. Human Exoskeleton Simulation in a Virtual Environment

The lower-limb exoskeleton model consists of a waist, left and right femurs, left and right tibias, and left and right feet. The human exoskeleton simulated model was created and controlled in a robot operating system (ROS) with average human anthropometric height. The waist functions as a base link to which the left and right femurs are connected. Consequently, other links are assembled to each other as parent and child links [23]. Figure 1a,b shows the hierarchical relation of the 3D exoskeleton model in an open-source dynamic simulator that imitates the behavior in a real environment.

To simulate the interaction between the exoskeleton and patient, a 3D human model was added in the Gazebo simulation environment. The foot, tibia, and femur of human are coupled to its related links. The configurations of coupling between human and exoskeleton model are shown in Figure 2.

To ensure the simulation accurately mimicked to the gait training task, parallel bars were added to avoid the human model from falling down by providing stability. Both hands of the human model are constrained to the parallel bars using a slide joint to allow forward and backward motion. Figure 3 demonstrates the gait training simulation of the human exoskeleton model walking with both hands holding the parallel bars.

3. Human Exoskeleton Dynamic Modeling

Figure 4 illustrates a free body diagram of the exoskeleton and human model. $O_1$, $O_2$, and $O_3$ are the center position of hip, knee, and ankle joint, respectively; ${\theta }_1$, ${\theta }_2$, and ${\theta }_3$ represent angles of hip, knee, and ankle, respectively; $l_1$, $l_2$, and $l_3$ are length of the femur, tibia, and foot, respectively; and $l_{Gi}$ is the center of gravity (CoG) of links 1, 2, and 3.

In this paper, each link of the exoskeleton and human is modeled as one DoF pendulum model, which replicates one joint of the exoskeleton and human, as shown in Figure 5. The summation of the physical effect of each lower link such as mass and inertia is considered to be $M_i$ [24].

In this study, the Lagrangian formulation of dynamic equation of the pendulum model is expressed as follows:

$$L_i=E_{ki}-E_{pi},$$

(1)

$$T_i=\frac{d}{dt}\left(\frac{d{L}_i}{d\dot{{\theta }_i}}\right)-\frac{d{L}_i}{d{\theta }_i},$$

(2)

where i represents number of the links; $L_i$ is the Lagrangian function; ${\theta }_i$ and $\dot{{\theta }_i}$ are the angular position and velocity, respectively; $T_i$ is the torque that applies to each joint; and $E_{ki}$ and $E_{pi}$ represent kinetic and potential energy [25], respectively, which are given as follows:

$$E_{ki}=\frac{1}{2}{I}_i{\dot{{\theta }_i}}^2+\frac{1}{2}{m}_i({\dot{x_{Gi}}}^2+{\dot{y_{Gi}}}^2)+\frac{1}{2}{M}_i({\dot{x_i}}^2+{\dot{y_i}}^2),$$

(3)

$$E_{pi}=m_{i}g{y}_{Gi}+M_{i}g{y}_i,$$

(4)

where $m_i$ and $I_i$ are the mass and inertia of the link, respectively; $M_i$ is the mass of lower links; g represents gravitational acceleration; $x_{Gi}$, $y_{Gi}$, $x_i$, $y_i$$\dot{x_{Gi}}$, $\dot{y_{Gi}}$$\dot{x_i}$, and $\dot{y_i}$, represent the position of the CoG, end point of the link, and their velocity, respectively, as follows:

$$x_{Gi}=l_{Gi}sin{\theta }_i{\dot{x}}_{Gi}=l_{Gi}\dot{{\theta }_i}cos{\theta }_i{x}_i=l_{i}sin{\theta }_i\dot{x_i}=l_i\dot{{\theta }_i}cos{\theta }_i,$$

(5)

$$y_{Gi}=l_{Gi}cos{\theta }_i{\dot{y}}_{Gi}=-l_{Gi}\dot{{\theta }_i}sin{\theta }_i{y}_i=l_{i}cos{\theta }_i\dot{y_i}=-l_i\dot{{\theta }_i}sin{\theta }_i,$$

(6)

where $l_{Gi}$ and $l_i$ are the length of the CoG and the link i, respectively. The total kinetic and potential energies are determined as follows:

$$E_{ki}=\frac{1}{2}{I}_i{\dot{\theta }}_i^2+\frac{1}{2}{m}_i{l_{Gi}}^2{\dot{{\theta }_i}}^2+\frac{1}{2}{M}_i{l}_i^2{\dot{{\theta }_i}}^2,$$

(7)

$$E_{pi}=gcos{\theta }_i(m_i{l}_{Gi}+M_i{l}_i).$$

(8)

The dynamic equation of the pendulum model is determined as follows:

$$T_i=(m_i{l_G}^2+M_i{l}_i^2+I_i)\ddot{{\theta }_i}+b\dot{{\theta }_i}+(l_i{M}_i+l_{Gi}{m}_i)gsin{\theta }_i,$$

(9)

where b is the friction coefficient and $sin{\theta }_i$ is approximated as ${\theta }_i$. The transfer function $G\left(s\right)$ is represented in Laplace as follows:   

$$G\left(s\right)=\frac{{\theta }_i\left(s\right)}{T_i\left(s\right)}=\frac{1}{s^2+\frac{b}{(m_i{l_{Gi}}^2+M_i{l}_i^2+I_i)}s+\frac{(l_i{M}_i+l_{Gi}{m}_i)g}{(m_i{l_{Gi}}^2+M_i{l}_i^2+I_i)}}.$$

(10)

Table 1. Physical features of the exoskeleton and human model.

Exoskeleton	Femur	Tibia	Foot
m (kg)	0.2943	0.2159	0.1150
l (m)	0.41	0.39	0.195
I (kg m2)	$5.7\times {10}^{-3}$	$4.3\times {10}^{-3}$	$3.9\times {10}^{-4}$
$l_G$ (m)	0.15	0.11	0.04
Human	Femur	Tibia	Foot
m (kg)	7.33	3.4503	1.075
l (m)	0.407	0.4334	0.275
I (kg m2)	0.1502	0.0505	0.0038
$l_G$ (m)	0.1763	0.1849	0.1179

lists the physical features of the exoskeleton based on computer-automated design (CAD). The human model’s features were obtained from References [26,27,28]. The mass and inertia of the human model were for a subject that weighed 74.2 kg and was 175.5 cm tall.

Based on the features, the transfer function of the exoskeleton ($G_L\left(s\right)$) and human ($G_H\left(s\right)$) can be determined as shown in

Table 2. Transfer functions of the human–exoskeleton model.

Link	Exoskeleton	Human Model
Femur	$\frac{1}{s^2+13.25s+25.96}$	$\frac{1}{s^2+27.27}$
Tibia	$\frac{1}{s^2+36.88s+27.58}$	$\frac{1}{s^2+29.24}$
Foot	$\frac{1}{s^2+1568s+78.62}$	$\frac{1}{s^2+66.34}$

.

4. Control System

The plant of the system $G\left(s\right)$ is defined as a combination of transfer functions of the human $G_H\left(s\right)$ and exoskeleton $G_L\left(s\right)$ in parallel. Figure 6 shows the control system with the feedback loop for the joint angle for each joint. Since the human model and exoskeleton are coupled together, it is assumed they are parallel to each other and the joints’ angles are equal. The PID controller employs the input as the error between yhe actual and desired trajectories, and the output is the torque as a plant [29]. The error is defined as follows:

$$e={\theta }_{act}-{\theta }_{des}.$$

(11)

The PID controller and plant model are represented as follows:

$$C\left(s\right)=K_p+\frac{K_i}{s}+K_{d}s,$$

(12)

$$G\left(s\right)=G_L\left(s\right)+G_H\left(s\right),$$

(13)

where $C\left(s\right)$ is the PID controller; $K_P$, $K_i$, and $K_d$ are the proportional, integral, and derivative gains to be tuned to reduce the error, respectively [30]; and $G_L\left(s\right)$ and $G_H\left(s\right)$ are the transfer function of the exoskeleton and lower limb of the human, respectively.

5. Optimization of the Controller’s Gains

The PID controller is tuned using an algorithm by a combination of the Z-N method and APSO. The PID parameters are obtained from Z-N, which is the classical observation method of tuning [31]. The output is set as the initial value, which provides the searching space to avoid local optima [32]. APSO increases the speed of convergence and the chance of finding the global optimum. PSO is based on population, inspired from observation of biological social phenomena such as a flock of birds, swarms of bees, and schools of birds [33].

An objective function $f_{obj}$ of the optimization problem is defined as the integral time absolute error (ITAE), which is the summation of the absolute error weighted by the elapsed time [34,35], as follows:

$$f_{obj}={\int }_0^{\infty }t\left|e_{obj}\left(t\right)\right|dt,$$

(14)

where the error of control system is given as follows:

$$e_{obj}\left(s\right)=1-\frac{C\left(s\right)G_i\left(s\right)}{1-C\left(s\right)G_i\left(s\right)},$$

(15)

where t is the elapsed time; $e_{obj}\left(t\right)$ and $e_{obj}\left(s\right)$ represent the error of the control system in the time and frequency domains, respectively. The tuning algorithm starts with Z-N, where $K_i$ and $K_d$ are set to zero while the $K_p$ value is established as $K_u$. The value of $K_u$ is added until the first oscillation is observed, then its frequency is calculated as $T_u$. The gains of the PID are determined as follows:

$$K_p=0.6\times K_u,$$

(16)

$$K_i=\frac{1.2·K_u}{T_u},$$

(17)

$$K_d=\frac{3·K_u·T_u}{40},$$

(18)

$$x_{zn}=\left[K_p{K}_i{K}_d\right],$$

(19)

where $x_{zn}$ is the output that is input into APSO, which contains particles as design variables for minimizing the objective function. Initialization is created by randomizing $x_{zn}$ as follows:

$$x_{1j}=rand[x_{zn}-ϵ,x_{zn}+ϵ],$$

(20)

where $x_{1j}$ represents the particles in the first iteration; $ϵ$ is a positive value, which is selected as a random value between 0 and 1.

$$ϵ=rand(0,1)$$

(21)

The particles of the next iteration are a summation of the positions of previous particles with their velocities:

$$x_{k,j}=x_{k-1,j}+v_{k,j},$$

(22)

where k and j represent the numbers of iterations and particles, respectively; $x_{k,j}$ and $x_{k-1,j}$ are the position of the current particle and the particle of the previous iteration, respectively; and $v_{k,j}$ is the velocity and direction of each particle toward the next iteration as follows:

$$v_{k,j}={\omega }_k·v_{k-1,j}+C_1·ran{d}_1·(p_{bes{t}_k-1}-x_{k,j})+C_2·ran{d}_2·(g_{bes{t}_k}-x_{k,j}),$$

(23)

where $ran{d}_1$ and $ran{d}_2$ are random values between 0 and 1; ${\omega }_i$ is a positive coefficient, called inertia weight, which reduces linearly with iteration, and is determined as follows:

$${\omega }_i={\omega }_{max}·b,$$

(24)

where ${\omega }_{max}$ is the maximum value for inertia weight; b is a gain that starts from 1 in the first iteration and end to 0, which allows ${\omega }_i$ starts from the maximum value to have wider searching space and decrease to be narrowed toward the global optima [36]; and b is a positive descending coefficient, starting from 1 in the first iteration. b is represented as follows:

$$b=1-a,$$

(25)

where a is a positive ascending coefficient, which starts from 0 and increases to 1 at the last iteration, as follows:

$$a=\frac{k}{k_{max}}k=1,2,...,k_{max},$$

(26)

where k and $k_{max}$ are the maximum number of iterations; $C_1$ and $C_2$ are non-negative coefficients of self and social recognition, defined as follows:

$$C_1=1.8·b+2.1,$$

(27)

$$C_2=1.8·a+0.1.$$

(28)

The value of $C_1$ is greater than $C_2$, and their summation remains at 4 in all iterations [36]. In each iteration, the objective function is determined for each particle, and the smallest objective function is selected as $p_{best}$. The value of $g_{best}$ will be the lowest value of $p_{best}$ when the algorithm ends. ASPO stops when the number of iterations meets the maximum number of iterations $k_{max}$. In this paper, the numbers of iterations and particles are 150 and 30, respectively, which were determined by experimentation to obtain the convergence of the objective function values to zero. Algorithm 1 and Figure 7 illustrate the pseudo code and flow chart of the APSO algorithm, respectively, where $K_{step}$ is a constant value for a step of $K_u$.

Algorithm 1 Combination of Ziegler–Nichols (Z-N) and adaptive particle swarm optimization (APSO).

1: Start  2:    3: Set $K_i=0$;  4:    5: Set $K_d=0$;  6:    7: Set $K_p=K_u$;  8:    9: Initialize $K_u$; 10:   11: Set step value for $K_u$ 12:   13: while Observe oscillation do; 14:   15:     Set $K_{step}$; 16:   17:     $K_u=K_u+K_{step}$; 18:   19: end while 20:   21: Calculate frequency of the oscillation as $T_u$; 22:   23: Calculate $K_p$, $K_i$, and $K_d$ 24:   25: Initialize population around the output of Z-N results; 26:   27: Evaluate initial population; 28:   29: Selecting the initial $p_{best}$; 30:   31: while Number of iteration less than $i_{max}$ do; 32:   33:     Create new iteration; 34:   35:     Select the $p_{best}$; 36:   37:     if $p_{best}$ less than the previous $g_{best}$ then; 38:   39:         Store it as $g_{best}$ 40:   41:     end if 42:   43: end while 44:   45: Select the particle of $g_{best}$ as the result 46:   47: End

Figure 8 shows the $p_{best}$ in each iteration for three joints of the human exoskeleton.

The changes in $C_1$ and $C_2$ for 150 iterations are shown in Figure 9 and the $\omega$ in each iteration.

Table 3. Gains of the proportional-integral-derivative (PID) controller.

	${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{d}}$
Hip	1656	4731.41	24.15
Knee	1314	3864.62	0.5
Ankle	101.90	615.57	4.26

demonstrates the optimal gains of the PID controller.

6. Results and Discussion

The desired trajectory of each joint is based on gait training to move the joints simultaneously. The trajectories of each joint are driven by two different functions, repeated periodically.

Table 4. Equation of desired trajectory for each joint.

Joints	$0<\mathit{t}<\frac{{\mathit{t}}^{\prime }}{2}$	$\frac{{\mathit{t}}^{\prime }}{2}\le \mathit{t}<{\mathit{t}}^{\prime }$
Right hip	$-0.24sin\left(\pi t\right)$	$0.4sin\left(\pi t\right)$
Right knee	0	$-0.48sin\left(\pi t\right)$
Right ankle	$-0.04sin\left(\pi t\right)$	$0.4sin\left(\pi t\right)$
Left hip	$0.4sin\left(\pi t\right)$	$-0.24sin\left(\pi t\right)$
Left knee	$-0.48sin\left(\pi t\right)$	0
Left ankle	$0.4sin\left(\pi t\right)$	$-0.04sin\left(\pi t\right)$

shows the desired trajectory.

t represents the elapsed time and $t^{\prime }$ is the time of one walking cycle. Figure 10 provides the frame-by-frame pictures of gait training with parallel bars.

The APSO algorithm was coded in Python programming language, and the virtual environment was Gazebo, integrated with ROS. All tests were performed using a PC with a Ubuntu 18.04 LTS operation system with 16 GB RAM and an Intel core-i3 Xeon processor. Figure 11 exhibits the performance of the simulated model in Gazebo. The graphs compare the trajectories of Z-N, PSO, and APSO with the desired trajectory. The Z-N measured the initial values for PSO and APSO. The difference between PSO and APSO is the method of setting the parameters. In PSO, $\omega$, $C_1$, and $C_2$ are adjusted to 1.2, 2, and 2, respectively.

Figure 12 shows the angular error trajectory for the hip, knee, and ankle. For all joints, the APSO error was the lowest, followed by PSO and Z-N. For example, the average error determined by APSO for the hip was 24% and 48% is less than PSO and Z-N, respectively. APSO showed less overshoot than the other methods because of its smaller peak, whereas PSO and Z-N experienced more shock when the trajectory changed.

Figure 13 demonstrates the output torque of each link. The maximum torques of the knee and ankle were less than 3 and 0.05 Nm, respectively. The maximum amplitude of torque for the hip was 7 Nm, which is the highest because the actuator of the hip tolerates a heavier mass compared to the knee and ankle.

Table 5. Numerical analysis.

Joint	Left				Right
	$\mathit{M}\mathit{E}$	$\mathit{A}\mathit{E}$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	${\mathit{R}}^{\mathbf{2}}$	$\mathit{M}\mathit{E}$	$\mathit{A}\mathit{E}$	$\mathit{R}\mathit{M}\mathit{S}$	${\mathit{R}}^{\mathbf{2}}$
Hip	0.0285	0.0080	0.0096	0.9983	0.0184	0.0077	0.0096	0.9983
Knee	0.0092	0.0021	0.0024	0.9998	0.0037	0.0019	0.0022	0.9998
ankle	0.023	0.0042	0.0047	0.9991	0.0088	0.0044	0.0049	0.9990

provides the numerical analysis of Z-N and APSO performance.

$ME$, $AE$, $RMSE$, and $R^2$ represent the maximum error, average error, root mean square error, and coefficient of determination, respectively. The maxima of $ME$ and $AE$ were less than 0.03 m, which are within the acceptable zone [37]. The maximum of $RMSE$ was 0.0130, while the $R^2$ values of all joints were greater than 99%, indicating the close fitted of all variability in the actual trajectories.

7. Conclusions

In this study, a PID controller was optimally tuned using an optimization algorithm that combines of Z-N and APSO. Parameter changes in each iteration enabled the finding of the optimal solution and avoided the local optima. Furthermore, the PID gains were established by Z-N and used for initializing the APSO. The proposed method was compared and tested in a virtual environment on a human exoskeleton model. The results demonstrated the acceptable performance of the APSO controller by statistical analysis.

The optimal PID controller minimized the trajectory error of a multi-joint human exoskeleton structure for gait training. For a new trajectory and different training conditions, the optimization process should be repeated, which is the limitation of this study. In future work, the model can be applied in different situations such as climbing stairs and the stand-to-sit movement.