Robust PID optimal tuning of a Delta parallel robot based on a hybrid optimization algorithm of particle swarm optimization and differential evolution

Abstract

In this paper, we propose an approach to tune optimal parameters of a robust PID controller by means of computed torque control (CTC) strategy for trajectory tracking of a Delta parallel robot, using a hybrid optimization algorithm of Particle Swarm Optimization (PSO) and differential evolution (DE). It differs from previous works that they propose robust PID controller parameters tuning based on conventional gradient-based optimization algorithms and apply them to process control. First, we reduce the tuning problem of a robust PID controller with CTC strategy satisfying requirements including robustness and disturbance attenuation to an optimization problem with nonlinear constraints by considering the nonlinear dynamic model of a Delta parallel robot. Second, we set up the design characteristics for the trajectory tracking of a Delta parallel robot. Then, we propose a robust PID controller in a way of obtaining the global optimization solution of the nonlinear optimization problem by running a PSO-DE hybrid optimization algorithm of finding the global optimal solution by maintaining the diversity of swarm during evolution based on the evolution of cognitive experience. Simulation and experimental results demonstrate that the proposed controller outperforms previous works with respect to robust performance and active disturbance attenuation when it is applied to tracking control of a Delta parallel robot.

1. Introduction

Recent years have witnessed ever-increasing applications of high-speed pick-and-place parallel robots such as Delta parallel robot in fields of electronics, food, pharmaceutical, packaging and other light industries [1–3]. For the parallel robots, it is difficult to design the control system with high control performance due to the presence of nonlinearity in the dynamic model, the high interaction and the kinematical constraints between the robot links [1].

To enhance the trajectory tracking performance of a robot manipulator, some traditional control strategies have been employed such as the independent PID control on each axis, the PID control with an acceleration feedforward component, the PID control with gravity compensation and nonlinear PID-based control [4, 5]. Since these strategies never or partially use the dynamic model of the Delta parallel robot, dynamic interaction between robot links could not be considered. Hence, one cannot expect the high control performance when high speed and accuracy of a parallel robot are of concern.

Otherwise, the computed torque control (CTC) strategy has been proposed which employs all the dynamic model of the manipulator [6]. However, the CTC strategy provides a good control performance only when the obtained model is accurate enough, and it does not consider model uncertainties and external disturbances that affect the control performances of the robot manipulator.

Robust control strategies such as H infinity control theory have been suggested to achieve the good control performances despite the presence of model uncertainties and external disturbances. But the controller design and implementation based on the approach is complex [7, 8].

Therefore, some researchers have investigated the robust control techniques using PID controller with simple and easy-to-implement structure. These techniques mainly focus on how to tune the parameters of the PID controller so that the closed-loop system guarantees robustness and disturbance attenuation. In refs. [9–12], the authors discussed the controller parameter tuning by using system identification method. In ref. [13], PID controller tuning for a planar parallel robot is stated as a nonlinear optimization problem which is solved by differential evolution (DE) algorithm. However, these techniques consider only one or two factors such as the stability and tracking performance to tune the controller parameters.

Recently, active researches on the fractional order PID control have been conducted. The fractional order PID controller provides the improvement of stability and performance of the system although model uncertainties and external disturbances appear. A typical method of the fractional order PID controller tuning is Monje’s method [14, 15]. It has an advantage of tuning a controller so that the closed-loop system meets several conditions regarding robustness to plant uncertainties, load disturbances and high-frequency noise. In Monje’s method, tuning is reduced to a nonlinear optimization problem with nonlinear constraints. However, it is unlikely to obtain the accurate global minimum. Normally, the optimization problem was solved out by using the function FMINCON in MATLAB Optimization Toolbox, which employs the traditional gradient-based paradigm [2, 14, 15].

For a nonlinear optimization problem, however, the gradient-based approaches reveal the drawbacks of the premature convergence to local minimum, high sensitivity to selection of the initial condition and the failure of convergence when discontinuities exist. Therefore, the obtained solution may not be a global optimal solution. It implies that it is necessary to overcome these issues. This motivates the study of this paper.

To the best of the authors’ knowledge, the robust PID controller design approaches based on Monje’s method are mainly applied to process control [16, 17]. In this paper, to enhance the convergence in Monje’s method, we propose a method based on a hybrid algorithm of particle swarm optimization (PSO) and DE rather than FMINCON used in previous works. This PSO-DE hybrid optimization algorithm has advantage of global search capability which was shown in refs. [7, 8].

The main contribution of this paper is to design a robust PID controller that meets some requirements such as robustness, high-frequency noise attenuation and disturbance rejection by using PSO-DE and to apply it to the trajectory tracking control of Delta parallel robot. In this paper, the design problem of the robust PID controller for the parallel robot has been investigated in the following stages:

• The design specifications are set concretely.

• The design problem is referred to an optimization problem with several constraints regarding design specifications using Monje’s method.

• To obtain global solution, this optimization problem is solved using the PSO-DE.

• For comparison, another PID controller is also tuned using the method in ref. [2].

Simulation and Experimental results show that the designed robust PID controller has a robust performance and active disturbance attenuation when it is applied to tracking control of the Delta parallel robot. The rest of this paper is organized as follows. Section 2 shows the structure, the kinematic and dynamic models of the developed Delta parallel robot. Section 3 presents the control structure of robust PID controller with CTC strategy, design specifications and PSO-DE-based tuning algorithm of the robust PID controller. Simulation and experimental results are discussed in Section 4. Finally, conclusions are given in Section 5.

2. Modeling of the Delta parallel robot

The Delta parallel robot is a 3-degree-of-freedom robot, which consists of a fixed platform, a mobile platform and three identical RUU (Revolution, Universal, Universal) legs between the fixed platform and the mobile one. The revolute joints are actuated by fixed rotational actuators. Since actuators are placed in the fixed platform, the Delta parallel robot can achieve high speed and acceleration [2, 18].

Figure 1 shows prototype of the Delta parallel robot on which experiments were carried out. Figures 2 and 3 show physical parameters, a base coordinate system ( $O - xyz$ ) and chain coordinate systems ( ${O_i} -{x_i}{y_i}{z_i}, i = 1,2,3$ ) established in the Delta parallel robot. In addition, the physical parameters of the Delta parallel robot are listed in Table I.

Figure 1. Prototype of the Delta parallel robot.

Figure 2. The base coordinate systems and physical parameters.

Figure 3. The structure after simplification.

As shown in Fig. 2, the robot has a base coordinate system $O - xyz$ located at the center of the fixed platform, in which z-axis is in the reverse direction of the gravity, x-axis is directed to the upper leg of the first kinematic chain and y-axis is chosen to be the right-handed coordinated system. Over each kinematic chain, a chain coordinate system ${O_i} -{x_i}{y_i}{z_i}, i = 1,2,3$ ) is located at distance $R$ from $O - xyz$ and rotated an angle ${\phi _i}(0^\circ, 120^\circ, 240^\circ )$ , respectively.

For simplicity, upper legs and lower legs of the three kinematic chains are translated to the origin of the base coordinate system by radius $r$ of the mobile platform, as shown in Fig. 3.

2.1. Kinematic model

In general, the kinematic model of the robot consists of the inverse kinematic model and the direct kinematic model. The kinematic model of delta robot has been widely researched [2, 18–20].

Given the coordinates $P ={[{x_p},\,{y_p},\,{z_p}]^{\rm{T}}}$ of the mobile platform in the base Cartesian coordinate system, the inverse kinematic model computes actuated revolute joint angles $\theta ={[{\theta _1}, {\theta _2}, {\theta _3}]^{\rm{T}}}$ of upper legs. The inverse kinematic model can be summarized as follows [2]. Given the coordinates $P ={[{x_p},\,{y_p},\,{z_p}]^{\rm{T}}}$ of the mobile platform on the $O - xyz$ , its position in each chain coordinate system ${O_i} -{x_i}{y_i}{z_i}\,\,(i = 1,2,3)$ is converted using a rotation transformation matrix as follows:

(1) \begin{equation} \left [{\begin{array}{*{20}{c}}{{x_{{p_i}}}}\\[5pt] {{y_{{p_i}}}}\\[5pt] {{z_{{p_i}}}} \end{array}} \right ] ={\left [{\begin{array}{*{20}{c}}{\cos{\phi _i}}\;\;\;\;\;&{ - \sin{\phi _i}}\;\;\;\;\;&0\\[5pt] {\sin{\phi _i}}\;\;\;\;\;&{\cos{\phi _i}}\;\;\;\;\;&0\\[5pt] 0\;\;\;\;\;&0\;\;\;\;\;&1 \end{array}} \right ]^{ - 1}}\left [{\begin{array}{*{20}{c}}{{x_p}}\\[5pt] {{y_p}}\\[5pt] {{z_p}} \end{array}} \right ],\,\,\,({\phi _i} = 0,\frac{{2\pi }}{3},\frac{{4\pi }}{3};\,\,i = 1,2,3) \end{equation}

Since the upper leg link $L_1$ can rotate only on the plane ${x_i}{O_i}{z_i}$ , the joint $J_i$ is located on the circular trajectory with the center at $A_i$ and a $L_1$ radius.

(2) \begin{equation}{({x_i} - d)^2} + z_i^2 = L_1^2,\,\,\,\,(i = 1,2,3) \end{equation}

Moreover, the joint $J_i$ is also located on a spherical surface centered in $B_i$ with radius $L_2$ .

(3) \begin{equation}{({x_i} -{x_{{p_i}}})^2} +{({y_i} -{y_{{p_i}}})^2} +{({z_i} -{z_{{p_i}}})^2} = L_2^2,\,\,\,(i = 1,2,3) \end{equation}

Then, $J_i$ is located on the intersection between (2) and (3) in the following plane

(4) \begin{equation}{y_i} = 0,\,\,\,(i = 1,2,3). \end{equation}

From (2) to (4), the position of joint $J_i$ in each chain coordinate system can be stated as

(5) \begin{equation}{x_{{J_i}}} = \frac{{(d -{A_i}{B_i}) \pm \sqrt{(1 + A_i^2)L_1^2 -{{({A_i}d +{B_i})}^2}} }}{{1 + A_i^2}},\;{y_{{J_i}}} = 0,\;{z_{{J_i}}} ={A_i}{x_{{J_i}}} +{B_i}, \end{equation}

where

\begin{equation*}{A_i} = - \frac {{d - {x_{{p_i}}}}}{{{z_{{p_i}}}}}, {B_i} = - \frac {{x_{{p_i}}^2 + y_{{p_i}}^2 + z_{{p_i}}^2 + L_1^2 - L_2^2 - {d^2}}}{{2{z_{{p_i}}}}}\end{equation*}

Table I. Physical parameters for the Delta parallel robot.

Considering [21, 22], we have the revolute joint angle of the upper leg as below:

(6) \begin{equation}{\theta _i} = \arctan 2\left ({{z_{{J_i}}},\;{x_{{J_i}}} - d} \right ),i = 1,2,3. \end{equation}

In the Delta parallel robot, the direct kinematic model is the reverse process of the inverse kinematic model, which calculates the Cartesian position $P ={[{x_p},{y_p},{z_p}]^{\rm{T}}}$ of the mobile platform in the base coordinate system, given the revolute joint angles $\theta ={[{\theta _1},{\theta _2},{\theta _3}]^{\rm{T}}}$ of the upper legs.

The position of joint $J_i$ in each chain coordinate system is obtained as follows:

(7) \begin{equation} x_{{J_i}}^\prime = d +{L_1}\cos{\theta _i},\;y_{{J_i}}^\prime = 0,\;z_{{J_i}}^\prime = -{L_1}\sin{\theta _i}. \end{equation}

Then, the position of $J_i$ in base coordinate system is obtained using the matrix of revolution transformation as follows:

(8) \begin{equation} \left [{\begin{array}{*{20}{c}}{{x_{{J_i}}}}\\[5pt] {{y_{{J_i}}}}\\[5pt] {{z_{{J_i}}}} \end{array}} \right ] ={\left [{\begin{array}{*{20}{c}}{\cos{\phi _i}}\;\;\;\;\; &{ - \sin{\phi _i}}\;\;\;\;\; &0\\[5pt] {\sin{\phi _i}}\;\;\;\;\; &{\cos{\phi _i}}\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &0\;\;\;\;\; &1 \end{array}} \right ]^{ - 1}}\left [{\begin{array}{*{20}{c}}{x_{{J_i}}^\prime }\\[5pt] {y_{{J_i}}^\prime }\\[5pt] {z_{{J_i}}^\prime } \end{array}} \right ],\;({\phi _i} = 0,\;\frac{{2\pi }}{3},\;\frac{{4\pi }}{3};\;i = 1,2,3). \end{equation}

Moreover, the joint $J_i$ is also located on a spherical surface centered in $p$ with radius $L_2$ .

(9) \begin{equation}{({x_{{J_i}}} -{x_p})^2} +{({y_{{J_i}}} -{y_p})^2} +{({z_{{J_i}}} -{z_p})^2} = L_2^2,\;(i = 1,2,3). \end{equation}

From (7) to (9),

(10) \begin{equation} \left \{{\begin{array}{*{20}{l}}{{{({x_{{J_1}}} -{x_p})}^2} + y_p^2 +{{({z_{{J_1}}} -{z_p})}^2} = L_2^2}\\[5pt] {({x_{{J_2}}} -{x_{{J_1}}}){x_p} +{y_{{J_2}}}{y_p} + ({z_{{J_2}}} -{z_{{J_1}}}){z_p} ={\dfrac{1}{2}}({\rho _2} -{\rho _1})}\\[12pt] {({x_{{J_3}}} -{x_{{J_1}}}){x_p} +{y_{{J_3}}}{y_p} + ({z_{{J_3}}} -{z_{{J_1}}}){z_p} ={\dfrac{1}{2}}({\rho _3} -{\rho _1})} \end{array}} \right . \end{equation}

where ${\rho _i} = x_{{J_i}}^2 + y_{{J_i}}^2 + z_{{J_i}}^2,\;(i = 1,2,3).$

From last two equations of (10),

(11) \begin{equation}{x_p} ={K_1}{z_p} +{H_1},\;{y_p} ={K_2}{z_p} +{H_2}, \end{equation}

where

${K_1} = \dfrac{{{y_{{J_3}}}({z_{{J_2}}} -{z_{{J_1}}}) -{y_{{J_2}}}({z_{{J_3}}} -{z_{{J_1}}})}}{{({x_{{J_3}}} -{x_{{J_1}}}){y_{{J_2}}} - ({x_{{J_2}}} -{x_{{J_1}}}){y_{{J_3}}}}},\;{H_1} = \dfrac{1}{2} \cdot \dfrac{{{y_{{J_2}}}({\rho _3} -{\rho _1}) -{y_{{J_3}}}({\rho _2} -{\rho _1})}}{{({x_{{J_3}}} -{x_{{J_1}}}){y_{{J_2}}} - ({x_{{J_2}}} -{x_{{J_1}}}){y_{{J_3}}}}}$

${K_2} = \dfrac{{({x_{{J_2}}} -{x_{{J_1}}})({z_{{J_3}}} -{z_{{J_1}}}) - ({x_{{J_3}}} -{x_{{J_1}}})({z_{{J_2}}} -{z_{{J_1}}})}}{{({x_{{J_3}}} -{x_{{J_1}}}){y_{{J_2}}} - ({x_{{J_2}}} -{x_{{J_1}}}){y_{{J_3}}}}},\;{H_2} = \dfrac{1}{2} \cdot \dfrac{{({x_{{J_3}}} -{x_{{J_1}}})({\rho _2} -{\rho _1}) - ({x_{{J_2}}} -{x_{{J_1}}})({\rho _3} -{\rho _1})}}{{({x_{{J_3}}} -{x_{{J_1}}}){y_{{J_2}}} - ({x_{{J_2}}} -{x_{{J_1}}}){y_{{J_3}}}}}$

From first equation of (10) and (11),

(12) \begin{equation} (1 + K_1^2 + K_2^2)z_p^2 + 2({K_1}{H_1} +{K_2}{H_2} -{K_1}{x_{{J_1}}} -{z_{{J_1}}}){z_p} +{({H_1} -{x_{{J_1}}})^2} + H_2^2 + z_{{J_1}}^2 - L_2^2 = 0. \end{equation}

Then, the position of $P$ is obtained as follows:

(13) \begin{equation}{x_p} ={K_1}{z_p} +{H_1},\;{y_p} ={K_2}{z_p} +{H_2},\;{z_p} = \frac{{ - b + \sqrt{{b^2} - 4ac} }}{{2a}}, \end{equation}

where

\begin{equation*}a = (1 + K_1^2 + K_2^2),\;b = 2({K_1}{H_1} + {K_2}{H_2} - {K_1}{x_{{J_1}}} - {z_{{J_1}}}),\;c = {({H_1} - {x_{{J_1}}})^2} + H_2^2 + z_{{J_1}}^2 - L_2^2.\end{equation*}

2.2. Dynamic model

The dynamic equation can be expressed as follows:

(14) \begin{equation} M(\theta )\ddot \theta + C(\theta,\dot \theta )\dot \theta + N(\theta ) = \tau, \end{equation}

where $M(\theta )$ is the inertia matrix, $C(\theta,\dot \theta )$ is the Coriolis matrix, $N(\theta )$ is the gravity vector, and $\tau ={[{\tau _1},{\tau _2},{\tau _3}]^{\rm{T}}}$ is the torque vector applied at the upper leg joints.

The inertia matrix $M(\theta )$ in (14) is defined as

(15) \begin{equation} M(\theta ) = \left [{\begin{array}{*{20}{c}}{{m_{11}}}\;\;\;\;\;&{{m_{12}}}\;\;\;\;\;&{{m_{13}}}\\[5pt] {{m_{21}}}\;\;\;\;\;&{{m_{22}}}\;\;\;\;\;&{{m_{23}}}\\[5pt] {{m_{31}}}\;\;\;\;\;&{{m_{32}}}\;\;\;\;\;&{{m_{33}}} \end{array}} \right ] ={I_{bt}} + ({m_3} +{m_4} + \frac{1}{3}{m_2}){J^T}J, \end{equation}

where $I_{bt}$ is the inertia matrix of the driving axis and is defined as follows:

\begin{equation*}{I_{bt}} = \left [ {\begin {array}{*{20}{c}} {{I_{bti}}}\;\;\;\;\; &0\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &{{I_{bti}}}\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &0\;\;\;\;\; &{{I_{bti}}} \end {array}} \right ],\,\,\,\,{I_{bti}} = {I_m} + L_1^2\left ( {\frac {{{m_1}}}{3} + {m_{1j}} + \frac {2}{3}{m_2}} \right ).\end{equation*}

And $J$ in (15) is Jacobian matrix defined as follows:

\begin{equation*}J = - {\left [ {\begin {array}{*{20}{c}} {s_1^T}\\[5pt] {s_2^T}\\[5pt] {s_3^T} \end {array}} \right ]^{ - 1}}\left [ {\begin {array}{*{20}{c}} {s_1^T{b_1}}\;\;\;\;\; &0\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &{s_2^T{b_2}}\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &0\;\;\;\;\; &{s_3^T{b_3}} \end {array}} \right ], {s_i} = \left [ {\begin {array}{*{20}{c}} x\\[5pt] y\\[5pt] z \end {array}} \right ] - {T_i}\left ( {\left [ {\begin {array}{*{20}{c}} d\\[5pt] 0\\[5pt] 0 \end {array}} \right ] + \left [ {\begin {array}{*{20}{c}} {{L_1}\cos {\theta _i}}\\[5pt] 0\\[5pt] { - {L_1}\sin {\theta _i}} \end {array}} \right ]} \right ),\end{equation*}

\begin{equation*}{b_i} = {T_i}\!\left [ {\begin {array}{*{20}{c}} {{L_1}\sin {\theta _i}}\\[5pt] 0\\[5pt] {{L_1}\cos {\theta _i}} \end {array}} \right ],\, {T_i} = \left [ {\begin {array}{*{20}{c}} {\cos {\phi _i}}\;\;\;\;\; &{ - \sin {\phi _i}}\;\;\;\;\; &0\\[5pt] {\sin {\phi _i}}\;\;\;\;\; &{\cos {\phi _i}}\;\;\;\;\; &0\\[5pt] 0\;\;\;\;\; &0\;\;\;\;\; &1 \end {array}} \right ],\,\,\,({\phi _i} = 0,\frac {{2\pi }}{3},\frac {{4\pi }}{3};\,\,i = 1,2,3).\end{equation*}

$X ={[x,\,y,\,z]^T} = f(\theta )$ : Position of the mobile platform center

The matrix $C(\theta,\dot \theta )$ in (14) is calculated as follows:

(16) \begin{equation} C(\theta,\dot \theta ) = -{J^T}({m_3} +{m_4} + \frac{1}{3}{m_2}){\left [{\begin{array}{*{20}{c}}{s_1^T}\\[5pt] {s_2^T}\\[5pt] {s_3^T} \end{array}} \right ]^{ - 1}}\left ({\left [{\begin{array}{*{20}{c}}{\dot s_1^T}\\[5pt] {\dot s_2^T}\\[5pt] {\dot s_3^T} \end{array}} \right ]J +{M_{sb}}} \right ), \end{equation}

\begin{align*}{M_{sb}} & = \left [ {\begin {array}{*{20}{c}} {\dot s_1^T{b_1} + s_1^T{{\dot b}_1}}\;\;\;\;&0\;\;\;\;&0\\[5pt] 0\;\;\;\;&{\dot s_2^T{b_2} + s_2^T{{\dot b}_2}}\;\;\;\;&0\\[5pt] 0\;\;\;\;&0\;\;\;\;&{\dot s_3^T{b_3} + s_3^T{{\dot b}_3}} \end {array}} \right ], {\dot s_i} = \dot X + {b_i}{\dot \theta _i} = J\dot \theta + {b_i}{\dot \theta _i}, \\ & {\dot b_i} = {T_i}\!\left [ {\begin {array}{*{20}{c}} {{L_1}\cos {\theta _i}}\\[5pt] 0\\[5pt] { - {L_1}\sin {\theta _i}} \end {array}} \right ]{\dot \theta _i},\;(i = 1,2,3).\end{align*}

In addition, $N(\theta )$ in (14) is calculated as follows:

(17) \begin{equation} N(\theta ) ={J^T}\!\left ({{m_3} +{m_4} + \frac{3}{2}{m_2}} \right )\left [{\begin{array}{*{20}{c}} 0\\[5pt] 0\\[5pt] g \end{array}} \right ] -{L_1}\!\left ({\frac{1}{2}{m_1} +{m_{1j}} + \frac{1}{2}{m_2}} \right )g\left [{\begin{array}{*{20}{c}}{\cos{\theta _1}}\\[5pt] {\cos{\theta _2}}\\[5pt] {\cos{\theta _3}} \end{array}} \right ] . \end{equation}

3. Robust PID control based on a PSO-DE hybrid optimization algorithm

3.1. Control structure for robust PID controller with CTC strategy

The CTC is a control strategy, which obtains a linearized and decoupled model of robot manipulator using a feedback linearization approach. Applying this strategy, the robot system can be controlled by the linear control techniques.

Consider the following linearization law:

(18) \begin{equation} \tau = \hat M(\theta )u + \hat C(\theta,\;\dot \theta )\dot \theta + \hat N(\theta ), \end{equation}

where $u$ is an input of new system and $\hat M(\theta )$ , $\hat C(\theta,\;\dot \theta )$ , $\hat N(\theta )$ are estimations of $M(\theta )$ , $C(\theta,\;\dot \theta )$ , $N(\theta )$ , respectively.

Assuming that Eq. (18) is exact enough to describe a dynamic model of the robotic system

\begin{equation*}M(\theta ) = \hat M(\theta ),\;C(\theta,\dot \theta ) = \hat C(\theta,\dot \theta ),\;N(\theta ) = \hat N(\theta ).\end{equation*}

Then, we get the Eq. (19) from (14) and (18).

(19) \begin{equation} \hat M(\theta )u = M(\theta )\ddot \theta . \end{equation}

If the inertia matrix $M(\theta )$ is invertible, the dynamic model (14) is linearized and decoupled.

(20) \begin{equation} \ddot \theta (t) = u(t). \end{equation}

The obtained system is a linearized and decoupled system for joint variables $\theta (t)$ .

In practice, however, model uncertainties and external disturbances should not be ignored, which affect the control performances of robotic system. These influences and errors can be attenuated by the feedback control for linear system (20).

For the linearized system (20), the following PID control strategy is proposed:

(21) \begin{equation} u(t) ={\ddot \theta _d}(t) +{k_p}e(t) +{k_i}\int{e(t)dt} +{k_d}\dot e(t), \end{equation}

where $e(t)=\theta _d(t)-\theta (t)$ is the joint error, $\dot{e}(t)$ is the joint velocity error, ${\ddot{\theta }}_d(t)$ is the desired joint acceleration, and $k_p$ , $k_i$ and $k_d$ are the controller parameters to be determined.

If we replace the ideal derivative controller in (21) with the real derivative controller, then the resulting PID controller can be expressed as follows:

(22) \begin{equation} C(s) ={K_p}\!\left ({1 + \frac{1}{{{T_i}s}} +{T_d}s} \right )\frac{1}{{{T_f}s + 1}}, \end{equation}

where $K_p$ is the proportional coefficient, $T_i$ and $T_d$ are the integral and derivative time constants, respectively. In addition, an output first-order filter with the time constant $T_f$ makes the PID controller proper, and it filters the high-frequency noise. $T_f$ is determined as in ref. [23].

(23) \begin{equation}{T_f} = \min \!\left \{{{T_i}/10,\;{T_d}/10} \right \}. \end{equation}

Figure 4 shows the control structure of the robust PID controller with CTC strategy for the Delta parallel robot. Most controllers for Delta parallel robot use joint sensors that allow us to easily measure joint angles. The trajectory is planned in Cartesian space, and the trajectory in joint space is obtained using inverse kinematic model. The controller uses the trajectory in joint space and measured joint position. Then, the position in Cartesian space is obtained using direct kinematic model.

Figure 4. Structure for the robust PID controller with CTC strategy for Delta parallel robot.

In Fig 4, $\theta$ is measured joint position and $\theta _d$ is desired joint position. The error in joint position is reduced by individual PID controller with first-order filter. And the torque vector $\tau$ is calculated using (18). This calculation is performed in low-level controller equipped with ARM processor.

3.2. Design specifications for the robust PID controller

In this section, we discuss how to design the robust PID controller with CTC strategy to simultaneously meet several specifications regarding frequency domain.

Robustness to gain variations of the plant is one of the most common specifications that attribute performance of the robust PID controller [24]. The robustness forces the phase of the open-loop system to be flat at the gain crossover frequency $\omega _{cg}$ and to be almost constant within an interval around $\omega _{cg}$ . Then, the system will be more robust to the gain variations and the overshoot of the transient response is almost constant within the range of gain variations.

Therefore, the robust PID controller design can be referred to solve an optimization problem, which minimizes the following cost function as follows for each possible solution [14].

(24) \begin{equation} J({K_p},{T_i},{T_d}) = \int \limits _{{\omega _i}}^{{\omega _f}}{e_{FD}^2(\omega )d\omega } = \int \limits _{{\omega _i}}^{{\omega _f}}{{{(\varphi (\omega ) - ({\varphi _m} - \pi ))}^2}d\omega }, \end{equation}

where error $e_{FD}(\omega )$ is calculated in a predefined interval between $\omega _i$ and $\omega _f$ around $\omega _{cg}$ . $\varphi (\omega )$ is the phase of the open-loop system at frequency $\omega$ and $\varphi _m$ is the phase margin.

The advantage of this approach is that the frequency band [ $\omega _i$ , $\omega _f$ ] in which the closed-loop system is robust to the gain variations can be selected according to user requirements. More specifications can be included in the cost function (24) with the following constraints [14, 15].

• The gain crossover frequency $\omega _{cg}$ must be within frequency interval [ $\omega _{cgl}$ , $\omega _{cgh}$ ], and the phase margin must be at least $\varphi _m$ . Then, these constraints can be described as follows:

  (25) \begin{equation}{\omega _{cgl}} \le{\omega _{cg}} \le{\omega _{cgh}}, \end{equation}
  (26) \begin{equation} \arg (C(j{\omega _{cg}})G(j{\omega _{cg}})) \ge - \pi +{\varphi _m}, \end{equation}
  where $C(j\omega _{cg})$ and $G(j\omega _{cg})$ are transfer functions of controller (22) and plant (20), respectively.

• To attenuate the high-frequency noise, the complementary sensitivity function needs to fulfill the following condition:

  (27) \begin{equation}{\left |{T(j\omega ) = \frac{{C(j\omega )G(j\omega )}}{{1 + C(j\omega )G(j\omega )}}} \right |_{{\rm{dB}}}} \le{M_t}\,{\rm{dB,}}\,\;\forall \omega \ge{\omega _t}\,{\rm{rad/s, }}\,|T(j{\omega _t}){|_{{\rm{dB}}}} ={M_t}\,{\rm{dB}}, \end{equation}
  where $M_t$ means the desired noise attenuation for frequencies $\omega \geq \omega _t$ rad/s.

• To ensure a good output disturbance rejection, the sensitivity function $S(j\omega )$ is constrained by the following condition:

  (28) \begin{equation}{\left |{S(j\omega ) = \frac{1}{{1 + C(j\omega )G(j\omega )}}} \right |_{{\rm{dB}}}} \le{M_s}\,{\rm{dB,}}\;\forall \omega \le{\omega _s}{\rm{rad/s, }}\,|S(j{\omega _s}){|_{{\rm{dB}}}} ={M_s}\,{\rm{dB}}, \end{equation}
  where $M_s$ is the desired value of the sensitivity function for frequencies $\omega \le \omega _s$ .

Therefore, the design problem of the robust PID controller can be referred to a nonlinear optimization problem with nonlinear constraints as follows:

(29) \begin{equation} \min \,\,J({K_p},\,\,{T_i},\,\,{T_d}) = \int \limits _{{\omega _i}}^{{\omega _f}}{{{(\varphi (\omega ) - ({\varphi _m} - \pi ))}^2}d\omega }. \end{equation}

Subject to (25)-(28).

In refs. [1, 15], the FOPID controller design is referred to as the similar nonlinear optimization problem, and it was solved using the function FMINCON of MATLAB. For the optimization problem (29), however, it is hard to find the global optimal solution with traditional approximation optimization methods such as gradient method since the cost function and constraints are complex nonlinear function on the design parameters ( $K_p$ , $T_i$ , $T_d$ ).

Therefore, in the next section, the nonlinear optimization problem with nonlinear constraints is to be solved effectively by a PSO-DE hybrid optimization algorithm, which is one of the swarm intelligent optimization techniques [7, 8].

3.3. PSO-DE-based tuning algorithm of the robust PID controller

The feature of the PSO algorithm is to achieve the optimal solution by sharing of information among particles in the swarm. It is assumed that the swarm consists of NP particles in $D$ -dimensional search space (i.e., solution space). Each particle is a $D$ -dimensional vector and its position represents a candidate solution in the solution space. Since it needs to identify three parameters in consideration of (29) in our case, the dimension of the particle is $D=3$ .

If $j$ is the current time step (generation), then the position of the i-th particle in the swarm can be represented as $\eta _j^i=(\eta _j^{i,1},\ \eta _j^{i,2},\ \ldots \, ,\eta _j^{i,D})$ , and it corresponds to

\begin{equation*}\eta _j^i=[K_p^i,\ T_i^i,\ T_d^i]\end{equation*}

considering (22) and (24)-(29), where the subscript j is omitted for description brevity. The velocity of each particle is denoted by $V_j^i=(V_j^{i,1},\ V_j^{i,2},\ \ldots \, , V_j^{i,D})$ . The best current position (i.e., the best personal position of the particle) of the i-th particle is denoted by $P_j^i=(p_j^{i,1},\ p_j^{i,2},\ \ldots \, , p_j^{i,D})$ , and it is known as the cognitive experience. The best particle in the swarm (that is, in minimization problems, the particle with the smallest cost function value) is indicated by $P_j^{best}=(p_j^{best,1},\ p_j^{best,2},\ \ldots \, , p_j^{best,D})$ , and it is known as the social experience.

We discuss the PSO version proposed by Clerc and Kennedy [25], which uses a new parameter $\chi$ , known as the constriction factor that improves the convergence. For each time step $j$ (generation), the velocities and the positions of each particle in the swarm are represented as follows:

(30) \begin{equation} V_{j+1}^i=\chi (V_j^i+c_1r_1(P_j^i-\eta _j^i)+c_2r_2(P_j^{best}-\eta _j^i)), \end{equation}

(31) \begin{equation} \eta _{j+1}^i=\eta _j^i+V_{j+1}^i, \end{equation}

for $i=1,\ 2,\ \ldots \,, NP$ , where $\chi$ is the above-mentioned constriction factor parameter, $c_1$ and $c_2$ are positive constants which are referred to as cognitive and social parameters, respectively, $r_1$ and $r_2$ denote random numbers that are uniformly distributed in [0, 1] and are randomly chosen in each generation. This scheme is typically utilized for the constant $\chi =0.72984$ and $c_1=c_2=2.05$ [25, 26].

The DE algorithm is a population-based stochastic parallel direct search method, and it requires few control parameters. However, it has better convergence behavior than other well-known algorithms such as Evolution Algorithm [26]. Like PSO, the population consists of $NP$ individuals that are the $D$ -dimensional vector and is initialized at random in the search space. The individuals are evolved over successive iterations (generations) to find the optimal solution of the cost function in the solution space. At each generation, the evolution of individuals is conducted through mutation, crossover and selection operations [27].

Nowadays, numerous PSO-DE hybrid approaches have been developed to improve the accuracy of the optimal solution and the convergence behavior of PSO and DE algorithms. The evolutionary algorithms are known to be still efficient, but it requires a great computational cost due to a large number of function evaluations. The incorporation of DE algorithm in each generation of PSO algorithm may cause an increase of the function evaluations. So, we adopt a PSO-DE hybrid algorithm based on the evolution of cognitive and social experiences in order to suggest an approach of the parameter tuning of the robust PID controller, since it costs smaller amount of computation along with better convergence behavior than other hybrid versions [7, 8].

In the tuning algorithm, we define a set of the best personal experience (i.e., the cognitive experiences) of each particle at the j-th time step (generation) as $S_j^p=(P_j^1,\ P_j^2,\ \ldots \, , P_j^{NP})$ . After each time step of PSO algorithm according to (30) and (31), DE procedure of one generation is applied to the elements $P_j^i,\ i=1,\ \ldots \,, NP$ , in the $S_j^p$ set (i.e., the cognitive experience of each particle is evolved). It implies that the best personal positions $P_j^i,\ i=1,\ \ldots \, , NP$ correspond to the individuals $x_j^i,\ i=1,\ \ldots \, , NP$ , of the DE algorithm. The DE mutation operation uses individuals from the $S_j^p$ set. In particular, three evolution steps of DE algorithm (mutation, crossover and selection) are applied only to the ‘promising’ cognitive experiences which have been changed (i.e., improved) during the previous PSO step. From the viewpoint of the computational cost, this turns out to be more effective than the procedure of evolving all particles in the swarm.

The algorithmic scheme of the proposed approach to tune the parameters of the robust PID controller for the Delta parallel robot is as follows:

1: Individuals in the swarm (i.e., the controller parameters $\eta _{j0}^i=[K_{p0}^i,\ T_{i0}^i,\ T_{d0}^i], i=1,\ \ldots \, , NP)$ are randomly initialized in the given range. The population size NP should be chosen in the range from 5 to 10 times $D$ . In order to improve the global convergence, it is set to $NP=50$ .

2: For each particle $\eta _{j0}^i=[K_{p0}^i,\ T_{i0}^i,\ T_{d0}^i] (i=1,\ \ldots \, , NP)$ , if it agrees with inequalities of the stability and constraints (25)-(28), initialize the particle again.

3: For each particle $\eta _{j0}^i=[K_{p0}^i,\ T_{i0}^i,\ T_{d0}^i] (i=1,\ \ldots \, , NP)$ , compute the cost function value $J(\eta _0^i)$ using (29) and let the best personal position $P_0^i$ be the value. Then, let the social experience $P_0^{best}$ be the minimum of $P_0^i,\;i = 1,\,2,\,\ldots \, , NP$ .

4:for each time step j do

5:for each particle i in the swarm do

6: Update the particle $\eta _j^i=[K_{pj}^i,\ T_{ij}^i,\ T_{dj}^i]$ (parameter vector of the controller) using (30) and (31). If the updated parameter vector $\eta _j^i$ disagrees with the inequalities constraints (25)-(28), then updating process should be repeated until it does.

7: Evaluate particle $\eta _j^i=[K_{pj}^i,\ T_{ij}^i,\ T_{dj}^i], i.e.$ , compute the cost function value for the current parameter vector $\eta _j^i=[K_{pj}^i,\ T_{ij}^i,\ T_{dj}^i]$ using (29).

8: If $J(\eta _j^i)\lt J(P_{j-1}^i)$ , then update the cognitive experience as $P_j^i=\eta _j^i$ , else $P_j^i=P_{j-1}^i$ (i.e., do not update value of the previous generation). Also, if $J(P_j^i)\lt J(P_{j-1}^{best})$ , then update the social experience as $P_j^{best}=P_j^i$ , else $P_j^{best}=P_{j-1}^{best}$ .

9: If $P_j^i$ has been changed then

10: /* Evolve $P_j^i$ utilizing one DE step in $S_j^p$ */

11: Mutate $P_j^i$ and generate the corresponding mutant vector $v_j^i$ through the mutation strategy.

12: Crossover the mutant vector $v_j^i$ and target vector $P_j^i$ , and generate the corresponding trial vector $u_j^i$ .

13: Evaluate the trial vector $u_j^i$ , that is, under the trial vector $u_j^i$ agrees with the inequalities constraints (25)-(28) (If it disagrees with the constraints, then mutation and crossover should be repeated until it does), compute the cost function value for $u_j^i$ using (29).

14: If $J(u_j^i)\lt J(P_j^i)$ , then update the best personal position $P_j^i$ . Also, if $\ J(P_j^i)\lt J(P_j^{best})$ , then update the best position $P_j^{best}$ .

15: end if

16: end for

17: end for

Table II. Parameter values of PSO-DE, PSO and DE algorithms.

Figure 5. Convergence behavior of the tuning algorithm based on PSO-DE.

4. Experimental results and analysis

4.1. Robust PID controller design for the Delta parallel robot

Considering the specifications of the Delta parallel robot and results in refs. [14, 23], the design specifications in the frequency domain mentioned in the inequalities constraints (25)-(28) of the optimization problem (29) are listed in (a)-(e):

1. (a) Minimum phase margin $\varphi _m$ is equal to $45^{\circ }$ .

2. (b) Gain crossover frequency $\omega _{cg}$ must be in the interval $[\omega _{cgl},\ \omega _{cgh}]=[50\,{\rm rad/s},150\,{\rm rad/s}]$ .

3. (c) An interval ( $\omega _f-\omega _i$ ) where the phase is intended to be flat about the gain crossover frequency is equal to $0.25\,{\rm dec}$ .

4. (d) High-frequency noise attenuation: $\omega _t=350\,{\rm rad/s},\ M_t=-15\,{\rm dB}$ that is, $\forall \omega \geq \omega _t,\ |T(j\omega ) |_{\rm dB}\le -15\,{\rm dB}$ .

5. (e) Disturbance rejection: $\omega _s=10\,{\rm rad/s},\ M_s=-15\,{\rm dB}, {\rm that\,is}, \forall \omega \le \omega _s,\ |S(j\omega ) |_{\rm dB}\le -15\,{\rm dB}$ .

In order to find the parameters $(K_p,\ T_i,\ T_d,\ T_f)$ of the robust PID controller, this problem is converted into an optimization problem (29). And the controller is designed by the algorithm mentioned in the previous section. All the parameters used in our PSO-DE-based tuning algorithm are listed in Table II. PSO algorithm and DE algorithm also use the same parameters in Table II.

Figure 6. Bode diagram of open-loop system ( $PM=54.9^{\circ }$ , ${\ \omega }_{cg}=76.6\,{\rm rad/s}$ ).

Figure 7. Frequency specification of the complementary sensitivity function $T(j\omega )$ ( $\omega _t=350\,{\rm rad/s},$ $|T(j\omega _t)|_{\rm dB}=-17.1\,{\rm dB}$ ).

Figure 8. Frequency specification of the sensitivity function $S(j\omega )$ ( $\omega _s=10\,{\rm rad/s},$ $|S(j\omega _s)|_{\rm dB}=-25.6\,{\rm dB}$ ).

Table III. Path point of control simulation trajectory for Delta robot.

Figure 9. Trajectories of x, y, z axes in the Cartesian space.

Figure 10. Desired trajectory $\theta _d(t)$ in joint space.

Figure 11. Tracking error in Cartesian space (Orange: method in ref. [2], blue: robust PID).

Figure 12. Tracking error in joint space (Orange: method in [2], blue: robust PID).

Figure 5 shows the convergence behaviors of the robust PID controller tuning algorithm for the Delta parallel robot based on the PSO-DE, DE and PSO. PSO-DE has better global search capability than PSO and DE.

After 50 generations, the parameters of the robust PID controller are obtained as follows:

(32) \begin{equation} K_p=1853.4,\ T_i=1073.7,\ T_d=0.0413,\ T_f=0.0041. \end{equation}

The bode diagram of the open-loop system of designed controller is shown in Fig. 6. Two diagrams demonstrate that the designed controller meets the design specifications (a), (b) and (c).

The frequency responses of the complementary sensitivity $T(j\omega )$ and the sensitivity function $S(j\omega )$ of the designed closed-loop system are shown in Figs. 7 and 8, respectively. As shown in these figures, the complementary sensitivity function meets the design specification (d) and the sensitivity function meets the design specification (e).

Consequently, we can make sure that the robust PID controller designed by the proposed approach satisfies the above-mentioned design specifications in frequency domain.

Table IV. Simulation results.

Table V. Real mass of the Delta robot.

Figure 13. Tracking error in Cartesian space in robustness analysis simulation (Orange: method in ref. [2], blue: robust PID).

Figure 14. Tracking error in joint space in robustness analysis simulation (Orange: method in ref. [2], blue: robust PID).

4.2. Simulation results

In order to analyze the control performance of the designed robust PID controller in the time domain, a trajectory planning was implemented using the modified trapezoidal acceleration curve method, in which acceleration and jerk is small for pick-and-place task [28, 29].

Table III shows the waypoints of the track for the control simulation of the Delta parallel robot. The desired trajectories on x, y, z axes and in joint space are shown in Figs. 9 and 10, respectively.

The desired trajectories on x, y, z axes are obtained from waypoints in Table III using modified trapezoidal acceleration curve method, and the trajectories in joint are obtained from trajectories on task space using inverse kinematic model.

The following waypoints are determined as considering pick-and-place tasks in workspace.

To analyze the performance of the robust PID controller designed in the previous section, the MATLAB/Simulink model is built using the physical parameters of the Delta parallel robot in Table I, the nonlinear dynamic model (14) and the control structure of Fig. 4.

Another robust PID controller with CTC [2] is implemented and used for comparative analysis. The design specifications of the closed-loop system are set as follows:

\begin{equation*}\zeta = 0.9,\,\,{t_s} = 0.3\,{\rm s},\,\,{\omega _n} = 14.2\,{\rm {rad/s}},\,\,\,\beta = 15.\end{equation*}

Then, the obtained parameters are

\begin{equation*}{K_p}=583.42, {T_i}=3012.1, {T_d}=40.5073, {T_f}=0.0069.\end{equation*}

Figures 11 and 12 show that the tracking errors of the controllers in this paper and [2] and Table IV show the RMS of tracking errors. As shown in these figures and table, the controller in this paper has a smaller tracking error for the desired trajectory.

To analyze the robustness of controllers, real mass of delta robot is setup in Table V, which is different from Table I. Moreover, we assume that random noises of [-0.1 deg, 0.1 deg] arise in feedback loop and disturbances are present in the torque driving each upper leg in interval [4 s, 8 s].

\begin{equation*}0.5\sin(10t),\ 0.5\sin(15t),\ 0.5\sin(20t).\end{equation*}

Figures 13 and 14 show tracking errors of controllers in this paper and in ref. [2] for robustness analysis. Table VI shows RMS of tracking errors. As shown in these figures and table, the controller proposed in this paper has better robustness and good disturbance attenuation.

Table VI. Simulation results for robustness.

Figure 15. Tracking error in Cartesian space when ${m_4}=100\,{\rm g}$ .

4.3. Experimental results

First, we built the Delta parallel robot to be controlled by the robust PID controller designed in Section 4.1. Then, the tracking performances for the desired trajectory in the Cartesian space and the joint space are shown in Figs. 15 and 16 when mass of load is 100 g. Figures 17 and 18 show tracking errors when mass of load is 300 g. Table VII shows RMS of tracking errors.

Figure 16. Tracking error in the joint space when ${m_4}=100\,{\rm g}$ .

Figure 17. Tracking error in Cartesian space when ${m_4}=300\,{\rm g}$ .

Figure 18. Tracking error in the joint space when ${m_4}=300\,{\rm g}$ .

In conclusion, the robust PID controller designed in this paper tracks the desired trajectory with a small error, and a slight tracking error occurs when the desired value changes. At the same time, our experimental results show a slight increase of the tracking error than that of simulation results. It follows from the fact that the control performance is affected by nonlinearities as backlash of reduction gears, which is connected with the driving motors. If the load mass is less than 100g, the position tracking error is less than ±3.5 mm.

5. Conclusions

In this paper, we propose an approach to tune an optimal robust PID controller with CTC strategy for trajectory tracking of a Delta parallel robot by adopting a PSO-DE hybrid optimization algorithm. Considering the nonlinear dynamic model of the Delta parallel robot, the design problem of a robust PID controller with CTC strategy that guarantees robustness and disturbance attenuation has been formulated as a nonlinear optimization problem with nonlinear constraints. To obtain a more accurate global optimum of the nonlinear optimization problem, an algorithm based on a PSO-DE hybrid optimization algorithm has been newly proposed. Simulation and experimental results carried out on the Delta parallel robot show that the designed robust PID controller has robust performance and active disturbance attenuation despite the presence of model uncertainties and external disturbances.

Table VII. Experiment results.

The proposed method needs much calculations, but it is solved by hardware that has good computing power. FOPID has 5 parameters, so FOPID has better performances such as robustness than PID. In future works, the FOPID and ADRC (active disturbance rejection control) are proposed and applicated in other parallel robots.