# Interpolatory Methods for Generic BizJet Gust Load Alleviation Function

C. Poussot-Vassal, P. Vuillemin, O. Cantinaud and F. Sève C. Poussot-Vassal
and P. Vuillemin are with Université de Toulouse, F-31055 Toulouse, France.
Contact<EMAIL_ADDRESS>O. Cantinaud and F. Sève are with
Dassault-Aviation, France.This work has been funded within the frame of the
Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology
Demonstrator platform AIRFRAME ITD (contract N. CSJU-CS2-GAM-AIR-2014-15-01
Annex 1, Issue B04, October 2nd, 2015) being part of the Horizon 2020 research
and Innovation framework program of the European Commission.

###### Abstract

The paper’s main contribution concerns the use of interpolatory methods to
solve end to end industrial control problems involving complex linear
dynamical systems. More in details, contributions show how the rational data
and function interpolation framework is a pivotal tool _(i)_ to construct
(frequency-limited) reduced order dynamical models appropriate for model-based
control design and _(ii)_ to accurately discretise controllers in view of on-
board computer-limited implementation. These contributions are illustrated
along the paper through the design of an active feedback gust load alleviation
function, applied on an industrial generic business jet aircraft use-case. The
closed-loop validation and performances evaluation are assessed through the
use of an industrial dedicated simulator and considering certification
objectives. Although application is centred on aircraft applications, the
method is not restrictive and can be extended to any linear dynamical systems.

##  1 Introduction

###  1.1 General aircraft context

Aircraft mobility plays an important role in our life style and societal
organisation. As this transportation mean is facing severe environmental and
societal challenges (_e.g._ global warming and CO2 emissions), it is the role
of researchers to provide innovative answers and industry-oriented tools to
address them. These solutions should fulfil the safety and design requirements
and be numerically efficient and simple to implement within the aircraft
industrial design value chain. Among others scientific disciplines, it is well
admitted that civil aircraft industry relies on dynamical systems theory,
linear algebra and computational sciences to address these issues. In light of
these statements, recent developments in these scientific communities may have
a major impact in the overall aircraft conception and exploitation enhancement
[2, 45]. Hopefully, it can also be a substantial lever in the _civil aviation
footprint reduction_. This is one of the main expectation of this work.

In order to reach these ambitious objectives, the paper highlights the pivotal
role and relevance of _interpolatory methods_ [4, 5] in the context of _linear
time invariant ( LTI) large-scale systems_ [34].

###  1.2 Generic business jet aircraft context

More specifically, through a complete industrial gust load alleviation control
design and validation problem applied to a generic business jet aircraft, we
show how rational function interpolation is an appealing ingredient for
engineers. It is used _(i)_ to construct _reduced order dynamical models_
appropriate for feedback control design and _(ii)_ to compute more accurate
_sampled-time_ gust load dynamical controller functions that can be
implemented in a constrained on-board computer. The complete closed-loop
stability and performance analysis is done through a dedicated _industrial
simulator_ to assess the approach.

###  1.3 Generic business gust load alleviation feedback control objective

An important aircraft design criteria concerns the _gust load envelope_ that
should not reach some limit to ensure structure integrity. To this aim, it is
standard considering vertical gusts wavelengths, modelled through the
following ”1-cosine” profile

$\mathbf{w}_{g}(t)=\left\\{\begin{array}[]{ll}\dfrac{W}{2}\left(1-\cos\left(\dfrac{\pi
V}{L}t\right)\right)&\text{ for $0\leq t\leq\dfrac{2L}{V}$}\\\ 0&\text{ for
$0>\dfrac{2L}{V}$}\end{array}\right.,$ (1)

where $W$ is the gust velocity (in feet), $L$ is the gust wavelength (in
meters) and $V$ is the aircraft true airspeed (in meter per seconds)111Typical
values for these parameters are provided by authorities and results in
hundreds of different gust configurations.. The _gust load envelope_ is simply
the worst case load responses along the wing span in reaction to the set of
time-domain vertical wind gust profiles eq. 1 affecting the aircraft.
Similarly to aerodynamical effects inducing vibrations [25], these gusts
represent typical phenomena the aircraft might face during its exploitation.
Prior any flight test, aircraft manufacturers should guarantee authorities
that they can be handled by the system.

In the preliminary conception step (before considering any control functions),
the aircraft is designed by experts so that the wings should support a given
nominal load envelope, dictated by physical considerations such as desired
aircraft manoeuvrability, gust, and many other manufacturing constraints. The
larger the supported loads are, the larger the structural stiffeners and mass
reinforcements should be. The aircraft mass is consequently bigger and its
consumption during flight increased. In this context, the _Gust Load
Alleviation_ (GLA) control function plays an important role in the aircraft
conception: it is aimed at lowering the loads envelope and thus at reducing
the aircraft overall mass, consumption and emissions [1, 14]. In this work,
the gain brought by the GLA, on the _worst upward gain_ is denoted by
$\mathcal{E}(x_{i})$, is computed as

$\mathcal{E}(x_{i})=\max_{\mathbf{w}_{g}\in\cal W\\\
}\left\lVert\dfrac{\mathbf{z}_{\text{loads}}(x_{i})-\mathbf{z}^{\text{(GLA)}}_{\text{loads}}(x_{i})}{\mathbf{z}_{\text{loads}}(x_{i})}\right\rVert_{\infty},$
(2)

where $x_{i}$ (in meters) is the wing location points where the loads are
computed (here five locations $x_{i}=[x_{1},x_{2},x_{3},x_{4},x_{5}]$m) and
where $\cal W$ is the set of gust profiles as in (1). Then,
$\mathbf{z}_{\text{loads}}(x_{i})$ denotes the load at location $x_{i}$ for
the baseline aircraft and $\mathbf{z}^{(GLA)}_{\text{loads}}(x_{i})$ denotes
the one when the _Gust Load Alleviation_ (GLA) control function is activated.
Both responses are obtained when the aircraft model is fed by eq. 1 (see also
section 2).

###  1.4 Generic business gust load alleviation feedback industrial
constraints

In an industrial value chain, the GLA function is traditionally designed after
the flight controller (focusing on handling qualities). Then, one requirement
is that GLA control should not impact the nominal (low frequencies) flight
performances and focus on the gust phenomena and load envelope only. In
addition, as the control functions should be implemented in a limited sampled-
time computer, with a constrained material architecture (sensors and actuators
limitations, sampling limitations, delays in the loop, _etc._ ), the
discretisation step should also be accurately taken into account before any
implementation and performance evaluation.

###  1.5 Main results and paper organisation

Given the above considerations and objectives, the rest of the paper is
focused on the central role of the _interpolation_ applied to the _gust load
alleviation control_ problem. The paper is organised as follows. The gust
load-oriented aircraft modelling is detailed in section 2, leading to a set of
medium-scale irrational models. As the main mathematical tool involved here,
reminders on the interpolation framework are recalled in section 3. Its
application for the construction of rational reduced order models set,
suitable for linear control design, is described in section 4. Then, after
briefly describing the continuous-time GLA controller synthesis and the
continuous and sampled-time interconnection framework, the controller
_discretisation_ through the interpolation framework is detailed in section 5.
This last point stands as a methodological contribution of the work. Finally,
conclusions are given in section 6, illustrating first the efficiency of the
proposed process to alleviate gust-driven loads in an industrial application
context, and second, to summarise the main contributions that can be applied
to any kind of linear dynamical models.

Although the paper is centred around the aircraft gust load function, the
proposed approach is readily applicable to any other aircraft control problems
(including aircraft flight quality, vibrations, _etc._ functions). It is also
valid to (m)any linear time invariant dynamical use-cases. Indeed, the
presented _interpolatory framework_ versatility allows to address a large
number of engineering applications.

###  1.6 Notations

Let us denote by $\mathbb{R}$ the set of real numbers, $\mathbb{C}$ the set of
complex numbers, $\mathbb{C}_{+}$ ($\mathbb{C}_{-}$) the open right (left)
half plane, $\mathcal{D}$ the open unit disk, $\partial\mathcal{D}$ its
boundary and $\overline{\mathcal{D}}$ the complementary of the closed unit
disk, respectively. The complex variable is given by $\imath=\sqrt{-1}$. Let
$\mathcal{L}_{2}(\mathcal{I})$ ($\mathcal{I}=\imath\mathbb{R}$ or
$\partial\mathcal{D}$) be the set of functions that are square integrable on
$\mathcal{I}$. Let $\mathcal{H}_{2}(\mathcal{\mathcal{D}})$ (resp.
$\mathcal{H}_{2}(\mathcal{\overline{\mathcal{D}}})$) be the subset of
$\mathcal{L}_{2}(\partial\mathcal{D})$ containing the functions analytic in
$\mathcal{D}$ (resp. $\overline{\mathcal{D}}$). Let
${\mathcal{H}_{2}}(\mathbb{C}_{+})$, shortly ${\mathcal{H}_{2}}$, (resp.
$\mathcal{H}_{2}(\mathbb{C}_{-})$) be the subset of
$\mathcal{L}_{2}(\imath\mathbb{R})$ containing the functions analytic in
$\mathbb{C}_{+}$ (resp. $\mathbb{C}_{-}$). Similarly, let
$\mathcal{L}_{\infty}(\mathcal{I})$ ($\mathcal{I}=\imath\mathbb{R}$ or
$\partial\mathcal{D}$) be the set of functions that are bounded on
$\mathcal{I}$. Let $\mathcal{H}_{\infty}(\mathcal{\mathcal{D}})$ (resp.
$\mathcal{H}_{\infty}(\mathcal{\overline{\mathcal{D}}})$) be the subset of
$\mathcal{L}_{\infty}(\partial\mathcal{D})$ containing the functions analytic
in $\mathcal{D}$ (resp. $\overline{\mathcal{D}}$) and
$\mathcal{H}_{\infty}(\mathbb{C}_{+})$, shortly ${\mathcal{H}_{\infty}}$, the
subset of $\mathcal{L}_{\infty}(\imath\mathbb{R})$ of functions analytic in
$\mathbb{C}_{+}$. The Fourier transform of a time-domain signal
$v\in\mathcal{L}_{2}(\mathbb{R})$ is denoted by $\overline{v}=\mathcal{F}(v)$.

##  2 Gust load oriented modelling

fig. 1 illustrates the considered feedback control loop architecture, where
the sampled-time GLA function denoted as ”GLA controller” is to be computed
and validated. One may note that the actuators, sensors, flight controller and
computational delay are considered as given and included in the considered
generic business jet (BizJet) aircraft model set, referred to as
$\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$, and detailed hereafter.

Actuators Sensors Computational delay $\tau_{m}$
$\mathbf{y}$$\mathbf{z}$$\mathbf{u}$$\mathbf{w}_{G}$ Flight controller Sampled
GLA controller $h$Generic BizJet aircraft $\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$

Figure 1: Closed-loop architecture of the GLA problem. The complete
aeroservoelastic dynamical aircraft models
$\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$ include the ”Flight controller”,
”Actuators”, ”Sensors” and ”Computational delay $\tau_{m}$”. The ”GLA
controller” is GLA function to be computed. Signals $\mathbf{w}$,
$\mathbf{u}$, $\mathbf{z}$ and $\mathbf{y}$ denote the exogenous inputs,
control inputs, performance outputs and measurements, respectively. Then $h$
denotes the sampling time usd in the GLA function.

In section 2.1, the large-scale aeroservoelastic models
$\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$, as provided by physics and loads teams,
is described. Its input-output simplification is then done in section 2.2,
leading to the the set of irrational models
$\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$.

###  2.1 Aeroservoelastic aircraft modelling

To guarantee the safety and performance levels, GLA function is developed
using very accurate environment and models, issued from aeroservoelastic
modelling teams [31, 33, 32, 20].

Following fig. 1, $\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$ denotes the gust load
oriented dynamical models family complex-valued functions. Each
$\mathbf{G}_{i}$ ($i=1,\dots,n_{s}$) represents a _linear time invariant
dynamical model_ of the aircraft, evaluated at a given flight and mass
condition. Its construction is a know-how of the aircraft manufacturer and
usually results from multiple steps performed by different experts from fluid
mechanics, structural, tests _etc._ teams. Each model
$\\{\mathbf{G}_{i}\\}_{i=1}^{n_{s}}$ is described by the following realisation
denoted by $\\{\mathcal{G}_{i}\\}_{i=1}^{n_{s}}$

$\\{\mathcal{G}_{i}\\}_{i=1}^{n_{s}}:\left\\{\begin{array}[]{rcl}\mathbf{\dot{x}}^{(i)}(t)&=&{A}^{(i)}\mathbf{x}^{(i)}(t)+{B}^{(i)}\left[\begin{array}[]{c}\mathbf{w}_{G}(t)\\\
\mathbf{u}(t)\end{array}\right]\\\ \left[\begin{array}[]{c}\mathbf{y}(t)\\\
\mathbf{z}(t)\end{array}\right]&=&{C}^{(i)}\mathbf{x}^{(i)}(t)+{C}_{d}^{(i)}\mathbf{x}^{(i)}(t-\tau_{m})\end{array}\right.$
(3)

where $\mathbf{x}^{(i)}(t)\in\mathbb{R}^{M_{i}}$,
$\mathbf{u}(t)\in\mathbb{R}^{n_{u}}$,
$\mathbf{w}_{G}(t)\in\mathbb{R}^{n_{wG}}$,
$\mathbf{y}(t)\in\mathbb{R}^{n_{y}}$ and $\mathbf{z}(t)\in\mathbb{R}^{n_{z}}$
are the internal variables, input, exogenous input, output and performance
output signals, respectively. The dimensions $n_{u}$, $n_{wG}$, $n_{y}$ and
$n_{z}$ are constant for all $n_{s}$ models, while $M_{i}$ depends on the
model configuration ($i=1,\dots,n_{s}$). Finally, $\tau_{m}\in\mathbb{R}_{+}$
represents the computational delay needed to capture the $\mathbf{y}$ signal
(and to compute the - to be designed - GLA control law).

In the considered use-case, $\mathbf{u}$ gathers three control inputs, the
horizontal tail, the inner and outer ailerons ($n_{u}=3$). Exogenous input
$\mathbf{w}_{G}$ gathers $\mathbf{\delta}_{mc}^{\star}$, the aircraft
horizontal tail deflection given by the pilot, the gust disturbance and its
first and second derivatives, applied at three different locations of the
aircraft fuselage ($n_{w}=7$). This vector reads

$\begin{array}[]{rcl}\mathbf{w}_{G}(t)&=&\bigg{[}\mathbf{\delta}_{mc}^{\star}(t),\underbrace{\mathbf{w}_{g}(t),\mathbf{\dot{w}}_{g}(t),\mathbf{\ddot{w}}_{g}(t)}_{\text{aircraft
front}},\underbrace{\mathbf{w}_{g}(t-\tau_{1}^{(i)}),\mathbf{\dot{w}}_{g}(t-\tau_{1}^{(i)}),\mathbf{\ddot{w}}_{g}(t-\tau_{1}^{(i)})}_{\text{aircraft
middle}},\dots\\\
&&\underbrace{\mathbf{w}_{g}(t-\tau_{2}^{(i)}),\mathbf{\dot{w}}_{g}(t-\tau_{2}^{(i)}),\mathbf{\ddot{w}}_{g}(t-\tau_{2}^{(i)})}_{\text{aircraft
rear}}\bigg{]}^{T},\end{array}$ (4)

where the same gust (position, velocity and acceleration) enters in the model
at the front, the middle and the rear, with a delay
$\tau_{1}^{(i)}\in\mathbb{R}_{+}$ and $\tau_{2}^{(i)}\in\mathbb{R}_{+}$,
depending on the aircraft velocity ($i=1,\dots,n_{s}$). The only considered
measurement $\mathbf{y}$ is the aircraft angle of attack at the nose level
($n_{y}=1$). Finally, the performances $\mathbf{z}$ gathers outputs being the
tracking signal error between the load factor without
($\mathbf{n}_{z}^{\star}$) and with ($\mathbf{n}_{z}$) GLA controller
$\mathbf{n}_{z}^{\star}-\mathbf{n}_{z}$ and the loads envelope
$\mathbf{z}_{\text{loads}}(x_{i})$ ($n_{z}=6$). The ${A}^{(i)}$, ${B}^{(i)}$,
${C}^{(i)}$ and ${C}_{d}^{(i)}$ matrices are real and of appropriate
dimensions.

###  2.2 Toward gust load control-oriented models

As signals in eq. 4 are repeated and linked together, it can problematic to
deal with in a control design setup. One solution is to merge them in order to
deal with
$\mathbf{w}(t)=\left[\mathbf{\delta}_{mc}^{\star}(t),\mathbf{w}_{G}(t)\right]^{T}$
instead of eq. 4. This can be obtained by applying the following simple
transformation

$\left[\begin{array}[]{l}\mathbf{\delta}_{mc}^{\star}(t)\\\
\mathbf{w}_{g}(t)\\\ \mathbf{\dot{w}}_{g}(t)\\\ \mathbf{\ddot{w}}_{g}(t)\\\
\mathbf{w}_{g}(t-\tau_{1}^{(i)})\\\ \mathbf{\dot{w}}_{g}(t-\tau_{1}^{(i)})\\\
\mathbf{\ddot{w}}_{g}(t-\tau_{1}^{(i)})\\\ \mathbf{w}_{g}(t-\tau_{2}^{(i)})\\\
\mathbf{\dot{w}}_{g}(t-\tau_{2}^{(i)})\\\
\mathbf{\ddot{w}}_{g}(t-\tau_{2}^{(i)})\\\
\end{array}\right]=\left[\begin{array}[]{cl}1&0\\\ 0&1\\\ 0&s\\\ 0&s^{2}\\\
0&e^{-\tau_{1}^{(i)}s}\\\ 0&se^{-\tau_{1}^{(i)}s}\\\
0&s^{2}e^{-\tau_{1}^{(i)}s}\\\ 0&e^{-\tau_{2}^{(i)}s}\\\
0&se^{-\tau_{2}^{(i)}s}\\\ 0&s^{2}e^{-\tau_{2}^{(i)}s}\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{\delta}_{mc}^{\star}(t)\\\
\mathbf{w}_{G}(t)\end{array}\right],$ (5)

where $s$ denotes the Laplace variable. By merging eq. 3 and eq. 5, one can
now re-construct each model in a realisation form
$\\{\mathcal{S}_{i}\\}_{i=1}^{n_{s}}$ as

$\\{\mathcal{S}_{i}\\}_{i=1}^{n_{s}}:\left\\{\begin{array}[]{rcl}{E}^{(i)}\mathbf{\dot{x}}^{(i)}(t)&=&{A}^{(i)}_{0}\mathbf{x}^{(i)}(t)+{A}_{1}^{(i)}\mathbf{x}^{(i)}(t-\tau_{1}^{(i)})+{A}_{2}^{(i)}\mathbf{x}^{(i)}(t-\tau_{2}^{(i)})\\\
&&+{B}_{g}^{(i)}\left[\begin{array}[]{c}\mathbf{w}(t)\\\
\mathbf{u}(t)\end{array}\right]\\\ \left[\begin{array}[]{c}\mathbf{y}(t)\\\
\mathbf{z}(t)\end{array}\right]&=&{C}^{(i)}_{0}\mathbf{x}^{(i)}(t)+{C}^{(i)}_{1}\mathbf{x}^{(i)}(t-\tau_{m})\end{array}\right.$
(6)

where $\mathbf{x}^{(i)}(t)\in\mathbb{R}^{N_{i}}$,
$\mathbf{u}(t)\in\mathbb{R}^{n_{u}}$, $\mathbf{w}(t)\in\mathbb{R}^{n_{w}}$,
$\mathbf{y}(t)\in\mathbb{R}^{n_{y}}$ and $\mathbf{z}(t)\in\mathbb{R}^{n_{z}}$
are the internal variables, input, exogenous input, output and performance
output signals, respectively. With this new form $n_{u}$, $n_{y}$ and $n_{z}$
are unaffected, while $n_{w}=2$ (instead of 10 in eq. 4) and $N_{i}=M_{i}+6$,
due to the double derivative and delay structure added. The ${E}^{(i)}$,
${A}_{0}^{(i)}$, ${A}_{1}^{(i)}$, ${A}_{2}^{(i)}$, ${B}_{g}^{(i)}$,
${C}_{0}^{(i)}$ and ${C}_{1}^{(i)}$ matrices are real and of appropriate
dimensions.

######  Remark 1 (About the internal delays $\tau_{1}^{(i)}$ and
$\tau_{2}^{(i)}$, and the ${E}^{(i)}$ matrix rank)

In the considered use-case, the fuselage is subdivided in three patches.
Consequently, the gust vertical displacement signal $\mathbf{w}_{g}(t)$ given
in eq. 1, enters the model at $t$, $t+\tau_{1}^{(i)}$ and $t+\tau_{2}^{(i)}$
($0<\tau_{1}^{(i)}<\tau_{2}^{(i)}$). To consider one single gust input signal
(instead of three delayed), internal delay are added in the model. Similarly,
to be able to accurately compute the loads along the wings, the model must
takes as gust input, its vertical displacement as in eq. 1, and its first and
second derivatives as in eq. 4 [33]. Once again, to limit the model number of
input disturbances to one instead of three (position, velocity and
acceleration), the first and second derivatives of eq. 1 are embedded in the
model, leading to a descriptor form as in eq. 6. More specifically, in our
case, the ${E}^{(i)}$ matrix rank is then equal to $N_{i}-6$: two rank loss
per fuselage patch (one for the first and one for the second derivative).

Following $\\{\mathcal{S}_{i}\\}_{i=1}^{n_{s}}$ given in equation eq. 6, the
gust load oriented model is a family of dynamical models which transfer
function $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$, from $\mathbf{u}$ and
$\mathbf{w}$ to $\mathbf{y}$ and $\mathbf{z}$, reads

$\mathbf{H}_{i}(s)=\left({C}^{(i)}_{0}+{C}^{(i)}_{1}e^{-\tau^{(i)}_{m}s}\right)\left(s{E}^{(i)}-{A}^{(i)}_{0}-{A}^{(i)}_{1}e^{\tau^{(i)}_{1}s}-{A}^{(i)}_{2}e^{\tau^{(i)}_{2}s}\right)^{-1}{B}^{(i)}.$
(7)

##  3 Preliminaries in rational interpolation and approximation

The entire proposal relies on a specific use of the model interpolation tools:
namely, the _Loewner_ and the _optimal (frequency-limited) ${\mathcal{H}_{2}}$
frameworks_. The Loewner framework is recalled in its general form, in section
3.1. The ${\mathcal{H}_{2}}$ model order dimension reduction methods is
recalled in section 3.2, section 3.3 and section 3.4. In the following
subsections, indices and exponents in eq. 6 and eq. 7 are dropped for
simplicity.

###  3.1 Reminder of the Loewner framework

The main elements of the Loewner framework are recalled thereafter in the
multi-input multi-output (MIMO) general square case. For a complete
description readers may refer to [5, 24], and to [3] for insight in the
rectangular case. Under mild considerations, the Loewner approach is a data-
driven method aimed at building a rational descriptor LTI dynamical model
$\mathbf{H}^{m}$ of dimension $m$ which interpolates given complex data, here
generated by a model $\mathbf{H}$. Let us be given the left or row data, and
the right or column data:

$\left.\begin{array}[]{c}(\mu_{j},\mathbf{l}_{j}^{H},\mathbf{v}_{j}^{H})\\\
\text{for $\displaystyle
j=1,\dots,m$}\end{array}\right\\}\text{~{}~{}and~{}~{}}\left\\{\begin{array}[]{c}(\lambda_{i},\mathbf{r}_{i},\mathbf{w}_{i})\\\
\text{for $\displaystyle i=1,\dots,m$}\end{array}\right.,$ (8)

where $\mathbf{v}_{j}^{H}=\mathbf{l}_{j}^{H}\mathbf{H}(\mu_{j})$ and
$\mathbf{w}_{i}=\mathbf{H}(\lambda_{i})\mathbf{r}_{i}$, with
$\mathbf{l}_{j}\in\mathbb{C}^{n_{y}\times 1}$,
$\mathbf{r}_{i}\in\mathbb{C}^{n_{u}\times 1}$,
$\mathbf{v}_{j}\in\mathbb{C}^{n_{u}\times 1}$ and
$\mathbf{w}_{i}\in\mathbb{C}^{n_{y}\times 1}$. The $\mathbf{w}_{i}$ and
$\mathbf{v}_{j}$ are the complex measurement data at the points $\mu_{j}$ and
$\lambda_{i}$, along the tangential directions $\mathbf{l}_{j}$ and
$\mathbf{r}_{i}$. In addition, let us assume that $\lambda_{i}\in\mathbb{C}$
and $\mu_{j}\in\mathbb{C}$ represent a set of distinct interpolation points
$\\{z_{k}\\}_{k=1}^{2m}\subset\mathbb{C}$ which is split in two equal subsets
as

$\\{z_{k}\\}_{k=1}^{2m}=\\{\mu_{j}\\}_{j=1}^{m}\cup\\{\lambda_{i}\\}_{i=1}^{m}.$
(9)

The method then consists in building the _Loewner_
${\mathds{L}}\in\mathbb{C}^{m\times m}$ and _shifted Loewner_
${\mathds{L}_{\sigma}}\in\mathbb{C}^{m\times m}$ matrices defined as follows,
for $i=1,\dots,m$ and $j=1,\dots,m$:

$\begin{array}[]{rclcl}[{\mathds{L}}]_{j,i}&=&\dfrac{\mathbf{v}_{j}^{H}\mathbf{r}_{i}-\mathbf{l}_{j}^{H}\mathbf{w}_{i}}{\mu_{j}-\lambda_{i}}&=&\dfrac{\mathbf{l}_{j}^{H}\big{(}\mathbf{H}(\mu_{j})-\mathbf{H}(\lambda_{i})\big{)}\mathbf{r}_{i}}{\mu_{j}-\lambda_{i}}\\\
\,[{\mathds{L}_{\sigma}}]_{j,i}&=&\dfrac{\mu_{j}\mathbf{v}_{j}^{H}\mathbf{r}_{i}-\lambda_{i}\mathbf{l}_{j}^{H}\mathbf{w}_{i}}{\mu_{j}-\lambda_{i}}&=&\dfrac{\mathbf{l}_{j}^{H}\big{(}\mu_{j}\mathbf{H}(\mu_{j})-\lambda_{i}\mathbf{H}(\lambda_{i})\big{)}\mathbf{r}_{i}}{\mu_{j}-\lambda_{i}}\end{array}.$
(10)

Then, the model $\mathbf{H}^{m}$ given by the following descriptor
realisation,

$\mathcal{S}^{m}:\left\\{\begin{array}[]{rcl}{E}^{m}\delta\left\\{\mathbf{x}(\cdot)\right\\}&=&{A}^{m}\mathbf{x}(\cdot)+{B}^{m}\mathbf{u}(\cdot)\\\
\mathbf{y}(\cdot)&=&{C}^{m}\mathbf{x}(\cdot)\end{array}\right.,$ (11)

where ${E}^{m}=-{\mathds{L}}$, ${A}^{m}=-{\mathds{L}_{\sigma}}$,
$[{B}^{m}]_{k}=\mathbf{v}_{k}^{H}$ and $[{C}^{m}]_{k}=\mathbf{w}_{k}$ (for
$k=1,\ldots,m)$, and whose related transfer function

$\mathbf{H}^{m}(\xi)={C}^{m}\left(\mathbf{\xi}{E}^{m}-{A}^{m}\right)^{-1}{B}^{m},$
(12)

interpolates $\mathbf{H}$ at the given driving frequencies and directions
defined in eq. 8, _i.e._ satisfies

$\begin{array}[]{rcl}\mathbf{l}_{j}^{H}\mathbf{H}^{m}(\mu_{j})&=&\mathbf{l}_{j}^{H}\mathbf{H}(\mu_{j})\\\
\mathbf{H}^{m}(\lambda_{i})\mathbf{r}_{i}&=&\mathbf{H}(\lambda_{i})\mathbf{r}_{i}\end{array}.$
(13)

######  Remark 2 (About $\delta\left\\{\cdot\right\\}$ and $\xi$ notations)

In eq. 11, “$(\cdot)$” denotes the considered time-domain variable; this
latter can either be “$(t)$” for continuous-time models ($t\in\mathbb{R}_{+}$)
or “$[q]$” for sampled-time models ($q\in\mathbb{Z}$). Similarly, in eq. 11
“$\delta\left\\{\cdot\right\\}$” stands as the shift operator being either
$\delta\\{\mathbf{x}(t)\\}=\mathbf{\dot{x}}(t)$ in the continuous-time case
and $\delta\\{\mathbf{x}(q)\\}=\mathbf{x}[q+1]$ in the sampled one. With
reference to eq. 12 $\mathbf{\xi}$ stands as the associated complex version
being the Laplace variable $\mathbf{\xi}=s$ in the continuous-time case and
the forward shift $\mathbf{\xi}=z$ in the sampled-time one.

Assuming that the number $2m$ of available data is large enough, then it has
been shown in [24] that a minimal model $\mathbf{H}^{n}$ of dimension $n<m$
that still interpolates the data can be built with a projection of eq. 11
provided that, for $k=1,\ldots,2m$, where $z_{k}$ are as in eq. 9,

$\operatorname*{\mathbf{rank}}\left(z_{k}{\mathds{L}}-{\mathds{L}_{\sigma}}\right)=\operatorname*{\mathbf{rank}}\left([{\mathds{L}},{\mathds{L}_{\sigma}}]\right)=\operatorname*{\mathbf{rank}}\left([\mathbb{L}^{H},{\mathds{L}_{\sigma}}^{H}]^{H}\right)=n.$
(14)

In that case, let us denote by $Y\in\mathbb{C}^{m\times n}$ the matrix
containing the first $n$ left singular vectors of
$[{\mathds{L}},{\mathds{L}_{\sigma}}]$ and $X\in\mathbb{C}^{m\times n}$ the
matrix containing the first $n$ right singular vectors of
$[{\mathds{L}}^{H},{\mathds{L}_{\sigma}}^{H}]^{H}$. Then,

${E}^{n}=Y^{H}{E}^{m}X,\,{A}^{n}=Y^{H}{A}^{m}X,\,{B}^{n}=Y^{H}{B}^{m},\,{C}^{n}={C}^{m}X,$
(15)

is a realisation of the model $\mathbf{H}^{n}$ given as,

$\mathbf{H}^{n}(\xi)={C}^{n}\left(\mathbf{\xi}{E}^{n}-{A}^{n}\right)^{-1}{B}^{n},$
(16)

with the same structure as eq. 12, embedding a _minimal McMillan degree_ equal
to $\operatorname*{\mathbf{rank}}\left({\mathds{L}}\right)$. The quadruple
$\mathcal{S}^{n}:({E}^{n},{A}^{n},{B}^{n},{C}^{n})$ is a descriptor
realisation of $\mathbf{H}^{n}$. Note that if $n$ in eq. 14 is greater than
$\operatorname*{\mathbf{rank}}({\mathds{L}})$ then $\mathbf{H}^{n}$ can either
have a direct-feedthrough or a polynomial part. Reader may note that the
number $n$ of singular vectors composing $Y$ and $X$ used to project the
system $\mathbf{H}^{n}$ in eq. 15 may be decreased even further than $n$ at
the cost of approximate interpolation of the data. This allows for a trade-off
between complexity of the resulting model and accuracy of the interpolation.
Reader may note that Loewner extensions to some nonlinear systems exist,
_e.g._ bilinear [6].

In this work we will always consider exact interpolation while the model
reduction step, recalled in section 3.2, will be performed by (frequency
limited) ${\mathcal{H}_{2}}$-oriented interpolatory methods instead. These
latter are detailed in a rather large framework in the following section 3.3
and section 3.4.

######  Remark 3 (Assumptions for model order reduction)

In what follows, we assume that $\mathbf{H}^{n}$ is stable, strictly proper
and has semi-simple poles. These assumptions can be theoretically removed but
at the cost of a more complicated developments, not appropriate in our gust
load setting. Indeed, any polynomial or constant terms, as highlighted in eq.
7, will be kept for simulation but discarded for control design, with a low if
not null impact on the process.

###  3.2 Model dimension reduction

Given the finite $n$-th order function $\mathbf{H}^{n}$ equipped with a
realisation defined by the quadruple
$\mathcal{S}^{n}:({E}^{n},{A}^{n},{B}^{n},{C}^{n})$ given in eq. 15, the model
reduction goal consists in constructing reduced $r$-th order ($r\ll n$) model

$\mathcal{\hat{S}}:\left\\{\begin{array}[]{rcl}{\hat{E}}\delta\left\\{\mathbf{\hat{x}}(\cdot)\right\\}&=&{\hat{A}}\mathbf{\hat{x}}(\cdot)+{\hat{B}}\mathbf{u}(\cdot)\\\
\mathbf{\hat{y}}(\cdot)&=&{\hat{C}}\mathbf{\hat{x}}(\cdot)\end{array}\right.,$
(17)

where, $\mathbf{\hat{x}}(\cdot)\in\mathbb{R}^{r}$ is the reduced internal
variables and $\mathbf{\hat{y}}(\cdot)\in\mathbb{R}^{n_{y}}$ is the
approximated output, and where

${\hat{E}}\in\mathbb{R}^{r\times r},{\hat{A}}\in\mathbb{R}^{r\times
r},{\hat{B}}\in\mathbb{R}^{r\times n_{u}}\text{ and
}{\hat{C}}\in\mathbb{R}^{n_{y}\times r}$ (18)

are constant matrices, such that the input output behaviour of
$\mathbf{\hat{H}}$ defined as

$\mathbf{\hat{H}}(\xi)={\hat{C}}\left(\xi{\hat{E}}-{\hat{A}}\right)^{-1}{\hat{B}},$
(19)

is similar to $\mathbf{H}^{n}$, _i.e._ for the same input $\mathbf{u}$, the
mismatch $\mathbf{y}-\mathbf{\hat{y}}$ is small in some sense.

###  3.3 Rational ${\mathcal{H}_{2}}$ approximation and reduction by
interpolation

One way to find $\mathbf{\hat{H}}$ as in eq. 19 is to solve the so-called
${\mathcal{H}_{2}}$ model reduction problem222Optimising through the
${\mathcal{H}_{2}}$ norm is relevant as a it provides an upper bound on the
worst case time-domain error in response to an input signal [18]., which aims
at solving the following problem:

$\mathbf{\hat{H}}=\arg\min_{\mathbf{G}\in{\mathcal{H}_{2}}}||\mathbf{H}^{n}-\mathbf{G}||_{{\mathcal{H}_{2}}}.$
(20)

The most common approach to obtain the solution of eq. 20 is to work with the
_first order necessary optimality conditions_ which were developed in a series
of papers, see _e.g._ [10, 37] for theoretical insight (see also [9] for
existence and general properties). The _interpolation-based_ approach
revisited by a sequence of contributions [15, 38, 18] and resulting in an
interpolatory problem is used here. In the case of semi-simple poles only (the
case of higher order poles is presented in [39]), if $\mathbf{\hat{H}}$,
equipped with a realisation
$\mathcal{\hat{S}}:({\hat{E}},{\hat{A}},{\hat{B}},{\hat{C}})$ is a solution of
the ${\mathcal{H}_{2}}$ approximation problem eq. 20, then

$\begin{array}[]{rcl}\mathbf{H}^{n}(\hat{\kappa}_{l})\hat{\mathbf{b}}_{l}&=&\mathbf{\hat{H}}(\hat{\kappa}_{l})\hat{\mathbf{b}}_{l}\\\
\hat{\mathbf{c}}_{l}^{H}\mathbf{H}^{n}(\hat{\kappa}_{l})&=&\hat{\mathbf{c}}_{l}^{H}\mathbf{\hat{H}}(\hat{\kappa}_{l})\\\
\hat{\mathbf{c}}_{i}^{H}\left.\dfrac{d\mathbf{H}^{n}}{d\xi}\right\rvert_{\xi=\hat{\kappa}_{l}}\hat{\mathbf{b}}_{l}&=&\hat{\mathbf{c}}_{l}^{H}\left.\dfrac{d\mathbf{\hat{H}}}{d\xi}\right\rvert_{\xi=\hat{\kappa}_{l}}\hat{\mathbf{b}}_{l}\end{array},$
(21)

where $[\hat{\mathbf{b}}_{1},\dots,\hat{\mathbf{b}}_{r}]^{H}=R{\hat{B}}$ and
$[\hat{\mathbf{c}}_{1},\dots,\hat{\mathbf{c}}_{r}]={\hat{C}}L$ and where
$L\in\mathbb{C}^{r\times r}$ and $R\in\mathbb{C}^{r\times r}$ are the left and
right eigenvectors associated to $\hat{\lambda}_{l}$, the eigenvalues of the
$({\hat{E}},{\hat{A}})$ pair. Then, in the continuous-time case,
$\hat{\kappa}_{l}=-\hat{\lambda}_{l}$ while in the discrete-time one,
$\hat{\kappa}_{l}=1/\hat{\lambda}_{l}$ ($l=1,\dots,r$) [18, 11].

In [18], authors derive a Petrov-Galerkin like projection-based approach [40],
accompanied with a fixed point algorithm, celebrated as MIMO Iterative
Rational Krylov Algorithm (MIMO IRKA), allowing reaching the optimality
conditions eq. 21. The solution results in a procedure with a fairly cheap
computational cost, embedding linear equation resolutions for which efficient
methods are available [35, 22, 36, 19]. Details are skipped here but readers
can also refer to the book [4] or monograph [29, Chap. 2] for practical
details and numerical applications.

###  3.4 Rational frequency-limited ${\mathcal{H}_{2}}$ approximation and
reduction by interpolation

Similarly to the above ${\mathcal{H}_{2}}$ model approximation problem, it
might be interesting approximating over a frequency limited range
$\Omega=[0,\omega]$ ($\omega\in\mathbb{R}_{+}$) only. This is the purpose of
the _frequency-limited ${\mathcal{H}_{2}}$_, shortly
${\mathcal{H}_{2,\Omega}}$, approximation problem defined as [44, 27]:

$\mathbf{\hat{H}}=\arg\min_{\mathbf{G}\in{\mathcal{H}_{\infty}}}||\mathbf{H}^{n}-\mathbf{G}||_{{\mathcal{H}_{2,\Omega}}}.$
(22)

Interestingly, as in the seminal works of [18] treating the
${\mathcal{H}_{2}}$ objective, [41, Chap. 8] shown that this problem can also
be recast as an interpolatory one.

However, unlike the ${\mathcal{H}_{2}}$ interpolation conditions, the
frequency-limited ones do not involve directly the transfer functions
$\mathbf{H}$ and $\mathbf{\hat{H}}$ but irrational functions
$\mathbf{T}^{n}_{\omega}(\mathbf{H})$ and
$\mathbf{\hat{T}}_{\omega}(\mathbf{\hat{H}})$, parametrised by $\mathbf{H}$
and $\mathbf{\hat{H}}$. Here again, these functions should match at images of
the poles of the reduced-order model (see details in [41, Chap. 8]).

As it, these interpolatory conditions are difficult to practically exploit.
Indeed, no Krylov-like subspace have been clearly identified yet.
Consequently, a MIMO IRKA like procedure as in [18] is not straightforward to
develop333Note that in [28], a similar procedure achieving first order
${\mathcal{H}_{2}}$ optimality with interpolatory conditions with an input-
output delayed LTI model has been developed.. So far, this
${\mathcal{H}_{2,\Omega}}$ problem has been attacked using the DARPO procedure
[41, Chap. 8], a descent algorithm, and FL-ISTIA procedure [43], involving
frequency-limited gramian, taking advantage of a subset of the interpolation
conditions eq. 21. In the remaining of the paper, we will consider the FL-
ISTIA procedure developed in [43].

##  4 Interpolatory-driven aero-servoelastic aircraft gust load oriented
reduced modelling

Back to our GLA problem, we will now invoke the interpolatory results of
section 3 to construct an appropriate control-oriented model for the GLA
control design, simulation and analysis. The following section 4.2 and section
4.3 successively describe the exact continuous-time rational approximation and
frequency-limited model order reduction of the aeroelastic business jet
aircraft models. The reason for these successive sections are justified in
section 4.1. Numerical results are presented in section 4.4.

###  4.1 Models limitations for control purpose

This BizJet use-case is described by models defined as eq. 6 in state-space
form $\\{\mathcal{S}_{i}\\}_{i=1}^{n_{s}}$ and eq. 7 in its transfer form
$\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$. These relations are not appropriate to
manipulate in a control design objective for the following reasons: _(i)_ the
internal vector dimension of each model $\mathcal{S}^{(i)}$ is large (_i.e._
$N_{i}\approx 300$), leading to medium-scale matrices and a computational
burden inappropriate with the control design; _(ii)_ the presence of internal
delays renders the model irrational and results in an infinite number of
eigenvalues, for which dedicated tools exist but are complex to manipulate in
an industrial context; _(iii)_ the ${E}^{(i)}$ matrix being rank deflective,
the computational tools are not always appropriate and transferable in an
applicative context.

###  4.2 BizJet aircraft rational and polynomial approximation

The first step of the process consists in replacing each infinite444These
models, embedding a state-space of order $N_{i}$ are of infinite dimension due
to the internal delays, leading to infinite number of eigenvalues. dimensional
model $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$ eq. 7 by a finite order rational
$n_{i}$-th order one, $\\{\mathbf{H}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$ ($i=1,\dots
n_{s}$). By considering each $n_{s}$ continuous-time model of the considered
GLA use-case described in section 1, one can apply the rational approximation
by interpolation using the method recalled in section 3.1. Then one obtains a
set of $n_{s}$ continuous-time rational approximated models
$\\{\mathbf{H}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$ as eq. 16 equipped with
realisation $\\{\mathcal{S}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$ with matrices as eq.
15. Note that each of the $n_{s}$ models may have a different dimension
$n_{i}$; this latter being automatically computed by the rank condition given
in eq. 14. Here, as the original models eq. 7 are given in continuous-time, it
is convenient to select the interpolation points along the imaginary axis. In
our case, for all configurations, the interpolating points eq. 9 have been
selected as follows :

$\left\\{z_{k}\right\\}_{k=1}^{2m}=\underbrace{\left\\{\imath\omega_{j},-\imath\omega_{j}\right\\}_{j=1}^{m/2}}_{\left\\{\mu_{j}\right\\}_{j=1}^{m}}\cup\underbrace{\left\\{\imath\nu_{i},-\imath\nu_{i}\right\\}_{i=1}^{m/2}}_{\left\\{\lambda_{i}\right\\}_{i=1}^{m}},$
(23)

where $\omega_{j},\nu_{i}\in\mathbb{R}_{+}$ are the pulsations at which one
evaluates each transfer $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$, with $m=1000$.
In our application $\omega_{j}$ and $\nu_{i}$ are selected to be
logarithmically spaced from $10^{-2}$ to $10^{3}$. This choice allows to focus
on the frequency range of interest. Indeed, in the case of irrational models,
the method is efficient for interpolation but not necessary for extrapolation.

In addition, as the stability of the obtained rational model
$\\{\mathbf{H}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$ is not guaranteed by the Loewner
interpolatory procedure, a _post stabilisation_ is performed using the
procedure presented in [21]. This latter consists in projecting the rational
models $\\{\mathbf{H}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$ onto its closest stable
subset, here using the ${\mathcal{H}_{\infty}}$-norm, leading to a stable
model of the same dimension. Mathematically, given a realisation $\mathcal{S}$
associated to $\mathbf{H}\in{\mathcal{L}_{\infty}}$, one aims at finding
$P_{\infty}(\mathbf{H})\in{\mathcal{H}_{\infty}}$ such that,

$P_{\infty}(\mathbf{H})=\arg\inf_{\mathbf{G}\in{\mathcal{H}_{\infty}}}||\mathbf{H}-\mathbf{G}||_{\mathcal{H}_{\infty}}.$
(24)

Technical details and assumption can be found in [21]. Each model now share a
rational structure as the one given in eq. 17-eq. 19, and is now stable and of
finite order $n_{i}$, for $i=1,\dots,n_{s}$. Note that here, the stability
enforcement is performed since one knows that original models are structurally
and physically stable. At this point, the delay difficulty is removed and
traded with rational order model $n_{i}$ automatically determined by eq. 14.

######  Remark 4 (About a Padé delay approximation)

It can be appealing to replace the delays with a Padé approximation, which
preserves the gain but modifies the phase. Indeed, while this is classically
performed in many applications, it is, to the authors experience, not the best
way to deal with delays, especially in flexible structure models where phase
vary a lot. In addition, the use of Padé will drastically increases the model
internal vector. Although it may not be a problem, Padé usually results in
poor accuracy, specifically in flexible structure, affecting the phase.
Instead, we prefer rational multi-point interpolation while Padé may be viewed
as an interpolation at higher order, at zero.

###  4.3 BizJet aircraft control-oriented model reduction

As a second step, and rooted on $\\{\mathbf{H}_{i}^{n_{i}}\\}_{i=1}^{n_{s}}$,
the obtained continuous-time rational models, one now invokes again the
_interpolatory framework_ for a dimension reduction. Now we follow the
approach detailed in section 3.4 to reduce the model over a frequency-limited
range. For the considered application, and in view of control design, it is
important to reduce as much as possible the model while staying
representative. Indeed, as the controller designed later is a solution of an
optimisation problem, the simpler the model is, the more accurate and faster
the optimisation is.

In addition, as the GLA function should act in a given frequency range without
altering the flight control laws, it is relevant to reduce the model over a
frequency limited range $[0,\omega]$, only (where $\omega\in\mathbb{R}_{+}$).
FL-ISTIA process [43] being also a fixed point procedure, its initialisation
is done by selecting starting interpolating points as follows:

$\begin{array}[]{rcll}\\{\hat{\kappa}^{(0)}\\}_{k=1}^{r}&=&\\{\imath\omega_{k},-\imath\omega_{k}\\}_{k=1}^{r/2}&\text{
if $r$ is even}\\\
\\{\hat{\kappa}^{(0)}\\}_{k=1}^{r}&=&\\{\imath\omega_{k},-\imath\omega_{k}\\}_{k=1}^{\left\lfloor{r/2}\right\rfloor}\cup\omega_{0}&\text{
if $r$ is odd}\end{array},$ (25)

where $0<\omega_{k}<\omega$ and $\omega_{0}<\omega$. The interpolatory
conditions presented in eq. 21 can also be partially ensured.

###  4.4 Application to the BizJet models

As an illustration, fig. 2 compares the frequency and phase responses from the
gust input to a wing root bending moment of the original irrational model
$\mathbf{H}_{7}$ with its rational approximations $\mathbf{H}_{7}^{n_{7}}$ and
frequency-limited reduced model $\mathbf{\hat{H}}_{7}$ for one single
configuration point. Sigma plot and mismatch errors are also reported in fig.
3. All other cases are similarly obtained.

Figure 2: Gain (top) and phase (bottom) frequency responses of one transfer
from gust to a wing root bending moment. The original model $\mathbf{H}_{7}$
(solid blue), the approximated and approximated plus projected stable models
$\mathbf{H}^{n_{7}}_{7}$ (dashed magenta dash dotted black, respectively) and
the reduced order model $\mathbf{\hat{H}}_{7}$ with dimension $r=30$. Upper
bound of the reduction frequency band is materialised by the vertical green
dashed line.

Figure 3: Gain (top) and gain error (bottom) of the singular frequency
responses of the different models.

With reference to fig. 2 and fig. 3, the following comments can be done.
First, the rational approximation, satisfying the tangential interpolatory
conditions eq. 13, allows to accurately capture the descriptor delayed model
eq. 7 gain, phase and sigma, with a simple rational model (here for the 7-th
use-case, $n_{7}=211$) of the form eq. 16. Note that the projection [21]
performed afterward to force the model to be stable let the matching accurate
enough. Second, the frequency-limited order reduction still preserves a good
matching in the considered frequency range, with a rational model as eq. 19 of
dimension $r=30$ only. Analysing the bottom frame of fig. 3, the
$\mathbf{H}_{7}-\mathbf{H}_{7}^{n7}$ (unstable) mismatch is close to machine
precision up to the upper frequency limit $\omega$; this is caused by the
tangential directions which may produce less accurate transfers. When
projecting to obtain a stable model, the $\mathbf{H}_{7}-\mathbf{H}_{7}^{n7}$
(stable) mismatch results in a loss of accuracy, traded with the stability
property. The reduced model $\mathbf{\hat{H}}_{7}$ clearly is accurate below
the upper frequency bound, and can be considered as simple but accurate enough
for GLA control design.

##  5 Interpolatory-driven sampled-time feedback gust load alleviation
function computation

Given the continuous-time reduced order models
$\\{\mathbf{\hat{H}}_{i}\\}_{i=1}^{n_{s}}$, we are now ready to design the GLA
control function. As this control law has to be implemented on a on-board
computer sampled with a fixed step time $h\in\mathbb{R}_{+}$, a sampled-time
(or discrete-time) controller is sought. In the rest of this section, we first
present the context of discrete-time control design in section 5.1. Then,
based on the continuous-time reduced models
$\\{\mathbf{\hat{H}}_{i}\\}_{i=1}^{n_{s}}$, the continuous-time synthesis of
$\mathbf{K}$, in the continuous-time ${\mathcal{H}_{\infty}}$-norm, is briefly
exposed in section 5.2. Then, section 5.3 exposes the discretisation problem,
which is then solved in the _new interpolatory framework_ in section 5.4,
leading to algorithm 1.

###  5.1 Preliminary words on discrete and continuous-time control design

In dynamical systems and control theory, the continuous-time and discrete-time
domains coexist. While most of the tools available in one domain have a
counterpart in the other, it is not unusual that a specific application
requires to switch from one domain to the other. In particular, in the
considered GLA control-oriented application, engineers start from a
continuous-time physical model set $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$,
simplified by a reduced continuous-time one
$\\{\mathbf{\hat{H}}_{i}\\}_{i=1}^{n_{s}}$. The latter is then used for the
design of the GLA control law. To this aim, as illustrated schematically in
fig. 4, three approaches can be conducted: _(i)_ one may discretise the analog
plants $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$ and then design a discrete
control-law $\mathbf{K}_{d}$ with the same sampling period $h$, _(ii)_
conversely, one may design an analog control-law $\mathbf{K}$ and then
discretise it, _(iii)_ or, using dedicated techniques from sampled-data
systems theory (see _e.g._ [13, Chap. 12-13]), one may directly synthesise
$\mathbf{K}_{d}$ from $\\{\mathbf{H}_{i}\\}_{i=1}^{n_{s}}$.

$\mathbf{H}$ Continous-time plant $\mathbf{H}_{d}$ Sampled-time plant Plant
discretisation (zoh) $\mathbf{K}_{d}$ Sampled-time controller Discrete
synthesis $\mathbf{K}$ Continuous-time controller Controller discretisation
(zoh, tustin, etc.) Analog synthesis Sampled-data synthesis Figure 4:
Different paths to design a discrete control-law $\mathbf{K}_{d}$ from an
analog plant $\mathbf{H}$.

While the direct nature of the latter method _(iii)_ is appealing, it requires
dedicated theoretical and numerical tools that are not as widespread as usual
ones, especially in the industry where this approach would require their whole
control design and analysis process to be rethought. For these reasons, one
focuses here on the indirect methods that remain of practical interest, and
more specifically on the continuous controller, followed by a discretisation.
Performing in the other ways, _i.e._ discretising the model and then
synthesising the controller would lead to an inaccurate model response due to
the low frequency sampling imposed by the on-board computer.

For both indirect approaches, a discretisation step is required. This may have
a detrimental impact with respect to the expected dynamical behaviour,
especially if hard computational constraints on the sampling period $h$ are
imposed by the technology, which may happen with critical systems. In that
context, the availability of an efficient discretisation method is of
particular interest. Here, an approach based on the _Loewner interpolatory
framework_ [24], recalled in section 3.1, is used. It offers an interesting
alternative to usual discretisation processes as it enables to reach a better
frequency and time-domain matching with the continuous-time description.
Before detailing this approach, for sake of completeness, section 5.2 briefly
presents the continuous-time GLA controller design.

###  5.2 Preliminary continuous-time GLA controller design

In the considered application, as the objective is to attenuate a worst case
amplification, the ${\mathcal{H}_{\infty}}$-norm minimisation framework is
clearly tailored. This optimisation method is then used, on the basis of the
reduced continuous-time model representation. Note that this framework, and
linear ones, are largely used and mastered in aeronautical applications [26,
30, 1], and is preferred here also. Being given the reduced continuous-time
models $\\{\mathbf{\hat{H}}_{i}\\}_{i=1}^{n_{s}}$, the continuous-time
controller $\mathbf{K}$, mapping outputs $\mathbf{y}$ to control inputs
$\mathbf{u}$, as interconnected on fig. 1, is obtained by solving the
following ${\mathcal{H}_{\infty}}$-norm problem, for $i=1,\dots n_{s}$,

${\mathbf{K}}=\arg\min_{\tilde{\mathbf{K}}\in\mathcal{K}\subseteq{\mathcal{H}_{2}}}||\mathcal{F}_{l}(\mathbf{G}(\mathbf{\hat{H}}_{i}),\tilde{\mathbf{K}})||_{\mathcal{H}_{\infty}},$
(26)

where $\mathbf{G}(\mathbf{\hat{H}}_{i})\in{\mathcal{H}_{\infty}}$ is the so-
called generalised plant involving the input and output weight functions (not
detailed here as out of the scope) and the reduced order models. Moreover
$\mathcal{F}_{l}(\cdot,\cdot)$ denotes the lower fractional operator. The
control objective thus consists in finding the locally-optimal continuous-time
controller $\mathbf{K}$, such that eq. 26 is solved, _i.e._ such that the
${\mathcal{H}_{\infty}}$-norm of the
$\mathcal{F}_{l}(\mathbf{G}(\mathbf{\hat{H}}_{i}),{\mathbf{K}})$ loop, mapping
exogenous inputs $\mathbf{w}$ to performance output $\mathbf{z}$ (here being
the loads and some tracking objectives), is minimised. In the considered case,
one seeks this controller such that it belongs to a subset
$\mathcal{K}\subseteq{\mathcal{H}_{2}}$, reducing the controller search space
to reduced order structures without direct feed-through, avoiding issues for
industrial implementation. The problem is solved by the routine developed by
the authors of [7] and results in a stable moderate $n_{c}$-th order LTI
continuous-time controller $\mathbf{K}$.

###  5.3 Controller discretisation problem statement and discretisation error
derivation

Being given the linear continuous-time controller $\mathbf{K}$, the objective
here is to determine, for a fixed time-step $h\in\mathbb{R}_{+}$, a discrete-
time controller represented by the recurrence state-space equation,

$\mathcal{C}_{d}:\left\\{\begin{array}[]{rcl}\mathbf{x}_{d}[q+1]&=&{A}_{d}\mathbf{x}_{d}[q]+{B}_{d}\mathbf{y}_{d}[q]\\\
\mathbf{u}_{d}[q]&=&{C}_{d}\mathbf{x}_{d}[q]+{D}_{d}\mathbf{y}_{d}[q]\end{array}\right.$
(27)

where $\mathbf{x}_{d}[q]\in\mathbb{R}^{n_{d}}$,
$\mathbf{u}_{d}[q]\in\mathbb{R}^{n_{u}}$,
$\mathbf{y}_{d}[q]\in\mathbb{R}^{n_{y}}$ are the internal variables, sampled
control and measurements signals respectively. One seeks for eq. 27 such that
_(i)_ its associated transfer function
$\mathbf{K}_{d}(z)={C}_{d}(zI_{n_{d}}-{A}_{d})^{-1}{B}_{d}+{D}_{d}$ is stable,
_i.e._ $\mathbf{K}_{d}\in\mathcal{H}_{\infty}(\overline{\mathcal{D}})$, and
_(ii)_ the input-output behaviour of $\mathbf{K}$ is well reproduced by
$\mathbf{K}_{d}$ up to the Nyquist frequency. While discretisation methods
such as the bilinear (Tustin) are clearly able to build a discrete-time model
satisfying _(i)_ the input-output behaviour may be quite far from the original
one when $h$ is too large. Here, a method is proposed to build such a model by
first using the _Loewner framework_ to interpolate a specific set of frequency
data and then, similarly to section 4.2, projecting the resulting model onto
the stable subspace $\mathcal{H}_{\infty}(\overline{\mathcal{D}})$ to enforce
sampled-time controller stability.

####  5.3.1 Measure of the discretisation error

Quantifying the error induced by the discretisation of $\mathbf{K}$ is not
trivial due to the incompatible nature with $\mathbf{K}_{d}$ that prevents
from interconnecting the two systems. Indeed, one cannot just consider
$\mathbf{K}-\mathbf{K}_{d}$. Digital to analog converters are required. For a
sampling time $h\in\mathbb{R}_{+}$, the latter are modelled here as the ideal
sampler $S$ and holder $H$. Then, the error system is given by the
interconnection of fig. 5 where discrete-time signals are represented by
dashed lines. Such an interconnection of continuous and discrete time models
is called a _Sampled-data System_ (SD).

The controller $\tilde{\mathbf{K}}=S\mathbf{K}_{d}H$ is no longer LTI but
$h$-periodic and consequently, so is the error $\mathbf{K}-\tilde{\mathbf{K}}$
for which usual system norms cannot directly be applied. This problem has been
addressed in the literature for direct SD synthesis (see _e.g._ [8, 13] and
references therein). The recurring idea is to use continuous lifting to
transform a periodic system into a discrete-time LTI model with infinite input
and output spaces on which equivalent ${\mathcal{H}_{2}}$ or
${\mathcal{H}_{\infty}}$ norms can be defined.

While this framework appears well suited to evaluate
$\mathbf{K}-\tilde{\mathbf{K}}$, the infinite dimensionality of the lifted
model makes it quite technical. That is why here, a more straightforward
approach based on the _frequency characterisation of the discretisation error_
formulated in [13, section 3.5] is considered instead. For that purpose,
models for the ideal sampler and holder are recalled in section 5.3.2 and the
frequency error is presented in section 5.3.3555Note that in [13], the authors
use the $\lambda$-transform while here, one considers the $z$-transform, the
sign of the exponential in the Fourier transforms of discrete signal is thus
modified..

$\mathbf{K}$ $\mathbf{K}_{d}$ $S$ $H$
$\mathbf{u}$$\mathbf{e}$$-$$+$$\mathbf{y}_{s}$$\mathbf{y}$

Figure 5: Interconnection for the measurement of the discretisation error
between $G$ and $G_{d}$.

####  5.3.2 Ideal sampler and holder

Let us consider a continuous-time signal $v(t)$ and the sampling period $h$.
The ideal sampler $S$ transforms $v(t)$ into a discrete sequence $v_{d}[q]$
such that,

$v_{d}[q]=v(qh),\,q\in\mathbb{Z}.$ (28)

As shown in [13, Lemma 3.3.1], the Fourier transforms
$\overline{v}_{d}=\mathcal{F}(v_{d})$ and $\overline{v}=\mathcal{F}(v)$ are
linked as follows,

$\overline{v}_{d}(e^{\imath\omega
h})=\frac{1}{h}\sum_{q\in\mathbb{Z}}\overline{v}(\imath\omega+\imath
q\omega_{s}),$ (29)

where $\omega_{s}=2\pi/h$ is the sampling frequency. Equation eq. 29
highlights the frequency aliasing phenomena since all the multiples of the
sampling frequency are indistinguishable in the output. Note that the sampling
operator $S$ is not bounded for any signal in $\mathcal{L}_{2}(\mathbb{R})$ as
shown in [12]. To be bounded, the input signal $v$ must be restricted to the
class of band-limited $\mathcal{L}_{2}(\mathbb{R})$ signals or it must be
filtered by a finite-dimensional stable and strictly causal system, _i.e._ an
anti-aliasing filter.

Similarly, the holder $H$ transforms a sequence $v_{d}[q]$ into a continuous-
time signal $v(t)$ such that,

$v(t)=v_{d}[q],\quad qh\leq t<(q+1)h.$ (30)

The impulse response of the holder can be defined as the difference of two
unit steps delayed by $h$. Let

$R(s)=\frac{1-e^{-sh}}{sh}$ (31)

be the associated transfer function. Then, as shown in [13, Lemma 3.3.2], the
Fourier transforms of $v$ and $v_{d}$ are linked as follows,

$\overline{v}(\imath\omega)=hR(\imath\omega)\overline{v}_{d}(e^{\imath\omega
h}).$ (32)

Back to fig. 5, coupling equations (29) and (32) enables to express the
frequency-domain relationship between $\mathbf{u}$ and $\mathbf{y}_{s}$ (for
$\left\lvert\omega\right\rvert<\left\lvert\omega_{s}\right\rvert$),

$\overline{\mathbf{y}}_{s}(\imath\omega)=R(\imath\omega)\mathbf{K}_{d}(e^{\imath\omega
h})\sum_{k\in\mathbb{Z}}\overline{\mathbf{u}}(\imath\omega+\imath
q\omega_{s}).$ (33)

####  5.3.3 A frequency domain error

Using equation (33), one can express the frequency-domain relationship between
$\mathbf{u}$ and the discretisation error
$\mathbf{e}=\mathbf{y}-\mathbf{y}_{s}$ from Figure 5 as

$\overline{\mathbf{e}}(\imath\omega)=\mathbf{K}(\imath\omega)\overline{\mathbf{u}}(\imath\omega)-R(\imath\omega)\mathbf{K}_{d}(e^{\imath\omega
h})\sum_{k\in\mathbb{Z}}\overline{\mathbf{u}}(\imath\omega+\imath
q\omega_{s}).$ (34)

Assuming that $\mathbf{u}$ is band-limited, _i.e._ that
$\overline{\mathbf{u}}(\imath\omega)=0$ for
$|\omega|>\omega_{s}/2=\omega_{N}$, then (34) becomes

$\overline{\mathbf{e}}(\imath\omega)=\left(\mathbf{K}(\imath\omega)-R(\imath\omega)\mathbf{K}_{d}(e^{\imath\omega
h})\right)\overline{\mathbf{u}}(\imath\omega).$ (35)

This readily suggests to consider

$e_{\infty}(\mathbf{K},\mathbf{K}_{d})=\max_{\omega<\omega_{N}}\left|\mathbf{K}(\imath\omega)-R(\imath\omega)\mathbf{K}_{d}(e^{\imath\omega
h})\right|,$ (36)

to quantify the discretisation error. This frequency-domain characterisation
of the error inspired the discretisation process presented in the next
section. For the MIMO case, the absolute value in (36) may be replaced by the
$2$-norm.

###  5.4 Application of the Loewner framework for discretisation

Now we are ready to describe the _interpolatory driven discretisation process_
allowing to get $\mathbf{K}_{d}$ from the continuous-time controller
$\mathbf{K}$.

####  5.4.1 Principle

Being given the considerations of section 5.3, let us consider the transfer
function of the controller $\mathbf{K}\in{\mathcal{H}_{\infty}}$ to be
discretised at the sampling time $h$. Building a $n_{c}$-th order discrete-
time model $\mathbf{K}_{d}$ that matches the input-output behaviour of
$\mathbf{K}$, the frequency-domain characterisation of the discretisation
error eq. 36 suggests that $\mathbf{K}_{d}$ must be such that

$R(\imath\omega)\mathbf{K}_{d}(e^{\imath\omega h})=\mathbf{K}(\imath\omega),$
(37)

for $|\omega|<\omega_{N}$. Equation eq. 37 represents an infinite number of
interpolation conditions that may be approximated by sampling the interval
$[0,\omega_{N}]$ such that for $k=1,\ldots,2m$,

$R(\imath\omega_{k})\mathbf{K}_{d}(e^{\imath\omega_{k}h})=\mathbf{K}(\imath\omega_{k}).$
(38)

Such a model can be built by applying the _Loewner interpolation_ framework
recalled in section 3.1 to the following set of frequency data,

$\left\\{e^{\imath\omega_{k}h},R(\imath\omega_{k})^{-1}\mathbf{K}(\imath\omega_{k})\right\\}_{k=1}^{2m},$
(39)

where $e^{\imath\omega_{k}h}$ are the new interpolating points and
$R(\imath\omega_{k})^{-1}\mathbf{K}(\imath\omega_{k})$ the new function to
evaluate. With reference to eq. 8 and eq. 9 $e^{\imath\omega_{k}h}$ are the
$z_{k}$. Then by splitting them in the same ways,
$R(\imath\omega_{k})^{-1}\mathbf{K}(\imath\omega_{k})$ represents the couple
$\mathbf{v}_{j}^{H}$ and $\mathbf{w}_{i}$. Following the notations of section
3.1, one obtains the $m$-th order $\mathbf{K}_{d}^{m}$ controller that can be
lowered without loss of interpolatory accuracy to $\mathbf{K}_{d}^{n}$
($n<m$). Obviously, should the order $n_{c}$ be lower than the McMillan order
$n$ of the exact interpolating model given by the Loewner approach, then the
interpolation conditions (38) will not be perfectly satisfied, but the
controller will keep its original dimension.

####  5.4.2 About the stability of the discretised model

The Loewner framework does not ensure the stability of the resulting
interpolaed model, here controller. Therefore, even if
$\mathbf{K}\in{\mathcal{H}_{2}}$ (or ${\mathcal{H}_{\infty}}$), the transfer
function $\mathbf{K}_{d}$ obtained via the process described above may not lie
in $\mathcal{H}_{\infty}(\overline{\mathcal{D}})$ which is a major drawback in
comparison to stability preserving discretisation schemes like the bilinear
(Tustin) of backward ones.

To overcome this issue, one can apply the same process as in [17] and as
performed in section 4.2 (in the continuous-time). It consists in projecting
the unstable model $\mathbf{K}_{d}$ onto
$\mathcal{H}_{\infty}(\overline{\mathcal{D}})$ so that the
$\mathcal{L}_{\infty}$ norm between $\mathbf{K}_{d}$ and its projection is
minimised. This is the so-called Nehari problem for which solutions have been
given in the continuous-time domain [16] and in the discrete-time domain, see
_e.g._ [23].

Let us denote by $P_{\infty}$ this projection operator so that if
$\mathbf{K}_{d}\in\mathcal{L}_{\infty}(\partial\mathcal{D})$ then
$P_{\infty}(\mathbf{K}_{d})\in\mathcal{H}_{\infty}(\overline{\mathcal{D}})$
minimises

$\|\mathbf{K}_{d}-P_{\infty}(\mathbf{K}_{d})\|_{\mathcal{L}_{\infty}}.$ (40)

Note that in the sampled-time case, the order of $P_{\infty}(\mathbf{K}_{d})$
depends on the number of unstable poles of $\mathbf{K}_{d}$ and is lower than
$n_{c}$ when $\mathbf{K}_{d}$ is unstable. In particular, if $\mathbf{K}_{d}$
has $n_{+}$ stable and $n_{-}$ unstable poles and that $m$ is the multiplicity
of the largest unstable Hankel singular value, then
$P_{\infty}(\mathbf{K}_{d})$ has order $n_{+}+n_{-}-m$ (see [23]). To avoid
this issue and for a numerically more robust approach, sub-optimal projection
methods may alternatively be considered (see _e.g._ [21]).

######  Remark 5 (Alternative approach)

A similar problem can be considered using the $\mathcal{L}_{2}$-norm and is
actually much easier to solve considering the decomposition
$\mathcal{L}_{2}(\partial\mathcal{D})=\mathcal{H}_{2}(\mathcal{D})\oplus\mathcal{H}_{2}(\overline{\mathcal{D}})$.
The solution is simply obtained by discarding the unstable part of the model.
However this solution has generally a greater impact on the frequency
behaviour of the model which is not desirable here.

####  5.4.3 Loewner-driven discretisation

Algorithm 1 Loewner-driven discretisation

0: A continuous-time model
$\mathbf{K}\in\mathcal{H}_{\infty}(\mathbb{C}_{+})$, a sampling time $h>0$, an
upper bound $\bar{n}$ of the desired order and a number $m$ of interpolation
points.

1: Sample the interval $]0,\omega_{N}[$ in $2m$ points $\omega_{k}$

2: Evaluate $R(\imath\omega_{k})^{-1}\mathbf{K}(\imath\omega_{k})$

3: Apply Loewner to the data-set (39) to get a first model
$\mathbf{K}_{d}^{m}$ that matches all the data and get the underlying minimal
McMillian order $n$

4: Set $r=\min(n,\bar{n})$

5: Reduce $\mathbf{K}_{d}^{n}$ to $\mathbf{K}_{d}^{r}$ using any interpolation
method

6: Project $\mathbf{K}_{d}^{r}$ onto a stable subspace, _i.e._ compute the
$n_{c}$-th order $\mathbf{K}_{d}=P_{\infty}(\mathbf{K}_{d}^{r})$

7: return the model
$\mathbf{K}_{d}\in\mathcal{H}_{\infty}(\overline{\mathcal{D}})$ of order
$n_{c}\leq r\leq\bar{n}$

The Loewner-driven discretisation process of the controller is summarised in
algorithm 1. The following comments can be added:

  * •

the reduction step from $\mathbf{K}_{d}^{n}$ to $\mathbf{K}_{d}^{r}$ (step
$5$) may also be achieved by standard model approximation techniques as the
one presented in section 3.3 or, more interestingly, using the frequency-
limited version of section 3.4. Then $\mathbf{K}_{d}^{r}$ should be projected
on $\mathcal{H}_{\infty}(\overline{\mathcal{D}})$ and finally a stability
preserving model reduction method could be used to obtain $\mathbf{K}_{d}$.

  * •

As both $\mathbf{K}_{d}^{r}$ and $\mathbf{K}_{d}$ are easily available during
the process, $\|\mathbf{K}_{d}^{r}-\mathbf{K}_{d}\|_{\mathcal{L}_{\infty}}$
can be evaluated to get an estimation of the discretisation error
$e_{\infty}(\mathbf{K},\mathbf{K}_{d})$ in (36).

  * •

Consequently, when $n_{c}=r$, computing
$\|\mathbf{K}_{d}^{r}-\mathbf{K}_{d}\|_{\mathcal{L}_{\infty}}$ gives insight
on whether the maximal allowed order $\bar{n}$ is enough to ensure a low
error.

  * •

If $\mathbf{K}_{d}^{r}$ is unstable, then step $6$ leads to a loss of order
(see section section 5.4.2). Therefore, when $n_{c}<r$ (which is likely to
happen666Indeed, one tries to approximate data coming from an infinite
dimensional model by a rational one. A large order $r$ is generally required
to achieve an exact interpolation.), $r$ should be increased above $\bar{n}$
so that $P_{\infty}(\mathbf{K}_{d})$ is of order $\bar{n}$ thus avoiding any
loss of accuracy.

It should also be noted that algorithm 1 does not exploit the state-space
structure of $\mathbf{K}$ and only requires frequency data from it. Therefore,
as presented in the preliminary work [42], the approach may be applied to a
wider class of models than just those described by a state-space realisation.
algorithm 1 is applied on the obtained continuous-time GLA controller. The
obtained sampled-time controller sigma plot is shown in fig. 6, comparing the
responses of the continuous-time controller $\mathbf{K}$ interconnected with a
holder, with different sampled controllers obtained with different methods
(but same order $n_{c}$). In addition fig. 7 reports the gain and phase
mismatches of the different sampled controllers.

Figure 6: Comparison of the sigma plot of different controllers: the
continuous-time one $\mathbf{K}$ with holder (solid blue), the discrete-time
obtained with backward method $\mathbf{K}_{d}^{zoh}$ (dotted black), the
discrete-time obtained with bilinear method $\mathbf{K}_{d}^{tus}$ (dash
dotted magenta) and the discrete-time obtained with the interpolatory method
method $\mathbf{K}_{d}$ (dashed red).

Figure 7: Gain (left) and phase (right) errors: comparison of the Bode plot of
different sampled controllers with the continuous-time one $\mathbf{K}$ with
holder. The discrete-time obtained with backward method $\mathbf{K}_{d}^{zoh}$
(dotted black), the discrete-time obtained with bilinear method
$\mathbf{K}_{d}^{tus}$ (dash dotted magenta) and the discrete-time obtained
with the interpolatory method method $\mathbf{K}_{d}$ (dashed red).

With reference to fig. 6 and fig. 7, one may observe that the obtained
$\mathbf{K}_{d}$ well reproduces the original controller $\mathbf{K}$
behaviour up to the Nyquist frequency, while the bilinear approach
$\mathbf{K}_{d}^{tus}$ tends to compress the responses, both bilinear and
backward importantly modify the phase. This accurate restitution obtained by
$\mathbf{K}_{d}$ is first justified by embedding the holder model in the
control law before discretisation. Second, as the complete controller is now
delayed, once again the interpolatory framework is adapted and thus
outperforms the classic ones. More important effects on the closed-loop will
be presented in the conclusion section 6.

##  6 Conclusions

In this paper, a complete end to end approach for the construction of a GLA
control function has been presented, applied and implemented on a real-life
generic BizJet aircraft use-case simulator. As a conclusion to this paper,
section 6.1 presents the obtained performances when plugging the discrete-time
controller in the complete industrial simulator including the time-delayed
continuous-time model (hundreds of simulations are performed). As the proposed
design scheme can be readily applied to any linear dynamical systems and not
only aircraft control, in section 6.2 we summarise the main methodological
contributions and potential improvement axes.

###  6.1 BizJet aircraft GLA function performances

The main objective of GLA is to reduce the so-called load envelope, _i.e._ the
load amplification along the wing span in response to a gust signal family
usually considered by aircraft manufacturers and authorities. This alleviation
feature is presented on fig. 8, for different discretised controllers, with
the same dimension and sampling-time value $h$.

Figure 8: Gust load attenuation envelope gain $\mathcal{E}(x_{i})$ resulting
to the interconnection presented in fig. 1 for different GLA controllers: the
continuous-time one $\mathbf{K}$ (red), the discrete-time obtained with
backward method $\mathbf{K}_{d}^{zoh}$ (purple), the discrete-time obtained
with bilinear method $\mathbf{K}_{d}^{tus}$ (cyan) and the discrete-time
obtained with the interpolatory method method $\mathbf{K}_{d}$ (green).

While, the “optimal” gust load alleviation level obtained by the continuous-
time controller $\mathbf{K}$ leads to attenuation between 6 to 13% (except at
the wing tip, where the loads are very small), the sampled-time versions with
a sampling value $h$ fixed by the available on-board computer clearly
deteriorate the performances. Still, depending on the discretisation method,
one clearly observes an important gain when using the proposed discretisation
approach of algorithm 1, with respect to the standard ones. It is to be noted
that every additional percent may represent hundreds of kilograms of mass
reduction, and thus of aircraft global consumption. In addition, as a standard
industrial requirement, the proposed load alleviation should not result in a
change of flight performances. This is ensured by the control design and
illustrated on fig. 9 where the normalised load factor variation in response
to an elevator request $\mathbf{\delta}_{mc}^{\star}(t)$, results in similar
behaviours with and without the GLA function, which fulfils the industrial
process and constraints.

Figure 9: Time-domain response of the interconnection fig. 1 without (solid
blue) and with (dashed red) the GLA function, to a pilot stick input (solid
black).

In this paper, results are shown for few (worst case) configurations, but
obviously, the entire process has been validated over all hundreds of cases
and gust waves. This contributes in addressing both industrial safety and
quality, as well as helping in the certification process. Moreover, all the
presented tools were made available in the industrial process, running
robustly in the workflow.

###  6.2 General contributions overview and improvements axes

The paper considers the challenge of GLA control function design and
implementation while handling industrial constraints. Although this
application is very challenging as it involves a family of complex medium-
scale delayed models and considers different control specifications, both in
time and frequency domains, authors believe that the process can be readily
applied to any kind of complex linear dynamical systems.

Indeed, the main message of the article is to place _interpolatory methods_ at
the center of the modelling, design and discretisation steps. The main
contributions of the paper can then be viewed from three angles. First, from a
_methodological point of view_ , the proposed process is flexible, numerically
affordable and tailored to any (finite and infinite dimensional) linear
dynamical models. Second, from a more _theoretical point of view_ , it
provides a complete view of interpolatory methods for dynamical systems,
applicable to both realisation-based and transfer function-based models.
Finally, algorithm 1 gives a complete process to perform a dynamical
model/controller discretisation resulting in an alternate solution to the
classical literature methods, standing then as a third, _algorithmic_
contribution.

## References

  * [1] M. Alama, M. Hromcikb, and T. Hanis, Active gust load alleviation system for flexible aircraft: Mixed feedforward/feedback approach, Aerospace Science and Technology, 41 (2015), pp. 122 – 133.
  * [2] D. Amsallem and C. Farhat, An online method for interpolating linear parametric reduced-order models, SIAM Journal of Scientific Computing, 33 (2011), pp. 2169–2198.
  * [3] A. Antoulas, The Loewner framework and transfer functions of singular/rectangular systems, Applied Mathematics Letters, 54 (2016), pp. 36–47.
  * [4] A. Antoulas, C. Beattie, and S. Gugercin, Interpolatory methods for model reduction, SIAM Computational Science and Engineering, Philadelphia, 2020\.
  * [5] A. Antoulas, S. Lefteriu, and A. Ionita, Model reduction and approximation theory and algorithms, SIAM, Philadelphia, P. Benner, A. Cohen, M. Ohlberger and K. Willcox Eds, 2016, ch. A tutorial introduction to the Loewner framework for model reduction.
  * [6] A. C. Antoulas, I. V. Gosea, and A. C. Ionita, Model reduction of bilinear systems in the Loewner framework, SIAM Journal on Scientific Computing, 38 (2016), pp. B889–B916.
  * [7] P. Apkarian and D. Noll, Nonsmooth $\mathcal{H}_{\infty}$ Synthesis, IEEE Transaction on Automatic Control, 51 (2006), pp. 71–86.
  * [8] B. Bamieh, J. Boyd Pearson, B. Francis, and A. Tannenbaum, A lifting technique for linear periodic systems with applications to sampled-data control, Systems & Control Letters, 17 (1991), pp. 79 – 88.
  * [9] L. Baratchart, Existence and Generic Properties of $L^{2}$ Approximants for Linear Systems, IMA Journal of Mathematical Control & Information, 3 (1986), pp. 89–101.
  * [10] L. Baratchart, M. Cardelli, and M. Olivi, Identification and rational $\mathcal{L}_{2}$ approximation: A gradient algorithm, Automatica, 27 (1991), pp. 413–417.
  * [11] A. Bunse-Gerstner, D. Kubalinska, G.Vossen, and D.Wilczek, $h_{2}$-norm optimal model reduction for large scale discrete dynamical MIMO systems, Journal of Computational and Applied Mathematics, 233 (2010), pp. 1202–1216.
  * [12] T. Chen and B. Francis, Input-output stability of sampled-data systems, IEEE Transactions on Automatic Control, 36 (1991), pp. 50–58.
  * [13] T. Chen and B. Francis, Optimal sampled-data control systems, Springer Science & Business Media, 1995.
  * [14] R. Gadient, E. Lavretsky, and K. Wise, Very Flexible Aircraft Control Challenge Problem, in Proceedings of the AIAA Guidance, Navigation, and Control Conference, Minneapolis, Minnesota, USA, August 2012.
  * [15] K. A. Gallivan, A. Vanderope, and P. Van Dooren, Model reduction of MIMO systems via tangential interpolation, SIAM Journal of Matrix Analysis and Application, 26 (2004), pp. 328–349.
  * [16] K. Glover, All Optimal Hankel Norm Approximation of Linear Multivariable Systems, and Their $\mathcal{L}_{\infty}$ error Bounds, International Journal Control, 39 (1984), pp. 1145–1193.
  * [17] I. Gosea and A. Antoulas, Stability preserving post-processing methods applied in the Loewner framework, in Workshop on Signal and Power Integrity, 2016, pp. 1–4.
  * [18] S. Gugercin, A. C. Antoulas, and C. A. Beattie, $\mathcal{H}_{2}$ Model Reduction for Large Scale Linear Dynamical Systems, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 609–638.
  * [19] N. Higham, Accuracy and stability of numerical algorithms (2nd edition), SIAM, 2002.
  * [20] T. Kier and G. Looye, Unifying manoeuvre and gust loads analysis models, in Proceedings of the International Forum on Aeroelasticity and Structural Dynamics, July 2009.
  * [21] M. Kohler, On the closest stable descriptor system in the respective spaces $\mathcal{RH}_{2}$ and $\mathcal{RH}_{\infty}$, Linear Algebra and its Applications, 443 (2014), pp. 34–49.
  * [22] R. B. Lehoucq and D. C. Sorensen, Deflation Techniques for an Implicitly Restarted Arnoldi Iteration, SIAM Journal on Matrix Analysis & Applications, 17 (1996), pp. 789–821.
  * [23] J. Mari, Modifications of rational transfer matrices to achieve positive realness, Signal Processing, 80 (2000), pp. 615–635.
  * [24] A. J. Mayo and A. C. Antoulas, A framework for the solution of the generalized realization problem, Linear Algebra and its Applications, 425 (2007), pp. 634–662.
  * [25] C. Meyer, G. Broux, J. Prodigue, O. Cantinaud, and C. Poussot-Vassal, Demonstration of innovative vibration control on a Falcon Business Jet, in Proceedings of the International Forum on Aeroelasticity and Structural Dynamics, Como, Italy, June 2017.
  * [26] D. Ossmann and C. Poussot-Vassal, Design and assessment of a two degree of freedom gust load alleviation system, in Euro GNC, Milan, Italy, April 2019.
  * [27] D. Petersson and J. Lofberg, Model reduction using a frequency-limited $\mathcal{H}_{2}$-cost, System & Control Letters, 67 (2014), pp. 32 – 39.
  * [28] I. Pontes Duff, C. Poussot-Vassal, and C. Seren, Optimal $\mathcal{H}_{2}$ model approximation based on multiple input/output delays systems, Systems & Control Letters, 117 (2018), pp. 60–67.
  * [29] C. Poussot-Vassal, Large-scale dynamical model approximation and its applications, HDR, habilitation thesis, Onera, INP Toulouse, Toulouse, France, July 2019.
  * [30] M. Pusch, Allocation of Distributed Flaps for Gust Load Alleviation, in Proceedings of the Conference on Control Technology and Application, Waimea, USA, July 2018, pp. 27.8–30.8.
  * [31] M. Pusch, D. Ossmann, J. Dillinger, T. Kier, M. Tang, and J. Lubker, Aeroelastic modeling and control of an experimental flexible wing, in Proceedings of the AIAA Scitech 2019 Forum, 2019, p. 0131.
  * [32] D. Quero, An Aeroelastic Reduced Order Model for Dynamic Response Prediction to Gust Encounters, PhD thesis, Technical University of Berlin, 2017.
  * [33] D. Quero, P. Vuillemin, and C. Poussot-Vassal, A Generalized State-Space Aeroservoelastic Model based on Tangential Interpolation, Aerospace, (2019).
  * [34] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2000.
  * [35] Y. Saad and M. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM Journal of Scientific and Statistical Computing, 7 (1986), pp. 856–869.
  * [36] G. Sleijpen and H. Van der Vorst, A Jacobi-Davidson iteration method for linear eigenvalue problems, SIAM Journal on Matrix Analysis and Applications, 17 (1996), pp. 401–425.
  * [37] J. Spanos, M. Milman, and D. Mingori, A new algorithm for $\mathcal{L}_{2}$ optimal model reduction, Automatica, 28 (1992), pp. 897–909.
  * [38] P. Van Dooren, K. A. Gallivan, and P. A. Absil, $\mathcal{H}_{2}$-optimal model reduction of MIMO systems, Applied Mathematics Letters, 21 (2008), pp. 53–62.
  * [39] P. Van Dooren, K. A. Gallivan, and P. A. Absil, $\mathcal{H}_{2}$-optimal model reduction with higher order poles, SIAM Journal on Matrix Analysis and Applications, 31 (2010), pp. 2738–2753.
  * [40] C. D. Villemagne and R. E. Skelton, Model Reduction using a Projection Formulation, International Journal Control, 40 (1987), pp. 2141–2169.
  * [41] P. Vuillemin, Frequency-limited approximation of large-scale LTI dynamical models, Ph.D. thesis, Onera, ISAE, Toulouse University, Toulouse, France, November 2014.
  * [42] P. Vuillemin and C. Poussot-Vassal, Discretisation of continuous-time linear dynamical model with the Loewner interpolation framework, arXiv 1907.10956, (2019).
  * [43] P. Vuillemin, C. Poussot-Vassal, and D. Alazard, $\mathcal{H}_{2}$ optimal and frequency limited approximation methods for large-scale LTI dynamical systems, in Proceedings of the 6th IFAC Symposium on Systems Structure and Control, Grenoble, France, February 2013, pp. 719–724.
  * [44] P. Vuillemin, C. Poussot-Vassal, and D. Alazard, Spectral expression for the Frequency-Limited $\mathcal{H}_{2}$-norm of LTI Dynamical Systems with High Order Poles, in Proceedings of the 13th European Control Conference, Strasbourg, France, June 2014, pp. 55–60.
  * [45] K. Willcox, J. Peraire, and J. White, An Arnoldi approach for generation of reduced-order models for turbomachinery, Journal of Computers & Fluids, 31 (2002), pp. 369–389.