# An Innovative Aerial Manipulator with Tandem Ducted Fans: Modeling, Control, and Simulation

Authors: Yibo Zhang; Wei Fan; Changle Xiang; Bin Xu; Tianfu Ai; Lei Yuan; Yang Liu
Journal: Complexity
 (2020)

---

## Abstract

This paper proposes an innovative ducted fan aerial manipulator, which is particularly suitable for the tasks in confined environment, where traditional multirotors and helicopters would be inaccessible. The dynamic model of the aerial manipulator is established by comprehensive mechanism and parametric frequency-domain identification. On this basis, a composite controller of the aerial platform is proposed. A basic static robust controller is designed via H-infinity synthesis to achieve basic performance, and an adaptive auxiliary loop is designed to estimate and compensate for the effect acting on the vehicle from the manipulator. The computer simulation analyses show good stability of the aerial vehicle under the manipulator motion and good tracking performance of the manipulator end effector, which verify the feasibility of the proposed aerial manipulator design and the effectiveness of the proposed controller, indicating that the system can meet the requirements of high precision operation tasks well.


---

## Body

## 1. Introduction

In recent years, the application of unmanned autonomous robots is increasingly diverse, and the interaction between the autonomous robot system and the environment is developed from information interaction (such as sound, light, picture) to physical interaction (replacing manpower to complete operation work) [1]. On the one hand, using autonomous robots for operation work can greatly save labor costs and significantly improve work efficiency; on the other hand, it can liberate people from heavy labor, especially the labor in dangerous and harmful environments. At present, some ground mobile robots have been applied to postearthquake rescue [2], some underwater robots have been applied to oceanic biological sample collection [3], and some space robots have been used in space exploration [4]. Although ground mobile robots, underwater robots, and space robots have been widely used, the application of aerial robots is still in its infancy all over the world. The current unmanned aerial vehicles play an important role in monitoring activities, such as aerial photography and high-voltage line inspection [5], but do not have the ability to physically interact with the environment. However, the aerial operation robot (aerial manipulator) has great application value in the following three aspects: (1) replacing manpower to complete dangerous tasks, such as urban antiterrorism and high-rise building firefighting [6]; (2) replacing manpower to improve efficiency, as in wide-area scientific examination and collection [7]; (3) replacing manpower to reduce costs, as in infrastructure maintenance and remote operation in complex environment [8].The configuration of the aerial robot commonly consists of an aerial platform and an operation manipulator, but their combination brings completely new features to the system. The main challenges in the aerial robot system design include two aspects. First, since the aerial robot usually works in confined spaces that require operation, it should have high traffic ability and contact ability with the complex environment. Under this premise, it also should have as large operating payload as possible. Second, there is a serious coupling effect between the aerial platform and the robot arm, which makes the system face enormous challenge in terms of stability and manipulation accuracy, and physical contact and manipulation further exacerbate it. The research of aerial robots has begun to attract worldwide attention since 2010, and the representative work includes AIRobots Project [9], ARCAS Project [10], AEROWORKS Project [11].In respect to structure, most of aerial robots use helicopter or multirotor as aerial platform. Yale University designed a helicopter with a single DOF (degree of freedom) underactuated gripper, as shown in Figure1(a), and studied the quasi-static compliance control problem during the grabbing process through PID control method [12]. University of Drexel also used helicopter as base platform to expand the 1-DOF gripper into a multi-DOF manipulator and completed the grasping and placement tasks of cylindrical objects [13]. National Taipei University of Technology designed a small quadrotor with a 2-DOF manipulator and carried out simulation analysis of kinematic control in 2D plane based on a simplified vector model [14]. University of Seville designed a dual-arm system based on quadrotor, as shown in Figure 1(b), which expands the maneuvering range and enables more complex operations [15]. However, because of the inherent characteristics of open rotor, the helicopter or multirotor structure is unable to interact with the environment closely, and the arm can only manipulate objects above [10] or below [12], which greatly limits the application scenarios. In order to grab the target side-on, Johns Hopkins University proposed a very long arm to avoid the rotor disc of the quadrotor platform [16], but the changes in the center of mass and moment of inertia caused by the long arm exert serious impacts on the stability and effective payload of the system. Compared to helicopter and multirotor, ducted fan has greater thrust in a more compact structure [17], having contact capability and increasing the payload of the system, which is more suitable as the platform of an aerial robot. Figure 2 shows the comparison of flight area requirements of different types of aerial robots near a wall (namely, our innovative ducted fan aerial robot, quadrotor, and helicopter). It can be seen that, under the same effective load, ducted fan aerial robot can be closer to the target and operate it from side-on with a smaller joint motion range and can pass narrower confined space and be safer. University of Bologna completed a series of studies on contact dynamics of the single ducted fan with vertical wall [18], but their platform cannot carry the manipulator due to the payload and controllability limitations.Figure 1
Examples of aerial robots. (a) Yale University. (b) University of Seville.
(a)(b)Figure 2
Comparison of flight area requirements of different types of aerial robots.With regard to control, because of the serious coupling between aerial platform and manipulator, on the one hand, the manipulation process will have three effects on the aerial platform: (1) the effects of the manipulator gravity, and the gravity moment caused by the noncoincidence of the center of mass between the aerial vehicle and the manipulator, (2) the inertia force and moment generated by manipulator dynamics, and (3) the impact of the external environment during contacting and operating; on the other hand, the drift of the aerial platform also affects the positioning accuracy of the end effector [19]. Some studies ignored the existence of manipulator and only considered the aerial vehicle [20], and some tried to improve the robustness of the aerial vehicle basic controller to ensure system stability [21–23], but both of them are only suitable for small size arm (e.g., lightweight 1-DOF gripper) and lightweight object. DLR in [20] proposed an impedance controller to stabilize the system in the presence of external forces, ignoring the coupling effects between the aircraft and the arm by limiting the motion of the arm. Shenyang Institute of Automation of CAS in [23] designed a linear LQR controller of a helicopter with 1-DOF arm. The simulation results show that the controller has good performance when the arm moves in a small range near the equilibrium point, but the LQR controller cannot stabilize the system when the swing range of the arm is relatively large. AGH University of Science and Technology considered the change in the system’s center of mass caused by the arm movement and compensated for the influence through variable-parameter PID control to achieve system’s stability [24], while the University of Pennsylvania considered the impact of the payload acting on the system’s center of mass [25]. The University of Seville designed a controller to compensate for the change of both the center of mass and the moment of inertia when the arm is in different positions, but the dynamic characteristics of the arm are not considered [26]. The University of Naples Federico II designed a Cartesian space impedance controller based on the integrated dynamic model of the aerial manipulator, which fully considered the coupling effect between the aerial vehicle and the robot arm [27]. However, the controller is highly dependent on the model and has a complicated structure, which is not easy in practice.Based on the discussion above, the contribution of this paper mainly includes two parts. First, an innovative aerial manipulator based on tandem ducted fans is proposed, which has both great trafficability and effective payload. The dynamic model of the aerial manipulator is established by comprehensive mechanism and parametric frequency domain identification. Owing to the small lateral size, the aerial manipulator can easily realize the omnidirectional manipulation of side-on and below and is particularly suitable for the tasks in confined environment, where traditional multirotor and helicopter would be inaccessible. The application aims of this novel design are canopy sampling in dense forests and insulator lubricating in dense high-voltage wires. Second, a composite controller of the aerial platform is proposed, considering and compensating for both the static and the dynamic disturbances of the manipulator on the vehicle. A basic static robust controller is designed via H-infinity synthesis for basic performance, and an additional adaptive loop is designed for disturbance estimation and compensation from the manipulator to improve the platform stability and the end effector tracking accuracy. The computer simulations verify the effectiveness of the proposed controller.The organization of this paper is as follows. In Section2, the structure of the innovative ducted fan aerial manipulator is introduced and the dynamic model is established. Then, based on the dynamic model, a control-oriented state space model of the aerial platform is derived based on parametric identification to facilitate the controller design. In Section 3, the basic robust controller and the adaptive auxiliary controller are designed and analyzed in detail. In Section 4, the simulation results and analysis of the proposed controller are carried out. Finally, some conclusions are drawn in Section 5.


## 2. Modeling of the Ducted Fan Aerial Manipulator

### 2.1. System Description

The configuration of the proposed aerial manipulator mainly includes two ducted fan systems with coaxial rotors, two sets of control rudders, a 3-DOF manipulator, the control unit, and the landing gears, as shown in Figure3. The pitch channel of the vehicle is controlled by the thrust difference between the front and the rear ducted fans, which is caused by the difference of rotor speed; the roll channel is controlled by the rudder systems setting below the duct; the yaw channel is controlled by the reaction torque difference between the upper and lower rotor disc. The total mass of the aerial vehicle is 4.6 kg, with an effective payload of 2 kg. The manipulator weight is 0.9 kg, with a max grasping weight of 0.5 kg. The parameters of the system are detailed in Table 1, where the mass and the structural dimension parameters of the platform are measured directly, and the moments of inertia and the effective resistance areas along three axes of the vehicle are estimated by the CATIA® 3D model. The airfoil used in the prototype is NACA 0012, and the airfoil parameters can be obtained from [28]. The manipulator is installed on the center of the vehicle body and has 3 DOF (one lumbar joint, one shoulder joint, and one elbow joint) to reach any position in 3D space. The end effector is a gripper, and, from the perspective of our application purpose, the posture of the gripper is not considered. The parameters of the arm are described using standard D-H (Denavit–Hartenberg) method [29], as shown in Table 2.Figure 3
The configuration of the novel aerial manipulator.Table 1
Structural parameters of the system.
ParameterPhysical descriptionValuembVehicle mass (including battery)4.6 kgmrManipulator mass0.9 kgmgraspMax grasping mass0.5 kgIxxInertia tensor of the vehicle aroundx-axis0.092 kg·m2IyyInertia tensor of the vehicle aroundy-axis0.283 kg·m2IzzInertia tensor of the vehicle aroundz-axis0.245 kg·m2pcdDistance between duct center and CG of vehicle0.32 mDDuct diameter0.33 mSxVehicle effective resistance area aroundx-axis0.034 m2SyVehicle effective resistance area aroundy-axis0.127 m2SzVehicle effective resistance area aroundz-axis0.133 m2nBlade number of each disc4cBlade chord length0.027 mθ0Attack angle at the root of blade35 degθrwTorsion rate of blade18 degσBlade tip clearance0.001 mTable 2
D-H parameters of the manipulator.
Linka (m)α (deg)d (m)θ range (deg)10900.08180 [−90, 90]20.1500150 [0, 150]30.1600150 [−150, 0]The hardware of control unit is shown in Figure4. The on-board controller is Emlid® Navio2 based on Raspberry Pi® 3, which integrates dual IMU module, GPS module, barometer module, and 14 PWM output channels. The rotors are driven by four GARTT® motors, and the motors are driven by HOBBYWING® Electronic Speed Controllers (ESC), which are controlled by the main controller via PWM signals. Each duct system has two sets of control rudders, which are driven synchronously by one KST® servo. The servos are also controlled by PWM signals from the main controller. The Dynamixel® XH430 servo is chosen as arm joint servo, since it has both position control mode and torque control mode and has feedback function of actual position, velocity, and torque of the joint. The on-board main controller communicates with the ground station through the 3DR® radio telemetry.Figure 4
The hardware of control unit.


### 2.2. Dynamic Model

The coordinate system of the aerial manipulator is introduced in Figure5. The aerial manipulator is modeled as a multibody system consisting of four interconnected rigid bodies. Let Σ be the earth-fixed Cartesian coordinate frame following the north-east-down rules, Σb be the vehicle body-fixed coordinate frame at the center of mass of the vehicle, Σ0 be the manipulator base-fixed coordinate frame, and Σi (i = 1, 2, 3) be the coordinate frame of each link of the manipulator following the D-H rules. Notice that Σ0 coincides with the origin of Σb, only rotated 90° around the Zb axis. In the following, the superscript i means that the variable is related to the coordinate system Σi.Figure 5
The coordinate system of the aerial manipulator.The dynamic formulation is established using iterative Newton–Euler method [29], which is widely used in dynamic modeling [30, 31]. First, give the definition of the motion states of the system. Let pb=x⋅y⋅zT and Φb=φ⋅θ⋅ψT be the position vector and Euler angle vector of the aerial vehicle platform in earth-fixed frame, and let qi⋅i=1⋅2⋅3 be the joint angle of the manipulator. Let vb and ωb be the velocity vector and angular rate vector of the aerial vehicle, and let vi and ωi be the velocity and angular rate of the origin of link i of the manipulator. Then, the acceleration of the multibody system can be calculated outward iteratively from the aerial vehicle platform to the manipulator end effector as(1)vbb=uvwT=p˙bb=RbTp˙b,ωbb=pqrT=RbTωb=RbTΤbΦ˙b=QbΦ˙b,v00=Rb0vbb,ω00=Rb0ωbb,ωii=Ri−1iωi−1i−1+q˙iezi,ω˙ii=Ri−1iω˙i−1i−1+Ri−1iωi−1i−1×q˙iezi+q¨iezi,v˙ii=Ri−1iω˙i−1i−1×pii−1+ωi−1i−1×ωi−1i−1×pii−1+v˙i−1i−1,aCii=ω˙ii×pCii+ωii×ωii×pCii+v˙ii,where the superscript b refers to the body-fixed coordinate system, the superscript i refers to the ith link coordinate system, and in particular i = 0 for the manipulator base-fixed (link 0) coordinate system, as shown in Figure 5. Specifically, vbb and ωbb are the velocity and angular rate of the aerial vehicle with respect to the body-fixed coordinate frame, vii and ωii are the velocity and angular rate of the origin of link i with respect to the ith link coordinate frame. aCii is the acceleration of the center of mass of link i. pii−1 is the position vector of the origin of ith link frame with respect to the i−1th link coordinate frame, pCii is the position vector of the center of mass of ith link with respect to the ith link coordinate frame, and ezi refers to the projection along the z-axis of link i. Rii−1 refers to the transformation matrix from the i−1th link coordinate frame to the ith link coordinate frame. In particular, Rii−1 refers to the transformation matrix from the body-fixed frame to the arm base-fixed (link 0) frame. Rb and Qb are the linear velocity and angular rate transformation matrix between the earth-fixed frame and body-fixed frame expressed in the form of Euler angles.Then, the Newton–Euler dynamics formulation can be derived as(2)Fbb=mbv˙bb+ωbb×vbb,Mbb=Ibω˙bb+ωbb×Ibωbb,Fii=miaCii,Mii=Iiω˙ii+ωii×Iiωii,where mb and mi are the mass of the aerial vehicle and link i, Ib, and Ii are the inertia matrix of the aerial vehicle and link i. Fbb and Fii are the total external force exerted on the aerial vehicle and link i, and Mbb and Mii are the total external moment exerted on the aerial vehicle and link i, which can be calculated inward iteratively from the manipulator end effector to the aerial vehicle platform as(3)fii=Fii+Ri+1ifi+1i+1−RiTmig,nii=Mii+Ri+1ini+1i+1+pi+1i×Ri+1ifi+1i+1+pCii×Fii−RiTmig,τi=niiTezi,Fbb=Faero1+Faero2+Frudder1+Frudder2+Ffus+RbTmbg−R0bR10f11,Mbb=Maero1+Maero2+pcd1b×Faero1+pcd2b×Faero2+pcr1b×Frudder1+pcr2b×Frudder2+Mgyro−R0bR10n11,where fii and nii are the force and moment exerted on link i by link i−1 and Ri is the transformation matrix from the ith link coordinate frame to the earth-fixed coordinate frame. τi is the joint torque of the manipulator.The total force of the aerial vehicle includes five parts: the aerodynamic force of the ducted fans, the resultant force generated by the rudders, the fuselage resistance, the gravity, and the reaction force of the manipulator. The total moment of the aerial vehicle also includes five parts: the aerodynamic moment of the ducted fans, the torque generated by the aerodynamic forces, the torque generated by the rudders, the reaction torque of the manipulator, and the gyro moment produced by rotors, as shown in Figure6. In equation (3), the subscript “aero” refers to the aerodynamic force and moment, subscript 1 denotes the front ducted fan system, and subscript 2 denotes the rear ducted fan system. Same as the above, the subscript rudder refers to the force of the front and the rear rudder system. fus refers to the fuselage resistance, and gyro refers to the gyro torque. pcdb is the position vector of the center of duct with respect to the center of the vehicle, and pcdb is the position vector of the aerodynamic center of rudder with respect to the center of the vehicle. These forces and moments belong to the ducted fan system dynamics, which provide the main thrust and attitude control moments for the aerial vehicle, and will be introduced in detail in the following section.Figure 6
Diagram of the dynamic model.


### 2.3. Ducted Fan System

As mentioned in the Introduction, most of the literature of aerial manipulator uses multirotor as aerial platform. Because of the simple, symmetrical, and decoupled open rotor structure, the rotor aerodynamics is usually simplified as a thrust and a reaction torque in quadratic relation to the rotor speed, ignoring any other effects of the aerodynamic characteristics (e.g., [22, 26, 27]). However, the aerodynamics of ducted fan is significantly different from the traditional open rotor [32]. First, the duct lip causes airflow deflection effect, which changes the inflow at the rotor disc, while additional thrust is generated by the suction flow of the duct lip; second, duct suppresses the rotor tip vortex and reduces the momentum loss; third, the duct exit prevents the airflow from contracting, which increases the outflow pressure of the rotor and thereby increases the thrust. Johnson and Turbe [17] propose a Blade Element Momentum Theory (BEMT) model for ducted-single-rotor system. By improving this method, the dynamic model of our ducted-coaxial-rotor system is proposed. The basic idea of BEMT is to establish a set of relationships between rotor aerodynamic force and induced velocity through momentum theorem and blade element theory, respectively, and then solve the simultaneous equations iteratively. The following is the derivation of the front ducted fan, and the rear one can be obtained in the same way.Momentum Theorem. The inflow process of the coaxial ducted fan under airflow deflection effect is shown in Figure 7. It can be divided into five parts: the free inflow, the deflected airflow whose direction is deflected by duct lip, the airflow at upper rotor disc, the airflow at lower rotor disc, and the outflow of duct exit. The free flow at the ducted fan is(4)Va=uavawaT=vbb−RbTvw,where vbb is the body velocity of the aerial vehicle and vw is the local wind velocity. Express va in the inflow plane; the inflow velocity v0 is calculated as(5)V0=V0cosαi+V0sinαj,V0=Va,α=arctanwaua2+va2.Figure 7
The inflow process model.Then, the deflected airflow, the airflow at upper and lower rotor disc, and the outflow can be calculated by(6)Vt=V0cosαti+V0sinαtj,Vr1=V0cosαti+V0sinαt+vi1j,Vr2=V0cosαti+V0sinαt+vi1+vi2j,Ve=V0cosαti+V0sinαt+v∞j,where vi1 and vi2 refer to the induced velocities of upper and lower rotor, v∞ refers to the overall impact of the ducted fan on free flow, and αt is the deflection angle caused by duct effect at duct lip, which can be expressed as a function of free inflow angle and airflow deflection factor as(7)αt=α+ktπ2−α.According to momentum theorem and kinetic energy theorem, considering the control body from the free inflow to the upper disc exit and from the upper disc exit to the duct exit and the entire duct, respectively, the following equations are obtained:(8)Tr1V0sinαt+vi1=12ρSDV0V0sinαt+vi12−V0sinα2,Tr2V0sinαt+vi1+vi2=12ρSDV0V0sinαt+v∞2−V0sinαt+vi12,T=1+kductTr1+Tr2=ρSDV0V0sinαt+v∞−V0sinα,Td=ρSDV0V0cosαt−V0cosα,where Tr1 and Tr2 are the thrusts of upper rotor and lower rotor in the vertical direction, Td is the momentum resistance in the horizontal plane, kduct is the thrust augment factor of duct, and SD is the rotor disc area.Blade Element Theory. Blade element theory analyzes the aerodynamics of each micro element of the blade, then integrates it along the radial direction of the blade, and averages it along the circumferential direction to solve the forces and moments of the entire rotor. According to equations (5) and (6), the airflow velocities at the blades of upper and lower rotors are(9)Vra1=ura1vra1wra1T=V0cosαtuaua2+va2V0cosαtvaua2+va2V0sinαtT,Vra2=ura2vra2wra2T=V0cosαtuaua2+va2V0cosαtvaua2+va2V0sinαt+vi1T.Since the blade element theory is a common method, the results after integration of the rotor thrust Tri, the reaction torque Qri, the pneumatic rolling torque Lri, and the pneumatic pitching torque Mri of the upper and lower rotor are directly given below. Notice that i = 1 for upper rotor and i = 2 for lower rotor:(10)Tri=n2π∫02π∫0RdTridψ0=nρcCaΩi2R3423+μui2+μvi2θ0−1+μui2+μvi2θrw2−λwi,Qri=n2π∫02π∫0RdQridψ0=nρcCaΩi2R4423θ0−12θrw−λwiλwi+121+μui2+μvi2Cd,Lri=n2π∫02π∫0RdLridψ0=nρcCaΩi2R44μvi23θ0−12λwi−12θrw,Mri=n2π∫02π∫0RdMridψ0=nρcCaΩi2R44μui23θ0−12λwi−12θrw.In the above formula,n is the blade number of each disc, R is the rotor radius, ψ0 is the blade azimuth angle, c is the blade chord length, Ca is the lift coefficient slope, Cd is the resistance coefficient, Ωi is the rotor speed, θ0 is the blade attack angle, and θrw is the blade torsion rate. μui, μvi, and λwi denote the dimensionless ratios of forward velocity, lateral velocity, and inflow velocity:(11)μui=uraiΩiR,μvi=vraiΩiR,λw=vii−wraiΩiR.Solve equations (8) and (10) simultaneously; the aero force and moment of the ducted fan in body-fixed frame can be derived as(12)Faero=−uaua2+va2Td−vaua2+va2Td−TT,Maero=Lr1−Lr2Mr1−Mr2Qr1−Qr2T.Rudder Dynamics. Same as the ducted fan system, the following gives the analysis of the front rudder, and the rear one can be obtained in the same way. The rudder is in the outflow of the duct exit, and its attack angle can be calculated by(13)αv=αe+δ,where αe is the angle of the duct outflow velocity Ve, and δ is the rudder control angle. Then the aero-lift and aero-drag of the rudder can be derived as(14)Lv=12ρVe2SvClv,Dv=12ρVe2SvCdv,where ρ is the airflow density, Sv is the rudder area, and Clv and Cdv are the lift coefficient and drag coefficient, which are the functions of αv. The force of the control rudder in body-fixed frame can then be derived. nv is the number of the rudders in one duct:(15)Frudder=0nvLvcosαe+Dvsinαenv−Lvsinαe+DvcosαeT.In summary, the expression of each component in equation (3) is obtained, and the comprehensive nonlinear mechanism model is established.


### 2.4. Model Identification of the Aerial Platform

The comprehensive nonlinear mechanism model established in the previous section can fully reflect the characteristics of the system in all-envelope range and be used for system performance prediction and comprehensive simulation analysis, but the complex mechanism is not suitable for the controller design process [33]. Since the working scenarios of the aerial vehicle in this paper are hovering or near-hovering conditions, in order to reduce the complexity of nonlinear identification and improve the identification accuracy under the noise of low-cost sensors, the parametric frequency domain identification is performed. First, a simplified parametric model can be derived from the nonlinear mechanism model of the aerial platform by small perturbation assumption under hover equilibrium point; then, a control-oriented state space model of the aerial platform is derived based on frequency domain identification method. The identification model is verified in both time and frequency domains.The nonlinear mechanism model of the aerial vehicle can be described as(16)x˙b=fxb,ub,xb=uvwpqrφθψT,ub=ualturolupituyawT,where ub is the normalization control input vector, referring to the altitude channel, roll channel, pitch channel, and yaw channel, and is mapped to the rotor speed difference and rudder angle (described in 2.1). Find the partial derivative of equation (16) at hover equilibrium point to obtain the system’s Jacobian matrix:(17)x˙b≈Axb−xhover+Bub−uhover,A=∂f∂xb,B=∂f∂ub,fxhover,uhover=0.Then, the simplified parametric matrix A and B to be identified is obtained by small perturbation assumption. Since xhover and uhover are constant and xhover = 0, in the rest of this paper, for the simplicity of the symbol, let xb ≡ xb − xhover, ub ≡ ub − uhover. Because the aerial vehicle platform is inherently unstable, the closed-loop identification is carried out following the procedure in [34]. A simple PID controller is used to ensure the stability during identification process. A set of Chirp sweep signals (introduced in [35]) is performed in each of the input channels, and the data is processed by CIFER® identification tool [35]. The model fitting process is an optimization problem as(18)minJ=∑i=1n20nω∑ω1ωnωWγWgGfi−Gi2+Wp∠Gfi−∠Gi2,where Gfi is the transfer function to be fitted corresponding to Gi, converting from the parametric state space equation, and n is the number of the transfer functions. |G| refers to the amplitude and ∠G refers to the phase. nω is the number of the selected frequency sampling points, and ω1 and ωnω refer to the initial and termination sampling frequency respectively. Wγ, Wg, and Wp are the weighting functions referring to coherence value, amplitude, and phase, respectively, and are chosen following the directions in CIFER®. The frequency response identification results are shown in Figures 8 and 9, and the identified model is(19)A=−0.0876000000−9.801000−0.087600009.801000−0.1178−0.1172−1.0415−0.0012−0.0209−0.05090000−0.680100000000.0940000−1.06990.01320000000−0.01220000000100000000010000000001000,B=00000000−20.931200.0547−0.001508.349400.02600−0.002510.0876−0.028100.00160.13373.3981000000000000.Figure 8
Frequency response identification results of linear velocities.
(a)(b)(c)(d)(e)(f)Figure 9
Frequency response identification results of angular rates.
(a)(b)(c)(d)(e)(f)In order to verify the accuracy of the identified model, a series of flight tests are carried out in time domain, as shown in Figure10. The low-pass filtered bipolar square wave signal, which is different from the identification excitation signal, is used as the excitation signal in each channel, and the evaluation of the model is realized by comparing the output of the identified model simulation and the actual test measurement. Notice that the angular rate errors in Figure 10 are large while other variables match well, because the dynamic characteristics of angular rate are ignored in high frequency in the linearized model and the angular rate channel has relatively large measurement noise. According to [35], Theil inequality coefficient (TIC) can be used to estimate the model fitting accuracy, which is calculated as(20)JTIC=1/n∑i=1nytest−ysimTytest−ysim1/n∑i=1nysimTysim+1/n∑i=1nytestTytest,where ysim and ytest are the output vector of identified model simulation and actual flight test, and n is the number of sampling points. A smaller value of TIC indicates that the accuracy of the identified model is higher, and a value less than 0.25 indicates that the model has good predictive performance [35]. JTIC of our model is 0.1596, indicating that the nominal model has satisfactory accuracy.Figure 10
Time-domain verification of the identified model.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 2.1. System Description

The configuration of the proposed aerial manipulator mainly includes two ducted fan systems with coaxial rotors, two sets of control rudders, a 3-DOF manipulator, the control unit, and the landing gears, as shown in Figure3. The pitch channel of the vehicle is controlled by the thrust difference between the front and the rear ducted fans, which is caused by the difference of rotor speed; the roll channel is controlled by the rudder systems setting below the duct; the yaw channel is controlled by the reaction torque difference between the upper and lower rotor disc. The total mass of the aerial vehicle is 4.6 kg, with an effective payload of 2 kg. The manipulator weight is 0.9 kg, with a max grasping weight of 0.5 kg. The parameters of the system are detailed in Table 1, where the mass and the structural dimension parameters of the platform are measured directly, and the moments of inertia and the effective resistance areas along three axes of the vehicle are estimated by the CATIA® 3D model. The airfoil used in the prototype is NACA 0012, and the airfoil parameters can be obtained from [28]. The manipulator is installed on the center of the vehicle body and has 3 DOF (one lumbar joint, one shoulder joint, and one elbow joint) to reach any position in 3D space. The end effector is a gripper, and, from the perspective of our application purpose, the posture of the gripper is not considered. The parameters of the arm are described using standard D-H (Denavit–Hartenberg) method [29], as shown in Table 2.Figure 3
The configuration of the novel aerial manipulator.Table 1
Structural parameters of the system.
ParameterPhysical descriptionValuembVehicle mass (including battery)4.6 kgmrManipulator mass0.9 kgmgraspMax grasping mass0.5 kgIxxInertia tensor of the vehicle aroundx-axis0.092 kg·m2IyyInertia tensor of the vehicle aroundy-axis0.283 kg·m2IzzInertia tensor of the vehicle aroundz-axis0.245 kg·m2pcdDistance between duct center and CG of vehicle0.32 mDDuct diameter0.33 mSxVehicle effective resistance area aroundx-axis0.034 m2SyVehicle effective resistance area aroundy-axis0.127 m2SzVehicle effective resistance area aroundz-axis0.133 m2nBlade number of each disc4cBlade chord length0.027 mθ0Attack angle at the root of blade35 degθrwTorsion rate of blade18 degσBlade tip clearance0.001 mTable 2
D-H parameters of the manipulator.
Linka (m)α (deg)d (m)θ range (deg)10900.08180 [−90, 90]20.1500150 [0, 150]30.1600150 [−150, 0]The hardware of control unit is shown in Figure4. The on-board controller is Emlid® Navio2 based on Raspberry Pi® 3, which integrates dual IMU module, GPS module, barometer module, and 14 PWM output channels. The rotors are driven by four GARTT® motors, and the motors are driven by HOBBYWING® Electronic Speed Controllers (ESC), which are controlled by the main controller via PWM signals. Each duct system has two sets of control rudders, which are driven synchronously by one KST® servo. The servos are also controlled by PWM signals from the main controller. The Dynamixel® XH430 servo is chosen as arm joint servo, since it has both position control mode and torque control mode and has feedback function of actual position, velocity, and torque of the joint. The on-board main controller communicates with the ground station through the 3DR® radio telemetry.Figure 4
The hardware of control unit.


## 2.2. Dynamic Model

The coordinate system of the aerial manipulator is introduced in Figure5. The aerial manipulator is modeled as a multibody system consisting of four interconnected rigid bodies. Let Σ be the earth-fixed Cartesian coordinate frame following the north-east-down rules, Σb be the vehicle body-fixed coordinate frame at the center of mass of the vehicle, Σ0 be the manipulator base-fixed coordinate frame, and Σi (i = 1, 2, 3) be the coordinate frame of each link of the manipulator following the D-H rules. Notice that Σ0 coincides with the origin of Σb, only rotated 90° around the Zb axis. In the following, the superscript i means that the variable is related to the coordinate system Σi.Figure 5
The coordinate system of the aerial manipulator.The dynamic formulation is established using iterative Newton–Euler method [29], which is widely used in dynamic modeling [30, 31]. First, give the definition of the motion states of the system. Let pb=x⋅y⋅zT and Φb=φ⋅θ⋅ψT be the position vector and Euler angle vector of the aerial vehicle platform in earth-fixed frame, and let qi⋅i=1⋅2⋅3 be the joint angle of the manipulator. Let vb and ωb be the velocity vector and angular rate vector of the aerial vehicle, and let vi and ωi be the velocity and angular rate of the origin of link i of the manipulator. Then, the acceleration of the multibody system can be calculated outward iteratively from the aerial vehicle platform to the manipulator end effector as(1)vbb=uvwT=p˙bb=RbTp˙b,ωbb=pqrT=RbTωb=RbTΤbΦ˙b=QbΦ˙b,v00=Rb0vbb,ω00=Rb0ωbb,ωii=Ri−1iωi−1i−1+q˙iezi,ω˙ii=Ri−1iω˙i−1i−1+Ri−1iωi−1i−1×q˙iezi+q¨iezi,v˙ii=Ri−1iω˙i−1i−1×pii−1+ωi−1i−1×ωi−1i−1×pii−1+v˙i−1i−1,aCii=ω˙ii×pCii+ωii×ωii×pCii+v˙ii,where the superscript b refers to the body-fixed coordinate system, the superscript i refers to the ith link coordinate system, and in particular i = 0 for the manipulator base-fixed (link 0) coordinate system, as shown in Figure 5. Specifically, vbb and ωbb are the velocity and angular rate of the aerial vehicle with respect to the body-fixed coordinate frame, vii and ωii are the velocity and angular rate of the origin of link i with respect to the ith link coordinate frame. aCii is the acceleration of the center of mass of link i. pii−1 is the position vector of the origin of ith link frame with respect to the i−1th link coordinate frame, pCii is the position vector of the center of mass of ith link with respect to the ith link coordinate frame, and ezi refers to the projection along the z-axis of link i. Rii−1 refers to the transformation matrix from the i−1th link coordinate frame to the ith link coordinate frame. In particular, Rii−1 refers to the transformation matrix from the body-fixed frame to the arm base-fixed (link 0) frame. Rb and Qb are the linear velocity and angular rate transformation matrix between the earth-fixed frame and body-fixed frame expressed in the form of Euler angles.Then, the Newton–Euler dynamics formulation can be derived as(2)Fbb=mbv˙bb+ωbb×vbb,Mbb=Ibω˙bb+ωbb×Ibωbb,Fii=miaCii,Mii=Iiω˙ii+ωii×Iiωii,where mb and mi are the mass of the aerial vehicle and link i, Ib, and Ii are the inertia matrix of the aerial vehicle and link i. Fbb and Fii are the total external force exerted on the aerial vehicle and link i, and Mbb and Mii are the total external moment exerted on the aerial vehicle and link i, which can be calculated inward iteratively from the manipulator end effector to the aerial vehicle platform as(3)fii=Fii+Ri+1ifi+1i+1−RiTmig,nii=Mii+Ri+1ini+1i+1+pi+1i×Ri+1ifi+1i+1+pCii×Fii−RiTmig,τi=niiTezi,Fbb=Faero1+Faero2+Frudder1+Frudder2+Ffus+RbTmbg−R0bR10f11,Mbb=Maero1+Maero2+pcd1b×Faero1+pcd2b×Faero2+pcr1b×Frudder1+pcr2b×Frudder2+Mgyro−R0bR10n11,where fii and nii are the force and moment exerted on link i by link i−1 and Ri is the transformation matrix from the ith link coordinate frame to the earth-fixed coordinate frame. τi is the joint torque of the manipulator.The total force of the aerial vehicle includes five parts: the aerodynamic force of the ducted fans, the resultant force generated by the rudders, the fuselage resistance, the gravity, and the reaction force of the manipulator. The total moment of the aerial vehicle also includes five parts: the aerodynamic moment of the ducted fans, the torque generated by the aerodynamic forces, the torque generated by the rudders, the reaction torque of the manipulator, and the gyro moment produced by rotors, as shown in Figure6. In equation (3), the subscript “aero” refers to the aerodynamic force and moment, subscript 1 denotes the front ducted fan system, and subscript 2 denotes the rear ducted fan system. Same as the above, the subscript rudder refers to the force of the front and the rear rudder system. fus refers to the fuselage resistance, and gyro refers to the gyro torque. pcdb is the position vector of the center of duct with respect to the center of the vehicle, and pcdb is the position vector of the aerodynamic center of rudder with respect to the center of the vehicle. These forces and moments belong to the ducted fan system dynamics, which provide the main thrust and attitude control moments for the aerial vehicle, and will be introduced in detail in the following section.Figure 6
Diagram of the dynamic model.


## 2.3. Ducted Fan System

As mentioned in the Introduction, most of the literature of aerial manipulator uses multirotor as aerial platform. Because of the simple, symmetrical, and decoupled open rotor structure, the rotor aerodynamics is usually simplified as a thrust and a reaction torque in quadratic relation to the rotor speed, ignoring any other effects of the aerodynamic characteristics (e.g., [22, 26, 27]). However, the aerodynamics of ducted fan is significantly different from the traditional open rotor [32]. First, the duct lip causes airflow deflection effect, which changes the inflow at the rotor disc, while additional thrust is generated by the suction flow of the duct lip; second, duct suppresses the rotor tip vortex and reduces the momentum loss; third, the duct exit prevents the airflow from contracting, which increases the outflow pressure of the rotor and thereby increases the thrust. Johnson and Turbe [17] propose a Blade Element Momentum Theory (BEMT) model for ducted-single-rotor system. By improving this method, the dynamic model of our ducted-coaxial-rotor system is proposed. The basic idea of BEMT is to establish a set of relationships between rotor aerodynamic force and induced velocity through momentum theorem and blade element theory, respectively, and then solve the simultaneous equations iteratively. The following is the derivation of the front ducted fan, and the rear one can be obtained in the same way.Momentum Theorem. The inflow process of the coaxial ducted fan under airflow deflection effect is shown in Figure 7. It can be divided into five parts: the free inflow, the deflected airflow whose direction is deflected by duct lip, the airflow at upper rotor disc, the airflow at lower rotor disc, and the outflow of duct exit. The free flow at the ducted fan is(4)Va=uavawaT=vbb−RbTvw,where vbb is the body velocity of the aerial vehicle and vw is the local wind velocity. Express va in the inflow plane; the inflow velocity v0 is calculated as(5)V0=V0cosαi+V0sinαj,V0=Va,α=arctanwaua2+va2.Figure 7
The inflow process model.Then, the deflected airflow, the airflow at upper and lower rotor disc, and the outflow can be calculated by(6)Vt=V0cosαti+V0sinαtj,Vr1=V0cosαti+V0sinαt+vi1j,Vr2=V0cosαti+V0sinαt+vi1+vi2j,Ve=V0cosαti+V0sinαt+v∞j,where vi1 and vi2 refer to the induced velocities of upper and lower rotor, v∞ refers to the overall impact of the ducted fan on free flow, and αt is the deflection angle caused by duct effect at duct lip, which can be expressed as a function of free inflow angle and airflow deflection factor as(7)αt=α+ktπ2−α.According to momentum theorem and kinetic energy theorem, considering the control body from the free inflow to the upper disc exit and from the upper disc exit to the duct exit and the entire duct, respectively, the following equations are obtained:(8)Tr1V0sinαt+vi1=12ρSDV0V0sinαt+vi12−V0sinα2,Tr2V0sinαt+vi1+vi2=12ρSDV0V0sinαt+v∞2−V0sinαt+vi12,T=1+kductTr1+Tr2=ρSDV0V0sinαt+v∞−V0sinα,Td=ρSDV0V0cosαt−V0cosα,where Tr1 and Tr2 are the thrusts of upper rotor and lower rotor in the vertical direction, Td is the momentum resistance in the horizontal plane, kduct is the thrust augment factor of duct, and SD is the rotor disc area.Blade Element Theory. Blade element theory analyzes the aerodynamics of each micro element of the blade, then integrates it along the radial direction of the blade, and averages it along the circumferential direction to solve the forces and moments of the entire rotor. According to equations (5) and (6), the airflow velocities at the blades of upper and lower rotors are(9)Vra1=ura1vra1wra1T=V0cosαtuaua2+va2V0cosαtvaua2+va2V0sinαtT,Vra2=ura2vra2wra2T=V0cosαtuaua2+va2V0cosαtvaua2+va2V0sinαt+vi1T.Since the blade element theory is a common method, the results after integration of the rotor thrust Tri, the reaction torque Qri, the pneumatic rolling torque Lri, and the pneumatic pitching torque Mri of the upper and lower rotor are directly given below. Notice that i = 1 for upper rotor and i = 2 for lower rotor:(10)Tri=n2π∫02π∫0RdTridψ0=nρcCaΩi2R3423+μui2+μvi2θ0−1+μui2+μvi2θrw2−λwi,Qri=n2π∫02π∫0RdQridψ0=nρcCaΩi2R4423θ0−12θrw−λwiλwi+121+μui2+μvi2Cd,Lri=n2π∫02π∫0RdLridψ0=nρcCaΩi2R44μvi23θ0−12λwi−12θrw,Mri=n2π∫02π∫0RdMridψ0=nρcCaΩi2R44μui23θ0−12λwi−12θrw.In the above formula,n is the blade number of each disc, R is the rotor radius, ψ0 is the blade azimuth angle, c is the blade chord length, Ca is the lift coefficient slope, Cd is the resistance coefficient, Ωi is the rotor speed, θ0 is the blade attack angle, and θrw is the blade torsion rate. μui, μvi, and λwi denote the dimensionless ratios of forward velocity, lateral velocity, and inflow velocity:(11)μui=uraiΩiR,μvi=vraiΩiR,λw=vii−wraiΩiR.Solve equations (8) and (10) simultaneously; the aero force and moment of the ducted fan in body-fixed frame can be derived as(12)Faero=−uaua2+va2Td−vaua2+va2Td−TT,Maero=Lr1−Lr2Mr1−Mr2Qr1−Qr2T.Rudder Dynamics. Same as the ducted fan system, the following gives the analysis of the front rudder, and the rear one can be obtained in the same way. The rudder is in the outflow of the duct exit, and its attack angle can be calculated by(13)αv=αe+δ,where αe is the angle of the duct outflow velocity Ve, and δ is the rudder control angle. Then the aero-lift and aero-drag of the rudder can be derived as(14)Lv=12ρVe2SvClv,Dv=12ρVe2SvCdv,where ρ is the airflow density, Sv is the rudder area, and Clv and Cdv are the lift coefficient and drag coefficient, which are the functions of αv. The force of the control rudder in body-fixed frame can then be derived. nv is the number of the rudders in one duct:(15)Frudder=0nvLvcosαe+Dvsinαenv−Lvsinαe+DvcosαeT.In summary, the expression of each component in equation (3) is obtained, and the comprehensive nonlinear mechanism model is established.


## 2.4. Model Identification of the Aerial Platform

The comprehensive nonlinear mechanism model established in the previous section can fully reflect the characteristics of the system in all-envelope range and be used for system performance prediction and comprehensive simulation analysis, but the complex mechanism is not suitable for the controller design process [33]. Since the working scenarios of the aerial vehicle in this paper are hovering or near-hovering conditions, in order to reduce the complexity of nonlinear identification and improve the identification accuracy under the noise of low-cost sensors, the parametric frequency domain identification is performed. First, a simplified parametric model can be derived from the nonlinear mechanism model of the aerial platform by small perturbation assumption under hover equilibrium point; then, a control-oriented state space model of the aerial platform is derived based on frequency domain identification method. The identification model is verified in both time and frequency domains.The nonlinear mechanism model of the aerial vehicle can be described as(16)x˙b=fxb,ub,xb=uvwpqrφθψT,ub=ualturolupituyawT,where ub is the normalization control input vector, referring to the altitude channel, roll channel, pitch channel, and yaw channel, and is mapped to the rotor speed difference and rudder angle (described in 2.1). Find the partial derivative of equation (16) at hover equilibrium point to obtain the system’s Jacobian matrix:(17)x˙b≈Axb−xhover+Bub−uhover,A=∂f∂xb,B=∂f∂ub,fxhover,uhover=0.Then, the simplified parametric matrix A and B to be identified is obtained by small perturbation assumption. Since xhover and uhover are constant and xhover = 0, in the rest of this paper, for the simplicity of the symbol, let xb ≡ xb − xhover, ub ≡ ub − uhover. Because the aerial vehicle platform is inherently unstable, the closed-loop identification is carried out following the procedure in [34]. A simple PID controller is used to ensure the stability during identification process. A set of Chirp sweep signals (introduced in [35]) is performed in each of the input channels, and the data is processed by CIFER® identification tool [35]. The model fitting process is an optimization problem as(18)minJ=∑i=1n20nω∑ω1ωnωWγWgGfi−Gi2+Wp∠Gfi−∠Gi2,where Gfi is the transfer function to be fitted corresponding to Gi, converting from the parametric state space equation, and n is the number of the transfer functions. |G| refers to the amplitude and ∠G refers to the phase. nω is the number of the selected frequency sampling points, and ω1 and ωnω refer to the initial and termination sampling frequency respectively. Wγ, Wg, and Wp are the weighting functions referring to coherence value, amplitude, and phase, respectively, and are chosen following the directions in CIFER®. The frequency response identification results are shown in Figures 8 and 9, and the identified model is(19)A=−0.0876000000−9.801000−0.087600009.801000−0.1178−0.1172−1.0415−0.0012−0.0209−0.05090000−0.680100000000.0940000−1.06990.01320000000−0.01220000000100000000010000000001000,B=00000000−20.931200.0547−0.001508.349400.02600−0.002510.0876−0.028100.00160.13373.3981000000000000.Figure 8
Frequency response identification results of linear velocities.
(a)(b)(c)(d)(e)(f)Figure 9
Frequency response identification results of angular rates.
(a)(b)(c)(d)(e)(f)In order to verify the accuracy of the identified model, a series of flight tests are carried out in time domain, as shown in Figure10. The low-pass filtered bipolar square wave signal, which is different from the identification excitation signal, is used as the excitation signal in each channel, and the evaluation of the model is realized by comparing the output of the identified model simulation and the actual test measurement. Notice that the angular rate errors in Figure 10 are large while other variables match well, because the dynamic characteristics of angular rate are ignored in high frequency in the linearized model and the angular rate channel has relatively large measurement noise. According to [35], Theil inequality coefficient (TIC) can be used to estimate the model fitting accuracy, which is calculated as(20)JTIC=1/n∑i=1nytest−ysimTytest−ysim1/n∑i=1nysimTysim+1/n∑i=1nytestTytest,where ysim and ytest are the output vector of identified model simulation and actual flight test, and n is the number of sampling points. A smaller value of TIC indicates that the accuracy of the identified model is higher, and a value less than 0.25 indicates that the model has good predictive performance [35]. JTIC of our model is 0.1596, indicating that the nominal model has satisfactory accuracy.Figure 10
Time-domain verification of the identified model.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3. Composite Controller Design

The overall architecture of the aerial manipulator control system is shown in Figure11. A typical application scenario of the aerial manipulator is that, given the reference position of the target object, the aerial vehicle flies near the target, and the manipulator moves, grabs the object, and then leaves. Since the process of arriving at the desired position of the aerial platform and completing the operation task of the manipulator is independent, and the control objectives, dynamic characteristics, and response speed of the aerial platform and the manipulator are different, the controllers of the aerial platform and the manipulator are designed separately for the purposes of reducing the controller complexity, reducing the hardware cost, and facilitating the engineering practice. In this framework, the most important challenge is the aerial platform controller design under the motion of manipulator and external interaction disturbances of target object.Figure 11
The architecture of the control system.This section mainly focuses on the composite controller design of the aerial platform under the existence of manipulator disturbances. First, a two-layer basic controller is designed to ensure the stability, decoupling, and tracking performance of the platform. On the one hand, the aerial vehicle is a strongly coupled, nonlinear, multi-input multi-output (MIMO) system; thus, traditional controllers (such as PID controller) cannot guarantee good performance of the system. On the other hand, based on the premise of guaranteeing performance, the controller should have low complexity, clear structure, and high solution rate to be easy in practice. Considering these two aspects, a robust static feedback controller based on H-infinite synthesis and nonsmooth optimization is designed as the basic inner loop controller, which guarantees good robustness and stability of the system while not increasing the system order. Second, considering the movement of the manipulator and the disturbances of the environment, an auxiliary adaptive loop is designed to estimate and compensate for the disturbances of the aerial platform. Both the static (gravity force and moment) and dynamic (inertia force and moment) disturbances of the manipulator are considered to improve the platform stability and the end effector tracking accuracy.In addition, in order to improve the control accuracy of the manipulator under the motion of the aerial platform, a computed torque PID controller, considering the dynamics of the aerial platform, is used for manipulator control. Since the computed torque PID control is a conventional method, and the manipulator dynamics is elaborated in Section2, the controller of the manipulator is not detailed in this paper.


### 3.1. Robust Basic Controller Design

The design objective of the basic controller is to achieve the stability and accurate position tracking of the aerial platform. The diagram of the basic controller is shown in Figure12, which includes two loops. The inner loop is an H-infinite static full-state feedback controller Kb to ensure the decoupling and stabilization of the system, and a prepositive static controller Kf is proposed for tracking performance of the inner loop. Since the inner loop has been well decoupled and stabilized, the outer loop (position loop) is designed as a series of simple PD controllers for position tracking. According to (17), the control-oriented dynamic model of the aerial platform is(21)x˙bt=Axbt+Bubt+ft,ybt=Cxbt,xb=uvwpqrφθψT,yb=uvwrT,where f is the disturbance vector. The reference command is rb = [uref vref wref rref]T. Accordingly, the outer loop reference command is pr = [xref yref zrefψref]T. The kinematics module in Figure 12 converts the variables between the body-fixed coordinate frame and earth-fixed coordinate frame according to equation (1).Figure 12
Basic controller architecture.Focusing on the inner loop controller, we can define the output error as(22)eb=rb−yb.Rewrite the system model in frequency domain and augment it, and the structured decoupling standard form of H-infinite synthesis problem [36] can be obtained as in Figure 13:(23)yaug=ybxaugT=Gauguaug,xaug=ebxbT,uaug=frbubT,Gaugω=CGf0CGb−CGfI−CGbGf0Gb,where the control law of the robust controller can be calculated as(24)ub=Kfeb−Kbxb=Kstxaug,Kst=Kf00−Kb.Figure 13
Standard structure of H-infinity synthesis.For the augmented system, the closed-loop transfer function from the inputm to the output n can be calculated by linear fractional transformation (LFT) as(25)Tm⟶ns=FlGaug,Kst.Then, based on H-infinite synthesis method, the controller tuning process can be converted into an optimization problem as(26)minKstJ=WΔTΔGaug,Kst∞subject toKststabilizesGauginternallyWΔTΔGaug,Kst∞<1,where WΔTΔ∞ are the performance constraint functions expressed by H-infinite norm and WΔ are the corresponding weighting functions. The solvability of the optimization problem guarantees the robust stability of the closed-loop system. For the aerial platform in this paper, the following functions are carried out for ensuring tracking performance, disturbance rejection performance, and input energy limitation, respectively:(27)WesTrb⟶ebs∞<1,WbsTrb⟶ybs∞<1,WfsTf⟶ybs∞<1,WusTrb⟶ubs∞<1.Solve equation (26) using nonsmooth optimization algorithm [37]; the controller parameters are as follows. In addition, after the inner loop controller has been solved, the outer loop PD controller can be tuned easily by optimization, and the results are also given here:(28)Kb=0.00560.00550.04980.0012−0.00010.00060.0086−0.00400.00010.0001−0.16480.00021.16580.01110.25149.3447−0.0118−0.00080.0751−0.0001−0.00010.02370.80120.03280.01487.3617−0.0003−0.00500.00780.00020.18880.02572.6449−0.8605−0.51680.0005,Kf=−0.04780.00090.0005−0.0019−0.00340.95420.00090.24730.00070.0049−0.75020.0351−0.0167−0.07000.05132.6541,PD=diag0.85+0.8s0.85+0.8s0.95+1.0s0.95+1.0s.Figure14 illustrates the robust stability margin of the closed-loop system, where the gain margin exceeds 6 dB, and the phase margin exceeds 45°. The results demonstrate that the proposed robust controller can guarantee the basic performance of the system.Figure 14
The robust stability margin with the basic controller.
(a)(b)


### 3.2. Adaptive Auxiliary Controller Design

Although the designed robust controller has achieved basic stability and tracking performance and has antidisturbance ability to some extent, it cannot deal with the aerial manipulator system well because the disturbance caused by manipulator operation far exceeds the stability margin of the robust controller, which makes the system performance deteriorate dramatically. Therefore, it is necessary to design an adaptive auxiliary loop to compensate for the impact of the manipulator. The architecture of the adaptive controller is shown in Figure15. The input of the adaptive loop is the output of the outer loop basic PD controller, and its output is the reference input of the inner closed-loop system. Because the inner loop robust controller is static, it does not increase or decrease the state variables and does not change the physical meanings of the system states; thus, the inner closed-loop system can be obtained easily from Section 3.1 as(29)x˙bt=Aadxbt+Baduadt+ft,ybt=Cxbt,which is clearly bounded-input bounded-state (BIBS) stable relying on the condition of (26) in the inner loop controller design process, and Aad is Hurwitz matrix. Besides, according to the previous description, all the state variables are measurable and xb0=0. f(t) is the time-varying disturbance, which is bounded and differentiable, and also the differential is bounded. The boundary of f(t) can be estimated by the manipulator dynamics and max grasping payload introduced in Section 2, with some redundancy on this basis. Thus, the following condition is given:(30)∃compact setF,∀t≥0,ft∈F,∃FanddF,∀t≥0,ft2≤F,ft2≤dF,Figure 15
Adaptive controller architecture.The adaptive loop consists of three parts: first, a state predictor is designed for state prediction; then, the adaptive law is designed for disturbance estimation based on the prediction error; finally, the control law is reconstructed to compensate for the disturbance and track the desired reference. In particular, inspired by the theory in [38], a filter matrix C(s) is applied to shape the control input, which decouples the adaptive rate from system robustness. Thus, the large adaptive gain can be adopted while avoiding adverse effects on the system states, which ensures the accurate estimation of disturbance and guarantees good transient and steady-state performance of the adaptive controller.


#### 3.2.1. Design of the Adaptive Auxiliary Controller

Consider the following state predictor of the system in equation (29):(31)x^˙bt=Aadx^bt+Baduadt+f^t,y^bt=Cx^bt,where the cap ^ denotes the estimation of the corresponding variable. Then, the estimation error dynamics can be calculated from (29) and (31) as(32)x˜˙bt=Aadx˜bt+f˜t,x˜bt=x^bt−xbt,f˜t=f^t−ft.In order to ensure the stability of estimation error, consider the following candidate Lyapunov function:(33)Vx˜bt,f˜t=x˜bTtPx˜bt+Γ−1f˜Ttf˜t,where Γ is the adaptive gain. Since Aad is Hurwitz, P is the solution of the following algebraic Lyapunov equation:(34)AadTP+PAad=−I.The derivation of (33) is(35)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t.Based on Lyapunov stability criterion, let the adaptive law be(36)f^˙t=ΓProjf^t,−x˜bTtPT,where Proj denotes the projection operator, which prevents the divergence and ensures the boundedness of the adaptive parameter, detailed in [39]. Thus, (35) can be further derived as(37)V˙t=−x˜bTtx˜bt≤0.Notice that in the adaptive law derivation in (35), it is assumed that the change rate of the disturbance f is slower than the convergence rate of the adaptive estimation; thus, its derivative is approximately 0. In the performance analysis of the next section, this assumption will be relaxed and more general situations will be considered. The estimation error stability under arbitrary disturbance including arbitrary fast time-varying disturbance will be proved in detail with the proposed adaptive law in (36).Then the control law in Laplace form can be calculated as follows (refer to Figure15):(38)uads=CsKgrbs−F1sf^s.where Kg=−CAad−1Bad−1,F1s=H1−1sCH0s,andH0s=sI−Aad−1,Hs=H0sBad,H1s=CHs. C(s) is strictly proper and stable low-pass filter matrix. Following the idea in [38], C(s) could decouple the adaptive rate from system robustness, and the selection of it should satisfy the following conditions:(39)Csis  strictly proper and stable withC0=I,∃ρr1,CsH1−1sL1<ρr1,∃ρr2,H0sI−BadCsF1sL1<ρr2,where ΔsL1 refers to the L1 norm of Δs and the last two conditions indicate that the corresponding L1 norm gain is finite.Thus, the design process of the adaptive controller is completed. To sum up, for the system in (29), choose the appropriate adaptive gain Γ and filter matrix C(s) with satisfying condition (39); the control input is calculated via (31), (36), and (38) as(40)x^˙bt=Aadx^bt+Baduadt+f^t,f^˙t=ΓProjf^t,−x˜bTtPT,uads=CsKgrbs−F1sf^s.


#### 3.2.2. Performance Analysis of the Adaptive Auxiliary Controller

Consider the ideal form of the adaptive controller, which is assumed to compensate for the disturbance perfectly. The closed-loop reference system is defined as(41)x˙reft=Aadxreft+Badureft+ft,yreft=Cxreft,urefs=CsKgrbs−F1sfs,and xref0=0. Using Laplace transformation, the closed-loop system (41) can be described as(42)xrefs=HsCsKgrbs+H0sI−BadCsF1sfs.The following lemma describes the L1 norm stability theorem [38].Lemma 1.
A continuous-time linear time-invariant (LTI) MIMO system y(s) = G(s)u(s) is BIBO stable if its L1 norm is bounded and for arbitrary bounded u(t) one has(43)yL∞≤GsL1uL∞,where ΔsL1 refers to the L1 norm of Δs, and ΔsL∞ refers to the L∞ norm of Δs.SinceAad is Hurwitz and the conditions in (30) and (39) are satisfied, for all t ≥ 0, the L1 norm of system (42) is finite; thus, ||xref||L∞ is uniformly bounded, the closed-loop reference system (41) is BIBS stable, and the reference system has the following performance bound:(44)xrefL∞≤HsCsKgL1rbL∞+H0sI−BadCsF1sL1fL∞.For the defined state predictor (31), consider the candidate Lyapunov function (33) again; its derivation under arbitrary fast time-varying disturbance with the adaptive law (40) can be derived as(45)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t−2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Projf^t,−x˜bTtPT+2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+4Γ−1dFF.Also, from (30) and (33) one has(46)maxt≥0f˜Ttf˜t≤4maxf∈Ff22=4F2,V0=Γ−1f˜T0f˜0≤4Γ−1F2.Let Ωm: = 4λmax(P)dFF + 4F2; if there exist some t > 0 such that V(t) > Γ−1Ωm, then, it follows from (33) and (46) that(47)x˜bTtPx˜bt>4Γ−1λmaxPdFF,and thus,(48)x˜bTtx˜bt≥λmax−1Px˜bTtPx˜bt>4Γ−1dFF.Then, from equations (45) and (48), one can deduce that(49)V˙t<0.According to equation (46), it is derived that V(0) < Γ−1Ωm. Based on the deduction above, it can be obtained that(50)∀t≥0,Vt≤Γ−1Ωm.Then, the bound of the state estimation error can be derived as(51)λminPx˜bt22≤x˜bTtPx˜bt≤Vt≤Γ−1Ωm⇒x˜bt2≤λmin−1PΓ−1Ωm,which proves that the estimation error is uniformly bounded and the state predictor is Lyapunov stable. It also can be drawn from equation (51) that the performance of the state predictor is inversely proportional to the square root of the adaptive gain Γ, which can be improved by increasing the adaptive gain.It follows from equations (29) and (40) that(52)xbs=HsCsKgrbs+H0s−I−BadCsF1sfs−HsCsF1sf˜s,uads=CsKgrbs−F1sfs−F1sf˜s.Consider equations (32), (42), and (52) simultaneously; the adaptive state error can be calculated as(53)xrefs−xbs=HsCsH1−1sCx˜bs.According to the condition in equation (39) and the bound in equation (51), for all t ≥ 0, the adaptive state error is uniformly bounded:(54)xref−xbL∞≤HsCsH1−1sCL1x˜bL∞≤HsCsH1−1sCL1λmin−1PΓ−1Ωm.The adaptive input error can be calculated by equations (41) and (52), and in a similar way to the state error, the bound of adaptive input error for all t ≥ 0 can be derived as(55)uref−uadL∞≤CsH1−1sCL1x˜bL∞≤CsH1−1sCL1λmin−1PΓ−1Ωm.Since the closed-loop reference system proved to be BIBS stable in (44), combined with (54) and (55), it can be derived that the closed-loop system with the proposed adaptive controller is BIBS stable, and the performance bound relies on the adaptive gain Γ.Because the system robustness is decoupled with the adaptive gain by the filter matrix C(s), a large Γ can be applied to guarantee good transient and steady-state performance. Theoretically, the choice of adaptive gain is limited only by the performance of the on-board processor hardware. In some literature works, the disturbance f is called unmatched uncertainty, and, for a general open-loop MIMO system, only its effect on the system output (not the system state) can be compensated for by the control input. However, our system for adaptive controller design is BIBS stable and the system states are equal to system outputs (C is an identity matrix), which are measurable with actual physical meanings, and the disturbance boundary can be obtained by the manipulator dynamics. Moreover, thanks to the filter matrix C(s) and the large adaptive gain, the good performance of the adaptive system can be achieved while not leading to adverse effect on the system states. Therefore, the disturbance f of our system can be estimated and compensated for well. Based on the above discussion, the parameters of the adaptive controller are given as(56)Cs=25s2+6s+25I,Γ=1000.


## 3.1. Robust Basic Controller Design

The design objective of the basic controller is to achieve the stability and accurate position tracking of the aerial platform. The diagram of the basic controller is shown in Figure12, which includes two loops. The inner loop is an H-infinite static full-state feedback controller Kb to ensure the decoupling and stabilization of the system, and a prepositive static controller Kf is proposed for tracking performance of the inner loop. Since the inner loop has been well decoupled and stabilized, the outer loop (position loop) is designed as a series of simple PD controllers for position tracking. According to (17), the control-oriented dynamic model of the aerial platform is(21)x˙bt=Axbt+Bubt+ft,ybt=Cxbt,xb=uvwpqrφθψT,yb=uvwrT,where f is the disturbance vector. The reference command is rb = [uref vref wref rref]T. Accordingly, the outer loop reference command is pr = [xref yref zrefψref]T. The kinematics module in Figure 12 converts the variables between the body-fixed coordinate frame and earth-fixed coordinate frame according to equation (1).Figure 12
Basic controller architecture.Focusing on the inner loop controller, we can define the output error as(22)eb=rb−yb.Rewrite the system model in frequency domain and augment it, and the structured decoupling standard form of H-infinite synthesis problem [36] can be obtained as in Figure 13:(23)yaug=ybxaugT=Gauguaug,xaug=ebxbT,uaug=frbubT,Gaugω=CGf0CGb−CGfI−CGbGf0Gb,where the control law of the robust controller can be calculated as(24)ub=Kfeb−Kbxb=Kstxaug,Kst=Kf00−Kb.Figure 13
Standard structure of H-infinity synthesis.For the augmented system, the closed-loop transfer function from the inputm to the output n can be calculated by linear fractional transformation (LFT) as(25)Tm⟶ns=FlGaug,Kst.Then, based on H-infinite synthesis method, the controller tuning process can be converted into an optimization problem as(26)minKstJ=WΔTΔGaug,Kst∞subject toKststabilizesGauginternallyWΔTΔGaug,Kst∞<1,where WΔTΔ∞ are the performance constraint functions expressed by H-infinite norm and WΔ are the corresponding weighting functions. The solvability of the optimization problem guarantees the robust stability of the closed-loop system. For the aerial platform in this paper, the following functions are carried out for ensuring tracking performance, disturbance rejection performance, and input energy limitation, respectively:(27)WesTrb⟶ebs∞<1,WbsTrb⟶ybs∞<1,WfsTf⟶ybs∞<1,WusTrb⟶ubs∞<1.Solve equation (26) using nonsmooth optimization algorithm [37]; the controller parameters are as follows. In addition, after the inner loop controller has been solved, the outer loop PD controller can be tuned easily by optimization, and the results are also given here:(28)Kb=0.00560.00550.04980.0012−0.00010.00060.0086−0.00400.00010.0001−0.16480.00021.16580.01110.25149.3447−0.0118−0.00080.0751−0.0001−0.00010.02370.80120.03280.01487.3617−0.0003−0.00500.00780.00020.18880.02572.6449−0.8605−0.51680.0005,Kf=−0.04780.00090.0005−0.0019−0.00340.95420.00090.24730.00070.0049−0.75020.0351−0.0167−0.07000.05132.6541,PD=diag0.85+0.8s0.85+0.8s0.95+1.0s0.95+1.0s.Figure14 illustrates the robust stability margin of the closed-loop system, where the gain margin exceeds 6 dB, and the phase margin exceeds 45°. The results demonstrate that the proposed robust controller can guarantee the basic performance of the system.Figure 14
The robust stability margin with the basic controller.
(a)(b)


## 3.2. Adaptive Auxiliary Controller Design

Although the designed robust controller has achieved basic stability and tracking performance and has antidisturbance ability to some extent, it cannot deal with the aerial manipulator system well because the disturbance caused by manipulator operation far exceeds the stability margin of the robust controller, which makes the system performance deteriorate dramatically. Therefore, it is necessary to design an adaptive auxiliary loop to compensate for the impact of the manipulator. The architecture of the adaptive controller is shown in Figure15. The input of the adaptive loop is the output of the outer loop basic PD controller, and its output is the reference input of the inner closed-loop system. Because the inner loop robust controller is static, it does not increase or decrease the state variables and does not change the physical meanings of the system states; thus, the inner closed-loop system can be obtained easily from Section 3.1 as(29)x˙bt=Aadxbt+Baduadt+ft,ybt=Cxbt,which is clearly bounded-input bounded-state (BIBS) stable relying on the condition of (26) in the inner loop controller design process, and Aad is Hurwitz matrix. Besides, according to the previous description, all the state variables are measurable and xb0=0. f(t) is the time-varying disturbance, which is bounded and differentiable, and also the differential is bounded. The boundary of f(t) can be estimated by the manipulator dynamics and max grasping payload introduced in Section 2, with some redundancy on this basis. Thus, the following condition is given:(30)∃compact setF,∀t≥0,ft∈F,∃FanddF,∀t≥0,ft2≤F,ft2≤dF,Figure 15
Adaptive controller architecture.The adaptive loop consists of three parts: first, a state predictor is designed for state prediction; then, the adaptive law is designed for disturbance estimation based on the prediction error; finally, the control law is reconstructed to compensate for the disturbance and track the desired reference. In particular, inspired by the theory in [38], a filter matrix C(s) is applied to shape the control input, which decouples the adaptive rate from system robustness. Thus, the large adaptive gain can be adopted while avoiding adverse effects on the system states, which ensures the accurate estimation of disturbance and guarantees good transient and steady-state performance of the adaptive controller.


### 3.2.1. Design of the Adaptive Auxiliary Controller

Consider the following state predictor of the system in equation (29):(31)x^˙bt=Aadx^bt+Baduadt+f^t,y^bt=Cx^bt,where the cap ^ denotes the estimation of the corresponding variable. Then, the estimation error dynamics can be calculated from (29) and (31) as(32)x˜˙bt=Aadx˜bt+f˜t,x˜bt=x^bt−xbt,f˜t=f^t−ft.In order to ensure the stability of estimation error, consider the following candidate Lyapunov function:(33)Vx˜bt,f˜t=x˜bTtPx˜bt+Γ−1f˜Ttf˜t,where Γ is the adaptive gain. Since Aad is Hurwitz, P is the solution of the following algebraic Lyapunov equation:(34)AadTP+PAad=−I.The derivation of (33) is(35)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t.Based on Lyapunov stability criterion, let the adaptive law be(36)f^˙t=ΓProjf^t,−x˜bTtPT,where Proj denotes the projection operator, which prevents the divergence and ensures the boundedness of the adaptive parameter, detailed in [39]. Thus, (35) can be further derived as(37)V˙t=−x˜bTtx˜bt≤0.Notice that in the adaptive law derivation in (35), it is assumed that the change rate of the disturbance f is slower than the convergence rate of the adaptive estimation; thus, its derivative is approximately 0. In the performance analysis of the next section, this assumption will be relaxed and more general situations will be considered. The estimation error stability under arbitrary disturbance including arbitrary fast time-varying disturbance will be proved in detail with the proposed adaptive law in (36).Then the control law in Laplace form can be calculated as follows (refer to Figure15):(38)uads=CsKgrbs−F1sf^s.where Kg=−CAad−1Bad−1,F1s=H1−1sCH0s,andH0s=sI−Aad−1,Hs=H0sBad,H1s=CHs. C(s) is strictly proper and stable low-pass filter matrix. Following the idea in [38], C(s) could decouple the adaptive rate from system robustness, and the selection of it should satisfy the following conditions:(39)Csis  strictly proper and stable withC0=I,∃ρr1,CsH1−1sL1<ρr1,∃ρr2,H0sI−BadCsF1sL1<ρr2,where ΔsL1 refers to the L1 norm of Δs and the last two conditions indicate that the corresponding L1 norm gain is finite.Thus, the design process of the adaptive controller is completed. To sum up, for the system in (29), choose the appropriate adaptive gain Γ and filter matrix C(s) with satisfying condition (39); the control input is calculated via (31), (36), and (38) as(40)x^˙bt=Aadx^bt+Baduadt+f^t,f^˙t=ΓProjf^t,−x˜bTtPT,uads=CsKgrbs−F1sf^s.


### 3.2.2. Performance Analysis of the Adaptive Auxiliary Controller

Consider the ideal form of the adaptive controller, which is assumed to compensate for the disturbance perfectly. The closed-loop reference system is defined as(41)x˙reft=Aadxreft+Badureft+ft,yreft=Cxreft,urefs=CsKgrbs−F1sfs,and xref0=0. Using Laplace transformation, the closed-loop system (41) can be described as(42)xrefs=HsCsKgrbs+H0sI−BadCsF1sfs.The following lemma describes the L1 norm stability theorem [38].Lemma 1.
A continuous-time linear time-invariant (LTI) MIMO system y(s) = G(s)u(s) is BIBO stable if its L1 norm is bounded and for arbitrary bounded u(t) one has(43)yL∞≤GsL1uL∞,where ΔsL1 refers to the L1 norm of Δs, and ΔsL∞ refers to the L∞ norm of Δs.SinceAad is Hurwitz and the conditions in (30) and (39) are satisfied, for all t ≥ 0, the L1 norm of system (42) is finite; thus, ||xref||L∞ is uniformly bounded, the closed-loop reference system (41) is BIBS stable, and the reference system has the following performance bound:(44)xrefL∞≤HsCsKgL1rbL∞+H0sI−BadCsF1sL1fL∞.For the defined state predictor (31), consider the candidate Lyapunov function (33) again; its derivation under arbitrary fast time-varying disturbance with the adaptive law (40) can be derived as(45)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t−2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Projf^t,−x˜bTtPT+2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+4Γ−1dFF.Also, from (30) and (33) one has(46)maxt≥0f˜Ttf˜t≤4maxf∈Ff22=4F2,V0=Γ−1f˜T0f˜0≤4Γ−1F2.Let Ωm: = 4λmax(P)dFF + 4F2; if there exist some t > 0 such that V(t) > Γ−1Ωm, then, it follows from (33) and (46) that(47)x˜bTtPx˜bt>4Γ−1λmaxPdFF,and thus,(48)x˜bTtx˜bt≥λmax−1Px˜bTtPx˜bt>4Γ−1dFF.Then, from equations (45) and (48), one can deduce that(49)V˙t<0.According to equation (46), it is derived that V(0) < Γ−1Ωm. Based on the deduction above, it can be obtained that(50)∀t≥0,Vt≤Γ−1Ωm.Then, the bound of the state estimation error can be derived as(51)λminPx˜bt22≤x˜bTtPx˜bt≤Vt≤Γ−1Ωm⇒x˜bt2≤λmin−1PΓ−1Ωm,which proves that the estimation error is uniformly bounded and the state predictor is Lyapunov stable. It also can be drawn from equation (51) that the performance of the state predictor is inversely proportional to the square root of the adaptive gain Γ, which can be improved by increasing the adaptive gain.It follows from equations (29) and (40) that(52)xbs=HsCsKgrbs+H0s−I−BadCsF1sfs−HsCsF1sf˜s,uads=CsKgrbs−F1sfs−F1sf˜s.Consider equations (32), (42), and (52) simultaneously; the adaptive state error can be calculated as(53)xrefs−xbs=HsCsH1−1sCx˜bs.According to the condition in equation (39) and the bound in equation (51), for all t ≥ 0, the adaptive state error is uniformly bounded:(54)xref−xbL∞≤HsCsH1−1sCL1x˜bL∞≤HsCsH1−1sCL1λmin−1PΓ−1Ωm.The adaptive input error can be calculated by equations (41) and (52), and in a similar way to the state error, the bound of adaptive input error for all t ≥ 0 can be derived as(55)uref−uadL∞≤CsH1−1sCL1x˜bL∞≤CsH1−1sCL1λmin−1PΓ−1Ωm.Since the closed-loop reference system proved to be BIBS stable in (44), combined with (54) and (55), it can be derived that the closed-loop system with the proposed adaptive controller is BIBS stable, and the performance bound relies on the adaptive gain Γ.Because the system robustness is decoupled with the adaptive gain by the filter matrix C(s), a large Γ can be applied to guarantee good transient and steady-state performance. Theoretically, the choice of adaptive gain is limited only by the performance of the on-board processor hardware. In some literature works, the disturbance f is called unmatched uncertainty, and, for a general open-loop MIMO system, only its effect on the system output (not the system state) can be compensated for by the control input. However, our system for adaptive controller design is BIBS stable and the system states are equal to system outputs (C is an identity matrix), which are measurable with actual physical meanings, and the disturbance boundary can be obtained by the manipulator dynamics. Moreover, thanks to the filter matrix C(s) and the large adaptive gain, the good performance of the adaptive system can be achieved while not leading to adverse effect on the system states. Therefore, the disturbance f of our system can be estimated and compensated for well. Based on the above discussion, the parameters of the adaptive controller are given as(56)Cs=25s2+6s+25I,Γ=1000.


## 3.2.1. Design of the Adaptive Auxiliary Controller

Consider the following state predictor of the system in equation (29):(31)x^˙bt=Aadx^bt+Baduadt+f^t,y^bt=Cx^bt,where the cap ^ denotes the estimation of the corresponding variable. Then, the estimation error dynamics can be calculated from (29) and (31) as(32)x˜˙bt=Aadx˜bt+f˜t,x˜bt=x^bt−xbt,f˜t=f^t−ft.In order to ensure the stability of estimation error, consider the following candidate Lyapunov function:(33)Vx˜bt,f˜t=x˜bTtPx˜bt+Γ−1f˜Ttf˜t,where Γ is the adaptive gain. Since Aad is Hurwitz, P is the solution of the following algebraic Lyapunov equation:(34)AadTP+PAad=−I.The derivation of (33) is(35)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t.Based on Lyapunov stability criterion, let the adaptive law be(36)f^˙t=ΓProjf^t,−x˜bTtPT,where Proj denotes the projection operator, which prevents the divergence and ensures the boundedness of the adaptive parameter, detailed in [39]. Thus, (35) can be further derived as(37)V˙t=−x˜bTtx˜bt≤0.Notice that in the adaptive law derivation in (35), it is assumed that the change rate of the disturbance f is slower than the convergence rate of the adaptive estimation; thus, its derivative is approximately 0. In the performance analysis of the next section, this assumption will be relaxed and more general situations will be considered. The estimation error stability under arbitrary disturbance including arbitrary fast time-varying disturbance will be proved in detail with the proposed adaptive law in (36).Then the control law in Laplace form can be calculated as follows (refer to Figure15):(38)uads=CsKgrbs−F1sf^s.where Kg=−CAad−1Bad−1,F1s=H1−1sCH0s,andH0s=sI−Aad−1,Hs=H0sBad,H1s=CHs. C(s) is strictly proper and stable low-pass filter matrix. Following the idea in [38], C(s) could decouple the adaptive rate from system robustness, and the selection of it should satisfy the following conditions:(39)Csis  strictly proper and stable withC0=I,∃ρr1,CsH1−1sL1<ρr1,∃ρr2,H0sI−BadCsF1sL1<ρr2,where ΔsL1 refers to the L1 norm of Δs and the last two conditions indicate that the corresponding L1 norm gain is finite.Thus, the design process of the adaptive controller is completed. To sum up, for the system in (29), choose the appropriate adaptive gain Γ and filter matrix C(s) with satisfying condition (39); the control input is calculated via (31), (36), and (38) as(40)x^˙bt=Aadx^bt+Baduadt+f^t,f^˙t=ΓProjf^t,−x˜bTtPT,uads=CsKgrbs−F1sf^s.


## 3.2.2. Performance Analysis of the Adaptive Auxiliary Controller

Consider the ideal form of the adaptive controller, which is assumed to compensate for the disturbance perfectly. The closed-loop reference system is defined as(41)x˙reft=Aadxreft+Badureft+ft,yreft=Cxreft,urefs=CsKgrbs−F1sfs,and xref0=0. Using Laplace transformation, the closed-loop system (41) can be described as(42)xrefs=HsCsKgrbs+H0sI−BadCsF1sfs.The following lemma describes the L1 norm stability theorem [38].Lemma 1.
A continuous-time linear time-invariant (LTI) MIMO system y(s) = G(s)u(s) is BIBO stable if its L1 norm is bounded and for arbitrary bounded u(t) one has(43)yL∞≤GsL1uL∞,where ΔsL1 refers to the L1 norm of Δs, and ΔsL∞ refers to the L∞ norm of Δs.SinceAad is Hurwitz and the conditions in (30) and (39) are satisfied, for all t ≥ 0, the L1 norm of system (42) is finite; thus, ||xref||L∞ is uniformly bounded, the closed-loop reference system (41) is BIBS stable, and the reference system has the following performance bound:(44)xrefL∞≤HsCsKgL1rbL∞+H0sI−BadCsF1sL1fL∞.For the defined state predictor (31), consider the candidate Lyapunov function (33) again; its derivation under arbitrary fast time-varying disturbance with the adaptive law (40) can be derived as(45)V˙t=x˜˙bTtPx˜bt+x˜bTtPx˜˙bt+Γ−1f˜˙Ttf˜t+f˜Ttf˜˙t=−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Γ−1f^˙t−2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+2f˜Ttx˜bTtPT+Projf^t,−x˜bTtPT+2Γ−1f˜Ttf˙t≤−x˜bTtx˜bt+4Γ−1dFF.Also, from (30) and (33) one has(46)maxt≥0f˜Ttf˜t≤4maxf∈Ff22=4F2,V0=Γ−1f˜T0f˜0≤4Γ−1F2.Let Ωm: = 4λmax(P)dFF + 4F2; if there exist some t > 0 such that V(t) > Γ−1Ωm, then, it follows from (33) and (46) that(47)x˜bTtPx˜bt>4Γ−1λmaxPdFF,and thus,(48)x˜bTtx˜bt≥λmax−1Px˜bTtPx˜bt>4Γ−1dFF.Then, from equations (45) and (48), one can deduce that(49)V˙t<0.According to equation (46), it is derived that V(0) < Γ−1Ωm. Based on the deduction above, it can be obtained that(50)∀t≥0,Vt≤Γ−1Ωm.Then, the bound of the state estimation error can be derived as(51)λminPx˜bt22≤x˜bTtPx˜bt≤Vt≤Γ−1Ωm⇒x˜bt2≤λmin−1PΓ−1Ωm,which proves that the estimation error is uniformly bounded and the state predictor is Lyapunov stable. It also can be drawn from equation (51) that the performance of the state predictor is inversely proportional to the square root of the adaptive gain Γ, which can be improved by increasing the adaptive gain.It follows from equations (29) and (40) that(52)xbs=HsCsKgrbs+H0s−I−BadCsF1sfs−HsCsF1sf˜s,uads=CsKgrbs−F1sfs−F1sf˜s.Consider equations (32), (42), and (52) simultaneously; the adaptive state error can be calculated as(53)xrefs−xbs=HsCsH1−1sCx˜bs.According to the condition in equation (39) and the bound in equation (51), for all t ≥ 0, the adaptive state error is uniformly bounded:(54)xref−xbL∞≤HsCsH1−1sCL1x˜bL∞≤HsCsH1−1sCL1λmin−1PΓ−1Ωm.The adaptive input error can be calculated by equations (41) and (52), and in a similar way to the state error, the bound of adaptive input error for all t ≥ 0 can be derived as(55)uref−uadL∞≤CsH1−1sCL1x˜bL∞≤CsH1−1sCL1λmin−1PΓ−1Ωm.Since the closed-loop reference system proved to be BIBS stable in (44), combined with (54) and (55), it can be derived that the closed-loop system with the proposed adaptive controller is BIBS stable, and the performance bound relies on the adaptive gain Γ.Because the system robustness is decoupled with the adaptive gain by the filter matrix C(s), a large Γ can be applied to guarantee good transient and steady-state performance. Theoretically, the choice of adaptive gain is limited only by the performance of the on-board processor hardware. In some literature works, the disturbance f is called unmatched uncertainty, and, for a general open-loop MIMO system, only its effect on the system output (not the system state) can be compensated for by the control input. However, our system for adaptive controller design is BIBS stable and the system states are equal to system outputs (C is an identity matrix), which are measurable with actual physical meanings, and the disturbance boundary can be obtained by the manipulator dynamics. Moreover, thanks to the filter matrix C(s) and the large adaptive gain, the good performance of the adaptive system can be achieved while not leading to adverse effect on the system states. Therefore, the disturbance f of our system can be estimated and compensated for well. Based on the above discussion, the parameters of the adaptive controller are given as(56)Cs=25s2+6s+25I,Γ=1000.


## 4. Simulation Verification and Analysis

The aim of the simulation tests is to verify that the proposed aerial manipulator with the composite controller can meet the requirements of high precision operation tasks well. Concretely, they include three aspects: first, verifying the feasibility of the proposed innovative aerial manipulator design with tandem ducted fans; second, verifying the effectiveness of the proposed basic controller for the basic performance; third, comparing the stability performance and tracking performance of the aerial manipulator with and without the adaptive auxiliary loop to verify the effectiveness of the adaptive loop design.There are three scenarios considered: (1) the basic test of the aerial platform without the manipulator moving to verify the basic performance; (2) the tracking test of the aerial platform under the random moving of the manipulator to verify the effectiveness of the adaptive auxiliary controller; (3) the three-dimensional space tracking test of the manipulator end effector with given ideal reference trajectory to verify the positioning and tracking accuracy of the gripper. All simulation parameters follow the real physical system described in Section2.1.Scenario 1. Let the manipulator be in its initial posture, where the manipulator is vertically below the aerial vehicle (the corresponding initial joint angles are 0°, 90°, and 0°, respectively). At this condition, the gravity center of the manipulator is located on the Zb-axis of the vehicle body-fixed frame, and the manipulator is symmetrical with respect to the Xb- and Yb-axis. Give the step references of the aerial platform position x, y, z and yaw angle ψ when the manipulator remains in its initial posture, and the system responses are shown in Figure 16, where BC refers to the basic controller and AC refers to the adaptive controller with the adaptive auxiliary loop.Figure 16
System responses in Scenario 1.
(a)(b)(c)(d)(e)(f)The results show that when the manipulator is not moving, the basic controller can track the desired reference well with no overshoot and steady-state error. The response times of the longitudinal channel and the lateral channel are about 4 s, the altitude channel is about 3 s, and the yaw channel is the shortest, which is 2.5 s. The results illustrate that the innovative structure design of the aerial manipulator with tandem ducted fans is feasible, and the basic controller can achieve good performance in basic condition. The result with the adaptive auxiliary loop is also given in Figure16, which has faster response time than the basic controller, especially in the longitudinal channel and the lateral channel, which are 2.5 s and 2 s, respectively. This is because the existence of the manipulator will produce gravity moment in the opposite direction when the vehicles roll or pitch to generate the control moment, which makes the actual control moment less than the ideal value, and the adaptive auxiliary controller can compensate for it so as to improve the response speed of the system. Since the initial posture of the manipulator coincides with the Zb-axis of the vehicle, it has little effect on the altitude channel and yaw channel, and thus the response time of the adaptive controller is close to that of the basic controller.Scenario 2. Let the manipulator joints move randomly; the tracking responses of the aerial platform position and yaw angle under the manipulator disturbances are shown in Figure 17. Similarly, BC refers to the basic controller and AC refers to the adaptive controller. The motions of the three joints of the manipulator are shown in Figure 18. Notice that the initial angles of the three joints are 0°, 90°, and 0°, respectively, and the maximum ranges of the motion follow the limits in Table 2.Figure 17
System responses in Scenario 2.
(a)(b)(c)(d)(e)(f)Figure 18
Random motions of the manipulator joints in Scenario 2.The results show that the system responses with the basic controller deviate significantly from the expected references and cannot meet the performance requirements under the random motion disturbances of the manipulator. The lateral channel has the largest fluctuation with a maximum close to 100%, the longitudinal channel takes the second place, and the influences of the altitude channel and the yaw channel are relatively small. This is consistent with the structural characteristics of the aerial manipulator. As shown in Figure3, the manipulator extends from the lateral side of the aerial platform, so it mainly affects the lateral channel of the platform, and the main challenge of the control system is also in the lateral channel. In addition, when joint 1 of the manipulator rotates away from its initial angle, it will also have some influence on the longitudinal channel. The adaptive controller can effectively estimate and compensate for the motion disturbance of the manipulator to ensure good disturbance rejection and tracking performance. Figure 19 also shows the disturbance estimation error of the adaptive controller, indicating that the controller has good estimation accuracy. In order to evaluate the performance of the adaptive controller more accurately, the root mean square error (RMSE) values of the system responses under the two controllers are given in Table 3. The results show that, in the lateral channel, which has the largest fluctuation, the proposed adaptive auxiliary loop reduces the tracking error by more than 40%. The two controllers have similar performance in the altitude and yaw channel due to the correspondingly small disturbance caused by the manipulator.Figure 19
Disturbance estimation performance in Scenario 2.Table 3
RMSE in Scenario 2.
State variablesx (m)y (m)z (m)ψ (deg)Basic controller0.09810.13590.07598.8462Adaptive controller0.07820.08100.07509.0389Scenario 3. Given the reference trajectory of the end effector of the manipulator in three-dimensional space, the responses of the end effector with and without the adaptive loop are shown in Figures 20–22. The reference trajectory starts from the manipulator initial posture, then passes through four reference points in the Y-Z plane, with each point staying for 3 seconds, and finally returns to the initial posture. In order to fully evaluate the performance of the controller, three conditions are considered, of which Scenario 3.1 is the nominal condition, Scenario 3.2 assumes that the end effector grabs a 100 g object, and Scenario 3.3 considers gust disturbances, simulated by the Dryden model [40].Figure 20
System responses in Scenario 3.1.
(a)(b)(c)Figure 21
System responses in Scenario 3.2.
(a)(b)(c)Figure 22
System responses in Scenario 3.3.
(a)(b)(c)In Scenario 3.1, the adaptive controller can achieve accurate trajectory tracking of the end effector with a tracking error of approximately ±0.01 m. The response with the basic controller has a serious deviation in theY direction, and the system cannot completely track the desired trajectory. Since the manipulator moves in the Y-Z plane, the fluctuations of the system in the X and Z directions are relatively small. In Scenario 3.2, the adaptive controller can also track the desired position, but, due to the influence of the grasping object, the overall tracking accuracy of the system decreases compared with Scenario 3.1, and its maximum fluctuation in the Y direction is about 0.02 m. A 0.04 m offset occurs in the altitude direction under the action of the object, and then the adaptive controller quickly realizes the estimation convergence and compensation. Similar to Scenario 3.1, the trajectory with the basic controller also has a serious deviation in the Y direction. Moreover, it diverges rapidly in the Z direction. In Scenario 3.3, the responses of the system in all three directions have greater fluctuations due to the influence of wind gusts. The adaptive controller can still basically achieve the tracking of the reference, and by contrast the responses with the basic controller experience severe fluctuations, which cannot meet the system requirements. In conclusion, the proposed adaptive auxiliary controller can realize precise positioning and trajectory tracking of the end effector of the aerial manipulator.


## 5. Conclusion

This paper proposes an innovative aerial manipulator based on tandem ducted fans, which can easily realize the omnidirectional manipulation of side-on and below and is particularly suitable for the tasks in confined environment. The comprehensive nonlinear dynamic model of the aerial manipulator is established, based on which a control-oriented state space model of the aerial platform is derived using frequency domain identification method. The identification model is verified in both time and frequency domain. A composite controller of the aerial platform is designed under the existence of manipulator. First, a robust H-infinite controller is designed to ensure the basic stability, decoupling, and robust performance of the platform. Then, an auxiliary adaptive controller is designed to estimate and compensate for the effect of the manipulator dynamics. The computer simulation tests verify the feasibility of the proposed innovative aerial manipulator design and the effectiveness of the proposed controller. The results show that, with the proposed composite adaptive controller, the aerial vehicle can achieve good performance under the manipulator motion, and the manipulator end effector can realize precise positioning and trajectory tracking with the grasping load and wind gusts.In future work, the indoor positioning system will be studied and applied, so as to lay the foundation for actual flight tests. In addition, the trajectory planning of the aerial manipulator will also be studied. Then the actual grasping test will be carried out. Furthermore, the vision-based target recognition and the interaction dynamics would be considered to achieve more complex operation tasks.



---

*Source: 7923539-2020-09-08.xml*