# Unconventional superconductivity after the BCS paradigm and empirical rules
for the exploration of high temperature superconductors

Hai-Hu Wen  Center for Superconducting Physics and Materials, National
Laboratory of Solid State Microstructures and Department of Physics, Nanjing
University, Nanjing 210093, China<EMAIL_ADDRESS>
###### Abstract

Superconducting state is achieved through quantum condensation of Cooper pairs
which are new types of charge carriers other than single electrons in normal
metals. The theory established by Bardeen-Cooper-Schrieffer (BCS) in 1957 can
successfully explain the phenomenon of superconductivity in many single-
element and alloy superconductors. Within the BCS scheme, the Cooper pairs are
formed by exchanging the virtual vibrations of lattice (phonons) between two
electrons with opposite momentum near the Fermi surface. The BCS theory has
dominated the field of superconductivity over 64 years. Many superconductors
discovered in past four decades, such as the heavy Fermion superconductors,
cuprates, iron pnictide/chalcogenide and nickelates seem, however, to strongly
violate the BCS picture. The most important issue is that, perhaps the BCS
picture based on electron-phonon coupling are the special case for
superconductivity, there are a lot of other reasons or routes for the Cooper
pairing and superconductivity. In this short overview paper, we will summarize
part of these progresses and try to guide readers to some new possible schemes
of superconductivity after the BCS paradigm. We also propose several empirical
rules for the exploration of high-temperature unconventional superconductors.

##  I Introduction

Superconductivity is a very interesting quantum state which exhibits zero
resistivity and expulsion of magnetic field (Meissner effect). It has been
well understood that this fascinating phenomenon is induced by the
condensation of electron pairs (termed as Cooper pairs) at a finite
temperature. Up to date, the superconducting state can only be achieved with
these electronic Cooper pairs although some other novel pictures concerning
other type of charged bosons have been proposed. Since its first discovery in
April 1911 by Kamerling Onnes et al., the field of superconductivity has
experienced glorious developments and achievements in past 111 years. Owing to
the deep and interesting physics involved by the phenomenon of
superconductivity, and the potential of large scale applications, the field
has been pushed forward step by step, most of time by the discovery of new
superconductors. We can expect more interesting phenomena and materials, which
will eventually lead to the revolution, both on the fundamental understanding
of condensed matter physics and bring in the dreamful industrial applications
of superconductors.

The development of the field of superconductivity is however not
straightforward, it is always mixed with some excitement and frustrations. In
Fig.1, we present the time dependence of the superconducting transition
temperatures Tc of some representative superconducting systems when they were
discovered. One can see that, the value of Tc ramps up with a small slope
before October 1986. After the tremendous efforts in 75 years, the Tc finally
reached about 23.2 K in $Nb_{3}Ge$ system. Most superconductors discovered
before this time were single element metals or alloys. The basic reason for
superconductivity, called as the superconducting mechanism, was unraveled in
1957 by Bardeen, Cooper and Schrieffer and named as the BCS theoryBCS . In
this short review, we omit the description about the BCS paradigm. A short
review can be found in our earlier review paperWenHHJPS2014 . Since the field
is rapidly moving forward, the understanding on the so-called unconventional
superconductivity have not come to the mature status yet. Thus this paper
serves only as a brief introduction and may be taken as a supplementary to
other beautiful reviews in this fascinating field. Nowadays, it is known that
either in conventional or unconventional superconductors, we have an energy
gap that protects the superconducting condensate from breaking Cooper pairs or
exciting so-called quasiparticles. In the BCS theory, the gap is given by

$\Delta_{s}=2\hbar\omega_{D}exp^{-1/N(0)V}.$ (1)

Here $\hbar\omega_{D}$ = $k_{B}\Theta_{D}$ is the Debye energy; N(0) is the
density of states (DOS) near the Fermi energy; $V$ is the interaction between
two electrons involved in the pair-scattering process. In the weak coupling
case, it is predicted that $\Delta_{s}=1.75k_{B}T_{c}$. Therefore the
superconducting state can be achieved at a finite temperature. One can see
from the above description that, one of the basic conditions in the
establishment of the BCS scheme is that $\Delta_{s}<<\hbar\omega_{D}<<E_{F}$.
Once we have a system with $\Delta_{s}$ comparable with the Fermi energy
$E_{F}$, this basic requirement is thus lost and we run into a new scheme of
superconductivity. For superconductors with conventional superconducting
pairing mechanism, the one with the highest transition temperature at ambient
pressure is $MgB_{2}$ which was discovered in 2001MgB2 . In some pressurized
hydridesHydride , it was reported that the superconducting transition
temperature rises to about 200K or beyond, but the pairing mechanism can still
be described by the phonon mediated BCS picture.

Figure 1: Correlation between the superconducting transition temperatures and
the year for them to be found. The three horizontal red dashed lines represent
the boiling temperatures of liquid helium, hydrogen, and nitrogen. The
circles, squares, diamonds, stars represent the conventional, cuprate, iron-
based and nickelate superconductors, respectively.

##  II Unconventional superconductivity in cuprates and iron
pnictides/chacogenides

In past four decades, many new superconductors were discovered. These systems
exhibit some features which are clearly at odds with the BCS theory. These
superconductors include the heavy FermionHeavyFermion , organicOrganic ,
cuprates, iron based superconductors and nickelates. Among them, the cuprate
is a typical one which was firstly discovered by K. A. Müller and J. G.
BednorzCuprate . They found that the superconductivity occurred at about 35 K
in a typical copper oxide $La_{2-x}Ba_{x}CuO_{4}$. This discovery ignited an
explosion in the field of superconductivity. In later 7-8 years, many other
new superconducting systems were discovered, for example the
$YBa_{2}Cu_{3}O_{7-\delta}$ ($T_{c}\approx 90K$)ChuCW ; ZhaoZX , the bismuth
family $Bi_{2}Sr_{2}CaCu_{2}O_{8}$ ($T_{c}\approx 90K$) and
$Bi_{2}Sr_{2}Ca_{2}Cu_{3}O_{10}$ ($T_{c}\approx 125K$)BiSystem , the thallium
familyTlSystem  and the mercury family HgSystem  etc. So far the highest
superconducting transition temperature in cuprates is Tc = 164 K, which was
discovered in the $HgBa_{2}Ca_{2}Cu_{3}O_{8+\delta}$ under a high pressure of
45 GPa.164K .

Figure 2: Sketch of the $CuO_{2}$ planes of the cuprate superconductors. The
small filled circles with arrows represent the Cu ions with the $3d^{9}$ as
the outermost orbital of electrons. The large circles represent the oxygen
ions with the 2p6 outermost orbitals.  Figure 3: A general electronic phase
diagram of the cuprate superconductors. The central yellow region shows the
parent phase with a long range ordered antiferromagnetic insulating phase. The
two dome like highlighted by red color show the regions for superconductivity.
In the electron doped side, the AF order may enter the superconducting dome,
showing a coexistence of superconductivity and AF order. The pairing symmetry
has been proved to be d-wave in the hole doped region, it may be non-monotonic
d-wave or s-wave pairing in the electron doped region. The solid and dashed
lines represent the pseudogap temperatures according to different definitions
or experimental measurements.

No sooner after the discovery by Müller and Bednorz, it was found that the
major players for superconductivity in the cuprates are the $CuO_{2}$ planes
in the middle of the pyramid constructed by the octahedra of Cu and oxygen
atoms. This plane is sketched in Fig.2. If we have a look at the formulae of
some parent phase of cuprate systems, such as $La_{2}CuO_{4}$, a natural
counting on the ionic state would lead to $Cu^{2+}$ with 3$d^{9}$, and
$O^{2-}$ with 2$p^{6}$ as the outermost shell of electron orbitals. Taking the
structural parameters of the system into consideration, a simple density-
functional-theory (DFT) calculation will find out that the material should be
a metal with $Cu-3d$ orbitals as the conduction band, forming a large electron
like Fermi surface surrounding the $\Gamma$ point, while later experiments
show that it is an insulator due to Mottness, thus it is called as a Mott
insulatorAndersonMott . Many transport properties indicate that the
resistivity close to the undoped case cannot be described by a band insulating
picture, rather it is described by a law $\rho(T)\propto log(1/T)$ which
manifests itself some exotic reasons. Almost simultaneously, researchers found
that the undoped phase, here called as the parent phase of the cuprate
superconductors, has a long range anti-ferromagnetic order with an ordering
temperature at around 300 K and the super-exchange energy $J\approx 150meV$.
Clearly the insulating behavior of the parent phase is due to the strong
correlation effect. A Hubbard model based picture, namely the Zhang-Rice
singletZhangRice  model, was proposed to interpret the electric conduction in
the doped cuprates. In the cuprate system, due to its crystalline field
effect, the $3d_{x^{2}-y^{2}}$ orbital is lifted to the highest energy and is
half filled. However, because of the strong Hubbard U (6-8 eV), this half-
filled band will be split further into lower and upper Hubbard bands with a
large U between them. In this case, the oxygen $2p$ band will hybridize with
the lower Hubbard band, thus the doped holes would accommodate in this Cu-3d
and O-$2p$ hybridized band. The Zhang-Rice singlet model postulates that each
doped hole on the oxygen site will couple to the $Cu^{2+}$ cation and form a
spin singlet state, it is this composed singlet that moves, leading to
electric conduction. Meanwhile there are many other models to interpret this
interplay of charge and spin in the doped cuprate, for example the t-J model
incorporating a SU(2) symmetryWenLee , and the phase string modelWengZY , etc.
All these different models can also interpret the finite coherent DOS at the
Fermi energy and a superconducting state with a diluted superfluid density
arising from a strongly correlated background. A typical electronic phase
diagram is shown in Fig.3. More interestingly, this insulating feature will
leave the way to a strange metallic behavior if we dope it with holes (hole
doping) or electrons (electron doping). Above a certain threshold of doping
($p\approx 0.05$ for example in $La_{2-x}Sr_{x}CuO_{4}$) the superconductivity
emerges at a finite temperature. The $T_{c}$ versus doping normally exhibits a
dome like shape with the optimal $T_{c}$ around a doping level of about
$p=0.16$ in the hole doped side. In the overdoped side the $T_{c}$ drops down
again and superconductivity vanishes at about $p\approx 0.25-0.30$. Phase
separation picture was proposed to explain the dropping down of the superfluid
density and the superconducting transition temperature in the overdoped
regionWenHHOverDope . Recently, mutual inductance measurements on the
systematically doped thin films illustrate also a rough linear correlation
between the superfluid density and $T_{c}$ in the overdoped region,
manifesting also an non-BCS mechanismBozovicNature .

The normal state of the cuprate superconductors is also drastically deviating
from the very base for building up the BCS paradigm, the Landau-Fermi liquid
picture. One of the exotic features here is the appearance of the so-called
pseudogap far above $T_{c}$, as shown in Fig.3 by the solid and dashed lines.
According to the different definitions, the pseudogap temperature $T^{*}$ is
also quite different. For example, the Knight shift in the NMR measurements
shows a decrease at a quite high temperature. This has been regarded as the
opening of a spin gapAlloul . The pseudogap can also be detected in the
temperature dependence of resistivity in many cuprate systems. Below the
pseudogap temperature, the resistivity shows a linear behavior with a kink at
the pseudogap temperature. Similar behavior occurs in the specific heat data
which exhibits also a decrease of the electronic specific heat coefficient
below the pseudogap temperatureLoram . While the Nernst effect measurements
reveal a strong signal which may correspond to the existence of the
simultaneous vortex-antivortex excitationsXuZA . The more direct evidence for
the pseudogap is coming from the angle resolved photoemission spectroscopy
(ARPES)ShenZX . One appealing picture for the normal state is the gapped
feature near the anti-nodal area due to the correlation effect, while the
electrons near the nodal point have still finite life-time recovering a
certain weight of quasiparticles. This leads to the so-called truncated “Fermi
surface” or Fermi arc. Although to reckon on the number of the conduction
electrons through the well known Luttinger theorem is very questionable, while
the superconducting condensation in the cuprates seem to be achieved by
building up a gap on the so-called Fermi arc near the nodal point. When
temperature is increased above $T_{c}$, the gap on the Fermi arc vanishes and
the “Fermi arc metal” is recovered as the normal stateWenLee . Connected with
this pseudogap feature, there are many other “intertwined orders/phases”
existing in this region, such as the charge density way (CDW) order, pair-
density wave (PDW) order, spin-glass, etc. A recent overview gives a nice
summary about this complexityKeimerNature . Thus we see many very exotic
features of the normal state of the cuprate superconductors. These abnormal
features in the normal state of cuprate superconductors are certainly shaking
the very base for the BCS picture. It is still under a hot debate how these
facts can be put together to reach a reconciled picture. Although some people
suggest that the BCS picture may be still applicable to interpret the
superconductivity in the cuprate system, while significant modifications have
to be undertaken, especially the focus is about whether the retarded weak
coupling picture, as the case for phonon mediated pairing, is still applicable
here.

The discovery of iron based superconductorsHosono  seems to help unraveling
the puzzle of unconventional superconductivity. Up to date, many FeAs or FeSe
based superconductors have been synthesized, which enriches the families of
the iron based superconducting systemsWenHHReview . A common knowledge now
about the iron based superconductors is that the major electric conduction and
the superconductivity are fulfilled by the FeAs or FeSe based planes.
Theoretical calculations indicate that the Fe has a cationic state of
$Fe^{2+}$ with 6 electrons in the five 3d orbitals, namely $d_{xy}$, $d_{xz}$,
$d_{yz}$, $d_{x^{2}-y^{2}}$ and $d_{3z^{2}-r^{2}}$. One can estimate the
crystalline field splitting of these five 3d orbitals in the system, which
gives a quite weak effect. Therefore in principle, all the five 3d orbitals
cross the Fermi energy leading to several sets of Fermi surfaces, two hole-
like cylindrical Fermi pockets and a small 3D like electron pocket in the
Brillouin zone center ($\Gamma$ point), and two more electron Fermi pockets in
the zone corner ($M$ point). For most FeAs based superconductors, the Fermi
surface contains a set of hole pockets centering at $\Gamma$ and a set of
electron pockets near M point, as sketched in Fig.4(b). Early but swift
theoretical calculations indicate that the phonon spectrum and the electron-
phonon interaction does not support to give superconductivity above 1 KBoeri .
Interestingly the parent phase, for example, RFeAsO (R= rare earth elements)
and BaFe2As2 all show an antiferromagnetic order at temperatures between
100-200 K. The wave vector of these AF phase measures exactly the momentum
displacement between the electron pockets located at M($\pm\pi,0$) or
($0,\pm\pi$) to $\Gamma$(0,0) in most families of iron based superconductors.
There are some preliminary evidence to show that the iron based
superconductors also have the pseudogap effectPseudogap ; Shimojima , although
it is quite rare compared with that in the cuprate superconductors. Meanwhile
the normal state, especially those in phosphorous doped systems (LaFePO, or
BaFe2As2-xPx) indeed indicates a super-linear temperature dependence of
resistivity near the possible quantum critical point where the AF
vanishesMatsudaQCP . Therefore we have the strong reason to believe that the
superconductivity in the iron based superconductors should be closely related
to the antiferromagnetic phaseImai , and may have the same origin as the
cuprates.

Recently, superconductivity up to 15 K was found in hole doped
$Nd_{1-x}Sr_{x}NiO_{2}$Ni112SC . Interests arose immediately after this
discovery since this infinite layer system may also have a $3d^{9}$ orbital as
the dominant one for the electric conduction. The discovery of
superconductivity has also been extended to other compounds, for example
$Pr_{1-x}Sr_{x}NiO_{2}$ and $La_{1-x}Ca_{x}NiO_{2}$. The calculations based on
DFT reveal that for the parent phase $Nd(La)NiO_{2}$, the major bands crossing
Fermi energy are $Ni-3d_{x^{2}-y^{2}}$, $Nd-5d_{3z^{2}-r^{2}}$, and
$Nd-5d_{xy}$ with some hybridizations. Thus it is clear that the major
contribution of superfluid density comes from the orbital of
$Ni-3d_{x^{2}-y^{2}}$, but the $Nd-5d$ bands also contribute partially to the
DOS at the Fermi levelNorman ; Thomalle . A common understanding is that the
$Ni-3d_{x^{2}-y^{2}}$ band constructs a large hole pocket around
$\Gamma(0,0,0)$ at the $k_{z}=0$, which changes into an electron pocket at the
$k_{z}$ = $\pi$ due to the dispersion along $k_{z}$ direction. The hole pocket
at $k_{z}$ = 0 looks like that of underdoped cuprate; while on the cutting
plane at $k_{z}$ = $\pi$, the electron pocket shrinks and becomes similar to
the Fermi pocket in overdoped cuprate. There are two more 3D electron pockets,
one is surrounding $\Gamma$ which is contributed by the hybridization of the
$Nd-5d_{3z^{2}-r^{2}}$ and $Ni-3d_{3z^{2}-r^{2}}$, another one surrounding
A($\pi,\pi,\pi$) is constructed by the hybridization of $Nd-5d_{xy}$ and
$Ni-3d_{3z^{2}-r^{2}}$. A spin excitation spectrum has been measured in the
undoped system $NdNiO_{2}$, which indicates a wide band spin excitations up to
an energy of 200 meVLeeWS . Based on the scheme of repulsive
interactionsThomalle , $d$-wave pairing as the dominant one has been predicted
on the central hole pocket. By measuring scanning tunneling microscopy, two
types of gaps are observedWenNi112 , one has a $d$-wave gap function with gap
maximum about 3.9 meV, another one exhibits a full gap about 2.35 meV. Some
time a mixture of the two types of gap features appears on one single
spectrum, which indicates that the $Nd_{1-x}Sr_{x}NiO_{2}$ thin film is a
multigap system and probably with a $d$-wave pairing tendency as the dominant
one. Due to the difficulties of synthesizing good samples, more results are
desired to unravel the mystery of superconductivity in the nickelate systems.

Figure 4: Cartoon picture for the pair-scattering process based on a retarded
weak coupling in (a) cuprates and (b) iron based superconductors. The open
arrows represent two electrons with opposite momenta and spins before the
pair-scattering; the filled arrows show the final states after the pair-
scattering. The dashed lines show the pairing channel in the pair-scattering
process through a retarded interaction. If the interaction V is positive, this
kind of pairing will induce sign change of the superconducting gaps in the
momentum space, for example d-wave in cuprates and s±-wave in iron-based
superconductors.

##  III Preliminary understanding of the unconventional superconductivity

If there exist strong Coulomb interactions $U$ among electrons, or the ratio
between the repulsive coulomb interaction and the band width $t$, i.e., $U/t$
is large, the free electron gas or the so-called Landau-Fermi liquid picture
immediately becomes questionable. In this case one cannot assume a negligible
potential energy in the electronic Hamiltonian, the latter has been adopted in
building up the BCS picture. With the presence of strong correlations, the
Hamiltonian has been considered based on the $t-U$, or $t-J$ or $t-J-U$
models. Since tremendous number of electrons need to be considered in the
system, this multibody problem is not solvable. However, this problem may be
at least inspected by assuming a low energy approach. In this case, a
microscopic Hubbard-type Hamiltonian with a four fermion interaction can still
be written asChubukov

$H=\sum_{k,\alpha}\varepsilon_{k}\psi^{+}_{k\alpha}\psi_{k,\alpha}+\sum_{k_{i},\alpha_{i}}U^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}}_{k_{1},k_{2},k_{3},k_{4}}\psi^{+}_{k_{1},\alpha_{1}}\psi^{+}_{k_{2},\alpha_{2}}\psi_{k_{3},\alpha_{3}}\psi_{k_{4},\alpha_{4}}.$
(2)

Here
$U^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}}_{k_{1},k_{2},k_{3},k_{4}}$
represents the four electron interactions, $\psi^{+}_{k_{i},\alpha_{i}}$
($\psi_{k_{i},\alpha_{i}}$) is the creation (annihilation) operator for
electrons with the spin $\alpha_{i}$ and momentum $k_{i}$ (i=1,2,3,4). For a
single band Hubbard model with local Coulomb interactions, we have

$U^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}}_{k_{1},k_{2},k_{3},k_{4}}=U\delta_{k_{1}+k_{2}-k_{3}-k_{4}}(\delta_{\alpha_{1}\alpha_{4}}\delta_{\alpha_{2}\alpha_{3}}-\delta_{\alpha_{1}\alpha_{3}}\delta_{\alpha_{2}\alpha_{4}}).$
(3)

If only the interactions in the low energy area are considered, the situation
can be treated by using the so-called random phase approximation (RPA)
methodScalapinoRPA . For a system with very good Fermi surface nesting, the
electronic system will be instable to some long range ordered state, such as
the charge-density-wave (CDW) or the spin-density-wave (SDW) state. In the
iron based systems, the perfect nesting between the hole and electron pockets
may just satisfy this situation leading to the formation of a long range SDW
order in the parent phase in the FeAs based families. There are some arguments
that the AF order is induced by the local super-exchange effect. By chemical
doping, or applying a pressure, the perfect nesting situation deteriorates,
and the long range ordered SDW state will be replaced by the state with strong
antiferromagnetic spin fluctuations (AF-SF). The elementary excitations of
this state are the so-called paramagnons. When an itinerant electron is
passing through a position, it will polarize the spin cloud around it, this
polarized spin cloud will be sensed by the second itinerant electron passing
by. This process will build up a pairing interaction between the two electrons
with a pairing interaction as

$\Gamma_{s}(k,k^{\prime})=\frac{3}{2}U^{2}\frac{\chi_{0}(q)}{1-U\chi_{0}(q)},$
(4)

where $U$ is the repulsive energy, $\chi_{0}$ is the bare spin susceptibility
in absence of interactions, the subscript $s$ here represents the singlet
pairing channel. One can see that $\Gamma_{s}$ is also peaked at the wave
vector $q$ but is positive. If we use this pairing interaction and follow the
BCS approach, the gap function can be written as

$\Delta_{k}=-\sum_{k^{\prime}}\Gamma_{s}(k,k^{\prime})\frac{\Delta_{k}^{\prime}}{2E_{k}^{\prime}}tanh\frac{E_{k}^{\prime}}{2T}$
(5)

Here
$E_{k^{\prime}}=\pm\sqrt{\varepsilon_{k^{\prime}}^{2}+\Delta_{k^{\prime}}^{2}}$.
In iron based superconductors, this interesting pairing picture was first
proposed by Mazin et al.Mazin , and later further formalized by several other
groups using different theoretical approachesKuroki ; WangZD ; Hirschfeld1 ;
LeeDH . More details about the pairing order parameter can be found in a
recent reviewHirschfeld2 . The basic idea of weak coupling picture is
illustrated in Fig.4(a) and (b) for cuprates and iron based superconductors,
respectively. It still remains unclear whether this weak coupling based
picture applies to the cuprates, since the pairing is very strong and Fermi
surface is incomplete (truncated). In iron based superconductors, especially
in FeAs based superconductors, the two electrons on the electron (hole) FSs
(marked with the thick arrows indicating opposite directions of spins) are
scattered to the hole (electron) FSs by exchanging the antiferromagnetic spin
fluctuations (AF-SFs). If the pairing is established through exchanging the
AF-SFs, according to above equations, the so-called $S^{\pm}$ pairing model is
expected, i.e., the DOS on individual FS are fully gapped (although the gap
magnitude has some momentum dependent variation), but the signs of the
superconducting gaps on the electron and hole pockets should be opposite. This
kind of s-wave gap has been supported by the scanning tunneling experiment
alreadyHanaguri . The recent experiment in observing the in-gap quasi-particle
states by non-magnetic Cu impurities in $NaFe_{0.97-x}Co_{0.03}Cu_{x}As$ may
provide a more decisive evidence for the $S^{\pm}$ pairingCuImpurity .
Furthermore, this sign reversal of gaps on the electron and hole pockets in
$FeTe_{1-x}Se_{x}$ has been well checked by the phase referenced quasiparticle
interference techniqueChenMYPRB . This RPA based picture may also be extended
to the case of cuprates for understanding the well known $d$-wave pairing
symmetry. As shown in Fig.4(a), since there is only one set of Fermi surface
in the cuprate, the pair-scattering would lead to a sign change on the
different segments of the Fermi surface, this naturally leads to a $d$-wave
gap function.

Figure 5:  A roughly linear relation between the neutron resonance or bosonic
mode energy and the superconducting transition temperature $T_{c}$ in many
iron based superconductors. The small sized filled symbols are data collected
in the inelastic neutron scattering measurements, the two red bigger square
and circle symbols are determined from the bosonic mode on the tunneling
spectrum. This figure is made by quoting to the data in Ref.[52]

One of the consequences of the electronic system with the pairing gap with
opposite signs is the appearance of a resonance peak in the imaginary part of
the spin susceptibility at the wave vector of the SDW. This resonance peak
with a typical energy $\Omega\leq 2\Delta_{s}$ can be detected with inelastic
neutron scattering experiment and has been observed in cupratesRPCuprate ,
heavy FermionRPHeavyFermion  and recently in the iron pnictide
superconductorsRPIron . This resonance peak may be understood as due to the
spin flipping scattering when exciting a quasiparticle from the $k$ to $k+q$
with the pairing gap with the opposite signs, since in this case the BCS
coherence factor is given byScalapinoRPA ; Hirschfeld2

$I_{res}\propto\frac{1}{2}[1-\frac{\Delta(k)\Delta(k+q)}{E(k)E(k+q)}].$ (6)

Here $I_{res}$ gives the intensity of the resonance peak, $E(k)$ is the
dispersion of the Bogoliubov quasiparticles in the superconducting state
$E(k)=\pm\sqrt{\varepsilon_{k}^{2}+\Delta^{2}(k)}$. One can see that, if the
superconducting gap has a sign reversal at momenta connected by vector $q$,
that will produce a peak at this particular vector. The RPA calculations tell
that this peak should locate at a energy of $\Omega\leq 2\Delta_{s}$ depending
on the correlation strength. Although the resonance peak is regarded as the
consequence of the sign-reversal gaps, sometimes it is called as the “egg” of
the superconducting pairing state, not the pairing glue for superconductivity,
we must emphasize that this “egg” can be naturally traced back to the scenario
of the magnetic pairing, namely the origin of pairing. Therefore, although the
resonance peak may be due to the spin-1 exciton mode of the sign reversal
gaps, we can conclude that it is only possible that the pairing is induced
through the electronic origin, either through exchanging the AF-SFs, or the
superexchange itself, but certainly not through some non-magnetic origins,
like exchanging phonons. In Fig.5, we present a collection of many resonant
energies versus the superconducting transition temperature $T_{c}$ in
different iron based superconductors. Here the filled symbols are data
collected in the inelastic neutron scattering measurements, the two open
symbols are determined from the bosonic mode on the tunneling spectrum. In
some iron based superconductors, the scanning tunneling spectroscopy
measurements reveal extra peaks (humps) outside the superconducting coherence
peaks, which is explained as the coupling of the quasiparticles with the spin
resonance mode. In Fig.5, a simple ratio $\Omega\approx(4.3\pm 0.5)k_{B}T_{c}$
is shownNatPhysWangZY . This strongly suggests that the neutron resonance is
indeed closely related to superconductivity. Similar linear relation between
$\Omega$ and $T_{c}$ was also found in cuprate superconductorsPNAS . A study
on the superconducting condensation energy measured through specific heat and
the intensity of the resonance peak in the
$Pr_{0.88}LaCe_{0.12}CuO_{4-\delta}$ finds that both quantities decay with the
applied magnetic field in the same wayPNAS , indicating a close relationship
between the resonance mode energy and superconductivity.

Figure 6: Evolution of the pairing from weak coupling to strong coupling,
depending on the ratio $U/t$. The left hand picture shows the case for
$U/t<<1$, the BCS scheme works well. Here the four blue circles represent the
distortion of the atoms due to the electron-phonon interaction. This
distortion will be “seen” by the second electron passing by in a retarded
duration. The right hand one represents the possible situation of strong
correlation effect where the local super-exchange may dominate the pairing. In
the middle, the $U/t\approx 1$ which may be applicable for the case of AF-SF
mediated pairing. Figure 7: Schematic camp for superconductivity. The BCS
theory with weak coupling occupies a small region in the left hand side. We
need a general or global camp to cover all possibilities for
superconductivity, from weak coupling to strong local pairing.

In above discussion, we have shown that the low energy approach, or called as
the picture of weak coupling limit can give a reasonably good understanding of
the superconductivity in many unconventional superconducting systems. In the
systems with moderate Coulomb repulsion energy $U$, this seems to work quite
well. In some systems with extremely large ratio $U/t$, it is questioned
whether this picture based on retarded scattering can still be valid. In this
case, some people argue that the local super-exchange $J$ dominates the
pairing process, leading to the simultaneous pairingLocalPair1 ; LocalPair2
with a very strong pairing strength. One of the extreme cases is the so-called
resonance-valence-bond (RVB) spin liquid state for the $s=1/2$ systems, for
example in the cuprates. In this case, the spin singlet pairs are naturally
formed yielding a highly fluctuating quantum RVB state.RVB  This is certainly
a very interesting picture which may get an indirect support from the strong
gap, very strong superconducting fluctuation, or residual Copper pairs in the
normal state far above $T_{c}$. The convincing evidence for this RVB state is
still lacking. In Fig.6 we present a schematic show for the superconducting
pairing for the cases ranging from the weak coupling (e.g., phonon mediated)
pairing to the extreme case of the simultaneous pairing through the local
super-exchange. The situation evolves with the ratio of $t/U$: (1) $t/U\gg 1$,
weak coupling induced pairing; (2) $t/U\approx 1$, intermediate correlation
induced pairing; (3) $t/U\ll 1$ strong correlation and local pairing could
occur. This picture follows a similar thought as that when considering the
subtle balance between the correlation effect and the effective Drude weight
for optimized superconductivity in a variety of compoundsBasov . In the middle
of Fig.6, we have a delicate balance between the two extreme cases: complete
itinerant ($U/t<<1$) and strong correlation ($U/t>>1$). The middle case may
apply to the iron based superconductors, which still needs the retarded
interaction as the pairing origin, but through exchanging the paramagnon or
AF-SFs, not phonons. At the moment or in the near future, debate will go on
about whether we need to have a more exotic picture which assumes the pairing
completely coming from the local strong superexchange. Based on this idea, in
Fig.7, we sketch a big camp for superconductivity. The BCS theory occupies
only a small area in the weak coupling limit. It remains to know how the
pairing can evolve smoothly from the weak coupling limit to the extremely
strong coupling limit. It is expectable that the condensation process for
superconductivity will also evolve from the BCS type to BEC type, the
superconducting transition temperature is determined by the gap in the BCS
case, but by the phase stiffness in the BEC caseKievelson ; Randeria .

##  IV Some empirical rules for the exploration of high temperature
superconductors

It is always enthusiastic to explore new high temperature superconductors.
However, beside those of pressurized hydrides, no theory can really predict
how to get success, thus it has been a quite frustrated task. Even so, some
empirical rules may be summarized here. Here we would like to propose several
key points which may be the necessary ingredients for exploring high
temperature superconductivity at ambient pressure.

###  IV.1 Subtle balance between Mottness and itinerancy of $d$ orbital
electrons

As illustrated by Qazilbash et al.Basov , the unconventional superconductors
should have a delicate balance between Mottness and itineracy which governs
the effective DOS at the Fermi energy. Strong Mottness may manifest a large
superexchange energy $J$ for the local pairing, but the ground state of that
may be an insulator. Thus we need to effectively separate the dual roles of
the $d$-orbital electrons. The 3$d$ orbital electrons possess by them-self a
possibility to balance between these two factors. Thus it may give more
chances in the compounds with the 3$d$ transition metal elements, like Cr, Mn,
Fe, Co, Ni, Cu. Some 4$d$ elements may also be considered, such as Ru, Pd, Rh
etc.

###  IV.2 Bad metal in vicinity of an AF order

Recalling both the iron based and cuprate superconductors, it is known that
the parent phases are either Mott insulators or bad metals with or near an AF
order. Most importantly, this AF order can be easily suppressed by chemical
doping or applying a pressure. To satisfy this condition, it is important to
have ions with lower magnetic moments, for example below 2$\mu_{B}$. The spin
density wave order due to the Fermi surface instability would help. To avoid
the large magnetic moments, one may request the low filling or even number
filling states of the 3$d$ or 4$d$ orbitals, such as the Cu-3$d^{9}$ or
Fe-3$d^{6}$ cationic states. This is actually the case in iron based
superconductors. In BaFe2As2, for example, the cationic state is Fe2+, which
yields a magnetic moment generally smaller than 2$\mu_{B}$. In the cuprate,
for example, each Cu ion in $La_{2}CuO_{4}$ has a magnetic moment of 0.5-1.0
$\mu_{B}$. The strong quantum fluctuation of the gauge field near the quantum
critical point (QCP) may also help. An odd filling number of 3, 5, 7 of the
$d$ orbital will induce a strong Hund’s coupling leading to strong magnetic
moments and localization of electrons.

###  IV.3 Shallow and flat band near the Fermi energy

Some shallow and flat bands, or saddle points near the Fermi energy will
result in the von Hove singularity of DOS, this is certainly helpful for the
pairing instability. In iron based superconductors, it is known that some
bands are very shallow, which yields a high DOS. In addition, a shallow band
possesses a small Fermi energy $E_{F}$, which leads to the small ratio
$E_{F}/\Delta_{s}$. This quantity measures the Cooper pair numbers within a
coherence length $\xi$ or volume $\xi^{3}$. If this value is small, one may
expect a BEC type superconducting transition. Therefore, to have both strong
pairing and high superfluid density, we probably need a system with both
shallow and wide bands and together with a strong interband scattering, the
shallow band produces a strong local pairing, the wide band contributes a high
superfluid density.

###  IV.4 Anion mediated superexchange may help

In cuprates and iron based superconductors, there are anions to mediate the
superexchange for the antiferromagnetic order or spin fluctuations. Thus in
the dominant conductive planes, some atomic chains like T-A-T are necessary,
where $T$ represents a transition metal with $d$ orbital electrons, $A$ is an
anion. For example, in cuprates the pairing through $J_{1}$ superexchange is
established by the Cu-O-Cu bonds. While in iron based superconductors, the
pairing may be established by the $J_{2}$ superexchange which is generated by
the Fe-As-Fe bonds. However, this requirement may not be that strict for high
$T_{c}$ superconductivity. In many binary alloy compounds, unconventional
superconductivity coexists with the magnetic excitations. Thus for
unconventional superconductivity, the anion bridge may not be the necessary
component, but may be helpful for stabilizing the structure.

## Acknowledgement

We acknowledge the useful discussions with Qianghua Wang, Ilya Eremin, Peter
Hirschfeld, Andrey Chubukov, Igor Mazin. This work was supported by NSF of
China with Grant Nos. 11927809, NSFC-DFG1206113001(ER463/14-1), and the
Ministry of Science and Technology of China.

## References

## References

  * (1) Bardeen J, Cooper L N and Schrieffer J R 1957 Phys. Rev. 108 1175
  * (2) Wen, H.-H 2014 Japn Phys. Soc. Conf. Proc.1 011005
  * (3) Nagaatsu J, Nakagawa N, Muranaka T et al. 2001 Nature 410 63-65
  * (4) Drozdov A P, Eremets M I, Troyan I A, Ksenofontov V and Shylin S I 2015 Nature 525 73¨C76
  * (5) Steglich F, Aarts J, Bredl C D, Lieke W, Meschede D, Franz W and H. Schäfer 1979 Phys. Rev. Lett. 43 1892
  * (6) Jerome D, Mazaud A, Ribault M and Bechgaard K 1980 J. Phys. Lett. 41 L95. For a recent review on organic superconductors, see Ardavan A et al 2012 J. Phys. Soc. Jpn. 81 011004
  * (7) Bednorz J G and Müller K A 1986 Z. Phys. B64 189
  * (8) Wu M K, Ashburn J R, Tirng C J, Meng R L, Gao L, Huang Z L, Wang Y Q and Chu C W 1987 Phys. Rev. Lett. 58 908
  * (9) Zhao Z X et al 1987 KeXue TongBao 32 412
  * (10) Maeda H, Tanaka Y, Fukutomi M and Asano T 1988 Jpn. J. Appl. Phys. 27 L209
  * (11) Sheng Z Z, Herman A M, Ali A, Almasan C, Estrada J, Datta T and Matsonet R J 1988 Phys. Rev. Lett. 60 937
  * (12) Putilin S N, Antipov E V, Chmaissem O and Marezio M 1993 Nature 362 226
  * (13) Gao L, Xue Y Y, Chen F, Xiong Q, Meng R L, Ramirez D, Chu C W, Eggert J H and Mao H K 1994 Phys. Rev. B50 4260(R)
  * (14) Anderson P W 1987 Science 235 1196
  * (15) Zhang F C and Rice T M 1988 Phys. Rev. B 37 3759
  * (16) Wen X G and Lee P A 1996 Phys. Rev. Lett. 76 50; Lee P A and Wen X G 1997 Phys. Rev. Lett. 78 4111
  * (17) Weng Z Y, Sheng D N, Chen Y C and Ting C S 1997 Phys. Rev. B 55 3894; Sheng D N, Chen Y C and Weng Z Y 1996 Phys. Rev. Lett. 77 5102
  * (18) Wen H H, Yang W L, Zhao Z X and Ni Y M 1999 Phys. Rev. Lett. 82 410
  * (19) Bozovic I, He X, Wu J and Bollinger A T 2016 Nature 536 309¨C311
  * (20) Newns D M and Tsuei C C 2007 Nature Phys. 3 184
  * (21) Alloul H, Ohno T and Mendels D 1989 Phys. Rev. Lett. 63 1700
  * (22) Loram J W, Mirza K A, Cooper J R and Liang W Y 1993 Phys. Rev. Lett. 71 1740
  * (23) Xu Z A, Ong N P, Wang Y, Kakeshita T and Uchida S 2000 Nature 406 486
  * (24) For a review of the ARPES study on the Fermi surface see Damascelli A, Hussain Z and Shen Z X 2003 Rev. Mod. Phys. 75 473
  * (25) Keimer B, Kivelson S A, Norman M R, Uchida S and Zaanen J 2015 Nature 518 179
  * (26) Kamihara Y, Watanabe T, Hirano M and Hosono H 2008 J. Am. Chem. Soc. 130 3296
  * (27) Wen H H and Li S L 2011 Annu. Rev. Condens. Matter Phys. 2 121
  * (28) Boeri L, Dolgov O V and Golubov A A 2008 Phy. Rev. Lett. 101 026403
  * (29) Dai Y M, Xu B, Shen B, Wen H H, Hu J P, Qiu X G and Lobo R P S M 2012 Phys. Rev. B. 86 100501(R)
  * (30) Shimojima T, Sonobe T, Malaeb W et al. et al 2014 Phys. Rev. B 89 045101
  * (31) Hashimoto K et al 2012 Science 336 1554
  * (32) Ning F L, Ahilan K, Imai T, Sefat A S, McGuire M A, Sales B C, Mandrus D, Cheng P, Shen B and Wen H H 2010 Phys. Rev. Lett. 104 037001
  * (33) Li D F, Lee K, Wang B Y, Osada M, Crossley S, Lee H R, Cui Y, Hikita Y and Hwang H Y 2019 Nature 572 624
  * (34) Botana A S and Norman M R 2020 Phys. Rev. X 10 011024
  * (35) Wu X X, Sante D D, Schwemmer T, Hanke W, Hwang H Y, Raghu S and Thomale R 2020 Phys. Rev. B101 060504(R)
  * (36) Lu H et al 2021 Science 373 213
  * (37) Gu Q Q, Li Y Y, Wan S Y, Li H Z, Guo W, Yang H, Li Q, Zhu X Y, Pan X Q, Nie Y F and Wen H H 2020 Nat. Commun. 11 6027
  * (38) Chubukov A V 2012 Annu. Rev. Condens. Matter Phys. 3 57
  * (39) Scalapino D J 2012 Rev. Mod. Phys. 84 1383
  * (40) Mazin I I, Singh D J, Johannes M D and Du M H 2008 Phys. Rev. Lett. 101 057003
  * (41) Kuroki K, Onari S and Arita R 2008 Phys. Rev. Lett. 101 087004
  * (42) Yao Z J, Li J X and Wang Z D 2009 New J. Phys. 11 025009
  * (43) Maier T, Graser S, Hirschfeld P J and Scalapino D J 2009 Phys. Rev. B79 224510
  * (44) Wang F and Lee D H 2009 Phys. Rev. Lett. 102 047005
  * (45) Hirschfeld P J, Korshunov M M and Mazin I I 2011 Rep. Prog. Phys. 74 124508
  * (46) Hanaguri T, Niitaka S, Kuroki K and Takagi K 2010 Science 328 474
  * (47) Yang H, Wang Z Y, Fang D L, Deng Q, Wang Q H, Xiang Y Y, Yang Y and Wen H H 2013 Nat. Commun. 4 2749
  * (48) Chen M Y, Tang Q K, Chen X Y, Gu Q Q, Yang H, Du Z Y, Zhu X Y, Wang E Y, Wang Q H and Wen H H 2019 Phys. Rev. B99 014507
  * (49) Rossat-Mignod J, Regnault L P, Vittier C, Bourges P, Burlet P and Bossy J 1991 Physica C185-189 86
  * (50) Stock C, Broholm C, Hudis J, Kang H J and Petrovic C 2008 Phys. Rev. Lett. 100 087001
  * (51) Christianson A D et al 2008 Nature 456 930
  * (52) Wang Z Y, Yang H, Fang D L, Shen B, Wang Q H, Shan L, Zhang C L, Dai P C and Wen H H 2013 Nature Phys. 9 42
  * (53) Wilson S D, Li S L, Zhao J, Mu G, Wen H H, Lynn J W, Freeman P G, Regnault L P, Habicht K and Dai P C 2007 Proc. Natl. Acad. Sci. U.S.A. 104 15259
  * (54) Si Q M and Abrahams E 2008 Phys. Rev. Lett. 101 076401
  * (55) Seo K, Bernevig B A and Hu J P 2008 Phys. Rev. Lett. 101 206404
  * (56) Anderson P W, Lee P A, Randeria M, Rice T M, Trivedi N and Zhang F C 2004 J. Phys. Cond. Matter 16 R755
  * (57) Qazilbash M M, Hamlin J J, Baumbach R E, Zhang L, Singh D J, Maple M B and Basov D N 2009 Nature Phys. 5 647
  * (58) Emery V J and Kivelson S A 1995 Nature 374 434
  * (59) Randeria M and Taylor E 2014 Annu. Rev. Condens. Matter Phys. 5 209


# Integrated guidance and control framework for the waypoint navigation of a
miniature aircraft with highly coupled longitudinal and lateral dynamics

K Harikumar, Jinraj V. Pushpangathan, Sidhant Dhall, and M. Seetharama Bhat K.
Harikumar is a Research Fellow @ Singapore University of Technology and
Design, Singapore, e-mail<EMAIL_ADDRESS>Jinraj V. Pushpangathan is
a Research Fellow @ Department of Aerospace Engineering, Indian Institute of
Science, Bangalore, India, e-mail<EMAIL_ADDRESS>Sidhant Dhall is a
Research Assistant @ Department of Aerospace Engineering, Indian Institute of
Science, Bangalore, India, email<EMAIL_ADDRESS>M. Seetharama Bhat
is a Professor @ Department of Aerospace Engineering, Indian Institute of
Science, Bangalore, India, e-mail<EMAIL_ADDRESS>
###### Abstract

A solution to the waypoint navigation problem for fixed wing micro air
vehicles (MAV) is addressed in this paper, in the framework of integrated
guidance and control (IGC). IGC yields a single step solution to the waypoint
navigation problem, unlike conventional multiple loop design. The pure
proportional navigation (PPN) guidance law is integrated with the MAV
dynamics. A multivariable static output feedback (SOF) controller is designed
for the linear state space model formulated in IGC framework. The waypoint
navigation algorithm handles the minimum turn radius constraint of the MAV.
The algorithm also evaluates the feasibility of reaching a waypoint. Extensive
non-linear simulations are performed on high fidelity 150 mm wingspan MAV
model to demonstrate the potential advantages of the proposed waypoint
navigation algorithm.

###### Index Terms:

Integrated guidance and control, coupled dynamics, micro air vehicle, static
output feedback, waypoint navigation.

##  I Introduction

Fixed wing micro air vehicles (MAV) are extremely agile, lightweight aircraft
with a maximum wingspan of 150 mm [1]. In any of the MAV applications, it has
to reach a set of waypoints, specified by the end user. In most of the
waypoint navigation problems, guidance and control are treated as two separate
problems. Control system is usually designed to give higher bandwidth to track
the commands generated by the guidance law. But when the guidance law is
integrated with the inner loop control system, the combined system often fails
to meet the required performance and stability specifications. The major
reasons are due to the assumption of lower order dynamics for the inner loop
control system, neglecting actuator bandwidth and assuming decoupled
longitudinal and lateral dynamics, especially for the case of MAVs. Many of
the guidance law in the literature assumes lower order autopilot for inner
loop control system [2]. The assumption of lower order inner loop autopilot
provides fictitious high amplitude gain margin for guidance law to operate.
When the combined system is not stable, the guidance law parameters are
retuned to ensure stability. The re-tuning of guidance law parameters often
lead to inadequate performance. The actuator bandwidth is of major concern for
MAVs. The short period mode and Dutch roll mode natural frequencies of MAVs
are comparable to actuator bandwidth [3]. The relatively lower actuator
bandwidth in MAVs leads to poor wind disturbance rejection. The wind
disturbances change angle of attack $(\alpha)$ and side slip angle $(\beta)$
considerably from the trim operating conditions. The additional forces and
moments generated by wind disturbances are often neglected while designing the
guidance law in the presence of wind [4]. The coupling between longitudinal
and lateral dynamics is severe for MAVs and cannot be neglected as in the case
of bigger UAVs [3],[5]. The guidance law design needs to consider the combined
longitudinal and lateral dynamics, failure of which leads to instability
during abrupt maneuvers.
The topic of integrated guidance and control (IGC) has gained considerable
interest in the area of missile guidance and control [6]-[11] due to the
following advantages. The synergy between the guidance and control subsystems
can be exploited to make the combined system more stable and effective than
the individual subsystems [6]. Despite the fact that separate guidance and
control has been applied successfully, still there remains scope for further
optimizing the performance and enhancing the stability [7]. The prime focus on
the missile IGC design is to achieve improved terminal miss distances. When
the missile approaches the target, the time scale separation between the
control loop and guidance loop is no longer valid [8]. The conventional
approach of designing separate guidance and control leads to a conservative
design with the requirement of actuators with higher bandwidth than necessary
[9].
The idea of IGC is also applied to fixed wing UAVs [12]-[16]. A linear time
invariant model in the framework of integrated guidance and control is
developed for a fixed wing UAV of wingspan 3.8 $m$ to track helical trim
trajectories in [12]. For a UAV flying along the trim trajectory, the net
forces and moments acting on the UAV are zero. The method developed in [12]
explains how to follow a given helical trajectory in which UAV remains in a
trimmed state. But finding such a trimmed trajectory connecting two waypoints
is very computationally intensive. Moreover, the method cannot be applied for
in-flight waypoint changes since it is difficult to compute the trim
trajectories online with the limited onboard computational power available
with MAV. The IGC design for a fixed wing UAV of 8.8 $m$ wingspan for
following a predefined path using sliding mode control is explained in [13].
The design assumes that the roll, pitch and yaw dynamics are decoupled. But
for MAVs, the roll, pitch and yaw dynamics are tightly coupled [5]. Nonlinear
dynamic inversion (NDI) based IGC framework is developed for UAVs in [14] and
[15]. NDI utilizes full state feedback structure. In the case of MAVs,
measurement of $\alpha$ and $\beta$ is tedious due to the lack of availability
of lightweight sensors. Moreover, the perturbations in $\alpha$ and $\beta$
are very high for the MAVs when compared to bigger UAVs under nominal wind
conditions. IGC framework developed for lateral dynamics of a fixed wing UAV
assumes decoupled longitudinal dynamics in [16]. But the assumptions of
decoupled longitudinal and lateral dynamics is not valid for 150 mm class of
MAVs [3],[5].
To overcome the above-mentioned limitations of existing IGC methodology, a
novel IGC framework is derived in this paper, combining the idea of pure
proportional navigation (PPN) with MAV dynamics. The proposed waypoint
navigation algorithm based on IGC framework is non-iterative and hence
computationally inexpensive. The issue of coupling between longitudinal and
lateral dynamics of the MAV is handled by the proposed IGC framework. Unlike
the existing algorithms based on IGC framework, the proposed waypoint
navigation algorithm evaluates the feasibility of reaching a waypoint based on
the acceleration constraints of the MAV. A $H_{\infty}$ SOF controller is
designed for the combined guidance and control loop using Linear Matrix
Inequality (LMI) techniques. Numerical simulation results are presented for
the straight line following and rectangular path following using high fidelity
nonlinear model of a MAV obtained from wind tunnel data.
The paper is organized as follows. The coupled linear state space model of the
MAV is presented in section II. The guidance law used in IGC framework is
given in section III. The integration of guidance law to the coupled dynamics
of MAV and the waypoint navigation algorithm is presented in section IV. The
non-linear six degrees of freedom simulations for straight line following and
rectangular path following is given in section V, followed by conclusions.

##  II Preliminaries

###  II-A Coupled linear state space model of the MAV

The specifications of the MAV considered in this paper is given in Table I,
[1], [5].

TABLE I: Specifications of the MAV Planform | Rectangular
---|---
Take-off weight | 53 grams
Cruise airspeed | 8 $m/s$
Stall airspeed | 6 $m/s$
Wing - chord length | 0.11 $m$
Wing - span length | 0.15 $m$
Airfoil | E387
Control surfaces | Elevator and rudder

The coupling between longitudinal and lateral dynamics is severe for the case
of MAVs. The coupled linear state space shows an unstable spiral mode, whereas
the decoupled model fails to capture the unstable spiral mode [5]. The effects
of asymmetric propeller wake, motor-propeller counter torque and gyroscopic
effects lead to unstable spiral mode dynamics [5], [17]. The coupled linear
state space model of the MAV is given in (1).

$\dot{X}=A_{c}X+B_{c}U$ (1)

where

$\displaystyle
X=[\widetilde{u},\widetilde{w},\widetilde{q},\widetilde{\theta},\widetilde{h},\widetilde{v},\widetilde{p},\widetilde{r},\widetilde{\phi}]^{T}$
(2)

and

$U=[\widetilde{\delta}_{es},\widetilde{\delta}_{rs},\widetilde{\delta_{t}}]^{T}$
(3)

The set ($\widetilde{u},\widetilde{v},\widetilde{w}$) represents the
linearized velocity along body axis,
($\widetilde{p},\widetilde{q},\widetilde{r}$) are the linearized angular
rates, $\widetilde{\phi}$ is the linearized roll angle, $\widetilde{\theta}$
is the linearized pitch angle and $\widetilde{h}$ is the linearized altitude.
The linearized elevator servo output $\widetilde{\delta}_{es}$, linearized
rudder servo output $\widetilde{\delta}_{rs}$ and the linearized thrust
generated by motor-propeller denoted by $\widetilde{\delta_{t}}$ are the
control inputs. $A_{c}$ is the system matrix for coupled state space model and
$B_{c}$ is the input matrix. The matrices $A_{c}$ and $B_{c}$ are given in [5]
for the case of 150 mm wingspan fixed wing MAV.

The states of longitudinal dynamics is given by

$X_{1}=[\widetilde{u},\,\,\widetilde{w},\,\,\widetilde{q},\,\,\widetilde{\theta},\,\,\widetilde{h}]^{T}$
(4)

Similarly, the states of lateral dynamics is given by

$X_{2}=[\widetilde{v},\,\widetilde{p},\,\,\widetilde{r},\,\,\widetilde{\phi}]^{T}$
(5)

The matrix $A_{c}$ in (1) can be written as

${A_{c}}=\left(\begin{array}[]{cc}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{array}\right)$ (6)

where $A_{11}$ denotes the state matrix for decoupled longitudinal dynamics
and $A_{22}$ the state matrix for decoupled lateral dynamics. The coupling
between longitudinal and lateral dynamics is introduced by the matrices
$A_{12}$ and $A_{21}$. For bigger UAVs, the eigenvalues of $A_{c}$ denoted by
$\bar{\lambda}(A_{c})$ = $\bar{\lambda}(A_{11})$ $\cup$
$\bar{\lambda}(A_{22})$. But for the case of 150 mm wingspan MAVs,
$\bar{\lambda}(A_{c})$ $\neq$ $\bar{\lambda}(A_{11})$ $\cup$
$\bar{\lambda}(A_{22})$ [5].
Remark: In the rest of this paper, a variable with $\tilde{.}$ symbol denotes
the linearized variable. For example, the variables ($u,\,,v,\,w$) denotes the
velocity components along the body axis and
($\widetilde{u},\,\widetilde{v},\,\widetilde{w}$) represents their linearized
version.

##  III Guidance logic for waypoint navigation in two dimensional space

The origin of the local inertial coordinate system ($X_{I}Y_{I}Z_{I}$) is the
home location of the MAV as shown in Fig. 1 with positive $X_{I}$ axis towards
the geographic north pole, positive $Y_{I}$ axis in the direction of $90^{o}$
clockwise rotation from positive $X_{I}$ axis and positive $Z_{I}$ axis points
towards the down direction. For the waypoint following in 2D space, altitude
hold is employed using altimeter feedback, so motion in the local inertial
$X_{I}Y_{I}$ plane is only considered. The guidance logic is explained first
in this section and the integration of the guidance logic with the control
system is explained later in section IV.
In Fig. 2, axes $X^{\prime}Y^{\prime}$ are parallel to $X_{I}Y_{I}$, $(x,y)$
is the current MAV position in the $X_{I}Y_{I}$ plane, $(x_{f},y_{f})$ is the
next waypoint to be followed, $r_{a}$ is the range, $\sigma$ is the angle
between the line joining $(x,y)$ , $(x_{f},y_{f})$ and the $Y^{\prime}$ axis,
$\chi$ is the angle made by the velocity vector with the $X^{\prime}$ axis.
The miss distance $d$ is the perpendicular distance between $(x_{f},y_{f})$
and the velocity vector. The applied acceleration is denoted by $a_{c}$ and is
perpendicular to the velocity vector and lies in the $X_{I}Y_{I}$ plane. The
units of $(x,y)$, $(x_{f},y_{f})$, $d$, $r_{a}$ are in meters (m), $V_{a}$ is
in meters/second ($m/s$), $a_{c}$ is in meters/second-square($m/s^{2}$) and
$\sigma$, $\chi$ are in radians(rad). From Fig. 2, we can write the following
relations given in (7) to (11).

Figure 1: Diagram showing inertial $X_{I}Y_{I}Z_{I}$ co-ordinate sytem Figure
2: Diagram explaining the guidance law

$\displaystyle\sigma=tan^{-1}(\frac{x_{f}-x}{y_{f}-y})$ (7)
$\displaystyle\rho=\frac{\pi}{2}-(\chi+\sigma)$ (8) $\displaystyle
r_{a}=\sqrt{(x_{f}-x)^{2}+(y_{f}-y)^{2}}$ (9) $\displaystyle d=r_{a}sin\rho$
(10) $\displaystyle a_{c}=k_{1}d$ (11)

In (11), $k_{1}$ is the proportional gain. The magnitude of applied
acceleration depends upon $k_{1}$. An expression for the gain $k_{1}$ is given
in the subsection III-A. The acceleration is applied perpendicular to the
velocity vector in such a way that it does not change the magnitude of
velocity but only the direction. This is analogues to Pure Proportional
Navigation (PPN) where the acceleration applied is perpendicular to the
velocity vector [19]. Here it is assumed that altitude and $V_{a}$ is
constant. When the miss distance is zero, from (10), we obtain either
$r_{a}=0$ or $\rho=0$. The condition $r_{a}=0$ implies that $x=x_{f}$ and
$y=y_{f}$. From Fig. 2 , $\rho=0$ implies that the velocity vector is pointing
towards the goal point. So when the miss distance is zero, the MAV has either
reached the desired waypoint or the velocity vector of the MAV is along the
line connecting the waypoint and current MAV position such that no further
acceleration has to be applied to change the direction of the MAV.

###  III-A Relation to PPN Guidance Law and choice of gain $k_{1}$

This subsection gives the relation between the proposed guidance law given in
(11) and the Pure Proportional Navigation (PPN) guidance law. Consider Fig. 3,
in which $a_{ppn}$ denotes the acceleration applied in PPN guidance law and
rest all variables are same as in Fig. 2.

Figure 3: Diagram used for explaining the relation between the proposed
guidance law and PPN guidance law

The expression for PPN guidance law [19] is given by

$\displaystyle a_{ppn}=-NV_{a}\dot{\sigma}$ (12)

where $\dot{\sigma}$ is the line of sight rate and $N$ is the navigation
constant.

$\displaystyle\dot{\sigma}=-\frac{V_{a}\,sin{\rho}}{r_{a}}$ (13)

Using (12) to (13),

$\displaystyle a_{ppn}=NV_{a}(\frac{V_{a}\,sin\rho}{r_{a}})$ (14)

From (10) and (11),

$\displaystyle a_{c}=k_{1}r_{a}sin\rho$ (15)

Substituting (15) into (14) gives,

$\displaystyle a_{ppn}=(\frac{N}{k_{1}}){(\frac{V_{a}}{r_{a}})^{2}}a_{c}$ (16)

If we choose $k_{1}$ as given in (17), the proposed guidance law is same as
the PPN guidance law and is optimal for N=3 as shown in [20].

$\displaystyle k_{1}=\frac{N(V_{a})^{2}}{{r_{a}}^{2}}$ (17)

Thus the guidance law given in (11) takes the form

$\displaystyle a_{c}=(\frac{N({V_{a}})^{2}}{{r_{a}}^{2}})d$ (18)

and is same as the PPN guidance law. The guidance law given in (18) will
ensure that the MAV will reach the waypoint for $N>=2$, [19].

###  III-B Constraint on minimum turn radius

Micro air vehicles have a constraint on the minimum allowable turn radius. The
minimum value of the turn radius denoted by $R_{min}$ depends upon the
airspeed, stall angle of attack and the maximum thrust available. Since the
applied acceleration is perpendicular to the velocity vector, the MAV follows
a circular path with an instantaneous radius of curvature inversely
proportional to the magnitude of the applied acceleration. The constraint on
the applied acceleration $a_{c}$ can be given as [21]

$\displaystyle a_{c}<\frac{{V_{a}}^{2}}{R_{min}}$ (19)

Using (18) and (10), (19) can be written as

$\displaystyle(\frac{N(V_{a})^{2}}{{r_{a}}^{2}})(r_{a}sin\rho)<\frac{{V_{a}}^{2}}{R_{min}}$
(20)

After simplifying, (20) becomes

$\displaystyle\frac{Nsin\rho}{r_{a}}<\frac{1}{R_{min}}$ (21)

The minimum value for $N=2$. Hence (21) becomes,

$\displaystyle r_{a}>2R_{min}sin\rho$ (22)

At any instant, the relation given by (22) ensures that the minimum turn
radius constraint is not violated. A waypoint is considered to be not feasible
if the condition given in (22) is violated.

##  IV Integrated guidance and control for waypoint navigation in two
dimensional space

In separate guidance and control, the input to the actuator is computed in two
sequential stages from the navigation data. Whereas in integrated guidance and
control, the input to the actuator is computed in a single step from the
navigation data. The block diagram of integrated guidance and control (IGC) is
given in Fig. 4. The integrated guidance and control block generates the
actuator input ($U_{I}$) based on the output ($Y_{I}$) and the next waypoint
position ($x_{n},y_{n})$. $U_{P}$ is the actuator output that deflects the
elevator, rudder control surfaces and also changes the thrust input.

ActuatorMAV DynamicsIGC$(x_{f},y_{f})$$U_{P}$$Y_{I}$$U_{I}$ Figure 4:
Architecture of integrated guidance and control for MAV

The measured output $Y_{I}$ is given by

$\displaystyle Y_{I}=[q,\theta,h,p,r,\phi,d]^{T}$ (23)

The yaw angle $\psi$ , the current position of the MAV ($x,y$) and the next
waypoint position ($x_{f},y_{f}$) is used for the computation of the miss
distance $d$ as given in (7) to (10). Here the heading angle $\chi\approx\psi$
assuming the angle of side slip ($\beta$) to be small. Let $x_{1}=d$, and
$\dot{x_{1}}=\dot{d}=x_{2}$. The state variables used in the linear state
space model of integrated guidance and control logic are

$\displaystyle X_{I}=[X,X_{g},X_{a}]^{T}$ (24)

where $X$ is given in (2), $X_{g}=[\widetilde{x_{1}},\widetilde{x_{2}}]$ are
the states corresponding to the guidance logic and
$X_{a}=[{\widetilde{\delta}}_{ep},{\widetilde{\delta}}_{ep1},{\widetilde{\delta}}_{er},{\widetilde{\delta}}_{er1}]$
corresponds to second order elevator and rudder servo actuator transfer
function. Thus

$\displaystyle
X_{I}=[\widetilde{u},\widetilde{w},\widetilde{q},\widetilde{\theta},\widetilde{h},\widetilde{v},\widetilde{p},\widetilde{r},\widetilde{\phi},\widetilde{x_{1}},\widetilde{x_{2}},{\widetilde{\delta}}_{ep},{\widetilde{\delta}}_{ep1},{\widetilde{\delta}}_{er},{\widetilde{\delta}}_{er1}]^{T}$
(25)

In (25), $\widetilde{x_{1}},\widetilde{x_{2}}$ are the variables obtained
through the linearization of $x_{1}$ and $x_{2}$. The first step in the
linearization process is to find a relation between $\dot{x_{2}}$($\ddot{d}$)
and the acceleration command $a_{c}$. The second step is to find relation
between $a_{c}$ and the turn rate of the vehicle which is a function of
$\dot{\psi}$ and $\dot{\beta}$. These two steps will lead to the integration
of $X_{g}=[\widetilde{x_{1}},\widetilde{x_{2}}]$ to the state vector $X$ to
form $X_{I}$, since $\dot{\psi}$ and $\dot{\beta}$ is a function of $X$ . As a
first step, the commanded acceleration is simplified and is valid for small
angles of $\rho$ as given in the below theorem.
Theorem I: For small angle $\rho$, the commanded acceleration can be written
as

$a_{c}=\frac{1}{N-1}\ddot{d}$ (26)

Proof: The approximations $sin\rho=\rho$ and $cos\rho=1$ are done for small
angle $\rho$. From (10),

$\displaystyle d=r_{a}\rho$ (27)

Differentiating the above equation gives the following,

$\displaystyle\dot{d}=\dot{r_{a}}\rho+r_{a}\dot{\rho}$ (28)

From (9),

$\displaystyle\dot{r_{a}}=-V_{a}$ (29)

From (8),

$\displaystyle\dot{\rho}=-\dot{\chi}-\dot{\sigma}$ (30)

Using the relation $a_{c}=V_{a}\dot{\chi}$ and from (18)

$\displaystyle\dot{\chi}=\frac{NV_{a}d}{{r_{a}}^{2}}$ (31)

From (13) and (10) we obtain,

$\displaystyle\dot{\sigma}=-\frac{V_{a}d}{{r_{a}}^{2}}$ (32)

Substituting (31) and (32) into (30) gives,

$\displaystyle\dot{\rho}=-(N-1)\frac{V_{a}d}{{r_{a}}^{2}}$ (33)

From (29), (33) and (28) we obtain,

$\displaystyle\dot{d}=-\frac{NV_{a}d}{r_{a}}$ (34)

Differentiating the above equation and using (10) gives,

$\displaystyle\ddot{d}=-NV_{a}\frac{d}{dt}(\rho)$ (35)

Using (31) and (32), (35) can be rewritten as

$\displaystyle\ddot{d}=(N-1)\frac{{{NV_{a}}^{2}}d}{{r_{a}}^{2}}$ (36)

Noting that $a_{c}=\frac{N{V_{a}}^{2}d}{{r_{a}}^{2}}$, the desired result is
obtained.

The approximation of the non-linear guidance law that is valid for small
angles of $\rho$ is given in (26). Now as the second step to enable the
linearization process, a relationship between $a_{c}$ and $\dot{\psi}$ ,
$\dot{\beta}$ is derived. In Fig. 5, $X^{\prime}Y^{\prime}$ is the body fixed
lateral inertial plane of the MAV with X’ pointing towards north and Y’
pointing towards east, $X_{B}$ denotes the body X axis, $\beta$ is the
sideslip angle, $\psi$ is the yaw angle. The origin of the body fixed lateral
inertial plane is the current position of the MAV denoted by $(x,y)$. When an
acceleration $a_{c}$ is applied perpendicular to the velocity vector in the
yaw plane as shown in the Fig. 5, the MAV will move in a circular path in the
body fixed lateral inertial plane with an angular velocity
$\dot{\psi}+\dot{\beta}$ . Then the applied acceleration can be equated to
$V_{a}(\dot{\psi}+\dot{\beta})$ as given in (37).

Figure 5:  Body fixed lateral inertial plane of the MAV

The heading angle $\chi=\psi+\beta$. From the Fig. 5, we obtain

$\displaystyle a_{c}=V_{a}(\dot{\psi}+\dot{\beta})$ (37)

The expression for $\dot{\psi}$ and $\dot{\beta}$ is given in (38) and (39)
respectively.

$\displaystyle\dot{\psi}=qsin\phi sec\theta+rcos\phi sec\theta$ (38)
$\displaystyle\dot{\beta}=\frac{1}{V_{a}}(pw-ru+gcos\theta
sin\phi+\frac{Y_{a}}{m})$ (39)

where $Y_{a}$ is the side-force given by

$\displaystyle Y_{a}=\bar{Q}SC_{y}$ (40)

where $\bar{Q}=0.5\rho_{a}{V_{a}}^{2}$ is the dynamic pressure in $N/m^{2}$
and $C_{y}$ is the side-force coefficient given by

$\displaystyle
C_{y}=(C_{y\beta}\beta+C_{yp}(\frac{0.5b}{V_{a}})p+C_{yr}(\frac{0.5b}{V_{a}})r+C_{y\delta
r}\delta_{r})$ (41)

where $C_{y\beta}$, $C_{yp}$, $C_{yr}$ and $C_{y\delta_{r}}$ denotes the
change in side-force coefficient due to $\beta$, $p$, $r$ and $\delta_{r}$
respectively. From (37), (38) and (39) we can see that the applied
acceleration $a_{c}$ is a function of the state variables and the rudder
input. Thus the state equations for the guidance logic are

$\displaystyle\dot{x_{1}}=x_{2}$ (42) $\displaystyle\dot{x_{2}}=(N-1)a_{c}$
(43)

Equations (42) and (43) are linearized and augmented with the MAV dynamic
model to obtain the state space model for integrated guidance and control. Let
the linearized variables for $x_{1}$ and $x_{2}$ be $\widetilde{x_{1}}$ and
$\widetilde{x_{2}}$ respectively. Then the linear state equation and output
equation for the integrated guidance and control for two dimensional space is
given by

$\displaystyle\dot{X_{I}}=AX_{I}+BU_{I}$ (44) $\displaystyle Y_{I}=CX_{I}$
(45)

where $X_{I}$ is given in (25) and

$U_{I}=[\widetilde{\delta_{e}},\widetilde{\delta_{r}},\widetilde{\delta_{t}}]^{T}$
(46)
$Y_{I}=[\widetilde{q},\widetilde{\theta},\widetilde{h},\widetilde{p},\widetilde{r},\widetilde{\phi},\widetilde{x_{1}}]^{T}$
(47)

###  IV-A Controller design for IGC framework

For the trim conditions at $V_{a}=8m/s$, $\dot{\psi}=0.8rad/s$ and climb rate
$\dot{h}=1m/s$, the poles of linear state space model corresponding to IGC
framework is given in Table II. The spiral mode is unstable for the open loop
system.

TABLE II: Poles of linear state space model - IGC framework Short period mode ($\omega_{sp}$, $\zeta_{sp}$) | 35.6 $rad/s$, 0.285
---|---
Phugoid mode ($\omega_{ph}$, $\zeta_{ph}$) | 1.15 $rad/s$, 0.219
Dutch roll mode, ($\omega_{dr}$, $\zeta_{dr}$) | 43.8 $rad/s$, 0.337
Roll subsidence mode | -2.44
Spiral mode | 1.51 (unstable)
Elevator servo ($\omega_{e}$, $\zeta_{e}$) | 48.7 $rad/s$, 0.742
Rudder servo ($\omega_{r}$, $\zeta_{r}$) | 48.7 $rad/s$, 0.742
Poles for $\widetilde{h}$, $\widetilde{x_{1}},\widetilde{x_{2}}$ | 0, 0, 0

The controller is synthesized in $H_{\infty}$ robust control framework as it
can handle uncertainties associated with the MAV dynamics. The architecture of
$H_{\infty}$ control for IGC framework is shown in Fig. 6. Here the reference
signal $Y_{r}=[\widetilde{h}_{r},\,\widetilde{x_{1}}_{r}]^{T}$, and the
measured outputs $Y_{1}=[\widetilde{h},\,\widetilde{x_{1}}]^{T}$,
$Y_{2}=[\widetilde{q},\widetilde{\theta},\widetilde{p},\widetilde{r},\widetilde{\phi}]^{T}$.
The weighting transfer functions $W_{1}$, $W_{2}$ and $W_{3}$ are low pass
filter, high pass filter and constant respectively [22]. $Z_{1}$, $Z_{2}$,
$Z_{3}$ are the performance outputs to be minimized.

$K_{c}$MAV$U_{I}$$W_{1}$$W_{3}$$W_{2}$$Z_{2}$$Z_{1}$$Z_{3}$$Y_{r}$$E$$Y_{I}$$-$$Y_{1}$$Y_{2}$
Figure 6: $H_{\infty}$ control architecture for IGC framework

The controller synthesis is performed in the discrete time domain using the
algorithm described in [5]. The controller sampling time is 0.02 seconds (the
standard for servo control). The algorithm is summarized as given below. The
discrete time equivalent linear state space model of the generalized plant
shown in Fig. 6 is given in (48) to (52).

$\displaystyle{X_{dg}}(k+1)=A_{dg}X_{dg}(k)+B_{du}{U_{dc}}(k)+B_{dw}{W_{dd}}(k)$
(48)
$\displaystyle{Z_{d1}}(k)=C_{d1}X_{dg}(k)+D_{d11}{U_{dc}}(k)+D_{d12}{W_{dd}}(k)$
(49)
$\displaystyle{Z_{d2}}(k)=C_{d2}X_{dg}(k)+D_{d21}{U_{dc}}(k)+D_{d22}{W_{dd}}(k)$
(50)
$\displaystyle{Z_{d3}}(k)=C_{d3}X_{dg}(k)+D_{d31}{U_{dc}}(k)+D_{d32}{W_{dd}}(k)$
(51) $\displaystyle Y_{d}(k)=C_{d}X_{dg}(k)$ (52)

The controller design objective is to synthesis a SOF control law,

${U_{dc}}(k)=F_{d}Y_{d}(k)$ (53)

that minimizes the weighted sensitivity, weighted complimentary sensitivity
and weighted control input sensitivity function. The gain matrix $F_{d}$ is
obtained by solving the LMI (linear matrix inequality) given in (56). The
matrix $N_{d}$ is selected using the genetic algorithm (GA) for minimizing a
performance index considering the closed loop damping ratio and stability of
spiral mode and the resultant LMI is solved for $P_{d}$ and $F_{d}$. The
matrix $P_{d}>0$ for a feasible solution $F_{d}$ that stabilizes the closed
loop plant.

$\displaystyle\left(\begin{array}[]{cc}-P_{d}&-(A_{dg}+B_{du}F_{d}C_{d})^{T}N_{d}^{T}\\\
-N_{d}(A_{dg}+B_{du}F_{d}C_{d})&P_{d}-N_{d}-N_{d}^{T}\end{array}\right)<0$
(56)

The controller gain matrix $F_{d}$ is given in (57).

$F_{d}=[F_{d1},\,F_{d2}]$ (57)
$F_{d1}=\left[\begin{array}[]{ccc}-0.0060&-0.3378&-0.8032\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr 0&0.0987&0\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr-1.152&36.54&0.0042\end{array}\right]$ (58)
$F_{d2}=\left[\begin{array}[]{cccc}0.0002&-0.0007&0&-0.1008\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr-0.0528&0.0008&-0.0937&-0.0002\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr-5.325&-19.1&0.0004&0.0091\end{array}\right]$ (59)

For the trim conditions at $V_{a}=8m/s$, $\dot{\psi}=0.8rad/s$ and climb rate
$\dot{h}=1m/s$, the dynamic modes of IGC closed loop linear state space model
in continuous time domain is given in Table III. From the Table III, we can
see that the spiral mode is stabilized and the damping ratio of other modes
have improved when compared to open loop plant. The three poles at the origin
corresponding to $\widetilde{h}$, $\widetilde{x_{1}}$ and $\widetilde{x_{2}}$
are placed in the left half of the s-plane.

TABLE III: Poles of closed loop IGC linear state space model Short period mode ($\omega_{sp}$, $\zeta_{sp}$) | 36.8 $rad/s$, 0.404
---|---
Phugoid mode ($\omega_{ph}$, $\zeta_{ph}$) | 2.25 $rad/s$, 0.323
Dutch roll mode, ($\omega_{dr}$, $\zeta_{dr}$) | 30.0 $rad/s$, 0.429
Roll subsidence mode | -15.5
Spiral mode | -3.24
Pole for $\widetilde{h}$ | -1.50
Poles for $\widetilde{x_{1}}$, $\widetilde{x_{2}}$ ($\omega$, $\zeta$) | 0.1 $rad/s$, 0.343
Elevator servo ($\omega_{e}$, $\zeta_{e}$) | 44.4 $rad/s$, 0.671
Rudder servo ($\omega_{r}$, $\zeta_{r}$) | 49.5 $rad/s$, 0.603

###  IV-B Waypoint navigation algorithm using IGC framework

Inputs to the algorithm are the current $(x,y)$ position of the MAV obtained
from GPS, the desired waypoint position $(x_{f},y_{f})$, altitude $h$ obtained
from altimeter, rate gyro outputs p, q, r and Euler angles $\phi$, $\theta$,
$\psi$. The steps are summarized in Algorithm I. For a given waypoint, the
angle $\rho$ is computed from (8) and a constant roll angle command is issued
if $|\rho|>20^{0}$. Here the small angle approximation $sin\rho=\rho$ and
$cos\rho=1$ is used for $|\rho|\leq 20^{0}$ and the expression given in (26)
is valid for the same. The roll angle control for 150 mm MAV given in [3] is
used here. The altitude is maintained constant while executing the roll angle
command. The MAV follows a constant roll angle command
$\phi_{ref}=sgn(\rho)\phi_{r}$ for $|\rho|>20^{0}$. Here $\phi_{r}$ is a
constant value less than the maximum allowed roll angle $\phi_{max}$ of the
MAV. When $|\rho|\leq 20^{0}$, the feasibility of reaching the waypoint is
verified using (19) and subsequently the control law given in (53) is employed
as the IGC framework is valid. When the range ($r_{a}$) is less than a pre-
defined threshold $r_{t}$, the waypoint is considered to be reached and the
entire process is repeated for a new waypoint.

Algorithm I:
STEP I: Compute $\rho$ from (8).
STEP II: else-if $|\rho|>20^{0}$ Issue a roll angle hold command
$\phi_{ref}=sgn(\rho)\phi_{r}$, where $0<\phi_{r}<\phi_{max}$ and goto STEP I.
else-if $r_{a}>2R_{min}sin\rho$ Use the controller structure given in (53) and
goto STEP III.
Current waypoint is not feasible, goto STEP IV.
STEP III: else-if $r_{a}<r_{t}$ Waypoint reached, goto STEP IV else goto STEP
I.
STEP IV: Switch to a new waypoint and goto STEP I.

##  V Numerical simulation results for waypoint navigation

The waypoint navigation algorithm is simulated for the non-linear six degrees
of freedom (6DOF) model of the 150 mm MAV given in [1]. Flight test results
are presented in [1], validating the non-linear model used for simulation. Two
simulations are performed in which the first one is to demonstrate straight
line path following and for the second one, the objective is to follow a
rectangular path. Constraints are put on the actuators while performing the
simulations. The range of elevator deflection is $[-35^{o},15^{o}]$, rudder
deflection is $[-25^{o},25^{o}]$ and motor thrust is 0.0-0.45 N. The second
order transfer function model is used for elevator and rudder servos. The
control input update is done at every 0.02 second. The $XY$ position update is
done at every 1 second. All other outputs are sampled at a rate of 0.02
second. The waypoint navigation logic is started when the MAV reach at an
altitude of 20 m from the ground level. Altitude feedback is used to maintain
the desired altitude of 20 m. For both the simulations, the initial conditions
are shown in Table IV.

TABLE IV: Initial conditions used for the non-linear simulation Variable |  Value
---|---
$V_{a}$, $\alpha$, $\beta$ |  8 m/s, 0.2286 rad, -0.0532 rad
$\delta_{e}$, $\delta_{r}$, $\delta_{t}$ |  -0.2592 rad, 0.1255 rad, 0.231 N
$\phi$, $\theta$, $\psi$ |  -0.0361 rad, 0.2304 rad, 0 rad
p, q, r |  0 rad/s, 0 rad/s, 0 rad/s
u, v, w |  7.78 m/s, -0.43 m/s, 1.81 m/s
$x$, $y$, $h$ |  0 m, 0 m, 20 m

The objective of the first simulation is to follow a straight line. The XY
(North-East) path of the MAV is shown in Fig. 7. In Fig. 7, red circles
represent the waypoints commanded and the solid blue line represents the path
followed by the MAV. The plot of miss distance and the angle $\rho$ is given
in Fig. 8. During the initial stage, the miss distance is high as the first
waypoint is not along the line-of-sight (LOS) of the MAV. After 25 seconds, we
can see that the value of $\rho$ and miss distance is closer to zero as the
MAV got aligned to the straight line. The small angle approximation used in
the derivation of the IGC framework is used for $\rho<20^{o}$. Initially when
$\rho>20^{o}$ till 1.4 seconds, a constant roll angle command of
$\phi_{r}=10^{o}$ is employed. A plot of Euler angles and angular rates are
given in Fig. 9. The roll angle $\phi$ goes to $10^{o}$ as shown in Fig. 9
accordingly. The plot of control inputs is given in Fig. 10. The peak
amplitude of elevator deflection is less than $-16^{o}$ and that of the rudder
is less than $8^{o}$. The peak value of elevator deflection is less than
$50\%$ of the saturation value of $-35^{o}$. The peak value of rudder
deflection is less than $33\%$ of the saturation value of $25^{o}$. Thus the
control inputs are well within the saturation limits. The variation airspeed
is between 7.7 m/s and 8.2 m/s as seen from Fig. 10 and is closer to the trim
velocity of 8.0 m/s.

Figure 7: XY plot for MAV straight line following Figure 8: Miss distance and
$\rho$ angle plots for MAV straight line following Figure 9: Euler angles and
body axis angular rates for MAV straight line following Figure 10: Control
inputs and airspeed for straight line following

The objective of the second simulation is to follow a rectangular path.
Following a rectangular path enables the MAV to perform surveillance over a
given area. The rectangular path considered is of size 100m x 200m. The path
followed by the MAV is shown in Fig. 11. In Fig. 11, small red circles denotes
the waypoints along the rectangular path and blue line represents the path
followed by the MAV. The starting point of the flight path is at (0,0). The
MAV takes a turn after reaching the first waypoint at (80,0). The plot of $d$
and $\rho$ given in Fig. 12 show oscillations whenever MAV switches to a new
waypoint. The miss distance goes to zero as the MAV approaches a waypoint. The
plot of control inputs and airspeed for rectangular path following is shown in
Fig. 13. The airspeed is between 7-8.5 $m/s$ for most of the time duration of
the MAV flight.

Figure 11: XY plot for rectangular path following Figure 12: Miss distance and
$\rho$ angle plots for rectangular path following Figure 13: Control inputs
and airspeed for rectangular path following

##  VI Conclusion

The paper presents a novel waypoint navigation algorithm in IGC framework. The
IGC methodology combines unstable MAV coupled dynamics with the pure
proportional navigation guidance law. The resulting linear closed loop system
is found to be stable with a static output feedback control law. The proposed
waypoint navigation algorithm handles the minimum turn radius constraint of
the MAV and determines the feasibility of reaching a given waypoint. The
waypoint navigation algorithm does not require any iterative methods to
compute the control inputs. The high fidelity nonlinear numerical simulations
for straight line and rectangular path following validates the utility of the
algorithm. The airspeed of the MAV is maintained well above the stall value
while performing the path following. The extension of this algorithm for
waypoint navigation in three-dimensional space would be interesting future
work.

## Acknowledgment

The authors would like to thank National Programme for Micro Air Vehicle
(NPMICAV) and Aeronautical Development and Research Board (ARDB), government
of India for funding the project.

## References

  * [1] K. Harikumar, J. V. Pushpangathan, T. Bera, S. Dhall, and M. S. Bhat. ”Modeling and closed loop flight testing of a fixed wing micro air vehicle,” Micromachines, Vol. 9, No. 3, p.111, 2018.
  * [2] D.R. Nelson, D.B. Barber, T.W. McLain and R.W. Beard, ”Vector field path following for miniature air vehicles,” IEEE Transactions on Robotics, Vol. 23, No. 3, pp. 519-529, 2007.
  * [3] K. Harikumar, Sidhant Dhall and M. Seetharama Bhat, ”Nonlinear modeling and control of coupled dynamics of a fixed wing micro air vehicle,” Indian Control Conference, Hyderabad, India, January 2016.
  * [4] Sikha Hota and Debasish Ghose, ”Time-Optimal Convergence to a Rectilinear Path in the Presence of Wind,” Journal of Intelligent Robotic Systems Vol. 74, pp. 791-815, 2014.
  * [5] K. Harikumar, Sidhant Dhall and M. Seetharama Bhat, ”Design and experimental validation of a Robust output feedback control for the coupled dynamics of a micro air vehicle,” Springer International Journal of Control, Automation and Systems, vol. 17, no. 1, pp. 155-167, 2019.
  * [6] H. Yan and H. Ji, ”Integrated guidance and control for dual-control missiles based on small-gain theorem,” Automatica Vol. 48, pp. 2686-2692, 2012.
  * [7] M. Xin, S. N. Balakrishnan and E. J. Ohlmeyer, ”Integrated Guidance and Control of Missiles With $\theta$\- D Method,” IEEE Transactions on Control Systems Technology, Vol. 14, No. 6, pp. 981-992, 2006.
  * [8] Y. Ke, T. Qingke, L. Siyuan, L. Qingdong and R. Zhang, ”Integrated guidance and control design based on motion tracking,” Proceedings of the 35th Chinese Control Conference, Chengdu, China, pp. 5721-5725, 2016.
  * [9] P.K. Menon and E.J. Ohlmeyer, ”Integrated design of agile missile guidance and autopilot systems,” Control Engineering Practice, Vol. 9, No. 10, pp. 1095-1106, 2001.
  * [10] C. Zhang and W. Yun-jie, ”Non-singular terminal dynamic surface control based integrated guidance and control design and simulation,” ISA transactions, Vol. 63, pp. 112-120, 2016.
  * [11] G. Jianguo, Y. Xiong, and J. Zhou, ”A new sliding mode control design for integrated missile guidance and control system,” Aerospace Science and Technology, Vol.78, pp. 54-61, 2018.
  * [12] I. Kaminer, A. Pascoal, E. Hallberg and C. Silvestre, ”Trajectory tracking for autonomous vehicles: an integrated approach to guidance and control,” Journal of Guidance, Control and Dynamics, Vol. 21, No. 1, pp. 29-38, 1998.
  * [13] T. Yamasaki, S.N. Balakrishnan and H. Takano, ”Separate channel integrated guidance and control for automatic path following,” Journal of Guidance, Control and Dynamics, Vol. 36, No. 1, pp. 25-34, 2013.
  * [14] R. Padhi, P. R. Rakesh and R. Venkataraman, ”Formation Flying with Nonlinear Partial Integrated Guidance and Control,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 50, No. 4, pp. 2847-2859, 2014.
  * [15] X. Liu, W. Huang, L. Du , P. Lan and Y. Sun, ”Three-Dimensional Integrated Guidance and Control for BTT Aircraft Constrained by Terminal Flight Angles, ” 27th Chinese Control and Decision Conference (CCDC), pp.107-112, 2015.
  * [16] M. Z. Babar, R. Samar, A. I. Bhatti and M. Baglietto, ”Robust Integrated Lateral Guidance and Control of UAVs,” 20th IEEE International Conference on Methods and Models in Automation and Robotics (MMAR), pp. 523-528, 2015.
  * [17] J. V. Pushpangathan, M. S. Bhat, and K. Harikumar, ”Effects of Gyroscopic Coupling and Countertorque in a Fixed-Wing Nano Air Vehicle,” Journal of Aircraft, vol. 55, no.1, pp. 239-250, 2018.
  * [18] J. Roskam, ”Airplane design,” Part 6, Roskam Aviation and Engineering Corporation, pp. 371-461, 1990.
  * [19] N.A. Shneydor, ”Missile Guidance and Pursuit: Kinematics, Dynamics and control,” Horwood publishing, England, pp. 109-117, 1998.
  * [20] R. Yanushevsky, ”Guidance of Unmanned Aerial Vehicles,” CRC Press, pp. 9-18, 2011\.
  * [21] J.L. Meriam, L.G. Kraige, ”Engineering mechanics - Dynamics,” Sixth edition, Wi- ley India, pp. 55-57, 2008.
  * [22] H. Jiankun, C. Bohn and H.R. Lu, ”Systematic $H_{\infty}$ weighting function selection and its application to the real-time control of a vertical take-off aircraft,” Control Engineering Practice, Vol. 8, No. 3, pp. 241-252, 2000.