Multi-objective-optimization-based control parameters auto-tuning for aerial manipulators

Abstract

The aerial manipulator has recently attracted much research attention due to its wide applications such as aerial cleaning, aerial transportation, and aerial manipulation. It is important to design a reliable controller for the aerial manipulator to robustly perform aerial tasks with different settings. However, current controllers still employ manual parameters tuning methods, which is mostly limited to a specific setting like a fixed aerial manipulator configuration or an unchanged environment. In fact, there could be diverse configurations of aerial manipulators and uncertain environments in practice, which requires the manual tuning process to be frequently repeated. This repetition is easy to be unavailable due to its significant cost of time and expensive involvements of control-tuning experts. To solve these problems, a novel multi-objective-optimization-based control parameters auto-tuning method is proposed for the aerial manipulator. Based on a conventional proportional–integral–derivative control structure, an evolutionary-algorithm-based optimization is used to automatically find optimal proportional–integral–derivative control parameters to satisfy conflicting objectives such as minimizing the integrated time square error and the control rate. Simulation results prove that the proposed method can achieve better control performances like smaller overshoots and faster stabilization time than manual tuning methods.

Introduction

Thanks to the physical interactions with surrounding environments, the capacities of unmanned aerial vehicles (UAVs) have been extended from “passive” tasks such as inspection, ¹ remote sensing, ² and surveillance ³ to “active” tasks like aerial manipulation. ^{4 –20} Although many tools including grippers, ^4,5 cables, ^{6 –9} screw drivers, ¹⁰ and brushes ¹¹ have been used for the important physical interactions, attaching a robotic arm to a UAV, named an aerial manipulator, ^{12 –20} has lately been considered as the most efficient and promising way. The aerial manipulator does not restrict the attitude of the payload or the reachable space of the end-effector as a stationary gripper-based platform does. It is also able to directly regulate the movement of the payload, of which the cable-based platform is not capable.

It is important to design a reliable controller for the aerial manipulator to robustly and accurately finish tasks. Some such controllers have been reported for different tasks. For example, a Cartesian impedance controller ¹² was proposed for a quadrotor equipped with a three-degree-of-freedom (3-DOF) robotic arm to execute dexterous manipulation tasks. A novel hierarchical motion control scheme, ¹³ including a top layer of an inverse kinematics algorithm and a bottom layer of a motion control algorithm, was introduced for a quadrotor with a 5-DOF robotic arm to follow circular helix trajectories. A sliding mode controller ¹⁴ was designed to transport objects using a quadrotor with a 2-DOF robotic arm. For the Aerial Robotics Cooperative Assembly System (ARCAS) project, a multilayer control architecture, ¹⁵ including a novel battery movement compensation, arm static compensation, and external generalized forces estimation and compensation, was proposed for an octorotor with a 6-DOF robotic arm. To allow outdoor operations, a stable backstepping-based controller and an admittance controller were developed for an octorotor and its 7-DOF robotic arm, respectively. ¹⁶ To fly with unknown objects, an augmented passivity-based controller ¹⁷ was designed for a hexacopter with a 2-DOF robotic arm by online estimating unknown parameters of the payload. The attitude control of a hexacopter equipped with a 2-DOF robotic arm was achieved by a quaternion-based backstepping method. ¹⁸ A variable parameter integral backstepping controller ¹⁹ was presented for an unmanned helicopter equipped with multi-link arms to perform structure assembly and manipulation tasks. An advanced controller, including an attitude and a rate command augmentation system, a feedforward torque compensation system, and an $L_1$ adaptive augmentation law, was designed to simultaneously address challenges of maximizing workspace in constrained spaces and achieving stable flights when lifting payloads with unknown masses. ²⁰

However, these controllers employed manual tuning of control parameters, which may severely hinder their applications in practice. One reason is that the robot configuration or the environment is not always unchanged in practice. As a well-tuned control parameter set is generally tuned for a specific dynamic model case, the practical time-varying dynamics could easily make the current controller no longer suitable. To address this problem, the manual tuning needs to be frequently repeated; however, it requires a large amount of time like days and expert knowledge. In addition, practical applications often require the controller to satisfy multiple but often conflicting goals like both a good task accuracy and a short completion time. This is much more difficult for the manual tuning compared with satisfying only one single objective. Even if the previous requirements are all achievable and acceptable, the manual tuning could be still difficult or even impossible due to human’s own limit, named bounded rationality, ^21,22 in handling complex systems. For example, humans need to tune 27 parameters for a 6-DOF quadrotor equipped with a 3-DOFs robotic arm using a proportional–integral–derivative (PID) controller; however, the bounded rationality argues that humans begin to struggle when four variables are involved in problem solving and more than five variables to be solved at the same time nearly impossible.

One way to solve the manual tuning problems is automatically tuning parameters based on multi-objective optimization (MOO), and its effectiveness has been proven by other practical applications such as automatic voltage regulator system, ²³ chemical titration and neutralization plant, ²⁴ and linear brushless direct current (DC) motor. ²⁵ Although the optimization of gains control is not a new idea, it has never been studied for aerial manipulators. As a useful tool for practical aerial tasks, it is meaningful to formulize the specific MOO-based auto-tuning method for the aerial manipulator and validate its feasibility since the aerial manipulator is much more complex than previous systems. ^{23 –25} Thus, we propose a new MOO-based control parameters auto-tuning method for the complex aerial manipulator system with specific problem formulizations. It is tested in two conventional parallel PID controllers that are designed for motion controls of a quadrotor and a 3-DOF robotic arm, respectively. Although the PID controller is basic and simple, it is proved to be useful with good performances in practical applications. For each PID controller, an evolutionary algorithm-based optimization, Non-Dominated Sorting in Genetic Algorithms (NSGA)-II, is used to automatically find optimal control parameters, avoiding the time-consuming manual tuning process. Two conflicting objectives, reliable static–dynamic performance and smooth control, are considered as an example of possibly conflicting performance criteria. The main contributions are considered as (1) the work provides a new MOO-based automatic parameters tuning to stably control the highly nonlinear and complex aerial manipulator system with possible dynamic changes in many different situations. (2) This method can be easily transferred to similar aerial manipulators with different UAVs and robotic arms using any controllers. This work does not need structure modifications or additional mechanical designs, ¹⁵ as the problem is fully solved through the control algorithms. (3) This work investigates meanings of multiple equivalent optimal solutions, as well as the consequences of adding new objectives. This could be used as a starting point for users to better understand how to select an appropriate set of control parameters according to their own requirements.

Modeling

As the aerial manipulator modeling is based on a fundamental but common Euler–Lagrange method, we only explain necessary kinematic and dynamic models, which have been studied in previous works, ^{12 –14,16 –18} for the following controller design.

Kinematics for the quadrotor–arm system

With the aerial manipulator system consisting of a quadrotor and a robotic arm shown in Figure 1, $O_w$ , $O_b$ , and $O_i$ represent the world frame (inertial frame), body frame, and link $i$ frame, where $i=1, 2, 3$ denotes the link number. If the inertial frame is arbitrarily chosen while the body frame and link frames are located at the center of gravity (CoG) of their respective rigid bodies, a vector including all the generalized coordinate variables can be defined as

Figure 1.

Coordinate configuration of the aerial manipulator system.

$\mathit{q}={\left[{\mathit{p}}^{\text{T}} {\mathit{\Omega }}^{\text{T}} {\mathit{\eta }}^{\text{T}}\right]}^{\text{T}}$

1

where $\mathit{p}={\left[x y z\right]}^{\text{T}}$ represents the position of the CoG of quadrotor in the world frame, $\Omega ={\left[\varphi \theta \psi \right]}^{\text{T}}$ indicates the Euler angles of the quadrotor, and $\mathit{\eta }={\left[{\eta }_1 {\eta }_2 {\eta }_3\right]}^{\text{T}}$ means the joint angles of the 3-DOF manipulator, which are defined about the positive $y_b$ axis.

The relationship between the velocities in different frames is given as follows

$\dot{\mathit{p}}=\mathit{R}{\dot{\mathit{p}}}^b$

2

$\mathit{\omega }=\mathit{R}{\mathit{\omega }}^b$

3

$\mathit{\omega }=\mathit{T}\dot{\Omega }$

4

${\mathit{\omega }}^b=\mathit{R}\mathit{T}\dot{\Omega }=\mathit{Q}\dot{\Omega }$

5

where $\dot{\mathit{p}}$ and ${\dot{\mathit{p}}}^b$ represent the translational velocity of the CoG of the quadrotor in the world and body-fixed frame, respectively. They are related by the rotation matrix $\mathit{R}\in SO\left(3\right)$ . $\mathit{\omega }=\left[{\omega }_x {\omega }_y {\omega }_z\right]$ and ${\mathit{\omega }}^b=\left[{\mathit{\omega }}_x^b {\mathit{\omega }}_y^b {\mathit{\omega }}_z^b\right]$ mean the angle velocity of the quadrotor in the world and body-fixed frame, respectively. They are also related by $\mathit{R}$ . Also, $\dot{\mathrm{\Omega }}$ can be mapped to $\mathit{\omega }$ by the transformation matrix $\mathit{T}$ . $\mathit{Q}=\mathit{R}\mathit{T}$ is the relationship matrix between ${\mathit{\omega }}^b$ and $\dot{\Omega }$ .

Let ${\mathit{p}}_i^b$ be the position of the CoG of link $i=1, 2, 3$ in the body-fixed frame $O_b$ . Then, ${\mathit{p}}_i$ , which is the position of the CoG of the link $i$ in the world frame $O_w$ , is related to ${\mathit{p}}_i^b$ by ${\mathit{p}}_i=\mathit{p}+\mathit{R}{\mathit{p}}_i^b$ . Also, if defining two Jacobian matrices ${\mathit{J}}_t\in {ℝ}^{3\times 3}$ and ${\mathit{J}}_r\in {ℝ}^{3\times 3}$ , which are used to relate the translational and angular velocity of each manipulator link with $\dot{\mathit{\eta }}$ by ${\dot{\mathit{p}}}_i^b={\mathit{J}}_t\dot{\mathit{\eta }}$ and ${\mathit{\omega }}_i^b={\mathit{J}}_r\dot{\mathit{\eta }}$ , then we can get the following equations

${\dot{\mathit{p}}}_i=\dot{\mathit{p}}+\dot{\mathit{R}}{\mathit{p}}_i^b+\mathit{R}{\dot{\mathit{p}}}_i^b$

6

${\dot{\mathit{p}}}_i=\dot{\mathit{p}}+{\widehat{\omega }}_b\mathit{R}{\mathit{p}}_i^b+\mathit{R}{\mathit{J}}_t\dot{\mathit{\eta }}$

7

${\mathit{\omega }}_i=\mathit{\omega }+\mathit{R}{\mathit{J}}_r\dot{\mathit{\eta }}$

8

For simplicity, equations (2), (4), (7), and (8) can be rewritten in the following matrix form ¹⁴

$\dot{\mathit{p}}=\left[{\mathit{I}}_{3\times 3} {\mathbf{0}}_{3\times 3} {\mathbf{0}}_{3\times 3}\right]\dot{\mathit{q}}={\mathit{M}}_{t,b}\dot{\mathit{q}}$

9

$\mathit{\omega }=\left[{\mathbf{0}}_{3\times 3} \mathit{T} {\mathbf{0}}_{3\times 3}\right]\dot{\mathit{q}}={\mathit{M}}_{r,b}\dot{\mathit{q}}$

10

${\dot{\mathit{p}}}_i=\left[{\mathit{I}}_{3\times 3}-{\left(\mathit{R}{\mathit{p}}_i^b\right)}^{\wedge }\mathit{T} \mathit{R}{\mathit{J}}_{t,i}\right]\dot{\mathit{q}}={\mathit{M}}_{t,i}\dot{\mathit{q}}$

11

${\mathit{\omega }}_i=\left[{\mathbf{0}}_{3\times 3} \mathit{T} \mathit{R}{\mathit{J}}_{r,i}\right]\dot{\mathit{q}}={\mathit{M}}_{r,i}\dot{\mathit{q}}$

12

where the subscript $\cdot \times \cdot$ represents the size of $\mathit{I}$ and $0$ . $\wedge$ is the operator that converts a vector into a skew-symmetric matrix. ${\mathit{J}}_{t,i}$ and ${\mathit{J}}_{r,i}$ are previously defined Jacobian matrices for link i. Through equations (9) to (12), $\dot{\mathit{q}}$ can be easily mapped into the translational and angular velocities of the quadrotor and each link in the inertial coordinate frame.

Dynamics for the quadrotor–arm system

Given the Euler–Lagrange method provides an easy way to derive the system dynamics, the following Lagrange-D’ Alembert equation is used

$\frac{\text{d}}{\text{d}t}\frac{\partial \mathcal{L}}{\partial \dot{\mathit{q}}}-\frac{\partial \mathcal{L}}{\partial \mathit{q}}=\mathit{u}+{\mathit{u}}_{\text{ext}}$

13

$\mathcal{L}=\mathcal{K}-\mathcal{U}$

14

where $\mathcal{K}$ and $\mathcal{U}$ are the total kinetic and potential energy of the combined system, $\mathit{u}$ represents the generalized force, and ${\mathit{u}}_{ext}$ indicates external disturbance applied to the system, such as wind and the interaction with an external object.

The total kinetic energy is contributed by the quadrotor and each individual link, so the total kinetic energy $\mathcal{K}$ is expressed as follows

$\mathcal{K}=\frac{1}{2}{\dot{\mathit{p}}}^{\text{T}}{m}_b\dot{\mathit{p}}+\frac{1}{2}{\dot{\Omega }}^{\text{T}}{\mathit{T}}^{\text{T}}\mathit{R}{\mathit{I}}_b{\mathit{R}}^{\text{T}}\mathit{T}\dot{\Omega }+{\displaystyle \sum }_{i=1}^3\frac{1}{2}{\dot{\mathit{p}}}_i^{\text{T}}{m}_i{\dot{\mathit{p}}}_i+\frac{1}{2}{\mathit{\omega }}_i^{\text{T}}\left(\mathit{R}{\mathit{R}}_i\right){\mathit{I}}_i{\left(\mathit{R}{\mathit{R}}_i\right)}^{\text{T}}{\mathit{\omega }}_i$

15

where m is the mass and $\mathit{I}$ is the inertia matrix. The subscripts b and i indicate the corresponding values with respect to $O_b$ and $O_i$ , respectively. ${\mathit{R}}_i$ is the coordinate transformation matrix from link i of the manipulator to the quadrotor body-fixed coordinates.

Similarly, the total potential energy is contributed by the potential energy of the quadrotor and each link, which is described by the equation below

$\mathcal{U}=m_{b}g{e}_3^{\text{T}}\mathit{p}+{\displaystyle \sum }_{i=1}^3{m}_{i}g{e}_3^{\text{T}}\left(\mathit{p}+\mathit{R}{\mathit{p}}_i^b\right)$

16

where the first and last terms are the potential energies of the quadrotor and link i, respectively. $e_3$ is the unit vector ${\left[0 0 1\right]}^{\text{T}}$ , and g is the gravity coefficient.

By computing equation (13) and considering the Christoffel symbols of the first type, ²⁶ the dynamics equation which includes all components as one system can be derived as the following

$\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}+\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}+\mathit{G}\left(\mathit{q}\right)=\mathit{u}+{\mathit{u}}_{ext}$

17

where $\mathit{M}\left(\mathit{q}\right)\in {ℝ}^{9\times 9}$ is the inertia matrix which is positive definite and symmetric. $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)$ is the Coriolis matrix, and $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}$ represents the Coriolis and Centrifugal forces. In addition, $\mathit{G}\left(\mathit{q}\right)$ includes gravity effects at each joint. ²⁷ $\mathit{u}$ is the vector of inputs, that is, $\mathit{u}={\left[{\mathit{u}}_f {\mathit{u}}_{\mu } {\mathit{u}}_{\tau }\right]}^{\text{T}}={\left[\mathit{R}{\mathit{f}}_b^b {\mathit{R}}^T\mathit{T}{\mathit{\mu }}_b^b \mathit{\tau }\right]}^{\text{T}}$ , $\mathit{\tau }\in {ℝ}^3$ is the vector of the manipulator joint torques, while ${\mathit{f}}_b^b\in {ℝ}^3$ and ${\mathit{\mu }}_b^b\in {ℝ}^3$ are the forces and torques generated by the four motors of the quadrotor, expressed in frame $O_b$ .

The input vectors ${\mathit{f}}_b^b$ and ${\mathit{\mu }}_b^b$ can be expressed as ${\mathit{f}}_b^b={\left[0 0 f_z\right]}^{\text{T}}$ and ${\mathit{\mu }}_b^b={\left[{\mu }_{\varphi } {\mu }_{\theta } {\mu }_{\psi }\right]}^{\text{T}}$ , where $f_z$ is the total thrust applied by the rotors along z-axis. Both $f_z$ and ${\mathit{\mu }}_b^b$ are related to the quadrotor’s four motor actuation forces $\mathit{f}$ via the following relation ²⁸

$\left[\begin{array}{c}{f}_z\\ {\mathit{\mu }}_b^b\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& l& 0& -l\\ -l& 0& l& 0\\ c& -c& c& -c\end{array}\right]\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\end{array}\right]=\Xi \mathit{f}$

18

where $l>0$ is the distance from each motor to the vehicle center of mass, $c={\gamma }_d/{\gamma }_t$ , and ${\gamma }_d$ , ${\gamma }_t$ are the drag and thrust coefficient, respectively, and $f_i$ is the force generated by rotor i.

The matrices $\mathit{M}\left(\mathit{q}\right)$ , $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)$ , $\mathit{G}\left(\mathit{q}\right)$ in equation (17) can be expressed by following equations. The derivation details can be found in previous works. ^12,13

$\begin{array}{c}\mathit{M}\left(\mathit{q}\right)={\mathit{M}}_{t,b}^{\text{T}}{m}_b{\mathit{M}}_{t,b}+{\mathit{M}}_{r,b}^{\text{T}}\mathit{R}{\mathit{I}}_b{\mathit{R}}^{\text{T}}{\mathit{M}}_{r,b}\\ +{\displaystyle \sum }_{i=1}^3{\mathit{M}}_{t,i}^{\text{T}}{m}_i{\mathit{M}}_{t,i}+{\mathit{M}}_{r,i}^{\text{T}}\left(\mathit{R}{\mathit{R}}_i\right){\mathit{I}}_i{\left(\mathit{R}{\mathit{R}}_i\right)}^{\text{T}}{\mathit{M}}_{r,i}\end{array}$

19

$c_{kj}={\displaystyle \sum }_{i=1}^9\frac{1}{2}\left\{\frac{\partial m_{kj}}{\partial {\mathit{q}}_i}+\frac{\partial m_{ki}}{\partial {\mathit{q}}_j}-\frac{\partial m_{ij}}{\partial {\mathit{q}}_k}\right\}$

20

$\mathit{G}\left(\mathit{q}\right)=\frac{\partial \mathcal{U}}{\partial \mathit{q}}$

21

where $m_{\alpha \beta }$ indicates the element $\alpha \beta$ of the inertia matrix $\mathit{M}\left(\mathit{q}\right)$ .

MOO-based PID control scheme

This section focuses on explaining the MOO-based auto-tuning method with conventional PID control designs. The combination of the PID controller and parameters auto-tuning method results in the MOO-based PID control method for aerial manipulators.

PID control structure for the aerial manipulator

Since we only discuss the stabilization of the aerial manipulator system in this article, the external disturbance/force in equation (17) is assumed to be zero. This assumption will not affect the validness of the proposed method because the external disturbance/force term can be moved to the left side of equation (17), which makes the design of control input $\mathit{u}$ the same. To globally linearize the closed-loop dynamics, the following control input $\mathit{u}$ can be considered

$\mathit{u}=\mathit{M}\left(\mathit{q}\right)\mathit{\sigma }+\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}+\mathit{G}\left(\mathit{q}\right)$

22

where the auxiliary input σ can be partitioned, according to equation (1), that is, $\mathit{\sigma }={\left[{\mathit{\sigma }}_{\mathit{p}}^{\text{T}} {\mathit{\sigma }}_{\Omega }^{\text{T}} {\mathit{\sigma }}_{\eta }^{\text{T}}\right]}^{\text{T}}$ , with ${\mathit{\sigma }}_{\Omega }={\left[{\sigma }_{\varphi }^{\text{T}} {\sigma }_{\theta }^{\text{T}} {\sigma }_{\psi }^{\text{T}}\right]}^{\text{T}}$ . Hence, the globally linearized closed-loop dynamics is given by $\ddot{\mathit{q}}=\mathit{\sigma }$ , that is

$\{\begin{array}{c}\ddot{\mathit{p}}={\mathit{\sigma }}_{\mathit{p}}\\ \ddot{\mathbf{\Omega }} ={\mathit{\sigma }}_{\Omega }\\ \ddot{\mathit{\eta }}={\mathit{\sigma }}_{\mathit{\eta }}\end{array}$

23

Because $\mathit{p}$ , $\Omega$ , and $\mathit{\eta }$ are vectors including more than one variable, a parallel PID controller structure ²⁹ designed for each specific variable can be considered. Take the translational position $x$ , which is an element in $\mathit{p}$ , for example. Its corresponding PID controller is

$\begin{array}{c}{\mathit{u}}_{{\mathit{p}}_x}={\ddot{x}}_{\text{des}}=k_{P\_x}\left(x-x_{\text{des}}\right)+\\ k_{\mathit{I}\_x}\underset{0}{\overset{t}{{\displaystyle \int }}}\left(x-x_{\text{des}}\right)\text{d}t+k_{D\_x}\left(\dot{x}-{\dot{x}}_{\text{des}}\right)\end{array}$

24

where $x_{\text{des}}$ is the desired quadrotor’s $x$ position, ${\dot{x}}_{\text{des}}$ is the rate of desired x position change, $k_{P\_x}$ is the proportional gain, $k_{\mathit{I}\_x}$ is integral gain, and $k_{D\_x}$ is the derivative gain. Other elements in $\mathit{p}$ , $\Omega$ , and $\mathit{\eta }$ are not explained in detail here because they can be achieved in a similar way.

MOO-based PID control parameters auto-tuning

As shown in Figure 2, the multi-objective evolutionary algorithm (MOEA) is used to automatically tune the PID control parameters. The integration of parameters auto-tuning process and the PID controller results in the MOO-based PID controller for aerial manipulators. This controller includes two PID control loops. One is for the quadrotor and another one is for the robotic arm. Each PID controller has a set of control parameters. Because the quadrotor has six specific variables, the quadrotor PID controller includes a set of (6 × 3) parameters. Similarly, the robotic arm PID controller contains a set of (3 × 3) parameters. Based on the difference between the desired and actual motions as well as the difference between current and previous control inputs, the MOEA finds optimal parameters to meet different criteria. Among all the MOEAs, we select the improved NSGA-II ³⁰ as the specific optimization algorithm due to its efficiency and the ease to be implemented in practical applications. While traditional optimization methods convert the MOO to single-objective problems by emphasizing one particular Pareto-optimal solution at a time and running many times, the NSGA-II algorithm is more efficient because its evolutionary mechanism can find multiple Pareto-optimal solutions in one single run. Compared with other popular evolutionary algorithms developed in the same time period such as Pareto-archived evolution strategy (PAES) ³¹ and strength-Pareto evolutionary algorithm (SPEA) ³² and the original NSGA, ³³ NSGA-II has an outperforming performance in terms of finding a diverse set of solutions and in converging near the true Pareto-optimal set because it has a better sorting algorithm and incorporates elitism. It has also been proved to be easily implemented in practical problems like generation expansion planning, ³⁴ water resources management, ³⁵ and distributed hydrologic modeling ³⁶ with good results. Due to the complexity of the aerial manipulator system, the auto-tuning process based on NSGA-II is off-line as the inside evolutionary scheme costs considerable time like minutes or hours to find optima. In this work, it takes 30 minutes to find the optimal PID parameters.

Figure 2.

MOO-based PID control structure for the quadrotor-arm system. MOO: multi-objective-optimization; PID: proportional–integral–derivative.

Two objectives are considered for the NSGA-II algorithm. One is the reliable static–dynamic performance, and the other is the smooth control. Both these two objectives are important for the aerial manipulator system to keep stability, but they cannot be achieved simultaneously because they conflict with each other. In general, seeking for some dynamic specifications usually causes the large variance of control law and significant oscillating of the actuators. The performance balance between these two objectives will be discussed in the next section. The reliable static–dynamic performance is represented by $G_1={\displaystyle \sum }_{k=1}^T{\displaystyle \sum }_{i=1}^9{e}_i^2\left(k\right)$ . Minimizing $G_1$ aims to provide good reference tracking and better disturbance rejection. The objective of smooth control is represented by $G_2={\displaystyle \sum }_{k=1}^T{\displaystyle \sum }_{j=1}^9\frac{1}{2}{\left(u_j\left(k\right)-u_j\left(k-1\right)\right)}^2$ so that minimizing $G_2$ aims to smooth the control signals and avoid significant oscillation of the actuators. Hence, the MOO-based PID control problem can be mathematically formulized as

$\underset{\left\{K_P, K_I, K_D\right\}}{\text{min}}{G}_1\left(\left\{K_P, K_I, K_D\right\}\right), G_2\left(\left\{K_P, K_I, K_D\right\}\right)$

25

$\text{s.t} \text{LB}\le \left\{K_P, K_I, K_D\right\}\le \text{UB}$

where $e_i\left(k\right)$ is the error between the desired and actual output, $T$ is the number of total indices of the time discretization. $\left\{K_P, K_I, K_D\right\}$ groups all the PID control parameters together with a vector-shaped lower bound $\text{LB}$ and upper bound $\text{UB}$ .

Numerical simulations

Simulation setup

In this work, all simulations are completed in MATLAB/Simulink [version R2015b]. The parameters of the aerial manipulator system ¹² are listed in Table 1.

Table 1.

Quadrotor simulation parameters.

Parameters	Values	Descriptions
${\mathit{I}}_{bx}$ , $\hspace{1em}{\mathit{I}}_{by}$	1.24 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around quadrotor’s X-, Y-axis
${\mathit{I}}_{bz}$	2.48 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around quadrotor’s Z-axis
$m_b$	2 kg	Mass of quadrotor
${\gamma }_t$	3.13 × 10−5	Thrust factor
${\gamma }_d$	7.5 × 10−7	Drag factor
l	0.1 m	Distance between the rotor center and quadrotor center
$m_1$ , $m_2$ , $m_3$	0.1 kg	Mass of the first, second, and third link of the arm
$L_1$ , $L_2$ , $L_3$	0.1 m	Length of the arm’s first, second, and third link
${\mathit{I}}_{x1}$ , $\hspace{1em}{\mathit{I}}_{y1}$	0.0025 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around first link’s X, Y-axis
${\mathit{I}}_{z1}$	0.0075 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around first link’s Z-axis
${\mathit{I}}_{x2}$ , $\hspace{1em}{\mathit{I}}_{y2}$	0.0025 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around second link’s X, Y-axis
${\mathit{I}}_{z2}$	0.0075 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around second link’s Z-axis
${\mathit{I}}_{x3}$ , $\hspace{1em}{\mathit{I}}_{y3}$	0.0025 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around third link’s X, Y-axis
${\mathit{I}}_{z3}$	0.0075 $\text{kg}\cdot {\text{m}}^2$	Moment of inertia around third link’s Z-axis

The parameters of NSGA-II algorithm are listed in Table 2. As mentioned before, the lower bound LB and upper bound UB should be vectors corresponding to each control variable. Because we consider all the boundary values to be the same, there is only one value in Table 2, for simplicity. In fact, the bounds of the gain parameters may be difficult to be determined beforehand in practical situations. They are assigned values here only for a convenient calculation and can be practically expanded to a larger value to include more possible solutions. Other NSGA-II-related parameters are empirically chosen.

Table 2.

Parameters used for NSGA-II.

Descriptions	Values
Population size	60
Number of generations	50
Crossover probability	0.9
Mutation probability	0.33
Lower bounds of the gain parameters	[0 0 0]
Upper bounds of the gain parameters	[50 50 50]

NSGA-II: Non-Dominated Sorting in Genetic Algorithms-II.

To compare the performance of different PID parameters, the quadrotor performs a basic but important action within the aerial manipulation. Specifically, the quadrotor is supposed to follow a planar and circular trajectory: $x=\text{sin}\left(0.2t\right)$ , $y=\text{cos}\left(0.2t\right)$ , $z=1$ . The initial configuration of the aerial manipulator system is $\mathit{p}=\left[0, 0, 3\right]$ (m), $\Omega =\left[0, 0, 0\right]$ (°), and $\mathit{\eta }=\left[-45, 45, 45\right]$ (°). The arm moves from configuration A to configuration B ( $\mathit{\eta }=\left[0, 45, 45\right]$ ). This movement is intentionally designed to include a stable tracking process in the whole simulation. The total simulation time is 32 s, corresponding to a complete circle, and the sampling time is 0.02 s. A visualization of the simulated process can be seen in Figure 3. The arm should move from initial configuration A (red solid lines) to desired configuration B (green dotted lines) with the quadrotor following the circular trajectory (purple dashed lines).

Figure 3.

Visualization of simulated movements. The quadrotor is designed to follow a circular trajectory. The arm is designed to move from configuration A to configuration B when the quadrotor starts flying.

In the next section, we will show a comparison between results using manual-tuned parameters (such as using Ziegler–Nichols method ³⁷ ) and auto-tuned parameters. To further discuss the performance of the automatic tuning method, we will also provide the simulation results of a state-of-the-art control method, model predictive control (MPC). MPC is an advanced process control method which calculates the future state of the system and organizes the control action accordingly. It predicts the future of the system and sends control signals in such a manner that it reduces a cost function defined as the error between the output and the desired tracking point over a particular prediction horizon. ³⁸ It has been proven in many complex systems ^{38 –41} to provide better performances, especially stronger robustness against uncertainties, than a majority of current control methods like PID control and sliding mode control. Two key parameters, the control horizon (M) and the prediction horizon (N), are both chosen to be 10 in this article.

Simulation results

The necessity of automatic tuning

The performances of several different sets of manual-tuned PID parameters are shown in Table 3. The dynamic responses for different PID parameter sets are totally different. For example, the overshoot of link 1 for parameter set 1 is 4.4° less than set 2, but the stabilization time of link 1 is 4.3 s longer than set 2. This makes sense because the overshoot usually behaves inversely with respect to the stabilization time. If we want the aerial manipulator system to stabilize quickly, we need to have a larger control output during the beginning phase, which leads to a larger overshoot. This makes it difficult to judge which parameter set is better because each application has a specific requirement. Just for an illustration, we will consider the less overshoot of links as a better performance in the remaining article. Due to a careful manual tuning by experts, the performances of set 3 and set 4 are much better than set 1 and set 2. Considering the link’s overshoot, set 4 is the best in all the four sets.

Table 3.

Manually tuned PID parameter values and corresponding performances.^a

		Set 1	Set 2	Set 3	Set 4 (best)
Parameters	[ $K_{P\_x}$ , $K_{\mathit{I}\_x}$ , $K_{D\_x}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_y}$ , $K_{\mathit{I}\_y}$ , $K_{D\_y}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_z}$ , $K_{\mathit{I}\_z}$ , $K_{D\_z}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_\varphi }$ , $K_{\mathit{I}\_\varphi }$ , $K_{D\_\varphi }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_\theta }$ , $K_{\mathit{I}\_\theta }$ , $K_{D\_\theta }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_\psi }$ , $K_{\mathit{I}\_\psi }$ , $K_{D\_\psi }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_{\eta }_1}$ , $K_{\mathit{I}\_{\eta }_1}$ , $K_{D\_{\eta }_1}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
	[ $K_{P\_{\eta }_2}$ , $K_{\mathit{I}\_{\eta }_2}$ , $K_{D\_{\eta }_2}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
	[ $K_{P\_{\eta }_3}$ , $K_{\mathit{I}\_{\eta }_3}$ , $K_{D\_{\eta }_3}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
Performances	Overshoot of link 1 (°)	20.1	24.5	4.8	4.6
	Stabilization time of link 1 (s)	10.6	6.3	5.0	3.5
	Overshoot of link 2 (°)	9.7	11.8	1.1	1.0
	Stabilization time of link 2 (s)	8.9	8.0	7.2	6.7
	Overshoot of link 3 (°)	10.4	12.5	1.1	1.1
	Stabilization time of link 3 (s)	9.0	8.0	7.3	6.8
	Maximum absolute deviation of quadrotor’s positions (m)	0.11	0.14	0.08	0.08
	Maximum absolute deviation of quadrotor’s angles (°)	6.8	8.3	2.8	2.4

PID: proportional–integral–derivative.

^aEach link length L = 0.1 m, each link mass m = 0.1 kg.

If we change the values of arm length and mass to be $L_1=L_2=L_3=0.05 \hspace{0.17em}\text{m}$ , $m_1=m_2=m_3=0.05 \hspace{0.17em}\text{kg}$ , the results, obtained by using the same parameter sets as in Table 3, are shown in Table 4. Although employing the same control parameters, the dynamic responses of the new aerial manipulator system have changed. For instance, the overshoot of link 1 for the new system is 8.9° when using parameter set 4. It is 4.3° larger than the corresponding value in Table 3. If there is a strict limitation on the link overshoot in practice like flying through a narrow hallway, set 4’s performance can be non-acceptable. If considering the link’s overshoot of all the four sets for the new arm, the best set changes to set 3. This change indicates that the control parameters may need to be tuned each time because the quadrotor may carry different robotic arms to execute different tasks. Unfortunately, this frequent parameter tuning is not always guaranteed with the manual tuning method due to the lack of experts and the high costs of time.

Table 4.

PID parameter values and corresponding performances for different arm settings.^b

		Set 1	Set 2	Set 3 (best)	Set 4
Parameters	[ $K_{P\_x}$ , $K_{\mathit{I}\_x}$ , $K_{D\_x}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_y}$ , $K_{\mathit{I}\_y}$ , $K_{D\_y}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_z}$ , $K_{\mathit{I}\_z}$ , $K_{D\_z}$ ]	[8, 0, 2]	[8, 0, 2]	[15, 0, 3]	[12, 0, 3]
	[ $K_{P\_\varphi }$ , $K_{\mathit{I}\_\varphi }$ , $K_{D\_\varphi }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_\theta }$ , $K_{\mathit{I}\_\theta }$ , $K_{D\_\theta }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_\psi }$ , $K_{\mathit{I}\_\psi }$ , $K_{D\_\psi }$ ]	[10, 0, 2]	[10, 0, 2]	[30, 0, 5]	[30, 0, 5]
	[ $K_{P\_{\eta }_1}$ , $K_{\mathit{I}\_{\eta }_1}$ , $K_{D\_{\eta }_1}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
	[ $K_{P\_{\eta }_2}$ , $K_{\mathit{I}\_{\eta }_2}$ , $K_{D\_{\eta }_2}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
	[ $K_{P\_{\eta }_3}$ , $K_{\mathit{I}\_{\eta }_3}$ , $K_{D\_{\eta }_3}$ ]	[4, 0, 2]	[16, 0, 4]	[16, 0, 8]	[20, 0, 10]
Performances	Overshoot of link 1 (°)	24.6	28.3	7.6	8.9
	Stabilization time of link 1 (s)	11.2	7.5	5.9	5.8
	Overshoot of link 2 (°)	12.2	13.6	2.8	5.0
	Stabilization time of link 2 (s)	9.8	8.5	6.8	7.6
	Overshoot of link 3 (°)	13.7	14.8	2.9	5.0
	Stabilization time of link 3 (s)	9.8	8.6	6.8	7.7
	Maximum absolute deviation of quadrotor’s positions (m)	0.18	0.22	0.08	0.11
	Maximum absolute deviation of quadrotor’s angles (°)	14.2	10.8	3.2	4.8

PID: proportional–integral–derivative.

^bEach link length L = 0.05 m, each link mass m = 0.05 kg.

The results from both Tables 3 and 4 show the importance and necessity of good tuning of PID control parameters. The well-tuned parameters lead to better control performances, generally satisfying the goals of each task. The major drawback of manual tuning is the time-consuming involvement of humans, which could be unacceptable in practical applications. Hence, it is necessary and convenient for the controller to have a set of automatically tuned parameters, adapting to different environments and systematical uncertainties.

Performance evaluation of MOO-based PID control

In this article, auto-tuning is achieved using MOO, NSGA-II. In fact, there are multiple sets of PID parameters that are obtained by NSGA-II. The set of these solutions, called the Pareto front, ⁴² can be seen in Figure 4. All these solutions are regarded as optimal because all of them satisfy the two objectives simultaneously and no objective can be improved without sacrificing the other objective. The differences between solutions indicate various trade-offs between the objectives.

Figure 4.

Pareto front for MOO in terms of integrated time square error and control rate. The solutions to the left of point A focus more on minimizing the integrated time square error while the solutions to the right of point B focus more on minimizing the control rate. The solutions between points A and B are more balanced between the integrated time square error and control rate. MOO: multi-objective-optimization.

The two axes of Figure 4 are two objectives, $G_1$ and $G_2$ . As shown in Figure 4, there are two significant points, A and B. The values of the two objectives $G_1$ and $G_2$ at point A are 90,000 and 0.0520 while the two values at point B are 113,000 and 0.0031. This difference is large enough to show a trade-off between objectives $G_1$ and $G_2$ at points A and B. The solution at point A focuses more on minimizing the integrated time square error ( $G_1$ ) while the solution at point B makes more efforts to keep the control rate ( $G_2$ ) minimal. In fact, the trade-offs can also be indicated by other solutions. The values of objective $G_2$ decrease quickly as objective $G_1$ increases to the left of point A. To the right of point B, the value of objective $G_2$ almost remains the same even while objective $G_1$ increases. Between points A and B, we have a smooth change of objective $G_2$ with respect to objective $G_1$ . The solutions between points A and B can be regarded as considering both objectives to be approximately equal, which is optimal for most situations.

Due to the trade-offs between two objectives in the Pareto front, the dynamic responses of each solution are totally different. For simplicity, we randomly pick two solutions between points A and B out of all the solutions to analyze their performances. The performance comparison can be seen in Table 5. Link 1’s overshoot using parameter set 6 is zero, which is better than set 5’s 1.6°. However, the stabilization time for link 1 using set 6 is 1.3 s longer than set 5. The conflicts and trade-offs can give users some suggestions on selecting the most suitable parameters. If users focus more on the control smoothness, then parameter sets close to point B in Figure 4 are recommended. If the integrated time square error is more important, then parameter sets close to point A are suggested.

Table 5.

Performances of auto-tuned PID parameter sets 5 and 6.

		Set 5	Set 6
Parameters	[ $K_{P\_x}$ , $K_{\mathit{I}\_x}$ , $K_{D\_x}$ ]	[13.15, 0.07, 2.94]	[15.84, 0.06, 3.02]
	[ $K_{P\_y}$ , $K_{\mathit{I}\_y}$ , $K_{D\_y}$ ]	[11.68, 0.03, 3.25]	[14.28, 0.04, 4.01]
	[ $K_{P\_z}$ , $K_{\mathit{I}\_z}$ , $K_{D\_z}$ ]	[18.46, 0.06, 3.99]	[20.06, 0.03, 4.49]
	[ $K_{P\_\varphi }$ , $K_{\mathit{I}\_\varphi }$ , $K_{D\_\varphi }$ ]	[26.82, 0.09, 5.01]	[25.53, 0.07, 3.41]
	[ $K_{P\_\theta }$ , $K_{\mathit{I}\_\theta }$ , $K_{D\_\theta }$ ]	[30.14, 0.11, 5.18]	[28.61, 0.06, 3.93]
	[ $K_{P\_\psi }$ , $K_{\mathit{I}\_\psi }$ , $K_{D\_\psi }$ ]	[32.86, 0.09, 6.20]	[29.84, 0.05, 4.07]
	[ $K_{P\_{\eta }_1}$ , $K_{\mathit{I}\_{\eta }_1}$ , $K_{D\_{\eta }_1}$ ]	[20.34, 0.01, 9.54]	[24.58, 0.01, 15.54]
	[ $K_{P\_{\eta }_2}$ , $K_{\mathit{I}\_{\eta }_2}$ , $K_{D\_{\eta }_2}$ ]	[19.42, 0.01, 8.98]	[22.05, 0.01, 13.45]
	[ $K_{P\_{\eta }_3}$ , $K_{\mathit{I}\_{\eta }_3}$ , $K_{D\_{\eta }_3}$ ]	[19.98, 0.01, 9.06]	[23.12, 0.01, 13.33]
Performances	Overshoot of link 1 (°)	1.6	0
	Stabilization time of link 1 (s)	2.1	3.4
	Overshoot of link 2 (°)	0.8	0.4
	Stabilization time of link 2 (s)	5.8	7.1
	Overshoot of link 3 (°)	0.9	0.5
	Stabilization time of link 3 (s)	5.8	7.0
	Maximum absolute deviation of quadrotor’s positions (m)	0.064	0.062
	Maximum absolute deviation of quadrotor’s angles (°)	2.5	2.6

Performance comparison between manual-tuning, auto-tuning, and MPC

For convenience, we directly use the auto-tuned PID parameter set 5 in Table 5 to compare its performances with those of other different control methods. The performance comparison between the auto-tuned parameter set 5 and manual-tuned parameter set (e.g. set 1 and set 4 in Table 3) is shown in Table 6. Also, as mentioned at the end of simulation setup, we add the MPC simulation results in Table 6 for comparison. Additionally, Figure 5 is presented to show the change of joint angles using PID parameter set 1, set 4, set 5, and MPC. Correspondingly, Figure 6 shows the quadrotor trajectories. Figure 7 uses a bar plot to present a more intuitive comparison of various system responses in Table 6. This case uses the original aerial manipulator system parameters in Table 1.

Figure 5.

Link performance comparison using PID parameter set 1 (manual-tuned), 4 (best manual-tuned), 5 (auto-tuned), and MPC. The auto-tuned performances are comparable to MPC because the differences between the red solid lines and green dashed–dotted lines are small. Although the best manual-tuned performances (purple dotted lines) are close to auto-tuned and MPC performances for links 2 and 3, the best manual-tuned link 1 performance is notably worse than the auto-tuned and MPC performance. The manual-tuned performances (blue dashed lines) are worst with significant oscillations. PID: proportional–integral–derivative; MPC: model predictive control.

Figure 6.

Quadrotor trajectories using PID parameter set 1 (manual-tuned), 4 (best manual-tuned), 5 (auto-tuned), and MPC. The auto-tuned parameter set 5 has comparable performances with MPC, but better performances than manual tuned parameter sets in terms of less deviations and faster stabilization time. PID: proportional–integral–derivative; MPC: model predictive control.

Figure 7.

Performance comparison of PID parameter set 1 (manual-tuned), 4 (best manual-tuned), 5 (auto-tuned), and MPC. All performances including the overshoot, stabilization time, quadrotor position deviation, and quadrotor angle deviation are considered better with smaller values. The auto-tuned and MPC performances are comparable while they are better than the manual-tuned performances. PID: proportional–integral–derivative; MPC: model predictive control.

Table 6.

Performance comparison of PID parameter set 1 (manual-tuned), 4 (best manual-tuned), 5 (auto-tuned), and MPC.

		Set 1 (manual- tuned)	Set 4 (best manual-tuned)	Set 5 (auto-tuned)	MPC
Parameters	[ $K_{P\_x}$ , $K_{\mathit{I}\_x}$ , $K_{D\_x}$ ]	[8, 0, 2]	[12, 0, 3]	[13.15, 0.07, 2.94]	M = 10; N = 10
	[ $K_{P\_y}$ , $K_{\mathit{I}\_y}$ , $K_{D\_y}$ ]	[8, 0, 2]	[12, 0, 3]	[11.68, 0.03, 3.25]
	[ $K_{P\_z}$ , $K_{\mathit{I}\_z}$ , $K_{D\_z}$ ]	[8, 0, 2]	[12, 0, 3]	[18.46, 0.06, 3.99]
	[ $K_{P\_\varphi }$ , $K_{\mathit{I}\_\varphi }$ , $K_{D\_\varphi }$ ]	[10, 0, 2]	[30, 0, 5]	[26.82, 0.09, 5.01]
	[ $K_{P\_\theta }$ , $K_{\mathit{I}\_\theta }$ , $K_{D\_\theta }$ ]	[10, 0, 2]	[30, 0, 5]	[30.14, 0.11, 5.18]
	[ $K_{P\_\psi }$ , $K_{\mathit{I}\_\psi }$ , $K_{D\_\psi }$ ]	[10, 0, 2]	[30, 0, 5]	[32.86, 0.09, 6.20]
	[ $K_{P\_{\eta }_1}$ , $K_{\mathit{I}\_{\eta }_1}$ , $K_{D\_{\eta }_1}$ ]	[4, 0, 2]	[20, 0, 10]	[20.34, 0.01, 9.54]
	[ $K_{P\_{\eta }_2}$ , $K_{\mathit{I}\_{\eta }_2}$ , $K_{D\_{\eta }_2}$ ]	[4, 0, 2]	[20, 0, 10]	[19.42, 0.01, 8.98]
	[ $K_{P\_{\eta }_3}$ , $K_{\mathit{I}\_{\eta }_3}$ , $K_{D\_{\eta }_3}$ ]	[4, 0, 2]	[20, 0, 10]	[19.98, 0.01, 9.06]
Performances	Overshoot of link 1 (°)	20.1	4.6	1.6	0
	Stabilization time of link 1 (s)	10.6	3.5	2.1	3.5
	Overshoot of link 2 (°)	9.7	1.0	0.8	0.2
	Stabilization time of link 2 (s)	8.9	6.7	5.8	5.7
	Overshoot of link 3 (°)	10.4	1.1	0.9	0.2
	Stabilization time of link 3 (s)	9.0	6.8	5.8	5.6
	Maximum absolute deviation of quadrotor’s positions (m)	0.11	0.08	0.064	0.060
	Maximum absolute deviation of quadrotor’s angles (°)	6.8	2.4	2.5	3.2

MPC: model predictive control; PID: proportional–integral–derivative.

It is worthy to notice that all the performances of parameter set 5 are better than set 1, as the optimization method is used. It can be seen in Figure 5 that there are less oscillations of parameter set 5 than set 1. This means the changes of control inputs are less and smoother, meeting objective $G_2$ . The integrated time square error of set 5 can also be clearly seen to be much less than set 1, satisfying objective $G_1$ . Although there is a relatively larger overshoot for link 1, the stabilization time of parameter set 5 is shorter than MPC. The performance of auto-tuned set 5 is also better than best manual-tuned set 4 although they are close. This indicates expert-tuned results may still be not optimal and can be improved. From Figure 6, all the quadrotor trajectories have an obvious deviation at the first several seconds. However, the trajectory using PID parameter set 5 and MPC quickly follows the desired trajectory around 6 s while the one using PID parameter set 1 exactly tracks around 11 s. Hence, the performance of parameter set 5 can be regarded as comparable to MPC. There is no doubt that the performance of MPC will be better, possibly outperforming our MOO-based PID method, if increasing the value of the control horizon (M) and the prediction horizon (N). However, it multiples the time used to find solutions, which is 15 min longer than the MOO-based PID method (the offline optimization takes 30 min) in current MPC setting. In terms of both time cost and performance, the MOO-based PID method is preferred.

Since there may be more than two objectives in practical applications, we also simulate a case with three objectives by adding a new one specifying no overshoot for the link 1 joint angle. Similar to the two-objective case, there is also a Pareto front with numerous solutions. For convenience, we picked out a solution that happens to be the same with parameter set 6 in Table 5. It does not have an overshoot for link 1’s joint angle, which is verified to meet the newly added objective. However, the three-objective optimization needs a longer tuning time for the whole simulation, which is about three times than the time used for the two-objective case. This indicates that we can add more and new objectives into the MOO-based PID control formulation, that is, equation (25), but the feasible number of objectives may have an upper limit due to the additional computational complexity. MPC needs much more time than MOO-based PID if adding new objectives due to its own high complexity, which sometimes limits its practical application. Also, the PID controller is most widely used controller in the industry, and it may be easier and more convenient to keep the PID structure unchanged rather than developing a new and complex MPC controller.

Conclusion

In this article, we present a MOO-based control parameters auto-tuning method for the aerial manipulator system to be better applied in practical applications. This method uses NSGA-II, an MOEA, to automatically find optimal parameters for different control goals that are related to the desired and actual measurements of the aerial manipulator system. The auto-tuning performances are comparable to a well-designed state-of-the-art model predictive controller while both are better than manual-tuned performances.

The reasons about the better performances of auto-tuned controllers can be concluded as (1) the flexibility to changing conditions and (2) inherent satisfaction of multiple objectives. The practical changing conditions (including both robot dynamics and task environments) easily make current manual-tuned controllers that are designed for specific configurations or tasks not suitable. In contrast, the auto-tuning method has more flexibilities to these changes due to the inside optimization mechanism. Even for the case that humans can tune parameters for each situation, although this is nearly impossible, it costs significantly more time than the auto-tuning method. Furthermore, it is sometimes extremely difficult for even experts to tune appropriate parameters for multiple control objectives. Humans may not be able to quantitatively and precisely balance between different control objectives due to human’s fuzzy logic and approximate reasoning. Instead, the auto-tuning method can use an exact mathematical form to accurately evaluate trade-offs between different objectives. Compared with the advanced MPC, the MOO-based auto-tuning method is easier to be implemented and extended due to less computational complexity and more flexibility to add additional objectives.

Although the proposed MOO-based auto-tuning method is an off-line method, it can still be applied to many practical applications that require relatively less time sensitivity. If the aerial manipulator is used to transport different loads in different tasks, for example, then this auto-tuning method could be more efficient than a manual-tuned aerial manipulator to adapt to a specific load. Moreover, there could be different combinations of UAVs and robotic arms to generate an aerial manipulator system. It will be easier to use the auto-tuning method to handle such a large variety.

Despite the advantages, the MOO-based auto-tuning method can still be improved. For example, it is worth investigating how to select the best solution from the Pareto-front. This may be solved using graphic methods to visualize the distribution of solutions in Pareto-front. It could be especially useful when the solution dimension is high. When the number of objectives becomes more than three, the complexity of MOO greatly increases and hence the efficiency correspondingly decreases. It may be necessary to find another efficient optimization algorithm for this problem. We also plan to conduct experiments on real aerial manipulators to see how the presented method works with unavoidable noises.