Practical System Identification and Incremental Control Design for a Subscale Fixed-Wing Aircraft

Abstract

An incremental differential proportional integral (iDPI) control law using eigenstructure assignment gain design is tested in flight on a subscale platform to validate its suitability for fixed-wing flight control. A kinematic relation for the aerodynamic side-slip angle rate is developed to apply a pseudo full state feedback. In order to perform the gain design and assessment, a plant model is estimated using flight test data from gyro, accelerometer, airspeed and surface deflection measurements during sine-sweep excitations. Transfer function models for the actuators and surface deflections are identified both in-flight and on the ground for several different actuators and control surfaces using hall sensor surface deflection measurements. The analysis reveals a large variation in bandwidth between the different types of servo motors. Flight test results are presented which demonstrates that the plant model estimates based on tests with good frequency excitation, high bandwidth actuators and surface deflection measurements can be used to reasonably predict the closed-loop dynamic behavior of the aircraft. The closed-loop flight test results of the iDPi control law show good performance and lays the groundwork for further development.

1. Introduction

This paper investigates through flight test the feasibility of using an incremental differential proportional integral (iDPI) control law as introduced in [1] for controlling the lateral-directional motion of a fixed-wing aircraft. The investigation covers identification of servo-motor dynamics, system identification of a plant model, setup of the control law structure, gain design, assessment of stability and robustness, and at last flight test results of the control law. The iDPI control concept was developed from the method of incremental nonlinear dynamic inversion (INDI) [2,3]. The iDPI control law is an implementation of a gain scheduled proportional-integral (PI) controller that avoids hidden coupling terms, and inherits the robust performance of INDI against model uncertainties. Avoiding hidden coupling terms enables the design and assessment of the controller using linear control theory, even when scheduling the controller gains makes the control law non-linear. This paper forms the starting point of the investigation of its usefulness for practical flight control. It provides validation of its usability in a real-world setting and lays the groundwork for further investigation.

System identification for aircraft is a wide topic, and many branches of the field are discussed in [4,5,6]. Recent advances in aircraft system identification include [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. Several methods above use frequency domain identification, which was also studied in [26,27,28,29]. Also, online identification has been used in e.g., [30,31,32,33].

In this paper, surface deflection measurements was used to both identify in-flight actuator dynamics and for identifying a flight dynamics model. The system consists of a hall sensor device near the surface hingeline, as developed and tested in [34,35,36]. In contrast, internal actuator measurements were used in [37,38].

Fixed wing flight control is also a wide field with many different approaches and objectives; an excerpt is discussed in [39,40,41,42,43,44,45,46,47,48,49,50,51]. The control method of this paper is based on incremental methods, which are an active field of research. Incremental control has been proposed to reduce the dependency of model information for control law design and has been studied recently in [52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85]. In this paper, the method of eigenstructure assignment is additionally used for gain design, a method that has been investigated in [86,87,88].

Flight testing of Unmanned Aerial Vehicles (UAVs) has increased with advancements in microcontroller hardware and sensors, and the ease of integration with tools such as MATLAB/Simulink. Recent UAV flight test campaigns include [89,90,91,92,93,94,95,96,97,98,99], which also include flight testing with incremental type control laws in [100,101,102,103,104]. Subscale flight testing and system identification was performed in [105,106,107,108,109]. In [110], both system identification and control law design were performed for the longitudinal motion of a fixed wing UAV. Complete autopilot systems for UAV have been built or described in, for example, [111,112,113,114,115,116] including the PX4 flight stack used in this study.

The goal of this paper is to validate the feasibility of iDPI in a practical setting. In addition, the goal is to obtain an actuator and plant model of a small fixed wing UAV of good enough quality to design and predict the closed-loop behavior with the controller. The contributions of this paper are:

• Frequency domain identification of servo dynamics of different actuators both in flight and on ground using surface deflection measurements.
• Presentation of experimental results for system identification of a subscale fixed-wing aircraft using surface deflection measurements.
• Verification that the model resulting from system identification can accurately predict the closed-loop control law performance.
• Verification that iDPI control laws can obtain good performance on a subscale aircraft.
• Application of eigenstructure assignment for an iDPI lateral/directional control law.
• Derivation of kinematic relations for flow angle rates to avoid derivatives of load factor or flow angle measurements applicable to differential-integral control laws.

2. Experimental Setup

The aircraft, denoted “The Albadrone”, used for flight testing is the X-UAV Mini Talon shown if Figure 1, with modifications based on [117]. The basic data are given in

Table 1. Albadrone basic data.

Wingspan	1.3 m
Length	0.83 m
Weight	1–2 kg
Material	Polystyrene Foam

. The aircraft has a V-tail configuration with pusher propeller. Pitot tubes are mounted one at each wingtip, with only the left one active for the following experiments. For accurate airspeed measurements, wingtip mounted pitot tubes are not ideal, but were set up as such to support another research project investigating the ability to estimate the wind gradient in banked flight. The flight control computer is a Pixhawk 4, with the PX4 flight stack [118]. The control laws are implemented with MATLAB/Simulink [119], and flashed onto the aircraft using the UAV Toolbox Support Package for PX4 Autopilots. The sensor setup is depicted in Figure 2. The inertial measurement unit (IMU) which includes gyros and accelerometers is placed almost in the center of gravity. The surface deflection measurements are sent with $I^{2}C$ to the Pixhawk, and a custom-made PX4 driver is used to relay the measurements to the Simulink Application. Two $I^{2}C$ buses are used to connect the sensors in order to have the required bandwidth and to avoid address clashing. The $I^{2}C$ packets are decoded in a custom-made PX4 application, which forwards the data to the PX4 Simulink Application using a uORB message. The logging is performed using a custom-made Simulink Bus packing block that allows for simple logging of all relevant internal signals stored in a Simulink Bus object. This allows for rapid changes of logged signals with minimum maintenance of the Simulink model. The signals with a maximum of 498 bytes are logged with 100 Hz, and saved to an onboard SD card. All maneuvers for system identification and controller validation are injected automatically via the flight control software, by using execution and mode switches on the remote controller. Two maneuver injection systems were used, a custom made system for rapid development, and a well tested system also used for manned aircraft [120]. Both systems automatically index the maneuvers, to ease the task of post processing and data analysis.

3. Identification of Servo+Surface Dynamics and Equivalent Time Delay

This section identifies a simplified servo+surface deflection dynamic model for ground and in-air operation at a fixed velocity. Different sets of actuators and linkages were tested and the analysis will uncover the different dynamic behavior. The following list is the objectives of the analysis:

• Identify servo+surface transfer function models.
• Compare in-flight versus ground performance.
• Compare servo type, servo-surface linkage type and surface type.

In order to obtain a simple model for control design, a second-order plus time delay model is chosen:

$$G_A\left(s\right)=\frac{K_A{\omega }_{0,A}^2}{s^2+2{\omega }_{0,A}{\zeta }_{A}s+{\omega }_{0,A}^2}{e}^{-T_{eq}s}$$

(1)

where $G_A$ denotes the actuator transfer function, ${\omega }_{0,A}$ is the actuator eigenfrequency, ${\zeta }_A$ is the relative damping, $K_A$ is the actuator static gain and $T_{eq}$ is the equivalent time delay. A second-order response was also observed in [35] for a similar class of servos.

3.1. Servo Motors

Several actuators were tested during the flight test campaign.

Table 2. Datasheet values for the tested servos.

	Torque (kgf$\mathit{·}$cm)	Speed ($\mathit{s}/60°$, °/s)	Voltage (V)
MG90S	1.8	0.1,    600	4.8
M5251H	3.3	0.04, 1500	7.4
M5252H	4.7	0.05, 1200	7.4
MS320	5.5	0.08,  750	7.4

gives an overview of the actuators and their primary datasheet values.

3.2. Excitation

A set of logarithmic sine-sweep maneuvers, as described in [4] (ch. 9), with length $T=12$ s, amplitude $A=10$ deg, starting frequency of $f_0=0.5$ Hz (3.1 rad/s) and end frequency of $f_1=18$ Hz (113 rad/s) was used. The frequencies were chosen to cover the expected range of operation and capability of the actuators. An example of a time domain plot is found in Figure 6.

3.3. Coherence

The actuators will have non-linear effects such as rate and acceleration saturations, dead zone and backlash effects. The magnitude squared coherence function, denoted ${\gamma }_{xy}^2$, can be used to indicate the degree of linearity as a function of frequency, and can be used to identify at which frequencies non-linearities have an effect. The control law should be designed such that there is good separation between the crossover frequencies and the frequency where the coherence starts dropping. All frequency plots for both actuators and system dynamics include a measure of coherence, using the MATLAB mscohere() function with default parameters.

3.4. Model Estimation

A continuous time transfer function model is estimated using the method presented in [5] [ch. 11], by solving an optimization problem reducing the error between the model and a set of measured frequency responses. The cost function, which will also be used for estimation of a plant model, is given by [5]:

$$J=\frac{1}{n_{TF}}\sum _{l=1}^{n_{TF}}\left[\frac{20}{n_{\omega }}\sum _{\omega ={\omega }_1}^{{\omega }_{n_{\omega }}}{W}_{\gamma }\left(W_g{\left(\left|{\widehat{G}}_l\left(i\omega \right)\right|-\left|G_l\left(i\omega \right)\right|\right)}^2+W_p{\left(\angle {\widehat{G}}_l\left(i\omega \right)-\angle G_l\left(i\omega \right)\right)}^2\right)\right],$$

(2)

where $n_{TF}$ is the number of transfer function matrix entries to estimate (in cases of MIMO systems) and ${\omega }_1$ to ${\omega }_{n_{\omega }}$ are the $n_{\omega }$ frequency points, chosen to be spaced with (almost) equal distance on a logarithmic scale, such that the model match will be distributed equally on a standard bode plot. ${\widehat{G}}_l$ is one of the estimated frequency response pairs l (e.g., aileron to roll rate) from flight test data, and $G_l$ is the evaluation of the model at that frequency. The weighting in the cost function is chosen based on recommendations from [5] as follows: $W_g=1$ and $W_p=0.01745$, which will set 1 dB magnitude error equal to $7.57$ deg phase error. The coherence weighting recommended in [5] is used and is given by

$$W_{\gamma }={\left(1.58\left(1-e^{-{\gamma }_{xy}^2}\right)\right)}^2,$$

(3)

which reduces the weighting of the data points as the coherence decreases; e.g., when ${\gamma }_{xy}^2=0.6$ the weight on the squared error is reduced by $50\%$. The frequency response estimation is performed with MATLAB tfestimate() to obtain ${\widehat{G}}_l$, using Welch’s estimation method and the $H_1$ estimator. The optimization problem for the cost function in Equation (2) is solved using the MATLAB fminsearch() function, which uses the Nelder–Mead simplex direct search method, a derivative-free, unconstrained multivariate optimization method. Equally logarithmically spaced frequency points were chosen on a range where the crossover frequencies of the final control law were expected to be placed with a good margin on either side. The expected range of the controller crossover frequencies is $\omega \in [1:45]$ rad/s, based on engineering judgement by testing hypothetical actuator performance values and performing the associated control design. At high frequency, several of the actuators deviated from second-order behavior, and this part was neglected in the estimation by reducing the highest estimation frequency. In Figures (Figure 3, Figure 4 and Figure 5), the frequency points for model estimation are indicated with black dots, the model response is indicated with a red line in the region of estimation, and a yellow line for the extrapolated part. Each figure contains a marking in the top right with relevant actuator and maneuver setup. “xiLi B” indicates aileron (xi) on the left side (L) using the inner (i) deflection measurement device. The linkage type is indicated with (B) for ball link and (SL) for standard plastic bearing. Also, the amplitude (A), start and end frequency (f0 and f1) and sine-sweep length (T) is indicated. The estimated model parameters are given in

Table 3. Comparison of estimated parameters for the actuators on ground and in flight (around 20 m/s ± 2 m/s). Delay is in seconds; eigenfrequency and bandwidth is in rad/s. ${\omega }_{BW}$ denotes the bandwidth, i.e., the frequency where the magnitude first drops 3 dB from the DC value, and ${\omega }_{\Delta 60}$ is the frequency where phase is first dropped by ${60}^{\circ }$.

On Ground	$T_{eq}$	${\omega }_{0,A}$	${\zeta }_A$	$K_A$	${\omega }_{BW}$	${\omega }_{\Delta 60}$	J
Aileron
MG90S	0.016	31.9	0.45	0.87	42.2	19.8	10.9
MS320 ball link	0.030	43.6	0.78	0.99	39.4	15.8	9.7
M5251H	0.032	95.9	1.00	0.97	71.1	21.2	1.8
M5252H ball link	0.039	92.9	1.00	0.93	92.6	20.8	1.9
Ruddervator
M5252H	0.028	92.4	0.80	0.93	79.9	22.4	7.8
M5252H ball link	0.030	108.2	0.94	0.93	76.3	23.0	10.5
In Flight	$T_{eq}$	${\omega }_{0,A}$	${\zeta }_A$	$K_A$	${\omega }_{BW}$	${\omega }_{\Delta 60}$	J
Aileron
MG90S	0.014	31.3	0.42	0.81	42.4	20.5	21.9
MS320 ball link	0.032	46.4	0.77	0.87	42.4	15.9	7.4
M5251H	0.028	95.4	0.98	0.89	63.5	22.0	3.3
M5252H ball link	0.028	87.9	0.73	0.85	85.5	23.3	3.8
Ruddervator
M5252H	0.028	88.1	0.75	0.93	82.7	23.0	7.4
M5252H ball link	0.028	85.3	0.73	0.82	82.5	22.8	5.3

.

3.5. Actuator Discussion

The estimated actuator parameters can be used to make a few statements on the use for flight control and the considerations necessary in the control design.

• A second-order system with time delay can be used to reasonably represent the dynamics for most of the actuators both on ground and in flight, at least up to approximately the bandwidth frequency. The faster actuators M5251H and M5252H have the best model fit. It may be considered to increase the model order for the slower MG90S and MS320 to obtain a better fit.
• The reduction in bandwidth from ground to in flight at the given conditions is at maximum ∼10%.
• There can be significant in-flight bias (DC offset) between actuator/surface command and surface position.
• There is noticeable reduction in gain ($K_A$) between ground and in-air tests.
• The M5252H actuator, which has both high torque and high speed indicated in the datasheet, has significantly higher bandwidth than the other actuators.
• The frequency response of the actuator in flight is reproducible with good accuracy.
• High coherence values extend to higher frequencies for M5252H and M5251H than for the MG90S.

As will be shown later, the fastest dynamic mode of the lateral/directional dynamics for this aircraft is around 32 rad/s (the roll mode). This is basically equal to the bandwidth of the MS320 and MG90S actuator, and leaves few options to control this dynamic mode. The bandwidth of the M5252H actuator is twice as fast, and is more suitable for this control task. It seems to be a general issue for control of subscale aircraft that the frequency separation between the actuators and the dynamic modes is small. For further flight test and control law testing, the M5251H and M5252H actuators will be used due to their higher performance.

4. Identification of Flight Dynamic Model

A linear state space model of the rigid, bare airframe dynamics is sought for control design. This can be obtained in many ways, but for the purpose of this paper, a single linear model is estimated. This model will be used for control design of the iDPI controller, in order to demonstrate the applicability of the concept.

4.1. Model Structure: Linear-Directional Lateral Dynamics

The longitudinal and lateral-directional model is assumed to be sufficiently decoupled. For a given trim point, the lateral-directional dynamics is assumed to be described via a linear state space model on the form (see [121]):

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\Delta {\dot{r}}_K\\ \Delta {\dot{\beta }}_K\\ \Delta {\dot{p}}_K\\ \Delta \dot{\varphi }\end{array}\right]& =\left[\begin{array}{cccc}{N}_r& N_{\beta }& N_p& 0\\ Y_r-cos{\alpha }_0& Y_{\beta }& Y_p+sin{\alpha }_0& \frac{g}{V_0}cos{\theta }_0\\ L_r& L_{\beta }& L_p& 0\\ tan{\theta }_0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}\Delta r_K\\ \Delta {\beta }_K\\ \Delta p_K\\ \Delta \varphi \end{array}\right]\hfill \\ \hfill & +\left[\begin{array}{cc}{N}_{\xi }& N_{\zeta }\\ Y_{\xi }& Y_{\zeta }\\ L_{\xi }& L_{\zeta }\\ 0& 0\end{array}\right]\left[\begin{array}{c}\Delta \xi \\ \Delta \zeta \end{array}\right]+\left[\begin{array}{ccc}-N_{\beta }/V_0& -N_p& -N_r\\ -Y_{\beta }/V_0& -Y_p& -Y_r\\ -L_{\beta }/V_0& -L_p& -L_r\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{v}_g\\ p_g\\ r_g\end{array}\right]\hfill \end{array}$$

(4)

with the output equation given by:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\Delta {\beta }_A\\ \Delta n_y\end{array}\right]& =\left[\begin{array}{cccc}0& 1& 0& 0\\ \left(V_0/g\right)Y_r& \left(V_0/g\right)Y_{\beta }& \left(V_0/g\right)Y_p& 0\end{array}\right]\left[\begin{array}{c}\Delta r_K\\ \Delta {\beta }_K\\ \Delta p_K\\ \Delta \varphi \end{array}\right]\hfill \\ \hfill & +\left[\begin{array}{cc}0& 0\\ \left(V_0/g\right)Y_{\xi }& \left(V_0/g\right)Y_{\zeta }\end{array}\right]\left[\begin{array}{c}\Delta \xi \\ \Delta \zeta \end{array}\right]\hfill \\ \hfill & +\left[\begin{array}{ccc}-1/V_0& 0& 0\\ -Y_{\beta }/g& -\left(V_0/g\right)Y_p& -\left(V_0/g\right)Y_r\end{array}\right]\left[\begin{array}{c}{v}_g\\ p_g\\ r_g\end{array}\right]\hfill \end{array}$$

(5)

where $V_0$ is the trim true airspeed, the states $r_K$, ${\beta }_K$, $p_K$ and $\varphi$ are the kinematic yaw rate, side slip angle, roll rate and bank angle, and the outputs ${\beta }_A$ and $n_y$ are the aerodynamic side slip angle and the lateral load-factor in the aerodynamic frame at the center of gravity. The lateral load-factor is equal to the specific force in the lateral direction denoted $f_y/g$. The inputs $\xi$ and $\zeta$ are the aileron and rudder deflections, and $v_g$, $p_g$ and $r_g$ are the side gust velocity and rotational gust components. $N_r$, $N_{\beta }$, $N_p$, $N_{\xi }$ and $N_{\zeta }$ are the yaw moment derivatives; $Y_r$, $Y_{\beta }$, $Y_p$, $Y_{\xi }$ and $Y_{\zeta }$ are the lateral specific force derivatives; $L_r$, $L_{\beta }$, $L_p$, $L_{\xi }$ and $L_{\zeta }$ are roll moment derivatives. The delta $(\Delta )$ denotes deviations from the trim state.

4.2. Experiment Design

In practice, all maneuvers have to be flown line of sight. To have proper control over the general flight path of the aircraft, allowing for sufficient space to maneuver at the end of each maneuver, a maximum maneuver (recording) length is given by $T_{rec}=14$ s, which includes a 1 s straight segment before the excitation starts, a 12 s logarithmic sine-sweep and then again a 1 s segment with no excitation (the sine-sweep could be improved by fading in/out the response to avoid the sharp beginning/end). The 1 s segments allow one to verify correct trim conditions before the beginning of the maneuver, and allow for the response to fade away after the maneuver. The frequency was varied from 0.5 Hz to 18 Hz, with an amplitude of ${10}^{\circ }$. This frequency range was found to cover the dynamic modes of the lateral dynamics and the effects of the actuators. The amplitude was found as a compromise between obtaining a sufficient signal-to-noise ratio and not forcing the aircraft too far away from the trim point. Figure 6 and Figure 7 show examples of the used excitation maneuver. The maneuver was repeated three times (in three different runs), and then concatenated afterwards to form a sequence with length $T_F=42$ s. The concatenation allowed a larger window size when estimating the frequency response functions with Welch’s estimate. This results in an increased frequency resolution with higher coherence. A relatively high frequency resolution is required to prevent that the peaks from the dutch roll damping become smeared out in the estimated frequency response, which would result in an overestimation of the damping ratio. The multisteps were performed with asymmetric excitation in order to avoid drift away of the aircraft after the maneuver.

4.3. Estimation Data

A single sine-sweep run of the lateral and the directional channel (flown separately) is shown on the left in Figure 6 and Figure 7, respectively. For the sine-sweep on the (virtual) rudder, the low dampened dutch-roll mode is clearly visible in the yaw rate and lateral load factor response. Also, the direct feed-through from the virtual rudder to lateral load-factor is clearly seen at the high frequency part of the sweep as a constant amplitude response in the lateral load-factor. The angular acceleration estimates, ${\widehat{\dot{p}}}_{pp}$ and ${\widehat{\dot{r}}}_{pp}$, were obtained via post processing the angular rates, with a zero-phase forward–backward derivative filter, and are shown for reference only; the data were not used for estimation. The aircraft was flown between 15 and 50 m above the ground for the system identification and control law verification. It was sought to perform the flight test in low turbulence and low wind conditions. For the data used in this paper, the wind speed was below 4 m/s.

4.4. Parameter Estimation

The model parameters to be estimated are listed in

Table 4. Estimated model parameters and an approximate 95% confidence interval based on the diagonal of the inverse of the Hessian matrix.

$N_r$	−1.825	±0.744	$Y_r$	0.016	±0.030	$L_r$	1.638	±6.003
$N_{\beta }$	5.366	±4.847	$Y_{\beta }$	−1.218	±0.415	$L_{\beta }$	−218.730	±52.031
$N_p$	−3.058	±0.653	$Y_p$	−0.081	±0.058	$L_p$	−31.440	±7.043
$N_{\xi }$	−18.301	±4.396	$Y_{\xi }$	−0.518	±0.532	$L_{\xi }$	−239.410	±46.642
$N_{\zeta }$	−20.024	±2.499	$Y_{\zeta }$	0.271	±0.040	$L_{\zeta }$	19.672	±16.865

, in addition to the input–output time delay $T_d$. The parameter estimation is performed using the cost function in Equation (2). The model containing the parameters to be estimated is given by Equations (4) and (5), connected with the sensor filters as described in Appendix C.1. This allows one to use the surface deflection measurements as input and the sensor measurements as outputs. For the remaining parameters, $V_0$ is taken as the average measured velocity for the test runs, for ${\theta }_0$ the average of the estimated initial pitch angle is used and for the trim angle of attack, the relation ${\alpha }_0={\theta }_0-{\gamma }_0$ is used, where ${\gamma }_0$ is taken as the initial inertial flight path angle estimated from the GPS. The frequency response estimation and optimization is solved in the same manner as the servo-dynamics estimation in Section 2, using the default options for the Welch estimate. As opposed to the SISO servo estimation, the estimation for the plant uses a MISO estimate (using MATLAB tfestimate(___,‘mimo’)) of the frequency response pairs; i.e., both virtual rudder and aileron were used as input, although only the estimation related to the input containing the sine-sweep was used. The other input is there to compensate for the fact that corrective aileron was used during some virtual rudder maneuvers and vice versa to keep the aircraft around the trim point. Five frequency response input–output pairs were chosen for the model estimation, $r/\xi$, $p/\xi$, $r/\zeta$, $n_y/\zeta$, $p/\zeta$. Aileron to lateral load factor were left out as there is almost no coupling, and the two bank angle pairs were left out, as the bank angle estimate arises from an internal estimation filter, rather than a direct measurement. Often this can be compensated by flight path reconstruction, but was left out in order to avoid introducing more hyperparameters. The frequency estimation region was chosen for the individual pairs according to the quality of the data, i.e., coherence and noise on the frequency response estimate.

4.5. Estimation Results

The estimated model parameters are given in Table 4. The input–output time delay was estimated to be 0 s, given by the lower boundary (note that this is the time delay between surface deflection and flight dynamic response, which does not contain the time delay from the actuator commands to the surface deflection response, which was estimated in Section 3). The minimum of the optimization was found at $J=30.9$, which is below the guideline of 50 from [5]. The fit of the model for the considered frequency response pairs is shown in Figure 8, Figure 9 and Figure 10. It is seen that the dutch roll is well identified in the virtual rudder response pairs. In the rudder to lateral load factor frequency response, it was chosen to not include the amplification increase at high frequency (above 50 rad/s), as indicated by the yellow line. As the controller crossover frequencies are well below this value, it was not deemed necessary to increase the model order to capture these dynamics. Also, the aileron to roll-rate frequency response has some discrepancy at high frequency. The flight test data for this estimation had different actuators for each aileron, namely M5252H and M5251H with different bandwidth, which may have led to asymmetric excitation at high frequency. These discrepancies motivate the additional use of non-parametric stability margin assessment considered in Section 5.6. The non-parametric assessment contains all the un-modelled dynamics, and hence provides a higher confidence in the stability and robustness of the control law for the envelope point of the test data. The correlation between the parameter estimates was calculated using a numerical approximation of the Hessian matrix, at the optimal solution for parameters. The estimated correlation matrix is given in Figure 11. It is seen that the correlation is low between most parameters. A general guideline is that a correlation above 0.9 for a given parameter pair can cause problems due to collinearity [4]. High correlations are present for the aileron effectiveness and roll rate parameters, in addition to the high correlation between the side-force coefficient due to side-slip and roll rate. As can be seen from the time domain multistep shown in Figure 6, a small oscillation at around 6 Hz is present at all times in the roll rate response. Further investigation will be conducted into the cause and effect of this, which could arise from wing structural vibration, turbulence or another source.

The identified model is also validated open loop in the time domain with a multistep response as shown in Figure 6 and Figure 7 to the right. The model data are an open-loop simulation with surface deflection measurements as input, and sensor model values as output. It is seen that the model represents reasonably well the dynamics of the aircraft also in the time domain.

The identified model is also validated closed-loop in the time domain with a multistep response on the lateral-load factor command as shown in Figure 25. The model data are a closed-loop simulation with the plant model, the actuator model and the developed controller, starting from the trim value. It is seen that the model reasonably well represents the closed-loop dynamics of the aircraft and controller also in closed loop, validating its usefulness for control law design.

5. Control Design and Evaluation

In the following section, the control design and control evaluation through flight test will be presented. The control technique will be based on the incremental-differential-proportional-integral (iDP) technique as introduced in [1]. As described in [1], the advantages of the iDPI control law structure over a classical PI control law are:

• Avoids hidden coupling terms when scheduling the controller gains;
• Bump less transfer between different gain groups;
• Enhanced integrator wind-up prevention;
• Robust performance with respect to variations in model parameters;
• Enhanced disturbance rejection;
• Less sensitive to unknown time delays at the input.

5.1. Control Law Structure

In order to control the lateral-directional motion of the aircraft, a MIMO control law is designed using eigenstructure assignment for the gain design. The control variables are chosen to be the bank angle $\varphi$ and the lateral load-factor in the body frame $n_{y,B}$. As shown later, a command filter is used to create a rate command attitude hold behavior for the roll channel such that the piloted command is roll rate $p_c$ or ${\dot{\varphi }}_c$. Steering commands from an outer-loop can be directly applied as bank angle commands ${\varphi }_c$. In order to apply the iDPI control law structure, the measurements will have to be differentiated, and the commands have to be applied to the incremental loop with ${\overline{\mathbf{u}}}_c$, as shown in Figure 12. The derivative filters before the controller gains cancel with the virtual integrator filter after the control gains in the fixed gain linear analysis, and remove the hidden coupling terms when analyzing the gain scheduled controller, as shown in [1]. Instead of a filtered derivative of the side-slip angle, filtering the kinematic relation for the aerodynamic side-slip angle rate is used as feedback, as derived in Section 5.2. This allows for a pseudo full state feedback of the flight dynamic variables of the lateral motion described by Equation (4). The control law block diagram is shown in Figure 12.

The bank angle command is given by ${\varphi }_c$ and the lateral load-factor command by $n_{y,c}$. ${\overline{\mathbf{u}}}_c$ is the unscaled input command, which is scaled to the actual actuator command $u_c$ via multiplication of the scale factor matrix ${\tilde{\mathbf{K}}}_A^{-1}$. The scale factor matrix is calculated such that the commands from ${\overline{\mathbf{u}}}_c$ to $\mathbf{u}$ have unit gain. The actual implementation of the incremental loop in this flight test campaign was carried out using the actuator models ${\widehat{G}}_A$ from Section 3, instead of the actual measurements, as shown in Figure 13. ${\widehat{\mathbf{G}}}_A$ is the transfer function matrix with the respective SISO actuator models on the diagonal.

In order to create a rate command attitude hold behavior, an integrator is used to transform a roll-rate $p_c$ (or ${\dot{\varphi }}_c$) command into the bank angle command for the control law, as given in Equation (6). Additionally, a zero is added to hide the spiral pole from the roll-rate command to bank angle response, in a similar fashion to [88]. The gain design moved the spiral pole to the left in order to increase the rejection of disturbances to the bank angle. The prefilter equation is given in Equation (6) and the block diagram in Figure 14, where ${\tau }_S$ is the closed-loop spiral pole.

$${\varphi }_c=\frac{\left(\frac{1}{{\tau }_S}s+1\right)}{s}{p}_c$$

(6)

5.2. Flow Angle Derivatives for State Feedback

For the lateral motion, classical controllers often employ feedback of either lateral load-factor or side-slip angle in order to modify the frequency of the dutch-roll mode. In both DPI and iDPI control law structures, the feedback signals must be differentiated. In order to avoid differentiating the lateral load-factor or side-slip angle (which may be be available through noisy flow angle measurement or through an uncertain model based estimate), the kinematic relation for the aerodynamic side-slip angle rate can be used, which does not contain any model parameters. The kinematic relations for the aerodynamic flow angles were derived in [70], by considering the translational equations of motion, inserting the relation between aerodynamic, kinematic and wind velocity, and then expanding using the rotation of the aerodynamic frame with respect to the body frame. The resulting equations are given by Equation (7):

$$\begin{array}{cc}\hfill {\left(\begin{array}{c}{\left({\dot{V}}_A^G\right)}^{EA}/{\left(V_A^G\right)}^E\\ {\dot{\beta }}_A\\ {\dot{\alpha }}_{A}cos{\beta }_A\end{array}\right)}_A=& \frac{g}{{\left(V_A^G\right)}^E}{\mathbf{M}}_{AB}\left(\frac{{\left(\mathbf{f}\right)}_B}{g}+{\mathbf{M}}_{BO}\frac{{\left(\mathbf{g}\right)}_O}{g}\right)\hfill \\ & -{\left(\begin{array}{c}0\\ r_{K}cos{\alpha }_A-p_{K}sin{\alpha }_A\\ p_{K}cos{\alpha }_{A}sin{\beta }_A-q_{K}cos{\beta }_A+r_{K}sin{\alpha }_{A}sin{\beta }_A\end{array}\right)}_A\hfill \\ & -\frac{1}{{\left(V_A^G\right)}^E}{\mathbf{M}}_{AB}{\mathbf{M}}_{BO}{\left({\dot{\mathbf{V}}}_W^G\right)}_O^{EO}\hfill \end{array}$$

(7)

where ${\left(V_A^G\right)}^E$ is the magnitude of the aerodynamic velocity (TAS), superscript G refers to the center of gravity of the aircraft, subscript A denotes that the velocity is the aerodynamic velocity and superscript E indicates that the velocities are taken with respect to the earth fixed earth centered reference frame. The scalar time derivative of the true airspeed in the aerodynamic frame is given by ${\left({\dot{V}}_A^G\right)}^{EA}$, where the second superscript A denotes the derivative of the velocity is taken in the A frame. ${\mathbf{M}}_{AB}$ is the rotation matrix between the body frame B and the aerodynamic reference frame A given in Equation (A18), and ${\mathbf{M}}_{BO}$ is the rotation matrix between the NED frame O and the body frame, given in Equation (A19). ${\left(\mathbf{f}\right)}_B$ is the (mass) specific forces denoted in the B frame, which is the quantity measured by the accelerometers. ${\left(\mathbf{g}\right)}_O$ is the gravity vector in the NED frame. ${\left({\dot{\mathbf{V}}}_W^G\right)}_O^{EO}$ is the change of the wind velocity in the NED frame with respect to time. In a constant wind field, this term is zero. If the wind gusts’ effects on the side-slip angle rate are not considered, the last term in Equation (7) can be neglected. The resulting side-slip rate estimation is given by:

$${\widehat{\dot{\beta }}}_A=\frac{g}{{\left(V_A^G\right)}^E}{\mathbf{e}}_y^T{\mathbf{M}}_{AB}\left(\frac{{\left(\mathbf{f}\right)}_B}{g}+{\mathbf{M}}_{BO}\frac{{\left(\mathbf{g}\right)}_O}{g}\right)-r_{K}cos{\alpha }_A+p_{K}sin{\alpha }_A$$

(8)

with ${\mathbf{e}}_y^T$ being the unit vector in the y direction of the associated frame (in this case in the A frame). For many aircraft, including small or subscale aircraft, a measurement or estimate of ${\alpha }_A$ and ${\beta }_A$, as, for example, investigated in [122], may not be feasible. An approximation may be obtained considering ${\mathbf{M}}_{AB}\approx \mathbf{I}$, $r_{K}cos{\alpha }_A\approx r_K$ and $p_{K}sin{\alpha }_A\approx 0$. This results in a side-slip angle rate estimate:

$${\widehat{\dot{\beta }}}_A=\frac{g}{{\left(V_A^G\right)}^E}\left({\left(n_y\right)}_B+cos\theta sin\varphi \right)-r_K$$

(9)

where ${\left(n_y\right)}_B=\frac{f_{y,B}}{g}$ and

$${\mathbf{M}}_{BO}\frac{{\left(\mathbf{g}\right)}_O}{g}={\left(\begin{array}{c}-sin\theta \\ cos\theta sin\varphi \\ cos\theta cos\varphi \end{array}\right)}_B.$$

(10)

A similar expression can be derived for the approximate angle of attack rate:

$$\widehat{\dot{\alpha }}=\frac{-g}{{\left(V_A^G\right)}^E}\left({\left(n_z\right)}_B-cos\theta cos\varphi \right)+q_K$$

(11)

It is important to note that the simplifications in Equations (9) and (11) lead to the correct steady condition (i.e., if ${\dot{\beta }}_A=0$, then also ${\widehat{\dot{\beta }}}_A=0$) in a non-accelerated straight flight. This observation is important for the implementation of iDPI, as this results in correct steady-state tracking of the control variable in non-accelerated straight flight. The expressions for Equations (9) and (11) are included in the gain design model, such that the effects of the simplifications are minimized.

5.3. Turn Coordination

Turn coordination is applied to the body angular rate measurements in order to account for the non-linear effects caused by the rolling and yawing motion during a roll and a subsequent turn. This will minimize side-slip excursions due to the non-linearities. The correction is applied in order to avoid the controller reacting to the non-linear contribution of the body angular rates caused by a roll/turn. The incremental angular rate feedback is given by:

$$\begin{array}{cc}\hfill \Delta p_{K,turn}=& p_K-p_{K,turn}\hfill \\ \hfill \Delta r_{K,turn}=& r_K-r_{K,turn}\hfill \end{array}$$

(12)

where the turn contributions are given by:

$$\begin{array}{cc}\hfill p_{K,turn}=& -\frac{g}{{\left(V_A^G\right)}^E}tan\varphi sin\theta \hfill \\ \hfill r_{K,turn}=& \frac{g}{{\left(V_A^G\right)}^E}cos\theta sin\varphi \hfill \end{array}$$

(13)

as derived in Appendix A.

5.4. Gain Design

The controller gain matrices ${\mathbf{k}}_{\varphi i}$, ${\mathbf{k}}_{n_{y}i}$, ${\mathbf{k}}_r$, ${\mathbf{k}}_{\beta }$, ${\mathbf{k}}_p$ and ${\mathbf{k}}_{\varphi }$ are calculated via eigenstructure assignment as in [88]. The perceived control law integrator behavior was removed in the command path by placing the zeros introduced by the feed forward gains ${\mathbf{h}}_{\varphi }$, ${\mathbf{h}}_{n_y}$ on the respective closed-loop integrator in the relevant command/control-variable input–output pair. The two sets of gain design parameters in eigenstructure assignment are (1) the location of the poles and (2) the corresponding eigenvector structure defined by a design matrix. Several requirements have to be fulfilled, and are verified using the design and assessment models. Not all requirements have a straightforward relationship with the gain design parameters, and an iterative process is applied, with variation in the gain design parameters in the eigenstructure assignment until all requirements are fulfilled.

The design model is constructed using MATLAB state space models for the plant, sensor, actuator, actuator reference model, digital signal processing and controller filters. The models are connected to a design model using the MATLAB connect() command to generate one design state space model with states:

$$\mathbf{x}=\left(\begin{array}{c}{\mathbf{x}}_{\mathit{fdm}}\\ {\mathbf{x}}_{\mathit{sens}/\mathit{dsp}}\\ {\mathbf{x}}_{\mathrm{act}}\\ {\mathbf{x}}_{\mathit{act},\mathit{mdl}}\\ {\mathbf{x}}_{\mathit{filters}}\end{array}\right),$$

(14)

with input ${\overline{\mathbf{u}}}_c$ and output being $(\varphi ,n_{y,B},\Delta r_{K,turn},{\widehat{\dot{\beta }}}_A,\Delta p_{K,turn},\varphi )$ filtered with the appropriate filters in Figure 12. Note that there are no dedicated integrator states, as these are inherently contained in the incremental loop.

The requirements to be fulfilled consist of:

• Eigenmode characteristics as given in

  Table 5. Open-loop and desired (subscript d) closed-loop eigenvalue location (omega in rad/s).

  	$\mathit{\omega }$	$\mathit{\zeta }$	${\mathit{\omega }}_{\mathit{d}}$	${\mathit{\zeta }}_{\mathit{d}}$
  Roll Mode	32.0	–	32.0	–
  Dutch-Roll Mode	5.6	0.19	4.0	0.85
  Spiral Mode	0.17	–	1.25	–
  Integrator $\varphi$	–	–	70.0	–
  Integrator $n_y$	–	–	75.0	–

  . The dutch-roll damping was increased and the spiral pole was further stabilized, as also shown in Figure 15 and Figure 16.
• Eigenmode response decoupling by assigning the desired $\varphi /\beta$ ratio as given in

  Table 6. Phi/beta ratio ($\left|\varphi \right|/\left|\beta \right|$).

  	$\left|\mathit{\varphi }\right|/\left|\mathit{\beta }\right|$
  Open-Loop	1.27
  Closed-Loop	0.1

  by modifying the desired eigenvector structure given by

  Table 7. Assignment of closed-loop eigenvector entries. $e_{\varphi }$ and $e_{n_y}$ are the control law (virtual) integrator states.

  . This can be compared to the open-loop eigenvectors given in

  Table 8. Open-loop eigenvector entry magnitudes.
• Integral behavior by assigning closed-loop integrator poles as given by Table 5.
• Gain margin (>6 dB) and phase margin (>${45}^{\circ }$). Additionally, it is sought to have the time delay margin above four samples (>0.04 s). The final design margins are shown in Figure 17.
• Limit the gain and phase crossover frequencies to below 45 rad/s in order to provide sufficient margin against model uncertainties at higher frequencies.
• Limitations on maximum actuator rate and position in relation to step commands and gusts as depicted in Figure 18 and Figure 19. The following limits are empirically found to work, but should be increased and verified via further flight test. For a step gust of $v_g=3$ m/s or command $n_{y,c}=0.2$, the actuator rate is sought to be less than $10\%$ of the unloaded maximum, and the actuator position less than $50\%$ of the maximum.

Several additional requirements could be added, including turbulence response, noise suppression, attenuation of structural modes, additional robustness metrics, bandwidth, handling qualities, etc.; however, these are not directly relevant for the current investigation, and could be added for further work.

5.5. Implementation

The controller was implemented in Simulink in order to create auto-generated code to be flashed onto the flight control computer. The continuous time control law is converted to a discrete time implementation. The filters depicted in Figure 12 with time constant ${\tau }_f$ were implemented with a forward Euler discrete integrator. As the time constant ${\tau }_f=0.125$ is very high compared to the sampling time $T_s=0.01$, the discretization effects on the frequency response are minor compared to the advantage of the simplicity of implementation, compared to, for example, a Tustin transform with prewarping. For the actuator models in ${\widehat{\mathbf{G}}}_A$ depicted in Figure 13, a zero-order hold transform is used. Referring to Figure 21, it is noted that the controller outputs actuator commands with a zero-order hold (ZOH) in between commands, and hence this transform is used to represent the dynamics of the continuous actuators at the sampled time instants.

5.6. Assessment

An automatized assessment is performed to ensure that the control law is safe to fly. The assessment is performed using three different models:

• Gain design model (GDM):
  Linear continuous time model, possibly with simplifications used to calculate the controller gains.
• Linear parametric assessment model (LPAM):
  Discrete time linearization of the controller implementation model connected with the zoh transform of the identified linear continuous time parametric models of the actuators, plant and sensors. The zero order hold transform is used because the controller outputs commands that are held between each time update of the controller, i.e., zero order hold of the commands; using the zoh transform for assessment correctly represents the dynamics of the sampled closed-loop continuous time system.
• Non-parametric assessment model (NPAM):
  Discrete time linearization of the controller implementation model connected with the non-parametric frequency response estimates from flight test.

The interconnection between the three models is sketched in Figure 20.

The gain design model (GDM) is described in Section 5.4.

The LPAM is created using the continuous time linear models for the actuators ${\mathbf{G}}_{\mathit{A}}\left(s\right)$, plant ${\mathbf{G}}_{\mathit{P}}\left(s\right)$ and sensors ${\mathbf{G}}_{\mathit{S}}\left(s\right)$, and discrete time linear models for the digital signal processing ${\mathbf{C}}_{\mathit{dsp}}\left(z\right)$ and the discrete time controller $\mathbf{C}\left(z\right)$. Even though the digital signal processing model contains the discrete time digital downsampling filters as described in Appendix C.1, these are considered continuous and embedded in the sensor model ${\mathbf{G}}_{\mathit{S}}\left(s\right)$ in the analysis, which may be reasonable as they are sampled 10× faster than the controller, as shown in Figure 21 with $T_{s,1}=0.001$ and $T_s=0.01$. An alternative would be to represent both the controller and dsp at the higher sampling frequency, i.e., upsampling the controller with no loss of accuracy. No additional dynamics of the sensor was found in the addition to the digital signal processing (it is common to have a filter in the sensor for the downsampling; however, the PX4 software moved the filtering to the flight control computer, which samples at high rate). The dynamics of the estimate of $\varphi$ was neglected. The continuous time actuator model, plant model and in this analysis also the combined sensor/dsp model are discretized with the zero order hold transform, using the same argument as in Section 5.5; i.e., the controller outputs zero order hold commands. This will produce correct evaluations when considering stability margins, but may lead to slight errors when considering frequency responses for gusts and turbulence, the inputs of which are in continuous time. However, since the sampling frequency is much higher than the frequency of the turbulence, this is also deemed acceptable for this analysis. The LPAM model creation is fully automatic using Simulink slLinearizer, which automatically linearizes and discretizes the controller implementation model together with the actuator model, plant model and sensor model.

The NPAM is created using the linearized controller implementation model, similar to the LPAM. However, instead of using the parametric models for the other subsystems, the frequency response functions (FRF) obtained directly from flight test are used. These FRFs contain all the effects on the real aircraft, and give the most accurate pointwise estimate of the stability and robustness characteristics, as depicted in Figure 20. The SISO loop breaks are calculated using Equations (17) and (18), which ensures only the loop break under consideration is open; i.e., the other loop break is closed. The structure is depicted in Figure 22. For each frequency ${\omega }_k$ in the frequency response estimate, assemble the frequency response estimate matrix ${\mathbf{G}}_{FRF}\left(i{\omega }_k\right)$:

$${\widehat{\mathbf{G}}}_{FRF}=\left[\begin{array}{cc}{\widehat{G}}_{p_K,{\xi }_c}& {\widehat{G}}_{p_K,{\zeta }_c}\\ {\widehat{G}}_{r_K,{\xi }_c}& {\widehat{G}}_{r_K,{\zeta }_c}\\ {\widehat{G}}_{n_{y,B},{\xi }_c}& {\widehat{G}}_{n_{y,B},{\zeta }_c}\\ {\widehat{G}}_{\varphi ,{\xi }_c}& {\widehat{G}}_{\varphi ,{\zeta }_c}\end{array}\right]$$

(15)

Assemble the helper matrix $\tilde{\mathbf{G}}\left(i{\omega }_k\right)$:

$$\tilde{\mathbf{G}}\left(i{\omega }_k\right)=\mathbf{C}\left(e^{i{\omega }_k{T}_s}\right){\widehat{\mathbf{G}}}_{FRF}\left(i{\omega }_k\right)$$

(16)

where the discrete controller transfer function is evaluated in $z=e^{i{\omega }_k{T}_s}$ using MATLAB’s freqresp(). The loop transfer function with the rudder ${\zeta }_c$ closed can be calculated by inspecting Figure 22 as shown in Appendix D.1:

$$L_{\xi }\left(i{\omega }_k\right)={\tilde{\mathbf{G}}}_{11}\left(i{\omega }_k\right)+{\tilde{\mathbf{G}}}_{12}\left(i{\omega }_k\right){\left(1-{\tilde{\mathbf{G}}}_{22}\left(i{\omega }_k\right)\right)}^{-1}{\tilde{\mathbf{G}}}_{21}\left(i{\omega }_k\right)$$

(17)

and similarly for the rudder ${\zeta }_c$:

$$L_{\zeta }\left(i{\omega }_k\right)={\tilde{\mathbf{G}}}_{22}\left(i{\omega }_k\right)+{\tilde{\mathbf{G}}}_{21}\left(i{\omega }_k\right){\left(1-{\tilde{\mathbf{G}}}_{11}\left(i{\omega }_k\right)\right)}^{-1}{\tilde{\mathbf{G}}}_{12}\left(i{\omega }_k\right)$$

(18)

The loop transfer function shapes are shown in Figure 23 and Figure 24 with results from GDM, LPAM and NPAM.

5.7. Lateral Flight Test Results

The control law was flown in flight test to verify the behavior and applicability to a subscale remote controlled aircraft. The closed-loop response was compared against the design model, and showed very accurate prediction of the closed-loop time domain response. A multistep test for the lateral load-factor is shown in Figure 25, where the expected design model response is plotted on top of the actual aircraft response. Note that the gains in the flight test slightly differed from those above because of changes in aircraft configuration between flights and an intermittent crash; however, the model comparison is made with a matching gain set. This is also the reason why there is no bank angle step response evaluation. The flight test gain set had a desired $\varphi /\beta$ ratio of 0.9, which causes the higher bank angle coupling. There is a slight difference between the bank angle response and the design model; however, the average trend is correct. This discrepancy may arise from the oscillations on the roll rate measurements as discussed in Section 4.5, which may translate into the bank angle estimate. It is also notable that both the virtual rudder and aileron commands have small high-frequency oscillations, roughly at the same frequency as the oscillations on the yaw rate and roll rate. These oscillations also exists in open loop; see Figure 6; hence, they seem to not originate due to the control law. Further investigation should be conducted to see whether notch filtering can be used to filter out the oscillations in the measurement, and avoid them propagating to the actuator commands, as done in [1].

6. Conclusions

This paper investigated the feasibility of practical system identification and incremental control design for a subscale fixed wing aircraft. A subscale demonstrator platform was equipped with sensors for system identification and control, in particular surface deflection measurements, which is uncommon for this class of aircraft. An iDPI control law was synthesized to control the lateral-directional motion, and it was shown that the control law performed well in flight test. It was shown that the identified model was able to accurately predict the closed-loop behavior, motivating the use of the presented process for further control law development. It was shown that eigenstructure assignment could successfully be applied together with iDPI to form a multiple input/multiple output incremental control law. The controller was shown to be easily implemented in discrete time and deployed on a flight control computer. The control law was flight tested on the subscale aircraft, and it was demonstrated that the nominal performance was as designed. In addition, several actuator models were identified using surface deflection measurements, both on ground and in flight. These models were used to support proper control design. Several improvements can be made for future development, including composite windowing techniques, increased sampling rates, estimation of non-dimensional aerodynamic/propulsion coefficients, full envelope gain scheduling, investigation of handling qualities, experimental estimation of actuator non-linearities, estimation of servo dynamic variation with airspeed, comparison of internal and external servo position measurements, optimization of input design, identification of structural couplings and notch filter design and incorporation.