# The Quantum Mechanical Problem of a Particle on a Ring with Delta Well

Raphael J.F. Berger  Division of Chemistry and Physics of Materials, Paris-
Lodron University Salzburg, Salzburg, A-5020, Austria
<EMAIL_ADDRESS>
(Date: August 28, 2024)

###### Abstract.

The problem of a spin-free electron with mass $m$, charge $e$ confined onto a
ring of radius $R_{0}$ and with an attractive Dirac delta potentential with
scaling factor (depth) $\kappa$ in non-relativistic theory has closed form
analytical solutions. The single bound state function is of the form of a
hyperbolic cosine that however contains a paramter $d>0$ which is the single
positive real solution of the transcentdental equation $\coth(d)=\lambda d$
for non zero real $\lambda=\frac{2}{\pi\kappa}$. The energy eigenvalue of the
bound state
$\varepsilon=-\frac{d^{2}}{2\pi^{2}}\approx\frac{qemR_{0}}{2\hbar^{2}}$. In
addition a discrete inifinty of unbounded solutions exists, formally they are
obtained from the terms for the bound solution by substituting $d\to id$
yielding $\cot(d)=\lambda d$ as characteristic equation with the correspondig
set of solutions $d_{k},k\in\mathbb{N}$, the respective state functions obtain
via $\cosh(x)\overset{x\to ix}{\longrightarrow}\cos(x)$ of the form of cosine
functions.

RJFB acknowledges funding from DFG (German Research Foundation) within the
Priority Program SPP1807 “Control of LD in molecular chemistry”, grant
BE4632/2-2, project no. 271386299)

##  1\. Introduction

There is only roughly a dozen of quantum mechanical (QM) systems with an
analytical solution.[Wikipedia contributors(2022)] QM problems with analytical
solutions are not only of great of didactical use for showin how these
problems can be solved but often serve also as physical toy model systems for
otherwise unsovable problems which however share the principal
characteristics.

One example herefore is the particle in the Dirac delta potential[Wikipedia
contributors(2022)] which not ony serves with its bound state as an one-
dimensional analog of the hydrogen atom but can also be interpreted as the
simplest model for electron scattering when one regards the unbound states and
allows for simple calculation of reflection and transmissoin rates on step-
and related potentials, which is for example of relevance for the theory of
scanning tunnel microscopy.

In the context of our research on symmetry breaking in rotationally invariant
systems[Berger and Viel(2020)] we came across the analogous problem where the
particle is but confined to an atomic scale ring. As it turned out that this
QM model system has not yet been described in the literature (to the best of
our knowledge) we report on our results in the following.

##  2\. Solution

###  2.1. Schrödinger equation

A spin-free electron (i.e. a particle with charge $-e$ and mass $m_{e}$) on in
a ring shaped space with radius $R_{0}$ and a $\delta$ function well of
amplitude $-qe$ corresponding to an attractive potential with an integrated
total charge of $+qe$ is regarded in non-relativistic quantum mechanic theory.

As is well known from textbooks the Hamiltonian for the particle on a ring of
radius $R_{0}$ in 2D polar coordinates
($r\in\mathbb{R}^{+},\vartheta\in[0,2\pi)$) in SI units is given by

(1) $\displaystyle\hat{H^{\prime}}$
$\displaystyle=-\frac{\hbar^{2}}{2mR_{0}^{2}}\frac{\partial^{2}}{\partial\vartheta^{2}}$

The ring shall contain an attractive potential $\hat{V}$ with respect to the
electron in the form of a Dirac-$\delta$ function111The dimension of the
argument of the $\delta$ function has to be coosen such as $\int\delta\;dx=1$
is dimensionless and the potential shall integrate over the whole space
($R_{0}\times[0,2\pi)$) to the product of electron and potential charge of
$-eq<0$:

(2) $\displaystyle\hat{V}$
$\displaystyle=-qe\delta((\vartheta-\vartheta_{0})l),$

were $\vartheta$ is given in units of $\mathrm{[}rad]$, thus formally
$l=1\mathrm{rad}^{-1}=1\frac{1m}{1m}=1$ and hence can be dropped in the
following. Without loss of generality we will set $\vartheta_{0}=0$ such that
the Schrödinger equation for the problem becomes

$\displaystyle\hat{H}\psi(\vartheta)$ $\displaystyle=E\psi(\vartheta)$
$\displaystyle(\hat{H^{\prime}}+\hat{V})\psi(\vartheta)$
$\displaystyle=E\psi(\vartheta)$ (3)
$\displaystyle\left(-\frac{\hbar^{2}}{2mR_{0}^{2}}\frac{\partial^{2}}{\partial\vartheta^{2}}-qe\delta(\vartheta)\right)\psi(\vartheta)$
$\displaystyle=E\psi(\vartheta)$
$\displaystyle\left(\frac{\partial^{2}}{\partial\vartheta^{2}}+\frac{2qemR_{0}^{2}}{\hbar^{2}}\delta(\vartheta)+E\frac{2mR_{0}^{2}}{\hbar^{2}}\right)\psi(\vartheta)$
$\displaystyle=0$

For simplicity we combine the constants in

(4) $\displaystyle\kappa=\frac{2qemR_{0}^{2}}{\hbar^{2}}$

and

(5) $\displaystyle\epsilon=\frac{2mR_{0}^{2}}{\hbar^{2}}E$

such that we obtain

(6) $\displaystyle\psi^{\prime\prime}+(\epsilon+\kappa\delta)\psi$
$\displaystyle=0$

###  2.2. Bound state ($E<0$)

The assumption of $(E-V)<0$ or (ignoring the singularity at the origin) $E<0$,
or $\epsilon<0$, respectively thus leads to the bound state solutions. For the
symmetry of the problem and the form of the differential equation (6) we chose
the Ansatz

(7)
$\displaystyle\psi(\vartheta;d)=N^{\prime}(e^{-d\vartheta}+e^{d(\vartheta-2\pi)})=N\cosh(d(x-\pi))$

Yielding the normalization constant

(8) $\displaystyle N=\sqrt{\frac{\sinh(2\pi d)}{2d}+\pi}^{-1}$

To determine the exponent $d$ one in principle has to insert (7) into (6) and
attempt to match $d$ to the boundary conditions. One boundary condition, the
symmtry of the system, was already accounted for choosing the same exponent
$d$ for both functions but with different sign. Since a $\delta$-function is
appearing one has to chose an appropriate strategy to evaluate the result of
combining (7) and (6). The strategy is to integrate the Schrödinger equation
in an $\epsilon$ ball around the origin of the $\delta$-function and to
perform the limit of $\epsilon\to 0$, this results in

$\displaystyle\lim_{\epsilon\to
0}\int_{0-\epsilon}^{0+\epsilon}\left(\frac{\partial^{2}}{\partial\vartheta^{2}}\psi(\vartheta)\right)+\kappa\delta(\vartheta)\psi(\vartheta)-E^{\prime}\psi(\vartheta)\;d\vartheta$
$\displaystyle=\lim_{\epsilon\to 0}\int_{0-\epsilon}^{0+\epsilon}0d\vartheta$
(9) $\displaystyle\psi^{\prime}(0^{+})-\psi^{\prime}(0^{-})+\kappa\psi(0)$
$\displaystyle=0$

Inserting (7) in (9) yields

$\displaystyle-2d+2de^{-2\pi d}+\kappa(1+e^{-2\pi d})$ $\displaystyle=0$ (10)
$\displaystyle\coth(\pi d)$ $\displaystyle=\frac{2}{\kappa}d$

using $d^{\prime}=\pi d$ and $\lambda=\frac{2}{\pi\kappa}$ we obtain

(11) $\displaystyle\coth(d^{\prime})$ $\displaystyle=\lambda d^{\prime}$

here for $\lambda>0$ and real $d^{\prime}$ exactly one pair of solutions
$d^{\prime}_{+}=-d^{\prime}_{-}$ exists and yielding the same function
$\psi(\vartheta;d)$ due to the axial symmetry of $\cosh$, hence we can drop
the $\pm$ indices in the following and only regard the positive solution
$d^{\prime}$, which is it the same time the only solution to (6).

Equation (11) has no symbolic closed form solution, but

(12) $\coth(d^{\prime})\approx 1$

holds with accuracy increasing in $d^{\prime}$ (see Fig. 1).

Figure 1. Realtive error in % obtained from using (12) in dependence on
$d^{\prime}$

Using (12) the solution of (10) can be approximated as

$\displaystyle d\approx d^{0}=\frac{\kappa}{2}=\frac{qemR_{0}^{2}}{\hbar^{2}}$

with deviations decreasing with charge and ring size. In the Table 1 some
exemplary values for $R_{0}$ set to 1 (all atomic units) are shown. At
$R_{0}=1$ bohr and $q=1e$ we have $\epsilon=-d^{2}=-\frac{d^{\prime
2}}{\pi^{2}}$ and $E=-\frac{d^{\prime
2}\hbar^{2}}{2m\pi^{2}R_{0}^{2}}\approx-0.0145$ Hartree (for comparision in
the approximation (12) it yields
$E^{0}=-\frac{{d^{0}}^{2}}{\pi^{2}}\approx-0.0127$ Hartree (further values are
listed in Table 1). In addition we note that with increasing radius $R_{0}$
the approximative solutions will be of increasing quality.

$q$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
${d^{0}}^{\prime}$ | 0.5 | 1.0 | 1.5 | 2.0 | 2.5
$d^{\prime}$ | 0.53575 | 1.00366 | 1.50024 | 2.00001 | 2.5
${E^{0}}^{\prime}$ | -0.01267 | -0.05066 | -0.11399 | -0.202642 | -0.31663
$E^{0}$ | -0.01454 | -0.05103 | -0.11402 | -0.202642 | -0.31663
Table 1. Numerically exact ($d^{\prime}$) and approximated solutions
(${d^{0}}^{\prime}$) of (10) at different total charges $q$ of the Dirac delta
potential and corresponding approximate (${E^{0}}^{\prime}$) and numerically
exact energies (${E^{0}}^{\prime}$) at $1$ Bohr radius.

A graphical display of the wave function of the bound state is shown in Fig.
2.

Figure 2. Stereographic projection of $\psi(\vartheta)$ for the bound state of
the partice in the ring with $\delta$ well.

###  2.3. Unbound states ($E>0$)

In the spirit of (7) the unbound states can be obtained from the Ansatz

(13)
$\displaystyle\psi(\vartheta;di)=C(e^{-di\vartheta}+e^{di(\vartheta-2\pi)})$

Yielding in analogy to (10)

$\displaystyle-2id+2ide^{-i2\pi d}+\kappa(1+e^{-i2\pi d})$ $\displaystyle=0$
(14) $\displaystyle\cot(\pi d)$ $\displaystyle=\frac{2}{\kappa}d$

In contrast to (10) this yields for finite positive $\kappa$ an infinite set
of solutions. As in (10) the system is strictly not analytically solvable in
terms of a finite closed expression. However, for sufficiently large $d$ or
large values of $\kappa$ the solutions are approaching

(15) $d_{n_{+/-}}\approx d^{(0)}_{n_{+/-}}=\pm n\frac{\kappa}{2}$

The first 5 unbound solutions for $d$ and the corresponding energies of the
sytem with charge $+1e$ and $R=1$ bohr are

$n$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$d^{(0)}_{n_{+}}$ | 0.5 | 1.0 | 2.0 | 3.0 | 4.0
$d_{n_{+}}$ | 0.34278 | 1.15979 | 2.09395 | 3.06518 | 4.04963
$E_{n}$ | 0.05875 | 0.67256 | 2.19231 | 4.69766 | 8.19976

Since $\psi(\vartheta;di)$ are not purely real functions, we first decompose
them into real and imaginary part

(16) $\displaystyle\Re{[\psi(\vartheta;di)]}$
$\displaystyle=C[\cos(d\vartheta)+\cos(d(2\pi-\vartheta))]$ (17)
$\displaystyle=\left(\frac{\sin(2\pi
d)}{2d}+\pi\right)^{-\frac{1}{2}}\cos[d(\pi-\vartheta^{\prime})]$ (18)
$\displaystyle\Im{[\psi(\vartheta;di)]}$
$\displaystyle=C[\sin(d\vartheta)+\sin(d(2\pi-\vartheta))]$

and we note that (18) in general is dicontinuous at the $\delta$ well, thus
must be rejected. (16) can be rewritten as a single cosine function
originating at the position opposing the origin $\vartheta^{\prime}_{o}=\pi$

$\displaystyle\Re{[\psi(\vartheta;di)]}$
$\displaystyle=2C\cos(d\pi)[\cos(d(\pi-\vartheta))]$ (19)
$\displaystyle=C^{\prime}\cos[d(\pi-\vartheta^{\prime})]$

with $-\pi\leq\vartheta^{\prime}<\pi$ and the normalisation constant

(20) $\displaystyle C^{\prime}=\left(\frac{\sin(2\pi
d)}{2d}+\pi\right)^{-\frac{1}{2}}$

In summary this yields the unbounded state functions

(21) $\displaystyle\psi_{n}$ $\displaystyle=\sqrt{\frac{\sin(4\pi
d_{n})}{4d_{n}}+\pi}^{-1}\cos[d_{n}(\pi-\vartheta)]$

for $-\pi\leq\vartheta<\pi$ and $n>0$ with energies

(22) $\displaystyle E_{n}=-\frac{\hbar^{2}}{2mR_{0}^{2}}d_{n}^{2}$

where $d_{n}$ are corresponding to the positive solutions of

(23) $\displaystyle\cot(\pi d_{n})$
$\displaystyle=\frac{2\pi\hbar^{2}}{qemR_{0}}d_{n}.$

which are approximated for large $q$, $R_{0}$ or $n$ by

(24) $d_{n}\approx d^{0}_{n}=n\frac{qR_{0}em}{2\pi\hbar^{2}}$

where we have removed the symmetry equivalent negative solutions, since $\cos$
is an even function, thus droped the sign indices, as compared to (15). Hereby
we note that in comparison to the particle in the ring (without additional
well potential) we have lost the two-fold degeneracy of the higher (non-
ground) states. Which is an obvious consequence of the symmetry breaking due
to the potential.

###  2.4. $\kappa$ dependence

The $\kappa$ dependence of the real solutions of (28) is illustrated in the
graph below. The approximative solutions (24) we have used are based on the
asymptotic approximation of the $\tan$ branches to $x=2n+1$ for
$n\in\mathbb{R}$.

Figure 3. Graphical representation of the solutions of the characteristic
equation (23) for the problem with $\kappa=2$

##  3\. Smooth wave function for the bound state

Independently of the physical significance, one can ask for conditions under
which the wave function becomes smooth at the position of the Dirac delta
function potential. Smoothness at $\vartheta=0$ requires $d=1/2$ in the
bounded state (7). Hence we obtain from (10), (4) and using atomic units

(25) $\displaystyle R_{0}=\sqrt{\frac{\tanh{\pi/2}}{2q}}$

yielding, e.g. $R_{0}=0.677183$ Bohr for $q=1$ at $E_{0}=-0.01267$ Hartree.

##  4\. Summary

Due to the relations $\cos(ix)=\cosh(x)$ and $\sin(ix)=\sinh(x)$ we can
combine the bound (7) and the unbound (21) solutions in one expression

(26) $\displaystyle\psi_{n}$ $\displaystyle=\sqrt{\frac{\sin(4\pi
a_{n})}{4a_{n}}+\pi}^{-1}\cos[a_{n}(\pi-\vartheta)]$

with $\kappa=qR_{0}\frac{em}{\pi\hbar^{2}}$, and the corresponding energies

(27) $\displaystyle E_{n}=\frac{\hbar^{2}}{2mR_{0}^{2}}a_{n}^{2},\;\;\forall
n\in\mathbb{N}_{0}.$

and where $a_{0}$ is the (single) purely imaginary solution and $a_{n}$ with
$n>0$ are the purely real solutions of the equation

(28) $\cot(\pi a)=\frac{2}{\kappa}a$

## References

  * [Wikipedia contributors(2022)] Wikipedia contributors, _List of quantum-mechanical systems with analytical solutions — Wikipedia, The Free Encyclopedia_ , https://en.wikipedia.org/w/index.php?title=List_of_quantum-mechanical_systems_with_analytical_solutions&oldid=1121884573, 2022, [Online; accessed 28-November-2022].
  * [Wikipedia contributors(2022)] Wikipedia contributors, _Delta potential — Wikipedia, The Free Encyclopedia_ , https://en.wikipedia.org/w/index.php?title=Delta_potential&oldid=1117753099, 2022, [Online; accessed 28-November-2022].
  * [Berger and Viel(2020)] R. J. F. Berger and A. Viel, _Zeitschrift für Naturforschung B_ , 2020, 75, 327–339.


# A Data-Driven Autopilot for Fixed-Wing Aircraft
Based on Model Predictive Control

Riley J. Richards, Juan A. Paredes, and Dennis S. Bernstein  Riley J.
Richards, Juan A. Paredes, and Dennis S. Bernstein are with the Department of
Aerospace Engineering, University of Michigan, Ann Arbor, MI, USA. {rileyric,
jparedes<EMAIL_ADDRESS>
###### Abstract

Autopilots for fixed-wing aircraft are typically designed based on linearized
aerodynamic models consisting of stability and control derivatives obtained
from wind-tunnel testing. The resulting local controllers are then pieced
together using gain scheduling. For applications in which the aerodynamics are
unmodeled, the present paper proposes an autopilot based on predictive cost
adaptive control (PCAC). As an indirect adaptive control extension of model
predictive control, PCAC uses recursive least squares (RLS) with variable-rate
forgetting for online, closed-loop system identification. At each time step,
RLS-based system identification updates the coefficients of an input-output
model whose order is a hyperparameter specified by the user. For MPC, the
receding-horizon optimization can be performed by either the backward-
propagating Riccati equation or quadratic programming. The present paper
investigates the performance of PCAC for fixed-wing aircraft without the use
of any aerodynamic modeling or offline/prior data collection.

##  I Introduction

A fundamental necessity for autonomous atmospheric flight vehicles is a
reliable autopilot for controlling the attitude and flight path. For a fixed-
wing vehicle, stability and control derivatives are typically determined
through wind-tunnel testing or computational modeling over a range of Mach
number, angle of attack, and sideslip angle. This modeling data is then used
to develop an autopilot based on gain scheduling, feedback linearization, or
dynamic inversion [1, 2]. In practice, however, the aerodynamics of an
aircraft may be too expensive to model with high accuracy or may change due to
atmospheric conditions, such as icing, as well as damage. This possibility
motivates the need to develop an adaptive autopilot that can learn and respond
to unknown, changing conditions [3, 4, 5].

The present paper introduces a novel adaptive autopilot for fixed-wing
aircraft based on model predictive control (MPC). MPC is widely viewed as the
most effective modern control technique, due to its ability to enforce state
and control constraints in both linear and nonlinear systems through receding-
horizon optimization [6, 7, 8, 9].

As its name suggests, MPC requires a model for optimization. In the absence of
a reliable model, data-driven techniques can be used to update the plant model
during operation [10, 11, 12]. In effect, data-driven techniques perform
closed-loop system identification, at least to a level of accuracy that is
sufficient for the feedback controller to achieve desired closed-loop
performance. The interplay between system identification and control is a
fundamental, longstanding problem in control theory [13, 14, 15, 16, 17].

To develop an adaptive autopilot for fixed-wing aircraft, the present paper
focuses on predictive cost adaptive control (PCAC) [18, 19], which is an
indirect adaptive control extension of MPC. For online, closed-loop system
identification, PCAC uses recursive least squares (RLS) with variable-rate
forgetting [20, 21, 22, 23, 24]. At each time step, RLS-based system
identification updates the coefficients of a SISO or MIMO input-output model,
where the model order is a hyperparameter specified by the user. For MPC, the
receding-horizon optimization can be performed by either the backward-
propagating Riccati equation [25, 6] or quadratic programming.

The objective of the present paper is to investigate the performance of PCAC
as a data-driven autopilot for an aircraft with unmodeled kinematics,
dynamics, and aerodynamics. In particular, PCAC is implemented as a cold-start
indirect adaptive controller, where the plant model order is specified as a
hyperparameter, but otherwise no plant model is assumed to be available. The
identified model updated by RLS is linear, and thus it is suitable for
modeling the aircraft dynamics near trim. In practice, an autopilot designed
to operate over a wide range of flight conditions depends on gain scheduling
of multiple linear controllers. The goal of this study is to investigate, via
numerical experiments, the viability and potential performance of PCAC under
conditions of high uncertainty, in effect, no prior modeling information,
without the need for gain scheduling.

For this study, we first consider scenarios involving the longitudinal and
lateral dynamics of a linearized fixed-wing aircraft model provided by Athena
Vortex Lattice (AVL). Next, we consider the longitudinal dynamics of a
nonlinear fixed-wing aircraft model provided in the MATLAB aerospace toolbox.

The contents of the paper are as follows. Section II introduces the predictive
control problem. Section III reviews the PCAC formulation and algorithm.
Section IV applies PCAC to a linear 6DOF aircraft model for longitudinal and
lateral control, and Section V applies PCAC to a nonlinear 3DOF aircraft model
for longitudinal control. Finally, Section VI presents conclusions and future
research directions.

Notation: ${\bf z}\in{\mathbb{C}}$ denotes the Z-transform variable. $x_{(i)}$
denotes the $i$th component of $x\in{\mathbb{R}}^{n}.$ ${\rm sprad}(A)$
denotes the spectral radius of $A\in{\mathbb{R}}^{n\times n}.$ The symmetric
matrix $P\in{\mathbb{R}}^{n\times n}$ is positive semidefinite (resp.,
positive definite) if all of its eigenvalues are nonnegative (resp.,
positive). $\operatorname{vec}X\in{\mathbb{R}}^{nm}$ denotes the vector formed
by stacking the columns of $X\in{\mathbb{R}}^{n\times m}$, and $\otimes$
denotes the Kronecker product. $I_{n}$ is the $n\times n$ identity matrix, and
$0_{n\times m}$ is the $n\times m$ zeros matrix and $\mathds{1}_{n\times m}$
is the $n\times m$ ones matrix.

##  II Statement of the Control Problem

To reflect the practical implementation of digital controllers for physical
systems, we consider continuous-time dynamics under sampled-data control using
discrete-time predictive controllers. In particular, we consider the control
architecture shown in Figure 1, where ${\mathcal{M}}$ is the target
continuous-time system, $t\geq 0$, $u(t)\in{\mathbb{R}}^{m}$ is the control,
and $y(t)\in{\mathbb{R}}^{p}$ is the output of ${\mathcal{M}},$ which is
sampled to produce the measurement $y_{k}\in{\mathbb{R}}^{p},$ which, for all
$k\geq 0,$ is given by

$y_{k}\stackrel{{\scriptstyle\triangle}}{{=}}y(kT_{\rm s}),$ (1)

where $T_{\rm s}>0$ is the sample time.

The predictive controller, which is updated at each step $k,$ is denoted by
$G_{{\rm c},k}$. For all $k\geq 0,$ let $y_{{\rm
t},k}\stackrel{{\scriptstyle\triangle}}{{=}}C_{{\rm
t},k}y_{k}\in{\mathbb{R}}^{p_{\rm t}}$ be the command following output, where
$C_{{\rm t},k}\in{\mathbb{R}}^{p_{\rm t}\times p},$ and let
$r_{k}\in{\mathbb{R}}^{p_{\rm t}}$ be the command. The inputs to $G_{{\rm
c},k}$ are $r_{k},$ $y_{k},$ and $y_{{\rm t},k},$ and the output is the
requested discrete-time control $u_{{\rm req},k}\in{\mathbb{R}}^{m}.$ The
predictive controller uses $y_{k}$ for system identification, and $r_{k}$ and
$y_{{\rm t},k}$ trajectory command following. Since the response of a real
actuator is subjected to hardware constraints, the implemented discrete-time
control is

$u_{k}\stackrel{{\scriptstyle\triangle}}{{=}}\gamma_{k}(\sigma(u_{{\rm
req},k})),$ (2)

where $\sigma\colon{\mathbb{R}}^{m}\to{\mathbb{R}}^{m}$ is the control
magnitude saturation function

$\sigma(u)\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}\bar{\sigma}(u_{(1)})\\\
\vdots\\\ \bar{\sigma}(u_{(m)})\end{bmatrix},$ (3)

where $\bar{\sigma}\colon{\mathbb{R}}\to{\mathbb{R}}$ is defined by

$\bar{\sigma}(u)\stackrel{{\scriptstyle\triangle}}{{=}}\begin{cases}u_{\max},&u>u_{\max},\\\
u,&u_{\min}\leq u\leq u_{\max},\\\ u_{\min},&u<u_{\min},\end{cases}$ (4)

and $u_{\min},u_{\max}\in{\mathbb{R}}$ are the lower and upper magnitude
saturation levels, respectively, and
$\gamma_{k}\colon{\mathbb{R}}^{m}\to{\mathbb{R}}^{m}$ is the control rate
(move-size) saturation function

$\gamma_{k}(u_{k})\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}\bar{\gamma}_{k}(u_{k(1)})\\\
\vdots\\\ \bar{\gamma}_{k}(u_{k(m)})\end{bmatrix},$ (5)

where $\bar{\gamma}_{k}\colon{\mathbb{R}}\to{\mathbb{R}}$ is defined by

$\bar{\gamma}_{k}(u_{k})\stackrel{{\scriptstyle\triangle}}{{=}}\begin{cases}u_{k-1}+\Delta
u_{\max},&u_{k}-u_{k-1}>\Delta u_{\max},\\\ u_{k},&\Delta u_{\min}\leq
u_{k}-u_{k-1}\leq\Delta u_{\max},\\\ u_{k-1}+\Delta
u_{\min},&u_{k}-u_{k-1}<\Delta u_{\min},\end{cases}$ (6)

and $\Delta u_{\min},\Delta u_{\max}\in{\mathbb{R}}$ are the lower and upper
move-size saturation levels, respectively. The continuous-time control signal
$u(t)$ applied to the structure is generated by applying a zero-order-hold
operation to $u_{k},$ that is, for all $k\geq 0,$ and, for all $t\in[kT_{\rm
s},(k+1)T_{\rm s}),$

$u(t)=u_{k}.$ (7)

The objective of the predictive controller is to yield an input signal that
minimizes the norm of the difference between the command following output and
the command, that is, yield $u_{k}$ such that $\sum_{k=0}^{\infty}||y_{{\rm
t},k}-r_{k}||_{2}$ is minimized.

$G_{{\rm c},k}$$\sigma$$\gamma_{k}$ZOH$C_{\rm t}$${\mathcal{M}}$$u_{{\rm
req},k}$$u_{k}$$u(t)$$T_{\rm s}$$y(t)$$y_{k}$$r_{k}$$y_{{\rm t},k}$

Figure 1: Sampled-data implementation of predictive controller for control of
continuous-time system ${\mathcal{M}}.$ All sample-and-hold operations are
synchronous. The predictive controller $G_{{\rm c},k}$ generates the requested
discrete-time control $u_{{\rm req},k}\in{\mathbb{R}}^{m}$ at each step $k$.
The implemented discrete-time control is $u_{k}=\gamma_{k}(\sigma(u_{{\rm
req},k}))$, where $\sigma\colon{\mathbb{R}}^{m}\to{\mathbb{R}}^{m}$ represents
control-magnitude saturation and
$\gamma_{k}\colon{\mathbb{R}}^{m}\to{\mathbb{R}}^{m}$ represents move-size
saturation. The resulting continuous-time control $u(t)$ is generated by
applying a zero-order-hold operation to $u_{k}$. For this work,
${\mathcal{M}}$ represents a fixed-wing aircraft simulation model.

##  III Review of Predictive Cost Adaptive Control

PCAC combines online identification with output-feedback MPC. The PCAC
algorithm is presented in this section. Subsection III-A describes the
technique used for online identification, namely, RLS with variable-rate
forgetting based on the F-test [24]. Subsection III-B presents the block
observable canonical form (BOCF), which is used to represent the input-output
dynamics model as a state space model whose state is given explicitly in terms
of inputs, outputs, and model-coefficient estimates. Subsection III-C reviews
the MPC technique for receding-horizon optimization.

###  III-A Online Identification Using Recursive Least Squares with Variable-
Rate Forgetting Based on the F-Test

Let $\hat{n}\geq 0$ and, for all $k\geq 0,$ let $F_{{\rm
m},1,k},\ldots,F_{{\rm m},\hat{n},k}\in{\mathbb{R}}^{p\times p}$ and $G_{{\rm
m},1,k},\ldots,G_{{\rm m},\hat{n},k}\in{\mathbb{R}}^{p\times m}$ be the
coefficient matrices to be estimated using RLS. Furthermore, let
$\hat{y}_{k}\in{\mathbb{R}}^{p}$ be an estimate of $y_{k}$ defined by

$\hat{y}_{k}\stackrel{{\scriptstyle\triangle}}{{=}}-\sum_{i=1}^{\hat{n}}F_{{\rm
m},i,k}y_{k-i}+\sum_{i=1}^{\hat{n}}{G}_{{\rm m},i,k}u_{k-i},$ (8)

where

$\displaystyle y_{-\hat{n}}=\cdots=y_{-1}=0,$ (9) $\displaystyle
u_{-\hat{n}}=\cdots=u_{-1}=0.$ (10)

Using the identity ${\rm vec}(XY)=(Y^{\mathsf{T}}\otimes I){\rm vec}X,$ it
follows from (8) that, for all $k\geq 0,$

$\hat{y}_{k}=\phi_{k}\theta_{k},$ (11)

where

$\displaystyle\theta_{k}\stackrel{{\scriptstyle\triangle}}{{=}}$
$\displaystyle\ \begin{bmatrix}\theta_{F_{\rm m},k}^{\rm T}&\theta_{G_{\rm
m},k}^{\rm T}\end{bmatrix}^{\rm T}\in{\mathbb{R}}^{\hat{n}p(m+p)},$ (12)
$\displaystyle\theta_{F_{\rm m},k}\stackrel{{\scriptstyle\triangle}}{{=}}$
$\displaystyle\ {\rm vec}\begin{bmatrix}F_{{\rm m},1,k}&\cdots&F_{{\rm
m},\hat{n},k}\end{bmatrix}\in{\mathbb{R}}^{\hat{n}p^{2}},$ (13)
$\displaystyle\theta_{G_{\rm m},k}\stackrel{{\scriptstyle\triangle}}{{=}}$
$\displaystyle\ {\rm vec}\begin{bmatrix}G_{{\rm m},1,k}&\cdots&G_{{\rm
m},\hat{n},k}\end{bmatrix}\in{\mathbb{R}}^{\hat{n}pm},$ (14)
$\displaystyle\phi_{k}\stackrel{{\scriptstyle\triangle}}{{=}}$
$\displaystyle\begin{bmatrix}-y_{k-1}^{\rm T}&\cdots&-y_{k-\hat{n}}^{\rm
T}&u_{k-1}^{\rm T}&\cdots&u_{k-\hat{n}}^{\rm T}\end{bmatrix}\otimes I_{p}$
$\displaystyle\in{\mathbb{R}}^{p\times\hat{n}p(m+p)}.$ (15)

To determine the update equations for $\theta_{k}$, for all $k\geq 0$, define
$e_{k}\colon{\mathbb{R}}^{\hat{n}p(m+p)}\to{\mathbb{R}}^{p}$ by

$e_{k}(\bar{\theta})\stackrel{{\scriptstyle\triangle}}{{=}}y_{k}-\phi_{k}\bar{\theta},$
(16)

where $\bar{\theta}\in{\mathbb{R}}^{\hat{n}p(m+p)}.$ Using (11), the
identification error at step $k$ is defined by

$e_{k}(\theta_{k})=y_{k}-\hat{y}_{k}.$ (17)

For all $k\geq 0$, the RLS cumulative cost
$J_{k}\colon{\mathbb{R}}^{\hat{n}p(m+p)}\to[0,\infty)$ is defined by [20]

$J_{k}(\bar{\theta})\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{i=0}^{k}\frac{\rho_{i}}{\rho_{k}}e_{i}^{\rm
T}(\bar{\theta})e_{i}(\bar{\theta})+\frac{1}{\rho_{k}}(\bar{\theta}-\theta_{0})^{\rm
T}\psi_{0}^{-1}(\bar{\theta}-\theta_{0}),$ (18)

where $\psi_{0}\in{\mathbb{R}}^{\hat{n}p(m+p)\times\hat{n}p(m+p)}$ is positive
definite, $\theta_{0}\in{\mathbb{R}}^{\hat{n}p(m+p)}$ is the initial estimate
of the coefficient vector, and, for all $i\geq 0,$

$\rho_{i}\stackrel{{\scriptstyle\triangle}}{{=}}\prod_{j=0}^{i}\lambda_{j}^{-1}.$
(19)

For all $j\geq 0$, the parameter $\lambda_{j}\in(0,1]$ is the forgetting
factor defined by
$\lambda_{j}\stackrel{{\scriptstyle\triangle}}{{=}}\beta_{j}^{-1}$, where

$\beta_{j}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{cases}1,&j<\tau_{\rm
d},\\\ 1+\eta\bar{\beta}_{j},&j\geq\tau_{\rm d},\end{cases}$ (20)
$\bar{\beta}_{j}\stackrel{{\scriptstyle\triangle}}{{=}}g(e_{j-\tau_{\rm
d}}(\theta_{j-\tau_{\rm
d}}),\ldots,e_{j}(\theta_{j}))\textbf{1}\big{(}g(e_{j-\tau_{\rm
d}}(\theta_{j-\tau_{\rm d}}),\ldots,e_{j}(\theta_{j}))\big{)},$ (21)

and $\tau_{\rm d}>p$, $\eta>0$, $\textbf{1}\colon{\mathbb{R}}\to\\{0,1\\}$ is
the unit step function, and $g$ is a function of past RLS identification
errors which is given by (10) in [24] in the case where $p=1$ and (13) in [24]
in the case where $p>1.$ Note that $g$ includes forgetting terms based on the
inverse cumulative distribution function of the F-distribution and depends on
$\tau_{\rm d},$ $\tau_{\rm n}\in[p,\tau_{\rm d}),$ and significance level
$\alpha_{F}\in(0,1].$

Finally, for all $k\geq 0$, the unique global minimizer of $J_{k}$ is given by
[20]

$\theta_{k+1}=\theta_{k}+\psi_{k+1}\phi_{k}^{\rm
T}(y_{k}-\phi_{k}\theta_{k}),$ (22)

where

$\displaystyle\psi_{k+1}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\beta_{k}\psi_{k}-\beta_{k}\psi_{k}\phi_{k}^{\rm
T}(\tfrac{1}{\beta_{k}}I_{p}+\phi_{k}\psi_{k}\phi_{k}^{\rm
T})^{-1}\phi_{k}\psi_{k},$ (23)

and $\psi_{0}$ is the performance-regularization weighting in (18). Additional
details concerning RLS with forgetting based on the F-distribution are given
in [24].

###  III-B Input-Output Model and the Block Observable Canonical Form

Considering the estimate $\hat{y}_{k}$ of $y_{k}$ given by (8), it follows
that, for all $k\geq 0,$

$y_{k}\approx-\sum_{i=1}^{\hat{n}}F_{{\rm
m},i,k}y_{k-i}+\sum_{i=1}^{\hat{n}}{G}_{{\rm m},i,k}u_{k-i}.$ (24)

Viewing (24) as an equality, it follows that, for all $k\geq 0,$ the BOCF
state-space realization of (24) is given by [26]

$\displaystyle x_{{\rm m},k+1}$ $\displaystyle=A_{{\rm m},k}x_{{\rm
m},k}+B_{{\rm m},k}u_{k},$ (25) $\displaystyle y_{k}$ $\displaystyle=C_{\rm
m}x_{{\rm m},k},$ (26)

where

$\displaystyle A_{{\rm
m},k}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}-F_{{\rm
m},1,k+1}&I_{p}&\cdots&\cdots&0_{p\times p}\\\ -F_{{\rm m},2,k+1}&0_{p\times
p}&\ddots&&\vdots\\\ \vdots&{\vdots}&\ddots&\ddots&0_{p\times p}\\\
\vdots&\vdots&&\ddots&I_{p}\\\ -F_{{\rm m},\hat{n},k+1}&0_{p\times
p}&\cdots&\cdots&0_{p\times
p}\end{bmatrix}\in{\mathbb{R}}^{\hat{n}p\times\hat{n}p},$ (27) $\displaystyle
B_{{\rm m},k}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}G_{{\rm
m},1,k+1}\\\ G_{{\rm m},2,k+1}\\\ \vdots\\\ G_{{\rm
m},\hat{n},k+1}\end{bmatrix}\in{\mathbb{R}}^{\hat{n}p\times m},$ (28)
$\displaystyle C_{\rm
m}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}I_{p}&0_{p\times
p}&\cdots&0_{p\times p}\end{bmatrix}\in{\mathbb{R}}^{p\times\hat{n}p},$ (29)

and

$x_{{\rm m},k}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}x_{{\rm
m},k(1)}\\\ \vdots\\\ x_{{\rm
m},k(\hat{n})}\end{bmatrix}\in{\mathbb{R}}^{\hat{n}p},$ (30)

where

$\displaystyle x_{{\rm m},k(1)}\stackrel{{\scriptstyle\triangle}}{{=}}y_{k},$
(31)

and, for all $j=2,\ldots,\hat{n},$

$\displaystyle x_{{\rm m},k(j)}\stackrel{{\scriptstyle\triangle}}{{=}}$
$\displaystyle-\sum_{i=1}^{\hat{n}-j+1}F_{{\rm m},i+j-1,k+1}y_{k-i}$
$\displaystyle+\sum_{i=1}^{\hat{n}-j+1}G_{{\rm m},i+j-1,k+1}u_{k-i}.$ (32)

Note that multiplying both sides of (25) by $C_{\rm m}$ and using (26)–(32)
implies that, for all $k\geq 0,$

$\displaystyle y_{k+1}$ $\displaystyle=C_{\rm m}x_{{\rm m},k+1}$
$\displaystyle=C_{\rm m}(A_{{\rm m},k}x_{{\rm m},k}+B_{{\rm m},k}u_{k})$
$\displaystyle=-F_{{\rm m},1,k+1}x_{{\rm m},k(1)}+x_{{\rm m},k(2)}+G_{{\rm
m},1,k+1}u_{k}$ $\displaystyle=-F_{{\rm
m},1,k+1}y_{k}-\sum_{i=1}^{\hat{n}-1}F_{{\rm m},i+1,k+1}y_{k-i}$
$\displaystyle+\sum_{i=1}^{\hat{n}-1}G_{{\rm m},i+1,k+1}u_{k-i}+G_{{\rm
m},1,k+1}u_{k}$ $\displaystyle=-\sum_{i=1}^{\hat{n}}F_{{\rm
m},i,k+1}y_{k+1-i}+\sum_{i=1}^{\hat{n}}G_{{\rm m},i,k+1}u_{k+1-i},$ (33)

which is approximately equivalent to (24) with $k$ in (24) replaced by $k+1$.

###  III-C Model Predictive Control (MPC)

Let $\ell\geq 1$ be the horizon and, for all $k\geq 0$ and all
$i=1,\ldots,\ell,$ let $x_{{\rm m},k|i}\in{\mathbb{R}}^{\hat{n}p}$ be the
$i$-step predicted state, $y_{{\rm m},k|i}\in{\mathbb{R}}^{p}$ be the $i$-step
predicted output, and $u_{k|i}\in{\mathbb{R}}^{m}$ be the $i$-step predicted
control. Then, the $\ell$-step predicted output of (26) for a sequence of
$\ell$ future controls is given by

$\displaystyle Y_{\ell,k|1}=\Gamma_{\ell,k}x_{{\rm
m},k|1}+T_{\ell,k}U_{\ell,k|1},$ (34)

where

$\displaystyle
Y_{\ell,k|1}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}y_{{\rm
m},k|1}\\\ \vdots\\\ y_{{\rm m},k|\ell}\end{bmatrix}\in{\mathbb{R}}^{\ell
p},\quad
U_{\ell,k|1}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}u_{k|1}\\\
\vdots\\\ u_{k|\ell}\end{bmatrix}\in{\mathbb{R}}^{\ell m},$ (35)
$\Gamma_{\ell,k}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}C_{\rm
m}^{\rm T}&(C_{\rm m}A_{{\rm m},k})^{\rm T}&\cdots&(C_{\rm m}A_{{\rm
m},k}^{\ell-1})^{\rm T}\end{bmatrix}^{\rm T}\in{\mathbb{R}}^{\ell
p\times\hat{n}p},$ (36) $\displaystyle T_{\ell,k}$
$\displaystyle\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}0_{p\times
m}&\cdots&\cdots&\cdots&\cdots&0_{p\times m}\\\ H_{k,1}&0_{p\times
m}&\cdots&\cdots&\cdots&0_{p\times m}\\\ H_{k,2}&H_{k,1}&0_{p\times
m}&\cdots&\cdots&0_{p\times m}\\\ H_{k,3}&H_{k,2}&H_{k,1}&0_{p\times
m}&\cdots&0_{p\times m}\\\ \vdots&\vdots&\vdots&\ddots&\ddots&0_{p\times m}\\\
H_{k,\ell-1}&H_{k,\ell-2}&H_{k,\ell-3}&\cdots&H_{k,1}&0_{p\times
m}\end{bmatrix}$ $\displaystyle\in{\mathbb{R}}^{\ell p\times\ell m},$ (37)

where $H_{k,i}\stackrel{{\scriptstyle\triangle}}{{=}}C_{\rm m}A_{{\rm
m},k}^{i-1}B_{{\rm m},k}\in{\mathbb{R}}^{p\times m}$ for all
$i=1,\ldots,\ell-1.$

Let
$\mathcal{R}_{\ell,k}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}r_{k+1}^{\rm
T}\ \cdots\ r_{k+\ell}^{\rm T}\end{bmatrix}^{\rm T}\in{\mathbb{R}}^{\ell
p_{\rm t}}$ be a vector composed of $\ell$ future commands, let $y_{{\rm
m},{\rm t},k|i}\stackrel{{\scriptstyle\triangle}}{{=}}C_{\rm t}y_{{\rm
m},k|i}\in{\mathbb{R}}^{p_{\rm t}}$ be the $i$-step predicted command-
following output, let $Y_{{\rm
t},\ell,k|1}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}y_{{\rm
m},{\rm t},k|1}^{\rm T}&\cdots&y_{{\rm m},{\rm t},k|\ell}^{\rm
T}\end{bmatrix}^{\rm T}=C_{{\rm t},\ell}Y_{\ell,k|1}\in{\mathbb{R}}^{\ell
p_{\rm t}},$ where $C_{{\rm
t},\ell}\stackrel{{\scriptstyle\triangle}}{{=}}I_{\ell}\otimes C_{\rm
t}\in{\mathbb{R}}^{\ell p_{\rm t}\times\ell p},$ and define

$\displaystyle\Delta
U_{\ell,k|1}\stackrel{{\scriptstyle\triangle}}{{=}}\begin{bmatrix}u_{k|1}-u_{k}\\\
u_{k|2}-u_{k|1}\\\ \vdots\\\
u_{k|\ell}-u_{k|\ell-1}\end{bmatrix}\in{\mathbb{R}}^{\ell m}.$ (38)

Then, the receding horizon optimization problem is given by

$\displaystyle\min_{U_{\ell,k|1}}\left(Y_{{\rm
t},\ell,k|1}-\mathcal{R}_{\ell,k}\right)^{\rm T}$ $\displaystyle
Q\left(Y_{{\rm t},\ell,k|1}-\mathcal{R}_{\ell,k}\right)$
$\displaystyle\quad+\Delta U_{\ell,k|1}^{\rm T}R\Delta U_{\ell,k|1},$ (39)

subject to

$\displaystyle U_{\min}\leq U_{\ell,k|1}$ $\displaystyle\leq U_{\max}$ (40)
$\displaystyle\Delta U_{\min}\leq\Delta U_{\ell,k|1}$ $\displaystyle\leq\Delta
U_{\max},$ (41)

where $Q\triangleq\begin{bmatrix}\bar{Q}&0_{p_{\text{t}}\times
p_{\text{t}}}\\\ 0_{p_{\text{t}}\times
p_{\text{t}}}&\bar{P}\end{bmatrix}\in\mathbb{R}^{\ell p_{\text{t}}\times\ell
p_{\text{t}}}$ is the positive-definite output weighting,
$\bar{Q}\in\mathbb{R}^{(\ell-1)p_{\text{t}}\times(\ell-1)p_{\text{t}}}$ is the
positive-definite cost-to-go output weighting,
$\bar{P}\in\mathbb{R}^{p_{\text{t}}\times p_{\text{t}}}$ is the positive-
definite terminal output weighting, $R\in{\mathbb{R}}^{\ell m\times\ell m}$ is
the positive definite control move-size weight,
$U_{\min}\stackrel{{\scriptstyle\triangle}}{{=}}1_{\ell}\otimes
u_{\min}\in{\mathbb{R}}^{\ell m}$,
$U_{\max}\stackrel{{\scriptstyle\triangle}}{{=}}1_{\ell}\otimes
u_{\max}\in{\mathbb{R}}^{\ell m}$, $\Delta
U_{\min}\stackrel{{\scriptstyle\triangle}}{{=}}1_{\ell}\otimes\Delta
u_{\min}\in{\mathbb{R}}^{\ell m}$, and $\Delta
U_{\max}\stackrel{{\scriptstyle\triangle}}{{=}}1_{\ell}\otimes\Delta
u_{\max}\in{\mathbb{R}}^{\ell m}$.

In summary, at each time step, online identification is performed to find
input-output model coefficients $\theta_{k+1}$, which are then used to create
a state space realization $\left(A_{{\rm m},k},B_{{\rm m},k},C_{\rm
m}\right)$. Then, the state-space realization is used in a receding horizon
optimization problem to solve for the $\ell$-step controls $U_{\ell,k|1}.$ The
control input for the next step is then given by $u_{k|1}$, and the rest of
the components of $U_{\ell,k|1}$ are discarded.

##  IV Adaptive Control of a Linearized 6DOF Aircraft Model

In this section, a linearized 6DOF fixed-wing RC aircraft model from Athena
Vortex Lattice (AVL) is considered for both longitudinal and lateral control.
The first example uses altitude, elevation angle, and bank angle measurements
for altitude and bank-angle command following using the elevator and aileron.
The second example uses altitude, elevation, bank, and azimuth measurements
for altitude and azimuth-angle command following using all three control
surfaces. Note that all numerical examples in this section use a sampling rate
of $T_{s}=0.05$ sec/step and the aircraft model is initialized at an altitude
$h_{0}=500$ m and airspeed $V_{0}=10$ m/s with all other initial model
parameters set to 0.

###  IV-A Altitude and Bank-Angle Command Following

For this example, two PCAC loops are used to follow altitude and bank-angle
commands. The longitudinal loop uses the altitude $h$ and elevation angle
$\Theta$ to identify the longitudinal dynamics and to specify the elevator
deflection $\delta e$. The lateral loop uses the bank angle $\Phi$ to identify
the lateral dynamics and to specify the aileron deflection $\delta a$. This
control architecture is illustrated in Figure 2.

Figure 2: Command-following block diagram for altitude $h$ and bank angle
$\Phi$. This architecture uses two PCAC controllers, namely, PCAC/long for the
longitudinal dynamics and PCAC/lat for the lateral dynamics. Note that
PCAC/long uses the elevation angle $\Theta$ for system identification.

PCAC/long is initialized with $\psi_{\rm long,0}=10^{5}I_{36}$, $\theta_{\rm
long,0}=0.01\textbf{1}_{36\times 1}$, $\hat{n}=6$, $\ell=40$,
$\bar{Q}=0.1I_{39}$, $\bar{P}=0.1$, $R=0.1$, $|\Delta\delta e|\leq\Delta\delta
e_{\rm max}=0.5$ deg/step, and $|\delta e|\leq\delta e_{\rm max}=10$ deg with
no forgetting. PCAC/lat is initialized with $\psi_{\rm lat,0}=10^{5}I_{36}$,
$\theta_{\rm lat,0}=0.01\textbf{1}_{36\times 1}$, $\hat{n}=6$, $\ell=40$,
$\bar{Q}=0.1I_{39}$, $\bar{P}=0.1$, $R=0.01$, $|\Delta\delta
a|\leq\Delta\delta a_{\rm max}=0.5$ deg/step, and $|\delta a|\leq\delta a_{\rm
max}=2$ deg with no forgetting. The resulting flight path data is shown in
Figure 3. An altitude ramp command is given for $t\in[25,50]$ s, followed by a
bank-angle ramp command for $t\in[50,60]$. Furthermore, altitude and bank-
angle ramp commands are given simultaneously for $t\in[75,85]$ s. As shown in
Figure 3, PCAC follows the altitude and bank-angle commands with a slight
overshoot in the altitude command following.

Figure 3: Altitude and bank-angle command following. (a) compares the command
$h_{\rm cmd}$ and the simulated $h$. (b) compares the command $\Phi_{\rm cmd}$
and the simulated $\Phi$. Note that learning takes place for the longitudinal
dynamics for $t\in[0,5]$ s causing the excitation seen in (a).

Figures 4 and 5 show the control-surface deflections $\delta e$ and $\delta a$
as well as the estimated coefficient vectors $\theta_{\rm long}$ and
$\theta_{\rm lat}$, where $\theta_{\rm long}$ is the estimated coefficient
vector of PCAC/long and $\theta_{\rm lat}$ is the estimated coefficient vector
of PCAC/lat. $\diamond$

Figure 4: Aileron and elevator deflections for altitude and bank-angle command
following. (a) shows the elevator deflection $\delta e$ (blue) with the
constraints on $|\delta e|$ (dashed black). (b) shows the aileron deflection
$\delta a$ (blue) with the constraints on $|\delta a|$ (dashed black).

Figure 5: Estimated coefficient vectors for the altitude command-following
loop $\theta_{\rm long}$ (a) and the bank-angle command-following loop
$\theta_{\rm lat}$ (b). No prior modeling information is assumed, where
$\theta_{\rm long,0}=\theta_{\rm lat,0}=0.01\textbf{1}_{36\times 1}$.

###  IV-B Altitude and Azimuth-Angle Command Following

For this example, two PCAC loops are used to follow altitude and azimuth-angle
commands. The longitudinal loop uses the altitude $h$ and elevation angle
$\Theta$ to identify the longitudinal dynamics and to specify the elevator
deflection $\delta e$. The lateral loop uses the bank angle $\Phi$ and azimuth
angle $\Psi$ to identify the lateral dynamics and to specify the aileron and
rudder deflections $\delta a$ and $\delta r$. This control architecture is
illustrated in Figure 6.

Figure 6: Command-following block diagram for altitude $h$ and azimuth angle
$\Psi$. This architecture uses two PCAC controllers, namely, PCAC/long for the
longitudinal dynamics and PCAC/lat for the lateral dynamics. Note that
PCAC/long uses the elevation angle $\Theta$ for system identification and
PCAC/lat uses the bank angle $\Phi$ for system identification.

The same PCAC/long parameters from Section IV-A were used to initialize the
PCAC/long loop in this example. PCAC/lat is initialized with $\psi_{\rm
lat,0}=2\times 10^{7}I_{48}$, $\theta_{\rm lat,0}=0.01\textbf{1}_{48\times
1}$, $\hat{n}=6$, $\ell=50$, $\bar{Q}=60I_{49}$, $\bar{P}=60$,
$R=1.5I_{2\times 2}$, $[|\Delta\delta a|,\ |\Delta\delta r|]^{\rm
T}\leq[\Delta\delta a_{\rm max},\ \Delta\delta r_{\rm max}]^{\rm T}=[0.2,\
0.2]^{\rm T}$ deg/step, and $[|\delta a|,\ |\delta r|]^{\rm T}\leq[\delta
a_{\rm max},\ \delta r_{\rm max}]^{\rm T}=[5,\ 5]^{\rm T}$ deg with no
forgetting. The resulting flight path data is shown in Figure 7. Altitude ramp
commands are given for $t\in[20,40]$ s and $t\in[60,80]$ s while azimuth-angle
ramp commands are given for $t\in[20,45]$ s and $t\in[55,80]$ s. As shown in
Figure 7, PCAC follows the altitude and azimuth-angle commands with slight
overshoot in the altitude command following.

Figure 7: Altitude and azimuth-angle command following. (a) compares the
command $h_{\rm cmd}$ and the simulated $h$. (b) compares the command
$\Psi_{\rm cmd}$ and the simulated $\Psi$. Note that learning takes place for
the longitudinal dynamics for $t\in[0,5]$ s causing the excitation seen in
(a), while learning for the lateral dynamics takes place for $t\in[20,25]$ s,
causing the excitation seen in (b).

Figures 8 and 9 show the control surface deflections $\delta e$, $\delta a$,
and $\delta r$ as well as the estimated coefficient vectors $\theta_{\rm
long}$ and $\theta_{\rm lat}$, where $\theta_{\rm long}$ is the estimated
coefficient vector of PCAC/long and $\theta_{\rm lat}$ is the estimated
coefficient vector of PCAC/lat. $\diamond$

Figure 8: Elevator, aileron, and rudder deflections for altitude and azimuth-
angle command following. (a) shows the elevator deflection $\delta e$ (blue)
with the constraint on $|\delta e|$ (dashed black). (b) shows the aileron
deflection $\delta a$ (blue) with the constraint on $|\delta a|$ (dashed
black). (c) shows the rudder deflection $\delta r$ (blue) with the constraint
$|\delta r|$ (dashed black).

Figure 9: Estimated coefficient vectors for the altitude command-following
loop $\theta_{\rm long}$ (a) and the azimuth-angle command-following loop
$\theta_{\rm lat}$ (b). No prior modeling information is assumed, where
$\theta_{\rm long,0}=0.01\textbf{1}_{39\times 1}$ and $\theta_{\rm
lat,0}=0.01\textbf{1}_{49\times 1}$.

##  V Adaptive Control of a Nonlinear 3DOF Aircraft Model

In this section, a nonlinear 3DOF fixed-wing passenger aircraft model provided
by Matlab’s aerospace toolbox and featured in [27] is considered for
longitudinal control. The first example uses altitude and angle of attack
measurements for altitude command following using the elevator. The second
example uses altitude, angle of attack, and airspeed measurements for altitude
and airspeed command following using the elevator and throttle. Note that all
simulations in this section use a sampling rate of $T_{s}=0.05$ sec/step and
the aircraft model is initialized at an altitude $h_{0}=2000$ m and airspeed
$V_{0}=85$ m/s with all other initial model parameters set to 0. Both examples
use forgetting for altitude command following with $\eta=0.025$, $\tau_{\rm
n}=40$, $\tau_{\rm d}=200$, and $\alpha_{F}=0.001$.

###  V-A Altitude Command Following

For this example, a single PCAC loop is used to follow altitude commands. The
loop uses the altitude $h$ and angle of attack $\alpha$ to identify the
longitudinal dynamics and generate an elevator deflection $\delta e$. This
control architecture is illustrated in Figure 10.

Figure 10: Command-following block diagram for altitude $h$. Note that the
loop uses the angle of attack $\alpha$ for system identification.

PCAC/long is initialized with $\psi_{\rm long,0}=10^{4}I_{60}$, $\theta_{\rm
long,0}=0.01\textbf{1}_{60\times 1}$, $\hat{n}=10$, $\ell=40$,
$\bar{Q}=0.01I_{39}$, $\bar{P}=500$, $R=0.01$, $|\Delta\delta
e|\leq\Delta\delta e_{\rm max}=1$ deg/step, and $|\delta e|\leq\delta e_{\rm
max}=5$ deg. The resulting flight path data is shown in Figure 11. Altitude
ramp commands are given for $t\in[30,50]$ s and $t\in[75,90]$ s. As shown in
Figure 11, PCAC follows the altitude commands with slight overshoot and
undershoot.

Figure 11: Altitude command following. The figure compares the command $h_{\rm
cmd}$ and the simulated $h$. Note that learning takes place for $t\in[0,5]$ s
causing the initial excitation seen in the plot.

Figure 12 shows the control surface deflection $\delta e$ and estimated
coefficient vector $\theta_{\rm long}$. $\diamond$

Figure 12: Elevator deflection and estimated coefficient vector for altitude
command following. (a) shows the elevator deflection $\delta e$ (blue) with
the constraint on $|\delta e|$ (dashed black). (b) shows the estimated
coefficient vector $\theta_{\rm long}$. No prior modeling information is
assumed, where $\theta_{\rm long,0}=0.01\textbf{1}_{60\times 1}$.

###  V-B Altitude and Airspeed Command Following

For this example, two PCAC loops are used to follow altitude and airspeed
commands. The longitudinal loop uses the altitude $h$ and angle of attack
$\alpha$ to identify the longitudinal dynamics and to specify the elevator
deflection $\delta e$. The airspeed loop uses the airspeed $V$ to identify the
throttle dynamics and generate a change in throttle $\delta T$. Note that the
throttle parameter $T\in[0,1]$ has a constant input of 0.5. This control
architecture is illustrated in Figure 13.

Figure 13: Command-following block diagram for altitude $h$ and airspeed $V$.
This architecture uses two PCAC controllers, namely, PCAC/long for the
longitudinal dynamics and PCAC/vel for the throttle dynamics. Note that
PCAC/long uses the angle of attack $\alpha$ for system identification.

The same PCAC/long parameters from Section V-A were used to initialize the
PCAC/long altitude command-following loop. PCAC/vel is initialized with
$\psi_{\rm vel,0}=10^{5}I_{8}$, $\theta_{\rm vel,0}=0.01\textbf{1}_{8\times
1}$, $\hat{n}=4$, $\ell=30$, $\bar{Q}=0.0001I_{29}$, $\bar{P}=0.01$, $R=0.01$,
$|\Delta\delta T|\leq\Delta\delta T_{\rm max}=0.1$ 1/step, and $|\delta
T|\leq\delta T_{\rm max}=0.25$ with no forgetting. The resulting flight path
data is shown in Figure 14. Airspeed ramp commands are given for $t\in[10,45]$
s followed by an altitude ramp command for $t\in[70,95]$ s. Furthermore,
altitude and airspeed commands are given simultaneously for $t\in[125,145]$ s.
As shown in Figure 14, PCAC follows the altitude and airspeed commands with
slight overshoot in the altitude command following.

Figure 14: Altitude and airspeed command following. (a) compares the command
$h_{\rm cmd}$ and the simulated $h$. (b) compares the command $V_{\rm cmd}$
and the simulated data $V$. Note that learning takes place for both the
longitudinal dynamics and throttle dynamics for $t\in[0,5]$ s causing the
initial excitations seen in the plots.

Figures 15 and 16 show the control surface deflection $\delta e$, throttle
changes $\delta T$, and estimated coefficient vectors $\theta_{\rm long}$ and
$\theta_{\rm vel}$, where $\theta_{\rm long}$ is the estimated coefficient
vector of PCAC/long and $\theta_{\rm vel}$ is the estimated coefficient vector
of PCAC/vel. $\diamond$

Figure 15: Elevator deflection and throttle change for altitude and airspeed
command following. (a) shows the elevator deflection $\delta e$ (blue) with
the constraint on $|\delta e|$ (dashed black). (b) shows the throttle change
$\delta T$ (blue) with the constraint on $|\delta T|$ (dashed black).

Figure 16: Estimated coefficient vectors for the altitude command-following
loop $\theta_{\rm long}$ (a) and the airspeed command-following loop
$\theta_{\rm vel}$ (b). No prior modeling information is assumed, where
$\theta_{\rm long,0}=0.01\textbf{1}_{60\times 1}$ and $\theta_{\rm
vel,0}=0.01\textbf{1}_{8\times 1}$.

##  VI Conclusions and future work

For autonomous control of aircraft with unmodeled aerodynamics, this paper
proposed an autopilot based on predictive cost adaptive control (PCAC), which
is an indirect adaptive control extension of model predictive control. PCAC
uses recursive least squares (RLS) with variable-rate forgetting for online,
closed-loop system identification, and receding-horizon optimization based on
quadratic programming. This technique was demonstrated numerically on a 6DOF
linearized aircraft model and a 3DOF nonlinear aircraft model. In both cases,
the adaptive autopilot was able to follow attitude, bank-angle, azimuth-angle,
and velocity commands over a range of operation.

The ability of PCAC to operate as an adaptive autopilot without aerodynamic
modeling has useful implications in practice. First and foremost, this
technique can mitigate the need for extensive wind tunnel testing. Second, it
can help avoid the need for gain scheduling and dynamic inversion. Finally,
PCAC can be used to accelerate the expensive and time-consuming
aircraft/autopilot design cycle. This potential will be investigated in future
research through embedded control of fixed-wing, autonomous aircraft.

## References

  * [1] A. Tewari, _Advanced Control of Aircraft, Spacecraft and Rockets_. Wiley, 2011.
  * [2] R. F. Stengel, _Flight Dynamics_ , 2nd ed. Princeton, 2022.
  * [3] N. Hovakimyan and C. Cao, _$L_{1}$ Adaptive Control Theory: Guaranteed Robustness with Fast Adaptation_. SIAM, 2010.
  * [4] E. Lavretsky and K. Wise, _Robust and Adaptive Control: With Aerospace Applications_. Springer, 2012.
  * [5] A. Ansari and D. S. Bernstein, “Retrospective Cost Adaptive Control of the Generic Transport Model Under Uncertainty and Failure,” _J. Aerospace Information Systems_ , vol. 14, no. 3, pp. 123–174, 2017.
  * [6] W. Kwon and S. Han, _Receding Horizon Control: Model Predictive Control for State Models_. Springer, 2006.
  * [7] E. F. Camacho and C. Bordons, _Model Predictive Control_ , 2nd ed. Springer, 2007.
  * [8] U. Eren, A. Prach, B. B. Koçer, S. V. Raković, E. Kayacan, and B. Açıkmeşe, “Model predictive control in aerospace systems: Current state and opportunities,” _Journal of Guidance, Control, and Dynamics_ , vol. 40, no. 7, pp. 1541–1566, 2017.
  * [9] J. B. Rawlings, D. Q. Mayne, and M. M. Diehl, _Model Predictive Control: Theory, Computation, and Design_. Nob Hill, 2020.
  * [10] H. J. van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: a new perspective on data-driven analysis and control,” _IEEE Trans. Autom. Contr._ , vol. 65, no. 11, pp. 4753–4768, 2020.
  * [11] S. A. U. Islam, T. Nguyen, I. Kolmanovsky, and D. S. Bernstein, “Data-Driven Retrospective Cost Adaptive Control for Flight Control Applications,” _J. Guid. Contr. Dyn._ , vol. 44, pp. 1732–1758, October 2021.
  * [12] R. Soloperto, M. A. Müller, and F. Allgöwer, “Guaranteed closed-loop learning in model predictive control,” _IEEE Trans. Autom. Contr._ , vol. 68, pp. 991–1006, 2022.
  * [13] A. A. Feldbaum, “Dual Control Theory, Parts I and II,” _Autom. Rem. Contr._ , vol. 21, pp. 874–880, 1033–1039, 1961.
  * [14] H. Hjalmarsson, M. Gevers, F. De Bruyne, and J. Leblond, “Identification for Control: Closing the Loop Gives More Accurate Controllers,” in _Proc. IEEE Conf. Dec. Contr._ , Lake Buena Vista, FL, December 1994, pp. 4150–4155.
  * [15] U. Forssell and L. Ljung, “Closed-Loop Identification Revisited,” _Automatica_ , vol. 35, no. 7, pp. 1215–1241, 1999.
  * [16] N. M. Filatov and H. Unbehauen, _Adaptive Dual Control: Theory and Applications_. Springer, 2004.
  * [17] A. Mesbah, “Stochastic model predictive control with active uncertainty learning: A Survey on dual control,” _Ann. Rev. Contr._ , vol. 45, pp. 107–117, 2018.
  * [18] T. W. Nguyen, S. A. U. Islam, D. S. Bernstein, and I. V. Kolmanovsky, “Predictive Cost Adaptive Control: A Numerical Investigation of Persistency, Consistency, and Exigency,” _IEEE Contr. Sys. Mag._ , vol. 41, pp. 64–96, December 2021.
  * [19] T. W. Nguyen, I. V. Kolmanovsky, and D. S. Bernstein, “Sampled-data output-feedback model predictive control of nonlinear plants using online linear system identification,” in _Proc. Amer. Contr. Conf._ , 2021, pp. 4682–4687.
  * [20] S. A. U. Islam and D. S. Bernstein, “Recursive least squares for real-time implementation,” _IEEE Contr. Syst. Mag._ , vol. 39, no. 3, pp. 82–85, 2019.
  * [21] A. L. Bruce, A. Goel, and D. S. Bernstein, “Convergence and consistency of recursive least squares with variable-rate forgetting,” _Automatica_ , vol. 119, p. 109052, 2020.
  * [22] A. Goel, A. L. Bruce, and D. S. Bernstein, “Recursive least squares with variable-direction forgetting: Compensating for the loss of persistency,” _IEEE Contr. Sys. Mag._ , vol. 40, no. 4, pp. 80–102, 2020.
  * [23] A. L. Bruce, A. Goel, and D. S. Bernstein, “Necessary and sufficient regressor conditions for the global asymptotic stability of recursive least squares,” _Sys. Contr. Lett._ , vol. 157, pp. 1–7, 2021, article 105005.
  * [24] N. Mohseni and D. S. Bernstein, “Recursive least squares with variable-rate forgetting based on the F-test,” in _Proc. Amer. Contr. Conf._ , 2022, pp. 3937–3942.
  * [25] W. H. Kwon and A. E. Pearson, “On feedback stabilization of time-varying discrete linear systems,” _IEEE Trans. Autom. Contr._ , vol. AC-23, no. 3, pp. 479–481, 1978.
  * [26] J. W. Polderman, “A state space approach to the problem of adaptive pole assignment,” _Mathematics of Control, Signals and Systems_ , vol. 2, no. 1, pp. 71–94, 1989.
  * [27] MathWorks Inc., “Lightweight airplane design,” Natick, Massachusetts, United States, 2023. [Online]. Available: https://www.mathworks.com/help/aeroblks/lightweight-airplane-design.html