NOMENCLATURE

Flight mechanics

CG

centre of gravity of the aircraft

g

gravitational acceleration ($\mathrm{m/s}^2$)

$I_{\cdot \cdot}$

aircraft moments and products of inertia with respect to body axis system ($\cdot = $ x, y or z) ($\mathrm{kg\,m}^2$)

$L, M, N$

components of the resultant moment about the CG in the body axis system (N m)

m

aircraft mass (kg)

$p, q, r$

components of the aircraft angular velocity vector in the body axis system (rad/s)

$u_K, v_K, w_K$

components of the flight-path velocity vector relative to the Earth in the body axis system (m/s)

$V_{\rm T\!AS}$

true airspeed (m/s)

$w_{\rm gust}$

vertical wind component from the gust (m/s)

$X, Y, Z$

components of the resultant force vector in the body axis system (N)

xyz

body axis system (fixed to the aircraft, with the CG as origin)

$x, y, z$

components of the inertial position vector (inertial displacements) (m)

$\alpha$

angle-of-attack (°)

$\delta_{e,{\rm cmd}}, \delta_{a,{\rm cmd}}^{\rm sym}$

commanded elevator and symmetric aileron deflections, respectively (°)

$\delta_{e,{\rm act}}, \delta_{a,{\rm act}}^{\rm sym}$

elevator and symmetric aileron deflections (output of the corresponding flight control actuator), respectively (°)

$\Phi, \Theta, \Psi$

Euler angles (bank, inclination and azimuth angles, respectively) (°)

Structural dynamics

BR1, BR4

bending moments at stations WR1 and WR4, respectively (N m)

$d_E$

total elastic deformation

$m_i, Q_i$

generalised mass (kg) and force (N) associated with the ith flexible mode, respectively

$n_z^{\rm HTR}$

normal load factor at station HTR

$n_z^{\rm pilot}$

normal load factor at pilot location

$\zeta_i, \omega_i$

modal damping (–) and natural frequency (rad/s), respectively

$\eta_i, \phi_i$

generalised coordinate and mode shape (eigenfunction), respectively

Control

$a_z^{\rm IMU}$

vertical acceleration measurements from the IMU ($\mathrm{m/s}^2$)

$\textbf{\textit{d}}, {\textbf{\textit{d}}}_p$

disturbance input and previewed disturbance input available for the control action, respectively

$q^{\rm IMU}$

pitch rate measurements from the IMU (rad/s)

h

preview horizon (continuous time) or preview length (discrete time)

$H_{\infty}, \| \cdot \|_{\infty}$

infinity norm of a system

$T_{zw}$

closed-loop transfer matrix from the exogenous input to the regulated output

${\textbf{\textit{u}}}, {\textbf{\textit{w}}}$

control and exogenous input vectors, respectively

${\textbf{\textit{x}}}$

state vector of the system

${\textbf{\textit{y}}}, {\textbf{\textit{z}}}$

measured and regulated output vectors, respectively

${\textbf{\textit{y}}}_{\textbf{\textit{d}}}$

preview measurements from the LIDAR

$\gamma$

$H_{\infty}$ performance value

$\boldsymbol{\xi}$

augmented state vector

2.0 DYNAMICS OF FLEXIBLE AIRCRAFT

Most modern aircraft feature high-aspect-ratio wings, a slender fuselage and a thin, lightweight structure, resulting in significant structural flexibility. There are two main reasons why aeroelasticity needs to be accounted for: first, the coupling between the frequencies of structural (or flexible) and rigid-body modes of the aircraft can, in the presence of non-linearities, introduce effects on the flight dynamics, leading to the degradation of performance and handling quality^{(Reference Baghdadi, Lowenberg and Isikveren28)}; and second, the structural modes spoil sensor measurements and can potentially reduce stability margins. Hence, the structural flexibility and structural modes have to be accounted for in the design of Flight Control Systems (FCSs). The problem of FCS–structural coupling and a proposed solution for it were discussed in more detail in Ref. [Reference Caldwell29]. The authors in Ref. [Reference Waszak and Schmidt30] developed an aeroelastic model and applied it to the flight dynamics analysis of the Rockwell B-1 bomber aircraft. The mathematical model used the free vibration modes of the aircraft and was represented in the mean-axis system^{(Reference Schmidt31)} to minimise inertial coupling, while the coupling was instead captured via the aerodynamic forces. The structure of this aeroelastic model, summarised in the following paragraphs, is used in this paper.

The rigid-body dynamics of an aircraft expressed in the body axis system, the rotational kinematics expressed based on Euler angles and the translational kinematics expressed in a geodetic frame (assuming a flat, fixed Earth) can be represented by Equations (1), (2) and (3), respectively.

(1) $\begin{equation}\begin{array}{rcl}X&=&m(\dot{u}_K+q\ w_K-r\ v_K)+m\ g\sin{\Theta}\\[3pt]Y&=&m(\dot{v}_K+r\ u_K-p\ w_K)-m\ g\cos{\Theta}\sin{\Phi}\\[3pt]Z&=&m(\dot{w}_K+p\ v_K-q\ u_K)-m\ g\cos{\Theta}\cos{\Phi}\\[3pt]L&=&I_{xx}\ \dot{p}-I_{xy}(\dot{q}-p\ r)-I_{xz}(\dot{r}+p\ q)-I_{yz}\!\left(q^2-r^2\right)+(I_{zz}-I_{yy})q\ r\\[3pt]M&=&I_{yy}\ \dot{q}-I_{xy}(\dot{p}+q\ r)-I_{yz}(\dot{r}-p\ q)+I_{xz}\!\left(p^2-r^2\right)+(I_{xx}-I_{zz})p\ r\\[3pt]N&=&I_{zz}\ \dot{r}-I_{xz}(\dot{p}-q\ r)-I_{yz}(\dot{q}+p\ r)-I_{xy}\!\left(p^2-q^2\right)+(I_{yy}-I_{xx})p\ q\end{array}\end{equation}$

(2) $\begin{equation}\hspace*{-137pt}\begin{array}{rcl}\dot{\Phi}&=&p+q(\sin{\Phi}\tan{\Theta})+r(\cos{\Phi}\tan{\Theta})\\[3pt]\dot{\Theta}&=&q\cos{\Phi}-r\sin{\Phi}\\[3pt]\dot{\Psi}&=&q(\sin{\Phi}\sec{\Theta})+r(\cos{\Phi}\sec{\Theta})\end{array}\end{equation}$

(3) $\begin{equation}\hspace*{-80pt}\begin{array}{rcl}\dot{x}&=&u_K(\cos{\Theta}\cos{\Psi})\\[3pt]&&+v_K(\sin{\Phi}\sin{\Theta}\cos{\Psi}-\cos{\Phi}\sin{\Psi})\\[3pt]&&+w_K(\cos{\Phi}\sin{\Theta}\cos{\Psi}+\sin{\Phi}\sin{\Psi})\\[3pt]\dot{y}&=&u_K(\cos{\Theta}\sin{\Psi})\\[3pt]&&+v_K(\sin{\Phi}\sin{\Theta}\sin{\Psi}+\cos{\Phi}\cos{\Psi})\\[3pt]&&+w_K(\cos{\Phi}\sin{\Theta}\sin{\Psi}-\sin{\Phi}\cos{\Psi})\\[3pt]\dot{z}&=&-u_K\sin{\Theta}+v_K(\sin{\Phi}\cos{\Theta})+w_K(\cos{\Phi}\cos{\Theta})\end{array}\end{equation}$

For the structural dynamics, the total elastic displacement of the structure can be expressed in terms of a modal expansion using n free vibration modes as given by Equation (4), whereas the n generalised coordinates associated with the respective n free vibration modes are governed by the n equations given in Equation (5). By using the mean-axis system, the aeroelastic model is constituted of Equations (1), (2), (3) and (5), which are 12+n in number, have 12+2n states and are non-linear coupled differential equations of first and second order.

(4) $\begin{equation}d_E(x,y,z,t)=\sum_{i=1}^n {\big(\,\phi_i(x,y,z)\ \eta_i(t)\,\big)}\end{equation}$

(5) $\begin{equation}\ddot{\eta}_i+2\ \zeta_i\ \omega_i\ \dot{\eta}_i+\omega_i^2\ \eta_i=\frac{Q_i}{m_i}\ ,\ \ i=1,2,\ldots, n\end{equation}$

Although the current work is primarily intended to be applied to large transport aircraft (Part 25 of the FAA airworthiness standards⁽¹⁾ or CS-25 of the EASA certification specifications⁽²⁾), an aeroelastic model of a sailplane is used in this paper. The sailplane is DLR’s Discus-2c, shown in Fig. 1. The general mass and geometry characteristics, as well as the modal characteristics of the sailplane can be found in Ref. [Reference Khalil32]. A flight test measurement system has been installed on the aircraft and allowed for the development of a non-linear aeroelastic model of the discussed structure (Equations (1), (2), (3) and (5)) in the work of Ref. [Reference Viana33] using system identification techniques. This allowed for the calculation of shear forces, as well as torsional and bending moments at seven load stations, of which only three (WR1, WR4 and HTR, shown in Fig. 2) are considered in this paper. The non-linear aeroelastic model accounts for the change in the local angle-of-attack of the horizontal tail due to the aircraft pitching motion as well as the time delay of the air flow reaching the horizontal tail (by delaying the downwash angle through a first-order approximation of the angle-of-attack); see Equations (5.16) and (5.17) of Ref. [Reference Viana33] or Equations (12) and (13) of Ref. [Reference Viana34] for more detail. A linear aeroelastic model to be used for the control synthesis has been obtained and validated with the non-linear aeroelastic model in Ref. [Reference Khalil32].

Figure 1. DLR’s Discus-2c sailplane in flight.

Figure 2. Load stations and distribution of flight test sensors.

At the time of flight tests used for obtaining the non-linear aeroelastic model, no actuators were installed on the aircraft. Instead, the control surfaces were mechanically connected to the pilot’s stick, and the pilot performed the flight manoeuvres by manually manipulating the stick. If however a GLA system is to be applied, control surface actuators and automatic manipulation are required. DLR has now integrated and tested elevator and aileron actuators. On-ground and also in-flight tests have been performed to identify the elevator and aileron actuators dynamics. Based on these tests, the actuators dynamics can be well represented by the first-order dynamics given by Equation (6).

(6) $\begin{equation}\begin{array}{c}G_{\delta_e}^{\rm act} = \dfrac{\delta_{e,{\rm act}}}{\delta_{e,{\rm cmd}}} = \dfrac{1.14048}{s+1.14048},\quad G_{\delta_a}^{\rm act} = \dfrac{\delta_{a,{\rm act}}^{\rm sym}}{\delta_{a,{\rm cmd}}^{\rm sym}} = \dfrac{7.7105}{s+7.7105}\end{array}\end{equation}$

3.0 PREVIEW CONTROL

In Section 1, it was argued that, by measuring the turbulence or gusts at an adequate distance ahead of the aircraft, the performance of GLA systems can be enhanced. To achieve this, it is first required to know how to measure the turbulence or gusts and how to integrate this measured signal into the control problem. A sensor technology that enables the measurement of turbulence or gusts is Doppler LIDAR. In this type of sensor, a laser beam with a previously known wavelength/frequency is transmitted and the back-scattered light (by aerosols and/or air molecules in the atmosphere) is then received. The velocity of the aerosols/molecules – which is correlated to the air velocity – can be calculated from the Doppler shift between the frequencies of the transmitted and received laser beams. Hence, the wind speed and direction can be determined from the measurements. This paper assumes that the wind field ahead of the aircraft is available to the GLA system. Hence, the focus is rather on the optimal exploitation of this measured signal rather than on how to obtain it.

From the discussion above, it can be concluded that the sensor measures the wind at the current time but at a location that the aircraft has not yet reached. Based on the remote wind information gathered ahead of the aircraft at the present time and in the past, and based on the aircraft motion, a ‘best guess’ of the future encountered wind can be made. The control problem can then be represented by the block diagram shown in Fig. 3. From this figure, it can be seen that the future encountered wind, which the controller ‘knows’ in advance, will affect the aircraft after some time delay h.

Figure 3. Control problem with previewed disturbance.

Figure 4. Standard form for control synthesis.

The problem of structural load alleviation and ride quality enhancement of an aircraft when encountering turbulence or gusts is by nature a disturbance rejection problem. In a disturbance rejection problem, the goal is to minimise the impact of the disturbance (e.g. gust or turbulence) on the regulated outputs, whose variations in response to the disturbance must be kept as small as possible (e.g. structural loads or accelerations at various locations in the aircraft cabin). For this purpose, an adequate measure of the transfer function, i.e. a system norm, must be selected and then minimised. Two of the most widely used system norms are the $H_2$ and $H_{\infty}$ norms. The $H_2$ norm can be interpreted as an average of the system gain taken across all frequencies, whereas the $H_{\infty}$ norm is a measure of the peak gain of the system (i.e. the largest amplification that the system may show to an input vector across all frequencies and all input directions, the worst-case scenario). As the gust input is neither fixed nor has a fixed power spectrum (i.e. the objective is to minimise the peak loads), the $H_{\infty}$ norm is better suited here. Hence, the considered control problem can be well expressed as an $H_{\infty}$-optimal control design problem, shown in Fig. 4. P is the weighted plant (i.e. the plant plus the weighting functions that the control designer selects to specify the control objectives) and K is the controller to be synthesised. The objective of the $H_{\infty}$ control design problem is to find the controller K that yields a stable closed-loop system $T_{{\textbf{\textit{zw}}}}$ and ensures that its infinity norm is less than a positive $H_{\infty}$ performance value $\gamma$ (i.e. $\|T_{{\textbf{\textit{zw}}}}\|_{\infty} < \gamma$). The plant P can be expressed in state-space form as given by Equation (7), where the matrices A, $B_1$, $B_2$, $C_1$, $C_2$, $D_{11}$, $D_{12}$, $D_{21}$ and $D_{22}$ are of appropriate dimension.

(7) $\begin{align} \dot{{\textbf{\textit{x}}}}(t) & = {\textbf{\textit{A}}}\ {\textbf{\textit{x}}}(t) + {{\textbf{\textit{B}}}_1}\ {\textbf{\textit{w}}}(t) + {{\textbf{\textit{B}}}_2}\ {\textbf{\textit{u}}}(t)\nonumber\\[3pt] P :=\ {\textbf{\textit{z}}}(t) & = {{\textbf{\textit{C}}}_1}\ {\textbf{\textit{x}}}(t) + {{\textbf{\textit{D}}}_{11}}\ {\textbf{\textit{w}}}(t) + {{\textbf{\textit{D}}}_{12}}\ {\textbf{\textit{u}}}(t)\\[3pt] {\textbf{\textit{y}}}(t) & = {{\textbf{\textit{C}}}_2}\ {\textbf{\textit{x}}}(t) + {{\textbf{\textit{D}}}_{21}}\ {\textbf{\textit{w}}}(t) + {{\textbf{\textit{D}}}_{22}}\ {\textbf{\textit{u}}}(t)\nonumber\end{align}$

The classical methods for synthesising $H_{\infty}$ $ controllers (by solving either algebraic Riccati equations^{(Reference Doyle, Glover, Khargonekar and Francis35)} or linear matrix inequalities^{(Reference Gahinet and Apkarian36)}) result in state-space controllers that have the same order as that of the weighted plant (i.e. full-order methods). To use these methods, all the system inputs and outputs should be grouped into a single channel, which can include unwanted cross-coupling terms in the overall transfer function. To overcome these disadvantages of full-order methods, researchers have investigated other methods, leading to a new $H_{\infty}$ $-optimal control method that relies on non-smooth optimisation techniques^{(Reference Apkarian and Noll37)}. The new method allows one to the choose the controller structure and order (i.e. fixed-structure method) as well as expressing the problem in a multi-channel formulation. The control problem can be illustrated in the case of two channels as shown in Fig. 5, for which an optimisation problem of the form given by Equation (8) is solved. The multi-channel formulation can be used for: first, achieving a trade-off between various $H_{\infty}$-norm criteria (in this case, the plants $P_1$, $P_2$, etc, are representing the same plant P, but with different inputs $w_i$ and outputs $z_i$; and the controllers $K_1$, $K_2$, etc, are the same controller K); and second, performing robust control design by simultaneously considering several models (e.g. corresponding to various operating points or with different values for some uncertain parameters). In this paper, both the full-order and fixed-structure methods will be used, where applicable^{Footnote 2} and compared with each other. This allows for the achievement of the required design flexibility by using the fixed-structure method while examining the optimality of the synthesised controller via comparison with that of the full-order methods.

(8) $\begin{equation}\begin{array}{l}\textrm{minimize } \quad \max\limits_{i}\left\{\left\|T_{z_{1}w_{1}}(P_1,K_1)\right\|_{\infty}\!, \left\|T_{z_{2}w_{2}}(P_2,K_2)\right\|_{\infty}\right\}\\[10pt] \textrm{subject to} \quad K_1\ \textrm{and}\ K_2\ \textrm{stabilise}\ P_1\ \textrm{and}\ P_2, \textrm{respectively} \end{array}\end{equation}$

Figure 5. Two-channel control synthesis.

Figure 6. Standard form for preview control synthesis.

The preview control problem can typically be represented by the block diagram shown in Fig. 6. The delay term can be accounted for in either the continuous^{(Reference Kojima and Ishijima38)} or discrete time^{(Reference Takaba39)} domains. Contrary to the continuous-time approach, the discrete-time approach provides a closed-form controller and allows the problem to be formulated in the context of full-order as well as fixed-structure methods. Hence the discrete-time approach is followed in this paper. However, in Ref. [Reference Takaba39], the problem was formulated as a full-information state-feedback control problem and then solved using a state augmentation technique. Since in our case the model dynamics are for flexible aircraft (i.e. they include flexible degrees of freedom that cannot be measured), the problem needs to be formulated in the output-feedback setting. The same state augmentation technique is still used. A complete derivation of the preview control problem can be found in Ref. [Reference Khalil40], which yields the discrete-time system given by Equation (9). This system is in the standard form of a $H_{\infty}$ control problem given by Equation (7), hence the characterisation of solvability relies exclusively on that of $H_{\infty}$ control theory.

(9) $\begin{align} \boldsymbol{\xi}(k+1) & = \overline{{\textbf{\textit{A}}}}\ \boldsymbol{\xi}(k) + \overline{{\textbf{\textit{B}}}}_{\boldsymbol{1}}\ {{{\textbf{\textit{d}}}_p}}(k) + \overline{{\textbf{\textit{B}}}}_{\boldsymbol{2}}\ {\textbf{\textit{u}}}(k)\nonumber\\[3pt] {\textbf{\textit{z}}}(k) & = \overline{{{\textbf{\textit{C}}}}}_{\boldsymbol{1}}\ \boldsymbol{\xi}(k) + \overline{{{\textbf{\textit{D}}}}}_{\boldsymbol{11}}\ {{{\textbf{\textit{d}}}_p}}(k) + \overline{{{\textbf{\textit{D}}}}}_{\boldsymbol{12}}\ {\textbf{\textit{u}}}(k)\\[3pt] \overline{{\textbf{\textit{y}}}}(k) & = \overline{{{\textbf{\textit{C}}}}}_{\boldsymbol{2}}\boldsymbol{\xi}(k) + \overline{{{\textbf{\textit{D}}}}}_{\boldsymbol{21}}\ {{{\textbf{\textit{d}}}_p}}(k) + \overline{{{\textbf{\textit{D}}}}}_{\boldsymbol{22}}\ {\textbf{\textit{u}}}(k)\nonumber\end{align}$

4.0 CONTROL DESIGN

4.2 Application

For step 1, the chosen trim flight condition is steady rectilinear flight defined by $V_{\rm T\!AS}=$ 160km/h and $h=$ 1,000m. In this paper, only the symmetric motion and loads, the short-period mode for the rigid-body motion, as well as the first and third symmetric modes (first and second wing vertical bending, respectively) for the flexible degrees of freedom are considered. The bare linear state-space model can be defined again as given by Equation (7), but without the disturbance input w(t).

For step 2, since only the symmetric motion and loads are considered here, the exogenous disturbance input w(t) is chosen to be a vertical gust or turbulence. Assuming small angles, this vertical gust might be subtracted from the vertical component of the inertial velocity vector of the aircraft in the body axis system, $w_K(t)$. This means that the ${{\textbf{\textit{B}}}_1}$ matrix in Equation (7) will be equal to the first column of the ${\textbf{\textit{A}}}$ matrix. Since there is no direct transfer from the gust or turbulence input to the regulated outputs and measurements, the matrices ${{\textbf{\textit{D}}}_{11}}$ and ${{\textbf{\textit{D}}}_{21}}$ are set to zero matrices of appropriate dimension. Hence, the system model is written as in Equation (7).

In step 3, the regulated outputs are normalised to facilitate the expression of the design preferences. A typical method for such normalisation is to divide them by their corresponding load envelope values, which are not available to the authors. So in this work, a discrete-gust design criterion from the airworthiness standards (Section 25.341 in Ref. [1]) or certification specifications (CS 25.341 in Ref. [2]) is followed for the normalisation. According to this criterion, a gust of $1-\cos$ shape, given by Equation (10), with a sufficient number of different gust gradient distances H in the range of 30–350ft (9–107m) must be investigated to find the critical response for each load quantity. $U_{\rm ds}$ is the design gust velocity in equivalent airspeed defined by Equation (11), s is the gust penetration distance, $U_{\rm ref}$ is the reference gust velocity and $F_g$ is the flight profile alleviation factor^(1,2) . According to the criterion, $U_{\rm ref}$ is to vary linearly with altitude with $H_{\rm ref} =$ 350ft (107m). $F_g$ is assumed to be equal to $0.5$. Since this criterion is for transport aircraft but the model aircraft used for the simulation is a sailplane, H and $U_{\rm ds}$ are scaled down by the ratio of the trim speed of DLR’s Discus-2c (160km/h) to that of a transport aircraft (assumed Mach 0.5). The $1-\cos$ vertical gust input is shown in Fig. 7.

(10) $\begin{equation}U = \left\{\begin{array}{c@{\quad}c} \dfrac{U_{\rm ds}}{2}\left[1-\cos\left(\dfrac{\pi s}{H}\right)\right] & 0 \leq s \leq 2H\\[15pt] 0 & s > 2H\end{array} \right.\end{equation}$

(11) $\begin{equation} \hspace*{-46pt} U_{\rm ds} = U_{\rm ref} F_g \left(\dfrac{H}{H_{\rm ref}}\right)^{1/6}\hspace*{40pt}\end{equation}$

From Fig. 7, note that, as the gust gradient distance increases (i.e. the frequency of the $1-\cos$ gust input decreases), the amplitude increases. To reflect this in the control synthesis, a frequency-shaping function must be added to the exogenous input of the plant (the gust input here) to penalise the low-frequency inputs more than the high-frequency inputs. The appropriate weighting function can be determined from the amplitude–frequency relation at each gust gradient distance in Fig. 7. By doing so, the frequency-shaping function is given by Equation (12).

(12) $\begin{equation}\begin{array}{c}G_{w}^{F} = \dfrac{25(s+3)}{(s+2.5)(s+30)}\end{array}\end{equation}$

Finally, to achieve acceptable control actuation, high-frequency control actions should also be penalised. This can be done by weighting the control actions with some functions, and then defining them as output performance channels to be included in the computed $H_{\infty}$ norm. The functions chosen here are given by Equation (13). The final obtained continuous-time system will have a state vector of length equal to 12 (2 rigid-body states, 4 flexible degrees of freedom, 2 states for actuator dynamics, 2 states for frequency-shaping function and 2 states for control-action weighting functions).

(13) $\begin{equation}\begin{array}{c}G_{\delta_e}^{F} = \dfrac{u_{\delta_e}^{F}}{\delta_{e,{\rm cmd}}} = \dfrac{0.5(s+0.01)}{s+1},\quad G_{\delta_a}^{F} = \dfrac{u_{\delta_a^{\rm sym}}^{F}}{\delta_{a,{\rm cmd}}^{\rm sym}} = \dfrac{s+0.01}{s+7}\end{array}\end{equation}$

Figure 7. $1-\cos$ vertical gust.

Figure 8. Singular values of the frequency response.

In step 4, an adequate sampling time (or frequency) is chosen for the plant discretisation. The minimum sampling frequency at which a system can be sampled without introducing errors (also known as the Nyquist rate) should be greater than twice the highest frequency of that system. For the considered aircraft and trim condition, the highest frequency is approximately 53rad/s or 8.4Hz. This means that the minimum sampling frequency is approximately 16.8Hz, or the maximum sampling time is approximately 0.059s. This is the minimum required sampling frequency; but, to accurately capture the shape of a signal, the sampling frequency should ideally be at least ten times the highest frequency, which is $\approx$84Hz for the case considered here. This means that the sampling time must be at most $\approx$0.011s. In Fig. 8, the singular values of the frequency response of the open-loop system to the gust input are compared for two sampling times: 0.05 and 0.01s. It can be seen that, for a sampling time of 0.05s, the frequency response of the discrete-time system matches that of the continuous-time system with good accuracy only up to frequencies of $\approx$20rad/s, which is lower than the highest frequency of the system. At higher frequencies, the matching is no longer good and cannot be accepted.

This result can be confirmed by looking at the time response of the outputs of the open-loop system (for instance, the change in the bending moment at station WR1 from its trim value) to a unit-step gust input for the same two sampling times, as shown in Fig. 9. It can be seen that, in case of the sampling time of 0.05s, the shape of the signal is not accurately represented. Consequently, in this work, the sampling time is chosen to be 0.01s. The final obtained discrete-time system will be as given by Equation (9). Its state vector is composed of the original system states (here 12 states), augmented with the additional states from the h unit delays (h states). Its measurement vector ${\textbf{\textit{y}}}$ is composed of the considered feedback measurements (here $q^{\rm IMU}$ and $a_z^{\rm IMU}$), augmented with the additional measurements from the preview signal ($h+1$ measurements).

Figure 9. Time response to unit step (1m/s) vertical gust input.

FUNDING SOURCES

This study was carried out during the research stay of the first author at the Institute of Flight Systems of DLR (German Aerospace Center) and was funded by the Egyptian Ministry of Higher Education and Scientific Research and the DAAD (German Academic Exchange Service, Deutscher Akademischer Austauschdienst) via the GERLS (German-Egyptian Research Longterm Scholarship) program. Special thanks go to them. The work of the second author has been funded within the frame of the Joint Technology Initiative JTI Clean Sky 2, AIRFRAME Integrated Technology Demonstrator platform ‘AIRFRAME ITD (contract no. CSJU-CS2-GAM-AIR-2014-15-01 Annex 1, Issue B04, 2 October 2015), being part of the Horizon 2020 research and Innovation framework programme of the European Commission. The methodology presented in this paper will be applied to large transport aircraft and business jets within the CleanSky2 NACOR project funded through, among many others, the aforementioned contract.