MIMO PID Controller Tuning Method for Quadrotor Based on LQR/LQG Theory

Abstract: In this work, a new pre-tuning multivariable PID (Proportional Integral Derivative) controllers method for quadrotors is put forward. A procedure based on LQR/LQG (Linear Quadratic Regulator/Gaussian) theory is proposed for attitude and altitude control, which suposes a considerable simplification of the design problem due to only one pretuning parameter being used. With the aim to analyze the performance and robustness of the proposed method, a non-linear mathematical model of the DJI-F450 quadrotor is employed, where rotors dynamics, together with sensors drift/bias properties and noise characteristics of low-cost commercial sensors typically used in this type of applications are considered. In order to estimate the state vector and compensate bias/drift effects in the measures, a combination of filtering and data fusion algorithms (Kalman filter and Madgwick algorithm for attitude estimation) are proposed and implemented. Performance and robustness analysis of the control system is carried out by employing numerical simulations, which take into account the presence of uncertainty in the plant model and external disturbances. The obtained results show the proposed controller design method for multivariable PID controller is robust with respect to: (a) parametric uncertainty in the plant model, (b) disturbances acting at the plant input, (c) sensors measurement and estimation errors.


Introduction
Small unmanned aerial vehicles (SUAV) are used in a great variety of application fields [1][2][3]. In order to achieve the required level of performance, the vehicle must include a proper automatic control system as well as electronics, sensors, communications, and signal processing algorithms. To develop and validate the control system is very useful to count with an accurate mathematical model of the vehicle, sensors and actuators dynamics [4][5][6][7][8][9][10][11][12][13]. The present work uses the model structure exposed in [6,12,14] as a starting point but, as a novel contribution, introduces a nonlinear mathematic model (obtained experimentally) of the rotors used in the DJI F-450 quadrotor (Figure 1), that is used to perform a more detailed and realistic analysis of the control system through numerical simulations.
One crucial concern to allow the indoor operation of a SUAV is the attitude and position estimation, typically based on the use of inertial measurement units (IMU) type MEMS (MicroElectroMechanical Systems) and cameras [6,8,[15][16][17][18][19][20][21]. To estimate the attitude, position and velocity state variables of the vehicle from the measurements provided by the sensors, a variety of methods are used, such as Additionally, the uncertainty affecting the measurements are taken into account through the errors and noise modelling of a LIDAR, an optical flow camera, and a MEMS type MARG (Magnetic, Angular Rate, Gravity) sensor [42][43][44][45]. The attitude of the vehicle is estimated using the sensor fusion algorithm proposed in [46], whereas the altitude and the vertical velocity are estimated by means of a Kalman filter. Thus, this study considers the effects of measurement errors, including the relatively large bias and noise errors inherent to the MEMS sensors technology, and includes all the required estimation and compensation algorithms necessary for implementation in a real system. A horizontal controller is also included to analyze the influence of an external control loop over the proposed modified LQR/LQG. Next, we study the performance and the robustness of the proposed controller design method through numeric simulations, considering a non-linear model of the DJI F-450 quadrotor and the rotors non-linear dynamics obtained experimentally. The influence of uncertainty and bounded perturbations on the robustness of the control system is also tested adding parametric uncertainty and significant external disturbances.
The organization of this paper is as follows: In Section 2, the mathematical model obtained experimentally, which is used to the analysis and evaluation of the control system, is presented, followed by a simplified model used for controller design. In Section 3, the proposed design methodology to obtain the attitude and altitude control laws, as well as the employed state estimation algorithms, are described. In Section 4 address the problem of the estimation and control of the horizontal position. Section 5 provides numeric simulation results. Finally, conclusions are drawn in Section 6.

Mathematical Model
To carry out the synthesis and evaluation of the control system two different models of the DJI F-450 quadrotor will be developed. The first is a non-linear model to analyze and evaluate (AEM model) the designed controllers. This model incorporates the effects of the rotors dynamics and measurement errors. The second is a non-linear model for control system design (CDM model). This model includes some simplifications compared with the previous one that will lead to a linear model. The linear model is used to design the controllers and in several estimation algorithms.

Analysis and Evaluation Model (AEM)
As it was said previously, the present work uses the model structure exposed in [6,12,14] as a starting point. Figure 2 shows the Earth and Body reference frames. The translation dynamics equations are: mu = (cos φ sin θ cos ψ + sin φ sin ψ)T b − K dx u|u| mv = (cos φ sin θ sin ψ − sin φ cos ψ)T b − K dy v|v| mẇ = (cos φ cos θ)T b − K dz w|w| − mg T b = T 1 + T 2 + T 3 + T 4 (1) where (x E , y E , z E ) is the position vector of the vehicle relative to the Earth frame, (u, v, w) is the linear velocity (u =ẋ E , v =ẏ E , w =ż E ), and T b is the total thrust force generated by the four rotors relative to the Body reference frame. On the other hand, the rotational dynamics is described by: where (φ, θ, ψ) represents the attitude using Euler angles, (p, q, r) is the angular velocity vector (measured in the Body frame), T i is the thrust force generated by the i rotor (relative to the Body frame), and N i is the rotation frequency of the i rotor (in RPM), for i ∈ {1, 2, 3, 4}.   (1) and (2), that have been obtained experimentally [47] through laboratory tests with the DJI F-450 quadrotor used in this work. The inputs of the subsystem composed by the electronic speed controllers (ESC), the brushless motors, and the propellers are the duty-cycles values of the PWM (Pulse Width Modulation) signal applied to each ESC, D i , and its outputs are the thrust forces and torques developed by each rotor. Despite the model structure being based on that proposed by [6,12,14], its parameters have been identified to obtain a trustworthy model of these parts of the DJI F-450 quadrotor. The identification has been carried out applying step and ramp signals to the inputs (D i ) and recording the rotation frequencies, thrusts, and torques. The following equipment has been used to accomplish the tests: a calibrated SEN-CC-FZ1439 load cell, a QDR1114 optical sensor functioning as a tachometer, an NI USB-6211 acquisition system, and an LPC4088-32 microcontroller development board. Figure 3 shows one of the experimental set-ups. The reference [47] provides further details. Each group ESC-motor-propeller is modeled as a first-order system with delay, whose dynamics are described by where τ is the system time constant, τ d is the delay, and K is the stationary gain. The values derived experimentally in the hover equilibrium condition were: τ = 0.038 s, τ d = 0.012 s y K = 82.3 RPM/%. The experimental relations that link the rotational frequency N (in RPM) with the thrust, T, and the torque, P, developed by each rotor are given by the Equations (4) and (5), respectively.
In order to carry out a realistic analysis and design, the errors and noise levels of conventional commercial sensors employed in these type of systems have been obtained from the manufacturer manuals and introduced in the model. These sensors are a PulsedLight LIDAR-LITE laser range sensor [42], an InvenSense MPU-9250 MARG [43,44], comprising a 3-Axis gyroscope, a 3-Axis accelerometer, and a 3-Axis magnetometer, and a PX4-FLOW optical flow camera [45].

Control System Design Model (CDM)
To ease the task of designing the controller, we will build a simplified model relative to the one presented in the previous section. The main differences between the two are the following:

•
The CDM model does not take into account the rotor dynamics.

•
The thrust and the torque are linear functions of the PWM signals duty-cycle.

•
The rotational dynamical model assumes relative low angular velocities and small enough Euler angles.
After applying the cited simplifications, the CDM model is defined by the system of non-linear Equations (6). Table 1 shows the parameter values of this model.
The relations between the components of the U and D vectors (the duty-cycles of the four PWM signals) is defined by [6]: where D is the vector of duty-cycles, and the components of U are the fictitious control signals U φ , U θ , U ψ y U z that allows split control of the attitude (φ, θ, ψ) and the altitude (z) variables. Thus, we can choose between considering the whole 4 × 4 MIMO problem or one 3 × 3 MIMO attitude problem plus one SISO altitude problem instead. The next section elaborates on the first option.

Attitude and Altitude Control System
The main objective of this system is the stabilization (hovering) of the vehicle at a certain altitude. Hovering is a manoeuvre in which the quadrotor is maintained in nearly motionless flight over a reference point at a constant altitude and on a constant heading. In order to control the vehicle, the complexity of the system and the strong coupling that exists between the manipulated variables and the controlled variables have been taken into account. Therefore, a multivariable control approach based on state space has been considered [40,48,49]. As a novel contribution of this work for multivariable PID tuning, it is proposed to use a modified LQR/LQG control law, in which feedback is performed considering the tracking error between the state vector and a setpoint or reference state vector. ISpecifically, for the design of the attitude and altitude controller the following state vector, X, is considered: The manipulated variables (MV) vector, U, is composed of the four components U φ , U θ , U ψ , and U z given in (6). The process variables (PV) vector, or controlled variables (CV) vector, is given by, An equilibrium condition for hovering at an altitude z eq is characterized by, Linearizing the eight non-linear differential equations of the quadrotor model (6) for hover condition, (X eq , U eq , Y eq ), and considering incremental (or deviation) variables with respect to equilibrium point, the linearized system dynamics is defined by (12).
In order to add integral action to the controller, the system model is modified including four new state variables corresponding to the integral of the tracking error of the controlled variables, resulting in the following augmented state vector, whereẋ I = r − y, and r is the incremental value of the set-point vector, r = Y SP − Y SPeq . Therefore, the augmented system matrices are given by, The conventional LQR control law with integral action (incremental vector value) is given by, and the controller output (CO) vector (U) is obtained by, For that, the feedback gain matrices, K c and K I , must be calculated solving the corresponding algebraic Riccati equation (AREc) [40,48,49], which results from the optimization problem of minimizing a cost function given by, where Q c and R c are the named weighting matrices.
With the LQR controller, the linearized closed-loop control system dynamics (which is an approximation of the nonlinear system for a neighborhood of the equilibrium point defined by the hover condition) is given by, A proposal made in this work is to modify the conventional LQR structure given in (15), where the new control law is defined by, where x SP represents the incremental reference (or set-point) state vector to follow, The proposed control law (19) gets the same poles for the closed-loop linearized system obtained with (15), but achieves to reduce the time response of the system to set-point changes due to (19) modifies the zeros of the closed-loop system reached with (15).

Pre-Tuning Method for Attitude and Altitude Controller
The aim of a pre-tuning procedure is to determine Q c and R c weighting matrices in (17) for which a suitable solution is obtained. For that, diagonal weighting matrices have been chosen, Taking into account experimental (computer simulations) observations of the closed-loop system, and physical meaning of the relative values used in the elements of the weighting matrices, the following relations are proposed as a starting point: (26) in this way, for a given q φ parameter value, the rest of the weighting parameters are calculated for the pre-tuning controller. In practice, due to non-linear dynamics and parameter uncertainty, it is necessary to make a fine-tuning to achieve adequate performance and robustness properties, which is made experimentally in an iterative way. In order to analyze the properties of the proposed controller, an illustrative example is shown below, where the following tuning parameters have been adjusted experimentally, With these weighting matrices, the following state feedback matrices, K c and K i , are obtained, 0.1846 0.0000 −0.0000 0.0000 0.0000 0.1846 0.0000 0.0000 0.0000 0.0000 0.0079 −0.0000 0.0000 −0.0000 −0.0000 1.2247 In order to analyze the effect of the modified LQR control law, in Figure 4 we show a comparison between the time response of the linearized system employing the conventional control law (15) and the proposed control law (19) in this paper. In both cases the initial state vector x T i (0) = (0.1, −0.3, π, −1, 0, 0, 0, 0), is expressed in the International System of Units. As can be seen, the control law (19) Figure 4. Comparative of closed-loop systems obtained with control laws (15) and (19).
The proposed controller can be interpreted as a multivariable PID controller for a MIMO (4 × 4) system, taking into account that the K c matrix can be expressed as 2.9394 0.0000 −0.0000 0.0000 0.0000 2.9408 0.0000 0.0000 0.0000 0.0000 0.5024 0.0000 0.0000 0.0000 0.0000 2.7333 0.3174 0.0000 0.0000 0.0000 0.0000 0.3401 0.0000 0.0000 0.0000 0.0000 0.1542 0.0000 0.0000 0.0000 0.0000 2.8867 where the respective K P , K D (29) and K I (30) matrices are practically diagonal matrices, and therefore they can be approximated by With these approximations, the PID control laws (using Laplace transform domain) for the components of the controller out vector are given by, where e φ , e θ , e ψ , and e z are the respective tracking error signals for the controlled variables, φ, θ, ψ and z E . Figure 5 shows block diagrams with the conventional LQR and the modified LQR control law proposed in ths work. Once the controller output vector U has been calculated, the duty cycle vector, D for the PWM signals applied to motors are obtained taking into account Equation (7), In order to implement the proposed algorithm, all state variables are needed for feedback. For that, and taking into account the characteristics of the system to control, low-cost commercial sensors characteristics are used in computer simulations based on data sheets given by manufacturers. A MARG (Magnetic, Angular Rate, and Gravity) sensors system (MPU-9250 model, with 3-axes accelerometer, rate-gyro, and magnetometer) is used for attitude and orientation measurements ( [43,44]). A LIDAR (Light Detection and Ranging) is used for altitude measurement [42]. Due to bias/drift effects acting on measurements given by inertial sensors, compensation algorithms must be considered as described in the next section.

Attitude Estimation
In order to achieve proper control of the system, it is very important to have accurate enough attitude measurements. The references [50,51] compare the performance of three sensor fusion algorithms that use the MARG measurements for attitude estimation. These are an extended Kalman filter (EFK), the non-linear complementary filter (NLCF) described in [25], and the error optimization algorithm detailed in [46] based on the gradient descent method (MHV). Table 2 shows the relative computation times on a typical embedded system as well as the RMS errors for the three algorithms. This work uses the MHV algorithm since it has a good compromise between computation cost and estimation errors.
In each sample interval, the MHV algorithm proceeds first finding the best alignment between the gravity and Earth magnetic vectors and those measured by the accelerometer and the magnetometer, respectively. To get that optimal alignment of these two pairs of vectors despite the presence of measurement errors the MHV algorithm uses the gradient descent method. The output of this stage is a quaternion ( B Eq w,t ) presenting the rate of change of the attitude. This quaternion is then combined with the obtained from the angular velocity measured by the gyroscope ( B Eq w,t ) to calculate a better estimation ( Figure 6) where β is an adjustable weighting factor representing the gyroscope measurement error expressed as the magnitude of a quaternion derivative. The algorithm also copes with the effects of magnetic disturbances and the bias present in the gyroscope. The angular velocity vector (p, q, r) is supposed to be taken directly from the gyroscope output and, for simulation purposes, is degraded with the noise levels stated in the datasheet [43].

Estimation of Linear Velocity and Position in Axis z E
The sensor LIDAR-LITE [42] gives a measurement from the quadrotor to ground, which can be used for estimating the vertical velocity and position. For that, taking into account the model given in (6), only the equations corresponding toż E andẇ are considered, and it is supposed that φ and θ take values next to zero (which is a good approximation for hover manoeuvre); with what turns out Equation (36).
If (36) is linearized for hover condition (w eq = 0, z eq = 0, U z eq ), a linear Kalman filter (KBF) can be used for w, z E estimation. Covariance matrices, (Q o , R o ) must be adjusted to obtain a satisfactory KBF. For that, an experimental iterative procedure has been applied, obtaining the following matrices, where R o takes into account the noise measurement gives by [42], The stationary observer gain matrix, K o , can be obtained solving the corresponding observer algebraic Riccati equation (AREo) [40,48,49]. In this way, results in an LQG controller, and the complete attitude-altitude controller is a combination of LQR (for attitude control) and LQG controller (for altitude control), called, in this work, LQR/LQG. A block diagram with the complete structure of the control system is shown in Figure 7.

Horizontal Position Controller
Although the purpose of this article is to design a control system for the attitude and altitude (AAC) of a quadrotor, a horizontal position control (HPC) loop is included below in order to test the complete system when the attitude set-point is generated by an additional control loop whose objective is to control the position with respect to the horizontal plane. For that two scalars (SISO) PID controllers are used. The tuning of these PIDs are out of the objectives of this paper, and they are only used for testing the AAC integrated with an external control loop. As it is shown in 7, the PID controllers for HPC generate as control signals the set-points for φ SP and θ SP for the AAC. The output magnitudes of the set-points generate by the PIDs of the HPC are limited to avoid unnecessary pitch and roll angles, as it is shown in Figure 7. The blocks indicated as Filters correspond to first order low-pass filters used for smoothing the set-point changes sent to the controllers. Equation (38) shows the control laws generated by the PIDs with filter to derivative action, employing the Laplace transform notation. Controllers parameters have been obtained experimentally, and are given in Table 3 (the negative sign is due to the sign criteria used for φ and y B , see Figure 2).
where e x B and e y B are the error components of the horizontal position vector (e H ) which are referenced to the body-axis reference system.  The horizontal positioning error can be referenced with respect to the body-axes reference system, e H (e x B , e y B ), as well as to the Earth-axes reference system, e H (e x E , e y E ). For small values of φ and θ angles, the relation between the horizontal coordinates is given by,

Horizontal Position Estimation
The horizontal position estimation of the quadrotor derives from a CMOS optical flow camera model PX4FLOW [45] pointing down. This camera computes a displacement vector utilizing the correlation between successive images, so its output is expressed in pixels units and not in real units. To compute the truth velocity is necessary to take into account the focal length of the camera lens as well as the distance from the camera to the floor, or altitude, that is already available from the LIDAR, as is shown in Equation (40). Also, the horizontal velocity has to be rotated (41) taking into account the ψ angle estimated previously to refer it to the Earth frame. Then, the horizontal velocity is integrated to get the final horizontal position in that frame ( where z is the distance to the floor and f is the focal distance of the camera lens. The manufacturer states that the average error in the camera output is zero, so the position measurement drift should be negligible. v

Results
This section presents the results obtained, in a realistic simulation scenario, when analyzing the LQR/LQG controller modified for the attitude and the altitude proposed in this paper. For this purpose, an analysis of the performance of this control system is carried out by applying it to the AEM model of the quadrotor. The non-linear AEM model, (1) and (2), also takes into account: rotor dynamics (3), measurement errors [42][43][44][45], the effect of estimation algorithms (MHV and KBF) and the effect of adding an external control loop for the horizontal position (38) of the quadrotor. This framework is illustrated in Figure 7, which shows the interaction between measurements provided by sensors, estimation algorithms, and controllers. It should be mentioned that in all the simulations presented the numerical integration method of Runge-Kutta of fourth order with an integration step of 0.1 ms has been used. The estimation (MHV and Kalman filter) and control algorithms are executed at a frequency of 400 Hz in simple arithmetic precision. The discretization of the controllers and the Kalman filter is carried out using the Euler backward method, and β = 0.01 is used in the MHV algorithm for fusion of sensor data. To conclude, a robustness analysis of the attitude and altitude control system to parametric uncertainty and to different types of external disturbances has been realized.
Before the attitude and altitude control performances are analysed, the MARG MPU-9250 sensor modelling is shown in Figure 8. As can be seen, all measurements are affected by noise and bias levels according to the values indicated by the manufacturer [43]. In addition, in the case of accelerometers, the measurement provided (ax m ) is affected by both gravity (ax g ) and the linear accelerations experienced by the quadrotor (ax l ). This consideration has a negative effect on the estimation of the attitude made by the MHV algorithm, which can be seen in Figure 9. As a comparison, this figure also shows the result of integrating the angular velocities provided by the MARG directly to estimate the attitude. This integration diverges with time due to the bias effect present in the measured angular velocities (p m , q m , r m ). On the other hand, the Kalman Filter for altitude and linear velocity w gets to estimate these state variables accurately, see Figure 10. This is in part due to the quality of the measurements provided by the sensor Lidar-Lite [42]. Once the estimated/measured state vector is presented, it can be pointed out that the most negative effect is present in the attitude estimation; the attitude controller must be robust to the coupling presented between the estimate of the angles (φ, θ) and the linear accelerations experienced by the quadrotor.   Using the state vector used for the attitude and altitude controller feedback presented here, the results obtained when applying this controller to the AEM mathematical model are shown. In Figure 11 it can be seen how the modified LQR/LQG control shows that the system follows the setpoints requested for (φ, θ, ψ, z) adequately, despite measurement errors (Figure 8), estimation errors (Figures 9 and 10) and uncertainty in the AEM model with respect to the linear design model. It should be noted that for both attitude and altitude, the error between the estimated variable and the requested setpoint is quite small (±0.2 degrees for φ est and θ est ). Therefore, if sensors and/or estimation algorithms of superior performance are available, the following of setpoints carried out by this controller would have a considerable improvement margin. The duty cycles (D) generated by the modified LQR/LQG controller to get the above tracking are shown in Figure 12. On the other hand, most high-level applications of SUAV systems require the vehicle to be able to track position commands in space or maintain static positions. Therefore, it is necessary to validate the attitude and altitude controller proposed in this paper by adding an external control loop (38) for the horizontal position of the system. As shown in Figure 13, the closed-loop system resultant gets a satisfactory tracking of the setpoints demanded by the positions (x E , y E , z E ). Even with the estimation errors committed by the MHV algorithm (due to the linear accelerations experienced by the system), the modified LQR/LQG control successfully follows the setpoints rapidly (φ sp and θ sp ) demanded by the horizontal position controller. Furthermore, the proposed control structure is valid for tracking complex trajectories. Figure 14 shows the position in the quadrotor space when tracking ascending circular paths. It can be seen that despite the coupling of the system and the errors present in the estimates and measurements, the results obtained can be considered satisfactory in terms of stabilization and setpoint tracking, although they can be improved by fine-tuning the controllers.

Robustness Analysis
In order to study the robustness of the modified LQR/LQG structure proposed in this work, for the pre-tuning of PID controllers applied to quadrotor systems, an analysis of attitude and altitude control is performed by adding: parametric uncertainty to the system inertia (which will determine the system response speed) and parametric uncertainty to the gain of the SISO model identified for the rotors. On the other hand, the robustness of the modified LQR/LQG control by exposing the quadrotor system to significant external disturbances (torques and forces) is also studied. Furthermore, in both cases, all measurement/estimation errors presented in the previous section are considered.

Parametric Uncertainty Robustness Analysis
In this analysis, the nominal inertia moments (I nom ) of the AEM model (Analysis and Evaluation Model, see Table 1) have been multiplied by a factor, I = I nom (1 + ∆I), which takes into account the level of uncertainty. The study focuses on this parameter because I x and I y are the ones that determine the fastest dynamics of the quadrotor and, therefore, those that affect in greater measure to the system when determining the performance of the attitude control. The study of the added uncertainty in the mass of the system has not been carried out, since this can be known with sufficient precision and does not vary with time. As result of the study, Table 4 shows the values obtained when calculating the IAE (Integrated Absolute Error) of all controlled variables (x, y, z, φ, θ, ψ) in function of the uncertainty added to the inertia moments (I add = ∆I × 100%). In addition, it records the IAu (Integrated Absolute) of the incremental variable u, which is defined as the sum of the incremental duty cycles (d) with respect to the equilibrium condition (D eq ). All the simulations carried out to obtain the values listed in Table 4 are reproductions of the one shown in Figure 13 and, moreover, the time interval used for its computation is equal to the sampling time (Tm = 2.5 ms). As visual support to the above, Figure 15 shows the variation of the IAE attitude-altitude indicators and the IAu control effort indicator.
As can be seen, the low variation of the indicators for uncertainty levels close to the nominal value ([−17.5, 17.5]%) shows that the modified LQR/LQG controller, together with the horizontal position controller, presents robust behaviour to the uncertainty added to this parameter. Furthermore, the control system achieves to stabilize the system with acceptable stationary error values (±0.5 o for attitude and ±0.04 m for position) in an uncertainty range of [−19.5, 80]%, see Figure 16. Secondly, it should be noted that the magnitude of the IAE indicators decreases considerably if these are calculated with respect to the estimated variables, see Table 5. This supports the previously mentioned, the results obtained when applying the control structure proposed in this work to the AEM model have a considerable margin of improvement if the quality of the measures/estimates is increased.    To conclude, Figure 17 shows the results obtained by adding a parametric uncertainty to: the inertia moments (I add = −10%) and the gains of the SISO models identified for the rotors, ∆K R = (4, −6, 3.2, −2.9)%. As can be seen, the modified LQR/LQG controller and the position controller follow the setpoints demanded in an acceptable way, although they present an increase in overshoot and tracking error than the one obtained in Figure 13 (under the absence of added parametric uncertainty). The negative value chosen for the uncertainty added to the moments of inertia has been fixed to obtain an unfavourable situation. Since, as shown by the indicators in Tables 4 and 5 and in Figure 15, lower inertia moments make the system faster and, therefore, if the control system is not adjusted, its performance decreases.

External Disturbance Robustness Analysis
In this section, the robustness to significant external disturbances of the modified LQR/LQG control proposed in this paper is studied. First, a scenario in which the AEM model is subjected to external torques that increase until a constant value is established, adding to these torques a stochastic component (simulating a continuous wind). Furthermore, in the same scenario, a significant increase in weight is applied to the quadrotor (one-third of its weight) in order to simulate a manoeuvre in which the quadrotor system would elevate a load outdoors for then transport it. The results obtained when carrying out the hover manoeuvre in this scenario are presented in Figure 18.
As shown in the figure above, external disturbances are successfully rejected by the attitude and altitude controller, although the rejection time could be improved if fine-tuning to the integral action of the controllers is carried out. It should be added that the stochastic component of the disturbances increases the oscillation in the controlled variables. Especially in the angles φ and θ, although in this scenario this oscillation can be considered admissible, since it is bounded between [−0.5, 0.5] degrees. The duty cycles generated by the modified LQR/LQG control system and the applied external disturbances are shown in Figure 19.  Figure 18. AEM model of the quadrotor carrying out manoeuvre hover while it is exposed to the set of external disturbances shown in Figure 19. Finally, another scenario in which the quadrotor system is exposed to occasional external disturbances (of significant magnitude) while performing the hover manoeuvre is considered. This scenario can occur in obstacle avoidance and path planning missions, in which SUAV systems are often used. Specifically, the external torques and forces to which the AEM model is exposed could be the result of punctual collisions with different obstacles, as a result of a fault in the avoidance system, see Figures 20 and 21. In addition, in order to create an unfavorable scenario, an uncertainty of I add = −10% has been added to the inertia moments of the AEM model and an uncertainty ∆K R = (−2.5, 2, −2.3, −2.2)% is also added to the gains of the SISO models identified for the rotors.  In this last scenario, the magnitude of the external punctual torques makes that the system deviates considerably (25 and −22.5 degrees in the angles φ and θ, respectively) from the equilibrium position hover, for which the modified LQR/LQG control has been designed. Moreover, the punctual torque on the z B axis causes an undesired turn of 45 degrees on the ψ angle. However, the attitude control system is able to correct the effect of these punctual disturbances (shown in Figure 21) and return the system to its equilibrium point (hover). It should be added that the altitude control also achieves to reject the force maintained for 16 s on the z E axis. To reject this set of disturbances, the controller generates the duty cycles for the PWM signals shown in Figure 21. Therefore, the robustness of the modified LQR/LQG controller has been checked under the presence of significant punctual disturbances and added parametric uncertainty.

First Pre-Tuning Tests on the Real System
In order to validate the pre-tuning method proposed in this work, the initial results obtained when implementing the set of algorithms presented in the real quadrotor DJI F-450 are shown. Note that it is not the object of this work to compare the proposed method with alternative methods; this will be studied in future works. In this section, we only want to show the validity of the proposed method through its implementation in a real system. The results obtained can be improved by fine-tuning. With the limitations of these real tests established, it should be added that the embedded system has been implemented on the PCB whose schematic with the main parts is shown in Figure 22.  Figure 23 shows the behaviour of the multirotor for setpoint changes in altitude (where z sp and z est denote setpoint and estimated variable, respectively). An acceptable average tracking error is obtained for a pre-tuning stage, and the maximum transient error is limited to ±0.07 m.  Figure 24 shows the attitude (φ, θ, ψ) of the quadrotor during takeoff (where the subindex SP and est denote setpoint and estimated variable, respectively). As it can be seen, the modified LQR controller needs a few seconds to get that the quadrotor follows the attitude setpoints demanded. After this transitory, the system follows the demanded setpoints with a tracking error bounded between ±2 degrees. These results are less favorable than those obtained with the simulated model but could be improved by fine-tuning the controller. This is due to the fact that there will always be differences between the real system and the simulation. These differences are mainly due to the uncertainty in the model used for design, as well as possible disturbances acting on the system. Finally, [52] presents a video in which the real quadrotor system is affected by significant external disturbances. As can be seen in this video, the pre-tuning method proposed in this work obtains a robust attitude and altitude controller to significant external disturbances.

Conclusions
A method for pre-tuning a multivariable PID controller for a quadrotor attitude and altitude control has been proposed, which is based on LQR/LQG theory, using a modified structure of the LQR, where the state feedback algorithm employs the tracking error vector instead of the state vector. The proposed method uses only one pre-tuning parameter, which easily facilitates the controller design. Realistic computer simulations have been carried out using a nonlinear mathematical model of the DJI-F450 quadrotor, considering parametric uncertainty and external disturbances acting at the plant input. Madgwick's algorithm and a Kalman filter have been used for data fusion and estimating state variables using data from MARG, LIDAR, and optical flow camera, which have been simulated using data sheets provided by manufactures. Simulation results show that the proposed method design gives MIMO PID controllers which are robust in the face of uncertainties and external disturbances acting at plant input.