Adaptive Flight Control in the Presence of Limits on Magnitude and Rate
Joseph E. Gaudio Anuradha M. Annaswamy Michael A. Bolender Eugene Lavretsky
math.OC cs.SY eess.SY













































































Adaptive Flight Control in the Presence of Limits on Magnitude and Rate
    Joseph E. Gaudio, Anuradha M. Annaswamy, Michael A. Bolender, and Eugene Lavretsky
J.E. Gaudio and A.M. Annaswamy are with the Department
of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge,
MA, 02139 USA (email: jegaudio@mit.edu and aanna@mit.edu).
M.A. Bolender is with the Air Force Research Laboratory, WPAFB, OH, 45433 USA.
E. Lavretsky is with The Boeing Company, Huntington Beach, CA, 92647 USA (email: eugene.lavretsky@boeing.com).


    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================







Input constraints as well as parametric uncertainties must be accounted for in the design of safe control systems. This paper presents an adaptive controller for multiple-input-multiple-output (MIMO) plants with input magnitude and rate saturation in the presence of parametric uncertainties. A filter is introduced in the control path to accommodate the presence of rate limits. An output feedback adaptive controller is designed to stabilize the closed loop system even in the presence of this filter. The overall control architecture includes adaptive laws that are modified to account for the magnitude and rate limits. Analytical guarantees of bounded solutions and satisfactory tracking are provided. Three flight control simulations with nonlinear models of the aircraft dynamics are provided to demonstrate the efficacy of the proposed adaptive controller for open loop stable and unstable systems in the presence of uncertainties in the dynamics as well as input magnitude and rate saturation.




Adaptive control, magnitude saturation, rate saturation, flight control





















§ INTRODUCTION


















Any advanced control system must be capable of incorporating realistic constraints on control inputs such as magnitude limits and rate limits. In flight control, such constraints are commonly present due to actuator limits. While magnitude saturation is often accounted for in the underlying control design, the nonlinearity arising from input rate saturation is often ignored. Actuator rate saturation can lead to Pilot Induced Oscillations (PIO), which expose the aircraft to the risk of failures, and in worst cases, to departing flight <cit.>. Crashes in the SAAB Grippen development are evidences of the latter <cit.>. Additionally, in the event of control surface damage, additional complexities may arise if rate saturation is dominant.

Over the past four decades, adaptive control of linear time-invariant plant models in the presence of parametric uncertainties, perturbations due to bounded disturbances, and unmodeled dynamics has been studied extensively <cit.>. Adaptive control in the presence of magnitude constraints has been addressed in <cit.>. With the first set of results on this topic reported in <cit.>, references <cit.> extended the analysis to the multiple-input state feedback case. Design of a state feedback magnitude saturation control architecture with a buffer region was shown in <cit.>, by modifying the reference model. Reference <cit.> presents an application dependent architecture. None of these references, however, address input rate saturation.

Rate saturation architectures have been considered in <cit.>. Reference <cit.> proposes directly differentiating the control signal in order to saturate the control rate before re-integrating the signal. References <cit.> also propose a non-adaptive, robustness based architectures. These non-adaptive architectures do not directly compensate for plant parametric uncertainty. Reference <cit.> proposes a state feedback indirect adaptive control architecture, while <cit.> proposes a direct model reference control architecture, although a matching condition is violated. The integrator anti-windup architecture presented in <cit.> is proposed for systems with input saturation, but does not include a proof of stability for an adaptive system and indeed may result in instability for certain adaptive systems <cit.>. The controller presented in this paper is significantly less restrictive than these papers in that it provides a solution for the MIMO case with guarantees of stability. It is assumed that the control input is subjected to a hard limit on its rate and therefore differs from our earlier work in <cit.> which imposed only a soft limit, i.e. the control rate was allowed to exceed the specified limits.

The problem addressed in this paper is the control and command tracking of plant dynamics in the presence of parametric uncertainties, using input-output measurements, even when the plant input is subjected to hard limits on its magnitude and rate. In order to introduce a rate limit on the control input without explicit differentiation, a filter is designed with hard saturation nonlinearities, similar to <cit.>. The introduction of such a filter however causes new challenges in the form of an additional lag. This in turn causes the underlying plant to have an increase in relative degree. This, and the fact that we only have outputs available for feedback, suggests the use of a MIMO adaptive controller that uses output feedback and that is applicable for a plant with higher relative degree. For this purpose, we propose an adaptive control approach that is along the lines of <cit.>. While the results of these papers guarantee stability, they do not consider either a magnitude or a rate-limit (and will be shown empirically to result in instability in this paper). The contribution of this paper is that a similar controller as in <cit.>, augmented to account for magnitude and rate limits, can still be shown to lead to bounded solutions even in the presence of these limits.

This paper proceeds as follows: Section <ref> presents mathematical preliminaries. Section <ref> describes the problem formulation of the output feedback control a plant with magnitude and rate saturation in the presence of parametric uncertainties. The adaptive control architecture is presented in Section <ref>. Stability analysis of the closed-loop system follows in Section <ref>. Section <ref> demonstrates the efficacy of the proposed controller in numerical simulations with fighter and hypersonic aircraft, with conclusions following in Section <ref>.














§ PRELIMINARIES

The following notation is used for a MIMO plant model with m inputs and p outputs: {A,B,C}:=C(sI-A)^-1B, where s=d/dt is the differential operator. The transmission zeros of this system are defined using Definition 1 of reference <cit.>. The columns of the matrix B may be partitioned as B=[b_1,b_2,⋯,b_m]. The input relative degree of the MIMO plant model is stated as follows.
<cit.>
	A linear square plant given by {A,B,C} has
	
		
  * input relative degree r=[r_1,r_2,...,r_m]^T∈ℕ^m× 1 if
		
			
  * ∀ j∈{1,...,m}, ∀ k∈{0,...,r_j-2}: CA^kb_j=0_p× 1, and
			
  * rank[CA^r_1-1b_1   CA^r_2-1b_2...CA^r_m-1b_m]=m
		
		
  * uniform input relative degree r∈ℕ if it has input relative degree r=[r_1,r_2,...,r_m]^T and r=r_1=r_2=...=r_m
	

Every MIMO plant model thus has an input relative degree. The following proposition describes a minimal realization where there are differentiators on the plant input.
<cit.>
	
	Given a linear system {A,B_2,C} with uniform input relative degree 2, the following two realizations are equivalent:
	
		
  * ẋ(t)=Ax(t)+B_2(a_1^1s+a_1^0)[u(t)],   y(t)=Cx(t)
		
  * ẋ'(t)=Ax'(t)+B_2^1u(t),   y(t)=Cx'(t)
	
where 

    x'(t)=x(t)-B_2a_1^1u(t)

	
    B_2^1 = AB_2a_1^1+B_2a_1^0


The following identity holds: (sI-A)^-1s=I+(sI-A)^-1A. Thus for a relative degree two input: C(sI-A)^-1B_2s=C(sI-A)^-1AB_2 as CB_2 = 0. In this paper elliptical saturation is considered.
<cit.>
An elliptical saturation function of a vector v(t) is defined as

    E_s(v(t),v_max) = 
    
    v(t)      ||v(t)||<g(v(t))
    v̅(t)      ||v(t)||>g(v(t))

where the function g(v(t)) is expressed as

    g(v(t))=(∑_i=1^m[ê_i/(v_max)_i]^2)^-1/2

where ê=v/||v|| and v̅=êg(v) (see Figure <ref>).

It can be noted that this saturation function can be alternatively implemented using the projection operator as in <cit.>.










§ PROBLEM FORMULATION

We consider a class of linear plants of the form

    ẋ_p(t)   =A_px_p(t)+B_pΛ^*[u_p(t)+Θ_p^*Tx_p(t)]
    
    		y_p(t)   =C_px_p(t)
    
    		z(t)   =C_pzx_p(t)+D_pzΛ^*[u_p(t)+Θ_p^*Tx_p(t)]

with known matrices of dimensions: A_p∈ℝ^n_p× n_p, B_p∈ℝ^n_p× m, C∈ℝ^p_p× n_p, C_pz∈ℝ^n_z× n_p, and D_pz∈ℝ^n_z× m, n_p≥ p_p≥ m≥ n_z. The diagonal matrix Λ^*∈ℝ^m× m, and Θ_p^*∈ℝ^n_p× m represent unknown constant matched uncertainty which enter the plant dynamics through the columns of B_p. These uncertainty locations arise in a variety of applications, including in aircraft dynamics <cit.>. x_p, u_p, y_p, and z denote the plant state, plant input, plant output, and regulated output respectively. The following assumptions are made of the plant in (<ref>):

	{A_p,B_p,C_p} is a minimal realization;


	All transmission zeros of {A_p,B_p,C_p} and {A_p,B_p,C_pz,D_pz} are stable;


	{A_p,B_p,C_p} is uniform relative degree one;


	The uncertainty Θ_p^* is bounded by a known value, i.e. ||Θ_p^*||<Θ_p,max;


	The uncertainty Λ^* is diagonal positive definite and bounded by a known value along with its inverse, i.e. ||Λ^*||<Λ_max, ||Λ^*-1||<Λ_inv,max.

Assumptions <ref> and <ref> are standard in output feedback adaptive control <cit.>, where the condition of stable transmission zeros of {A_p,B_p,C_pz,D_pz} is employed to ensure controllability of an extended plant in the presence of integral tracking (to be defined) <cit.>. The relative degree statement in Assumption <ref> is commonly satisfied for aircraft control systems. The additive uncertainty Θ_p^* is assumed to be bounded by a known value in Assumption <ref>. Assumption <ref> states that the uncertain control effectiveness of each input path are independent of one another, upper bounded by a known value and bounded away from zero by a known value. Known bounds on the norm of the system uncertainty as in Assumptions <ref> and <ref> are required given that constraints will be imposed on the plant input.

The goal is to design u_p so that z tracks a bounded command z_cmd. To ensure a small tracking error, an integral error state x_e is generated as <cit.>

    ẋ_e(t)=z(t)-z_cmd(t).




 §.§ Rate Saturation

The goal of this paper is to design an adaptive controller that will accommodate the parametric uncertainties in (<ref>) and carry out the desired tracking using a control input which is rate limited, and possibly magnitude limited as well. If the time derivative of the computed control input u were to be available, a filtered control rate u_r could be generated as

    τu̇_r(t) + u_r(t) = u̇(t)

where u is the output of an adaptive controller. Equation (<ref>) implies that for a small enough filter time constant τ, u_r(t)≈u̇(t). One could then recover the actual output of the controller u by simply integrating u_r and setting it to be equal to the control input u_p as

    u̇_p(t)=u_r(t).

It should be noted that u_r in (<ref>) can be realized without explicit differentiation of u, as

    u_r(t)=1/τu(t)-1/τ[1/τ s+1]u(t)

where s denotes the operator d/dt. Equations (<ref>), (<ref>) provide us with a method for generating a control input u_p(t) that is rate-limited by simply replacing u_r with a saturation function of u_r, and generate u_p(t) as

    τu̇_r(t) + E_s(u_r(t),u_r,max)    = u̇(t)
    u̇_p(t)   =E_s(u_r(t),u_r,max)

instead of equations (<ref>) and (<ref>), where u_r,max is the desired rate limit on the control input. As in (<ref>), we note that u_r can be realized without using explicit differentiation as

    u_r(t)=1/τ(u(t)-u_p(t))

which can be used to generate the approximate derivative u_r rather than the differential form (<ref>).




 §.§ Magnitude Saturation

In addition to rate saturation, it is easy to ensure that the control input u_p is magnitude limited as well. For this purpose, rather than (<ref>), we generate u_r as

    u_r(t)=1/τ(E_s(u(t),u_max)-u_p(t))

where u_max denotes the magnitude limit of u. Equations (<ref>) and (<ref>) lead to a realization of a control input u_p into the plant that is magnitude limited and rate limited, starting from the output u from the adaptive controller (see Figure <ref> for a schematic).








 §.§ Saturation Effects as Disturbances

The saturation functions in magnitude and rate can be viewed as two nonlinearities. As in <cit.>, we accommodate these nonlinearities by treating their impact as additive known disturbances. In particular defining two known disturbance terms Δ u and Δ u_r as

    Δ u(t)   =E_s(u(t),u_max)-u(t)
    Δ u_r(t)   =E_s(u_r(t),u_r,max)-u_r(t)

it is easy to see that if u does not reach its magnitude saturation limit u_max, then Δ u(t)≡ 0. Similarly, if the input rate u_r does not reach its rate saturation limit u_r, max, then Δ u_r(t)≡ 0, that is, these known disturbance terms become non-zero only if the magnitude or rate limits are exceeded.

Combining (<ref>), (<ref>), and (<ref>) results in the following compact relation between the plant input and the control input:

    u̇_p(t)=-1/τu_p(t)+1/τu(t)+1/τΔ u_2(t)

where Δ u_2(t)=(Δ u(t)+τΔ u_r(t)) represents the combined effects of magnitude and rate saturation (see Figure <ref> for a schematic). That is, equation (<ref>) determines the relation between u_p, the plant input, and u, the computed control input in a compact manner.




 §.§ Full Plant Equations

We now assemble the complete plant model that includes the plant dynamics, the integral of the tracking error, and the effects of rate and magnitude saturation, which are given by equations (<ref>), (<ref>), and (<ref>), respectively. This is given by

    2[ ẋ_p(t); ẇ_u(t); ẋ_e(t) ]_ẋ(t)   =      [   A_p   B_p     0;     0 -1/τI     0;  C_pz  D_pz     0 ]_A[ x_p(t); w_u(t); x_e(t) ]_x(t)
             +
    [  B_p;    0; D_pz;      ]_B_1Λ^*Θ^*T_px_p(t)
    +
    [    0; 1/τI;    0 ]_B_2Λ^* u(t)
    
             +
    [    0; 1/τI;    0 ]Λ^*Δ u_2(t)
    +
    [  0;  0; -I;    ]_B_z
    z_cmd(t)
    
    
    y(t)   =       [ C_p   0   0;   0   0   I ]_C[ x_p(t); w_u(t); x_e(t) ]
    
    z(t)   =       [ C_pz D_pz    0 ]_C_z[ x_p(t); w_u(t); x_e(t) ]
    +
    D_pzΛ^*Θ^*T_px_p(t).

where w_u:=Λ^*u_p. Equation (<ref>) can be expressed in a compact form as

    2ẋ(t)   =       Ax(t)+B_1Ψ^*T_1x(t)+B_2Λ^*u(t)
             +B_2Λ^*Δ u_2(t)+B_zz_cmd(t)
    
    y(t)   =       Cx(t)
    
    z(t)   =       C_zx(t)+D_pzΨ^*T_1x(t)

where A∈ℝ^n× n, B_1∈ℝ^n× m, B_2∈ℝ^n× m, B_z∈ℝ^n× n_z, C∈ℝ^p× n, C_z∈ℝ^n_z× n. The matrix of additive uncertainty is represented as Ψ^*T_1=
[ Λ^*Θ^*T_p         0         0 ] and the measured output is y.

Based on the structure of the full plant model in (<ref>) it can be noted that {A,B_2,C} is a minimal realization given Assumptions <ref> and <ref>. In particular Assumption <ref> ensures that the inclusion of integral tracking (<ref>) in the extended dynamics (<ref>) preserves controllability given that {A_p,B_p,C_pz,D_pz} does not have a transmission zero at the origin <cit.>. Observability of the extended dynamics (<ref>) is preserved (and may also be enhanced) with the inclusion of integral tracking (<ref>) as the integral tracking state x_e is measured. Furthermore it can be demonstrated that the inclusion of the filter dynamics with state w_u in (<ref>) maintains controllability and observability of the extended system, given the form of a low-pass filter in series with a controllable and observable plant with states x_p and x_e.

Additionally, given Assumptions <ref> and <ref>, all transmission zeros of {A,B_2,C} are stable and {A,B_2,C} is uniform relative degree two. It can be noticed that B_1 can be spanned by a linear combination of B_2 and AB_2 as B_1=τ AB_2+B_2. Furthermore Ψ^*_1 satisfies Ψ^*T_1B_2=0 (Ψ^*_1 is not in the same input path as the computed control input) and is bounded by a known value from Assumptions <ref> and <ref>, i.e. ||Ψ^*_1||<Ψ_max=Θ_p,maxΛ_max.



§ ADAPTIVE CONTROL DESIGN

This section presents an adaptive controller for the plant model of uniform relative degree two in equation (<ref>) whose structure is similar to that in <cit.>. As will be shown, modifications in the underlying adaptive laws are however needed to account for the effects of input magnitude and rate saturation. Section <ref> presents the control architecture with gains designed in Section <ref>. The construction of an underlying error model that is strictly positive real is presented in Section <ref>.




 §.§ Control Architecture

As in any adaptive control design, we begin the control architecture discussion with the introduction of a reference model, which denotes the desired dynamics from the plant when there is no parametric uncertainty and no saturation disturbance. This corresponds to the plant model in (<ref>) with Λ^*=I, Ψ_1^*=0, and Δ u_2(t)=0. A closed-loop term is added, as in recent investigations <cit.>, which leads to a closed-loop reference model (CRM) in order to ensure smooth control inputs. Given that this problem is in output feedback, the reference model additionally serves the dual purpose as an observer <cit.>. The use of x_m to denote the state of the CRM as compared to the common observer notation of x̂ is solely a notation choice, where both representations have appeared in the adaptive control literature. Based on the plant structure (<ref>) in Section <ref>, a CRM may be chosen as

    ẋ_m(t)   =Ax_m(t)+B_2u_bl(t)+L(y(t)-y_m(t))+B_zz_cmd(t)
    
    y_m(t)   =Cx_m(t)
    
    z_m(t)   =C_zx_m(t)+D_pzu_bl(t)

where the design of L is discussed in Section <ref>. The input u_bl is chosen so as to ensure a stable CRM, and is of the form

    u_bl(t)=-Kx_m(t)

where the matrix K∈ℝ^n× m can be designed using linear quadratic regulator (LQR) techniques <cit.> to provide for a desired reference model dynamics matrix A_m=A-B_2K. It can be noticed that in the presence of perfect output tracking of the reference model output (y(t)-y_m(t)=0), the CRM in (<ref>) has the form of the standard adaptive control open loop reference model with ẋ_m(t)=A_mx_m(t)+B_zz_cmd(t).

Given that the plant in (<ref>) contains parametric uncertainties, an adaptive control input will be used to counter these uncertainties. The difficulty however is that the introduction of the filter in (<ref>) causes the underlying plant dynamics to have a relative degree of two, given that the model in (<ref>) has a relative degree of one as noted in Assumption <ref>. This causes the corresponding error model, derived by subtracting (<ref>) from (<ref>), to have a relative degree two property as well. The relative degree two property prohibits the use of a typical adaptive control input of the form of u(t)=Ω^T(t)ξ(x_m(t)), where Ω is an estimate of unknown parameters and ξ is the regressor. In order to provide for a strictly positive real (SPR) error model, an extra zero must be added as discussed in <cit.> Chapter 5. For this purpose, we now design the computed control input u as <cit.>

    u(t)=(a^1_1s+a^0_1)Ω^T(t)ξ̅(t)






where the variables a^1_1>0 and a^0_1>0 can be selected as desired to place the extra zero. The adaptive parameter matrix Ω and regressor vector ξ̅ are defined as

    Ω(t)   =[Λ^T(t), Ψ_1^T(t), Ψ_2^T(t)]^T
    ξ̅(t)   =[u̅_bl^T(t), -x_m^T(t), -x̅_m^T(t)]^T

where Ω is an estimate of Ω^*=[Λ^*-1, Ψ̅^*T_1, Ψ̅^*T_2]^T, and Ψ̅_1^*=
[ τ/a_1^1Θ_p^*T             0             0 ]^T and Ψ̅_2^*=
[ (I-τa_1^0/a_1^1I)Θ_p^*T                       0                       0 ]^T.
The following filtered signals are defined as

    u̅_bl(t)   =1/a_1^1s+a_1^0u_bl(t)
    Δu̅_2(t)   =1/a_1^1s+a_1^0Δ u_2(t)
    x̅_m(t)   =1/a_1^1s+a_1^0x_m(t)

Since the elements of ξ̅ are filtered signals, it is easy to show that (<ref>) can be realized without explicit differentiation, provided the derivative of Ω can be directly synthesized. We will show in the following that this is indeed the case.








 §.§ The Design of L

We now address the design of L, which is accomplished so as to ensure an underlying transfer function to be SPR. Since we begin with a non-square plant, a lemma is needed to square-up {A,B_2,C} in order to apply the KYP lemma (see <cit.> Lemma 2.5). See <cit.> for a proof and design procedure of the contents of Lemma <ref>.
<cit.>
	
	For plant models of the form of (<ref>) which satisfy Assumptions <ref> to <ref>, there exists a matrix B_s1∈ℝ^n×(p-m) such that {A,B̅_2,C}, where B̅_2=[B_2,B_s1], has stable transmission zeros and nonuniform input relative degree r_i=2 for i=1,2,...,m and r_i=1 for i=m+1,m+2,...,p.

The matrix B_s1 is used only in the design of L as follows:

    B̅_2^1=[B_2^1,B_s1]


    S^T=(CB̅_2^1)=[S_2^T,S_1^T]


    C̅=SC


    R^-1(ε)=(C̅B̅_2^1)^-1[C̅AB̅_2^1+(C̅AB̅_2^1)^T](C̅B̅_2^1)^-1+ε I


    ε>ε_max(A,B,C̅,a_1^1,a_1^0,Λ_max,Ψ_max)


    L=B̅_2^1R^-1(ε)S


    A^*_L=(A+B_1Ψ^*T_1-LC).

where ε>ε_max is selected large enough to guarantee that A_L^* is Hurwitz and the realization {A^*_L,B,SC} is SPR. In summary, the closed-loop gain L is chosen so as to guarantee that {A^*_L,B,SC} is SPR. This property in turn will be used to derive the adaptive laws for adjusting Ω, which is shown in the next section.








 §.§ SPR Error Model

This section will derive an error model and propose adaptive update laws from its SPR properties. The model tracking error is defined as: e_x(t)=x(t)-x_m(t). Applying (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>), the following error model may be derived:

    2ė_x(t)   =       A^*_Le_x(t)+B_2^1Λ^*Ψ̅^*T_1x_m(t)
             +B_2Λ^*(a^1_1s+a^0_1)Ψ̅^*T_2x̅_m(t)
             -B_2Λ^*(a^1_1s+a^0_1)Λ^*-1u̅_bl(t)
             +B_2Λ^*(a^1_1s+a^0_1)Ω^T(t)ξ̅(t)
             +B_2Λ^*(a^1_1s+a^0_1)Δu̅_2(t).

Given the undesirable differentiators in the right hand side of this equation, Proposition <ref> may be applied with the modified model tracking error variable of the form

    2
    e_mx(t)   =       e_x(t)-B_2Λ^*a_1^1[Ψ̅^*T_2x̅_m(t)+Ω^T(t)ξ̅(t).
             .-Λ^*-1u̅_bl(t)+Δu̅_2(t)].

This results in a modified model tracking error model as

    ė_mx(t)   =A^*_Le_mx(t)+B_2^1Λ^*Ω̃^T(t)ξ̅(t)+ B_2^1Λ^*Δu̅_2(t)
    
    e_y(t)   =Ce_mx(t)=Ce_x(t)

where Ω̃(t)=Ω(t)-Ω^*.

It should be noted that equation (<ref>) cannot be used directly to determine the rules for adjusting the adaptive parameters Ω(t). Unlike the approaches in <cit.>, there are additional terms Δu̅_2(t) present here due to the presence of magnitude and rate saturation. This introduces a significant departure from the procedure in <cit.> and requires additional tools. It should also be noted that the methods proposed in <cit.> are not sufficient in determining adaptive laws for the stability of the overall system. How these difficulties are overcome in this MIMO setting correspond to the main contributions of this paper and are elaborated on below.

In order to eliminate the effects of Δu̅_2, an implementable auxiliary signal e_Δ is introduced as

    ė_Δ(t)   =A_Le_Δ(t)+B_2^1Ω_Δ^T(t)ξ̅_Δ(t)
    
    e_y,Δ(t)   =Ce_Δ(t).

The matrix Ω_Δ and vector ξ̅_Δ are

    Ω_Δ(t)   =[Λ̂^T(t), Ψ̂_1^T(t), Ψ̂_2^T(t)]^T
    ξ̅_Δ(t)   =[Δu̅_2^T(t), e_Δ^T(t), e̅_Δ^T(t)]^T

where Ω_Δ is an estimate of Ω_Δ^*=[Λ^*, Λ^*Ψ̅^*T_1, Λ^*Ψ̅^*T_2]^T, and where the following filtered signal is defined:

    e̅_Δ(t)=1/a_1^1s+a_1^0e_Δ(t).

The following is an equivalent representation of (<ref>), which is obtained using (<ref>) and (<ref>):

    2ė_Δ(t)   =       A^*_Le_Δ(t)-B_2^1Λ^*Ψ̅_1^*Te_Δ(t)
             -B_2(a_1^1s+a_1^0)Λ^*Ψ̅_2^*Te̅_Δ(t)+B_2^1Ω_Δ^T(t)ξ̅_Δ(t).

To eliminate the differentiator on the right hand side of this equation, Proposition <ref> can be once again applied with a modified auxiliary signal of the form

    e_mΔ(t)=e_Δ(t)+B_2a_1^1[Λ^*Ψ̅_2^*Te̅_Δ(t)].

This results in modified dynamics of the form

    ė_mΔ(t)   =A^*_Le_mΔ(t)+B_2^1Ω̃_Δ^T(t)ξ̅_Δ(t)+B_2^1Λ^*Δu̅_2(t)
    
    e_y,Δ(t)   =Ce_mΔ(t)=Ce_Δ(t)

where Ω̃_Δ(t)=Ω_Δ(t)-Ω_Δ^*.
To remove the effects of Δu̅_2 the following modified augmented error (e_mu(t)=e_mx(t)-e_mΔ(t)) model is introduced:

    ė_mu(t)   =A^*_Le_mu(t)+B_2^1Λ^*Ω̃(t)ξ̅(t)-B_2^1Ω̃_Δ^T(t)ξ̅_Δ(t)
    
    e_y,u(t)   =e_y(t)-e_y,Δ(t)

We note that e_y,u(t) is available at each t since e_y(t) is measurable and e_y,Δ(t) is a known signal that can be computed at each time t.
<cit.>
	
	Given L∈ℝ^n× p in (<ref>) and S∈ℝ^p× p in (<ref>), ε in (<ref>), and with Assumptions <ref> to <ref>, for plant models of the form of (<ref>), the transfer function {A^*_L,B̅_2^1,SC}  is strictly positive real.

Lemma <ref> also implies that {A^*_L,B_2^1,S_2C} is SPR given the partitions in (<ref>) and (<ref>). Thus given the structure of the modified augmented error model in (<ref>), the adaptive parameters Ω and Ω_Δ can be updated as

    Ω̇(t)=-Γ_Ωξ̅(t)e_y,u^T(t)S_2^T


    Ω̇_Δ(t)=Γ_Ω_Δξ̅_Δ(t)e_y,u^T(t)S_2^T

where Γ_Ω>0 and Γ_Ω_Δ>0 are adaptive update gains.

A few comments regarding the two update laws are in order. Equation (<ref>) represents a standard update to address parametric uncertainties in Ω^*. Equation (<ref>) represents an update due to the effects of saturation, since once the update law and e_Δ(0) are initialized at zero, they will only change when a saturation limit is reached. The adaptive update laws in (<ref>) and (<ref>) along with the fact that the realization {(A^*-LC,B_2^1,S_2C)} is SPR provides the foundation for the closed-loop system stability. This is established in the next section.









§ STABILITY ANALYSIS


Before proceeding to state the main stability result, a discussion of boundedness of the adaptive parameters is presented in Section <ref>. The boundedness of all remaining states is addressed in Section <ref>. It can be noted that due to the presence of saturation, not all signals can be tracked and thus results are local in nature. Bounds on the local stability results will be shown to be proportional to the level of saturation and uncertainty.



 §.§ Boundedness of Adaptive Parameters


We consider the following Lyapunov function candidate:

    2
    V(e_mu(t), Ω̃(t), Ω̃_Δ(t))   =       e_mu^T(t)Pe_mu(t)
             +Tr[Ω̃^T(t)Γ^-1_ΩΩ̃(t)|Λ^*|]
             +Tr[Ω̃_Δ^T(t)Γ^-1_Ω_ΔΩ̃_Δ(t)]

where P=P^T>0 is a positive definite matrix which can be used to guarantee the SPR properties of {(A^*_L,B̅_2^1,SC)} and satisfy

    A^*T_LP+PA^*_L   =-Q<0
    
    PB̅_2^1   =C^TS^T

for a positive definite matrix Q=Q^T>0. The following partition may be used:

    P[B_2^1,B_s1]=C^T[S_2^T,S_1^T].

Taking a derivative of (<ref>) with respect to time:

    2V̇(e_mu(t), Ω̃(t), Ω̃_Δ(t))   =       e_mu^T(t)(A^*T_LP+PA^*_L)e_mu(t)
             +2e_mu^T(t)PB_2^1Λ^*Ω̃(t)ξ̅(t)
             -2e_mu^T(t)PB_2^1Ω̃_Δ(t)ξ̅_Δ(t)
             +2Tr[Ω̃^T(t)Γ^-1_ΩΩ̇(t)|Λ^*|]
             +2Tr[Ω̃_Δ^T(t)Γ^-1_Ω_ΔΩ̇_Δ(t)].

Applying (<ref>), (<ref>), (<ref>), and (<ref>):

    V̇(e_mu(t), Ω̃(t), Ω̃_Δ(t))=-e_mu^T(t)Qe_mu(t)≤0.

Thus all of (e_mu, Ω̃, Ω̃_Δ) are bounded. This does not show that the model tracking error e_x is bounded, however the following relation for the tracking error can be obtained:

    ||e_x(t)||=𝒪[sup_ζ≤ t||Δ u_2(ζ)||].

Such a relation is similar to that derived in <cit.>, where the constants in the big 𝒪 notation are system parameters.









 §.§ Boundedness of the State


For clarity of presentation, the constant matrices, constant scalars and time-varying scalars found in this section are defined in Appendix <ref>. The following closed-loop system representation can be obtained through combining (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>):

    2χ̇(t)   =       A_clχ(t)-B_ΩΩ̃^T(t)C_ξ̅χ(t)
             +C_Δu̅_21/a_1^1Δ u_2(t)+B_Zz_cmd(t)

where the full state of the system to be shown bounded is:

    χ(t)=
    [  x̅_m^T(t)   x_m^T(t)  e_mx^T(t) Δu̅_2^T(t) ]^T.

It can be noted from the block upper triangular structure of A_cl in Appendix <ref>, that its eigenvalues can be placed arbitrarily stable. This is due to the controllability properties of the underlying dynamics and the use of a closed-loop reference model (<ref>). Showing boundedness of the state χ will result in boundedness of all of the signals in the closed-loop system given that e_mu, Ω̃, and Ω̃_Δ were shown to be bounded in Section <ref>. The main challenges that arise here are due to the various combinations of scenarios of magnitude and rate saturations that can occur, which causes Δu̅_2 to be non-zero. This in turn necessitates the use of multiple sub-cases that are discussed in the proof of Theorem <ref> in Appendix <ref>. Before proceeding to a discussion of the main stability theorem, given that the system is input constrained, the following reference command bound is required.

The reference command z_cmd is bounded as ||z_cmd(t)||≤ z_cmd,max and is chosen such that ρχ_min<χ_max, κ_3>0, κ_8>0, κ_5^2>4κ_4κ_6, κ_5^2>4κ_9κ_10.

Assumption <ref> implies that the magnitude of the reference command is upper bounded with respect to the level of saturation and system uncertainty. The larger the level of saturation and smaller the system uncertainty, the larger the allowable command. The bound ρχ_min<χ_max can always be achieved for a provably bounded command, as for example if z_cmd,max=0: χ_min=0 and χ_max is finite, with κ_3>0,κ_8>0 and κ_5^2>4κ_4κ_6, κ_5^2>4κ_9κ_10. Thus χ_min can be seen to represent a shift in the equilibrium of the system due to tracking of z_cmd. We now state the main theorem of stability.


	
	Under Assumption <ref> for the system in (<ref>) and (<ref>), control input (<ref>), known control disturbances (<ref>), adaptive laws in (<ref>), (<ref>), and Lyapunov function (<ref>), the closed-loop system has bounded trajectories for all t≥0 if the following two conditions are satisfied: 
	
		
  * |χ(0)|<1/ρχ_max
		
  * √(V(0))<√(λ_min/γ_max)Ω̅_max
	
	
	Further
	
	
    ||χ(t)||<χ_max,  ∀ t>0

	And
	
    ||e_x(t)||= 𝒪[sup_ζ≤ t||Δ u_2(ζ)||]



Given that the dynamics matrix A_cl can be chosen by design to be strictly stable as in Appendix <ref>, it satisfies the following:

    A_cl^TP_cl+P_clA_cl<-Q_cl

where the matrix P_cl may be computed using a positive definite matrix Q_cl. Define the following candidate Lyapunov function of the closed-loop dynamics:

    W(χ(t))=χ^T(t)P_clχ(t).

Define the following level set, ℬ, of W:

    ℬ:{χ(t)|W(χ(t))=p_minχ_max^2}.

The following annulus region is defined:

    𝒜:{χ(t)|χ_min≤ ||χ(t)||≤χ_max}.

The proof of boundedness of the full system state follows from two steps. The first step will show that ℬ⊂𝒜 using condition (ii) of Theorem 1. Step 2 will show that Ẇ(χ(t))<0, ∀χ∈𝒜. Condition (i) in Theorem 1 implies the following:

    W(χ(0))<W(ℬ).

The two steps thus result in the following:

    W(χ(t))≤ W(χ(0))  ∀ t≥0.

Theorem 1 follows from these two steps.

Proof Step 1: This step shows that ℬ⊂𝒜. Using condition (ii) from Theorem 1, it can be seen that Ω̃_max<Ω̅_max. Additionally by Assumption <ref>:

    ρχ_min<χ_max.

Equation (<ref>) can be used to show that W(χ(t)) is bounded from below by p_min||χ(t)||^2≤ W(χ(t)). This implies the following:

    ||χ(t)||≤χ_max  ∀χ∈ℬ.

In a similar manner from equation (<ref>), W(χ(t)) can be bounded from above by W(χ(t))≤ p_max||χ(t)||^2. These relations imply the following:

    χ_min<1/ρχ_max≤||χ(t)||,   ∀χ∈ℬ.

The definition of the annulus region 𝒜 can then be used to conclude that ℬ⊂𝒜.

Proof Step 2: It can then be shown that Ẇ(χ(t))< 0 ∀χ∈𝒜. Three cases are considered. Case 1 considers the system not in magnitude nor rate saturation (Δ u=0,Δ u_r=0). Case 2 considers the system in magnitude saturation (Δ u≠0,Δ u_r=0) and case 3 considers the system in rate saturation (Δ u=0,Δ u_r≠0). The proof of this step is located in Appendix <ref>.

Boundedness of the state χ in equation (<ref>) as well as boundedness of e_mu, Ω̃, and Ω̃_Δ in Section <ref> is sufficient to show that all of the signals in the closed-loop system remain bounded, thus concluding the proof of Theorem <ref>.




	Global stability cannot be achieved in the presence of magnitude and rate saturation for general (possibly open loop unstable) plant models. There always exists initial conditions that will cause the system to have unbounded trajectories regardless of the design of the controller. The results are thus dependent on formulations of regions of attraction, dependent on the saturation levels and amount of uncertainty, such as those presented in Theorem <ref>.


	If the plant dynamics are open-loop stable, equation (<ref>) is bounded-input, bounded-output (BIBO) stable and equations (<ref>) and (<ref>) results in a bounded control input. Thus the state trajectory χ is bounded regardless of initial condition.


	This architecture incorporates the effects of saturation in the auxiliary signal in (<ref>). This auxiliary signal changes in the presence of saturation and thus alters the control input through the adaptive update laws. The effect of saturation may also be incorporated in the reference model as done in the μ-mod architecture <cit.>. A similar error model can be formulated for the μ-mod architecture, thus the MIMO output feedback magnitude and rate saturation architecture presented this paper can be extended to the μ-mod architecture.




§ FLIGHT CONTROL NUMERICAL SIMULATIONS



This section presents numerical simulation results for an open loop stable nonlinear F-16 vehicle model for single-input-single-output (SISO) longitudinal dynamics and multiple-input-multiple-output (MIMO) lateral-directional dynamics, as well a numerical simulation of a nonlinear hypersonic vehicle model for open loop unstable longitudinal dynamics. For the F-16 vehicle, the proposed adaptive controller maintains stability in the presence of input saturation and demonstrates increased performance compared to a non-adaptive controller. The benefit of the proposed adaptive controller is even more apparent for the open loop unstable hypersonic vehicle dynamics, where the proposed adaptive controller provides for command tracking, whereas non-adaptive and standard adaptive responses are unacceptable.



 §.§ F-16 Numerical Simulation


This section applies the adaptive controller with limiter presented in this paper to two nonlinear F-16 simulations adapted from <cit.>. The aircraft simulations feature nonlinear equations of motion with aerodynamic forces and moments calculated using aerodynamic parameters scheduled from aerodynamic look-up tables based on the flight condition. A trim point for this nonlinear vehicle model was obtained at a straight and level flying condition at a velocity of 500 ft/s with an altitude of 15,000 ft. The vehicle model was linearized about this trim point in order to obtain linearized dynamics for control design as in (<ref>). The actuator time constants and saturation levels from the simulation documentation <cit.> were employed in this simulation, thus providing for physically consistent input constraints. In the following numerical simulations, the thrust input to the aircraft was held constant (at the trim value). It is common in the aerospace industry to separate the design of the full control system of fighter aircraft into longitudinal and lateral-directional dynamics <cit.>.

The longitudinal dynamics of an aircraft describe the short timescale motion in the pitch plane. The longitudinal variables for control design are given by

    x_p=
    [ α q ]^T
     
    u_p=δ_e
     
    y_p=q
     
    z_p=q

where the longitudinal state is composed of the vehicle's angle of attack α and pitch rate q. Pitch rate is both the measured and regulated variable. An elevator deflection δ_e represents the input to the dynamics. The elevator gain parameter τ and saturation levels are as follows:

    τ=0.0495
     
    u_max=± 25 deg 
    u_r,max=± 60 dps.


The lateral-directional dynamics govern motion in the roll and yaw axes and were simulated to demonstrate the MIMO case. The lateral-directional subsystem is given by

    x_p=
    [ β p r ]^T
     
    u_p=
    [ δ_a δ_r ]^T
     
    y_p=
    [ p r ]^T
     
    z_p=p

where the lateral-directional state is composed of the vehicle's angle of sideslip β, roll rate p, and yaw rate r. Roll rate and yaw rate are both measured while roll rate is the regulated variable. Aileron δ_a and rudder δ_r deflections represents the inputs to the dynamics. The aileron gain parameter τ and saturation levels are as follows:

    τ=0.0495
     
    u_max=± 21.5 deg 
    u_r,max=± 80 dps.

The rudder gain parameter τ and saturation levels are as follows:

    τ=0.0495
     
    u_max=± 30 deg 
    u_r,max=± 120 dps.


Separate controllers were designed for each of the two subsystems, longitudinal and lateral-directional. Integral command tracking (<ref>) was included for both subsystems in addition to the explicit inclusion of the dynamics of each magnitude and rate saturated actuator as in (<ref>), resulting in complete plant models as in (<ref>), (<ref>). The plant models for each subsystem satisfy Assumptions <ref>-<ref>, as required for the adaptive control design in this paper.

The first numerical simulation in Figures <ref> and <ref> show the longitudinal dynamics of the F-16 aircraft in the presence of a 75% decrease in the total pitching moment coefficient C_m. In this simulation, the pitch rate command is stepped every two seconds. The significant amount of uncertainty in C_m, results in large degradation of performance for the non-adaptive baseline-only controller (<ref>), where it can be seen in Figure <ref> that the baseline controller does not settle before the step command changes. This is prohibitive in the design of guidance controllers for aircraft which rely on step responses close to the desired response. This motivates the use of an adaptive control input in order to recover the desired performance, even in the presence of uncertainty.

To this end, the adaptive control input presented in this paper was designed as in (<ref>). In order to counter the uncertainty and recover the desired performance, the adaptive controller is seen to require large control rates, which surpass 200 dps and violate the physical actuator limits as seen in Figure <ref>. If the input limits are enforced for the adaptive control input as in <cit.>, which does not include the modified error model (<ref>) and adaptive update law (<ref>), the closed loop system can be seen to be unstable. This is the motivation for the explicit accounting for saturation nonlinearities in the design of the adaptive controller presented in this paper. Figure <ref> shows a stable response for the proposed adaptive controller using the modified error model (<ref>) and saturation update law (<ref>), as described in this paper. In particular, the adaptive controller presented in this paper has a significantly improved pitch rate response as compared to the baseline controller. This improved response is due to the faster response of the control input as seen by the input rate profile in Figure <ref>. The adaptive controller proposed in this paper allows for significant performance recovery in the presence of the input limits. Figure <ref> shows the integrated absolute value of the difference in pitch rate responses between the ideal not-limited adaptive control response and each of the baseline response and limited adaptive control response. In this figure the adaptive controller can be seen to provide for approximately half the integrated error. It should be noted that some accumulated error is expected given the constraints on the control input. Stability is maintained in the adaptive controller presented in this paper by explicitly taking into account extra disturbance due to rate saturation as seen in Figure <ref>.





In order to demonstrate the effectiveness of the proposed adaptive controller for a MIMO system, with both magnitude and rate limits encountered, similar simulations were carried out for the lateral-directional subsystem of the numerical nonlinear F-16 simulation in Figures <ref> and <ref>. In this simulation the aircraft experiences a 50% reduction in rolling moment coefficient C_ℓ. The uncertainty results in a degradation in the performance of the roll rate response for the non-adaptive baseline controller as compared to the ideal adaptive controller for the plant model without input limits, as seen in Figure <ref>. The adaptive controller designed without the error model and update laws modified for the presence of input rate and magnitude limits can be seen to have a severely degraded response. In comparison, the proposed adaptive controller presented in this paper more closely recovers the ideal unconstrained response (up to the performance limits due to the input magnitude and rate constraints). In this simulation the effects of rate limits is even more apparent, where in Figure <ref> it can be more clearly noticed that the limited adaptive control aileron angle response evolves linearly when an aileron rate limit is encountered. A comparison of the responses between the limited baseline and proposed adaptive controllers, along with a plot of the saturation disturbance Δ u_2 for the proposed adaptive controller is shown in Figure <ref>.



 §.§ Hypersonic Vehicle Numerical Simulation



This section provides numerical simulation results for an open loop unstable nonlinear hypersonic vehicle model from <cit.> in order to further demonstrate the need for an adaptive controller capable of handling both rate limits and vehicle uncertainties. This nonlinear vehicle model also has aerodynamic forces and moments calculated from aerodynamic parameters based on flight condition. In order to obtain linearized dynamics for control design as in (<ref>), a trim point was obtained at a straight and level flying condition at Mach 6 at an altitude of 80,000 ft. The control design for the hypersonic vehicle was separated into velocity, longitudinal and lateral-directional dynamics. In order to demonstrate uncertainty in the longitudinal dynamics, separate controllers were designed for the velocity and lateral-directional subsystems to maintain the trim condition in their respective subsystems. As in <cit.>, the longitudinal subsystem is restated as

    x_p=
    [ α q ]^T
     
    u_p=δ_e
     
    y_p=q
     
    z_p=q

where the longitudinal state is composed of the vehicle's angle of attack α and pitch rate q. Pitch rate is both the measured and regulated variable. An elevator deflection δ_e represents the input to the dynamics. The elevator gain parameter τ and saturation levels are as follows:

    τ=0.02
     
    u_max=± 30 deg 
    u_r,max=± 100 dps.






In contrast to the F-16 simulation in the previous subsection, the longitudinal dynamics of this hypersonic vehicle numerical simulation are open loop unstable. Once again, integral command tracking (<ref>) was included in addition to the explicit inclusion of the dynamics of the magnitude and rate saturated actuator as in (<ref>), resulting in complete plant models as in (<ref>), (<ref>) which satisfy Assumptions <ref>-<ref>, as required.

Figures <ref> and <ref> show the response of the hypersonic vehicle in the presence of a 50% decrease in the control effectiveness of all control surfaces and a shift in the center of gravity 0.1 ft rearward, with the pitch rate command stepped every two seconds. The uncertainty present in this simulation significantly degrades the ability of the non-adaptive baseline controller in (<ref>) to track the desired pitch rate command. The angle of attack response of the baseline controller as seen in Figure <ref> additionally becomes unacceptably large for a hypersonic vehicle. This is another example of the motivation to employ an adaptive controller which is capable of countering the effects of uncertainty and restore the desired system response. In the absence of input limits, the adaptive control input, as in <cit.>, and designed in (<ref>) can be seen to result satisfactory command tracking in Figure <ref> by the use of large control rates initially as seen in Figure <ref>. However, in the presence of input magnitude and rate limits, this controller can be seen to result in instability. The proposed adaptive controller which includes the modified error model (<ref>), adaptive update law (<ref>) and explicitly accounts for the presence of input magnitude and rate limits can be seen to result in a stable response which recovers the desired command tracking performance. Figure <ref> shows the integrated absolute value of the error between the ideal unconstrained adaptive control response and each of the baseline response and the proposed adaptive control response. The baseline response can be seen to perform significantly worse as compared to the proposed adaptive control response which explicitly accounts for input limits. The open loop unstable nature of the longitudinal dynamics of the hypersonic vehicle can be seen to result in a greater difference in the performance of the adaptive input limiter controller as compared to the baseline controller.







§ CONCLUSION

This paper presented the first MIMO adaptive architecture for controlling a plant model in the presence of parametric uncertainties in the presence of hard limits on both input magnitude and rate. The plant input is limited in both magnitude and rate through the use of a filter with hard saturation nonlinearities placed in the control path. Augmented adaptive update laws which explicitly account for the presence of the input limits were derived to stabilize the system while in saturation. The proof of stability of signals in the closed-loop system contained two parts: boundedness of the adaptive parameters and boundedness of the state. The main stability result greatly extends the current state-of-the-art in MIMO adaptive control in the presence of rate limits, as it provides a solution in output feedback form with hard saturation nonlinearities employed. The proposed controller was applied to a nonlinear numerical simulation of an open loop stable F-16 aircraft as well as a nonlinear numerical simulation of an open loop unstable hypersonic vehicle. The MIMO simulations demonstrate successful input limiting, stability, and tracking performance of the adaptive controller presented in this paper.





















§ CONSTANTS AND TIME-VARYING SCALARS

Throughout this paper the norm symbol ||X|| denotes the 2-norm of element X.


 §.§ Constant Scalars

Define P_B such that: ||χ^T(t)P_clB_Ω||≤ P_B||χ(t)||

Define P_C such that: ||χ^T(t)P_clC_Δu̅_2||≤ P_C||χ(t)||

Define P_Z such that: ||χ^T(t)P_clB_z||≤ P_Z||χ(t)||

    γ_max=max[eig(Γ_Ω),eig(Γ_Ω_Δ)],  λ_min=min(eig(Λ^*))


    u_min=min_i(u_max,i),   u_r,min=min_i(u_r,max,i)


    p_min=min(eig(P_cl)),   p_max=max(eig(P_cl))


    ρ=√(p_max/p_min),   q_0=min(eig(Q_cl))


    Ω̃_max=sup_t||Ω̃(t)||


    ||B_ξ,Ω||=(||Ω^*||+Ω̃_max)||B_ξ||


    ||K_ξ,Ω||=(||Ω^*||+Ω̃_max)||K_ξ||


    ||   K_u_p||=
       ||C_u_p||(||B_2Λ^*a_1^1||(Ω̃_max||C_ξ̅||+(||Ω^*||+1))+2)


    κ_1=2P_C||S_C||Ω̃_max||Γ_C||


    κ_2=|-q_0+2Ω̃_maxP_B||C_ξ̅||+2P_C1/a_1^1||K_ξ,Ω|||


    κ_3=α u_min-(2P_Z+2P_C1/a_1^1||B_ξ,Ω||)z_cmd,max


    κ_4=Ω̃_max^2a_1^1P_B||S_C||·||Γ_C||·||C_ξ̅||/||K_ξ,Ω||


    κ_5=q_0-3Ω̃_maxP_B||C_ξ̅||


    κ_6=(2P_Z+Ω̃_maxP_B||B_ξ,Ω||·||C_ξ̅||/||K_ξ,Ω||)z_cmd,max


    κ_7=|-q_0+2Ω̃_maxP_B||C_ξ̅||+2P_C1/a_1^1(||K_ξ,Ω||+||K_u_p||)|


    κ_8=β u_r,min-(2P_Z+2P_C1/a_1^1||B_ξ,Ω||)z_cmd,max


    κ_9=Ω̃_max^2a_1^1P_B||S_C||·||Γ_C||·||C_ξ̅||/||K_ξ,Ω||+||K_u_p||


    κ_10=(2P_Z+Ω̃_maxP_B||B_ξ,Ω||·||C_ξ̅||/||K_ξ,Ω||+||K_u_p||)z_cmd,max


    κ_11=2P_Zz_cmd,max


    κ_12=q_0-2Ω̃_maxP_B||C_ξ̅||


    χ_min=max   (κ_5-√(κ_5^2-4κ_4κ_6)/2κ_4,.
       .κ_5-√(κ_5^2-4κ_9κ_10)/2κ_9,κ_11/κ_12)


    χ_max=min   (√(κ_2^2+4κ_1κ_3)-κ_2/2κ_1,κ_5+√(κ_5^2-4κ_4κ_6)/2κ_4,.
       .√(κ_7^2+4κ_1κ_8)-κ_7/2κ_1,κ_5+√(κ_5^2-4κ_9κ_10)/2κ_9)


    Ω̅_max=q_0/3P_B||C_ξ̅||

Note that Ω̃_max<√(V(0)γ_max/λ_min). Additionally, the ratio q_0/p_max is maximized with a choice of Q_cl=I as seen in reference <cit.>. A bound on Ω̃ may be enforced using the projection operator as is common in adaptive control <cit.>.




 §.§ Constant Matrices


    A_cl=
    [ -a_1^0/a_1^1I      1/a_1^1I             0             0;             0        A-B_2K            LC             0;             0             0         A^*_L      B_2^1Λ^*;             0             0             0 -a_1^0/a_1^1I ]


    B_Ω^T=
    [            0            0 (B_2^1Λ^*)^T            0 ]


    B_Z^T=
    [     0 B_z^T     0     0 ],   
    B_ξ^T=
    [          0 a_1^1B_z^T          0 ]


    C_Δu̅_2^T=
    [ 0 0 0 I ],  
    C_u_p=
    [     0 Λ^*-1     0 ]


    C_ξ̅=
    [ K 0 0 0; 0 I 0 0; I 0 0 0 ]


    K_ξ=
    [                      0                      K                      0                      0;                      0 (a_1^0I+a_1^1(A-B_2K))                a_1^1LC                      0;                      0                      I                      0                      0 ]


    S_C=S_2C,  Γ_C=C_ξ̅^TΓ_ΩC_ξ̅






 §.§ Time-Varying Scalars

For ease of exposition, the following time varying scalars which represent the ratio of the saturated input magnitude and rate respectively to the unsaturated input magnitude and rate:

    𝕌=||E_s(u(t),u_max)||/||u(t)||,  𝕌_r=||E_s(u_r(t),u_r,max)||/||u_r(t)||

The magnitude of the computed control input (<ref>) may be expressed as follows:

    2
    ||u(t)||   =       ||(a^1_1s+a^0_1)Ω^T(t)ξ̅(t)||
       =       ||Ω^T(t)(a_1^1s+a_1^0)ξ̅(t)+a_1^1Ω̇^T(t)ξ̅(t)||
       =       ||-Ω^T(t)K_ξχ(t)-Ω^T(t)B_ξz_cmd(t)
              -a_1^1S_2Ce_mu(t)χ^T(t)Γ_Cχ(t)||
       ≤       ||K_ξ,Ω||·||χ(t)||+||B_ξ,Ω||z_cmd,max
             +a_1^1||S_C||Ω̃_max||Γ_C||·||χ(t)||^2

The magnitude of the plant input may be expressed as follows, through applying equation (<ref>):

    2
    ||u_p(t)||   =       ||C_u_p(e_x(t)+x_m(t))||
       =       ||C_u_p(B_2Λ^*a_1^1[Ψ̅^*T_2x̅_m(t)+Ω^T(t)ξ̅(t).
             .-Λ^*-1u̅_bl(t)+Δu̅_2(t)]+e_mx(t)+x_m(t))||
       ≤       ||K_u_p||·||χ(t)||

Additionally, from equation (<ref>) the following inequality holds:

    ||u_r(t)||≤1/τ(||u(t)||+||u_p(t)||)




§ PROOF OF THEOREM <REF>, STEP 2



 §.§ Δ u=0,Δ u_r=0

In this case (<ref>) can be simplified as follows:

    χ̇(t)=A_clχ(t)-B_ΩΩ̃^T(t)C_ξ̅χ(t)+B_Zz_cmd(t)

This leads to the following time derivative of the Lyapunov function in equation (<ref>):

    2Ẇ(χ(t))   =       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clB_Zz_cmd(t)

The right hand side may be bounded as follows:

    2Ẇ(χ(t))   ≤       -(q_0-2Ω̃_maxP_B||C_ξ̅||)||χ(t)||^2
             +(2P_Zz_cmd,max)||χ(t)||

Using the definition of Ω̃_max and condition (ii) of Theorem <ref>, implies the following:

    Ẇ(χ(t))<0,   ||χ(t)||>κ_11/κ_12

Choosing χ_min from equation (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 1











 §.§ Δ u≠0,Δ u_r=0

In this case (<ref>) can be rewritten as:

    2χ̇(t)   =       A_clχ(t)-B_ΩΩ̃^T(t)C_ξ̅χ(t)
             +C_Δu̅_21/a_1^1(𝕌-1)u(t)+B_Zz_cmd(t)

This leads to the following time derivative of the Lyapunov function in equation (<ref>):

    2Ẇ(χ(t))   =       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1(𝕌-1)u(t)
             +2χ^T(t)P_clB_Zz_cmd(t)

Three sub-cases are considered in this section.



  §.§.§ Sub-case i


    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)<-α u_min||χ(t)||

The time derivative of the Lyapunov function may be bounded as follows using the condition of this sub-case and (<ref>):

    Ẇ(χ(t))≤(2P_C||S_C||Ω̃_max||Γ_C||)||χ(t)||^3
       +(-q_0+2Ω̃_maxP_B||C_ξ̅||+2P_C1/a_1^1||K_ξ,Ω||)||χ(t)||^2
       -(α u_min-(2P_Z+2P_C1/a_1^1||B_ξ,Ω||)z_cmd,max)||χ(t)||

Thus:

    Ẇ(χ(t))<0,   ||χ(t)||<√(κ_2^2+4κ_1κ_3)-κ_2/2κ_1

Choosing χ_max from equation (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 2, sub-case (i)











  §.§.§ Sub-case ii


    0>2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)>-α u_min||χ(t)||

The condition of this sub-case implies the following:

    0<2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)+α u_min||χ(t)||

Apply u_min/||E_s(u(t),u_max)||≤ 1:

    0<2χ^T(t)P_clC_Δu̅_21/a_1^1u(t)+α||u(t)||·||χ(t)||

Add terms to create Ẇ(χ(t)) on the left hand side:

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1(𝕌-1)u(t)
             +2χ^T(t)P_clB_Zz_cmd(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1u(t)
             +α||u(t)||·||χ(t)||

Apply case inequality:

    0>2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)

Thus:

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +α||u(t)||·||χ(t)||+2χ^T(t)P_clB_Zz_cmd(t)

Apply (<ref>):

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +α ||K_ξ,Ω||·||χ(t)||^2
             +α ||B_ξ,Ω||z_cmd,max||χ(t)||
             +α a_1^1||S_C||Ω̃_max||Γ_C||·||χ(t)||^3
             +2χ^T(t)P_clB_Zz_cmd(t)

Maximize the right hand side:

    2Ẇ(χ(t))   <       (α a_1^1||S_C||Ω̃_max||Γ_C||)||χ(t)||^3
             -(q_0-2P_BΩ̃_max||C_ξ̅||-α ||K_ξ,Ω||)||χ(t)||^2
             +((2P_Z+α ||B_ξ,Ω||)z_cmd,max)||χ(t)||

Apply α=P_BΩ̃_max||C_ξ̅||/||K_ξ,Ω||:

    2Ẇ(χ(t))   <       (Ω̃_max^2a_1^1P_B||S_C||·||Γ_C||·||C_ξ̅||/||K_ξ,Ω||)||χ(t)||^3
             -(q_0-3Ω̃_maxP_B||C_ξ̅||)||χ(t)||^2
             +((2P_Z+Ω̃_maxP_B||B_ξ,Ω||·||C_ξ̅||/||K_ξ,Ω||).
             .× z_cmd,max)||χ(t)||

Using the definition of Ω̃_max and condition (ii) of Theorem <ref>, implies the following:

    Ẇ(χ(t))<0,  for
       κ_5-√(κ_5^2-4κ_4κ_6)/2κ_4<||χ(t)||<κ_5+√(κ_5^2-4κ_4κ_6)/2κ_4

Choosing χ_min,χ_max from equations (<ref>) and (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 2, sub-case (ii)





  §.§.§ Sub-case iii


    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)>0

The condition of this sub-case implies the following inequality:

    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌u(t)<2χ^T(t)P_clC_Δu̅_21/a_1^1u(t)

Add terms to create Ẇ(χ(t)) on the left hand side:

    2Ẇ(χ(t))   ≤       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clB_Zz_cmd(t)

The right hand side may be maximized as:

    2Ẇ(χ(t))   ≤       -(q_0-2Ω̃_maxP_B||C_ξ̅||)||χ(t)||^2
             +(2P_Zz_cmd,max)||χ(t)||

Using the definition of Ω̃_max and condition (ii) of Theorem <ref>, implies the following:

    Ẇ(χ(t))<0,   ||χ(t)||>κ_11/κ_12

Choosing χ_min from equation (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 2, sub-case (iii)






 §.§ Δ u=0,Δ u_r≠0

In this case (<ref>) can be rewritten as:

    2χ̇(t)   =       A_clχ(t)-B_ΩΩ̃^T(t)C_ξ̅χ(t)
             +C_Δu̅_21/a_1^1(𝕌_r-1)τ u_r(t)+B_Zz_cmd(t)

The following is the time derivative of the Lyapunov function in equation (<ref>):

    2Ẇ(χ(t))   =       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1(𝕌_r-1)τ u_r(t)
             +2χ^T(t)P_clB_Zz_cmd(t)

Three sub-cases are considered in this section.















  §.§.§ Sub-case i


    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)<-β u_r,min||χ(t)||

The time derivative of the Lyapunov function may be bounded as follows using the condition of this sub-case and (<ref>):

    Ẇ(χ(t))≤(2P_C||S_C||Ω̃_max||Γ_C||)||χ(t)||^3
       +(-q_0+2Ω̃_maxP_B||C_ξ̅||.
       .+2P_C1/a_1^1(||K_ξ,Ω||+||K_u_p||))||χ(t)||^2
       -(β u_r,min-(2P_Z+2P_C1/a_1^1||B_ξ,Ω||)z_cmd,max)||χ(t)||

Thus:

    Ẇ(χ(t))<0,   ||χ(t)||<√(κ_7^2+4κ_1κ_8)-κ_7/2κ_1

Choosing χ_max from equation (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 3, sub-case (i)








  §.§.§ Sub-case ii


    0>2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)>-β u_r,min||χ(t)||

The condition of this sub-case implies the following:

    0<2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)+β u_r,min||χ(t)||

Apply u_r,min/||E_s(u_r(t),u_r,max)||≤ 1:

    0<2χ^T(t)P_clC_Δu̅_21/a_1^1τ u_r(t)+β||u_r(t)||·||χ(t)||

Add terms to create Ẇ(χ(t)) on the left hand side:

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1(𝕌_r-1)τ u_r(t)
             +2χ^T(t)P_clB_Zz_cmd(t)
             +2χ^T(t)P_clC_Δu̅_21/a_1^1τ u_r(t)
             +β||u_r(t)||·||χ(t)||

Apply case inequality:

    0>2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)

Thus:

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +β||u_r(t)||·||χ(t)||+2χ^T(t)P_clB_Zz_cmd(t)

Apply (<ref>), (<ref>), and (<ref>):

    2Ẇ(χ(t))   <       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +β/τ(||K_ξ,Ω||+||K_u_p||)||χ(t)||^2
             +β/τ ||B_ξ,Ω||z_cmd,max||χ(t)||
             +β/τ a_1^1||S_C||Ω̃_max||Γ_C||·||χ(t)||^3
             +2χ^T(t)P_clB_Zz_cmd(t)

Maximize the right hand side:

    2Ẇ(χ(t))   <       (β/τ a_1^1||S_C||Ω̃_max||Γ_C||)||χ(t)||^3
             +(-q_0+2P_BΩ̃_max||C_ξ̅||.
             .+β/τ(||K_ξ,Ω||+||K_u_p||))||χ(t)||^2
             +((2P_Z+β/τ||B_ξ,Ω||)z_cmd,max)||χ(t)||

Apply β=τ P_BΩ̃_max||C_ξ̅||/||K_ξ,Ω||+||K_u_p||:

    2Ẇ(χ(t))   <       (Ω̃_max^2a_1^1P_B||S_C||·||Γ_C||·||C_ξ̅||/||K_ξ,Ω||+||K_u_p||)||χ(t)||^3
             -(q_0-3Ω̃_maxP_B||C_ξ̅||)||χ(t)||^2
             +((2P_Z+Ω̃_maxP_B||B_ξ,Ω||·||C_ξ̅||/||K_ξ,Ω||+||K_u_p||).
             .× z_cmd,max)||χ(t)||

Using the definition of Ω̃_max and condition (ii) of Theorem <ref>, implies the following:

    Ẇ(χ(t))<0,  for
       κ_5-√(κ_5^2-4κ_9κ_10)/2κ_9<||χ(t)||<κ_5+√(κ_5^2-4κ_9κ_10)/2κ_9

Choosing χ_min,χ_max from equations (<ref>) and (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 3, sub-case (ii)










  §.§.§ Sub-case iii


    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)>0

The condition of this sub-case implies the following inequality:

    2χ^T(t)P_clC_Δu̅_21/a_1^1𝕌_rτ u_r(t)<2χ^T(t)P_clC_Δu̅_21/a_1^1τ u_r(t)

Add terms to create Ẇ(χ(t)) on the left hand side:

    2Ẇ(χ(t))   ≤       -χ^T(t)Q_clχ(t)-2χ^T(t)P_clB_ΩΩ̃^T(t)C_ξ̅χ(t)
             +2χ^T(t)P_clB_Zz_cmd(t)

The right hand side may be maximized as:

    2Ẇ(χ(t))   ≤       -(q_0-2Ω̃_maxP_B||C_ξ̅||)||χ(t)||^2
             +(2P_Zz_cmd,max)||χ(t)||

Using the definition of Ω̃_max and condition (ii) of Theorem <ref>, implies the following:

    Ẇ(χ(t))<0,   ||χ(t)||>κ_11/κ_12

Choosing χ_min from equation (<ref>) implies the following:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜 in case 3, sub-case (iii)

From equations (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) it can be concluded that:

    Ẇ(χ(t))<0,  ∀χ(t)∈𝒜











§ ACKNOWLEDGMENT



This work was supported by the Air Force Research Laboratory, Collaborative Research and Development for Innovative Aerospace Leadership (CRDInAL), Thrust 3 - Control Automation and Mechanization grant FA 8650-16-C-2642 and the Boeing Strategic University Initiative. The authors acknowledge Dr. Daniel Wiese for providing useful hypersonic vehicle numerical simulation scripts.
























IEEEtran























































