Flight Control Design for a Tailless Aircraft Using Eigenstructure Assignment

Abstract

We apply eigenstructure assignment to the design of a flight control system for a wind tunnel model of a tailless aircraft. The aircraft, known as the innovative control effectors (ICEs) aircraft, has unconventional control surfaces plus pitch and yaw thrust vectoring. We linearize the aircraft in straight and level flight at an altitude of 15,000 feet and Mach number 0.4. Then, we separately design flight control systems for the longitudinal and lateral dynamics. We use a control allocation scheme with weights so that the lateral pseudoinputs are yaw and roll moment, and the longitudinal pseudoinput is pitching moment. In contrast to previous eigenstructure assignment designs for the ICE aircraft, we consider the phugoid mode, thrust vectoring, and stability margins. We show how to simultaneously stabilize the phugoid mode, satisfy MIL-F-8785C mode specifications, and satisfy MIL-F-9490D phase and gain margin specifications. We also use a cstar command system that is preferable to earlier pitch-rate command systems. Finally, we present simulation results of the combined longitudinal/lateral flight control system using a full 6DOF nonlinear simulation with approximately 20,000 values for the aerodynamic coefficients. Our simulation includes limiters on actuator deflections, deflection rates, and control system integrators.

1. Introduction

We consider the design of a flight control system using eigenstructure assignment for a wind tunnel model of the innovative control effectors (ICEs) aircraft. This tailless aircraft program was first described by Dorsett and Mehl [1] and by Dorsett et al. [2]. The ICE aircraft has many unconventional control surfaces plus pitch and yaw thrust vectoring. Several authors have proposed flight control system designs for the ICE aircraft. Ngo et al. [3] use dynamic inversion with structured singular value synthesis. However, the authors remove the bank angle equation from the model which causes an unstable complex mode to be replaced with an unstable real mode. This occurs because the ICE aircraft does not exhibit the conventional real spiral mode. Sparks [4] uses linear parameter-varying control. Schumacher and Johnson [5] use dynamic inversion with adaptation for self reconfiguring. Shtessel et al. [6] propose reconfigurable sliding mode control with direct adaptation. Hess et al. [7] use sliding mode control with asymptotic observers. However, the results are based only upon linear simulation. The only previous design using eigenstructure assignment was proposed by Jones et al. [8]. The main emphasis of [8] is a method for gain scheduling. Thus, no details of the eigenstructure assignment design are shown. Furthermore, [8] ignores the phugoid mode, only uses pitch flap for the longitudinal controller, uses a pitch-rate command system, and only uses elevon for the lateral control system.

Eigenstructure assignment is an excellent method for incorporating classical specifications on damping, settling time, and mode decoupling into a modern multivariable control framework. Andry et al. [9] applied eigenstructure assignment to the design of a constant gain output feedback aircraft flight control system. Andry et al. [9] proposed a choice for the desired eigenvectors for the lateral dynamics of an aircraft based upon mode decoupling. The desired eigenvectors in [9] are based upon decoupling yaw rate and sideslip angle from roll rate and bank angle. Later, Davidson and Andrisani [10] proposed another choice for desired eigenvectors. The desired eigenvectors in [10] are based upon the observations that bank angle is the integral of roll rate and bank angle should induce a yaw rate. The desired eigenvectors in [10] have only two specified entries which is fewer than the number of specified entries in [9]. This is an important point because the lateral dynamics has only two independent control directions which are yawing moment and rolling moment.

In this paper, we apply eigenstructure assignment to the design of a flight control system for the ICE aircraft at one flight condition. Significant changes in the flight condition of the aircraft would require repeating our design method at each of the flight conditions. Then, the controller would be gain scheduled as a function of altitude, Mach number, and possibly angle of attack and sideslip. We linearize the aircraft in straight and level flight at an altitude of 15,000 feet and Mach number 0.4. Then, we separately design flight control systems for the longitudinal and lateral dynamics. Our longitudinal controller is a cstar [11] command system using symmetric pitch flaps and pitch thrust vectoring. The cstar quantity is a blend of normal acceleration at the pilot station and pitch rate. Our lateral controller is a stability axis roll-rate command system using right and left elevons, right and left all moving tips, and yaw thrust vectoring. We use a control allocation scheme with weights, so that the lateral pseudoinputs are yaw and roll moment, and the longitudinal pseudoinput is pitching moment. In contrast to the design of [8], we consider the phugoid mode, thrust vectoring, and stability margins. We show how to modify the desired eigenvalues to simultaneously achieve acceptable responses to pilot commands, stabilize the phugoid mode, satisfy MIL-F-8785C [12] mode specifications, and satisfy MIL-F-9490D [13] specifications on phase and gain margins.

We show time responses for our combined longitudinal/lateral flight control system using a 6DOF nonlinear simulation of the ICE aircraft. We start with the aircraft control toolbox [14] for use with MATLAB. This toolbox implements the nonlinear 6DOF simulation of an F-16 aircraft that is described by Stevens and Lewis [15]. We modify the aircraft control toolbox [14] to include look-up tables with interpolation for approximately 20,000 values of the ICE aircraft aerodynamic coefficients [16]. Then, we modify the force and moment equations to include yaw and pitch thrust vectoring. We also modify the trim program cost function to penalize nonzero trim bank angle and non-1g trim normal acceleration. Finally, we add limiters on actuator deflections, deflection rates, and control system integrators. Our nonlinear simulation results verify that we have achieved desirable performance to pilot cstar and stability axis roll-rate commands.

2. Eigenstructure Assignment

Consider an aircraft modeled by the linear time invariant matrix differential equation described by 𝑥̇=𝐴𝑥+𝐵𝑢,𝑦=𝐶𝑥, (1) where 𝑥 is the state vector (𝑛×1), 𝑢 is the control vector (𝑚×1) and 𝑦 is the output vector (𝑟×1). It is assumed that the 𝑚 inputs and the 𝑟 outputs are independent. Also, as is usually the case with aircraft, it is asssumed that 𝑚<𝑟<𝑛. For our design, 𝑚 will be the number of pseudocontrols. If there are no pilot commands, the feedback control vector 𝑢 equals a matrix times the output vector 𝑦, 𝑢=−𝐹𝑦. (2)

The feedback problem can be stated as follows: given a set of desired eigenvalues (𝜆^𝑑_𝑖), 𝑖=1,2,…,𝑟 and a corresponding set of desired eigenvectors, (𝑣^𝑑_𝑖), 𝑖=1,2,…,𝑟, find the real 𝑚×𝑟 matrix 𝐹 such that the eigenvalues of A-BFC contain (𝜆^𝑑_𝑖) as a subset, and the corresponding eigenvectors of A-BFC are close to the respective members of the set (𝑣^𝑑_𝑖).

Srinathkumar [17] has shown that if (𝐴,𝐵) is a controllable pair, then the feedback gain 𝐹 will exactly assign 𝑟 eigenvalues. It will also assign the corresponding eigenvectors, provided that (𝑣^𝑑_𝑖) is chosen to be in the subspace spanned by the columns of (𝜆_𝑖𝐼−𝐴)⁻¹𝐵 for 𝑖=1,2,…,𝑟. This subspace is of dimension 𝑚, which is the number of independent control variables. A numerically superior representation for this subspace is described by Kautsky et al. [18], who showed that the subspace in which the eigenvector (𝑣^𝑑_𝑖) must reside is also given by the nullspace of 𝑈^𝑇₁(𝜆_𝑖𝐼−𝐴). Let 𝐿_𝑖 be a matrix whose columns are a basis for this nullspace. The matrix 𝑈₁ is obtained from the singular value decomposition of 𝐵, which is given by 𝐵=𝑈₀𝑈₁Σ𝑉^𝑇0. (3)

If we choose an eigenvector 𝑣_𝑖 which lies precisely in the subspace spanned by the columns of 𝐿_𝑖, it will be achieved exactly. In general, however, a desired eigenvector 𝑣^𝑑_𝑖 will not reside in the prescribed subspace and hence cannot be achieved. Andry et al. [9] describe a way to find the “best possible choice” for an achievable eigenvector. This best possible eigenvector is the projection of 𝑣^𝑑_𝑖 onto the subspace spanned by the columns of 𝐿_𝑖 (in the least square sense). In many practical situations, complete specification of 𝑣_𝑖 is neither required nor known, but rather the designer is interested only in certain elements of the eigenvector. Thus, assume that 𝑣^𝑑_𝑖 has the following structure: 𝑣^𝑑_𝑖=𝑣_𝑖1,𝑥,𝑥,𝑥,𝑥,𝑣_𝑖𝑗,𝑥,𝑥,𝑣_𝑖𝑛^𝑇, (4) where 𝑣_𝑖𝑗 are designer-specified components, and 𝑥 is an unspecified component. Define, as shown by Andry et al. [9], a reordering operation {}^𝑅𝑖 such that 𝑣^𝑑_𝑖^𝑅𝑖=ℓ_𝑖𝑑_𝑖, (5) where ℓ_𝑖 is a vector of specified components of 𝑣^𝑑_𝑖, and 𝑑_𝑖 is a vector of unspecified components of 𝑣^𝑑_𝑖. The rows of the matrix 𝐿_𝑖 are also reordered to conform with the reordered components of 𝑣^𝑑_𝑖. Thus, 𝐿_𝑖^𝑅𝑖=𝐿_𝑖𝐷_𝑖. (6)

Then, the achievable eigenvector 𝑣^𝑎_𝑖 is given by 𝑣^𝑎_𝑖=𝐿_𝑖𝑧_𝑖, (7) where 𝑧_𝑖=𝐿^†_𝑖ℓ_𝑖 and where (⋅)^† denotes the appropriate pseudoinverse of (⋅).

The feedback gain matrix 𝐹 is given by 𝐹=−𝑉_𝑏Σ_𝑏⁻¹𝑈^𝑇_𝑏0(𝑉Λ−𝐴𝑉)𝑉_𝑟Σ_𝑟⁻¹𝑈^𝑇_𝑟0, (8) where the singular value decomposition of the matrices 𝐵 and 𝐶𝑉 is given by 𝐵=𝑈_𝑏0𝑈_𝑏1Σ_𝑏𝑉^𝑇_𝑏0,𝐶𝑉=𝑈_𝑟0𝑈_𝑟1Σ_𝑟𝑉^𝑇_𝑟0, (9) and where Λ is an 𝑟×𝑟 diagonal matrix with entries 𝜆_𝑖, 𝑖=1,2,…,𝑟.

3. Flight Control System Design

3.1. ICE Linearized Model

We linearize the dynamics of the ICE aircraft at a speed of Mach 0.4 and an altitude of 15,000 feet. The state variables are velocity 𝑉_𝑇 (ft/s), angle of attack 𝛼 (rad), pitch angle 𝜃 (rad), pitch rate 𝑞 (rad/s), engine power level, sideslip angle 𝛽 (rad), bank angle 𝜙 (rad), roll rate 𝑝 (rad/s), and yaw rate 𝑟 (rad/s). The control effectors are throttle 𝛿_th (0-1), symmetric pitch flap 𝛿_pflap (deg), left elevon 𝛿_el (deg), right elevon 𝛿_er (deg), left all moving tip 𝛿_amtl (deg), right all moving tip 𝛿_amtr (deg), pitch thrust vectoring 𝛿_ptv (deg), and yaw thrust vectoring 𝛿_ytv (deg). The deflection limits are |𝛿_pflap|≤30^∘, |𝛿_elevon|≤30^∘, −30^∘≤𝛿_amt≤60^∘, |𝛿_ptv|≤15^∘, and |𝛿_ytv|≤15^∘. The deflection rate limits are |𝛿̇_pflap|≤50 deg/s, |𝛿̇_elevon|≤150 deg/s, |𝛿̇_amt|≤150 deg/s, |𝛿̇_ptv|≤60 deg/s, and |𝛿̇_ytv|≤60 deg/s.

The dynamics are decoupled into longitudinal and lateral design models. The longitudinal states are 𝑉_𝑇 (ft/s), 𝛼 (rad), 𝜃 (rad), 𝑞 (rad/s), and integral of cstar 𝑐^∗. Here, cstar [11, 15] is defined as 𝑐^∗=𝑛_𝑧𝑝+12.4𝑞, where 𝑛_𝑧𝑝 is the normal acceleration at the pilot's station. The control effectors are symmetric pitch flap 𝛿_pflap (deg) and pitch thrust vectoring 𝛿_ptv (deg). The sensor measurements are 𝛼 (deg), 𝑞 (deg), and 𝑛_𝑧𝑝(𝑔). The lateral states are 𝛽 (rad), 𝜙 (rad), 𝑝 (rad/s), 𝑟 (rad/s), and a yaw rate washout filter state 𝑥_wo with a one second-time constant. The control effectors are left elevon 𝛿_el (deg), right elevon 𝛿_er (deg), left all moving tip 𝛿_amtl (deg), right all moving tip 𝛿_amtr (deg), and yaw thrust vectoring 𝛿_ytv (deg). The sensor measurements are 𝛽 (deg), roll rate 𝑝 (deg/s), and yaw rate 𝑟 (deg/s). In addition, we use longitudinal and lateral simulation models which are the design models augmented with first-order actuator dynamics with time constant 0.049505 seconds.

3.2. Lateral Flight Control System

The linearized lateral states are 𝛽, 𝜙, 𝑝, and 𝑟. The open-loop lateral eigenvalues and eigenvectors are shown in Table 1. We observe that the aircraft has two complex lateral modes with one of the modes being unstable. This is in sharp contrast to a conventional aircraft that exhibits a complex dutch roll mode, a real roll subsidence mode, and a real spiral mode. Furthermore, the eigenvectors indicate that both modes are strongly coupled and that neither mode can be identified as a dutch roll.

Table 1 

Eigenvalues and eigenvectors for the lateral dynamics^a.

We transform the lateral dynamics from body axis to stability axis, so that we can decouple stability axis roll rate 𝑝_s from sideslip angle. This is consistent with our choice of 𝑝_s as the pilot command. This transformation requires the value of trim alpha which is 0.1569 rad. The states are now sideslip angle 𝛽, bank angle 𝜙, stability axis roll rate 𝑝_s, stability axis yaw rate 𝑟_s, and a washout filter state 𝑥_wo. We choose the lateral feedbacks to be 𝛽, 𝑝_s, and washed-out stability axis yaw rate (𝑟_s)_wo which allows us to assign three eigenvalues. The military specifications [12, 19] for level 1, category A, class IV flight require (1) a maximum roll mode time constant of one second, (2) dutch roll minimum 𝜁_dr of 0.4, minimum 𝜁_dr𝜔_𝑛 of 0.4 r/s, minimum 𝜔_𝑛=1.0 r/s, and (3) stable spiral mode or unstable spiral mode with minimum time to double amplitude 𝑇₂>12 s. The desired dutch roll and roll subsidence eigenvalues and eigenvectors are those suggested by Davidson and Andrisani [10] and are shown in Table 1. We note that the desired eigenvalues satisfy the military specifications. Furthermore, we note that the dutch roll mode is a (𝛽,𝑟_s) mode that is decoupled from 𝜙, whereas the roll subsidence mode is a 𝑝_s mode that is decoupled from 𝛽.

We show the singular values of the control distribution matrix 𝐵_stab in Table 2. Observe that there are only two nonzero singular values which indicate that there are only two independent control directions. We use the first two columns of 𝑈_stab as the pseudocontrol distribution matrix. Column 1 corresponds to rolling moment 𝛿_roll, and column 2 corresponds to yawing moment 𝛿_yaw. We can generalize this method by adding weights, to the five lateral control effectors. This will allow larger control deflections for the effectors with larger deflection limits and vice versa. Mathematically, let 𝑊 contain the weights, and consider the singular value decomposition 𝐵_stab𝑊=𝑈_stabΣ_stab𝑉^𝑇_stab, (10) then 𝑥̇=𝐴𝑥+𝐵𝛿, (11) where 𝐵=𝑈_stab(),𝛿∶,1∶2𝛿=−𝐾𝑦=_roll,𝛿_yaw^𝑇. (12) Then, the inverse mapping to obtain the control 𝑢(𝑡) is 𝑢(𝑡)=𝑉_stabΣ⁻¹_stab𝐾𝑊𝑦(𝑡). (13)

Table 2 

SVD of lateral control distribution matrix^a.

The feedback gain matrix, achievable eigenvalues, and achievable eigenvectors are shown in Table 3. We observe that we have achieved the desired eigenvalues and eigenvectors. We also observe that the thrust vector gains are larger than the elevon and amt gains. A linear simulation, that includes actuator dynamics, shows that the thrust vector deflection is larger than the elevon and amt deflections. This is undesirable because the thrust vector deflection limit is only 15 degrees, whereas the elevon and amt deflection limits are 30 degrees. To remedy this problem, we reduce the thrust vector weight in (10) while maintaining the other weights at unity. Table 3 shows the results for a thrust vector weight of 0.25 where we observe that the thrust vector gains are now smaller than the elevon and amt gains. A linear simulation shows that the maximum thrust vector deflection to a stability axis roll rate pulse is approximately 40% and 66% of the elevon and amt maximum deflections, respectively. We also observe from Table 3 that the spiral mode is unstable with eigenvalue 𝜆_spiral=0.0077. The time to double amplitude for the unstable spiral mode is given by [19] 𝑇₂=0.693=unstableroot0.6930.0077=90sec, (14) which satisfies the specification of 𝑇₂>12 seconds.

Table 3 

Lateral eigenstructure assignment design^a.

We analyze the stability robustness of the design using both open-loop phase and gain margins and closed-loop mu analysis. We use the linearized simulation model that includes actuator dynamics. First, we use function loop margin from the MATLAB robust control toolbox [20] to compute the disk margins. The disk margins represent the largest region for each channel, such that the closed loop system is stable for all gain and phase variations inside the region. The disk margin analysis allows for simultaneous gain and phase variations, one loop at a time. The unstable open-loop complex mode poses a problem for gain and phase margin calculation. For a conventional aircraft with a real unstable spiral mode, we would use a similarity transformation to decouple the spiral mode. Then, we would obtain gain and phase margins for the reduced-order model that does not include the spiral mode. Unfortunately, the best that can be done here is to remove the bank angle equation from the open-loop aircraft. This results in closed-loop eigenvalues that are approximately the same for the model with the bank angle equation except for the absence of the spiral mode. We believe that this is acceptable especially because we use disk margins for simultaneous gain and phase variation, whereas MIL-F-9490D [13] specifies either gain or phase variation. Thus, our analysis is conservative with respect to the specifications. The disk margins are shown in Table 4 where we observe that the specifications are achieved for all loops except the sideslip loop. We found that increasing the dutch roll damping improved the sideslip loop margins. We choose the dutch roll damping to be 𝜁_dr=0.7 which yields a sideslip loop phase margin of ±45.31 degrees as compared with the specification of ±45 degrees. The feedback gains are shown in Table 3 where we observe that the gain magnitudes are larger for 𝜁_dr=0.7 as compared to the gain magnitudes for 𝜁_dr=0.4.

Table 4 

Lateral disk gain and phase margins^a.

A block diagram of the lateral flight control system is shown in Figure 1. We have added a gain of 4/3 to the pilot stick for the purpose of achieving zero steady-state error to a 𝑝_s command. The linearized lateral time responses to a one-degree initial sideslip are shown in Figure 2. We observe that the bank angle is approximately 0.1 degrees. This yields |𝜙/𝛽|≈0.1 which indicates excellent decoupling from sideslip to bank angle.

Figure 1 

Block diagram of lateral flight control system.

(a)

(b)

(c)

(d)

Figure 2 

Linear time responses to a 1 deg initial sideslip.

We use function robuststab from the MATLAB robust control toolbox [20] to compute the margin of stability for real uncertainty in the stability and control derivatives. This analysis is based upon the closed-loop aircraft which exhibits an unstable real spiral mode. We use a similarity transformation to decouple this unstable spiral mode. This yields a lower-order closed-loop model without the spiral mode that is used for the robustness analysis. A stability robustness margin greater than one means that the uncertain system is stable for all values of its modeled uncertainty. The function robuststab computes upper and lower bounds on the exact stability margin using mu analysis. The uncertainty in the stability derivatives is given by [3] 𝑌_𝛽=15%, 𝐿_𝛽=10%, 𝐿_𝑝=30%, 𝐿_𝑟=20%, 𝑁_𝛽=30%, 𝑁_𝑝=50%, and 𝑁_𝑟=15%. The uncertainty is 15% in all control derivatives. To obtain a reliable estimate of the robustness margin, we add 5% uncertain dynamics to the real parameter uncertainty [20] and reduce the real uncertainty by 5%. The upper and lower mu bounds are shown in Figure 3 where we observe that the closed loop is robustly stable for up to 108% of the modeled uncertainty. We analyze performance robustness to the real uncertainty in the stability and control derivatives by performing a monte carlo simulation. We use the linearized lateral simulation model that includes first-order actuator dynamics. We perform 100 runs of the linearized lateral dynamics over the real uncertainty in the stability and control derivatives.This is accomplished by using functions ureal and usample from the MATLAB robust control toolbox[20].The uncertainty in each parameter is uniformly distributed and independent of the other parameters. The responses to a one deg/s stability axis roll-rate pulse command are shown in Figure 4. We observe that the responses are well behaved for all uncertainty and that deviations from nominal indicate a robust design.

Figure 3 

Stability margin bounds for lateral dynamics.

(a)

(b)

(c)

(d)

Figure 4 

Monte carlo simulation for lateral dynamics.

3.3. Longitudinal Flight Control System

The linearized longitudinal states are 𝑉_𝑇, 𝛼, 𝜃, 𝑞, and integral of cstar ∫𝑐^∗. The open-loop short period and phugoid modes exhibit natural frequencies and damping ratios given by (𝜔_nsp,𝜁_sp)=(1.10,0.57) and (𝜔_nph,𝜁_ph)=(0.107,0.11), respectively. The damping ratio specifications for the short period and phugoid modes for level 1, category A flight are [12] 0.35≤𝜁_sp≤1.30 and 𝜁_ph≥0.04, respectively. For small angles, the steady-state load factor per steady-state angle of attack is given by [21] 𝑛𝛼=𝑉_trim𝐿_𝛼𝑔≈−𝑈₀𝑍_𝑤𝑔. (15) For our linearized longitudinal dynamics, the trim velocity 𝑉_trim=422.96 ft/s, 𝑍_𝑤=−0.5972 which yields 𝑛/𝛼=7.85. Using MIL-F-8785C [12] with 𝑛/𝛼=7.85, we find that the specification for the short-period natural frequency is 1.6≤𝜔_nsp≤5.5. We note that the open-loop aircraft does not meet the specification on 𝜔_nsp.

We show the singular values of the control distribution matrix 𝐵_long in Table 5. Observe that although there are two nonzero singular values, the second singular value is ten times smaller than the first singular value. This indicates that there is only one strongly independent control direction which corresponds to pitch moment denoted by 𝛿_pitch. We use the first column of 𝑈_long as the pseudocontrol distribution matrix. Therefore, this is a single-input problem which allows eigenvalue but not eigenvector assignment. We add the weights (1, 0.65), so that 𝛿_pflap will have maximum deflection of approximately twice the 𝛿_ptv deflection. This choice is because the pitch flap has a deflection limit that is twice the thrust vector deflection limit. Mathematically, let 𝑊 contain the weights, and consider the singular value decomposition 𝐵_long𝑊=𝑈_longΣ_long𝑉^𝑇_long, (16) then, 𝑥̇=𝐴𝑥+𝐵𝛿, (17) where 𝐵=𝑈_long(),𝛿∶,1_pitch=−𝐾𝑦. (18) Then, the inverse mapping to obtain the gain 𝐾₀ for the control 𝑢(𝑡)=−𝐾₀𝑦(𝑡) is 𝐾₀=𝑉_longΣ⁻¹_long𝐾𝑊. (19)

Table 5 

SVD of longitudinal control distribution matrix^a.

We choose the feedbacks to be 𝑐^∗, 𝑞, and ∫𝑐^∗. This causes the aircraft longitudinal dynamics to have an output vector 𝑦=𝐶𝑥+𝐷𝑢 with a nonzero matrix 𝐷. To apply eigenstructure assignment to these dynamics, we first compute 𝐾₀ from (19) with 𝐷=0. Then, it can be easily shown that the feedback gain matrix 𝐾_𝐷 for nonzero 𝐷 is given by 𝐾_𝐷=𝐼−𝐾₀𝐷⁻¹𝐾₀. (20) We can assign three eigenvalues because we have three outputs. We begin our design with the short-period values that are used byYechout et al.[19] which are 𝜁_sp=0.7 and 𝜔_nsp=3.98. The cstar response is slow, so we change the short-period natural frequency and damping to 𝜁_sp=0.6 and 𝜔_nsp=5.00. This yields an improved cstar response. For the integrator mode, we begin with 𝜆_int=−1 which is approximately the value used by Stevens and Lewis [15].

Velocity and bank angle are not usually fed back for inner-loop control such as a stability augmentation system (SAS) or command augmentation system (CAS). We find that it is difficult to simultaneously achieve short-period and phugoid specifications for the ICE aircraft when feeding back 𝑐^∗, 𝑞, and ∫𝑐^∗. We found that to stabilize the phugoid mode, it is necessary but not sufficient to feed back the angle of attack. However, we do not want to add a fourth output to the eigenstructure assignment because we do not have a fourth eigenvalue to assign. The modes that we are assigning are the short period and the integrator. So we chose to replace 𝑐^∗ with 𝑐^∗+𝛼. We can separate out the 𝑐^∗ and 𝛼 gains from the gains computed by eigenstructure assignment. We further found that the phugoid was unstable even with the 𝛼 feedback when the integrator eigenvalue was assigned to 𝜆_int=−1. We moved the integrator eigenvalue closer to the origin, and the phugoid mode was stabilized with 𝜁_ph=0.053 when 𝜆_int=−0.38. This achieves the phugoid damping specifications, but the gain and phase margins were too small in the 𝛿_pflap, 𝛿_ptv, 𝛼, and ∫𝑐^∗ loops. Finally, we achieved the short-period, phugoid, and gain/phase-margin specifications by assigning the integrator mode at 𝜆_int=−0.25. The phugoid mode exhibits a natural frequency and damping ratio given by (𝜔_nph,𝜁_ph)=(0.01,0.7). The feedback gain matrix and the disk gain/phase margins are shown in Tables 6 and 7, respectively. However, the cstar response was not inside the cstar response envelope [11]. This was because the integrator mode caused the response to be too slow. To remedy this problem, we add a lead compensator in series with the pitch stick. We found that a good cstar response was obtained with the compensator 5(𝑠+0.5)/(𝑠+2). A block diagram of the longitudinal flight control system is shown in Figure 5.

Table 6 

Longitudinal feedback gain matrix 𝐾_𝐷.

Table 7 

Longitudinal gain/phase margins^a.

Figure 5 

Block diagram of longitudinal flight control system.

We use function robuststab from the MATLAB robust control toolbox [20] to compute the margin of stability for real uncertainty in the stability and control derivatives. The uncertainty in the stability derivatives is given by [3] 𝑍_𝛼=2% and 𝑀_𝛼=4%. The uncertainty is 15% in all control derivatives. The uncertain stability and control derivatives also appear in the output equation for 𝑛_𝑧𝑝.

The normal acceleration equation from the numerical linearization is given by 𝑛_𝑧𝑝=0.0048𝑉_𝑇+7.604𝛼−0.3417𝑞−0.0094𝛿_pflap−0.0200𝛿_ptv. (21) However, (21) does not show the dependence of 𝑛_𝑧𝑝 on 𝑍_𝛼, 𝑀_𝛼, 𝑍_pflap, 𝑀_pflap, 𝑍_ptv, or 𝑀_ptv. Thus, we consider the algebraic equation for 𝑛_𝑧𝑝 given by 𝑛_𝑧𝑝≈𝑉_trim𝑔𝛼̇+ℓ𝑞−𝑔𝑞,̇ (22) where 𝑉_trim=422.96 ft/s, 𝑔=32.17 ft/s², and ℓ=17 ft. Upon substituting the 𝛼̇ and 𝑞̇ equations into (22), we obtain 𝑛_𝑧𝑝≈0.00537𝑉_𝑇+−𝑉_trim𝑍_𝛼+ℓ𝑀_𝛼𝑔+−𝑉𝛼−0.3417𝑞_trim𝑍_pflap+ℓ𝑀_pflap𝑔𝛿_pflap+−𝑉_trim𝑍_ptv+ℓ𝑀_ptv𝑔𝛿_ptv. (23) We evaluate (23) at the nominal values of 𝑍_𝛼, 𝑀_𝛼, 𝑍_pflap, 𝑀_pflap, 𝑍_ptv, and 𝑀_ptv which yields 𝑛_𝑧𝑝≈0.00537𝑉_𝑇+7.3988𝛼−0.3417𝑞−0.009814𝛿_pflap−0.01472𝛿_ptv. (24) We observe that the nominal coefficients in (21) and (24) are slightly different.We add the difference in the coefficients to (23), so that the nominal value of the coefficients will be those shown in (21). Then, the equation for 𝑛_𝑧𝑝 that is used in function robuststab is given by 𝑛_𝑧𝑝≈0.0048𝑉_𝑇+−𝑉_trim𝑍_𝛼+ℓ𝑀_𝛼𝑔+−𝑉+0.2052𝛼−0.3417𝑞_trim𝑍_pflap+ℓ𝑀_pflap𝑔+0.000364𝛿_pflap+−𝑉_trim𝑍_ptv+ℓ𝑀_ptv𝑔−0.00528𝛿_ptv. (25)

To obtain a reliable estimate of the robustness margin, we add 2% uncertain dynamics to the real parameter uncertainty [20] and reduce the real uncertainty by 2%. The upper and lower mu bounds are shown in Figure 6 where we observe that the closed loop is robustly stable for up to 178% of the modeled uncertainty.

Figure 6 

Stability margin bounds for longitudinal dynamics.

We analyze performance robustness to the real uncertainty in the stability and control derivatives by performing a monte carlo simulation. We use the linearized longitudinal simulation model that includes first-order actuator dynamics. We perform 100 runs of the linearized longitudinal dynamics over the real uncertainty in the stability and control derivatives. The responses to a one g step command in cstar are shown in Figure 7. We observe that the responses are well behaved for all uncertainty and that deviations from nominal indicate a robust design.

(a)

(b)

(c)

Figure 7 

Monte carlo simulation for longitudinal dynamics.

4. Nonlinear Simulation Results

We perform time responses using the linearized lateral simulation model that includes actuator dynamics and the full 6DOF nonlinear model with deflection and deflection rate limits. The nonlinear simulation includes both the lateral and longitudinal controllers. Figure 8 shows the lateral linearized (dashed line) and the nonlinear (solid line) responses to a 1 deg/s stability axis roll-rate pulse for 0≤𝑡≤2 sec. We observe that (1) sideslip is less than 0.03 deg, (2) steady-state bank angle is 2.3 deg, and (3) steady-state stability axis yaw rate is 0.17 deg/s. We also observe that the linear and nonlinear responses are almost identical. Figure 9 shows the control effector deflections. We observe that (1) maximum elevon deflection is 0.7 deg, (2) maximum amt deflection is 0.3 deg, and (3) maximum yaw thrust vector deflection is 0.2 deg. We conclude that the nonlinear responses are excellent. Next, we perform time responses using the linearized longitudinal simulation model that includes actuator dynamics and the full 6DOF nonlinear model with deflection and deflection rate limits. The nonlinear simulation includes both the lateral and longitudinal controllers. Figure 10 shows the longitudinal linearized (dashed line) and the nonlinear (solid line) responses to a 1 g cstar step. Normal acceleration and pitch rate are shown for 0≤𝑡≤10 sec, whereas normalized cstar is shown for 0≤𝑡≤3 sec. This is because the normalized cstar envelope is shown in [11] only for 0≤𝑡≤3 sec. We observe that (1) the normalized cstar response is within the category I flight envelope of [11] and (2) the linear and nonlinear responses are similar until approximately 5 seconds when the nonlinear normal acceleration tends to its one g equilibrium value. The maximum pitch flap and pitch thrust vector deflections are 4.5 deg and 2.2 deg, respectively, for a maximum of 1.4 g normal acceleration. These deflections are well below the maximum deflections and will allow for much larger values of normal acceleration.

(a)

(b)

(c)

(d)

Figure 8 

Nonlinear time responses to a 1 deg/s 𝑝s pulse.

(a)

(b)

(c)

(d)

(e)

Figure 9 

Nonlinear control responses to a 1 deg/s 𝑝s pulse.

(a)

(b)

Figure 10 

Nonlinear time responses to a 1 g cstar step.

5. Summary

We have designed a flight control system for the ICE tailless aircraft that includes thrust vectoring. Our longitudinal controller is a cstar command system, and our lateral controller is a stability axis roll-rate command system. Our design achieves MIL-F-8785C mode specifications on the short period, phugoid, dutch roll, roll subsidence, and spiral modes and achieves MIL-F-9490D specifications on phase and gain margins. The robustness of the design is further demonstrated by closed-loop margins based upon mu values and monte carlo simulations. The performance of the flight control system is evaluated by using a 6DOF nonlinear simulation of the ICE aircraft. Nonlinear simulation results show good performance to cstar and stability axis roll-rate commands.

Acknowledgment

The work of the first author was partially supported by The Alliance for Graduate Education and the Professoriate (AGEP) and the New York City Louis Stokes Alliance for Minority Participation (NYC-LSAMP).