Frequency Response Estimation for Multiple Aircraft Control Loops Using Orthogonal Phase-Optimized Multisine Inputs

Abstract

The latest advances made at NASA for simultaneous excitation in multiple axes and the identification of frequency responses for aircraft flight systems are discussed in this paper. These techniques are extended in the prospect of identifying multiple dynamics and control loops for multiple axes. Recent applications with flight test data and simulation data are also presented, along with a discussion of the practical aspects of the approach. A demonstration is also provided, using a simulation model for the X-59 airplane in which frequency responses for the bare-airframe, closed-loop, and broken-loop (both for the actuator and the sensor) dynamics were identified from a single 60 s maneuver. The results indicate that this approach can significantly shorten the duration of flight tests and Monte Carlo simulations to save time and costs, and can produce results in real time.

1. Introduction

Aircraft system identification is the determination of models that describe an aircraft and its components from experimental flight test data. For example, common products of an aircraft system identification analysis include estimated aerodynamic models for the bare-airframe dynamics, stability margins other empirical frequency responses, input–output transfer functions, and physics-based models in state-space form. (In this paper, “bare-airframe” dynamics refer to the motion of the aircraft resulting from control effector settings, and are usually written in terms of stability and control derivatives or their nonlinear generalizations [1,2,3]).

These products of system identification analysis have many important and practical uses, as summarized in refs. [2,4,5,6]. Identified models can be used for designing or tuning feedback control laws, performing various flight dynamics analyses, training pilots, rehearsing mission profiles, etc. The results can also be used to verify control laws, empirically assess closed-loop robustness via stability margins or evaluate compliance with specifications for handling qualities [5,7]. Reference [8] discusses more applications and specific examples, and notes that as aircraft become increasingly complicated and higher performance control systems are demanded, the required fidelity of the aircraft model also increases.

Ideally, flight tests for identifying these models are conducted in an efficient manner to minimize testing, which in turn reduces time and costs. In traditional tests for aircraft system identification, one input of interest is excited at a time in a sequential manner, while the output measurements are observed and recorded. Examples include the Rockwell-MBB X-31A [9,10], the Boeing X-48B [11,12,13], the Grumman X-29A [14], the Bell XV-15 [15], and the Learjet LJ-25D [16]. For instance, when identifying the closed-loop dynamics, the pilot longitudinal stick could be used for one maneuver, then the pilot lateral stick could be used for a second maneuver, and lastly the rudder pedals could be used for a third maneuver. When other control loops are to be identified, such as for computing stability margins, then other inputs are also applied, again one at a time, in other parts of the flight control system using separate maneuvers. This process is repeated for many combinations of speed, altitude, and vehicle weight over the flight envelope. As the number of control effectors and sensors increases, such as for aeroelastic aircraft [17,18,19,20] and urban air mobility (UAM)-type aircraft [21,22], the total duration of the flight test needed for system identification also increases. Due to these reasons and other factors, flight testing for system identification can be an expensive and time-consuming process.

Over the last decade at NASA, an approach has evolved for efficiently testing aircraft with many inputs for system identification analysis. The approach is based on the use of orthogonal phase-optimized multisine excitations on all inputs of interest, simultaneously. The relevant frequency responses are then computed as ratios of Fourier transform data. This formulation has led to several successful flight tests and wind tunnel tests in which frequency responses for different loops were identified from a single maneuver, which saved significant time and costs over using sequential single-input maneuvers in the traditional manner. This approach has also enabled real-time estimates of the frequency responses, which are useful for fault detection or the monitoring of stability margins. The approach has also been applied to numerical analyses, including the computation of stability margins from nonlinear simulations within a Monte Carlo simulation environment and simplifying complicated computational fluid dynamics (CFD) solvers with many inputs to reduced-order models.

In this paper, the latest advances that have been made at NASA in conducting efficient aircraft flight tests for frequency response estimation are presented. The techniques are explained and several examples are shown. In addition, these techniques are extended to the prospect of simultaneously identifying multiple control loops over multiple axes using a single maneuver. This is demonstrated through an example using the X-59 Quiet SuperSonic Technology (QueSST) vehicle model. Such techniques can result in a considerable reduction in test duration, thereby reducing time and costs.

The remainder of this paper is organized as follows. Section 2 presents the technical material of the methods that have been developed. The problem statement is discussed first, followed by the multisine input design, Fourier transforms, frequency response approaches with examples, real-time implementation, and the simultaneous estimation of multiple control loops. Section 3 provides an example using a simulation model of the X-59 airplane, in which the closed-loop and bare-airframe dynamics and broken-loop paths (opened at both the mixer input and the sensor output) were identified from a single maneuver. Section 4 then concludes the paper.

Several analyses and designs described in this paper used routines from the MATLAB^®-based software package called System IDentification Programs for AirCraft (SIDPAC), which is associated with ref. [2] and available through the NASA Software Catalog [23].

2. Methods

2.1. Problem Statement

Figure 1 shows a generic block diagram of a modern aircraft. The pilot interacts with the aircraft through the control inceptors, such as the stick and rudder pedals. Based on these inputs and feedback from the measured aircraft responses, the flight control system (FCS) produces a set of force and moment commands. In many aircraft, a mixer is then used to distribute these commands amongst multiple control effectors. Afterward, commands are sent to individual actuators, which move the control effectors and impart forces and moments on the aircraft bare airframe. The resulting motion is observed by sensors and fed back to the flight control system.

For most system identification applications, up to three types of loops are of interest. The first loop is for the bare-airframe dynamics of the aircraft, which describe how the aircraft responds to the control effectors. As the physics of the problem are generally well understood [1], the ultimate goal of this analysis is usually to estimate stability and control the derivatives that quantify the aerodynamic forces and moments for small perturbations from a reference condition.

The second loop of interest is for the closed-loop dynamics, which describe the response of the aircraft to the pilot inputs. These are helpful for verifying a control design, obtaining a low-order equivalent system of the aircraft [2,5,24], and assessing compliance with specifications for handling qualities [25].

The third type of loop is for the broken-loop dynamics. For a multiple input multiple output (MIMO) system, such as an aircraft, these loops are obtained by opening a single input single output (SISO) open-loop path at a particular point while keeping the remaining loops closed. Loops are traditionally opened at the inputs to the mixer [5], in part because that is where the loops are most similar to a set of decoupled SISO systems. However, it has become common practice at NASA to also check loop openings at the sensor outputs. This is also used for rotorcraft with attitude-hold modes to assess the disturbance rejection bandwidth [5,26].

In addition to these three types of loops, other loops are sometimes identified from flight test data. Although actuators are typically characterized from ground-based tests, these models can be updated using system identification when the commands and responses are measured during flight. Similarly, the programming of the flight control system can be verified (or identified for adaptive control systems) using recorded commands and responses.

To identify the dynamics of interest, appropriately designed excitation inputs [2,5] are applied at specific points in the system. The excitation points used in this work are depicted in Figure 1 by the inputs labeled $\mathbf{r}$. The subscripts denote the loop for which the excitation is intended: ${\mathbf{r}}_{ba}$ is for the bare-airframe dynamics; ${\mathbf{r}}_{cl}$ is for the closed-loop dynamics; and ${\mathbf{r}}_{mb}$ and ${\mathbf{r}}_{sb}$ are for the broken-loop dynamics opened at the mixer input and sensor output, respectively. Figure 1 also shows the relevant loop inputs $\mathbf{u}$ and outputs $\mathbf{y}$ with a similar scheme for the color and subscripts. This information is summarized in

Table 1. A summary of the excitation, input, and output signals for different control loops.

Loop/Dynamics	Excitation	Input Description	Output Description
Bare airframe	${\mathbf{r}}_{ba}$	Control effector positions, ${\mathbf{u}}_{ba}$	Sensor measurements, ${\mathbf{y}}_{ba}$
Closed loop	${\mathbf{r}}_{cl}$	FCS inputs from pilot station, ${\mathbf{u}}_{cl}$	Sensor measurements, ${\mathbf{y}}_{cl}$
Mixer broken loop	${\mathbf{r}}_{mb}$	Mixer inputs, ${\mathbf{u}}_{mb}$	FCS output commands, ${\mathbf{y}}_{mb}$
Sensor broken loop	${\mathbf{r}}_{sb}$	FCS inputs from feedback, ${\mathbf{u}}_{sb}$	Sensor measurements, ${\mathbf{y}}_{sb}$

. Note that the sensor outputs were used for the identification of several loops.

A major issue with identification of aircraft dynamics is the correlation of signals [2,5]. In parametric modeling, inputs must exhibit low correlation coefficients, whereas with frequency response estimation, inputs must also exhibit low correlations in frequency. This is one reason that traditional testing is performed using single inputs. As is discussed in the following sections, testing can be performed with multiple inputs for identifying the frequency responses of multiple control loops by using (1) multiple excitation points, (2) orthogonal inputs at different harmonic frequencies, and (3) appropriate methods for frequency response estimation.

2.2. Multisine Excitation Inputs

The excitation inputs $\mathbf{r}$ used in this work were sets of orthogonal phase-optimized multisines. These multisines were introduced in refs. [27,28] and are summarized in this section. More details on the design and implementation of multisine inputs are given in refs. [2,29,30,31]. Notable applications of multisines have included hypersonic aircraft [28], off-nominal flight conditions [29], aeroelastic aircraft [17,19,20,32], pilot tracking tasks [33], rotorcraft and multicopters [34,35], and others.

Each of the multisine inputs has the form:

$$r_j\left(t\right)=a_j\sum _{k\in K_j}\sqrt{{\mathrm{\Phi }}_k}sin\left({\omega }_{k}t+{\varphi }_k\right)$$

(1)

for excitations $j=1,2,\cdots ,n_r$. The design of a set of multisines typically begins by first picking the time duration T over which the excitation is applied. This selection defines the fundamental frequency of the excitation ($1/T$ Hz) and the associated harmonic frequencies (in rad/s):

$${\omega }_k=\frac{2\pi k}{T}$$

(2)

for each harmonic k. The time duration is typically dictated by the constraints of the flight space, the time afforded to a given maneuver or the lowest frequency of interest. Consequently, there is a trade-off between frequency resolution, the number of cycles observed at a particular frequency (and hence, the accuracy of the frequency response estimate [36,37]), and the number of inputs that can be excited during a maneuver. In Equation (1), $a_j$ is the overall gain applied to the multisine, ${\mathrm{\Phi }}_k$ is the normalized power spectrum, and ${\varphi }_k$ is the phase angle spectrum, which are discussed below.

The harmonics k that are included in the multisines are selected to cover the bandwidth(s) of interest for the identification, as per Equation (2). Harmonics pertaining to frequencies outside the bandwidth of interest are simply skipped during the summation in Equation (1), so that the excitations only contain power within a specific bandwidth. Furthermore, the power is only at the known and discrete harmonic frequencies, which is in contrast to the continuous frequency content characteristic of frequency sweeps and other inputs. An approximate logarithmic spacing, so that resulting Bode plots appear evenly drawn over frequency, can be chosen if desired [37].

When designing multisines for exciting multiple inputs, the harmonics are collected into the $n_r$ sets $K_j$. This is typically carried out by interleaving the harmonics amongst the inputs in an alternating manner so that each input contains a wide band of excitation. For example, the harmonics $k=4,7,10,\cdots$ may be assigned to the elevator, while $k=5,8,11,\cdots$ are assigned to the aileron and $k=6,9,12,\cdots$ to the rudder. This strategy provides an even distribution of frequencies for each input; however, this is not necessary for any of the analyses discussed in this paper.

Multisines are mutually orthogonal, which is a key property that is needed for the frequency response estimation of multiple loops and multiple axes. As explained in refs. [2,27,28], taking the inner product of any two sinusoids with harmonics $k_1\ne k_2$ over the excitation duration evaluates as:

$$\begin{array}{}\hfill & {\int }_0^T{a}_{1}sin\left(\frac{2\pi k_1}{T}t+{\varphi }_1\right)a_{2}sin\left(\frac{2\pi k_2}{T}t+{\varphi }_2\right)dt\hfill \hfill & =\frac{a_1{a}_2}{2}{\int }_0^{T}cos\left[\frac{2\pi (k_1-k_2)}{T}t+{\varphi }_1-{\varphi }_2\right]-cos\left[\frac{2\pi (k_1+k_2)}{T}t+{\varphi }_1+{\varphi }_2\right]dt\hfill & \mathrm{(3a)}\hfill & =\frac{a_1{a}_2}{2}\left\{\frac{T}{2\pi (k_1-k_2)}{\left[sin\left(\frac{2\pi (k_1-k_2)}{T}t+{\varphi }_1-{\varphi }_2\right)\right]}_0^T\right.\hfill \hfill & \left.-\frac{T}{2\pi (k_1+k_2)}{\left[sin\left(\frac{2\pi (k_1+k_2)}{T}t+{\varphi }_1+{\varphi }_2\right)\right]}_0^T\right\}\hfill & \mathrm{(3b)}\hfill & =0\hfill & \mathrm{(3c)}\end{array}$$

regardless of the amplitudes or phase angles. This means that any two harmonic sinusoids are orthogonal at the end of the excitation. This property extends to the multisine excitations as well because multisines are linear combinations of harmonic sinusoids, which is why the harmonics can be arbitrarily and uniquely distributed to different inputs. Furthermore, this result holds in both the time and frequency domains due to the linearity of the Fourier transform operator. Reference [30] showed that multisines are also orthogonal at times other than T, depending on the frequencies used in each multisine, and have pairwise correlations that decrease over time. This orthogonal property of the multisines was leveraged for the analyses discussed in the following sections for distinguishing the input–output pairs amongst the different control loops. Orthogonality has also been exploited in other fields for efficient testing, such as with orthogonal square waves [38,39]. However, sinusoids were used here, instead of square waves, to achieve the desired frequency bandwidth, encourage the development of the steady-state solution after the transients had decayed, and comply with the theoretical definition of the frequency response.

After selecting the time duration and allocating harmonics over the bandwidth of interest to potentially multiple inputs, the next step is to design the normalized power spectra for each input. The power spectrum can be arbitrarily chosen to amplify or diminish the relative amplitude of each excitation frequency, subject to the constraint that for each input:

$$\sum _{k\in K_j}{\mathrm{\Phi }}_k=1$$

(4)

In the simplest approach, the uniform distribution:

$${\mathrm{\Phi }}_k=\frac{1}{n_{f_j}}$$

(5)

typically works well to robustly cover a broad range of frequencies [2,40,41]. This spectrum also has the advantage that a single power spectrum can be designed for an entire flight test campaign, over which the aircraft modes and resonances can significantly shift with the flight condition [20]. Alternatively, a uniform distribution can be the starting point for further refinement and optimization when the models of interest are well known. Power can be reduced to diminish the level of response over a certain bandwidth, similar to using a notch filter but without the associated phase shifts. Another example is limiting the low-frequency elevator power to reduce changes in airspeed and dynamic pressure and therefore, flight condition.

Figure 2, which is adapted from ref. [29], shows an example design using different power spectra. This example was for AirSTAR flight tests with the T-2 airplane. These multisines were designed for a 10 s duration with power of $0.2$ to $2.2$ Hz. The aileron input (shown in red) was provided a uniform distribution, whereas the elevator and rudder inputs (shown in blue and orange, respectively) had increased power of around 1 Hz, which was where the short period and dutch roll modes were expected to reside. For other applications, power spectra could also be tailored to mimic the effects of measurement noise or turbulence (standard models or empirically determined), as in refs. [42,43].

In general, when sinusoids are combined into a multisine input, the phase angles can create large peaks in the multisine. For aircraft, such inputs could take the aircraft off the intended flight condition and reduce the accuracy of identified models based on linear modeling techniques, such as frequency responses. To reduce these large peak amplitudes, phase angles are sometimes randomly drawn until a suitable waveform is attained. Alternatively, ref. [44] proposed an analytical technique (now called Schroeder sweeps) for selecting phase angles that result in multisines with relatively low peak factors. In a similar vein, direct numerical searches are sometimes performed to find phase angles that minimize a metric related to the size of the peak amplitudes of the resulting multisine. For example, ref. [37] provides two algorithms for crest factor minimization. The multisines in this work were optimized for the minimum relative peak factor [2,27,28,29,30,31]:

$$\mathrm{RPF}\left[r_j\left(t\right)\right]=\frac{\max \left[r_j\left(t\right)\right]-\min \left[r_j\left(t\right)\right]}{2\sqrt{2}\mathrm{rms}\left[r_j\left(t\right)\right]}$$

(6)

which quantifies the efficiency of an input. The theoretical minimum RPF value is unity, and larger RPF values reflect larger peaks in the input. For example, the input design in Figure 2 (adapted from ref. [29]) resulted in RPF values of $1.03$, $1.14$, and $1.15$.

After the phase angles have been optimized for a multisine, the individual sinusoids are collectively shifted in time until the resulting waveform begins and ends at zero. This process removes an effective step response from the beginning and end of the multisine, which may excite transient responses, frequencies outside the bandwidth of interest, and nonlinearities, which all decrease the accuracy of the results analyzed using frequency responses.

The process of constructing a set of multisine inputs from specified design parameters is implemented in the SIDPAC function mkmsswp.m [2,23]. The optimization of the phase angles was performed using the Nelder–Mead simplex search method.

An example of the effect of phase angle optimization for the elevator power spectrum of Figure 2 is shown in Figure 3. Four different inputs are shown, which correspond to phase angles selected as 0 deg (drawn in blue), random phase angles (drawn in red), Schroeder phase angles (drawn in purple), and phase angles optimized using the SIDPAC function (drawn in green). For the waveforms shown, the optimized input had the lowest RPF value of 1.01.

Lastly, the overall gains $a_j$ are designed. For aircraft, this step is often carried out in two parts. First, each of the multisines are scaled relative to one another. For example, the elevator and rudder multisines could be twice as large as the antisymmetric aileron multisine. This is to account for the relatively small roll inertia that is typical of conventional aircraft. Sometimes the sign of the amplitude for a particular multisine is switched to force an initial nose-down moment to avoid reaching high angles of attack. Second, the family of all multisines are collectively scaled to meet the desired signal-to-noise ratios. Past flight tests at NASA were planned with several discrete amplitudes (corresponding to “small,” “medium,” and “large”) at each flight condition. During the flights, the amplitudes were increased from small to large until sufficient signal-to-noise ratios and modeling results were obtained.

2.3. Frequency Transforms

When estimating frequency responses, the measured data are transformed into the frequency domain for further processing. Analytically, this can be carried out using the finite–time Fourier transform, which is defined for a generic signal $z\left(t\right)$ as:

$$z\left(j\omega \right)={\int }_0^{T}z\left(t\right)e^{-j\omega t}dt$$

(7)

As shown in Equation (7), the argument of a variable is used in this paper to distinguish between the time history and Fourier transform. In practice, the Fourier transform is implemented numerically with measured data. One approach for doing this is the chirp Z-transform. Unlike the fast Fourier transform (FFT), the chirp Z-transform can be implemented with an arbitrary frequency resolution that is independent of the data record length. The SIDPAC function fint.m implements the chirp Z-transform with an option to apply cubic interpolation to improve the accuracy for deterministic signals [2,23,45]. More discussion on the Fourier transforms of flight data is provided in refs. [2,46].

For frequency response estimation when multisines are used, Fourier transforms only need to be computed at the harmonic frequencies ${\omega }_k$ in Equation (2). This is because for linear systems, which are assumed in frequency response analyses, the steady-state responses only contain power at frequencies in the input. As a result, only a relatively small number of data points in the frequency domain need to be processed, which saves computation time and effort. These data points, however, have relatively high signal-to-noise ratios because the excitation inputs have power at precisely those frequencies. This is again in contrast to frequency sweeps and other inputs, in which the power and signal-to-noise ratios vary with frequency.

Before the measured data are transformed into the frequency domain, trends are removed from the data. This helps to minimize spectral leakage [2,46], particularly for signals that have non-zero steady-state values, such as the angle of attack or vertical accelerometer measurements. This can sometimes be accomplished by subtracting the first point or the average of the first few samples of a time series or by removing the line of best fit from the data. Other approaches for doing this, such as removing a higher-order polynomial, are discussed in ref. [46].

2.4. Direct Approach for Estimating Frequency Responses

A frequency response quantifies the steady-state response of a linear time-invariant system to a sinusoidal input when all initial conditions are zero. Using multisine inputs and Fourier transforms, as described in Section 2.2 and Section 2.3, an estimate of the frequency response from each measured input $u_j$ to each measured output $y_i$ can be computed by this definition as:

$${\widehat{G}}_{ij}\left(j{\omega }_k\right)=\frac{y_i\left(j{\omega }_k\right)}{u_j\left(j{\omega }_k\right)}$$

(8)

where $i=1,2,\cdots ,n_y$ and $k\in K_j$ [2,46,47,48]. Note that, similar to the Fourier transforms, the frequency responses are only computed at the harmonic frequencies of the respective input because those are the only frequencies at which the input–output data have sufficient information. In matrix/vector notation, where the inputs and outputs are vectors and the frequency response is a matrix, the scalar version in Equation (8) is rewritten as:

$$\widehat{\mathbf{G}}\left(j{\omega }_k\right)=\frac{\mathbf{y}\left(j{\omega }_k\right)}{\mathbf{u}\left(j{\omega }_k\right)}$$

(9)

The estimates can be displayed in a variety of ways. In many applications, a Bode plot is desired in which the magnitude and phase components:

$$\begin{array}{cc}\hfill \parallel G_{ij}\left(j{\omega }_k\right)\parallel & =20{\log }_{10}\sqrt{\Re {\left[G_{ij}\left(j{\omega }_k\right)\right]}^2+\Im {\left[G_{ij}\left(j{\omega }_k\right)\right]}^2}\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill \angle G_{ij}\left(j{\omega }_k\right)& =\frac{180}{\pi }arctan\left(\frac{\Im \left[G_{ij}\left(j{\omega }_k\right)\right]}{\Re \left[G_{ij}\left(j{\omega }_k\right)\right]}\right)\hfill \end{array}$$

(10b)

in units of dB and deg, respectively, are plotted against radian frequency on a logarithmic scale. For examining stability margins, the estimated frequency responses are often plotted as a Nichols chart with magnitude in dB against phase in deg. In other applications, results are presented as a polar plot or Nyquist diagram with the real part vs. the imaginary part of the frequency response.

In general, Equation (9) is an unbiased estimate of the frequency response, which has a variance that is inversely proportional to the squared signal-to-noise ratio for steady-state data over the integer number of cycles [36]. As shown in ref. [49], Gaussian noise on the output measurements creates a proportional Gaussian distribution in the real and imaginary parts of the frequency response, which is important when fitting parametric models to frequency responses. Nonlinearities, unmodeled inputs, time-varying dynamics, and transient responses degrade the accuracy of frequency response estimates.

The approach outlined in this section for frequency response estimation using multisine inputs is analogous to the original technique of sine dwell testing, in which a single input was oscillated at a single frequency, and magnitude and phase differences could be determined from the steady-state response data. Running an experiment with multisines is similar, except many frequencies are used on many inputs at the same time, which effectively runs many tests concurrently to enable savings in time and costs. This method is also similar to the Fourier approach, in which the frequency response is computed from the ratios of Fourier transform data. Again, multisines are used here, rather than pulses, steps, and doublets, so that steady-state modeling data were achieved and the analysis was conducted at a relatively low number of harmonic frequencies with high signal-to-noise ratios. This removed engineering judgement on the part of the analysis and could enable real-time estimation, as discussed later in Section 2.6.

The direct frequency response estimates from Equation (9) are applicable in three situations. The first case is for identifying open-loop systems, such as an airplane flying without feedback control. The second case is for identifying input–output frequency responses for closed-loop systems, such as handling qualities analysis. The third case is for identifying the frequency response of a broken loop with all other loops closed, as performed to compute stability margins.

Figure 4 shows three examples in which this approach was recently used at NASA. The first example, shown in Figure 4a, displays the identified bare-airframe frequency response, from elevator to vertical accelerometer output, as a Bode plot for the T-2 airplane used in the NASA AirSTAR flight tests. For this maneuver, which was performed during Flight 33, the multisines from Figure 2 were simultaneously applied to the elevator, ailerons, and rudder with no feedback control. The black dots are the frequency response estimates obtained from the direct approach, as presented in Equation (9). The solid red lines correspond to the frequency response of a parametric model identified from time-domain data using the output-error approach [2,23]. Other input–output frequency responses from these data, not shown here, can be found in ref. [49].

The second example, shown in Figure 4b, was taken from ref. [50]. Frequency responses were used to efficiently determine a reduced-order model for the computational fluid dynamics (CFD) solver analysis of a half-span flexible wind tunnel article called the Integrated Adaptive Wing Technology Maturation (IAWTM). In total, 196 frequency responses were identified between the 14 inputs (11 prescribed structural mode displacements and 3 control surface deflections) and the 14 outputs (the corresponding generalized aerodynamic forces (GAFs)) over the frequency range of $0.25$ to 80 Hz using this approach. Using multisines significantly reduced the computation time needed to compute the full frequency response matrix compared to using single-input methods. The single frequency response shown in Figure 4b corresponds to the GAF for the first structural mode output resulting from forcing the second structural mode input. The Mach number for these data was $0.8$, the dynamic pressure was 276 lbf/ft ${}^2$, and some nonlinearities in the flow were present, which distorted the frequency response estimate. The frequency response is shown as a polar plot, as is typical in the aeroelasticity literature. The black dots are again the frequency response estimates obtained using the direct approach in Equation (9). The colored lines are the other estimates, discussed further in ref. [50], which were obtained using the ZAERO^™ panel code software (shown in blue) and the linearized frequency domain (LFD) approach (shown in red). The remaining 195 frequency responses from the analysis are presented in ref. [50].

The third example, shown in Figure 4c, is a real-time screenshot from the control room at the NASA Armstrong Flight Research Center (AFRC) during the testing of the X-56A Multi-Utility Technology Testbed (MUTT) aeroelastic demonstrator. This Bode plot shows the broken-loop frequency response for the roll axis loop, which was opened at the mixer input, using a classical take-off and landing control law during Flight 8. The black triangles are the estimates obtained using telemetered flight data and Equation (9), and the blue lines are predictions based on the pre-flight models. The 16 dB gain margin was identified during the flight and corroborated with pre-flight models, which contributed to the efficient progression to the next scheduled test point.

2.5. Joint Input–Output Approach for Estimating Frequency Responses

In many cases, it is desirable to identify an open-loop system from closed-loop data. This is the case when identifying the bare-airframe dynamics when the aircraft is flown with feedback, particularly when the aircraft is nominally unstable.

When applying the direct approach to this problem, incorrect results are obtained [51,52]. Although the multisines are orthogonal and distinct harmonic frequencies are applied to the input command path, the feedback puts all of the harmonic frequencies on all of the inputs that include the feedback. This correlates the inputs in frequency, although pairwise correlation coefficients remain low. Applying Equation (9) ignores any additional off-axis terms that are not negligible for this case. This problem does not impact the results of parametric modeling approaches, however, such as equation error and output error [2,36,52].

For low correlations in frequency, the direct approach can sometimes still be used to estimate open-loop frequency responses from closed-loop data, albeit with some error. For larger amounts of correlation, an indirect identification can be used wherein frequency responses are identified for the closed-loop system and then transfer function algebra and a model for the control law are used to solve for the frequency responses of interest [36].

Another approach, which also works with high amounts of correlation, is called the joint input–output (JIO) method [36]. Two intermediate frequency responses are estimated using the direct approach: one from the excitation to the system input and another from the excitation to the system output. The frequency responses of interest are then computed as:

$$\widehat{\mathbf{G}}\left(j{\omega }_k\right)=\frac{\mathbf{y}\left(j{\omega }_k\right)}{\mathbf{u}\left(j{\omega }_k\right)}=\left[\frac{\mathbf{y}\left(j{\omega }_k\right)}{\mathbf{r}\left(j{\omega }_k\right)}\right]{\left[\frac{\mathbf{u}\left(j{\omega }_k\right)}{\mathbf{r}\left(j{\omega }_k\right)}\right]}^{-1}$$

(11)

Note that the number of excitations must equal the number of system inputs to be able to compute the matrix inverse in Equation (11). This process produces an unbiased estimate of the frequency response for an open-loop system from closed-loop data. In the case of open-loop data, the frequency response from the excitation to the input is approximately unity and the direct approach in Equation (9) is recovered.

Reference [51] suggested using the JIO approach for aircraft system identification with highly correlated inputs and demonstrated results from flight test data for the Calspan Variable Stability System (VSS) Learjet LJ-25D using frequency sweeps. It was also suggested that the JIO approach could be used with multisine inputs when interpolation is applied to obtain frequency responses at all harmonic frequencies and not just at the harmonics in the corresponding excitation. Otherwise, the matrix inversion in Equation (11) cannot be performed. A similar approach was developed and demonstrated in refs. [52,53].

Using the JIO method is more complicated than the direct method, but it is necessary in some instances. As it depends on two sets of direct frequency response estimates, it requires more computation and is more sensitive to errors, particularly in $\mathbf{u}\left(j\omega \right)/\mathbf{r}\left(j\omega \right)$. However, when using multisines, these increased computations are not usually significant. As interpolation is used, the multisines must also be designed with sufficient frequency resolution when estimating lightly damped bare-airframe modes, as discussed in ref. [52].

The JIO approach is demonstrated here using multisine inputs and linear interpolation with two examples. In the first example, a simulation model of the T-2 airplane was excited using two multisines on the inner and outer elevator pairs, from which the pitch rate measurement was fed back with a proportional gain to the inner elevator pair, as in refs. [52,53]. This feedback resulted in a pairwise correlation coefficient of $0.4$ between the inner and outer elevator deflections. Figure 5a shows the resulting outer elevator to pitch rate frequency response. The solid markers correspond to the JIO estimates, computed with Equation (11), and match the true system dynamics, which are drawn with a solid red line. The open markers were computed with the direct approach in Equation (9) and were incorrect due to the correlated elevator deflections.

In the second example, five multisines were used to excite the 10 control surfaces on the X-56A in symmetric pairs while a feedback control system and mixer were active. Each multisine excitation contained 25 harmonics between $0.2$ and $6.25$ Hz. Figure 5b shows the frequency response estimates from the symmetric outermost control surface pairs to a pitch rate gyroscope located near the aircraft nose, taken from a maneuver from Flight 11 at NASA AFRC. The flight condition for this was straight and level flight at 48% fuel and 63% flutter speed. The 5.8 rad/s resonance was the short period mode and the 18.7 rad/s resonance was the first symmetric wing bending mode. The solid dots are again the frequency response estimates obtained using the JIO method, whereas the open dots are the results obtained using the direct method. The solid line is the frequency response of a parametric model that was estimated using output error with Fourier transform data [19]. The output error and JIO estimates agree. However, at frequencies below about 20 rad/s, where the feedback control law is operating, the results from the direct method differ.

2.6. Real-Time Estimation

The estimation approaches discussed in Section 2.4 and Section 2.5 were presented for batch computation after all of the data from a maneuver were recorded and processed. Many times, however, it can be beneficial to process the data in real time. Some examples include monitoring stability margins during the flight of new and unstable aircraft, efficient envelope expansion flight tests, and fault detection. The frequency response methods can be altered to enable real-time computation.

The Fourier transform, discussed in Section 2.3, can also be implemented by the summation:

$$z\left(j\omega \right)=\Delta t\sum _{i=1}^{N}z\left(t_i\right)e^{-j\omega t_i}$$

(12)

where $\Delta t$ is the sampling period, $t_i$ is a discrete representation of the time index, and N is the total number of samples collected. Equation (12) is an Euler approximation to the Fourier integral in Equation (7) and is sufficient for low frequencies compared to the Nyquist frequency.

By expanding the summation in Equation (12), the Fourier transform of the measured signal up to any arbitrary time $t_i$ can be written as the recursive update [2,46]:

$$z_i\left(j\omega \right)=z_{i-1}\left(j\omega \right)+z\left(t_i\right)e^{-j\omega t_i}\Delta t$$

(13)

starting from the initial condition:

$$z_0\left(j\omega \right)=0$$

(14)

In this way, the Fourier transforms are updated after each measurement. Relatively few computations are needed to perform this operation, several of which may be precomputed. Computational demands can be further decreased by updating Fourier transforms at a lower rate. A forgetting factor or sliding window can also be used to remove data from memory and adapt more quickly to changing conditions [2,40,46].

Afterward, the frequency responses are updated using either the direct method or the JIO method. As these are typically used for displays or to compute stability margins in a real-time setting, frequency response estimation calculations are usually performed at a slower rate than the Fourier transforms.

The real-time estimation of frequency responses is enabled by the combination of multisine inputs, recursive Fourier transform data, and the ratios of the Fourier transform data. Other approaches, such as those based on spectral density estimates, require the entire duration of data and cannot generally be implemented in real time during a maneuver.

2.7. Multiple Loop Estimation

The examples shown so far have illustrated the use of multiple inputs but only estimating frequency responses for a single loop. Such cases are more efficient than using single-input maneuvers to obtain the same information. The efficiency of testing can be further improved by applying sets of multisine inputs at different points throughout the aircraft system, as shown in the block diagram of Figure 1, and then proceeding with the same analysis.

The drawback of simultaneously estimating multiple loops is that for a given time duration, the frequency resolution on any particular frequency response becomes more coarse. To compensate, the multisine time duration should be increased until all estimated frequency responses have a sufficient frequency resolution.

Multiple control loops were identified from flight test data of the X-56A for computing real-time stability margins in flight. One set of maneuvers included a multisine being applied to the elevator command before the mixer and a second multisine being applied to the pitch rate measurement from the forward gyroscope. Figure 6 shows the resulting Bode plots for a maneuver from Flight 15 at NASA AFRC. The numbers have been removed from these plots because these data are restricted by ITAR. The results for the elevator loop, shown in Figure 6a, suffered from low signal-to-noise ratios at higher frequencies. This is evident from the low values and high scatter in the magnitude plot after the roll-off, as well as the scattered phase angle estimates at those frequencies. However, the sensor loop had adequate signal-to-noise, as shown in Figure 6b. Nonetheless, these results matched well enough to the pre-flight predictions (not shown) and were useful in establishing sufficient margins to continue flight testing and envelope expansion.

3. Example

In this section, an example is presented for identifying multiple control loops from a single maneuver using frequency responses and multisine inputs injected at different points throughout the aircraft system. This demonstration involved a simulation model of the X-59 Quiet SuperSonic Technology (QueSST) vehicle that is used in the NASA Low-Boom Flight Demonstrator (LBFD) project [54]. The aircraft was unstable in pitch and required feedback control, which correlated the signals. Frequency responses were identified for the longitudinal bare-airframe and closed-loop dynamics, as well as the broken-loop paths that were opened at the mixer input and the pitch rate sensor to estimate stability margins.

The X-59 airplane was designed to create a significantly quieter sonic boom. Figure 7 shows a rendering of the aircraft, which has a 99.6 ft length and 29.5 ft wingspan. Powered by a single engine, the X-59 was designed for research at a nominal Mach 1.4 operating condition. The aircraft has eight independent control surfaces: left and right stabilators; left and right ailerons; left and right flaps; rudder; and the “T-tail” horizontal surface on the top of the tail section. For primary flight control, the stabilators are moved symmetrically and the ailerons are moved anti-symmetrically. The canards mounted on the sides of the fuselage just in front of the pilot station are fixed.

A nonlinear, six degree of freedom flight dynamics simulation was constructed in MATLAB^® and Simulink^® by Lockheed Martin with input from NASA. The aerodynamic database was based on wind tunnel and CFD data. For the results shown below, this nonlinear simulation was numerically linearized at a flight condition with low altitude, slow speed, and medium weight in an up-and-away configuration.

Table 2. The geometry and mass properties of the X-59 in flight condition 7220.

Parameter	Value	Unit
b	$\mathrm{29.50}$	ft
$\overline{c}$	$21.96$	ft
S	$486.16$	ft ${}^2$
m	$637.16$	slug
$I_{xx}$	9438	slug - ft ${}^2$
$I_{yy}$	153,650	slug - ft ${}^2$
$I_{zz}$	160,110	slug - ft ${}^2$
$I_{xy}$	$-87$	slug - ft ${}^2$
$I_{xz}$	4296	slug - ft ${}^2$
$I_{yz}$	7	slug - ft ${}^2$
$x_{cm}$	841	in

presents the geometry and mass distribution parameters for this flight condition, and

Table 3. The trimmed flight variables of the X-59 in flight condition 7220.

Parameter	Value	Unit
h	20,000	ft
$\alpha$	$3.8$	deg
V	509	ft/s
KEAS	220	knots
M	$0.45$	—
${\delta }_{placmd}$	54	deg
${\delta }_{a_l}$	$1.7$	deg
${\delta }_{a_r}$	$1.3$	deg
${\delta }_s$	$5.8$	deg
${\delta }_r$	$-0.3$	deg
${\delta }_f$	$2.5$	deg

lists the trimmed flight variables.

Table 4. The nondimensional stability and control derivatives of the X-59 longitudinal bare-airframe dynamics in flight condition 7220.

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
$C_{D_V}$	$0.065$	$C_{D_{\alpha }}$	$0.28$	$C_{D_q}$	$0.081$	$C_{D_{{\delta }_s}}$	$0.12$
$C_{L_V}$	$0.50$	$C_{L_{\alpha }}$	$2.5$	$C_{L_q}$	$2.4$	$C_{L_{{\delta }_s}}$	$0.50$
$C_{m_V}$	$-0.0032$	$C_{m_{\alpha }}$	$0.52$	$C_{m_q}$	$-2.3$	$C_{m_{{\delta }_s}}$	$-0.42$

gives the longitudinal stability and control derivatives in nondimensional form for the flight condition. These derivatives were defined [2] as:

$$\begin{array}{cccccccc}\hfill C_{D_V}& =\frac{\partial C_D}{\partial \frac{V}{V_0}}\hfill & \hfill C_{D_{\alpha }}& =\frac{\partial C_D}{\partial \alpha }\hfill & \hfill C_{D_q}& =\frac{\partial C_D}{\partial \frac{q\overline{c}}{2V_0}}\hfill & \hfill C_{D_{{\delta }_s}}& =\frac{\partial C_D}{\partial {\delta }_s}\hfill \end{array}$$

(15a)

$$\begin{array}{cccccccc}\hfill C_{L_V}& =\frac{\partial C_L}{\partial \frac{V}{V_0}}\hfill & \hfill C_{L_{\alpha }}& =\frac{\partial C_L}{\partial \alpha }\hfill & \hfill C_{L_q}& =\frac{\partial C_L}{\partial \frac{q\overline{c}}{2V_0}}\hfill & \hfill C_{L_{{\delta }_s}}& =\frac{\partial C_L}{\partial {\delta }_s}\hfill \end{array}$$

(15b)

$$\begin{array}{cccccccc}\hfill C_{m_V}& =\frac{\partial C_m}{\partial \frac{V}{V_0}}\hfill & \hfill C_{m_{\alpha }}& =\frac{\partial C_m}{\partial \alpha }\hfill & \hfill C_{m_q}& =\frac{\partial C_m}{\partial \frac{q\overline{c}}{2V_0}}\hfill & \hfill C_{m_{{\delta }_s}}& =\frac{\partial C_m}{\partial {\delta }_s}\hfill \end{array}$$

(15c)

and were numerically evaluated at the reference condition. At the flight condition, the bare-airframe model was significantly unstable in the pitch axes due to the positive value of $C_{m_{\alpha }}$ and it had the poles $[0.23,0.092]$, $(-1.9)$, and $\left(3.0\right)$. These poles are written in the shorthand notation of ref. [1], in which first- and second-order polynomial factors are represented by $s-a\to \left(a\right)$ and $s^2+2\zeta {\omega }_{n}s+{\omega }_n^2\to [\zeta ,{\omega }_n]$, respectively. Although the first mode was characteristic of a classical oscillator phugoid mode, the complex pole pair that typically forms the short period mode was attributed to the two real poles: one stable and one unstable. This combination of poles is common for fighter aircraft with relaxed static stability in the longitudinal axis.

Four multisine excitation inputs were designed. As shown in Figure 1, these multisines were added at the pilot longitudinal stick input, pitch command at the mixer input, stabilator command at the actuator input, and pitch rate gyroscope at the sensor output. Additional multisines could have been simultaneously added to inputs for the other axes (such as the pilot lateral stick, roll command, aileron commands, and roll rate gyroscope) and the analysis would have proceeded in the same fashion. However, only these four inputs were used to keep the example relatively simple.

Summarized parameters for these multisine designs are listed in

Table 5. A summary of multisine design parameters for the X-59 example ($T=60$ s).

Excitation	Harmonics	Frequencies, Hz	Amplitude	RPF
${\mathbf{r}}_{ba}$	$6,10,\cdots ,118$	$0.1000$ to $1.9667$	1 deg	$1.14$
${\mathbf{r}}_{cl}$	$7,11,\cdots ,95$	$0.1167$ to $1.5833$	7.5%	$1.21$
${\mathbf{r}}_{mb}$	$4,8,\cdots ,240$	$0.0667$ to $4.0000$	1 deg	$1.16$
${\mathbf{r}}_{sb}$	$5,9,\cdots ,237$	$0.0833$ to $3.9500$	1 deg/s	$1.37$

. The inputs were designed for 60 s of excitation, which resulted in a fundamental frequency and frequency resolution of $0.0167$ Hz. For reference, the frequency sweeps applied to the X-29 airplane to estimate the longitudinal margins (broken at the mixer input) were also roughly 60 s in duration [14]. Each multisine was designed with a uniform power spectrum for simplicity. The actuator multisine was designed to span the bare-airframe modes, whereas the pilot multisine was selected to cover the frequency range that is typical of handling qualities analyses [25] and the mixer and sensor multisines were selected to cover the ranges of the crossover frequencies. The phase angles were optimized for the minimum RPF using the Nelder–Mead simplex optimization. The amplitudes were selected to be small but to also produce adequate signal-to-noise ratios for the frequency response modeling.

Figure 8 shows the time histories of the multisine excitation inputs, which are colored to match the inputs in Figure 1. An extra 5 s of time without any excitation was added to the beginning and end of the maneuver.

The excitations were applied to the closed-loop simulation and the data were sampled at 100 Hz. Measurement noise sequences were added to all sensor measurements. Figure 9 and Figure 10 show the inputs and outputs that were needed to identify the bare-airframe, closed-loop system, and broken-loop frequency responses. These are again color coded to match Figure 1. The only exception is the second output, which was common for all of the frequency responses except for the open-loop stabilator and is shown in black. Evidence of the measurement noise was present in all signals except for $u_{cl}$. Note that the designed excitations kept the inputs and outputs relatively small in amplitude but with good signal-to-noise ratios, the lowest being $7.6$ for the FCS output $y_{mb}$.

The time histories were transformed into the frequency domain using the high-accuracy Fourier transform mentioned in Section 2.3. The bare-airframe frequency response was computed using the JIO approach discussed in Section 2.5 with linear interpolation over the excitation frequencies. The closed-loop and broken-loop frequency responses were computed using the direct method from Section 2.4 at the appropriate frequencies listed in Table 5.

Figure 11 shows the resulting frequency responses as Bode plots. The black lines indicate the true frequency responses of the X-59 linear simulation model, whereas the colored dots indicate the estimates from the entire record of data and are color coded to match Figure 1. The estimates fitted the true frequency responses very closely. Some sensitivity to the measurement noise can be observed in the stabilator broken-loop results in Figure 11c at the lowest and highest frequencies. This is due to the lower signal-to-noise ratio of $y_{mb}$ and is similar to the experimental results seen during the X-56A flight testing presented in Figure 6a. Although these errors were small, results could be improved by increasing the excitation duration or amplitude or by increasing the power at these frequencies.

Figure 12 shows the real-time estimates of the stability margins for the stabilator broken loop. The black lines are the true margins for the linear model and the dots are the estimates. In this case, the recursive Fourier transform described in Section 2.3 was applied. Fourier transforms were updated with each new measurement at 100 Hz and frequency responses and stability margins were calculated at 1 Hz once the excitation began. The trim values used to detrend the time histories were taken as the measurement at the time that the excitation started. On a typical laptop computer running MATLAB^®, computing the Fourier transforms used about 20% of the period associated with the 100 Hz frame rate. When the frequency responses were updated and the stability margins calculated, about 70% of the 100 Hz frame was used. Although there were two gain margins ($-7.4$ dB at $1.4$ rad/s and $6.6$ dB at $16.1$ rad/s) and one phase margin (48 deg at $5.8$ rad/s) for the true simulation model, several margins were sometimes calculated, particularly near the start of the maneuver, because of noise and low information in the evolving frequency response estimates. However, after about 20 s into the maneuver, all stability margins were closely matched.

4. Conclusions

This paper presented the latest advances from NASA in efficiently estimating frequency responses for different aircraft control loops from flight test data using orthogonal phase-optimized multisine inputs, which were added at multiple points throughout the aircraft system. Relatively simple techniques were used to compute Fourier transforms and frequency responses. Several examples were shown, including flight test data, CFD experiments, and a demonstration using the X-59 simulation model.

The main points of this paper can be summarized as the following:

1.

This approach can be used to estimate the bare-airframe dynamics from closed-loop data, the closed-loop dynamics, the broken-loop dynamics (opened at the mixer input or the sensor output) or other dynamics of interest;

2.

These frequency response estimates can be computed in real time;

3.

Frequency responses can be estimated for a single loop and multiple axes, multiple loops over a single axis or multiple loops and multiple axes. The simultaneous identification of multiple loops and multiple axes attains the most efficient test results and the most savings in time and costs.