Hugo Fournier
email: hugo.h.fournier@airbus.com

Paolo Massioni

Minh Tu Pham

Laurent Bako

Robin Vernay

Robust Gust Load Alleviation of Flexible Aircraft equipped with Lidar

This work concerns the problem of gust load alleviation of a flexible aircraft by focusing on Airbus's XRF1 aircraft concept, with a fully actuated wing. The aircraft is equipped with a lidar sensor, which measures the wind velocity ahead of it, together with standard sensors.

Based on the available measurements, controllers are then designed by H ∞ and µ syntheses, with emphasis put on the multiple-input multiple-output robustness in order to reduce a selected set of loads due to the wind, hence potentially saving mass in the aircraft design. A statespace model of the flexible aircraft is obtained by means of an aeroelastic computation and system identification from frequency data. The controllers' performance is evaluated through their capability to reduce shear force, bending and torsion moments on different locations of the aircraft in response to different types of discrete gusts and continuous turbulence; the constraints of the sensors and the actuators are taken into account. The gain in performance due to the use of lidar is assessed, and a trade-off is discussed concerning the optimal measurement distance. Finally, the closed-loop robustness is assessed by simulations where different types of uncertainties are applied to the system. V ertical and lateral winds are known to cause loads in an aircraft by putting strain on the structures that interact the most with air flows, such as the wings and the tailplane. By creating high instantaneous loads and fatigue in the aircraft structure through its life cycle, gusts and turbulence require the addition of structural reinforcements.

Alleviating these loads can save a significant amount of mass during the design phase of the aircraft. This can be achieved through passive load alleviation by an adequate design of the aircraft geometry, and through active load alleviation. The latter has been successfully implemented since the 1970s with the Lockheed L-1011 with a mass-saving objective, and more recently on the Airbus A380 and the Boeing 787 [START_REF] Regan | Survey of Applications of Active Control Technology for Gust Alleviation and New Challenges for Lighter-weight Aircraft[END_REF]. By deflecting the ailerons symmetrically upwards when a vertical gust is detected, the lift is locally reduced and concentrated in the innermost part of the wings, which reduces the integrated loads at the wing root. Another strategy for active gust load alleviation (GLA) would consist in turning the plane towards the wind direction to reduce its influence, by use of the elevators.

Active gust load alleviation will most likely gain importance in the future generations of aircraft. The autonomy of the aircraft electric flight control system (EFCS) at the expense of the pilot's freedom to maneuver in flight will allow more efficient load alleviation control laws by reducing the risk of a detrimental pilot input and by coordinating the flight control and gust alleviation systems. Furthermore, many new aircraft concepts choose to increase the wing aspect ratio and to reduce its structural weight to enhance aircraft energy efficiency. This increases airframe flexibility, leading to aeroelastic modes frequencies that approach the flight dynamics frequencies, creating flight dynamics-aeroelastic interactions. The recent development of lidar sensors [START_REF] Herbst | Design of a monolithic Michelson interferometer for fringe imaging in a near-field, UV, direct-detection Doppler wind lidar[END_REF][START_REF] Rabadan | Airborne Lidar for Automatic Feedforward Control of Turbulent In-Flight Phenomena[END_REF] that can measure the vertical and lateral wind several hundreds of meters ahead of the aircraft allows new perspectives for gust load alleviation, by anticipating the response of the aircraft to the incoming wind perturbation. With the availability of computational capabilities, this motivates the development of more complex gust alleviation techniques, based on a full automation of the aircraft and taking advantage of its complete aeroelastic dynamics.

The first step in the development of gust load alleviation control laws is the aircraft modeling. This is of particular importance for flexible aircraft, where aeroelastic and flight dynamics modes interact due to low frequency separation [START_REF] Waszak | Flight dynamics of aeroelastic vehicles[END_REF]. Historically, the first models considered the six-degree-of-freedom equations of motions for the rigid aircraft, corrected to include quasi-static aeroelastic effects [START_REF] Bisplinghoff | Aeroelasticity[END_REF][START_REF] Etkin | Dynamics of Flight -Stability and Control[END_REF]. With the availability of computers, models that include n additional degrees of freedom related to the structural dynamics were developed in order to fully capture the interaction between aeroelastic and flight dynamics. They generally rely on Lagranges's equations which can be simplified by an accurate choice of reference axis. Mean axes, with respect to which the linear and angular momenta associated with structural vibrations are zero, have been particularly used in the literature [START_REF] Waszak | Flight dynamics of aeroelastic vehicles[END_REF][START_REF] Schmidt | Modeling and Simulation of Flexible Flight Vehicles[END_REF][START_REF] Reschke | Flight loads analysis with inertially coupled equations of motion[END_REF]. In [START_REF] Neto | Formulation of the Flight Dynamics of Flexible Aircraft Using General Body Axes[END_REF] a formulation with an arbitrary choice of body axes has been proposed. In their general formulation, Lagrange's equations lead to nonlinear models that include inertial coupling and can be used for simulations of a large class of aircraft including very flexible aircraft, in which large aeroelastic deflections can exist [START_REF] Patil | Flight Dynamics of Highly Flexible Flying Wings[END_REF][START_REF] Su | Nonlinear Aeroelasticity of a Very Flexible Blended-Wing-Body Aircraft[END_REF]. On the contrary, when the objective is to design a controller, a linear model is preferred in which the inertial effects are neglected. This is the case of the present work, where a moderately flexible aircraft is studied and small deflections are considered. Linear aeroelastic models making use of finite element models (FEM) have been developed [START_REF] Rodden | Equations of motion of a quasisteady flight vehicle utilizing restrained static aeroelastic characteristics[END_REF][START_REF] Dykman | Structural Dynamics and Quasistatic Aeroelastic Equations of Motion[END_REF][START_REF] Saltari | Finite-Element-Based Modeling for Flight Dynamics and Aeroelasticity of Flexible Aircraft[END_REF]. They can be coupled with linear unsteady aerodynamic methods such as the doublet lattice method (DLM) [START_REF] Albano | A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows[END_REF][START_REF] Kalman | Application of the Doublet-Lattice Method to Nonplanar Configurations in Subsonic Flow[END_REF] which compute the generalized aerodynamic forces (GAF) at different frequencies. A linear aeroelastic model can then be obtained by rational function approximation (RFA) of these forces by use of Roger [START_REF] Roger | Airplane Math Modeling Methods for Active Control Design[END_REF] or Karpel [START_REF] Karpel | Design for Active Flutter Suppression and Gust Alleviation Using State-Space Aeroelastic Modeling[END_REF] methods for instance, and improved by use of reliable external aeroelastic data [START_REF] Baldelli | Unified Aeroelastic and Flight Dynamic Formulation via Rational Function Approximations[END_REF]. However, such techniques can lead to high-order models when many flexible modes are included.

In [START_REF] Hesse | Reduced-Order Aeroelastic Models for Dynamics of Maneuvering Flexible Aircraft[END_REF] a reduction is performed by filtering out unnecessary high-frequency modes, and using balanced realizations.

In the present work, the frequency response of the aeroelastic system is computed by means of a full FEM model of the aircraft and an unsteady aerodynamic model obtained by DLM, with aerodynamic corrections coming from computational fluid dynamics (CFD) to account for DLM's inaccuracy at low frequencies. A state-space model can then be obtained by system identification based on these frequency data. The Loewner interpolation method [START_REF] Lefteriu | A New Approach to Modeling Multiport Systems From Frequency-Domain Data[END_REF][START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF] is used in this work for its simplicity and its capability to identify multiple-input multiple-output (MIMO) state-space models of relatively low order with good accuracy.

The most straightforward method for active gust load alleviation of flexible aircraft is feedback control [START_REF] Fonte | Gust Load Alleviation for a Regional Aircraft Through a Static Output Feedback[END_REF][START_REF] Che | Integrated Flight/Structural Mode Control for Very Flexible Aircraft Using L1 Adaptive Output Feedback Controller[END_REF][START_REF] Stanford | Gradient-Based Aeroservoelastic Optimization with Static Output Feedback[END_REF][START_REF] Mokrani | Multiple-Input Multiple-Output Experimental Aeroelastic Control Using a Receptance-Based Method[END_REF] also broadly used in flight control [START_REF] Luspay | Flight control design for a highly flexible flutter demonstrator[END_REF][START_REF] Ustinov | Control Law Design for Flexible Aircraft: Comparison of the H-Infinity-based and Classical Methods[END_REF]. It is based on the aircraft measured output, typically the pitch rate, yaw rate and vertical load factor obtained from gyroscopes and accelerometers. Feedback control allows a direct design of the closed loop poles, ensuring stability and robustness to the system uncertainties, without depending on the exact knowledge of the aircraft system. In order to improve the loads reduction, a feedforward controller that directly uses the measurement of the wind perturbation can be used in addition to feedback control [START_REF] Alam | Active gust load alleviation system for flexible aircraft: Mixed feedforward/feedback approach[END_REF][START_REF] Wildschek | Gust Load Alleviation on a Large Blended Wing Body Airliner[END_REF][START_REF] Magar | Active Gust Alleviation using Artificial Hair Sensors and Feedforward Control[END_REF] or it can be used alone [START_REF] Zhao | Gust Load Alleviation on a Large Transport Airplane[END_REF][START_REF] Li | Flight Testing an Adaptive Feedforward Controller for Gust Loads Alleviation on a Flexible Aircraft[END_REF][START_REF] Wildschek | Active Wing Load Alleviation with an Adaptive Feed-forward Control Algorithm[END_REF]. A perfect feedforward controller can theoretically be obtained by direct inversion of the system's dynamics, but this is generally impossible due to actuator limitations, and the resulting controller would be very sensitive to system uncertainties. In [START_REF] Caverly | Gust-Load Alleviation of a Flexible Aircraft using a Disturbance Observer[END_REF] such a technique is used in addition to a linear controller and a disturbance observer for the gust load alleviation of a flexible aircraft. In [START_REF] Wang | Flexible Aircraft Gust Load Alleviation with Incremental Nonlinear Dynamic Inversion[END_REF], a method based on dynamics inversion is used with state estimation while ensuring robustness to system's uncertainties. Both of these methods work with a nonlinear model of the aircraft. A promising alternative is the employment of adaptive control, in which the system's dynamics are constantly estimated and the control is adapted accordingly. It generally relies on a finite-impulse response (FIR) model of the controller [START_REF] Zhao | Gust Load Alleviation on a Large Transport Airplane[END_REF][START_REF] Li | Flight Testing an Adaptive Feedforward Controller for Gust Loads Alleviation on a Flexible Aircraft[END_REF][START_REF] Wildschek | Active Wing Load Alleviation with an Adaptive Feed-forward Control Algorithm[END_REF][START_REF] Zeng | Adaptive Feedforward Control for Gust Load Alleviation[END_REF][START_REF] Wildschek | Wind Tunnel Testing of an Adaptive Control System for Vibration Suppression on Aircraft[END_REF][START_REF] Wang | Adaptive Feedforward Control Design for Gust Loads Alleviation of Highly Flexible Aircraft[END_REF].

Developments in the last decades [START_REF] Herbst | Design of a monolithic Michelson interferometer for fringe imaging in a near-field, UV, direct-detection Doppler wind lidar[END_REF][START_REF] Rabadan | Airborne Lidar for Automatic Feedforward Control of Turbulent In-Flight Phenomena[END_REF] such as Airbus AWIATOR program [START_REF] Schmitt | The AWIATOR airborne LIDAR turbulence sensor[END_REF] or DLR's DELICAT program [START_REF] Vrancken | Airborne forward-pointing UV Rayleigh lidar for remote clear air turbulence detection: system design and performance[END_REF],

thanks to which airborne lidar sensors emit a LASER and measure the light scattered by air molecules, create new opportunities in gust load alleviation by not only diminishing the effects of wind but also preventing them. While simple controller structures are commonly used for feedback control, finding a controller structure that takes the wind predictions as input is not straightforward. Model predictive control (MPC) methods that perform an online optimization can be used, they have been applied to gust load alleviation [START_REF] Haghighat | Model-Predictive Gust Load Alleviation Controller for a Highly Flexible Aircraft[END_REF][START_REF] Giesseler | Model Predictive Control for Gust Load Alleviation[END_REF] and maneuver load alleviation [START_REF] De Freitas | Model Predictive Control Architectures for Maneuver Load Alleviation in Very Flexible Aircraft[END_REF][START_REF] Gaulocher | Aircraft Load Alleviation During Maneuvers Using Optimal Control Surface Combinations[END_REF]. More complex strategies based on wavelet decomposition of the lidar signal have been used in order to map the signal onto the different actuators depending on the frequency and the amplitude of the gust [START_REF] Fezans | Combined Feedback and LIDAR-Based Feedforward Active Load Alleviation[END_REF][START_REF] Fezans | In-flight remote sensing and identification of gusts, turbulence, and wake vortices using a Doppler LIDAR[END_REF]. The H ∞ and µ syntheses [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Skogestad | Multivariable feedback control : analysis and design[END_REF][START_REF] Vidyasagar | Robust controllers for uncertain linear multivariable systems[END_REF] provide a frame in which one can perform a multi-objective optimization based on multiple sensors measurements including the lidar, while ensuring compliance with robustness constraints and system's limitations. Unlike MPC, closed-loop stability and robustness can be imposed with these methods, and the on-line computational cost is much lower. Already used for robust GLA with classical sensors [START_REF] Zeng | Experimental Model-Based Aeroelastic Control for Flutter Suppression and Gust-Load Alleviation[END_REF][START_REF] Cook | Robust Gust Alleviation and Stabilization of Very Flexible Aircraft[END_REF][START_REF] Ferrier | Active Gust Load Alleviation of High-Aspect Ratio Flexible Wing Aircraft[END_REF] and for active flutter suppression of flexible aircraft [START_REF] Waitman | H-infinity Control Design for Active Flutter Suppression of Flexible-Wing Unmanned Aerial Vehicle Demonstrator[END_REF][START_REF] Pusch | Structured Control Design for a Highly Flexible Flutter Demonstrator[END_REF], the application of H ∞ to lidar-based GLA has been proposed in [START_REF] Khalil | Performance Enhancement of Gust Load Alleviation Systems for Flexible Aircraft using H∞ Optimal Control with Preview[END_REF] where the additional loads reduction obtained with a noiseless lidar was assessed on a simple sailplane model including two elastic modes, without robustness constraints.

This work can be considered as fitting into the same framework of what it described above, aiming at developing robust control synthesis techniques adapted to the GLA of a general class of flexible aircraft, with the following contributions:

• The loads reduction due to the lidar is evaluated in presence of measurement noise. This is achieved with realistic assumptions on the system constraints with guaranteed design margins. The optimal measurement distance of the lidar is investigated for two lidar configurations, which differ according to the number of measurements at intermediate distances they can perform.

• Robustness of the obtained controllers to realistic system uncertainties is assessed through different uncertainty cases. The multiple-input multiple-output robustness in both stability and performance is obtained thanks to the introduction of so-called disk uncertainties [START_REF] Blight | Practical control law design for aircraft using multivariable techniques[END_REF][START_REF] Seiler | An Introduction to Disk Margins [Lecture Notes[END_REF] in the controller design, which apply to several inputs and outputs simultaneously and can be defined in a way compatible with the H ∞ framework.

The GLA problem and the strategy employed to solve it are first formally stated in section II along with the system description. In section III the procedure that leads to the reduced order state-space model of the flexible aircraft is described. In section IV the H ∞ and µ syntheses are presented, and the detail of how to obtain a robust and efficient controller for the GLA of a flexible aircraft equipped with lidar is explained. The results obtained by application of these techniques are presented in section V.

II. Gust Load Alleviation of a Flexible Aircraft: Problem Statement and Strategy

This section presents the characteristics of the aircraft relevant to the design of a gust load alleviation control, as well as the lidar system with two plausible variants for future use in industry. The models used for the wind description, mainly based on official certification documents [START_REF]Certification Specifications for Large Aeroplanes, Amendment 3[END_REF], are then described before stating the problem studied in this work, which is the design of a robust control for alleviating the loads caused by gusts and turbulence in a flexible aircraft. The approach presented in this work is kept as general as possible and can be applied to the GLA of any flexible aircraft model. Table 1 Properties of the XRF1 aircraft concept

A. Aircraft Description

The models used in this work come from the XRF1 aircraft concept, a platform developed by Airbus for research purposes. A view of its finite element model is shown in Fig. 1. Some elements of description of the aircraft and its mass configuration are shown in Table 1. The position of the center of gravity (c.g.) of the aircraft is given as a percentage of the mean aerodynamic chord, backwards from the reference leading edge of the wing. It corresponds to one mass configuration for a given fuel tanks filling and a given payload; other configurations will be studied in the robustness analysis of section V.D.

The aircraft model possesses several control surfaces that can be used for gust load alleviation: one pair of elevator and one rudder on the tailplane, and twelve pairs of ailerons on the wing trailing edge. This last feature is a particularity of the model considered here. With many control surfaces employed to create local forces and moments onto the wing structure, a higher load reduction is expected. In addition to the predicted wind velocity measured by the lidar sensor, four measured variables are considered in this work: the pitch angle and rate, the vertical load factor and the angle of and alpha-probe, whose noise is low enough not to have any impact on the response.

B. Lidar Sensor Description

The aircraft studied in this work is equipped with a lidar sensor that can remotely measure the wind velocity ahead of the aircraft. This is achieved [START_REF] Herbst | Design of a monolithic Michelson interferometer for fringe imaging in a near-field, UV, direct-detection Doppler wind lidar[END_REF] by a LASER emitting light forward from the nose of the aircraft. The light is then scattered by air molecules through the quasi-static molecular scattering Rayleigh-Brillouin effect, a sensor at the nose of the aircraft measures by Doppler effect the frequency shift of the reflected light and the wind velocity can be recovered. By emitting light in different directions (four in the considered case) and measuring the projection of the relative velocity of the aircraft with respect to the air onto the LASER beams, the full velocity vector is obtained as seen in Fig. 2.

Two different types of lidar will be studied in this work. The first kind measures the wind velocity at one given distance ahead of the aircraft; by keeping in buffer memory the velocity, the information of the wind profile the aircraft is going to encounter can be used in the GLA control design. In the second type of lidar, intermediate measurements are performed at different distances of the aircraft, directly providing the wind profile. Because the lidar precision decreases with the distance [START_REF] Herbst | Design of a monolithic Michelson interferometer for fringe imaging in a near-field, UV, direct-detection Doppler wind lidar[END_REF], the closest lidar measurements are thought to provide more accurate information. The model of the lidar measurement noise will be studied in section IV.B. In Fig. 2, the first lidar configuration corresponds to the case in which only the measurement located at the last dot is performed, while the intermediate measurements located at the blue dots are also available in the second configuration.

C. Gust and Turbulence Model

High velocity wind can cause important loads in the aircraft structure, especially on the wing and tailplane, but also on the fuselage. For this reason, simulations and tests must be performed along the aircraft design phases in order to verify that maximum loads that can occur due to high winds during the aircraft's life cycle do not overcome the maximum design loads withstandable by the structure. The design of several parts of the aircraft is driven by the worst cases of winds. but also has a dynamic influence resulting from the aeroelastic interaction between the structure and the unsteady aerodynamics.

A conservative yet flexible design can be obtained by choosing a not large number of dynamical cases of wind profiles, representative of the different worst cases that can occur at different locations in the aircraft. The most common approach is to distinguish between:

• discrete gusts, short events of high wind velocity lasting at maximum a few tenths of seconds

• continuous turbulence, longer events of random wind velocity with less abrupt changes.

These two models correspond to the main cases described in official certification documents, such as the European CS-25 [START_REF]Certification Specifications for Large Aeroplanes, Amendment 3[END_REF] that states the requirements for large aeroplanes. Discrete gusts are defined in such a way as to target one frequency of the aeroelastic system, while the continuous turbulence excites a range of frequencies. In particular, the article 25.341 provides a mathematical definition of discrete gusts and continuous turbulence that must be used for certification. The vertical discrete gust v z is described in the time domain by the 1cos expression:

v z (t) =            U gu st 2 1 -cos πtV L gu s t if 0 ≤ t ≤ 2L gu s t V 0 otherwise ( 1
)
where t is the time in s, V the true aircraft speed in ft/s. L gust is called the gradient distance in ft and defines one half of the gust duration. U gu st is the gust amplitude in ft/s, which depends on L gust and other flight parameters. The range of gradient distances is generally comprised between 30 ft and 350 ft. In Fig. 3 time simulations of three discrete gusts of various gradient distances are shown for a true aircraft speed of 840 ft/s. Short gusts have a lower amplitude but can excite higher frequency modes: due to the aircraft inertia and actuators limitations, they are more difficult to reject. The continuous turbulence is a stochastic process, generally defined by its power spectral density (PSD). In this work, the empirical von Kármán PSD for vertical turbulence [START_REF] Diederich | Effect of Spanwise Variations in Gust Intensity on the Lift Due to Atmospheric Turbulence[END_REF] is used:

Φ(ω) = U 2 t ur b L t ur b πV 1 + 8 3 (1.339 ω L t ur b V ) 2 [1 + (1.339 ω L t ur b V ) 2 ] 11/6 ( 2
)
where ω is the angular frequency in rad/s, L t ur b is the scale of turbulence in ft and U t ur b is the turbulence standard deviation in ft/s. A von Kármán PSD corresponding to a 2500 ft scale length and a 85 ft/s standard deviation is represented in Fig. 4.

D. Problem Statement and Control Strategy

With the elements of description of the aircraft system and of the wind perturbations exposed above, the problem can be stated as follows. Control laws minimizing the loads caused by the different cases of gusts and turbulence shall be defined. Maximum values of the control surfaces deflection angles and deflection rates must be defined based on the actuators capabilities. Closed-loop stability must be ensured, and given robustness margins must be satisfied. The models used for the design and in the simulations shall account for the sensors and actuators limitations, and represent accurately the aircraft response to the wind input. The loads consist of the shear force, bending moment and torsional moment at different stations along the wing span. They are not directly measured.

To do so, state-space controllers are developed in this work, by use of the H ∞ framework described in section IV.

With more degrees of freedom than classical controllers based on gains and simple filters, the state-space controllers are thought to better exploit the full frequency content of the complex aeroelastic system. They also allow the natural inclusion of the multiple measured output (sensors and different values coming from the lidar measurements) and multiple control input in a unified formulation. The choice of the loads to minimize is of particular importance: indeed, trying to minimize all loads simultaneously would lead to poor overall performance since the control authority of each actuator is not the same with respect to all loads: for example, the loads on the wing are better controlled by the ailerons than the loads on the tail-plane. Furthermore, reducing one can sometimes increase another. For example, an upward aileron deflection locally decreases the lift, hence reducing the shear force and bending moment, but at the same time decreases the wing camber, creating a pitching moment that increases the wing torsion. In this work, priority is put on the alleviation of the bending moment at the wing root, as it could avoid strong structural reinforcements hence saving a significant amount of mass. The reduction of the shear force, bending and torsional moments on the horizontal tail plane (HTP) is also included in the control objectives.

III. Flexible Aircraft Model

The design of a controller requires modeling the flexible aircraft as a dynamical system. First an aeroelastic computation is performed, which accurately models the frequency response of the aircraft to surface controls and wind input thanks to FEM and DLM models. Aerodynamic corrections based on reliable CFD data are applied to account for the DLM's inaccuracy at low frequencies. Such response is obtained at different frequency values and a continuous state-space model is obtained from the frequency response for the purpose of being used by H ∞ and µ syntheses.

A. Frequency-Domain Modeling

In GLA the full dynamical interaction between the aircraft flexible structure and the unsteady aerodynamics must be accounted for. The structure of the aircraft is modeled by FEM. The user chooses the number of elastic modes they want to include, in addition to the flight dynamics modes. The aerodynamics are modeled by DLM [START_REF] Albano | A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows[END_REF][START_REF] Kalman | Application of the Doublet-Lattice Method to Nonplanar Configurations in Subsonic Flow[END_REF]. Aerodynamic corrections based on reliable CFD data are applied to account for the DLM's inaccuracy at low frequencies. From a meshing of the aircraft external surfaces, different from the FEM meshing, aerodynamic influence coefficient (AIC) matrices of external control surfaces are computed based on a linearized unsteady aerodynamic potential theory. From the structural, inertial and aerodynamic models, an aeroelastic solver such as sol145 of Nastran [START_REF][END_REF] computes the different elements of the equations of motion expressed in the modal base:

M hh s 2 + B hh s + K hh h = q Q hh (M, k)h + Q vh (M, k) v V (3)
where h are the generalized modal coordinates, including the flight dynamics and elastic degrees of freedom and the control surface deflection angles. s = jω is the Laplace operator, q = 1 2 ρV 2 is the dynamic pressure with ρ the air density and V the true aircraft speed, M hh is the modal mass matrix, B hh is the modal damping matrix, K hh the modal stiffness matrix, Q hh and Q vh are the modal and turbulent forces, v is the wind velocity. The computation is performed at different Mach numbers M and reduced frequencies k = ωb V with b the reference semi-chord and ω the angular frequency. The outputs of interest y are obtained directly from the modal displacements by:

y = H y (s)h (4)
where H y (s) is a known transfer function matrix. The output variable y contains the loads and the sensors measurements:

y =            y l oads y sens            (5)
Note that because of the dependence of the generalized forces Q hh and Q vh on the angular frequency, the frequency response is only computed at a finite number of values of ω. By taking the Fourier transform of Eq. ( 3) one obtains transfer functions from surface controls deflection angle input u included in h and wind input v to the output of interest y, at different frequencies. The frequency data response (FRD) are noted H i ∈ C N y × N u associated with the frequencies (ω i ) i=1... N f at which the response is computed, where N y is the number of output variables, including measured variables and loads, N u is the number of inputs, including the different control surfaces deflection angles and the wind velocity. Note that unlike techniques based on time data from experiments [START_REF] Halder | Signal Interpolation Augmented Linear Nonintrusive Reduced-Order Model for Aeroelastic Applications[END_REF], the frequency data obtained with the aeroelastic computation are noise-free.

B. Time-domain Modeling: System Identification from Frequency Response Data

Given the FRD (ω i , H i ) i=1... N f , a state-space model of the following form is sought:

           ẋ = Ax + B u u + B v v y = C x + D u u + D v v (6)
where x is the state vector, y is the measured output, u is the control input, v is the wind velocity, A, B u , B v , C, D u and D v are matrices to be identified. Note that y is the same as in Eq. ( 4). On the contrary, the state vector x in the state-space model does not correspond to the modal coordinates h of Eq. ( 3), and the input vector u was included in h.

With twelve pairs of ailerons in addition to the elevators and wind inputs, and the loads output at different locations of the aircraft structure, the model to identify is high-dimensional. In order to account for the different aeroelastic modes, the identified model must also be high-order. These particularities make the identification process numerically intensive, and generally requiring high amounts of computer memory. Furthermore, in order to be used in H ∞ and µ syntheses, the identified model must be of low order, hence the need of dedicated system identification techniques.

One method consists in taking Eq. ( 3) and performing a RFA of Q hh and Q vh , for example with Roger [START_REF] Roger | Airplane Math Modeling Methods for Active Control Design[END_REF] or Karpel [START_REF] Karpel | Design for Active Flutter Suppression and Gust Alleviation Using State-Space Aeroelastic Modeling[END_REF] methods. This approach leads to very high order models for a given required accuracy level, and low numerical conditioning.

A more direct approach can be performed by working directly on the frequency response data H and identifying a state-space system. The Loewner method [START_REF] Lefteriu | A New Approach to Modeling Multiport Systems From Frequency-Domain Data[END_REF][START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF] is particularly adapted to MIMO systems with noise-free data such as the one studied in this work. For the sake of conciseness only the general procedure is outlined here, the interested reader can refer to the references for the detail of the equations. In the Loewner method, the data response H is first split in two subsets of data that are used to create an augmented state-space model of the form given in Eq. ( 6). This high-order model is then reduced to a selected order by singular value decomposition (SVD). The advantages of this method include a high accuracy for relatively low model order, a good numerical conditioning of the final system, the possibility to work with high amount of data and the simplicity of the procedure. The main drawback of the Loewner method is the impossibility to impose stability in the obtained state-space model. When unstable modes appear in the identified model, they are simply removed, which is an optimal stable approximation in the sense of the H 2 norm according to [START_REF] Köhler | On the closest stable descriptor system in the respective spaces RH2 and RH-infinity[END_REF]. As the aeroelastic system to identify is stable, this procedure does not deteriorate the modeling accuracy in a significant way.

In order to reduce the order of the final state-space model, only the frequency range that is relevant to the GLA problem is included in the input data. The frequency content of the wind rolls off at high frequency both for the discrete gust and continuous turbulence, and their impact on the aircraft response can be neglected after 5 Hz.

IV. Controller Synthesis by H ∞ and µ Techniques

Robust H ∞ and µ control syntheses are the main techniques used in this work. After a brief formulation of the theoretical framework shared by these two techniques, their application to the particular case of the flexible aircraft with a lidar sensor is performed. The linear dynamical models of actuators and sensors used in the synthesis are presented, and the methodology for obtaining performance and robustness in presence of wind disturbance and sensors noise is exposed.

A. Theoretical Framework

In H ∞ and µ analyses [START_REF] Doyle | Feedback Control Theory[END_REF][START_REF] Skogestad | Multivariable feedback control : analysis and design[END_REF], a general MIMO plant P :

      w u       →       z y      
is controlled by a feedback controller K : y → u as seen in Fig. 5. u is the control signal, y is the output variable composed of measured variables and non-measured loads, and w and z are the exogenous input and output respectively. The H ∞ and µ analyses focus on the frequency-dependant transfer function T w →z ( jω), and in particular on its H ∞ norm. This matrix norm is formally defined as the least upper bound over all frequencies of the maximum singular values of the transfer function:

T w →z ∞ = sup ω σ(T w →z ( jω)) (7)
where ω is the angular frequency and σ( A) is the maximum singular value of a generic matrix A. The signals u and y correspond to the actual signals used by the control system. The signals w and z do not necessarily correspond to physical quantities, the meaning of the H ∞ norm will depend on the choice of these variables. The H ∞ synthesis looks for a controller K : y → u that stabilizes the system and minimizes the H ∞ norm T w →z ∞ under certain constraints. Practically, by an adequate choice of P described in section IV.B, this minimization problem is turned into the following sub-optimal form: find a controller K that ensures

P

T w →z (K, P) ∞ ≤ 1 (8)
In Eq. ( 8), both the objective function and the constraints are written in a compact form, relying on a single scalar inequality. Exogenous input are typically noise, disturbance and output tracking reference, while exogenous output are typically the control effort and the plant output to minimize. From a given physical system, by an adequate choice of exogenous input and output and with additional gains and filters added to shape the desired closed-loop transfer functions, an augmented plant P is obtained. A controller K that satisfies the inequality (8) can be obtained by convex optimization using linear matrix inequalities or by solving the Riccati equations. It will be of the same order as the augmented plant P when H-infinity synthesis is used, and of higher order with mu-synthesis, and referred to as "full-order controller". A controller of reduced order set by the user can also be obtained, in which case the optimization In the H ∞ framework exposed above, the augmented plant P is assumed perfectly known. As it will be further explained in section IV.B, robustness to system's uncertainties can be obtained by an accurate choice of exogenous input and output. However, only unstructured uncertainties defined by full-order, upper norm-bounded dynamical systems, can be handled by the H ∞ framework. Imposing robustness to unstructured uncertainties is a restrictive constraint that can deteriorate the performance of the controller. Furthermore, when several uncertainties are considered in the H ∞ framework, they are treated as independent from one another [START_REF] Skogestad | Multivariable feedback control : analysis and design[END_REF]. The µ-synthesis framework [START_REF] Vidyasagar | Robust controllers for uncertain linear multivariable systems[END_REF] aims at reducing the same H ∞ norm as in the H ∞ framework, but in the presence of system uncertainties that can admit a fixed structure and that are applied at different inputs and outputs simultaneously [START_REF] Skogestad | Multivariable feedback control : analysis and design[END_REF]. The iterative D-K process [START_REF] Vidyasagar | Robust controllers for uncertain linear multivariable systems[END_REF] alternates K steps where H ∞ synthesis is performed, and D steps where the closed-loop robustness is assessed and the H ∞ norm to be reduced by the next K step is scaled according to the robustness performance. Although convergence is not ensured, this algorithm is numerically efficient.

B. H ∞ Controller Synthesis for the Flexible Aircraft equipped with Lidar

The H ∞ and µ syntheses frameworks described in section IV.A are applied to the flexible aircraft equipped with a lidar. The detail of the augmented plant P from Fig. 5 is shown in Fig. 6. By reducing the H ∞ norm of the full system, the influence of each exogenous input to the exogenous output is reduced. The augmented plant P shown in Fig. 5 consists of

• the aircraft aeroelastic model G defined in section III

• the actuators and sensors dynamics

• the lidar sensor noted as pure delay z i i=[0:N ] .

• the exogenous input w and output z and their associated tunable filters.

The uncertainties are written in a general form, they can be of different types as explained in section IV.B.4. The aircraft model G is normalized in input and output. In addition to improving the numerical convergence of the algorithms, this allows the comparison of the relative values of the entries of the performance output z y . The details of the different components of the augmented plant are shown in Fig. 6. The filters and exogenous inputs and outputs are chosen in such a way as to define the GLA problem defined in section II in terms of a performance criterion and associated constraints, as an H ∞ problem of the form [START_REF] Reschke | Flight loads analysis with inertially coupled equations of motion[END_REF]. The next sections explain their design in more details.

Gust Load Alleviation performance

The main objective of the controller synthesis as explained in section II is the gust load attenuation. It can be summarized as reducing the influence of the wind velocity v, taken here as the vertical wind alone, towards the loads.

Contrarily to y sens the output y l oads is not measured by the aircraft, hence it cannot be used by the controller. However, by including y l oads in the exogenous output z y , the controller is designed in such a way as to minimize the loads with the knowledge of the sensors output y sens only. The synthesis of a gust load alleviating controller is then performed as follows. The transfer function H v of Fig. 6 takes the frequency content of the worst case vertical wind. A third order approximation of a von Kármán turbulence with worst case amplitude is considered here. The numerical values can be found in section V.A. With an appropriate output normalization of variable y, H z y can be taken as the constant diagonal transfer function, with ones only at the output to minimize. As explained in section II.D, only a few selected loads are to be minimized. With this choice of filters H v and H z y , the minimization of the norm T w v →z y ∞ can be achieved by H ∞ synthesis. With the full delay z -N , the exogenous input w v corresponds, after multiplication by its shaping filter H v , to the predicted wind.

Actuators and Output Sensors

The actuators must be modeled accurately in order to account for the phase delay and the constraints they create.

The transfer function matrix H act of Fig. 6 represents the actuators dynamics and are modeled by a second order low-pass filter in the Laplace domain

H act (s) = 1 s 2 + 2ξω 0 s + ω 2 0 I N u (9)
where I N u is the identity matrix of size N u , the number of actuators, ω 0 and ξ are the actuators' cutoff frequency and damping ratio respectively. Together with the effect of the aircraft inertia, this limits the capability of the controller to attenuate the high frequency response of the aircraft.

In addition to the actuators dynamics, a saturation function is added in their model. The control surfaces cannot exceed a certain deflection angle threshold, and their derivative cannot exceed a certain deflection rate threshold. As it will be seen in the results section, the deflection rate saturation limit is the most constraining in the control synthesis problem. In the H ∞ framework, the maximum deflection angle and deflection rate are imposed through the filter

H zu (s) =          1 u ma x I N u s umax ( s+1) I N u          (10)
where s ( s+1) is a pseudo derivative with 1, and u ma x and umax are the maximum deflection angle and deflection rate respectively. Note that it has been assumed for simplicity that all the control surfaces have the same saturation limits. With an adequate input normalization, this ensures that when the H ∞ norm of the system is lower than one the saturation constraints are respected. In particular, T w v →z u ∞ must be lower than one with z u of dimension 2N u . Note that these filters are first written in the Laplace domain, and then converted into the z domain.

The sensors are modeled by H sens as a pure delay, accounting for the sensors measure rate and the on-board computer refresh rate:

H sens (s) = e -τ se n s s (11)
where τ sens is a diagonal matrix of size N y with delays associated to each sensor on the diagonal. When written in the z-transform domain, this transfer function becomes rational. The tunable filter H ny (see Fig. 6) associated with the measurement noise is taken as the diagonal matrix with expected standard deviation of each sensor's noise on the diagonal.

Lidar Sensor and Noise

As explained in section II, two types of lidar are considered in this work: the first type only performs a measurement of the vertical velocity at a fixed distance ahead of the aircraft. The measurement is then delayed several times and buffered as in Fig. 6, providing the information of the disturbance at all time steps ahead in time over a certain time horizon. The second type of lidar records also intermediate measurements. In the augmented plant, the intermediate measurements are performed at distances i δx corresponding to each time step i, where δx = V T s is the spatial step between two lidar measurements as shown in Fig. 2 and T s is the sampling time. In case of a perfect measure (no noise), this makes no difference. The main difference regards the measurement noise, which is assumed to increase linearly with the distance. This assumption is supported by experiments with different existing lidars [START_REF] Herbst | Design of a monolithic Michelson interferometer for fringe imaging in a near-field, UV, direct-detection Doppler wind lidar[END_REF]. With the first type of lidar (single measurement), the noise standard deviation at the i-th step ahead (time or spatial step equivalently) can be written as

σ i v = k n L ( 12
)
where L is the measure distance; k n is the noise per unit distance, which depends mainly on the lidar accuracy, and to a lesser extent on the the aircraft velocity and the atmosphere turbulence, as explained in section II. Note that the noise does not depend on i. In the second type of lidar (multiple measurements), the noise standard deviation can be written as

σ i v = i k n δx (13)
In this case, the closest measures give less time to react to the incoming wind, but they are more accurate. The standard deviation is then multiplied by a unit white noise. The knowledge of the noise amplitude is included in the H ∞ framework by selecting the lidar noise tunable filters as the constant H i nv = σ i v . The exogenous input w i nv of Fig. 6 can then be considered as white noises.

With the single-measurement lidar, increasing the length of measurement can become detrimental when the noise becomes too high and starts to overcome the performance benefit obtained by increasing the prediction horizon. This tradeoff will be studied in detail in section V.B.

Robustness with H ∞ and µ syntheses

The main aim of this work is the study of robustness in stability and performance. The former kind of robustness is the capability of the closed-loop to remain stable when the controlled system changes, whereas the latter concerns the guarantee to maintain performance under system variations, by keeping the H ∞ norm less than one. Ensuring robustness with a state space controller with several input and output is more challenging than with classical controllers consisting of gains and filters as used in aeronautics. The controller must be robust to MIMO uncertainties, which apply simultaneously to the various inputs and outputs.

One way to apply robustness requirements in the H ∞ framework is to consider input and output reference exogenous input, noted as w r e f ,u and w r e f ,y in Fig. 6. Those are the primary input of aircraft flight control, but their use in GLA is purely mathematical. The input and output sensitivity functions are first defined as the transfer functions S y = T y r e f → y and S u = T u r e f → u respectively. When considering inverse multiplicative system uncertainties in output G(s) = (I + ∆ y (s)) -1 G nom (s) and in input G(s) = G nom (s)(I + ∆ u (s)) -1 where G nom is the nominal system, ∆ y and ∆ u are unstructured uncertainties (upper-bounded uncertain dynamical systems) then it can be shown [START_REF] Doyle | Feedback Control Theory[END_REF] that if there exist scalars β y and β u such that

S y ∞ ≤ 1 β y ; S u ∞ ≤ 1 β u (14)
then the closed-loop remains stable for all uncertain systems ∆ y and ∆ u that verify

∆ y ∞ ≤ β y ; ∆ u ∞ ≤ β u ( 15
)
The parameters β y and β u are referred to as output and input modulus margins respectively. The infinity norm of the sensitivity functions S y and S u can be controlled directly by the tunable filter H r e f ,y and H r e f ,u of Fig. 6 respectively, taken as constants, and by the already mentioned H z y . This procedure defines the MIMO margins that can be handled by an appropriate choice of the tunable filters. Note that including the reference exogenous input for robustness will enable new interactions between exogenous input and output. For example, the norm of T y r e f →z u is included in the H ∞ synthesis, which is not directly wanted given the problem as stated. The robustness condition ( 15) is also practically difficult to ensure by a controller, and it can be too restrictive. Furthermore, the uncertainties hence defined do not apply simultaneously in input and output.

The other approach for ensuring robustness is µ-synthesis. A broad range of uncertainty structures can be directly defined in the augmented plant. The reference exogenous input w r e f ,u and w r e f ,y can then be ignored. In this work, a conservative design based on disk margins [START_REF] Blight | Practical control law design for aircraft using multivariable techniques[END_REF][START_REF] Seiler | An Introduction to Disk Margins [Lecture Notes[END_REF] is presented. Output symmetric disk uncertainties are multiplied by the nominal plant from the left. They are defined as

∆ y di sk =               ∆ 1 0 0 0 . . . 0 0 0 ∆ N y               (16)
with

∆ i = 1 + ηδ i 1 -ηδ i i = 1 . . . N y ( 17
)
where δ i are uncertain complex parameters lying inside the unit disk, N y is the number of outputs and η is a positive parameter defining the uncertainty level. Similar input uncertainties ∆ u di sk are also defined, multiplying from the right and with dimension N u . Note that for η = 0, the multiplicative uncertainty is equal to identity, meaning that the system is nominal. The minimum value of η that leads to unstable closed-loop will define the disk margin. These uncertainties are less restrictive than the general dynamical uncertainties defined in Eq. [START_REF] Albano | A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows[END_REF] of this type can be found. An advantage of this formulation is the fact that the disk margins can be used to compute guaranteed gain and phase margins, as explained in [START_REF] Seiler | An Introduction to Disk Margins [Lecture Notes[END_REF]. Disk uncertainties provide a more comprehensive assessment of robust stability as they account for simultaneous gain and phase uncertain variations, at different inputs and outputs at the same time. Moreover, µ-synthesis implies both stability and performance robustness which is of particular interest for this work.

In the following, SISO margins will refer to the classical gain and phase margins applied at one input or output at a time, while MIMO margins will refer to the robustness to simultaneous gain and phase variations at different inputs and outputs, guaranteed by disk margins.

V. Results: Robust Gust Load Alleviation of the Flexible Aircraft Equipped with Lidar under System Constraints

The methodology described in the previous sections is applied in this section. After a description of the numerical values used for the control design and the simulations in section V.A, section V.B investigates the alleviation potential for different lidar configurations and measurement distances. The alleviation of structural loads along the wing and in the HTP is then assessed with controllers of different orders in section V.C. Finally, a complete robustness analysis is performed in section V.D by the computation of disk margins and by studying the performance capabilities of the controllers when the system is subject to different types of uncertainties.

A. Simulations Parameters and Hardware considerations

The numerical values used in the simulations are presented in this section and summarized in table 2. The aircraft models are obtained from the system identification technique described in section III up to 5 Hz, due to the fact that the aircraft response to discrete gusts and continuous turbulence is negligible beyond this frequency. Some design filters are given in this work by their Laplace transform expression, but are eventually defined in the z-domain after time discretization of the system with time step T s equal to 0.02 s. This time step allows the simulation of the shortest gust, which lasts approximately 0.05 s. In the control design, the wind is modeled as a 3-rd order approximation of the von Kármán turbulence treated as a transfer function:

H v (s) = 0.02s 3 + 0.3398(L t ur b /V ) 2 s 2 + 2.7478(L t ur b /V )s + 1 0.1539(L t ur b /V ) 3 s 3 + 1.9754(L t ur b /V ) 2 s 2 + 2.9958(L t ur b /V )s + 1 ( 18
)
The true transfer function associated with the von Kármán PSD ( 2) is of high order and would increase the computational cost of the controller synthesis, so it is only used in simulation. During the controller synthesis, the output of the aircraft system are normalized by the standard deviation of the open loop response to continuous turbulence. All output variables shown in simulation are increments with respect to the trim equilibrium of the aircraft. Computations have been performed with a 8-core CPU and 32 GB of random access memory, in parallel mode when needed. The design is performed on Matlab thanks to the musyn function, and the simulations are performed on Simulink.

The measured output of the model are the pitch angle, pitch rate, vertical load factor and angle of attack from sensors.

The model also includes loads outputs that are not measured by the aircraft, composed of the shear force, bending moment and moment of torsion at 31 locations regularly spaced along the wing span starting at about 10%. According to the strategy exposed in section II, the control designs aims at minimizing the wing root bending moment, as well as the wing root shear force and moment of torsion but with a weight five times lower. The shear force, bending moment and torsion on the tail plane are also reduced, with the same weights as the corresponding loads on the wing. The other variables of interest are only presented in the results, but not included in the controller design. The twenty-four ailerons are used, and since only the longitudinal dynamics are considered, they are deflected symmetrically and grouped in twelve pairs.

B. Loads Alleviation with µ-synthesis for different lidar configurations

The influences of the lidar prediction length and noise are studied in this section. These simulations have been obtained for a 209th order robust controller (full order), corresponding to the order of the augmented plant P and additional states created by the µ-synthesis for robustness. Time simulations of the responses to a 150 ft discrete gust and to a continuous turbulence are shown in Fig. 7, where a lidar with multiple measurements up to 400 ft distance is included.

Actuator deflection rate

Inner Aileron Medium Aileron

Outer Aileron Elevator Fig. 7 Time response to a 150 ft discrete gust and to a von Kármán continuous turbulence of 2500 ft scale length and 34 ft/s amplitude with the full order robust controller the bending moment is strongly attenuated.

In Fig. 8 robust controllers have been obtained by µ-synthesis with a lidar that can measure the wind velocity at different intermediary lengths, one per time step. As described in section IV.B, each intermediary measure is noisier than its predecessor. The displayed performance corresponds to the reduction with respect to the open loop of the highest value of WRM in the discrete gust simulations, and of the standard deviation in the continuous turbulence case.

The standard deviation is computed by numerical integration of the frequency response. Each point corresponds to a different design, where the block modeling the lidar sensor is changed and the tuning parameters are adjusted in order to obtain an adequate design. In order to compare the performance, the parameters have been tuned in such a way as to obtain constant control efforts for all control designs. Stability robustness is also set as a constraint, by imposing MIMO margins of 6 dB gain and 30°phase. Without the saturation and robustness constraints, the loads could theoretically be totally canceled, and the maximum performance is highly influenced by this limit. Imposing such constraints on the control effort is needed for comparing the performance, as they are directly influenced by the maximum deflection rates allowed by the design. The µ-synthesis has been found to lead to better robustness than the H ∞ , and the stability margins can be more easily tuned thanks to the disk margins described in section IV.B. Imposing robustness constraints and a constant control effort in these results requires different iterations with varying design parameters to be performed, and has been done with parallel computing. Fig. 8 shows the benefit of using a lidar, the gain in performance between the 0 ft wind prediction and the asymptote is clear. It shows that beyond a 300 ft lidar measurement distance, adding information of the wind prediction does not lead to performance increase. Higher distances involve lower frequencies, and beyond a certain point, the lidar provides enough information to cover all the dynamical modes of the aircraft.

Moreover, the measurements at a high distance ahead of the aircraft are more noisy, hence they provide information that is less likely to improve the performance. It can be seen that the response to short gust is more difficult to alleviate than for long gusts. This is because high frequencies of the aircraft are excited, and due to the aircraft inertia and the controller limits, the high frequency response is more difficult to control. The same results have been obtained without ailerons, but are not shown for conciseness. They show that when only the elevators are used, the optimal measurement distance is about 700 ft, and the maximum alleviation is 19%. This demonstrates the advantage of including control surfaces creating direct lift such as the ailerons. Similar results have been obtained with the lidar that performs only one measure at a time, and can be seen in Fig. 9. Here, one measure of the wind vertical velocity is performed at the distance in the x-axis, and by conserving each measure the controller can predict the incoming wind profile in a similar way as with the lidar with intermediate measurements. The main difference lies in the lidar measurement noise as defined in section IV.B. Indeed in this case, for a given measurement distance L, each predicted wind velocity will suffer from the same noise, which is proportional to L. This time, increasing the prediction length has the disadvantage of reducing the accuracy of the predictions that are the closest in time. This can directly be seen on Fig. 9, where instead of a horizontal asymptote, a slow asymptotic reduction of the performance follows a performance peak at around 300 ft. One can note that the asymptotic decreasing rate is smaller for the 350 ft and the turbulent simulations than for the short discrete gusts. The short gusts are fully covered by a low lidar measurement distance, while longer gusts and turbulence have a higher spatial extension, needing higher lidar measurement distances to cover. Having a high number of noisy measurements representing the wind profile, the controller will more easily filter the noise out than with short gusts.

C. Multi-Load Alleviation with H ∞ and µ syntheses

In this section, the lidar measure distance is fixed to 400 ft and the lidar can perform intermediate measurements.

The focus is now on the multi-loads alleviation of the aircraft and the technique used for control synthesis. As explained in sections II and IV.B, minimizing all loads at the same time is a difficult task, and better overall results are obtained c.g.

Wing root

Horizontal tail plane

n z T z M x M y T z M x
M y Full order with lidar -20% -55% -57% -20% -28% -46% +21% Full order without lidar +19% -18% -22% +9% +15% -5% -13%

2nd order with lidar -18% -35% -36% -25% -33% -43% -27% Full order with lidar (no margin) -34% -64% -72% -23% -75% -69% +31% Table 3 Attenuation of the vertical load factor n z , shear force T z , bending moment M x and torsion torque M y at the wing root and horizontal tail plane due to a von Kármán continuous turbulence of 2500 ft scale length and 34 ft/s. Negative signs mean a reduction of the load when minimizing a few number of loads that have a particular interest, here mainly the bending moment at wing root but also the shear force and moment of torsion at wing root.

The robustness and loads reduction of the full-order robust controller are compared to other designs. A full-order robust controller without lidar is also designed for comparison. A reduced second-order controller is designed by structured H ∞ synthesis. Note that in the case of a low-order controller, robustness is easier to obtain as the controller does not have the possibility to perfectly fit the system to control as in the full-order case. Low order controllers are simpler and easier to implement on current generations of aircraft. The comparative performance at wing root and on the HTP obtained with the described controllers are shown in Tab. 3. It can be seen that the bending moment can be strongly alleviated (57% alleviation) whilst ensuring stability margins and system constraints. The lidar clearly helps alleviating the loads, the alleviation is about twice times higher than without lidar. The 2 nd -order reduced controller attenuates the bending moment by 36%. The bending moment on the tail plane is also successfully attenuated in presence of a lidar, both by the full-order and the reduced order robust controllers. A full-order controller is designed with H ∞ with no imposed robustness in its design and the results are shown in the last row of the table. With a reduction of the WRM of 72% compared to the 57% when robustness is imposed, this gives an idea of the maximum performance that can be achieved by a stabilizing linear controller, and shows the cost of the robustness constraints in terms of loads reduction.

The shear force T z , bending moment M x and moment of torsion M y along the wing span at 31 locations regularly spaced along the wing span are shown in Fig. 10. In the discrete gust case (first column) the y axis corresponds to the maximum values obtained during the time simulations, while in the continuous turbulence case (second column) it corresponds to the standard deviation computed by numerical integration of the frequency response. It can be seen that all loads decrease along the span, due to the accumulation effect of shear (linear) and moments (quadratic). It is clear that reducing the loads at the wing root directly reduces those further in the wing. Notice that a high reduction of the bending moment is obtained, and that the shear force directly follows the bending moment and is well alleviated too.

The moment of torsion is not always reduced and it can even increase in the discrete gust case. As explained in the section II, not all loads can be controlled at the same time. The aileron upward deflection creates a direct negative lift increment, which alleviates the shear force and bending moment, but also reduces the wing camber hence creating a pitch moment that can increase the moment of torsion. It can be seen that the 2 nd -order controller leads to good results for the continuous turbulence.

D. Robustness analysis

The robust controllers have been obtained by µ synthesis, including disk uncertainties in the augmented plant P as described in section IV.B. This means that an H ∞ criterion has been satisfied even in the worst case of the disk 17). Now, different tools will be used to validate the controller's robustness, in terms of both stability and performance.

The robustness in stability is assessed using disk uncertainties too. In the controller design, the range of simultaneous gain and phase variations at different inputs and outputs of the system was imposed, while in the validation process, the minimum variations leading to instability of the closed loop is calculated. The disk margin is defined from Eq. ( 17) as the minimum value of the parameter η leading to instability, and associated MIMO gain and phase margins can be computed. Tab. 4 shows the MIMO margins of the previously described controllers. Classical SISO margins, corresponding the minimum one-at-a-time gain and phase variations at a single input or output leading to instability, are also calculated. The SISO margins shown in Tab. 4 corresponds to the lowest among all the inputs and outputs.

It is first clear from the last row that when no robustness is required in the controller's design, the closed loop shows very low robustness margins. All the robust controllers present very high SISO margins (more than 11 dB and 76°g ain and phase margins respectively). All the MIMO margins are above 6 dB and 30°in gain and phase, meaning that any combination of uncertainties less than a factor 2 and with phase shift lower than 30°occurring simultaneously at different inputs and outputs of the system will lead to a stable closed loop.

The robustness in performance is then assessed in Fig. 11 by applying MIMO uncertainties of 4 dB gain and 30°to

the system and looking at the worst case performance. 40 realizations of these uncertainties are randomly generated and shown in green. A realization of the uncertain system corresponding to the highest H ∞ norm of the transfer function from the wind velocity to the WRM is constructed by the matlab function wcgain in a similar way as done in [START_REF] Patartics | Construction of an Uncertainty to Maximize the Gain at Multiple Frequencies[END_REF] and shown in red in the Fig. 11. Note that the maximum in the frequency-domain does not necessarily correspond to the maximum in the time-domain. Robustness to such uncertainties are the most difficult to obtain, as a particular worst case of simultaneous uncertainties may cause instability or loss of performance. It can be seen on the figure that the system remains stable in all cases of simultaneous uncertainties, which was expected since the uncertainties are lower than the closed-loop robustness margins. Moreover, very low variations in the performance are obtained. This robustness in performance ensures that the load alleviation remains good for small variations of the system and shows that the controller is not only valid for a single precise configuration.

Robustness in performance and stability to SISO uncertainties are studied through the particular case of actuators The uncertainties used in the previous paragraphs are defined by the user to cover different cases of the system's variations. More realistic uncertainties are also studied for the sake of completeness of the robustness analysis. In Fig. 13, the controller is simultaneously tested on different aeroelastic models of the aircraft, which correspond to different mass configurations. As explained in section II, very important variations in the aircraft model can exist along a flight, due to fuel depletion in the tanks and depending on the payload, so different mass cases used for calculation of extreme loads are considered. The total mass of the fully loaded aircraft is about twice that of the empty configuration. These results show realistic extreme cases of variations of the system that can occur during the flight, and are then of particular interest for certification. The maximum load with these variable mass cases is very close to the one obtained with the nominal mass case and the system remains stable, confirming the ability of the designed controller to alleviate the loads

VI. Conclusion

This work establishes a framework for gust load alleviation of a flexible aircraft equipped with a lidar sensors, applied to the XRF1 model developed at Airbus. This model features a high number of ailerons, and a relatively flexible wing. The first step is the identification of the high-dimensional, high-order aeroelastic system with the Loewner technique, which can deal with a high amount of data and leads to models of good accuracy and limited complexity.

The design of robust GLA laws taking advantage of the lidar measurements and of the high number of control surfaces is then performed with H ∞ and µ syntheses. The optimal lidar configuration has been studied, and a robustness analysis has been performed.

It results from the different test cases that an optimal lidar distance exists for the GLA of a flexible aircraft; this distance is about 300 ft for the case study of this paper. Beyond this distance, the measurements are too noisy to improve the performance when the lidar performs intermediate measurements, and deteriorate it when one measurement is performed. When only elevators were used, an optimal distance of 700 ft has been found instead. With the methodology developed in this work, the WRM can be alleviated up to 57% with a lidar, against 22% without. The response to short Open loop with nominal system Closed loop with different mass cases Closed loop with nominal system Fig. 13 Robustness to dynamic uncertainties: wing root bending moment time response to a 150 ft discrete gust for 8 different mass configurations of the aircraft gusts is more difficult to alleviate due to the higher frequencies involved, the percentage of reduction ranges from 20% for the 30 ft gust to 33% for the 350 ft one. The loads on the horizontal tail plane are also reduced, such as the bending moment which is divided by two with the proposed controller with respect to the open-loop case. The controller order can be drastically reduced (from about 200 down to 2) with dedicated structured synthesis, at the expense of the load alleviation that goes to 36%. The loads on the HTP are also successfully reduced, the bending moment being almost halved with respect to the open-loop case. These percentages of load reduction must be seen in the light of the fact that following European certification document CS-25 [START_REF]Certification Specifications for Large Aeroplanes, Amendment 3[END_REF], the aircraft structure should be able to withstand two third of the ultimate loads when the GLA system fails. This means that any alleviation beyond one third may not be accounted for in the aircraft design. Finally, the cost of imposing robustness can be evaluated, as a non-robust controller shows performance up to 72% load alleviation. The robust controller shows very low loss of performance when a sensor or an actuator fails, or when realistic variations of the mass configurations are applied. The most critical sensor appears to be the incidence probe, that measures the angle of attack at the nose of the aircraft.

The assumption of constant altitude and velocity has been made in this work, whereas most modern flight control laws can adapt in real time during a flight as the aerodynamics of the aircraft change in a significant way. The robustness studied in this work might not be enough to guarantee good performance during a whole flight. Future work will focus on finding solutions to obtain a robust stabilizing GLA controller performing in different flight conditions.

U

  gu st , L gu st = Gust amplitude and gradient distance U t ur b , L t ur b = Turbulence standard deviation and scale M hh , B hh , K hh = Modal mass, damping and stiffness matrices Q hh ,Q vh = Generalized aerodynamic modal and turbulent forces N u , N y = Number of aeroelastic inputs and outputs H = Frequency data response A, B u , B v ,C, D u , D v = State-space model matrices x,u, y = State, input and output vectors K, P = Controller and Augmented plant T α →β = Transfer function from generic vector α to β H act , H sens = Actuators and sensors transfer functions ω 0 , ξ = Actuators' cutoff frequency and damping ratio u ma x , umax = Maximum actuators deflection angle and rate τ sens = Sensors delays L = Lidar measurement distance k n = Lidar noise coefficient H zu , H z y , H v , H r e f ,u , H r e f ,y , H i nv , H ny = H ∞ tunable filters z u , z y = Exogenous output w v , w r e f ,u , w r e f ,y , w ny , w nv = Exogenous input S u , S y = Input and output sensitivity functions β u , β y = Input and output modulus margins ∆ u , ∆ y = Input and output uncertainties I. Introduction

Fig. 1

 1 Fig. 1 Finite element model of the XRF1 aircraft concept

Fig. 2

 2 Fig. 2 Remote wind velocity measurement with a lidar

Fig. 3

 3 Fig. 3 Discrete gust profiles of different gradient distances

Fig. 4 von

 4 Fig. 4 von Kármán turbulence profile of scale length 2500 ft and amplitude 85 ft/s

Fig. 5

 5 Fig. 5 General H ∞ formulation

Fig. 6

 6 Fig.6Augmented system used for the H ∞ and µ syntheses

Fig. 7 -

 7 b shows the wing root bending moment (WRM), whereas Fig. 7-c and Fig. 7-d show the deflection angle and rates of three ailerons out of the twelve present on a wing and the elevator. Positive values of the control surfaces angles are associated with downward deflections. The controller successfully attenuates the impact of the gust, the damped oscillations of the bending moment are typical for a flexible aircraft. Similarly, the impact of turbulence on 20

Fig. 8

 8 Fig. 8 Vertical load factor alleviation for a lidar with multiple intermediate measurements for varying length of measure

Fig. 9

 9 Fig. 9 Vertical load factor alleviation for a lidar with one measure for varying length of measure

Fig. 10

 10 Fig. 10 Shear force (a), bending moment (b) and moment of torsion (c) at different stations along the wing span, caused a 150 ft discrete gust and a 2500 ft von Kármán turbulence

Fig. 11

 11 Fig. 11 Robustness to MIMO uncertainties: wing root bending moment time response to a 150 ft discrete gust with 6dB gain and 30°phase uncertainties, applied simultaneously on all inputs and outputs with 40 random samples

  loop with one faulty sensor Closed loop with nominal system

Fig. 12

 12 Fig. 12 Robustness to SISO uncertainties: wing root bending moment time response to a 150 ft discrete gust with failure of a single actuator (a) or sensor (b) at a time

  The wind not only has a static influence on the aircraft by applying static aerodynamic loads,

									30ft gust
		25							150ft gust
									350ft gust
	Vertical Wind Velocity(ft/s)	10 15 20
		5
		0
		0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
						Time [s]

Table 2

 2 Numerical values used in section V

	Quantity	Symbol	Value
	Mach number	M	0.85
	Altitude	-	30.000 ft
	Actuators cutoff frequency	ω 0	4 Hz
	Actuators damping ratio	ξ	0.85
	Maximum deflection angle	u ma x	30°M
	aximum deflection rate	umax	40 °/s
	Sensors delays	τ sens	40 ms for each sensor
	Lidar noise coefficient	k n	0.02 s -1
	Discrete gusts gradient	L gust	30 ft, 150 ft and 350 ft
	Scale of turbulence	L t ur b	2500 ft
	Turbulence standard deviation	U σ	34 ft/s
	Model order	-	50
	Time step	T s	0.02 s

and they can apply simultaneously on the various input and outputs. With µ-synthesis, a controller that ensures robustness to simultaneous MIMO uncertainties

Acknowledgments

This work has been funded by Airbus Operations SAS and the French National Agency for Technological Research (ANRT) under grant number 2019/1611. This support is gratefully appreciated.


Adrian Ion
email: aion@prip.tuwien.ac.at

Gabriel Peyré
email: bgabriel.peyre@ceremade.dauphine.fr

Yll Haxhimusa

Samuel Peltier

Walter G Kropatsch

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Shape Matching Using the Geodesic Eccentricity Transform -A Study 1)

This paper makes use of the continuous eccentricity transform to perform shape matching. The eccentricity transform has already been proven useful in a discrete graph-theoretic setting. We show how these ideas extend naturally to the continuous setting thus bringing a higher geometrical fidelity. The continuous eccentricity transform is used to compute multiscale descriptors for shapes. These descriptors are defined as histograms of the eccentricity transform of a scale-space representation of the shape. These multiscale descriptors are naturally invariant to euclidean motion and bending. They show promising results for shape discrimination.

Introduction

Shape recognition is a central topic in computer vision. It requires to set up a signature that characterizes the properties of interest for the recognition [START_REF] Veltkamp | Properties and performance of shape similarity measures[END_REF]. The invariance of this signature to local deformations such as articulations is important for the identification of 2D shapes. Matching can then be carried out over this reduced space of signatures.

Most shape descriptors are computed over a transformed domain that amplifies the important features of the object while throwing away ambiguities such as translation, rotation or local deformations. The Fourier transform of the boundary curve [START_REF] Zahn | Fourier descriptors for plane closed curves[END_REF] is an example of such transformed domain descriptor adapted to smooth shapes. Shape transformations computed with geodesic distances [START_REF] Bronstein | Matching two-dimensional articulated shapes using generalized multidimensional scaling[END_REF] lead to signatures invariant to isometric deformations such as bending or articulation. To capture salient features of objects, local quantities such as curvature [START_REF] Mokhtarian | A theory of multiscale, curvature-based shape representation for planar curves[END_REF] or shape contexts [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] can be computed. They can be extended to bending invariant signatures using geodesic distances [START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF]. More global features include the Laplace spectra [START_REF] Reuter | Laplace-spectra as fingerprints for shape matching[END_REF] and the skeleton [START_REF] Siddiqi | Shock graphs and shape matching[END_REF]. Some transformations involve the computation of a 2D function defined on the shape, for instance the solution to a linear partial differential equation [START_REF] Gorelick | Shape representation and classification using the poisson equation[END_REF] or geometric quantities [START_REF] Osada | Shape distributions[END_REF]. One can also use geodesic distance information such as the mean-geodesic transform [START_REF] Hamza | Geodesic object representation and recognition[END_REF]. In this paper, we propose to use the eccentricity transform [START_REF] Kropatsch | The eccentricity transform (of a digital shape)[END_REF], computed using continuous geodesic distances [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], for shape matching. In Sec. 2 the eccentricity transform is presented. Sec. 3 presents our eccentricity histogram descriptors. Experimental results are shown in Sec. 4.

A shape is a planar, connected, bounded and closed set S ⊂ R 2 with a piecewise-smooth boundary ∂S. In practice we consider a discretized version of S which can be represented using an image f S of n pixels where f S is the indicator of the shape: f S (x) = 1 for x ∈ S and f S (x) = 0 otherwise. Such a discretized representation is usually acquired from some digitized continuous shape.

The geodesic inner distance between two points x, y ∈ S is defined as

d S (x, y) def. = min γ∈P(x,y) L(γ)
where L(γ)

def. = d S (x 0 , x) to some starting point x 0 ∈ S can be computed efficiently as the solution of the non-linear Eikonal equation:

∀ x ∈ S, ||∇U(x)|| = 1 and U(0) = 0. ( 1
)
We note that this resolution is restricted to S, which makes the function U different from the traditional euclidean distance function x → ||x 0 -x|| for non-convex shapes. Fig. 1 shows a comparison of the geodesic and euclidean distances. The Eikonal equation ( 1) can be solved in O(n log(n)) operations for a grid of n points using the Fast Marching algorithm [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. The eccentricity transform of a shape S, assigns to each point p ∈ S the shortest geodesic distance to the point of S farthest away from it. The eccentricity of the shape S is defined as

ECC S (x) def. = max y∈S d S (x, y) = max y∈∂S d S (x, y). ( 2
)
This is a continuous and piecewise smooth function. A point y that reaches the global maximum in (2) is called eccentric. We denote E(S) the set of eccentric points. An important property due to [START_REF] Kropatsch | The eccentricity transform (of a digital shape)[END_REF] is that E(S) is included in the boundary of S: E(S) ⊂ ∂S. The set of eccentric points allows to define a segmentation of S into eccentric regions

S = x∈E(S) A x where A x = {y \ ECC S (y) = d S (y, x)} .
The eccentricity is computed by performing a Fast Marching propagation from each x ∈ ∂S in order to compute the set of distances {d(x, y)} x∈∂S . This method requires O(|∂S|n log(n)) operations where |∂S| is the number of pixels on the boundary of S. Fig. 2 shows the eccentricity transform of two shapes.

The eccentricity transform has already been studied on a discrete graph [START_REF] Kropatsch | The eccentricity transform (of a digital shape)[END_REF]. The set of values {ECC S (x)} x∈S is invariant under rigid motion and isometric transform of S, which includes bending. It is also nearly-invariant under articulations, see [START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF][START_REF] Kropatsch | The eccentricity transform (of a digital shape)[END_REF].

Another property of ECC S is that it is robust to salt and pepper noise that might create holes in S, because a hole or a segmentation error of size ε only modifies the eccentricity ECC S by no more than ε. It is very different from local descriptors such as the curvature [START_REF] Mokhtarian | A theory of multiscale, curvature-based shape representation for planar curves[END_REF] or even global ones such as the structure of the skeleton [START_REF] Siddiqi | Shock graphs and shape matching[END_REF] which are not robust to this kind of noise.

Eccentricity Histogram Matching

In order to match two shapes from two binary images we first create a shape descriptor for each of them and then match these descriptors to obtain a similarity measure. This paper proposes two approaches based on a single and a multiscale descriptor.

Mono-scale descriptor. The basic building block for our shape descriptor is the histogram h S of the eccentricity transform ECC S of the shape S. We use m bins to estimate the histogram and in numerical applications, we use m = 200. The histogram descriptor is then the vector h S ∈ R m defined by where |S| is the number of pixels in S. Fig. 3 shows examples of histograms for shapes with different geometric features. We note that the histogram h S is invariant under euclidean transformations, scaling and isometric bending of S.

∀ i = 1, . . . , m, h S (i) = 1 |S| # x ∈ S \ i -1 m ECC S (x) -min(ECC S ) max(ECC S ) -min(ECC S ) < i m ,
Multiscale descriptor. In order to capture more geometric information about a shape S, we compute a non-linear scale-space of S and extract the histograms of the eccentricity over a scale-space domain.

In order to smooth the shape, we perform the following mean curvature evolution [13] 1)

∂γ t ∂t (u) = κ t (u) n s (u) with γ 0 (u) = ∂S(u),
where κ t (u) is the curvature of the curve γ t (u) and n t (u) is the normal vector to the curve. The curve γ s is thus a smoothed version of the boundary of S after a diffusion during a time t. Fig. 4, left, shows the process of shape smoothing.

This curve smoothing allows to define a set of shapes {S k } K k=1 where S k is the shape whose boundary is ∂S k def.

= γ t k with dyadic time steps t k = τ 2 k . The number of shapes is set to K = 3 in our numerical tests and τ is of the order of 5 pixels.

Our shape descriptor is composed of the histograms of the eccentricity transform over the scale space domain

D S def. = {h S 1 , . . . , h S k }.
Fig. 4 shows an example of descriptor. During smoothing, elongated parts tend to disappear, which affects the histograms.

Comparison of histograms.

In order to match two descriptors D S and D S of two shapes S and S, we need to compute a meaningful distance between histograms. Let h ∈ R m and h ∈ R m be two histograms. We propose to use the simple L 2 -norm defined by

δ(h, h) def. = m i=1 (h(i) -h(i)) 2 .
1) Note that more elaborate methods, like the one presented in [START_REF] Mokhtarian | A theory of multiscale, curvature-based shape representation for planar curves[END_REF], could be used. One could use more elaborate metrics such as the χ 2 metric, the earth mover distance or those defined in [START_REF] Osada | Shape distributions[END_REF], but we found in numerical experiments that all these metrics give results similar to δ.

We can compute the distance Δ(S, S) between two shapes S, S as the weighted sum of the distance of their histogram descriptors D S = {h S 1 , . . . , h S K } and D S = {h S1 , . . . , h SK } Δ(S, S)

def. = K k=1 w k δ(h S k , h Sk ).
In our experiments we used

w k = 1 2 k-1 .

Experiments

For the experiments we have used three shape databases: Kimia 25 [START_REF] Sharvit | symmetry-based indexing of image databases[END_REF], Kimia 99 [START_REF] Sebastian | Recognition of shapes by editing their shock graphs[END_REF] and MPEG7 CE-Shape-1 [START_REF] Latecki | Shape descriptors for non-rigid shapes with a single closed contour[END_REF]. For the Kimia 25 and Kimia 99 databases, the eccentricity ECC S is computed directly on the images as shown in section 2. For the MPEG 7 database, which contains 1400 images with varying sizes, we rescale the images to 256 × 256 pixels.

A shape database is composed of q shapes {S i } q i=1 and each shape S i has a label (i) ∈ {1, . . . , max }. Each label value 1 max defines a class of shapes C def.

= {S i \ (i) = }. The left column of the three blocks of Fig. 5 shows some shapes from the Kimia database, ordered by classes (such as fish, planes, rabbits, etc.). Any shape matching algorithm α assigns to each shape S i a vector of best matches Φ i where Φ i (1) is the shape the more similar to S i , Φ i (2) is the second hit, and so on.

For the Kimia 99 database, max = 9 and q = 99. We measure the efficiency of various matching algorithms on Kimia databases by the number of correct matches for each ranking position k

Match k (Φ) def. = q i=1 1 (i)= (Φ i (k)) q.
Tables 1 and2 give the value of Match k for various shape matching algorithms.

In the case of the MPEG7 database, which contains max = 70 classes with 20 images each, so q = 70 × 20, the efficiency of matching algorithms is computed using the standard Bullseye

def. = 1 20q k 40 q i=1 1 (i)= (Φ i (k)) = 1 20q k 40 Match k (Φ) 1.
For each image there can be at most 20 correct hits i.e. a maximum of 20 × 1400 hits can be obtained and thus Bullseye(Φ) 1. Table 3 gives the value of Bullseye for various shape matching algorithms.

The overall results over both Kimia 25 and Kimia 99, and over MPEG 7 are slightly bellow the state of the art. It is important to consider that the proposed methods are simple since they only require the computation of geodesic distances and histograms. In contrast, the most efficient shape matching algorithms [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF][START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF] are more complicated and require extraction of salient features and local signatures that need to be aligned or registered.

Case Study -Kimia 25 Fig. 5 shows the retrieval results for Kimia 25 using the single scale descriptor. The first column shows the 25 shapes S i . The following set of shapes forms an array where the shape at row i and column k is Φ i (k), the rank-k shape associated to S i .

The airplanes and greebles (the shapes that do not look like anything we know). Two questions immediately rise when looking at these results:

1. Why are the greebles considered to be more similar to the hands than to other greebles? 2. Why does a rabbit appear in so many cases when the matching has failed?

For the first question, look at the histograms of the greebles and the unoccluded hands (see Fig. 6.a). The histograms are very similar even though the shapes are different, e.g. the histogram of the first greeble (Fig. 6.a top-left) looks more similar to the hands, than the second and third greeble. This is due to the abstraction of a 2D shape to a 1D histogram, which disregards the structure of distances/paths.

For the second question, we consider the shapes in Fig. Geometrical properties of the shapes are well captured by our low-dimensional descriptors. For instance, elongated shapes are well separated from more compact shapes. However, more advanced geometrical features, such as intricate structural properties are thrown away by our signature extraction. This is for instance why the class "greebles" is not separated enough from the class "hands". In this paper we have presented the continuous eccentricity transform. We studied it's application to shape matching using descriptors computed with the continuous eccentricity transform. These descriptors are defined as histograms of the eccentricity transform of a scale-space representation of the shape. Results of the proposed method on three well known databases are presented and compared with the state of the art methods. We plan to study decomposition of articulated shapes using the eccentricity transform and extend the current approach to part based shape matching.

= 1 0

 1 |γ (t)|dt, where P(x, y) is the set of paths γ(t) ⊂ S joining x to y: γ(0) = x and γ(1) = y. The computation of the distance function U(x) def.

Figure 1 :

 1 Figure 1: (a,c) Two examples of euclidean distance functions ||x 0 -x|| restricted to S. (b,d) The geodesic distances d(x 0 , x) to some starting point x 0 .

Figure 2 :

 2 Figure 2: Left: the binary shape indicator function f S for two different shapes. Center: eccentricity transform ECC S , the values range from 0 to max x ECC S (x). The eccentric center C(S) = argmin x E S (x) is displayed as a point. Right: eccentric points E(S) together with the segmentation in eccentric regions A x for x ∈ E(S).

Figure 3 :

 3 Figure 3: Top: eccentricity transform ECC S for some shapes S used during experiments. Bottom: corresponding histograms h S .

Figure 4 :

 4 Figure 4: Upper left : smoothed shapes S k for 0 k 5. Bottom left : corresponding eccentricity transform ECC S k . Right : overlay of the histograms h S k .

Figure 5 :

 5 Figure 5: Retrieval results for the single scale descriptor on the Kimia 25 database.

  [START_REF] Latecki | Shape descriptors for non-rigid shapes with a single closed contour[END_REF].b (a rabbit -S 19 , and two tools -S 25 and S 22 ), and the results, Φ 25 , in the bottom row of Fig.5. When matching S 25 , the rabbit has a better score than S 22 , even though one might say that the histograms of S 25 and S 22 reveal more similar distance characteristics than the histogram of S 19 (see Fig.6.b) i.e. both S 25 and S 22 have more long distances than medium, and short, while S 19 has a peak in the medium. This is due to histogram matching methods, inherently low level, failing to capture the high level context of the task. Besides the L 2 -norm we have also looked at the χ 2 statistic[START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF] and diffusion distance[START_REF] Ling | Diffusion distance for histogram comparison[END_REF] and both produce similar results.

Figure 6 :

 6 Figure 6: (a) Histograms for: top: greebles, and bottom: unoccluded hands. (b) top: Three shapes from the Kimia 25 database, and bottom: their eccentricity histograms.

Table 1 :

 1 The value of Match k (Φ) for

	Algorithm α	k=1 2	3
	Sharvit et. al	23	21 20
	ECC-Hist	25	20 16
	ECC-M.-Hist	25	21 15
	Gdalyahu and Weinshall	25	21 19
	Belongie et. al	25	24 22
	ID-Shape Context [1]	25	24 25

various algorithms on the Kimia 25 database. See [7] for a description of these algorithms.

	Algorithm α	k=1 2	3	4	5	6	7	8	9 10
	ECC-Hist	99	87 74 67 64 49 52 45 38 33
	ECC-M.-Hist	99	88 79 68 72 60 55 44 43 40
	Shape Context	97	91 88 85 84 77 75 66 56 37
	Gen. Model	99	97 99 98 96 96 94 83 75 48
	Shock Edit	99	99 99 98 98 97 96 95 93 82
	ID-Shape Context [1]	99	99 99 98 98 97 97 98 94 79

Table 2 : The value of Match

 2

k (Φ) for

various algorithms on the Kimia 99 database. See [7] for a description of these algorithms

test. This test counts the number of correct hits (same class) in the first 40 hits Bullseye(Φ)

Table 3 : The value of Bullseye(Φ) for various algorithms on the MPEG 7 databases.

 3

[START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] Partially supported by the Austrian Science Fund under grants S9103-N04 and P18716-N13.