Investigation of Multi-Input Multi-Output Robust Control Methods to Handle Parametric Uncertainties in Autopilot Design

Some level of uncertainty is unavoidable in acquiring the mass, geometry parameters and stability derivatives of an aerial vehicle. In certain instances tiny perturbations of these could potentially cause considerable variations in flight characteristics. This research considers the impact of varying these parameters altogether. This is a generalization of examining the effects of particular parameters on selected modes present in existing literature. Conventional autopilot designs commonly assume that each flight channel is independent and develop single-input single-output (SISO) controllers for every one, that are utilized in parallel for actual flight. It is demonstrated that an attitude controller built like this can function flawlessly on separate nominal cases, but can become unstable with a perturbation no more than 2%. Two robust multi-input multi-output (MIMO) design strategies, specifically loop-shaping and μ-synthesis are outlined as potential substitutes and are observed to handle large parametric changes of 30% while preserving decent performance. Duplicating the loop-shaping procedure for the outer loop, a complete flight control system is formed. It is confirmed through software-in-the-loop (SIL) verifications utilizing blade element theory (BET) that the autopilot is capable of navigation and landing exposed to high parametric variations and powerful winds.


Introduction
Quite a few applications require automatic flight control techniques that help or even replace the human pilot [1,2]. Autopilot technologies have advanced significantly over time to present-day autopilots that achieve path following and landing rather proficiently [3][4][5]. Autopilot technologies are no longer reserved to military fighters and large-airlines; coupled with the increased availability of unmanned aerial vehicles (UAVs) to the common researcher, they have found application in many smaller-scale civilian projects [6][7][8]. Procedures utilized for autopilot design comprise of dynamic inversion [9], nonlinear optimal predictive control [10], reconfigurable flight control laws [11], robust nonlinear control [12], Lyapunov vector fields [13], command filtered backstepping [14], sliding mode control [15], multiple model adaptive control [16], invariant manifolds [17], fault tolerant H 1 control [18] and geometric control [19].
Conventional autopilot models usually approximate every flight channel as independent and create single-input single-output (SISO) controllers [2]. To regulate the banking angle for instance, a SISO model from aileron deflection to roll angle is formulated, on which a SISO control technique (e.g. PID) is employed [4]. The operation is continued for every aspect being controlled and all the SISO designs are finally used concurrently in the ultimate design. This method disregards the couplings amongst the channels. However, these side effects are usually not powerful in most fixed-winged airplanes so it continues to remain a dominant approach. Multi-input multi-output (MIMO) methods are usually preferred for missiles [20][21][22][23], helicopters [24][25][26] and multirotor vehicles [27][28][29][30] where dynamical couplings are dominant.
Regrettably, numerous research work hint that there exist threats linked to the standard SISO design method for even fixed-winged aerial vehicles. High crosswinds, aerobatics, structural damage and actuator failures could lead to the aircraft getting forced away from its normal behaviour [31,32]. These kinds of scenarios break the regular weak-coupling suppositions. It is furthermore common in conventional SISO designs to presume that the airplane variables like geometry, mass, inertias and stability derivatives are known completely. In actuality it is hard to come by accurate estimates to these variables. In reports assessing popular methods such as Advanced Aircraft Analysis (AAA), Athena Vortex Lattice (AVL), wind tunnel based modelling and real flight dependent modelling, it has been discovered that one method could possibly produce parameter values quite different from an alternative approach [33][34][35][36]. These parameters get inserted inside the system model during linearization of the nonlinear dynamics about a particular operating condition. A controller constructed on this sort of a model is then utilized for the whole flight envelope. This often works nicely for the nominal case yet performance degrades and destabilization can take place for perturbed conditions [37].
The two primary cases for which departure occurs from the nominal circumstances are leaving the vicinity of the operating point and the existence of uncertainties in the aircraft parameters. Although many results are documented with regards to the former [38,39], it is hard to find a single study on the latter for the complete parameter set on the whole aircraft (i.e. not just particular parameters on particular modes [40]). In this paper a practical oriented study is given to investigate this issue. Via an example on a general aviation airplane, it is proven that the common technique of designing SISO controllers for each channel may have bad performance and even lose stability under modest parametric uncertainties. We then look at two robust MIMO design methods, namely loop-shaping and μ-synthesis, which can endure large parametric variations and at the same time maintain decent performance. Nonlinear simulations along with software-in-the-loop (SIL) verifications utilizing blade element theory (BET) [41], [42] are executed. BET is commonly perceived as an extremely reliable way of carrying out aerospace simulations in which all surfaces of the aircraft are represented as "blade elements" for force and moment calculations.
The scientific contribution of this work can be summarized under three items: 1. Analysis of the critical effects of varying all aircraft parameters simultaneously: In this work, all the parameters related to mass, geometry and the dynamical behavior of the aircraft are varied simultaneously to first uncover the possibility of fatal results with variations as little as 2% under traditional control approaches. These traditional approaches are presently embedded in many of the existing autopilot systems used in civilian and military aircraft.
2. Systematic construction of autopilots to remedy the risks associated with parameter variations: It is demonstrated that MIMO robust control approaches can endure large parametric variations and at the same time maintain good performance for aircrafts.
Step by step construction using two such methods, namely loop-shaping and μ-synthesis, is outlined. It is seen that as much as 30% changes on the full parameter set is acceptable without significant sacrifice in stability and performance.

Separation of design-models and verification-models:
A clear separation between the aircraft model used in design and that used for verification is made for this study. The autopilot designs utilize a mathematical model based on the theory of aircraft dynamics. This model consists of nonlinear differential equations for which the forces and moments are computed using coefficients called stability derivatives. Once the controllers are established, the final closed-loop verifications are done with SIL tests which are based on different and much more accurate models of the aircraft. These models use BET in which all surfaces of the aircraft are subdivided into small regions called blade elements for force and moment calculations. BET flight simulations are widely regarded to produce the most realistic results. Moreover, they are not based on ordinary differential equation models of the aircraft such as the one we use here for designing the controllers. Avoiding the same type of model in the testing phase is advantageous to increase the validity of the proposed methods.
The remainder of the paper is structured as follows: Section Methodology sets out the structure for the study together with background information on necessary tools. Section Results performs the analyses on a general aviation aircraft and presents the outcomes. Section Conclusions wraps up the article with conclusions and future research ideas.

Mathematical Model
The first step is the derivation of the mathematical model on which numerical analysis and controller design will be carried out. The model of the aircraft dynamics is acquired from rigid body force and moment equations: The coefficients of the matrix I are the moments and products of inertia of the rigid body and they are constant for a frame of reference fixed to the aircraft. A rearrangement of Eqs (1) and (2) yields After some manipulations the following non-linear state space system is obtained in terms of which the state space equations can be derived as: x e ¼ ½u cos y þ ðv sin 0 þ w cos 0Þ sin y cos c À ðv cos 0 À w sin 0Þ sin c _ y e ¼ ½u cos y þ ðv sin 0 þ w cos 0Þ sin y sin c þ ðv cos 0 À w sin 0Þ cos c where _ z e ¼ u sin y À ðv sin 0 þ w cos 0Þ cos y and P pp , P pq , P pr , P qq , P qr , P rr , P l , P m , P n , Q pp , Q pq , Q pr , Q qq , Q qr , Q rr , Q l , Q m , Q n , R pp , R pq , R pr , R qq , R qr , R rr , R l , R m , R n are values determined by the inertia values. The equations given listed here are in the East, North, Up (ENU) reference frame. The sign of _ z e can be reversed if North, East, Down (NED) is desired. For solving the abovementioned differential equations one needs to get the force and moment values F = [F x F y F z ] T and M = [L M N] T . These depend on a variety of mass and geometry parameters, the thrust mechanism, along with the control commands. These forces and moments are handily depicted with respect to stability derivatives, which capture the impact of various important variables on a given force or moment value. For example, the longitudinal aerodynamical force can be written as where , C X ad f are the stability derivatives capturing the effect of the term which they multiply on F x . Expressions for F y , F z , L, M, N may be expressed in a similar fashion in terms of their matching stability derivatives [43]. For the implementation of the mathematical model we utilize the Flight Dynamics & Control (FDC) Toolbox for MATLAB [43]. The simulation results produced by this package were tested on real aircraft with successful results. While a different aircraft was used in the FDC Toolbox, actual values for mass, geometry and stability derivatives for Cessna 172 are available in various reliable sources including reports form Cessna Aircraft Company itself [44], making it possible to obtain a valid and reliable Cessna 172 model.
It is difficult to perform a direct validation of the aircraft model based on real flight data since we do not currently have access to a real Cessna 172 and a flight recorder. Nor were we unable to locate any work in literature where such data is openly available. The closest is the study of Neuhart et al. [45] where actual flight data from a Cessna 172 was collected using a custom-built data acquisition system. While the comprehensive data set is not provided numerically, a scenario involving a stall-test is presented through some plots. In this scenario the aircraft was trimmed in a wings-level attitude for 30°flap with a fixed throttle setting. The pilot increased elevator input until stall was achieved. This test matches JAA test 2c8 and FAA test 2c9 [46,47]. Although numerical data for this test is not available, looking at the figures one can eyeball that the aircraft stalls when the elevator input is around 15°.
For validation purposes we make an attempt to replicate the scenario above with our mathematical model as follows: We trim the Cessna 172 model with 30°flap for straight and level flight, after which the elevator input is increased gradually until stall happens. The results can be seen in Fig 1 where it is observed that the aircraft stalls around roughly t = 20 s, when the elevator is at 0.25 rad = 14.324°. This is quite close to the 15°observed in actual flight [45]. This serves as an additional proof of our model's acceptable accuracy, as well as the feasibility and effectiveness of the analysis and control design to be carried out on it.

Trimming and Linearization
Trimming describes the procedure for identifying an operating point for a provided flight condition. In our design, the aim is to solve the following equations where V 0 is the chosen airspeed, z e, 0 is the desired attitude and the time derivatives are calculated from the nonlinear dynamical equations provided in Section Mathematical Model. The trimming procedure may be presented as the subsequent general optimization problem as min x;u f ðx; uÞ ð16Þ subject to Ideally f(x, u) = 0 yet because of numerical implementation and round-off errors f(x, u) < 10 −3 is appropriate in reality. When the problem is formed as an optimization problem as above, it can be solved by means of effective numerical methods like Sequential Quadratic Programming (SQP) [48,49]. The trim point solving the optimization problem in Eqs (16) and (17) is expressed as (x 0 , u 0 ) where x 0 is vector of the aircraft states at the operating condition and u 0 is the vector of control inputs to be applied at the trim condition. The nonlinear aircraft model could then be linearized around the operating conditions (x 0 , u 0 ), which produces a linear state-space system G of the form to be utilized in controller design G : ð18Þ The vectors fields f(x, u) and h(x, u) consist of respectively the equations for derivatives of the states and the outputs to be controlled. The linearized system may also be written in transfer function matrix form where I is the identity matrix.
The first MIMO control strategy utilized in this study is H 1 loop-shaping. In loop-shaping control design, the desired specifications are commonly expressed as where " s and " s denote minimum and maximum singular values respectively. Here, S(s) is the sensitivity function defined as T(s) is the complementary sensitivity function defined as

L(s) is the loop transfer function matrix
LðsÞ ¼ GðsÞKðsÞ; ð26Þ jW À 1 1 ð joÞj is the desired disturbance attenuation factor and |W 3 (jω)| is the largest anticipated uncertainty of the plant expressed as a multiplicative perturbation. Note that the singular values of S(jω) determine the disturbance attenuation since S(s) is actually the closed-loop transfer function from an output disturbance d to plant output y. Note also that T(s) is indeed the closed-loop transfer function on the whole system. A stabilizing H 1 controller K is computed for plant G to make the sigma plot of the loop transfer function GK have desired loop shape G d with accuracy γ. The specifications on disturbance attenuation and multiplicative stability margin in Eqs (22) and (23) can be written in terms of singular values of the loop transfer function since one can make the following approximation for " sðLðsÞÞ ) 1 and the one below for " sðLðsÞÞ ( 1 Therefore if ω c is the 0 dB crossover frequency of the singular values plot of G d (jω), the specifications can be stated as Thus, high tracking performance is achieved at low frequencies where the system model is more accurate, and high robustness is achieved at high frequencies where the system model is less accurate and noise effects are stronger. A stable minimum-phase loop-shaping, squaring-down prefilter W is computed using greatest common divisor (GCD) formulas [50] such that the shaped plant G s = GW is square and that the desired shape G d is achieved with good accuracy in a desired frequency range (ω min , ω max ) by the shaped plant; i.e., Normalized coprime factor synthesis theory is then used to compute an optimal loop-shaping controller for the shaped plant. If the shaped planet is factored as then any perturbed plant can be written as where Δ M and Δ N are stable and unknown transfer functions that represent uncertainties in the nominal plant. The objective of the robust controller design is to stabilize by a controller K, not only the nominal plant but also the family of perturbed plants defined as For robust stability, internal stability must be achieved for the nominal and perturbed plant. If there exist a K such that (M, N, K, ε) is robustly stable, then (M, N, ε) is said to be robustly stabilizable with stability margin ε [51]. For robust stability the following must be satisfied where the infimum is taken over all stabilizing controllers. The H 1 optimization problem allows ε −1 being chosen as small as possible. For actual implementation, the robust stabilization problem can be converted to a more suitable formulation. Let P ¼ : Then (Eq 35) can be seen to be equivalent to where K is gain chosen over all stabilizing controllers and P is a plant of standard form for H 1 optimization problem [52]. The final controller to be used is then computed as μ-synthesis using D-K Iteration Another MIMO control strategy utilized in this study is μ-synthesis, whose goal is to achieve robust performance in the presence of uncertainties. Consider the control system configuration in Fig 2. In the figure the nominal open-loop interconnected transfer function matrix is denoted P(s). This is the aircraft model without any uncertainties. The uncertainties in the parameters are represented by the transfer function matrix Δ(s) and the controller is denoted as K(s). The signals forming the interconnections are named d, v, w, z, u and y as shown in the figure. Based on these interconnections, P(s) can be partitioned as The signals y are the feedback signals to the controller (i.e. tracking errors for our case) and u are the control signals generated by the controller. (i.e. the thrust and the surface deflections for our case.) The blocks P(s) and K(s) can be composed into a single block M as follows: The goal of μ-synthesis it to find a stabilising controller K such that where μ is the structured singular value defined by An iterative approach was introduced in [53] to solve (Eq 44). The method is called the D-K iteration μ-synthesis method, and is based on solving the following optimization problem, for a stabilising controller K and a diagonal constant scaling matrix D, Referring to (Eq 42), a stabilising controller is to be found such that The D-K iteration method minimizes (Eq 45), i.e. reduces the left-hand-side value of (Eq 46) for K and D in turn, while keeping the other one fixed. For a given matrix D, (Eq 45) is a standard H 1 optimisation problem inf KðsÞ k DMðP; KÞD À 1 k 1 ð47Þ that can be written as compatible with the partition of P. For a fixed K(s), inf D " s½DMðP; KÞD À 1 ð joÞ is a convex optimization problem at each frequency ω. After minimization on a frequency range of interest, the resultant diagonal matrices can be approximated using curve fitting, by a stable minimum phase, rational transfer function matrix D(s). This is then used in the next iteration for K. The steps of the D-K iterative μ-synthesis algorithm can be summarized as follows: 1. Start with an initial guess for D, usually D = I.
2. Fix D and solve the H 1 -optimization for K, 3. Fix K and solve the following convex optimization problem for D at each frequency over a desired frequency range, 4. Curve fit D(jω) to get a stable, minimum-phase D(s). Go to Step 2 and repeat, until a desired convergence tolerance is met, (Eq 46) is achieved, or a prespecified maximum iteration count is reached.
It is well known that the solution to the H-infinity optimization problem (Eq 50) is not unique except in the scalar case and there are no analytic formulae for the solutions in general [54]. In practice one usually seeks a suboptimal solution close enough to the actual one, i.e. one tries to find a controller K such that for a small enough value of γ > 0. Several established methods exist for solving (Eq 52) among which we utilize the two-Riccati formulae, the mathematical details of which can be found in [55]. A slight extension is employed here in the sense that multiple iterations are performed with successively smaller values of γ. Starting with a conservative (high) and tight (low) bound guess for γ, a bisection algorithm is applied to approach the optimal γ value. At each step, the problem (Eq 52) may or may not be feasible depending on how small γ is. The algorithm terminates and returns the last feasible solution obtained when the relative difference between the last γ value that failed and the last γ value that succeeded is less than a specified tolerance (0.01 for this work). The D-K iteration approach used in this study is implemented using the numerical computing package MATLAB. It will be employed in the succeeding sections to design a MIMO controller for the aircraft model with a prescribed amount of parameter variation in an attempt to achieve robust performance over the entire uncertainty range.
Remark. At this point it worth recapping why loop-shaping and μ-synthesis were chosen as the control design methods above others. While various tools exists for building robust controllers, loop-shaping and μ-synthesis were preferred since their design procedures can be linked (directly or indirectly) to the control of MIMO systems under parametric variations. In loopshaping, the procedure is based on finding a controller to make the loop transfer function match a desired loop shape. The properties of the desired loop shape such as low/high frequency gains, bandwidth and crossover slope can be used to specify performance and robustness margins. These will in turn determine the behaviour under undesired circumstances including parameter uncertainties. As to the μ-synthesis based on D-K iteration technique, this is a numerical method where the amount of anticipated uncertainty can be specified directly as a constraint for optimization. Once the procedure converges to a solution, the resulting controller is guaranteed work for the modelled uncertainty.

Blade Element Simulation
After nonlinear dynamical simulations, the final test for the control system is software-in-theloop (SIL) verifications based on blade element theory (BET). The surfaces of the aircraft (e.g. propellers, wings, stabilizers) are divided into several sections, the lift/drag forces acting on each section is computed separately and the composite effect is applied to the entire aerial vehicle (Fig 3). This approach contrasts traditional flight simulations relying on empirical data (e.g. stability derivatives) in predefined lookup tables and is widely accepted to be more realistic albeit computationally expensive. It is also a good choice for testing the control design presented here as the mathematical model utilized is based on stability derivatives. Hence a different (and more accurate) flight simulation technique that does not rely on stability derivatives serves as a better test. The main idea of BET can be summarized on a propeller blade shown in Fig 4. The blade is divided into N elements, each of which experiences a slightly different flow. Lift and drag coefficients (C L , C D ) are readily available for numerous airfoil shapes from wind tunnel tests. Using relative velocities, the flow over each element can be related to these tests. The flow is slightly turned passing over the airfoil so inlet and exit flow conditions are averaged to improve accuracy. Carrying out the necessary computations yields where dL, dD, dF x , dF θ are respectively the lift, drag, axial and tangential forces, β is the relative flow angle, ρ is the air density, V is the flow velocity, a is the axial induction factor, r is the radius, σ 0 is the local solidity, B is the number of blades and c is chord length [41]. The procedure is carried out on the entire aircraft to compute all the forces, using which the flight dynamics can be simulated.

Results
In this part first the issues linked to the conventional independent SISO control strategy is highlighted on the attitude control (inner loop) of a popular general aviation aircraft, namely the Cessna 172. The mass, geometry parameters, performance data and stability derivatives of this aircraft are given in Tables 1 and 2. The objective is to use the four usual inputs, namely, thrust F X , elevator deflection δ e , aileron deflection δ a and rudder deflection δ r to manipulate the four outputs airspeed V, pitch angle θ, roll angle ϕ and sideslip angle β. The outcomes of the are conventional individual-channel SISO approach are compared to the MIMO control design methods loop-shaping and μ-synthesis, which drastically enhance the robustness to parameter changes. Next the functionality of the robust MIMO control approach is additionally validated via the inclusion of an outer loop for path following.

Independent SISO Control Design
We consider a typical flight condition for the Cessna 172, namely straight and level flight at V 0 = 65 m/s and elevation z e, 0 = 1000 m. A trim point for this condition is computed as representing the local behavior around V 0 = 65 m/s and z e, 0 = 1000 m. For individual SISO design, since there are four inputs, one can select four of these G ij (s) and design separate controllers. It makes sense to pick the diagonal entries G 11 , G 22 , G 33 , G 44 since for straight and level flight, V is mostly influenced by F x , θ by δ e , ϕ by δ a and β by δ r . For a different scenario an alternate choice may be preferable. For instance during landing many pilots think of regulating airspeed by elevator, altitude by thrust, keeping the aircraft level with ailerons, and aligning with the runway using the rudder. A standard SISO control design is then performed on each diagonal entry, hoping that the off-diagonal dynamics remain well-behaved in closed-loop. To design the individual controllers, several standard automated tuning methods such as Ziegler-Nichols PID, internal model control (IMC), linear quadratic Gaussian (LQG) and optimization based approaches were tested and the best results were obtained for the IMC design method [56,57] with a time constant of τ = 1/3 s % 0.3333 s, resulting in well-damped responses settling in about t s = 4τ % 1.3333 s as seen in Fig 6. While the controllers seem to perform well on their individual nominal models, the real test is whether they will regulate the aircraft attitude successfully when used simultaneously under parameter uncertainties in nonlinear simulations. A sample scenario is presented in Fig 7 where the aircraft is commanded to do the following: 1. Increase airspeed V by 1 m/s at t = 5 s and decrease it by the same amount at t = 20 s.   One hundred Monte-Carlo simulations on the nonlinear model were performed, allowing up to a mere 2% perturbation in each aircraft parameter. The resulting aircraft states are plotted in Fig 7. It can be seen that some of the responses provide acceptable tracking of the references; these are for parameters close to nominal values listed in Table 1. However, it is clear from the figure that certain parameter combinations have poor performance with significant oscillations, some of which even grow unboundedly. In practice this would imply the destabilization and potential loss of the aircraft. It is true that these scenarios are unlikely and constitute a small portion of all the runs. However any risk of losing an aircraft and the lives of those onboard justifies the need to consider alternate autopilot strategies. One such strategy is investigated in the next section.

Loop-shaping MIMO Control Design
For loop-shaping control design the system is first linearized around the nominal flight condition V 0 = 65 m/s and z e, 0 = 1000 m for the parameters given in Tables 1 and 2. The procedure is identical to that described in previous section until (Eq 59). Simple actuator dynamics are also augmented to G(s) for a more realistic model, i.e. G nom ðsÞ ¼ GðsÞG act ðsÞ ð60Þ and diag stands for diagonal matrix. G act captures the fact that the real Cessna 172 reacts slower to throttle command, somewhat faster to elevator and rudder commands, and the fastest to aileron commands. Therefore the controller must not generate commands beyond these bandwidth or else they will not be effective. The loop-shaping is controller designed following the procedure in Section H 1 Loop-shaping with desired loop-shape  Investigation of MIMO Robust Control Methods to Handle Parametric Uncertainties in Autopilot Design so the closed loop system will be able to track all references successfully with minimal overshoot and a settling time of approximately t s ¼ 5t ¼ 5 1 2 ¼ 2:5 seconds. Also the off-diagonal entries of T(s) are roughly zero, which indicates that the coupling between different commandresponse pairs are eliminated.
The results for the loop-shaping controller in closed-loop for the nonlinear aircraft with the nominal parameter values are shown in Fig 9. The nominal performance is seen to be quite good with accurate tracking and very little cross-coupling. The inputs generated by the controller are also shown in Fig 10. The thrust and surface deflections remain within reasonable limits at all times and the control inputs do not contain significant power at frequencies higher than about 0.5 Hz. This ensures that the control does not cause any sharp thrust changes or wild oscillations in control surfaces. The responses of the loop-shaping controller to perturbed nonlinear models are shown in Fig 11. One hundred Monte-Carlo simulations on the nonlinear model were performed, allowing up to 20% perturbation in each aircraft parameter. It is seen that the stability is never lost, and acceptable tracking performance is achieved for all cases. The inputs to the aircraft are also plotted in Fig 12 showing that the inputs always remain within reasonable amplitude and frequency ranges.

μ-synthesis MIMO Control Design
For μ-synthesis control design the system must be expressed in the form described in Section μ-synthesis using D-K Iteration. For this purpose, the nominal aircraft model G nom (s) is defined based on G(s) in (Eq 59), i.e. the linearization around V 0 = 65 m/s and z e, 0 = 1000 m for the parameters given in Tables 1 and 2. We also augment simple actuator dynamics for a more realistic model, i.e.

G nom ðsÞ ¼ GðsÞG act ðsÞ ð66Þ
where and diag stands for diagonal matrix. G act captures the fact that the real Cessna 172 reacts slower to throttle command, somewhat faster to elevator and rudder commands, and the fastest to aileron commands. Therefore the controller must not generate commands beyond these bandwidth or else they will not be effective. Parameter variations are expressed as an input multiplicative uncertainty Δ M (s) so that the real aircraft model is The goal is set to design a controller that can tolerate up to 20% uncertainty on the aircraft parameters; this figure is consistent with the typical variation in parameter values estimated by AAA, AVL, wind-tunnel tests and actual flight data as described in the Introduction. To compute Δ M (s) the parameters are randomly perturbed up to 20%, a linearized model is obtained for each perturbed parameter set, and the procedure is repeated 100 times to collect enough data. The variations between the nominal model G nom and all 100 models are computed and an uncertainty Δ M (s) is selected to cover the maximum variation over the frequency range of interest, which is chosen as ω 2 (0.01, 100) rad/s for this application. Due to the large data set, it would be time consuming and computationally intensive to estimate individual covers for all 16 channels so we utilize a scalar third order Δ M (s) to cover the uncertainties over all channels:  Fig 14 is constructed to carry out μ-synthesis control design. The main goal is to match the response of the system to a desired response contained in the block G des (s), i.e. to minimize the error e des . Similar to the independent SISO design case, we ask for the diagonal channels to be well damped with a time constant of τ = 1/3 s = 0.333 s. Unlike the independent SISO design however, we have the opportunity to specify the desired response for the off-diagonal entries, which would ideally be zero for ð69Þ W e and W u are weighting filters for the error and inputs to select certain frequencies during optimization for μ-synthesis. It is natural to emphasize low frequencies for tracking purposes so a first order filter with DC-gain 50, bandwidth 3 rad/s, and high frequency gain 0.02 is selected The inputs generated by the controller must adhere to the frequency limitations of the actuators so we pick W u = G act to represent this criterion. Establishing correspondence between Figs 14 and 2, Δ is Δ M , K is the block labeled "Controller" and P contains everything else. The external inputs are w = [r n] T , where r are the reference commands, and n is an output disturbance representing sensor noise, atmospheric effects and so on. The overall outputs of the system are z = [e filt u filt ] T , which are the weighted error and weighted input signals. The system is therefore in suitable for for D-K iteration described in Section μ-synthesis using D-K Iteration so the procedure is executed numerically using MATLAB to design the controller. The procedure returned a K satisfying sup ω μ[M(P, K)(jω)] < 0.9174 indicating a successful design in accordance with (Eq 42).
The results for the μ-synthesized controller in closed-loop for the nonlinear aircraft with the nominal parameter values is shown in Fig 15. The nominal performance is seen to be quite good with accurate tracking and very little cross-coupling. The inputs generated by the controller are also shown in Fig 16. The thrust and surface deflections remain within reasonable limits at all times and the control inputs do not contain significant power at frequencies higher than about 0.5 Hz. This ensures that the control does not cause any sharp thrust changes or wild oscillations in control surfaces. The responses of the μ-synthesized controller to perturbed nonlinear models are shown in Fig 17. One hundred Monte-Carlo simulations on the nonlinear model were performed, allowing up to 20% perturbation in each aircraft parameter. It is seen that the stability is never lost, and acceptable tracking performance is achieved for all cases. The inputs to the aircraft are also plotted in Fig 18 showing that the inputs always remain within reasonable amplitude and frequency ranges.
Remark. While the results above confirm robust performance and stability to the designed amount of uncertainty (20%), the controller in practice is able to handle much higher perturbations. We performed over 1000 simulations where we have seen that the performance remains acceptable to around 34%, and stability is not lost up to around 43% parameter variation. This is consistent with the knowledge that μ-synthesis control design via D-K iteration can be suboptimal and produce conservative results.

Outer Loop and Blade Element Simulation Tests
As a final test, an outer-loop controller is wrapped around the attitude controller to generate attitude references (r in Fig 14) from outer-loop references r outer = [r v r z e r χ r β ] T . These are the  airspeed command, altitude command, heading angle command and sideslip angle command respectively. The outer and inner loop collectively form a complete flight control system capable of navigation and landing. In this section the results are presented for both loops designed with loop-shaping but similar outcomes were also obtained using μ-synthesis.
The initial test is performed using pulse-type references both on the nominal model and on perturbed models with up to 30% uncertainty. Plots in Figs 19-22 indicate good performance for this situation. Finally, SIL simulations were performed with the flight controller implemented in MATLAB/Simulink and the blade element simulations carried out by the flight simulator X-Plane. The former sends thrust, elevator, aileron, rudder commands to the latter and receives flight simulation results in real-time every 25 ms through user datagram protocol (UDP) packets (Fig 5). Numerous scenarios with up to 30% perturbation were studied in this configuration with success and one is presented here as an example. In this case a perturbed Cessna 172 cruising at 65 m/s and 1000 m receives appropriate commands for maneuvers required to navigate to a target airfield, line-up with the runway, descend, flare and touchdown. r β is always zero for turn coordination and flight comfort. For the first 500 seconds the weather remains calm but after t = 500 s, the wind magnitude is steadily increased, reaching 15 m/s during landing phase with the crosswind component reaching 7 m/s. Some oscillation in the states are unavoidable for high winds; nevertheless the flight control system responds well to the commands as seen in Fig 23. The control inputs applied to the aircraft and their frequency spectra is presented in Fig 24, which are within practical limits. The 3D trajectory is also visualized in Fig 25. Overall the flight control system maintains stable flight, responds well to commands received and executes the landing, even under unfavorable atmospheric conditions.
Remark. The careful reader might notice that the thrust F x drops to zero around t = 800 s in Fig 24. This corresponds to the throttle being cut, which is a normal event within the landing phase. When the aircraft is higher than desired during landing, reducing/cutting the throttle is typical procedure; conversely the throttle is boosted if the is the aircraft is too low. Also recall that strong winds are present during the landing phase. Thus the aircraft can be suddenly pushed upward or downward by wind gusts, calling for a rapid correction from the autopilot. A similar reasoning can be made for the other control surfaces. It should be noted however that the long duration for the simulation (about 1000 s) compresses the plotted data and gives the false illusion that the fluctuations are quite wild. Observing the frequency responses in the bottom row of the figure reveals that the control inputs do not contain significant power at frequencies higher than about 0.5 Hz. This means that the autopilot does not cause any sharp thrust changes or wild oscillations in the control surfaces.

Conclusions
Two multi-input multi-output (MIMO) control design approaches were investigated to handle parametric uncertainties in autopilot design for aircrafts. In real-life it is impossible to perfectly determine geometry, mass and stability derivative parameters so some level of uncertainty is inevitable. The study was carried out on the aircraft dynamics as a whole and not on individual modes, which is essential to capture all the coupling effects. It was revealed that an attitude controller working perfectly on individual nominal models can lose stability with a perturbation as small as 2%. Robust MIMO design using loop-shaping and μ-synthesis were presented as remedies, which were seen to withstand high parametric variations of 30%, while retaining good performance. As a final test, the attitude controller was augmented with an outer loop controller designed also using loop-shaping, forming a complete flight control system. It was confirmed through software-in-the-loop (SIL) verifications using blade element theory (BET) that the autopilot is capable of performing navigation and landing under high parametric variations and strong winds.