Application of Structured Robust Synthesis for Flexible Aircraft Flutter Suppression

This article presents a method for structured robust control design for systems with a mixture of parametric and dynamic uncertainty. The proposed method alternates between an analysis step and a synthesis step. Samples of the parametric uncertainty are computed during the analysis steps, thus yielding an array of uncertain systems containing only dynamic uncertainty. The controller is then synthesized on this array of uncertain models. This synthesis step itself involves an alternation between constructing a D-scale for each of the uncertain systems and tuning a single controller for the entire collection of scaled plants. The controller tuning is performed using structured control design techniques. The proposed method is utilized to design a flutter suppression controller for a flexible aircraft. The aircraft dynamics are described by both a high-fidelity and a reduced-order model. The design objectives for flutter suppression are to achieve robust stabilization in the presence of mixed uncertainty. The proposed structured design method yields a single, low-order, linear time-invariant (LTI) controller, which increases the flutter speed by 15%. Additional robustness analyses and high-fidelity simulations are provided to assess the controller performance.

Abstract-This article presents a method for structured robust control design for systems with a mixture of parametric and dynamic uncertainty.The proposed method alternates between an analysis step and a synthesis step.Samples of the parametric uncertainty are computed during the analysis steps, thus yielding an array of uncertain systems containing only dynamic uncertainty.The controller is then synthesized on this array of uncertain models.This synthesis step itself involves an alternation between constructing a D-scale for each of the uncertain systems and tuning a single controller for the entire collection of scaled plants.The controller tuning is performed using structured control design techniques.The proposed method is utilized to design a flutter suppression controller for a flexible aircraft.The aircraft dynamics are described by both a high-fidelity and a reduced-order model.The design objectives for flutter suppression are to achieve robust stabilization in the presence of mixed uncertainty.The proposed structured design method yields a single, low-order, linear time-invariant (LTI) controller, which increases the flutter speed by 15%.Additional robustness analyses and high-fidelity simulations are provided to assess the controller performance.

I. INTRODUCTION
T HERE are a variety of optimal control synthesis methods available in the literature, e.g., H 2 and H ∞ methods [1], [2].However, the practical application of these methods can be limited because they provide full-structure state-space controllers.Industrial applications often prefer or require controllers with specific structure for implementation.Examples include the Vega launch vehicle [3], automated landing of a civilian aircraft [4], and chemical process control [5].
The H ∞ synthesis conditions can be expressed as convex constraints when the plant is nominal (no uncertainty) and the controller is allowed to be full-order and unstructured [6].Nominal synthesis becomes nonconvex, in most cases, once structure or reduced order is imposed on the controller.One notable exception is the quadratic invariance condition [7], which allows for convex synthesis.As a consequence, most research on structured control synthesis has focused on algorithms that are reliable and computationally efficient on typical engineering problems (but not necessarily with proofs of global optimality).
Two notable algorithms for structured control synthesis include the H ∞ /H 2 fixed order optimization (HIFOO) in [8] and structured H ∞ design in [9].This article focuses on the algorithm in [9] due to its ease of use in MATLAB [10] as the hinfstruct function.The approach has also been extended to handle multiple frequency-domain design requirements in the fixed structured design as implemented in the systune function [11].
This article builds on a previous conference paper [12] in proposing a structured synthesis method for systems with mixed (i.e. parametric and dynamic) uncertainty.The chosen objective of the design is the minimization of the worst case gain, i. e. the maximal closed-loop H ∞ norm over the set of allowable uncertainties.The algorithm utilizes an iteration in which controller synthesis is followed by analysis.This iteration continuously updates a set of "bad" samples for the parametric uncertainty.This yields a corresponding set of sampled design plants.The synthesis step itself is an iterative process, which resembles the traditional D-K iteration [13].The K-step is performed as a structured design (using systune) on the collection of sampled plants.The D-step solves for the scalings associated with each sampled plant.In the analysis step, the worst case gain of the closed loop is calculated along with a worst case uncertainty.A sample of parametric uncertainty is extracted from this worst case uncertainty.The subsequent synthesis step is repeated using this updated set of sampled plants.
The proposed algorithm resembles the hybrid relaxation approach in [14, Algorithm 2].The novelty of this article lies in the control synthesis step.In particular, the hybrid relaxation given in [14] has a synthesis step that involves parameterizing a single dynamic D-scale for all uncertainty sample.It then jointly optimizes over the tunable parameters in both the scaling and controller using structured synthesis techniques [9], [10].This approach is effective for many 1063-6536 © 2021 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
See https://www.ieee.org/publications/rights/index.html for more information.problems, but the restriction to a single dynamic scaling might be conservative on certain examples, e.g., see Section IV-B in [12].As commented in [14], it is possible to extend the hybrid relaxation to parameterize one scaling for each parametric uncertainty sample at the expense of a larger number of optimization variables.In contrast, the approach given in this article utilizes dynamic scalings that depend on the parametric uncertainty sample.The use of parameter-dependent scalings reduces conservatism at the expense of requiring a coordinatewise D-K iteration for the synthesis.This coordinatewise iteration alternates between optimizing the D-scaling while holding the controller fixed and vice versa.This tends to be less effective than the joint scaling and controller optimization used for the synthesis in the hybrid relaxation approach.If the D-scales depend strongly on the value of the uncertain parameters, however, the individual D-scale construction has the advantage.This is illustrated with a concrete example in [12].Another minor distinction is that we use a branch-and-bound implementation for the worst case gain analysis.The branch-and-bounding provides tighter analysis bounds.MATLAB implementation of the algorithm is available online [15].
To demonstrate the effectiveness of this synthesis method, and as an extension of our previous work in [12], a real-life flutter suppression control example is presented for the demonstrator aircraft of the FLEXOP project [16].The aeroelastic flutter is an adverse interaction between the aerodynamics and the structural dynamics of the aircraft, which causes undamped oscillations that can result in catastrophic structural failure.There are a number of active control solutions to mitigate flutter.The robust control framework was successfully applied in [17] and extended to the linear parameter-varying (LPV) case in [18].Optimal blending of the inputs and outputs is utilized in [19] to maximize control effectiveness for the flutter modes and minimize the excitation of the remaining modes.
Our approach is similar to [20] with additional focus on robustness.The design model is an uncertain control-oriented model of the aircraft with both parametric and dynamic uncertainties.The control problem is decomposed into the synthesis of two independent single-input-single-output (SISO) controllers with the performance criteria of robust stabilization and control effort minimization.The robustness of the flutter control loop is evaluated using the high-fidelity model of the aircraft.Loop margins are computed and time-domain simulations are conducted to investigate the extension of the safe flight envelope.
The remainder of this article is laid out as follows.The synthesis method is described in detail by Section II.Section III presents the nonlinear high-fidelity model of the aircraft, the bottom-up modeling that results in the control-oriented model, the construction of the uncertain model, and the control design.The analysis of the closed flutter loop is given in Section IV.Finally, conclusions are drawn in Section V.

II. STRUCTURED ROBUST CONTROL DESIGN ALGORITHM AGAINST MIXED UNCERTAINTY
The majority of this section is the restatement of the method in [12].The mathematical formulation of the design problem is given in Section II-A.Then, Section II-B provides an overview of the algorithm.The construction of the D-scales, and the synthesis and analysis steps are elaborated in Section II-C, II-D, and II-E, respectively.
A. Problem Formulation 1) Notation : Let N ∈ C r N ×c N and K ∈ C r K ×c K be given matrices with r K < c N and c K < r N .Partition N such that the lower right block is c K × r K and The lower linear fractional transformation (LFT) is Similarly, the upper LFT of N and ∈ C r ×c is defined as assuming compatible dimensions.The same definitions for both the lower and upper LFTs hold if N, K , and are linear time-invariant (LTI) systems.Finally, diag{ 1 , 2 , . . ., N } denotes the block diagonal concatenation of 1 , 2 , …, N .

2) Structured Robust Synthesis :
The robust synthesis problem is formulated using the feedback interconnection shown in Fig. 1.The uncertain plant P( ) = F U (M, ) consists of an LTI system M and structured LTI uncertainty .The objective is to design an LTI controller with fixed structure, K (κ), where κ ∈ R n κ is a vector of tunable design parameters.
The uncertainty is assumed to be block-structured consisting of unit norm-bounded parametric and dynamic uncertainty [21], [22].Specifically, define the following sets of parametric and dynamic uncertainties: where ||•|| ∞ denotes the H ∞ norm.The dynamic blocks in d are assumed to be square.This assumption can be relaxed with only notational changes.The plant uncertainty is given as ∈ where Note that || || ∞ ≤ 1 for all ∈ .We also define a set of D-scales D d associated with the set of dynamic uncertainties d .The set D d consists of stable, minimum phase, and invertible LTI systems that commute with the elements of This set of scalings will be used in the structured robust synthesis algorithm.
Define the worst case gain of the closed loop in Fig. 1 as This definition assumes that F L (P( ), K (κ)) is robustly stable, i.e., stable for all ∈ .If the closed loop is not robustly stable, then J (κ) is defined to be +∞.The synthesis   objective is to tune the parameters of the structured controller to minimize the closed-loop worst case gain, that is According to our engineering experience, optimizing the worst case gain is a better performance objective than the usual robust performance [13].If robust stability is reached, the design objective is to improve the performance without extending the stability margin any further.Structured robust synthesis is, in general, a nonconvex optimization problem.Sections II-B and II-E describe an iterative method to compute suboptimal controller parameters.

B. Overview of the Algorithm
This section introduces the proposed algorithm for computing a structured controller to minimize the closed-loop worst case gain.An overview summary of the proposed approach is given in Algorithm 1.The algorithm is iterative and utilizes the sample sets for the parametric uncertainty p,s and D-scales D d,s .The set p,s is initialized with the nominal value of the real parametric uncertainty (Step 2).These sets are updated throughout the iteration, as described further in the following.The algorithm also expects that a controller K (κ init ) is available for initialization.It is further assumed that this controller robustly stabilizes the system for the dynamic uncertainty ( d ) and the nominal value of the parametric uncertainty ( p = 0).In some problems, the open loop is robustly stable, and hence, K (κ init ) = 0 can be used for initialization.Otherwise, an initial robustly stabilizing controller can be computed with a related stability margin maximization algorithm.
The first key step of Algorithm 1 is to synthesize a structured robust controller using an iteration (Steps 4-7) that is similar to the (unstructured, full-order) D-K iteration [13].This Algorithm 1 Structured Robust Synthesis alternates between a D-step (Step 5) that computes a set of scalings D d,s and a K-step (Step 6) that computes a structured robust controller K (κ).This iteration terminates when γd does not decrease significantly compared to the previous iteration or when the maximum number of iterations is reached.Our specific implementation terminates the D-K inner loop if γd fails to decrease by more than 1% or the maximum number of 15 iterations is reached.
The D-step (Step 5) first constructs a set of closed loops with the structured controller K (κ), and each sample of the parametric uncertainty p ∈ p,s obtained in Step 8.Each closed-loop sample still has the remaining dynamic uncertainty d .Worst case gain analyses are performed to compute a scaling D d ∈ d for each closed loop with the corresponding sample p ∈ p,s .Thus, D d,s has the same number of elements as p,s .Details on this step are given in Section II-C.
The K-step (Step 6) forms a collection of scaled plants from the samples of parametric uncertainty and D-scales.Specifically, a scaled plant is constructed for each sample ( p , D d ) ∈ p,s ×D d,s , as shown in Fig. 2.Then, systune is used to compute a structured controller K (κ) to minimize the worst case gain on this collection of sampled, scaled systems.Details on this step are provided in Section II-D.Finally, an analysis step (Step 8) is performed on the resulting closed loop F L (P( ), K (κ)).In the analysis, the entire uncertainty is considered using a version of wcgain with additional branch-and-bounding.This generates lower and upper bounds on the worst case gain.It also computes a worst case uncertainty that achieves the lower bound.The parametric block of this worst case uncertainty is added to the sample set p,s (Step 9).Details on the analysis are given in Section II-E.
The algorithm outer-loop terminates (Step 3) if the worst case synthesis gain is within the interval of the analysis bounds, i.e., γ s ∈ [γ a , γa ], and γa fails to decrease significantly (more than 1%) over the last iteration.The loop also terminates if the maximum number of iterations counts ( 30) is reached.
The outline of Algorithm 1 is inspired by [14,Algorithm 2].The key difference between the two approaches is the control synthesis in Steps 4-7.Specifically, in Algorithm 1, a separate D d ∈ D d,s corresponds to the individual samples in p,s , instead of a common D d .This solution reduces conservatism and decreases the number of decision variables for the structured synthesis.This comes at the expense of having to use an iterative process for the controller synthesis, however.Another distinction between the two approaches is that we use branchand-bounding for the worst case gain computation, which results in tighter bounds.
The proposed algorithm is designed to handle mixed uncertainty.If the system only contains dynamic uncertainty ( = d ), then only Steps 4-7 are performed.On the other hand, if the uncertainty is purely parametric ( = p ), then the samples are LTI systems without uncertainty.In this case, Steps 4-7 are replaced by a single controller synthesis step, where we find a controller that simultaneously stabilizes all the samples and minimizes performance at the same time.This is done using hinfstruct [9], [10].

C. Parameter-Dependent D-Scales
Let a parametric uncertainty sample p ∈ p and controller K (κ) be given.The corresponding closed loop T ( p , κ), shown in Fig. 2, is defined as The worst case gain for the closed loop over the remaining dynamic uncertainty is given by It follows from standard μ analysis results [22], [23] that an upper bound on this worst case gain over d is given by min γ subject to: Fig. 2 shows the sampled closed loop with D-scales and performance gain scaling.For the remainder of the derivation, the notational dependence of T on ( p , κ) is dropped for Fig. 3. Gain of the scaled system with D-scales obtained using different methods.
simplicity.Form the partitioning T =: where T 2 has the row dimension of the error signal e. Recall that, for any complex matrix L, it follows that σ (L) ≤ 1 if and only if I − L * L 0. Here, 0 means positive semidefiniteness of the left-hand side, and L * is the conjugate transpose of L. This can be used to express the constraint in (10) as a (frequency-dependent) linear matrix inequality (LMI).By the Schur complement lemma, this constraint is equivalent to ∀ω, where X( j ω) , with each d i as a stable, minimum-phase, invertible SISO system.Thus, X has the structure X = diag x 1 I, . . ., x N d I , where Express each x i as x i = ψ * i ψ, where ψ is a (fixed) column vector of stable, minimum phase basis functions, and i = * i are decision variables to be optimized.
For the subsequent synthesis step, it is advantageous to find a D-scale that does not only minimize the peak of γ but minimizes it over a broad frequency range.To illustrate this, consider Fig. 3.If we solve (10) frequency-by-frequency while allowing the D-scales to vary arbitrarily between each point, then we obtain the theoretical worst case gain lower bound.This is the dashed blue curve in Fig. 3.If we search for a realizable D d (s) by enforcing (11) only at the frequency of the peak, then γ (ω) we get is in general only accurate at the peak and can be far from the theoretical value at other frequencies.See the red dashed-dotted line in Fig. 3. Using such a D-scale leaves little room for the control synthesis to improve the performance further.Therefore, it is better to enforce (11) at several frequencies at the same time.In practice, this calculation usually leads to some inaccuracy at the peak, but it provides lower γ values at the rest of the frequencies, as illustrated by the continuous yellow curve in Fig. 3.
Therefore, the constraint in ( 11) is enforced on a frequency grid {ω k } N ω k=1 .The cost function to be minimized is turned into Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.
where γ k is the gain at ω k and γ is the peak gain over all frequencies, that is is added to the constraints.
To ensure that x i ≥ 0 for all ω, additional Kalman-Yakubovich-Popov (KYP) lemma constraints are added [24].This guarantees that the solution X obtained from the optimization can be spectral factorized to get a stable, minimum phase scaling for all ω.By the KYP lemma, this condition is equivalent to the existence of The optimization in ( 10) is reformulated as the minimization of the cost function in (12) with constraints (11) for each ω k and γ k , ( 13) and ( 14).This is a convex SemiDefinite Program (SDP) in the variables γ , Finally, this optimization is solved to compute dynamic scaling D d for each uncertainty sample p ∈ p,s .
An alternative formulation of this optimization is possible in which the minimization of the cost function in ( 12) is decomposed into two steps.First, only the peak γ is minimized.Second, (13), where ε is a positive tolerance parameter.This approach may yield better results for some problems.Specifically, this alternative algorithm performed better on 13 of the 31 test examples (with no difference in 8 cases).However, the computation time is roughly doubled for the alternative algorithm since two optimization steps are performed.

D. Synthesis Step
The output of the D-step is a set of dynamic scalings D d,s each computed for a corresponding element of the parameter uncertainty set p,s .The synthesis step optimizes the controller K (κ) to minimize the worst case gain over the set of scaled, closed-loop samples.This is formulated as the following optimization: min κ∈R nκ γ >0 γ subject to: This control synthesis is directly addressable using the systune command in MATLAB [11].If K (κ init ) is not robustly stabilizing for the dynamic uncertainty, as assumed in Algorithm 1, the D-scales are set to identity in Step 5.Under these conditions, the optimization in ( 15) may still yield a solution (and in fact usually does).If there is no solution however, a different algorithm is invoked, which maximizes the stability margin of the system, thus providing a robustly stabilizing controller if possible.This algorithm is not detailed in this article for lack of space.

E. Analysis Step
The analysis step computes the upper and lower bounds on the worst case gain of the closed-loop given in (6) for a fixed controller K (κ).This is performed using the wcgain function in MATLAB but modified to incorporate branch-and-bound on the parametric uncertainty.The lower bound γ a is computed using a combination of a (skewed-μ) power iteration for the dynamic uncertainty blocks and a Hamiltonian-based coordinatewise iteration on the parametric uncertainty.This returns a worst case uncertainty that achieves the lower bound.The upper bound γa is computed using a frequency-gridded SDP with D-scales for the dynamic uncertainty and (D,G)-scales for the parametric uncertainty.Thus, instead of using the D-scales obtained during the controller synthesis like the algorithm in [14], we recompute the D-scales on a frequency grid for the analysis step.This is less conservative since the state-space realization of the D-scales is not required for the calculation of the worst case gain upper bound.
Branch-and-bounding of the real uncertainty is used to reduce the gap between the upper and lower bounds.Specifically, the real parameter uncertainty set p is described by a normalized hypercube.This hypercube is split and the upper/lower bounds are computed on each subcube.The splitting of subcubes continues until the gap between the overall upper and lower bounds is below some relative error or a maximum number of cube splits is reached.Details on this analysis step are given in [25].
When the worst case gain upper bound γa is infinite, the worst case uncertainty is chosen, which corresponds to the upper bound of the stability margin of the closed loop.This is computed using the robstab function in MATLAB [26].

III. FLUTTER SUPPRESSION CONTROL DESIGN FOR THE FLEXIBLE AIRCRAFT
The demonstrator of the FLEXOP project, shown in Fig. 4, was built specifically for flutter control experimentation [16].Fig. 5. Positions of the sensors and control surfaces used for the flutter suppression control.The control inputs and measurements are marked at the corresponding control surface or sensor.The signals q L , q R , and q denote angular rates along the y-axis, and u f,L and u f,R are deflection commands for the actuators.
The single-engined aircraft features a wingspan of 7 m, aspect ratio of 20, and takeoff weight between 55 and 66 kg.This section is dedicated to the application of the structured synthesis method to the flutter suppression control design for this platform.In Section III-A, the high-fidelity model is given, which is used in the analysis of the closed loop.The reduced-order control-oriented model is described in Section III-B.Finally, the formulation of the uncertain design model and the control problem along with the synthesis results are in Section III-C.

A. Nonlinear High-Fidelity Model of the Flexible Aircraft
Each wing of the FLEXOP aircraft is equipped with four control surfaces [27] with the outermost pair dedicated to flutter suppression, as shown in Fig. 5.To actuate these, a custom-made direct-drive actuator is developed with bandwidth greater than the flutter frequencies.Based on system identification, the direct-drive actuator has 0.1-ms delay and transfer function In addition to the GPS and air data probe, the aircraft features inertial measurement units (IMUs) at the center of gravity and in the wings, as shown in Fig. 5.The IMUs provide acceleration and angular rate measurements along all three body axes.
The construction of the aeroservoelastic (ASE) model is based on a subsystem approach, which involves the integration of aerodynamics, structural dynamics, and flight dynamics [28], [29], as shown in Fig. 6.The components in Fig. 6 are developed separately and then combined to form the ASE model.
The structural dynamics are modeled by the high-fidelity finite-element method [30] and are then condensed with Guyan reduction [31].The aerodynamics model of the aircraft is based on the vortex-lattice method (VLM) and the doublet-lattice method (DLM), and it is complemented by results from computational fluid dynamics (CFD) simulations.The nonlinear equations of motions are derived based on a mean axes reference frame [32].The mean axes approach Fig. 6.ASE subsystem interconnection.F aero represents the aerodynamic forces acting on the rigid body dynamics and F external represents external, i.e., the propulsion and gravitational forces; the rest of the variables are defined in the remainder of the section.
describes the dynamics of the flexible body by a set of equations that decouple the rigid body modes from the vibrational modes.The mean axes coordinates ensure that the coupling is restricted to external forcing terms only [32].
The nonlinear ASE model of the FLEXOP aircraft consists of 12 rigid body states, 100 states corresponding to flexible modes, and 1040 aerodynamic lag states in addition to the actuator dynamics.This model is considered as the high-fidelity model.Based on this model, the aircraft has two unstable aeroelastic modes.The symmetric flutter mode becomes unstable at 52 m/s and 50.2 rad/s, and the asymmetric mode becomes unstable at 55 m/s and 45.8 rad/s.

B. Bottom-Up Model Reduction
The high-fidelity model in III-A has over 1000 states.Control design for such a high-dimensional model is not practical; therefore, an appropriately low-order and numerically tractable control-oriented model is required [17].Because the size of this model poses considerable difficulty to the LPV reduction techniques [33], [34], the "bottom-up" modeling approach in [35] and [30] is applied in this article.The key idea is to reduce the components of the system in Fig. 6 separately.Since the structural and aerodynamics subsystems have a simpler structure than the combined ASE model, it is possible to simplify them using more straightforward reduction techniques.This approach leads to a sufficiently low-order ASE control-oriented model.
The outermost control surface pair is used for flutter suppression, as shown in Fig. 5.The actuating signal is the control surface deflection command received by the actuator denoted as u f,L and u f,R for the left and right wings, respectively.The sensors used for flutter control are three IMUs: one in the center of gravity and two at the 90% of the half wingspan on each wing.Only the angular rate measurements along the y-axis of the three IMUs are used (see Fig. 5).The signals of the IMUs on the left and right wings are denoted by q L and q R , respectively, while q denotes the pitch rate in the center of gravity.Only these inputs and outputs are considered in the reduction.
The frequency range of interest where the control-oriented model is expected to provide a good approximation of the high-fidelity model is [0, 100] rad/s; 100 rad/s is roughly twice the flutter frequencies 50.2 and 45.8 rad/s, which ensures the accurate representation of the flutter behavior.It is possible to retain the accuracy of the control-oriented model over a wider frequency range at the expense of additional dynamics.
The objective of the reduction is to decrease the number of states of the subsystems in Fig. 6 while maintaining an acceptably low ν-gap between the high-fidelity and the control-oriented model in the frequency range of interest ([0, 100] rad/s).The ν-gap metric δ ν (•, •) is used since it considers the feedback control objective.It assumes values between zero (for identical systems) and one.If a feedback controller stabilizes the system P 1 (s) with stability margin ε, then it also stabilizes P 2 (s) if δ ν (P 1 (s), P 2 (s)) < ε.A plant at a distance greater than ε from P 1 , on the other hand, will in general not be stabilized by the same controller [36].The ν-gap between P 1 (s) and P 2 (s) at the fixed frequency ω is where σ[•] denotes the largest singular value.The aerodynamic lag terms assume the state-space form where V TAS is the true airspeed (TAS), c is the reference chord, x rigid is the rigid body state, η is the state of the structural dynamics, and δ cs is the control surface deflection.Balancing transformation is applied for the aerodynamics model given by A lag , B lag , and C lag in (18).The order reduction is achieved by residualizing the states with the smallest Hankel singular values.Keeping two lag states results in acceptable accuracy.Coefficient C Q 3η 1 for the modal generalized force (see Chapter 7 of [32] for more details) affecting the symmetric flutter mode and coefficient C l p affecting the asymmetric flutter mode of the control-oriented model was scaled by 0.65 and 1.6, respectively.By this heuristic modification, the resulting control-oriented model matches the flutter speeds and frequencies of the high-fidelity model better.The effect of this modification on the rest of the dynamics is negligible.
The structural dynamics of the aircraft are of the form where M, C, and K are the modal mass, damping, and stiffness matrices, respectively, and F modal is the external excitation in modal coordinates.The structural dynamics model is an LTI system; therefore, state truncation is applied.Along with the first six structural modes, modes 19, 20, and 21 are retained since their removal results in a large increase in the νgap between the high-fidelity and the control-oriented model.In this way, the 100th order structural dynamics model is reduced to 18 states.It is assumed that the structural dynamics model has parametric uncertainty.Specifically, the first six modes of the control-oriented model have ±1% uncertainty in the natural frequency and ±10% in their damping.The resulting bottom-up control-oriented model has 56 states that consist of 12 rigid body states, 18 structural dynamics states, two aerodynamic lag states, and 24 actuator dynamics states.The design model is obtained by trimming and linearizing this model for straight and level flight at 36 equidistant points of the airspeed in the interval [30,65] m/s and for each combination of the perturbed parameters in the structural dynamics.The ν-gap between the high-fidelity and the control-oriented model for different airspeed values and the nominal structural dynamics is shown in Fig. 7.The ν-gap is calculated for the inputs and outputs used for flutter control (see Fig. 5).
The pole migration of the high-fidelity and the control-oriented model is compared in Fig. 8.The full-order model predicts flutter at 52 and 55 m/s at frequencies of 50.2 and 45.8 rad/s, respectively.In comparison, flutter occurs in the control-oriented model at 52 and 56.5 m/s at 50.3 and 46 rad/s, respectively.The flutter speed and frequency accuracy of the control-oriented model is deemed sufficient for control design.

C. Flutter Suppression Control Synthesis
In this section, the flutter suppression control design for the flexible aircraft in Section III-A is discussed.First, two uncertain SISO models are obtained from the reduced-order  control-oriented model detailed in Section III-B.Then, the performance specification used for the optimal control design is described.Finally, the results of the synthesis method in Section II are given.
1) Design Model : Selecting the control surfaces and sensor signals as described in Section III-B (and in Fig. 5) allows us to separate the symmetric and asymmetric flutter modes of the aircraft using the combination of these variables shown in Fig. 9. Two SISO systems are obtained this way: Gsym (s) and Gasym (s).The input and output of Gsym (s) are u f,L + u f,R and q L + q R − 2q, respectively.The states consist of w (vertical velocity), q (pitch rate), the modal coordinates that correspond to the symmetric deformations, the lag states, and the actuator states.For Gasym (s), the input and output are u f,L − u f,R and q L − q R , respectively.The states are v (horizontal velocity), p (roll rate), r (yaw rate), the modal coordinates corresponding to the asymmetric deformations, and lag and actuator states.
Both parametric and dynamic uncertainty are introduced using Gsym (s) and Gasym (s).As the result of the model reduction technique in Section III-B, the state-space matrices of these systems are given on a parameter grid.The grid consists of 36 equidistant points of the TAS between 30 and 65 m/s, five points of the natural frequency in the structural dynamics between ±1% of the nominal value, and five points of the damping in the structural dynamics between ±10% of the nominal value.These two arrays of LTI systems are used to introduce parametric uncertainty in the system.The dependence of the dynamics on the airspeed is expressed via the uncertain parameter δ V .The uncertainty of the natural frequency and the damping in the structural dynamics correspond to uncertain parameters δ ω 0 and δ ξ , respectively.The parameters are normalized uncertainties, i.e., |δ V | ≤ 1, δ ω 0 ≤ 1, and δ ξ ≤ 1. Least-squares fitting is performed to obtain the uncertain state-space matrices of the form , and {C i } 2 i=0 are constant matrices.The elements of A δ , B δ , and C δ are assumed to be a second-order polynomial in δ V based on Chapter 5 in [37].Only A δ depends on δ ω 0 and δ ξ since the perturbation in the damping and natural  frequency influences the position of the poles.Also, this form of A δ , B δ , and C δ provides low error when compared to Gasym and Gsym in the parameter grid points.
Dynamic uncertainty is added to account for the model reduction in Section III-B.As shown in Fig. 10, input multiplicative uncertainty structure is chosen, i.e., both uncertain SISO plants have the form where d (s) is the stable SISO uncertainty block for which The Bode magnitude plot of both of the weights is shown in Fig. 11.These weighting functions serve to add moderate uncertainty at low frequencies and significant uncertainty at high frequencies.The level of uncertainty starts to rise considerably at the flutter frequencies and it is above one (100%) after 100 rad/s.This is in accordance with the ν-gap results of the model reduction in Section III-B (see Fig. 7).
The resulting uncertain systems are written in the LFT form as where (s) is the structured uncertainty block containing the three uncertain parameters δ V , δ ω 0 , and δ ξ and the dynamic uncertainty block d (s).Therefore, for G sym (s) and G asym (s), respectively, the uncertainty sets sym and asym are defined as in ( 5) with   The number of repetitions is especially high for δ V and δ ω 0 .The LFR-toolbox for MATLAB could be used to reduce the size of (s) [38], but since the control synthesis method in Section II samples the uncertain parameters, this is unnecessary.
The poles of G sym (s) and G asym (s) migrate with the change of the uncertain parameters.The domain of migration of the flutter modes is shown in Fig. 12.The introduction of δ ω 0 and δ ξ G sym asym (s) captures the uncertainty and damping ratio in addition to the variation with airspeed.Note that the nominal systems are The gain of the two systems changes due to both types of In Fig. 13, the Bode magnitude plot the nominal systems is depicted along with random samples.In the frequency range below the flutter frequencies, only a moderate variation is observable, while for frequencies above 100 rad/s, the gain of the systems is highly uncertain.Since the nominal systems are stable, the gain at the flutter frequencies are finite for both G sym (s) and G asym (s).For some values of the uncertainty, the gain of these systems is significantly higher or even unbounded.
2) Performance Specification : The purpose of the flutter controller is to robustly stabilize the undamped vibrations of the aircraft.There are further considerations to be considered.The effect of the flutter controller on the rigid body behavior of the aircraft ought to be minimal so that the flutter controller does not interfere with the baseline controller governing the rigid body motion of the aircraft.The controller is to be implemented on an embedded computer whose sampling frequency is 200 Hz.Beyond this constraint, the bandwidth of the flutter controller must not breach the limitation posed by the actuator.Also, the closed loop must remain stable even in the presence of 15-ms output delay introduced by the autopilot hardware.
The generalized plant interconnection in Fig. 14 incorporates all the requirements above for both G sym (s) and G asym (s).The 15-ms delay is represented by τ (s) that is the fourth-order Padé approximation of e −0.015s .The controller is augmented with the filter to ensure appropriate bandwidth.In this way, the sampling constraint is met and the excitation of high-frequency dynamics is avoided.The Bode magnitude plot of F(s) along with the performance constraints is shown in Fig. 15.The task of robust stabilization is expressed as sensitivity minimization.The sensitivity function of both closed loops is S(s) = (1/1 + τ (s)G(s)K (s)F(s)).For any stable SISO loop, the minimal distance between the open-loop Nyquist curve and the −1 point is the inverse of the peak of |S( j ω)|, i.e., [1].It is also known that the gain margin of the loop is at least and the phase margin is at least 2 sin −1 ((1/2||S(s)|| ∞ )) [1].Therefore, we want to achieve ||S(s)|| ∞ ≤ 2 since this provides at least 6-dB gain margin and close to 30 • phase margin.To that end, the weight of the tracking error (i.e., the function) is chosen as W e (s) = (1/2).To limit actuator effort, the weight of the control input is W u (s) = (1/10 • ) = 5.78.
3) Synthesis : The structure of the flutter suppression controller is shown in Fig. 16.The SISO controllers are parameterized as general second-and fourth-order transfer functions, that is where {κ k } 14 k=1 are tunable design parameters.For the synthesis of K asym (s), a frequency grid of 30 points and five basis functions were used.For K sym (s), the number of frequency points and basis functions is 90 and 9, respectively.As shown in Fig. 16, the flutter controller is where Therefore, the final state order of K f (s) is 10.
The two SISO loops are tuned separately.As the result of the synthesis, the tunable components of K f (s) are K asym (s) = 0.0106s2 + 0.5476s + 53.6039 s 2 + 40.4171s + 904.595K sym (s) = 0.0126s 4 + 3.7765s 3 +1479.85s 2 + 48443.9s+ 444523 s 4 + 114.354s 3 +59809.2s2 + 538987s + 1.4737 • 10 6 .It takes ten iterations (i.e., ten samples of p,asym ) to obtain K asym (s) and nine iterations (nine samples of p,sym ) are required for K sym (s). 1 The singular values of the flutter controller are shown in Fig. 17.The controller has enough control authority at the frequencies, and it is sufficiently rolled at the sampling frequency.The performance of K f (s) is analyzed in closed loop with the high-fidelity model linearized at certain speed values in Section IV.

IV. EVALUATION OF THE FLUTTER CONTROLLER USING
THE NONLINEAR HIGH-FIDELITY MODEL This section details the evaluation of the flutter controller described in Section III.First, the baseline control architecture is presented that governs the rigid body motion of the aircraft.Then, the stability, performance, and robustness analysis are described with the baseline and flutter suppression control loops closed.These results are also contrasted with a design based on a single D-scale instead of the D-K iteration.Finally, two time-domain simulations are conducted to confirm the frequency-domain results.All tests are performed using the linearized versions of the nonlinear high-fidelity model in Section III-A.

A. Baseline Control Architecture
The rigid body motion of this aircraft is described by a standard nonlinear six-degree-of-freedom flight mechanics model (e.g., [39]) in terms of translational velocities u, v, and w and angular velocities p (roll), q (pitch), and r (yaw) in the body-fixed frame.Orientation in the earth-fixed reference frame is described in terms of Euler angles (bank), (pitch), and (heading).The angles between the body-fixed frame and the wind axes are angle of attack α and side-slip angle β.The flight path is described with respect to the earth by the path angle γ and the course angle χ.
For the actuation of the rigid body motion, four ruddervators on the V-tail of the aircraft are used, two on the left (u rv,L1 , u rv,L2 ) and two on the right side (u rv,R1 , u rv,R2 ), as shown in Fig. 18.These ruddervators are combining the functionalities of classical rudders and elevators.The symmetric deflections of the ruddervator correspond to classical elevator deflections, whereas asymmetric deflections exhibit rudder deflections.In addition, four control surfaces are available on each wing.As discussed in Section III, the outermost pair (u f,L , u f,R ) is used for flutter control, whereas the innermost pair is used as high lift devices during takeoff and landing.The remaining two   pairs (u a,L2 , u a,R2 , u a,L3 , and u a,R3 ) are utilized in the baseline control law as ailerons to actuate the roll motion of the aircraft.The structure of the baseline controller is a classical cascaded setup shown in Fig. 19.The lateral-directional control problem is necessarily multivariable and requires the coordinated use of aileron command u a and rudder command u r .The innermost loop features roll-attitude ( ) tracking, roll-damping augmentation via the roll rate ( p), and coordinated turn capabilities, i.e., turns without side slip, via feedback of the side-slip angle (β).The outer loop establishes control of the course angle (χ).All controllers are scheduled with velocity to increase the performance over the velocity range.Within the fully automated flight mode, the reference signals for the velocity (V ref ), altitude (H ref ), and course angle (χ ref ) are provided by a dedicated navigation algorithm.It uses the GPS longitudinal and lateral position of the aircraft (x a and y a ) as well as the current course angle (χ) to provide the commands.More details on the algorithm can be found in [40].
Structurewise, the control loops use scheduled elements of proportional-integral-derivative (PID) controllers with additional roll-off filters in the inner loops to ensure that no aeroelastic mode is excited by the baseline controller.
A scheduling in dependence of the TAS V TAS is used to ensure an adequate performance over the velocity range from 32 to 70 m/s.For the scheduling, a first-or second-order polynomial in V TAS is applied.The free parameters are directly included in a structured controller optimization problem.A comprehensive summary of the used controller structures for each cascaded loop is provided in Table I, including the channel description in the controller architecture and the implemented scheduling.
Note that these controller outputs u e , u a , and u r defer from the actual surface inputs to ease the control design task.Thus, they need to be transformed to physical actuator commands via an adequate control allocation.The commands to the actuators of the two aileron pairs are determined by u a,L2 = u a,L3 = 0.5u a u a,R2 = u a,R3 = −0.5u a to generate the required differential aileron deflections for roll motion control.For the ruddervators, superposition of the elevator command u e and the rudder command u r is applied by u rv,L1 = u rv,L2 = u e + 0.5u r u rv,R1 = u rv,R2 = u e − 0.5u r .
Thus, symmetric deflections on the left and right of the ruddervators correspond to elevator commands, whereas differential deflections establish rudder commands.The free parameters of the control laws are tuned to satisfy certain performance and robustness criteria as explained in [40].

B. Closed-Loop Stability and Performance Evaluation
The pole trajectories of the closed loop (with both the flutter and baseline controller) are shown in Fig. 20   whereas the symmetric flutter mode crosses over to the right half-plane at 68 m/s.Since the closed loop is unstable at 68 m/s, this is called the absolute flutter speed.Several other modes are influenced by the two controllers but none becomes unstable, and therefore, no adverse interaction between the flutter and baseline controller is observed.
Let us verify the result of the optimization in Section III in terms of the sensitivity functions.To this end, only the flutter controller is applied and the closed loop is opened loop-at-atime to get the sensitivity function of the two SISO loops separately.The magnitudes from 30 to 60 m/s are shown in Fig. 21.The regions where the sensitivity functions are greater than one are concentrated to low frequency due to the bandwidth of the flutter controller.The sensitivity functions are also reasonably flat.The peaks are 2.41 and 2.08 for the symmetric and asymmetric loop, respectively.This means that the objective of pushing the sensitivity functions below two was not completely accomplished.The resulting stability margins are computed in Section IV-C.
To demonstrate the increase in damping, we compare the gain of the system with and without flutter control.Specifically, Fig. 22 shows the gain from an additive disturbance on u f,L to the vertical acceleration measured by the IMU on the left wing (see Fig. 5).The figure shows a significant decrease in damping around the flutter frequencies.Comparing Fig. 22 to Fig. 21, we observe that the gain increased in the low-frequency range, where the sensitivity values are high, and remained close to the original in the high-frequency range where the sensitivity values are close to one.

C. Robustness Analysis
The robustness loop is analyzed using define margins, interconnection in Fig. 23.Here, k ∈ are complex scalar of α = (1 δ k δ k where δ k C, k = 1, . . ., 5. We differentiate between loop-at-a-time and multiloop disk margins.The loop-at-a-time disk margin is the largest simultaneous gain and phase variation in a single channel for the closed loop remains stable.It is the largest value of m l that the closed loop in Fig. 23 is well posed and stable for |δ l | < m l , while δ k = 0 for k = l.For multiloop margins, the in several loops are to vary simultaneously and independently.The classical gain and phase margin interpretation of the loop-at-a-time margins for each input and output channel in Fig. 23 is shown in Fig. 24.For four channels, the margins the same.the q channel is by the two control loops, it shows greater robustness.The loop-ata-time margins have acceptable values up to 60 m/s for all channels, the gain and margins are 5.5 dB and 34 • , respectively.The margins also degrade slowly after 60 m/s; therefore, there is a sufficiently wide safety region around 60 m/s.In view of these facts, the robust flutter speed (RFS), i.e., the speed at which the demonstrator may fly safely, is determined to be 60 m/s.The loop-at-a-time margins became zero at 68 m/s, which is the absolute flutter speed already established in Section IV-B.The multiloop margins of the closed loop are shown in Fig. 25.These margins also start to significantly degrade after 60 m/s.The most notable is the simultaneous input-output margin for which all five uncertainties are allowed to vary at the same time.At the RFS, the simultaneous input-output margin in 2.2 dB and 14 • .The rate of degradation for these margins is similar to the loop-at-a-time margins.This confirms the RFS to be 60 m/s.
The critical frequency corresponding to a loop margin is the frequency where the Nyquist curve of the open loop touches the critical point (−1 or 0).The critical frequencies for all margins and all channels are less than 65 rad/s.This ensures that these robustness measures are meaningful because this shows that the perturbation causes instability in the low-frequency range where the accuracy of the modeling is acceptable.At higher frequencies, any model is more uncertain, and therefore, larger gain and phase margin are required.
We designed a second flutter controller using a variation on Algorithm 1 for comparison.In this version, instead of the D-K iteration in Steps 4-7, a single parameterized D-scale is tuned together with the controller, similar to the hybrid relaxation method in [14].The D-scales are defined to be state-space systems with state order equal to the number of basis functions in Section III-C3 (i.e., 5 for the synthesis of K asym (s) and 9 for K sym (s)).The "companion" parameterization option of the tunableSS function is chosen for both.The performance setup and controller structure is identical to what is described in Section III-C2.The resulting controller guarantees 2-dB simultaneous input-output gain margin and 13 • phase margin up to 60 m/s.Comparing these values to 2.2 dB and 14 • obtained when using Algorithm 1 reveals the advantage of having separate D-scales for the samples in the synthesis.a controller designed using by values the uncertain parameters our design model can damping the flutter modes at speeds.We unable to beyond the open-loop speed (OLFS) using a controller designed for a single LTI system.

D. Time Domain Simulations
Two time-domain simulations are conducted to verify the frequency-domain results in Sections IV-B and IV-C.In both cases, the baseline controller described in Section IV-A and the flutter controller from Section III-C3 are operating at the same time.The setup for both simulations is shown in Fig. 26.The controllers are implemented in discrete time (with 200-Hz sampling frequency), and a 15-s output delay is included to simulate the delay caused by the sensors.Gaussian measurement noise is added to the plant output.The only reference command input of the baseline controller used is the speed reference to demonstrate the behavior of the closed loop at different speed values.
In the first simulation, the speed is increased starting at 46 m/s and following a staircase reference.No wind disturbance is used in this simulation.The vertical acceleration from the two IMUs on the wing in Fig. 5 is recorded to demonstrate the load caused by the oscillations.These are denoted by a z,L and a z,R for the left and right wings, respectively.The results are shown in Fig. 27.In the graph on the top, the airspeed is depicted as a function of time, while the graph on the bottom shows the vertical acceleration.If the flutter suppression controller is not activated and the aircraft passes the OLFS, the acceleration of the wingtip becomes greater than g indicating structural failure.With the flutter controller however, the demonstrator safely passes the OLFS.The wingtip  acceleration values are below critical at the RFS as well.The oscillations become damaging at 63 m/s, i.e., 3 m/s above the RFS.This is achieved with u f,L < 0.5 • and u f,R < 0.5 • .
In the second simulation, we test whether it is indeed possible to fly at the RFS even in the presence of continuous Dryden wind gusts.The aircraft is accelerated from 43 to 60 m/s, and then, it follows a constant 60-m/s speed reference command.This is illustrated in the top diagram of Fig. 28.Without flutter control, the wing tip accelerations become critical a few seconds after the OLFS was passed.When the flutter controller is used, the RFS is reached and acceptable acceleration values are maintained even if the wind gusts push speed beyond the demonstrates the robustness of the flutter controller discussed in Section IV-C.The control input is shown in Fig. 29.The control surface deflections increase when the open loop and the RFS are passed.For the entire simulation, the control input is within ±1.5 • bounds.

V. CONCLUSION
A control design algorithm is to minimize the worst case of the closed loop by tuning a fixed structure controller in the presence of mixed uncertainty.The method is iterative.In the analysis step, the worst case sample of the uncertain parameters is computed.In the synthesis step, D-scales are constructed by solving a convex optimization problem for each sample separately to account for the dynamic uncertainty.The controller is tuned for the collection of scaled samples.The applicability of this method for the flutter suppression control design of the flexible demonstrator aircraft of the FLEXOP project is also demonstrated.The problem is articulated as the minimization of the sensitivity function and control effort for two SISO loops.Disk margins are computed using the baseline controller and the high-fidelity model of the aircraft.Based on these margins and on time-domain simulations, the OLFS of 52 m/s is extended to 60 m/s, which is a 15% increase of the flight envelope.The controller is low order and LTI and therefore easy to implement on the onboard computer of the aircraft.
Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.
Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.

Fig. 1 .
Fig. 1.Closed loop interconnection for the control design.
Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.

Fig. 7 .Fig. 8 .
Fig. 7. ν-gap values as a function of frequency between the high-fidelity and the control-oriented model for the inputs and outputs used for the flutter suppression control design.
Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.

Fig. 9 .
Fig. 9. Structure of the control loop using two SISO controllers to stabilize the symmetric and asymmetric flutter modes separately.

Fig. 10 .
Fig. 10.Uncertain plant interconnection with the dynamic uncertainty d (s) and A δ , B δ , and C δ matrices depending on the parametric uncertainty.

Fig. 12 .
Fig. 12. Domain of the migration of the flutter mode poles due to the parametric uncertainty.

Fig. 13 .
Fig. 13.Bode magnitude plot of the nominal value and random samples of the uncertain systems G sym (s) and G asym (s).

Fig. 15 .
Fig.15.Bode magnitude plot of F(s) and the performance constraints.The actuator dynamics G act (s) are given in(16).

Fig. 17 .
Fig. 17.Singular values of the flutter controller along with the frequency-domain constraints.

Fig. 18 .
Fig. 18.Control surface configuration of the baseline control architecture.The control inputs are marked at the corresponding control surface.
Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.

Fig. 19 .
Fig. 19.Control architecture for fully automated flight and augmented flight (indicated in orange).

Fig. 20 .
Fig. 20.Change of pole trajectories due to the flutter and baseline controller.
. The damping of both flutter modes is increased significantly.The asymmetric flutter mode is stabilized up to 70 m/s, Authorized licensed use limited to: Institute for Computer Science and Control.Downloaded on January 10,2022 at 10:26:25 UTC from IEEE Xplore.Restrictions apply.

Fig. 22 .
Fig. 22. Gain from input disturbance to wingtip acceleration in open and closed loops.

TABLE I SUMMARY
OF THE CONTROL LOOPS OF THE FLEXOP BASELINE FLIGHT CONTROL SYSTEM WITH THE INNER LOOP FUNCTIONS (FIRST PART) AND AUTOPILOT FUNCTIONS (SECOND PART)