Quadrotor Trajectory Tracking Using Model Reference Adaptive Control, Neural Network-Based Parameter Uncertainty Compensator, and Different Plant Parameterizations

: A quadrotor trajectory tracking problem is addressed via the design of a model reference adaptive control (MRAC) system. As for real-world applications, the entire quadrotor dynamics is typically unknown. To take that into account, we consider a plant model, which contains uncertain nonlinear terms resulting from aerodynamic friction, blade ﬂapping, and the fact that the mass and inertia moments of the quadrotor may change from their nominal values. Unlike many known studies, the explicit equations of the parameter uncertainty for the position control loop are derived in two different ways using the differential ﬂatness approach: the control signals are ( i ) used and ( ii ) not used in the parametric uncertainty parameterization. After analysis, the neural network (NN) is chosen for both cases as a compensator of such uncertainty, and the set of NN input signals is justiﬁed for each of them. Unlike many known MRAC systems with NN for quadrotors, in this study, we use the k x x + k r r baseline controller, which follows from the control system derivation, with both time-invariant (parameterization ( i )) and adjustable (parameterization ( ii )) parameters instead of an arbitrarily chosen non-tunable PI/PD/PID-like one. Adaptive laws are derived to adjust the parameters of NN uncertainty compensator for both parameterizations. As a result, the position controller ensures the asymptotic stability of the tracking error for both cases under the assumption of perfect attitude loop tracking, which is ensured in the system previously developed by the authors. The results of the numerical experiments support the theoretical conclusions and provide a comparison of the effectiveness of the derived parameterizations. They also allow us to make conclusions on the necessity of the baseline controller adjustment.


Introduction
This study is an extended version of a previously published conference paper [1]. The problem of quadrotor control quality improvement under the condition of parameter uncertainty has gained considerable attention from the control community in the recent years [2] because of the numerous outdoor and indoor civilian applications of such devices in the following: agriculture, traffic monitoring, delivery, mapping coverage, etc.
The main advantages of the quadrotors are their small size, light weight, vertical take-off and landing, and dynamical maneuverability. However, on the other hand, their dynamics is non-linear with couplings between translational and rotational motions; we have only four control signals (thrust and rotational torques) for six degrees of freedom. Moreover, some parameters of such a plant can not be known precisely: aerodynamical coefficients, its mass, and inertia moments. For instance, if the quadrotor is spraying pesticides on crops, its mass is getting lower like a smooth function, while in the case of cargo delivery, such variation happens in a step-like manner. It should be noted that the inertia moments change their values in both cases.
Despite all the above-mentioned difficulties, the quadrotor position control is required to be of high quality, especially if the indoor missions are considered. To solve such a complicated problem, the control system is usually implemented in accordance with the cascade principle [3] with the inner attitude and outer position loops. The plants of both have their own parameter uncertainties to be compensated. The scope of this research is the position loop only, but interested readers are referred to [4], in which a short review of the existing methods of attitude adaptive control is given, and an effective solution of such problem is proposed.
As far as the trajectory tracking is concerned, the following different techniques have been applied to solve the control problem under consideration: the conventional PIDcontrol [5], LQ-regulator [5,6], backstepping [7], and feedback linearization [8]. However, all of them require the plant dynamics (model parameters) to be known a priori.
To relax such a strict requirement, two solutions should be mentioned. The first is an active disturbance rejection control [9], which is based on the idea of the disturbance (uncertainty) estimation with the help of the extended state observer (ESO). It is a linear high-gain observer, in which the large gain may induce the well-known peaking phenomenon. So, such a scheme is not considered further.
Secondly, the adaptive control techniques can be applied. In particular, the model reference adaptive control (MRAC) with the conventional adaptive law of the form of −Γe re f PBx is used in [10]. Unfortunately, it is not explained why the baseline controller is chosen as the PD instead of k x x + k r r, as well as the absence of the σ or e robust modification in the adaptive law to compensate for the parameter uncertainty, given that its basis functions are approximated using the set of radial-basis ones. The presented stability analysis of such a scheme holds only in case the uncertainty parameters are time invariant.
As the parameter uncertainty is both time varying and nonlinear, it is worth applying the universal approximators such as neural networks (NN) to compensate for it. NN is used to implement compensators in the MRAC-like robust scheme with a PD-based baseline controller [11] and the schemes of backstepping in [12,13]. The solution from [11] requires the known upper bounds of some constituent parts of the parameter uncertainty, while the application of backstepping results in complex controllers of a high dynamic order. Concerning NNs in [11][12][13], the input vector components are chosen without analysis of what signals form the uncertainty to be compensated, and only the output layer parameters are adjusted. So, the higher number of hidden neurons is to be chosen to obtain the required approximation error [14]. More detailed analyses of NN applications to solve the problem under consideration can be found in recent studies [15,16]. As for [15,16], the proposed solutions use the PI/PID controller with time-invariant parameters as a baseline, which contradicts the classical principles of MRAC system synthesis.
All in all, despite a great number of studies on the application of MRAC + NN to control quadrotor trajectory, all of them contains at least some of the following disadvantages: (d1) A simplified model of the quadrotor with time-invariant parameters of mass and inertia moments is used. (d2) The explicit equation of the parameter uncertainty to be compensated is not presented, as well as the numerical simulation results, which compare the outputs of the proposed compensator with the ones of such equation. This makes it difficult to understand to what extent the uncertainty has been compensated. (d3) The baseline controller is chosen arbitrarily but not as a result of the synthesis procedure. (d4) Baseline controller has time-invariant parameters, and plant parameterization, under which it can also be adjustable, is not considered.
This study fills this gap and designs the direct MRAC schemes with the NN-based compensator of matched parameter uncertainty for a quadrotor, which are free from the above-mentioned disadvantages.
The following contributions of this study are that simultaneously: (c1) The explicit equations of the parameter uncertainty are derived for the quadrotor trajectory tracking problem for two cases, when control signals are (i) used and (ii) not used in the parametric uncertainty parameterization; (c2) Using (c1), the application of NN-based uncertainty compensator and signals included into its input vector are justified; (c3) The MRAC-based schemes are implemented with the baseline controller of a k x x + k r r type (in two variants, with time-invariant and adjustable k x and k r ) and the NN-based compensator of the paramtric uncertainty, in which the parameters of both output and hidden layers are adjusted in real time.

Mathematical Model of Quadrotor
Schematic representation of quadrotor is shown in Figure 1. The equations of the quadrotor model [8] are written as follows (the brackets and the time argument are omitted whenever it is clear from the context and causes no confusion): where m is a plant mass, ψ is roll, θ stands for pitch, and φ denotes yaw, ξ = ψ, θ, φ T . The center of gravity of the quadrotor is the origin of a body-fixed frame (BFF). Considering the inertial frame (IF), p = (X, Y, Z) is the coordinate vector, V = ẊẎŻ T stands for the linear velocity, and w = ψθφ T is the plant angular velocity. As far as the BFF is concerned, the same vectors are denoted as p b = (x, y, z), V b = u, v, w T , and ω b = p, q, r T , respectively. F ext = f x , f y , f z T and τ ext = τ ψ , τ θ , τ φ T are the external forces and moments in BFF. τ ψ , τ θ , τ φ are calculated via ψ, θ, and φ controllers. R and J are transformation matrices, which are expressed as follows: Hereafter, s (.) stands for sin(.), c (.) is cos(.), and t (.) denotes tan(.).

Assumption 1.
The inertia matrix I b = diag(J x , J y , J z ) is diagonal as the quadrotor is symmetric in xB − zB and yB − zB planes.
Taking into account the drag force [8], Equation (1) are represented as follows: where K c , C x , C y , C z stand for aerodynamic coefficients; k ψ , k θ are the air resistance parameters; ρ is the air density; S x , S y , S z are the quadrotor surface areas corresponding to the respective axis; l is the distance between the gravity center and each motor rotor; J r stands for the inertia moment of the motor; and ω r is the total motors speed. The output of the altitude controller is T. Lift forces T 1 , T 2 , T 3 , T 4 and rotational torques M 1 , M 2 , M 3 , M 4 are calculated from T, τ ψ , τ θ , τ φ , using (13) from [8].

Trajectory Tracking Control Problem
The parameters values of the quadrotor model are unknown. p, q, r and u, v, w are measurable. Then, all other states from (3)-(5) can be obtained.

Assumption 3.
The inner attitude control system that includes ψ, θ, φ loops has already been designed using the approach from [4]. That means that despite (a) J x , J y , J z being unknown and timevarying and (b) having ∆ pqr in (4) and J in the right-hand part of (3), the NN-based compensator of the the inner loop proposed in [4] ensures ψ, θ, φ control quality defined using the attitude reference model.

Goal.
To solve the quadrotor trajectory tracking problem via development of MRAC system (green blocks in Figure 2), which is capable of the parameter uncertainty (see (5)) compensation and ensures asymptotic convergence of the tracking error to a compact set.

Uncertainty Parameterization
The left equation in (3) is differentiated with regards to time. The substitution of (5) into the obtained result yields the following: The actual quadrotor mass m is unknown and time varying. However, we know its nominal constant value m.
The aim is to ensure trajectory tracking independent of the value of the yaw angle, i.e., X d , Y d , Z d are to be followed with the required quality at any value of φ. To achieve this, the differential flatness approach [17] is applied: This equation is actually an implementation of the converter block in Figure 2. Further parameterization depends on whether u x , u y , u z is directly used in the uncertainty parameterization or not.

Case I: Control Signals Are Directly Used in Uncertainty Parameterization
Equation (6) is rewritten as follows: The only constant in (8) is m. Equation (7) is substituted into (8) to obtain the following: is the matched parameter uncertainty to be compensated. It can be rewritten as a linear regression with unknown m (time varying), ρ, C x , C y , C z , S x , S y , S z and measurable nonlinear regressor, which depends on the twelve known signals u x , u y , u z , ψ, θ, φ, u, v, w, p, q, r. The obtained parameterization has two main features: (i) the input gain matrix is known and identified, and (ii) control signals u x , u y , u z are included into uncertainty parameterization.

Case II: Control Signals Are Not Used in Uncertainty Parameterization
Equation (6) is used, and (7) is substituted into it as follows: In such a case, the matched parametric uncertainty again can be rewritten as a linear regression, but its regressor does not include u x , u y , u z . Furthermore, the system input gain matrix is unknown and represented as a diagonal matrix with m m on its main diagonal.

Representation of Plant in State-Space Form
The state x p = X,Ẋ, Y,Ẏ, Z,Ż T and control vectors u = u x , u y , u z T are introduced. Then, (9) and (10) are written as follows: Considering (9), we have the following: while it is written for (10) as follows: Λ is unknown. So, in the course of the following design, we will need u base with time-invariant parameters for Case I and adjustable ones for Case II.
The uncertainty ∆ is time varying and has a nonlinear regressor for both cases, so NN is chosen to compensate for it.

Reference Model
The reference model is introduced as follows: The following definitions are used:

Plant Representation: Neural Network Description
According to the proved approximation theorem [18], ∆ is expressed as NN with the sigmoid activation function σ in the hidden layer ( Figure 3) as follows: On the basis of the obtained parameter uncertainty parameterization, it is concluded that the number of the input neurons is N 1 = 12 for Case I (x nn = {u x , u y , u z , ψ, θ, φ, u, v, w, p, q, r}) and N 1 = 9 (x nn = {ψ, θ, φ, u, v, w, p, q, r}) for Case II, the number of output neurons is N 3 = 3 for both cases. ε is the approximation error, of which its upper bound ( ε ≤ ε N ) can be estimated using the results from [14]. It depends on the number of hidden neurons N 2 , the higher N 2 , and the lower ε N .
As ideal weights W and V are unknown, NN with adjustable parametersV andŴ is introduced as follows [19]: The adaptive laws forV andŴ will be derived further for both cases under consideration.

MRAC System Design
The approximation and estimation errors are introduced as follows: The dependence of e from V is nonlinear due to the sigmoid function, so it is linearized via the Taylor series expansion as follows: Introducing e re f = x p − x re f , the error equation is obtained as follows: The summands of u(t) are chosen as follows: As far as Case I is concerned, k x and k r can be directly calculated, while in Case II, this is impossible as B p is unknown. In such a case, the baseline control law with adjustable parameters should be introduced as follows: u base =k x x p +k r r =θ T ω, As for Case I, the substitution of (15) and (17) into (16) yields the following: Considering Case II, we have the following: We are in a position to present the main result of this research. If the parametric uncertainty parameterization includes the control signals u x , u y , u z , the following theorem holds. Theorem 1. Let the laws to adjust NN parameters be defined as follows:V where Γ v and Γ w are adaptive gains, σ v and σ w are σ-modification parameters, P = P T > 0 stands for the solution of the Lyapunov equation A re f T P + PA re f = −Q, the augmented error ζ = e T re f vec T W vec T Ṽ T is uniformly ultimately bounded (ζ ∈ UUB) with the ultimate bound R, and e re f asymptotically converges to a compact set with the bound R 1 .
Proof. The proof of this theorem, as well as the definitions of R and R 1 are presented in Appendix A.
If the uncertainty parameterization does not include control signals u x , u y , u z , the following theorem holds. It is noted that in such a case, we can write B p = m m I 3×3 ⊗ 0 1 .

Theorem 2.
Let the laws to adjust the baseline controller and NN parameters be defined as follows:θ where Γ θ , Γ v , and Γ w are adaptive gains, σ θ , σ v , and σ w are σ-modification parameters, P = P T > 0 stands for the solution of the Lyapunov equation A re f T P + PA re f = −Q, the augmented error is uniformly ultimately bounded (ζ ∈ UUB) with the ultimate bound R 2 , and e re f asymptotically converges to a compact set with the bound R 3 .
Proof. The proof of this theorem, as well as the definitions of R 2 and R 3 , are verbatim to the one of Theorem 1 up to one additional term in the Lyapunov function for the baseline controller adjustable parameters.

Numerical Experiments and Discussion
On the basis of Equations (1)-(5), the mathematical model of the quadrotor has been implemented in Matlab/Simulink. The values of its parameters were taken from the Parrot Mambo model [4]. The Euler solver with fixed step of 0.001 s was used for modeling.
The same flight plan was chosen for all experiments as follows: (1) to take-off to reach 1.1 m height (from start to fifth second), (2) to track Gerono lemniscate (figure eight) trajectory from the 5th to 65th seconds: X d = cos(2πt) and Y d = sin(4πt).
The parameters of MRAC system were picked as follows: ω n = 5, γ = 1.1, N 2 = 50, The known nominal values of inertia matrix and mass were chosen as follows: m = 0.063 2 kg, I b = diag ( 5.8 2 ) × 10 −5 , ( 7.17 2 ) × 10 −5 , ( 1 2 ) × 10 −4 kg·m 2 . The quadrotor model parameters were unknown to the control system. In the course of experiments, I b and m were changed as follows: (1) from start to 35th second (2) from 35th to 65th second I b = diag (5.8 × 10 −5 , 7.17 × 10 −5 , 10 −4 kg·m 2 , m = 0.063 kg. The attitude control was implemented in accordance with [4], while the system proposed in this study was used for the position control. The yaw setpoint was picked as φ d = π 4 . Here, we outline some comments on how the parameters for simulation were chosen. All quadrotor (and environment) parameters were taken from [4]. The initial values of the controllers parameters were picked as zeros. The reference model parameters were calculated on the basis of the analysis of the results from [4]. The values of k x and k r were calculated using (17). The value of φ d could be chosen arbitrarily. The values of the adaptive gains and σ-modification parameters were picked as a result of trial and error. The trade-off of such a choice was as follows. The only way to reduce the tracking error was to increase the adaptive gain. However, this increased the L 2 norm of the control signal, which resulted in oscillations [20]. Therefore, the trade-off was between the values of the tracking error and the oscillations. The values of the σ-modification parameters were in relation with the corresponding adaptive gain. Their lower values resulted in higher controller parameters and lower tracking error, and vice versa (please see proof of Theorem 1 in Appendix A).
The experiments were conducted as follows: (i) with NN-based compensator designed for Case I, (ii) with NN-based compensator designed for Case II, and (iii) without compensator.
First of all, the compensator for Case I was compared with the system without uncertainty compensation. Figure 4 presents the behavior of the Z coordinate. When the control system included only the baseline controller, a steady-state error occurred due to the difference between m and m. Such an error was compensated for in the neural network. The behavior of the coordinates X and Y is shown in Figure 5. The sum of squared errors (the difference between the reference figure-eight trajectory and the real one) was reduced via the NN-compensator by 16.62%. The obtained equation for ∆ was used to calculate the ideal value of the parameter uncertainty for each time instant of the numerical experiment. Figure 6 illustrates that the NN-based compensator approximated such ideal ∆ and verifies the theoretical result. The fact that m value was changed at 35th second did not cause the peaking phenomenon. As the attitude control loop synthesized according to [4] also included its own NNbased compensator, Figures 7 and 8 are shown to demonstrate that the parameter uncertainty of the ψ, θ, φ control loops and J x , J y , J z switch were approximated as well.  The same results were obtained using different values of φ d . Therefore, the validity of the converter synthesis procedure for Case I was corroborated.
The next step was to test the adaptive control system designed for Case II. As mentioned above, the initial values of the baseline controller adjustable parameters were zeros. The results obtained for X, Y, Z coordinates are shown in Figures 9 and 10. Figure 9. Transients of Z coordinate for control system with NN-based compensator for Case II using baseline controller with zero initial conditions. They demonstrate that the adaptive control system stabilized the quadrotor and ensured the asymptotic convergence of the tracking error to a compact set, but the steadystate error is too large for X and Y. This was mainly caused by the fact that both baseline and adaptive controllers started from zero initial conditions, which did not provide stabilization. Figures with curves to assess the quality of uncertainty approximation are not provided, as in Case II, both the NN and baseline controller have adjustable parameters and try to compensate for the whole uncertainty together.  Having compared the curves in Figures 4 and 11 and Figures 5 and 12, it was concluded that the design in Case I provided better transients quality for X, Y, and Z coordinates, so the additional adjustment of the baseline controller did not allow us to obtain any advantages. Therefore, according to the conducted experiments, Case I parameterization is a better solution to the stated problem. Additionally, the system on basis of the PID-controllers of altitude, attitude and position from the "ParrotMinidroneHover" project was also applied to the implemented quadrotor model. The comparison of such system with the proposed solution is considered to be fair, as we did not adjust PID-controllers by ourselves but used the best values of their parameters obtained by the developers of the "ParrotMinidroneHover" project, which were successfully applied not only for simulations but also for a real minidrone. The obtained behavior of the X, Y, and Z coordinates are shown in Figures 13 and 14. It was concluded that compared to the MRAC system with NN, the PID-controllers were not capable of the parameter uncertainty compensation.   Furthermore, finally, we compared the proposed solution with the ones in which the PI/PD/PID controllers with time-invariant parameters are used as the baseline. To do so, the scheme with PID-controllers was augmented with the NN-based compensator for the X, Y, and Z coordinates. The obtained results are presented in Figures 15 and 16. Comparing Figures 5, 14 and 16, it was concluded that the simple combination of some PID-controllers and NN-based compensators did not allow us to obtain results similar to the proposed solution. So, as in [15,16], the PID-controllers has to be adjusted via a trial and error procedure. Therefore, in [15,16], such a procedure is much more difficult, as both baseline controller parameters and adaptive controller hyperparameters have to be chosen manually. As far as the proposed solution (Case I) is concerned, the baseline controller parameters are calculated using mathematically-sound equations, and only adaptive gains and σ-modification parameters for NN are chosen via trial and error.  All in all, the obtained experimental results can be summarized as follows: (r1) Despite the fact that systems designed on the basis of Case I and Case II parameterizations had the same theoretical properties, the MRAC system with the NN-compesator with time-invariant baseline controller parameters allowed us to obtain better results in comparison with the adjustable baseline controller.
(r2) The system on the basis of PID-controllers was not able to fully compensate for the parametric uncertainty. (r3) The simple combination of the PID-based control system with the NN-based compensator did not allow us to obtain the same results as the proposed approach. So, the baseline controller should be derived on the basis of the MRAC design procedure and have the form k x x + k r r. Moreover, as far as Case I parameterization is considered, the values of k x and k r can be directly computed.

Conclusions
Considering the trajectory tracking problem, the equations of the quadrotor parameter uncertainty were obtained for two cases: control signals are (i) used and (ii) not used in the parametric uncertainty parameterization. Their analysis allowed us to choose the following: (1) the neural networks to implement the compensators of Case I and Case II uncertainties and (2) the signals to form their input layers. For Case I, the adaptive laws of NN parameters were derived. For Case II, the adaptive laws of both NN and baseline controller parameters were obtained, which ensured that the augmented tracking error was UUB. The conducted experiments demonstrated that, as far as transient quality is concerned, Case I parameterization seems to be a better solution.
The main problem of the proposed solution is the above-mentioned trade-off between the value of the augmented error and the transient quality, particularly, oscillations. The scope of further research is to (i) derive a parameterization of the quadrotor to apply an identification-based adaptive control and some extension and mixing schemes to ensure some transient quality, (ii) guarantee exponential convergence instead of the assymptotic one, and (iii) consider cases when the quadrotor dynamics is affected by the wind and gravity center displacement due to load.  Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Proof. Let a Lyapunov function candidate be chosen as follows: L e re f ,Ṽ,W = e T re f Pe re f + tr Ṽ T Γ −1 vṼ + tr W T Γ −1 wW , A T re f P + PA re f = −Q, Its derivative is written as follows: L = −e T re f Qe re f + 2tr W T − σ V T x nn −σ V T x nn V T x nn × e T re f PB p + Γ −1 wẆ + 2tr Ṽ T −x nn e T re f PB pŴ Tσ V T x nn + Γ −1 vV +2e T re f PB p δ.
(A2) Substituting (21) into (A2), it is obtained as follows: L ≤ −λ min Q e re f 2 + 2λ max (P) B p e re f δ 0 − 2 B p σ v Ṽ 2 The first two terms of the right hand side of (A3) are added to full square. Then, after routine operations, the upper bound ofL is rewritten as follows: Applying the comparison lemma, the inequality is solved as follows: So, ζ(t) is uniformly ultimately bounded. To prove the asymptotic convergence of the tracking error e re f , we use (A5) and take into account the definitions of γ 2 and t → ∞. As a result, the upper bound of the tracking error is obtained as follows: e re f ≤ 2λ max (P) λ min (Q)λ min (P) This completes the proof.