# Identifying and Explaining the Resilience of Ecological Networks

Cailan Jeynes-Smith  School of Mathematical Sciences, Queensland University of
Technology, Brisbane, 4000, Australia<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Michael Bode  School of Mathematical Sciences, Queensland
University of Technology, Brisbane, 4000, Australia.  Safeguarding Antartica’s
Environmental Future, Queensland University of Technology, Brisbane, 4000,
Australia.  Robyn P. Araujo  School of Mathematical Sciences, Queensland
University of Technology, Brisbane, 4000, Australia.
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

Resilient ecological systems will be better able to maintain their structure
and function in the emerging Anthropocene. Estimating the resilience of
different systems will therefore provide valuable insight for conservation
decision-makers, and is a priority goal of resilience theory. Current
estimation methods rely on the accurate parameterisation of ecosystem models,
or the identification of important motifs in the structure of the ecological
system network. However, both of these methods face significant empirical and
theoretical challenges. In this paper, we adapt tools developed for the
analysis of biochemical regulatory networks to prove that a form of resilience
- robust perfect adaptation - is a property of particular ecological networks,
and to explain the specific process by which the ecosystem maintains its
resilience. We undertake an exhaustive search for robust perfect adaptation
across all possible three-species ecological networks, under a generalised
Lotka-Volterra framework. From over 20,000 possible network structures, we
identify 23 network structures that are capable of robust perfect adaptation.
The resilient properties of these networks provide important insights into the
potential mechanisms that could promote resilience in ecosystems, and suggest
new avenues for measuring and understanding the property of ecological
resilience in larger, more realistic socioecological networks.

###### keywords:

resilience, adaptation, socioecological, ecosystem, management

Keywords: resilience, adaptation, socioecological, ecosystem, management

Author contributions statement

CJS conceptualised the idea, RPA contributed to the modelling framework, CJS
implemented methods, analysed results, created visualisations and produced the
first draft of the manuscript, all authors contributed to the review and
editing of the manuscript, and MB, and RPA supervised the project.

Acknowledgements:

Computational resources and services used in this work were provided by the
eResearch Office, Queensland University of Technology, Brisbane, Australia.

Funding: Robyn P. Araujo is the recipient of an Australian Research Council
(ARC) Future Fellowship (project number FT190100645) funded by the Australian
Government, Cailan Jeynes-Smith is supported by an Australian Government
Research Training Program Scholarship.

Conflict of Interest Statement:

We declare we have no competing interests.

Data Statement

Data sharing is not applicable to this article as no new data were created or
analyzed in this study. All code used to generate the results and figures in
this article can be found on Github:
https://github.com/JeynesSmith/PerfectResilience.git

## Introduction

As the Anthropocene drives accelerating global change, resilience is an
important and desirable characteristic of ecological and socioecological
systems [1, 2, 3]. Resilience has a range of definitions, but it primarily
refers to the ability of an ecosystem to return to a particular set of
ecological states following changes in environmental or social (policy)
pressures (Figure 1). Resilience is generally used to describe a property
exhibited by a particular ecological or socioecological system [4, 5, 6, 7, 8,
9, 10, 11], or as a measure of how well the population will recover from
perturbation [12, 13, 14, 15, 16]. Understanding whether, and to what extent,
an ecological or socioecological system is resilient helps us understand
whether an ecosystem is at risk of collapse, or how far it can be pushed by
environmental or social changes [1, 2, 3].

Figure 1: An example of the perfect resilience behaviour, where a species’
abundance is able to consistently return exactly to its original abundance
following a disturbance. This behaviour is equivalent to robust perfect
adaptation in biochemical reaction networks. Resilience, or imperfect
adaptation are similar behaviours in which the abundance does not strictly
return to the exact original abundance, but instead returns to within some
surrounding region. This figure is generated by simulating a version of the
network in Figure 3(i). We simulate the system for one hundred time steps
before doubling a stimulus ($S$) which disturbs the networks, specifically
$S=\\{1,2,4\\}$. Parameters: $r_{I}=0.61$, $d_{O}=-0.07$, $a_{13}=-0.37$,
$a_{21}=0.27$, $a_{22}=-0.81$, $a_{31}=0.32$, $a_{32}=-0.69$.

A major hurdle to the use of resilience theory in conservation decision-making
is its measurement. A primary approach is to create and parameterise a
dynamical systems model of the ecological system, and then to simulate its
response to perturbations [4, 5]. However, the complexity of ecological
systems, paired with sparse and noisy data sets, make the process of parameter
estimation incredibly difficult [17, 18]. If resilience could be defined as a
property of the structure of the ecological system, then parameter
identifiability problems could be overcome.

Previous studies into resilience of socioecological systems have proposed
theoretical structures - network motifs - that underpin resilience [19, 20].
However, this research is relatively atheoretical, with resilient network
motifs being proposed on the basis of intuition, and justified by statistical
association with observations of resilient dynamics. This limits the insights
that these methods can offer.

By contrast, resilience has been studied extensively in biochemical reaction
networks, where it is called adaptation. Resilience/adaptation is a relatively
common dynamical property in cellular systems, where it has evolved to
maintain the processes of life in a stochastic environment where external
stimuli are constantly perturbing the abundances of molecules in the system.
It has been observed in applications scaling from chemotaxis in single-celled
organisms [21, 22, 23, 24, 25, 26, 27, 28, 29] to complex sensory systems [30,
31, 32, 33], while the loss of adaptation has been linked to cancer
progression, and substance abuse [34, 35, 36]. Resilience/adaptation is
further categorised as either imperfect adaptation, where a particular element
of the system returns within some accepted tolerance of its original (pre-
stimulus) abundance [37, 38, 39, 40, 41], or as robust perfect adaptation,
where that target element returns precisely to the exact same (pre-stimulus)
abundance [42, 43, 44, 45, 46, 47, 48, 49, 50, 51] (see Figure 1).
Importantly, robust perfect adaptation is a property of the network
interaction structure, as opposed to a property of a specific parameterisation
of the network [42].

In this study we apply analytical tools from biochemical reaction network
theory to understand the process of resilience in ecological networks, and to
assess whether ecological systems can exhibit robust perfect adaptation, and
under what circumstances. To more closely match similar concepts in the
existing ecological literature, we herein refer to ‘perfect resilience’ as a
mathematically equivalent property to robust perfect adaptation. We study
ecological systems at the ‘operational layer’ [19], where the effect of
policies directly impact populations in an ecosystem. We perform an extensive
search of all three-species ecosystems using a novel application of the
generalised Lotka-Volterra equations - a commonly used mechanistic framework
in ecological modelling. We ask, are there network motifs that promote perfect
resilience under an ecological framework, and will the corresponding networks
still be representative of in situ ecological systems? If possible, these
specific motifs could be identified in ecosystems instead of current, more
generalised motifs [19].

## Methods

The structures capable of robust perfect adaptation (RPA) in biochemical
reaction networks have recently been identified in full generality [42, 44].
Ma et al. [37] performed extensive numerical simulations of three-node
chemical reaction (signalling) networks with Michaelis-Menten kinetics, and
identified that there were two motifs structures that support RPA (at least
approximately, as determined by tight numerical thresholds): a negative
feedback loop with buffer node (NFBLB), and an incoherent feedforward loop
with proportioner node (IFFLP). Araujo and Liotta [42] later determined that
for networks of arbitrary size and complexity, and for arbitrary interaction
kinetics, two-well defined subnetwork structures, or ‘modules’, constitute a
topological basis for RPA in any network: Opposer modules, which are
generalisations of three-node NFBLBs, and are feedback structured subnetworks;
and Balancer modules, which are generalisations of three-node IFFLPs, and are
feedforward-structured subnetworks. All RPA-capable networks, no matter how
large or complex, and no matter the ‘kinetics’ of the interacting elements,
are necessarily decomposable into these two special modular subnetwork
structures. See [42] for a detailed description of these mechanisms. More
recently, the intricate biochemical reaction structures that are compatible
with these overarching RPA-conferring mechanisms have also been determined
[44].

The biochemical networks that are capable of adaptation are often built from
enzyme-mediated reactions [37, 42, 40, 39, 45], where an enzyme combines with
a substrate molecule to form an intermediate complex. This complex can either
dissociate into the original molecules, or convert the substrate into a
modified (product) form, releasing the enzyme unmodified. Importantly, there
is no equivalent to these reactions in ecological systems, so it is unclear
how, or if, ecological networks could generate adaptive behaviours without
these fundamental reactions.

We provide an extensive study of all three-species interaction networks under
a generalised Lotka-Volterra framework containing the three arbitrary species
(or functional groups [52]), $I$, $M$, and $O$, and a stimulus (cause of the
change), $S$. We check for network structures in which the species, $O$, has
the perfect resilience property. The stimulus represents an external influence
on the system, which could be an intervention such as a harvest, or an
environmental factor such as a heatwave, and can affect any combination of
species in the system.

The full set of interactions in a three-species network modelled with the
generalised Lotka-Volterra equations are given by,

$\displaystyle\frac{\text{d}I}{\text{dt}}$
$\displaystyle=\overbrace{r_{I}I}^{\text{Intrinsic
Growth}}-\overbrace{a_{II}I^{2}}^{\text{Intraspecific
Competition}}+\overbrace{a_{IM}IM+a_{IO}IO}^{\text{Interspecific
Interactions}}+\overbrace{d_{I}SI}^{\text{Stimulus}},$ (1)
$\displaystyle\frac{\text{d}M}{\text{dt}}$
$\displaystyle=r_{M}M-a_{MM}M^{2}+a_{MI}IM+a_{MO}MO+d_{M}SM,$ (2)
$\displaystyle\frac{\text{d}O}{\text{dt}}$
$\displaystyle=r_{O}O-a_{OO}O^{2}+a_{OI}IO+a_{OM}MO+d_{O}SO,$ (3)

where: $I$, $M$, and $O$ are the abundances of the interacting species (see
Figure 2(a)); $r_{i}$ is the intrinsic growth rate of species $i$; $a_{ij}$ is
the (per-capita) interaction constant for how species $i$ is affected by
species $j$; and $d_{i}$ is the interaction constant for how the stimulus, S,
affects species $i$ ($d_{i}=0$ when species $i$ is not affected by the
stimulus).

The generalised Lotka-Volterra equations consist of three key terms: intrinsic
growth, intra-specific competition (self-regulation), and inter-specific
interactions. The intrinsic growth term models positive influences on a
populations growth. This is often included for lower trophic-level species
such as vegetation, where influences like rainfall and soil nutrients are
implicitly modelled. Intra-specific competition represents a limiting or
dampening term, whereby competition for the same resources within a population
will ultimately limit its growth. Lastly, inter-specific interactions are all
of the interactions between species $i$ and $j$. These can be positive or
negative, where the signs of a pair of terms, $(a_{ij},a_{ji})$, are
indicative of the type of interaction between species. Some common examples
include: competition, $a_{ij},a_{ji}<0$; mutualism, $a_{ij},a_{ji}>0$; and
predator-prey, $a_{ij}>0$, $a_{ji}<0$ where population $i$ is the predator. We
illustrate an application of the generalised Lotka-Volterra equations and the
graphical representation of a specific network in Figure 2.

Figure 2: (a) The graphical representation of a three-species network with the
species, $I$, $M$, and $O$, and a stimulus $S$, and (b) the associated
generalised Lotka-Volterra equations. This network depicts a system in which
the species $O$ predates on both $I$ and $M$. Species $I$ is affected by the
stimulus, and $O$ is the output of the network. In the network diagram (a) we
represent a positive interaction by a pointed arrow and a negative interaction
by a flat-ended arrow. The interactions and constants have been coloured based
on the type of interaction: intraspecific competition, intrinsic growth rates,
interspecific interactions (coloured light to dark blue respectively), and
stimulus interactions (coloured red). The green arrow indicates that species
$O$ is the output of the network, however this notation is dropped in later
networks to simplify diagrams.

We systematically examine the capacity for perfect resilience in networks
constructed from every possible combination of terms in Equations (1)-(3)
based on a full factorial design. In the case where interactions are present
or absent, i.e. not considering the sign of interactions, there are a total of
$5\times 2^{12}=20,480$ possible network structures, since there are twelve
parameters and five combinations in which the stimulus affects the network.
When accounting for the sign of interactions, i.e. considering interactions as
either absent, positive, or negative, there is a total of $20\times
2^{6}\times 3^{6}=933120$ possible network structures. We therefore require
efficient methods to test for the capacity of these networks to exhibit
perfect resilience.

Ma et al. [37] determined the capacity for imperfect adaptation in over 16,000
networks using metrics based on the simulated behaviour of 20,000 parameter
sets for each network. This method is computationally expensive and depends
heavily on an extensive sampling of parameter space for every individual
network motif. In the current study, we automate algebraic methods recently
developed by Araujo and Liotta [44] which can definitively and accurately
determine whether a given network motif has the capacity for RPA while
handling all parameters symbolically.

The cornerstone of this method is the recognition that any RPA-capable system,
with dynamical rate equations $f_{1},\ldots,f_{n}$ (assumed to be polynomial
functions of the interacting elements, $x_{1},\ldots,x_{n}$), is characterised
by an ideal whose two-variable geometric projection assumes the special form
of an RPA polynomial (see [44] for full details). In this case, RPA, and thus
perfect resilience, requires the existence of three polynomials,
$(p_{1},p_{2},p_{3})\subset\mathbb{R}[I,M,O]$, such that

$\displaystyle
p_{1}\frac{dI}{dt}+p_{2}\frac{dM}{dt}+p_{3}\frac{dO}{dt}=f(S,O)(O-k),$ (4)

where $k$ is a rational function of parameters, and $f(S,O)$ is a polynomial
(known as the ‘pairing function’ [44]) that is non-vanishing on a suitably
extensive region of the positive orthant. The right-hand side of Equation (4)
has the form of an RPA polynomial. If $(p_{1},\ p_{2},\ p_{3})$ can be found
that satisfy Equation (4) for any given network, then the network in question
has the capacity for perfect resilience, with the output steady-state value of
$k$ for all disturbances and all parameter choices. For the small network
structures considered here, the existence of suitable $p_{1},\ p_{2},\ p_{3}$
can be determined algorithmically by computation of a Gröbner basis [44] with
an elimination monomial ordering, and with variables $S$ and $O$ ordered last.
We automate this process in Matlab using the lexicographic monomial ordering,
and then use Matlab’s symbolic toolkit to automatically determine if a non-
zero projection onto $S$ and $O$ exists, and if the projection can be
factorised into an RPA polynomial (Equation (4)) based on the monomials. We
reject any networks for which there is no non-trivial two-variable projection,
or for which its two-variable projection is not an RPA polynomial. All code
developed for this study is provided at
https://github.com/JeynesSmith/PerfectResilience.git.

We provide two examples in which a network does, and does not, have the
capacity for perfect resilience based on the above projection test. The
network in Figure 3(i) (without $r_{O}$) is capable of achieving perfect
resilience. By calculating the Gröbner basis for this network, we obtain the
following projection,

$p_{1}\frac{dI}{dt}+p_{2}\frac{dM}{dt}+p_{3}\frac{dO}{dt}=c_{1}\mathbf{O^{2}S^{2}}+c_{2}\mathbf{O^{2}S}+c_{3}\mathbf{O^{2}}+c_{4}\mathbf{OS^{2}}+c_{5}\mathbf{OS}+c_{6}\mathbf{O},$

where $c_{i}$ is some function of the network parameters. The projection can
be factorised into a form which matches the right-hand-side of Equation (4),

$\displaystyle
c_{1}\mathbf{O^{2}S^{2}}+c_{2}\mathbf{O^{2}S}+c_{3}\mathbf{O^{2}}+c_{4}\mathbf{OS^{2}}+c_{5}\mathbf{OS}+c_{6}\mathbf{O}=O(S^{2}+c_{7}S+c_{8})(O-k),$
(5)

where $f(O,S)=O(S^{2}+c_{7}S+c_{8})$ to match Equation (4). Since the
factorisation exists, this network has the capacity for perfect resilience.

As a second example, the network from Figure 2 has the following Gröbner basis
projection,

$\displaystyle p_{1}\frac{dI}{dt}$
$\displaystyle+p_{2}\frac{dM}{dt}+p_{3}\frac{dO}{dt}=$ $\displaystyle
c_{1}\mathbf{O^{5}}+c_{2}\mathbf{O^{4}S}+c_{3}\mathbf{O^{4}}+c_{4}\mathbf{O^{3}S^{2}}+c_{5}\mathbf{O^{3}S}+c_{6}\mathbf{O^{3}}+c_{7}\mathbf{O^{2}S^{2}}+c_{8}\mathbf{O^{2}S}+c_{9}\mathbf{O^{2}}+c_{10}\mathbf{OS^{2}}+c_{11}\mathbf{OS}+c_{12}\mathbf{O}.$

There is no factorisation in which the right-hand-side of this equation will
match the form of Equation (4), even when based on the monomials ($S$ and $O$)
alone, and therefore this network does not have the capacity for perfect
resilience.

While the above condition is necessary for perfect resilience, we must still
determine whether the setpoint, $O=k$, is a feasible and stable steady state.
Feasibility ensures that all species have a positive abundance at steady
state, while stability ensures that our ecosystem can move towards that steady
state. Since the projection test has significantly reduced the number of
networks (20480 to 1072 networks, see Supplementary Figure 3), we test for
stability and feasibility using a sampling approach similar to Ma et al. [37]
in which we randomly select $4\times 10^{3}$ parameter sets. Parameters are
selected from uniform distributions, $x\in(-1,1)$ for interspecific or
stimulus interaction constants, $x\in(0,1)$ for intrinsic growth rates, and
$x\in(-1,0)$ for intraspecific competition constants. We calculate the steady
states of the network then substitute in random parameter values and check for
feasibility i.e. there exists a steady state in which every species has a
positive abundance. If the steady state is feasible, then we check stability
using Lyaponuv stability i.e. negative real components for all eigenvalues of
the Jacobian matrix [53]. If both stability and feasibility conditions are
met, then the parameter set is saved as successful and we continue for the
remaining $4\times 10^{3}$ parameter sets or until we obtain ten successful
parameter sets. We define that a network is capable of perfect resilience if
it has any successful parameter sets. Note, we aim to ensure that networks
with perfect resilience do not require strict constraints on the parameters,
and therefore do not require an intensive parameter search for stability and
feasibility.

We lastly ensure that successful systems have perfect resilience, and not a
trivial form in which the output species has no reaction to changes in
stimulus [37]. We simulate each system and check that the output reacts to a
change in stimulus by at least $1\%$ of its pre-stimulus abundance and then
returns to within $1\%$ of its pre-stimulus abundance at steady state. For
networks which pass all of the above tests, we generate $2000$ parameter sets
that enable feasible and stable steady states, and use these for further
analysis. See Supplementary Figure 3 for a graphical representation of this
process and the number of networks which proceed after each test.

## Results

In the following sections we identify all network configurations which have
the perfect resilience property under the generalised Lotka-Volterra
framework. We examine the structural constraints on ecosystem networks capable
of perfect resilience, and how these constraints and the associated transient
dynamics affect the possibility of observing perfect resilience in in situ
ecosystems.

### Network Structures with Perfect Resilience

Our extensive analysis of networks required testing a total of 20,480 possible
network structures. In Supplementary Figure 1, we provide an extensive list of
all 23 networks which are capable of perfect resilience before considering the
sign of interactions between species. This is approximately $0.1\%$ of
possible networks which are capable of perfect resilience. In Figure 3 we
present the 23 networks as ten unique, condensed network motifs which
illustrate the general trends in these structures.

Figure 3: Unique network motifs which enable perfect resilience in a
generalised Lotka-Volterra three-species system. Below each motif is the
number of variations on the motif which are capable of perfect resilience when
the sign structure is accounted for - a total of 82 networks. The full list of
23 networks (without assigned sign structure) which are capable of perfect
resilience can be found in Supplementary Figure 1. The 23 networks are
condensed into ten motifs by including optionally present growth rates, where
a growth rate may be present or absent and still have perfect resilience in
the associated network. For example, (f) could be split into two networks with
or without the $r_{O}$ growth rate. The motif in (c) is an exception where
perfect resilience is present in networks in which: only $r_{O}$ is present;
$r_{O}$ and $r_{I}$ are present; and when $r_{O}$, $r_{I}$ and $r_{M}$ are
present. The motifs in (f)-(i) have non-unique variations in which species $I$
and $M$ are switched. The interaction signs are determined by 2000 parameter
sets which permitted stable, feasible steady states, as discussed in the
Methods section.

After identifying the capacity for perfect resilience, we determine random
parameter regimes which enable feasible and stable steady states for these
networks (see Methods). Some parameters, like the growth rates and
intraspecific competition constants, can only take on one sign (strictly
positive or strictly negative respectively). To ensure feasibility and
stability, some interaction and stimulus constants must take a specific sign,
however most motifs had interactions that could take on positive or negative
values (arrows with rounded heads in Figure 3(a,c,e-h,j)). In these motifs, we
studied the correlation between parameter sets (associated with feasibility
and stability) and found that there is either a high correlation between the
sign of interactions, or no correlation. This can then be used to further
classify the motifs based on the sign of interactions - the sign structure.
For example, the motif in Figure 3(a) has three constants which can take on
either sign: $a_{IM}$, $d_{I}$, and $a_{OI}$. These three constants have
strong correlation which we use to classify the motif into two possible sign
structures (Figure 4(a)). This can also be observed for the motif in Figure
3(g), but one of the three interactions is weakly correlated with the others
and therefore we can define four sign structures (Figure 4(b)). The number of
sign structures are identified at the bottom of each motif in Figure 3, and
the specific sign structures can be found in Supplementary Figure 2. In total,
when we specify the sign of every interaction, we identified 82 networks which
had the perfect resilience property. These 82 networks represent less than
$0.01\%$ of the possible network structures (when interaction sign is
included).

Figure 4: Two examples of how the correlation between parameters can be used
to specify the sign structure of networks capable of perfect resilience. The
motifs in (a) and (b) are from Figure 3(a) and (g) respectively. The
correlation of parameter sets is represented in the grid where a value of $-1$
(black) represents a strong, negative correlation, and a value of $1$ (white)
represents a strong, positive correlation. Strong correlations between
parameters indicate a dependent relationship between the sign of those
parameters, whereas a weak correlation (around $0$) indicates independence.
These can then be used to create the associated tables of possible sign
structure combinations.

While further classifying motifs based on sign structure can be used to
provide more detailed motifs, we found that this did not provide
distinguishable insights into the structure or dynamics of the ten unique
network motifs in Figure 3 without the sign structure fully specified. In the
following section, we draw conclusions on the structural requirements of these
motifs and relate these to conditions we would expect from real ecosystems.

### Topological Constraints and Ecosystem Implications

In this extensive study, we examined the five possible configurations for how
the stimulus can interact with populations in the network (excluding the sign
and strength of those interactions), however only three were found to have
motifs capable of perfect resilience when the stimulus affects: $I$, $O$, and
$I$ and $O$ (Figure 3, red arrows). Crucially there was no capacity for
perfect resilience if the stimulus affects $I$ and $M$, or when it affects all
three species. This is particularly important in the latter case as a stimulus
representing environmental events, such as heatwaves or cyclones, are likely
to directly affect all species in the network. This indicates that perfect
resilience is never possible under these crucial types of stimuli which are
increasing in frequency with climate change [54].

When we translate the perfect resilience motifs into the mechanisms which
obtain robust perfect adaptation in biochemical reactions, we identified that
our motifs only used opposer mechanisms to generate perfect resilience - no
balancer mechanisms. Araujo and Liotta [42] determined that networks cannot
have balancer mechanisms if the stimulus directly affects the output, because
it requires multiple paths connecting the stimulus to the output to
effectively ‘balance out’ a change in stimulus. We do not observe multiple
paths in any of the possible networks where the stimulus does not directly
affect the output, $O$ (Figure 3(a)-(e)).

The structure of all ten motifs are dependent on a sequence of one-way
interactions between species (Figure 3, dark blue arrows). For example, the
network in Figure 3(a) has a one-way interaction connecting $M$ to $I$
followed by $I$ to $O$. In ecological systems having a one-way interaction is
possible, an example being that an orchid on a tree benefits from the tree,
but the tree has no significant benefit or harm from the orchid [55]. However,
having multiple one-way interactions, in sequence, is highly unlikely in real
ecological systems.

Lastly, in all of the motifs we observe either intrinsic growth or self-
regulation terms in only a subset of the species i.e. not for every species.
It is not uncommon to exclude intrinsic growth terms for some species since,
by definition, these terms are included to represent implicit increases for
that population. However, self-regulation terms are more often included for
all species as this ensures a carrying capacity for that species [56, 57, 58].
In our perfect resilience motifs, only one species in the network ever has a
self-regulation term (Figure 3, light blue looping arrow), resulting in the
potential for unbounded growth for two species. Self-regulation terms also
play an important role in the stability of the steady state, and dampening
oscillations of a species following a perturbation [56, 57, 58]. In the
following section we explore the dynamics of these networks and any
implications that these have on perfect resilience.

### Large Oscillations and Transient Dynamic Limitations

In all ten motifs we observed a particular variety of integral control -
constrained integral control [44, 59] \- in which perfect resilience is
present in the output species, $O$, as long as it does not go extinct. This
occurs when there is an isolated factor of $O$ in the projection (Equation
(4)), which can be observed in the example in Equation (5). While this is
realistic for an ecosystem, it is a stark contrast to chemical reaction
networks where a molecule can be created from reactions independent of its
abundance. Ecosystem networks are therefore limited by the strength of
perturbations which the output is able to recover from while avoiding
extinction.

Generalised Lotka-Volterra models are prone to generating transient dynamics
which are highly oscillatory [57, 58] and the networks which we identified as
having perfect resilience are not exempt from this behaviour (Figure 5(a)).
When the stimulus changes, these networks rapidly oscillate in response and
can take significant time to return to their pre-stimulus abundance. If
another perturbation occurs within this oscillatory period it is possible that
the abundance of the species can be perturbed to extinction and perfect
resilience will be lost (Figure 5(b)). Moreover, in reality having
oscillations which repeatedly bring the output close to extinction have a
higher risk of a stochastic event killing off the population.

Figure 5: Two examples of the oscillatory behaviours present in ecosystem
networks capable of perfect resilience. In (a) we demonstrate a representative
example of the highly oscillatory behaviour from the motif in Figure 3(g)
which is present in all motifs. Note that the oscillations have an increasing
frequency and amplitude with repeated perturbations. In (b) we demonstrate how
repeated perturbations before the network has returned to steady state results
in extinction (motif from Figure 3(h)). In both cases we simulate the system
for one hundred time steps before doubling the stimulus, specifically
$S=\\{1,2,4,8\\}$. Parameters: (a) $r_{I}=0.6$, $r_{M}=0.46$, $r_{O}=0.5$,
$d_{O}=-0.74$, $a_{13}=-0.02$, $a_{21}=0.46$, $a_{22}=-1$, $a_{31}=0.85$,
$a_{32}=-0.18$, (b) $r_{M}=0.04$, $r_{O}=0.09$, $d_{O}=-0.01$, $a_{12}=-0.18$,
$a_{13}=0.83$, $a_{22}=-0.26$, $a_{23}=0.64$, $a_{31}=-0.42$.

## Discussion

In this study we present a novel attempt to find perfect resilience in
ecological systems. Moreover we have developed fully automated methods for
identifying perfect resilience which are widely applicable to networks with
any number of species, or functional groups, and under alternate modelling
frameworks. We have identified that there is potential for ecological systems
to obtain perfect resilience - albeit in only 23 network structures out of a
possible 20,480 configurations in a three-species system. The motifs that we
have identified can provide an important insight into ecological management
decisions, particularly since this only requires an understanding of the
structure of a network - an easier task compared to determining the
interaction strengths in a network. By understanding the motif structure, it
can then assist in highlighting interventions to ensure perfect resilience
either does or does not occur for a target species or ecosystem function [1,
2, 3].

The 23 networks that we identified as being capable of perfect resilience were
all based on opposer mechanisms [42], instead of the mixture of opposers and
balancers which can promote robust perfect adaptation in biochemical reaction
networks. This matches observations in biochemical literature, since the
oscillatory behaviours of our networks are only possible under the feedback
structures found in opposer mechanisms [43]. In biochemical networks there is
a significantly better understanding of balancer mechanisms [43, 26, 60, 61,
62, 63], compared to opposer mechanisms [38, 41, 40, 39, 43, 64, 45, 37, 46,
51, 47, 48, 49, 65, 50]. There are two key drivers underlying opposer
mechanisms - ultrasensitivity and antithetic integral control [43, 45, 37, 46,
47, 48, 49, 65, 50, 51, 38, 39, 40, 41]. These are dependent, respectively: on
molecules rapidly converting to an active (‘on’) or inactive (‘off’) state,
like a switch, to regulate an output; or the annihilation of excess output
molecules, which get transformed into a combined (complexed) form before being
removed. Neither of these drivers are examples of realistic ecosystem
behaviours. Moreover, both drivers utilise the formation of enzyme-substrate
complexes, but there is no equivalent to this set of reactions in species
interactions.

Many of the distinctions between ecological systems and biochemical reaction
networks arise from structural differences between the two types of network. A
clear feature in all ten network motifs, and likely pivotal to the opposer
mechanism, is the presence of a sequence of one-way interactions between
species. One-way interactions, where only one population is affected by an
interaction, can be found in ecological networks, such as epiphytes living on
trees [55], but observing multiple one-way interactions in sequence is
unlikely. Chemical reaction networks can include one-way interactions,
dependent on the modelling framework used to describe the behaviour. For
example, Michaelis-Menten equations readily display one-way interactions with
enzymes which are able to alter molecular concentrations without having any
effect on their own abundance [37]. In all of the perfect resilience motifs we
also identified that only a subset of species have intraspecific competition
terms. These terms play a crucial role in generalised Lotka-Volterra equations
to dampen oscillations, and increase stability of the system [56, 57, 58]. We
observe that all of our motifs are associated with highly oscillatory
behaviours which make them more susceptible to successive perturbations and at
a higher risk of extinction from stochastic disturbances. Similar behaviours
have been observed in coral reefs which are impacted by successive
disturbances [66]. Another crucial finding of this work is that there were no
networks in which changes that affect all populations directly, can possibly
have the perfect resilience property. This implies that disturbances,
particularly environmental disturbances such as heatwaves, cannot be moderated
by perfect resilience.

In this study, we have focused primarily on the commonly used generalised
Lotka-Volterra equations and a stimulus which represents the direct effect of
policy and environmental changes. An important problem with these equations
are parameter identifiability [17, 18]. Generalised Lotka-Volterra equations
simplify complex dimensions such as stochasticity, spatial and genetic
heterogeneity, and complex interactions (e.g. prey switching). This makes it
challenging to measure, or infer, the net per capita strength of an
interaction between any two species. While perfect resilience overcomes this
problem by treating parameters as undefined variables, a more pressing concern
is that there can be difficulties in even identifying the sign of the net
interaction between species [67]. This poses a big problem with identifying
these motifs in real ecosystems. While these problems are likely not limited
to generalised Lotka-Volterra equations, it is important to note that the
methods used here for identifying perfect resilience can readily be applied to
other modelling frameworks, with only some constraints required for
calculating the Gröbner basis [44, 68].

The network motifs identified in this study have important implications for
conservation management and socioecological systems. While these motifs can be
used to identify resilient ecological networks, a pivotal problem that
managers will face is attempting to design, or alter, ecosystems. In
biochemical reactions identifying robust perfect adaptation can be used to
suggest targeted treatments which can later undergo in vitro or clinical
testing. However, for managers, ecosystem engineering faces a slew of problems
which cannot be tested before implementation. Removing species through culling
can be controversial for native species [69, 70, 71, 72, 73], and the
introduction of non-native species to an ecosystem has repeatedly had adverse
effects [74, 75, 76]. In this study, we only performed an extensive search of
three-species ecosystems. An existing study from biochemical literature [42]
identifies how smaller mechanisms can be embedded within larger networks while
ensuring that robust perfect adaptation, or perfect resilience, still obtains.
There is still a benefit in searching for perfect resilience in larger
networks to identify possible mechanisms that may only be possible with a
larger number of species (potentially balancer mechanisms), and the methods
developed here constitute an important foundation for further work in this
area. Future work may also be to identify the frequency with which these
motifs occur in situ [19].

## References

  * [1] Folke, C. _et al._ Regime shifts, resilience, and biodiversity in ecosystem management.  _Annu. Rev. Ecol. Evol. Syst._ 35, 557–581 (2004).
  * [2] May, R. M.  _Stability and complexity in model ecosystems_ (Princeton university press, 2019).
  * [3] Parrott, L. & Meyer, W. S.  Future landscapes: managing within complexity.  _Frontiers in Ecology and the Environment_ 10, 382–389 (2012).
  * [4] Mumby, P. J., Chollett, I., Bozec, Y.-M. & Wolff, N. H.  Ecological resilience, robustness and vulnerability: how do these concepts benefit ecosystem management?  _Current Opinion in Environmental Sustainability_ 7, 22–27 (2014).
  * [5] Urruty, N., Tailliez-Lefebvre, D. & Huyghe, C.  Stability, robustness, vulnerability and resilience of agricultural systems. a review.  _Agronomy for sustainable development_ 36, 1–15 (2016).
  * [6] Evans, D. M., Pocock, M. J. & Memmott, J.  The robustness of a network of ecological networks to habitat loss.  _Ecology letters_ 16, 844–852 (2013).
  * [7] DiGiano, M. L. & Racelis, A. E.  Robustness, adaptation and innovation: Forest communities in the wake of hurricane dean.  _Applied Geography_ 33, 151–158 (2012).
  * [8] Lake, P. S.  Resistance, resilience and restoration.  _Ecological management & restoration_ 14, 20–24 (2013).
  * [9] Hoover, D. L., Knapp, A. K. & Smith, M. D.  Resistance and resilience of a grassland ecosystem to climate extremes.  _Ecology_ 95, 2646–2656 (2014).
  * [10] MacGillivray, C., Grime, J. & Team, T. I. S. P. I.  Testing predictions of the resistance and resilience of vegetation subjected to extreme events.  _Functional Ecology_ 640–649 (1995).
  * [11] Holling, C. S.  Engineering resilience versus ecological resilience.  _Engineering within ecological constraints_ 31, 32 (1996).
  * [12] Carpenter, S., Walker, B., Anderies, J. M. & Abel, N.  From metaphor to measurement: resilience of what to what?  _Ecosystems_ 4, 765–781 (2001).
  * [13] Lavorel, S.  Ecological diversity and resilience of mediterranean vegetation to disturbance.  _Diversity and distributions_ 5, 3–13 (1999).
  * [14] Holling, C. S.  Resilience and stability of ecological systems.  _Annual review of ecology and systematics_ 4, 1–23 (1973).
  * [15] Neubert, M. G. & Caswell, H.  Alternatives to resilience for measuring the responses of ecological systems to perturbations.  _Ecology_ 78, 653–665 (1997).
  * [16] Meyer, K.  A mathematical review of resilience in ecology.  _Natural Resource Modeling_ 29, 339–352 (2016).
  * [17] Adams, M. P. _et al._ Informing management decisions for ecological networks, using dynamic models calibrated to noisy time-series data.  _Ecology Letters_ 23, 607–619 (2020).
  * [18] Remien, C. H., Eckwright, M. J. & Ridenhour, B. J.  Structural identifiability of the generalized lotka–volterra model for microbiome studies.  _Royal Society Open Science_ 8, 201378 (2021).
  * [19] Barnes, M. L. _et al._ The social structural foundations of adaptation and transformation in social–ecological systems.  _Ecology and Society_ 22 (2017).
  * [20] Janssen, M. A. _et al._ Toward a network perspective of the study of resilience in social-ecological systems.  _Ecology and society_ 11 (2006).
  * [21] Rao, C. V., Kirby, J. R. & Arkin, A. P.  Design and diversity in bacterial chemotaxis: a comparative study in escherichia coli and bacillus subtilis.  _PLoS biology_ 2, e49 (2004).
  * [22] Parent, C. A. & Devreotes, P. N.  A cell’s sense of direction.  _Science_ 284, 765–770 (1999).
  * [23] Macnab, R. M. & Koshland Jr, D. E.  The gradient-sensing mechanism in bacterial chemotaxis.  _Proceedings of the National Academy of Sciences_ 69, 2509–2512 (1972).
  * [24] Levchenko, A. & Iglesias, P. A.  Models of eukaryotic gradient sensing: application to chemotaxis of amoebae and neutrophils.  _Biophysical journal_ 82, 50–63 (2002).
  * [25] Kollmann, M., Løvdok, L., Bartholomé, K., Timmer, J. & Sourjik, V.  Design principles of a bacterial signalling network.  _Nature_ 438, 504–507 (2005).
  * [26] Goy, M. F., Springer, M. S. & Adler, J.  Sensory transduction in escherichia coli: role of a protein methylation reaction in sensory adaptation.  _Proceedings of the National Academy of Sciences_ 74, 4964–4968 (1977).
  * [27] Berg, H. C. & Brown, D. A.  Chemotaxis in escherichia coli analysed by three-dimensional tracking.  _Nature_ 239, 500–504 (1972).
  * [28] Alon, U., Surette, M. G., Barkai, N. & Leibler, S.  Robustness in bacterial chemotaxis.  _Nature_ 397, 168–171 (1999).
  * [29] Yi, T.-M., Huang, Y., Simon, M. I. & Doyle, J.  Robust perfect adaptation in bacterial chemotaxis through integral feedback control.  _Proceedings of the National Academy of Sciences_ 97, 4649–4653 (2000).
  * [30] Yau, K.-W. & Hardie, R. C.  Phototransduction motifs and variations.  _Cell_ 139, 246–264 (2009).
  * [31] Reisert, J. & Matthews, H. R.  Response properties of isolated mouse olfactory receptor cells.  _The Journal of physiology_ 530, 113–122 (2001).
  * [32] Matthews, H. R. & Reisert, J.  Calcium, the two-faced messenger of olfactory transduction and adaptation.  _Current opinion in neurobiology_ 13, 469–475 (2003).
  * [33] Kaupp, U. B.  Olfactory signalling in vertebrates and insects: differences and commonalities.  _Nature Reviews Neuroscience_ 11, 188–200 (2010).
  * [34] Medina-Gomez, G. _et al._ Adaptation and failure of pancreatic $\beta$ cells in murine models with different degrees of metabolic syndrome.  _Disease models & mechanisms_ 2, 582–592 (2009).
  * [35] Araujo, R. P., Petricoin, E. F. & Liotta, L. A.  Mathematical modeling of the cancer cell’s control circuitry: paving the way to individualized therapeutic strategies.  _Current Signal Transduction Therapy_ 2, 145–155 (2007).
  * [36] Fodale, V., Pierobon, M., Liotta, L. & Petricoin, E.  Mechanism of cell adaptation: when and how do cancer cells develop chemoresistance?  _Cancer journal (Sudbury, Mass.)_ 17, 89 (2011).
  * [37] Ma, W., Trusina, A., El-Samad, H., Lim, W. A. & Tang, C.  Defining network topologies that can achieve biochemical adaptation.  _Cell_ 138, 760–773 (2009).
  * [38] Frei, T., Chang, C.-H., Filo, M., Arampatzis, A. & Khammash, M.  A genetic mammalian proportional–integral feedback control circuit for robust and precise gene regulation.  _Proceedings of the National Academy of Sciences_ 119, e2122132119 (2022).
  * [39] Aoki, S. K. _et al._ A universal biomolecular integral feedback controller for robust perfect adaptation.  _Nature_ 570, 533–537 (2019).
  * [40] Briat, C., Gupta, A. & Khammash, M.  Antithetic integral feedback ensures robust perfect adaptation in noisy biomolecular networks.  _Cell systems_ 2, 15–26 (2016).
  * [41] Khammash, M. H.  Perfect adaptation in biology.  _Cell Systems_ 12, 509–521 (2021).
  * [42] Araujo, R. P. & Liotta, L. A.  The topological requirements for robust perfect adaptation in networks of any size.  _Nature communications_ 9, 1757 (2018).
  * [43] Araujo, R. P. & Liotta, L. A.  Design principles underlying robust adaptation of complex biochemical networks.  In _Computational Modeling of Signaling Networks_ , 3–32 (Springer, 2023).
  * [44] Araujo, R. P. & Liotta, L. A.  Universal structures for adaptation in biochemical reaction networks.  _Nature Communications_ 14, 2251 (2023).
  * [45] Jeynes-Smith, C. & Araujo, R.  Protein–protein complexes can undermine ultrasensitivity-dependent biological adaptation.  _Journal of the Royal Society Interface_ 20, 20220553 (2023).
  * [46] Barkai, N. & Leibler, S.  Robustness in simple biochemical networks.  _Nature_ 387, 913–917 (1997).
  * [47] Sontag, E. D.  Adaptation and regulation with signal detection implies internal model.  _Systems & control letters_ 50, 119–126 (2003).
  * [48] Ang, J., Bagh, S., Ingalls, B. P. & McMillen, D. R.  Considerations for using integral feedback control to construct a perfectly adapting synthetic gene network.  _Journal of theoretical biology_ 266, 723–738 (2010).
  * [49] Ang, J. & McMillen, D. R.  Physical constraints on biological integral control design for homeostasis and sensory adaptation.  _Biophysical journal_ 104, 505–515 (2013).
  * [50] Ferrell, J. E.  Perfect and near-perfect adaptation in cell signaling.  _Cell systems_ 2, 62–67 (2016).
  * [51] Koshland Jr, D. E., Goldbeter, A. & Stock, J. B.  Amplification and adaptation in regulatory and sensory systems.  _Science_ 217, 220–225 (1982).
  * [52] Zheng, D. W., Bengtsson, J. & Agren, G. I.  Soil food webs and ecosystem processes: decomposition in donor-control and lotka-volterra systems.  _The American Naturalist_ 149, 125–148 (1997).
  * [53] Lyapunov, A. M.  The general problem of the stability of motion.  _International journal of control_ 55, 531–534 (1992).
  * [54] Stott, P.  How climate change affects extreme weather events.  _Science_ 352, 1517–1518 (2016).
  * [55] Rasmussen, H. N. & Rasmussen, F. N.  The epiphytic habitat on a living host: reflections on the orchid–tree relationship.  _Botanical Journal of the Linnean Society_ 186, 456–472 (2018).
  * [56] Hening, A. & Nguyen, D. H.  Persistence in stochastic lotka–volterra food chains with intraspecific competition.  _Bulletin of mathematical biology_ 80, 2527–2560 (2018).
  * [57] Wangersky, P. J.  Lotka-volterra population models.  _Annual Review of Ecology and Systematics_ 9, 189–218 (1978).
  * [58] Anisiu, M.-C.  Lotka, volterra and their model.  _Didáctica mathematica_ 32 (2014).
  * [59] Xiao, F. & Doyle, J. C.  Robust perfect adaptation in biomolecular reaction networks.  In _2018 IEEE conference on decision and control (CDC)_ , 4345–4352 (IEEE, 2018).
  * [60] Mangan, S. & Alon, U.  Structure and function of the feed-forward loop network motif.  _Proceedings of the National Academy of Sciences_ 100, 11980–11985 (2003).
  * [61] Goentoro, L., Shoval, O., Kirschner, M. W. & Alon, U.  The incoherent feedforward loop can provide fold-change detection in gene regulation.  _Molecular cell_ 36, 894–899 (2009).
  * [62] Skataric, M., Nikolaev, E. V. & Sontag, E. D.  Fundamental limitation of the instantaneous approximation in fold-change detection models.  _IET systems biology_ 9, 1–15 (2015).
  * [63] Shinar, G. & Feinberg, M.  Structural sources of robustness in biochemical reaction networks.  _Science_ 327, 1389–1391 (2010).
  * [64] Jeynes-Smith, C. & Araujo, R. P.  Ultrasensitivity and bistability in covalent-modification cycles with positive autoregulation.  _Proceedings of the Royal Society A_ 477, 20210069 (2021).
  * [65] Drengstig, T. _et al._ A basic set of homeostatic controller motifs.  _Biophysical journal_ 103, 2000–2010 (2012).
  * [66] Ortiz, J.-C. _et al._ Impaired recovery of the great barrier reef under cumulative stress.  _Science advances_ 4, eaar6127 (2018).
  * [67] Peterson, K. A. _et al._ Reconstructing lost ecosystems: A risk analysis framework for planning multispecies reintroductions under severe uncertainty.  _Journal of Applied Ecology_ 58, 2171–2184 (2021).
  * [68] Buchberger, B.  Gröbner bases: A short introduction for systems theorists.  In _Computer Aided Systems Theory—EUROCAST 2001: A Selection of Papers from the 8th International Workshop on Computer Aided Systems Theory Las Palmas de Gran Canaria, Spain, February 19–23, 2001 Revised Papers_ , 1–19 (Springer, 2002).
  * [69] Read, J. L. _et al._ Improving kangaroo management: a joint statement (2021).
  * [70] Drijfhout, M., Kendal, D. & Green, P.  Mind the gap: Comparing expert and public opinions on managing overabundant koalas.  _Journal of Environmental Management_ 308, 114621 (2022).
  * [71] Drijfhout, M., Kendal, D. & Green, P. T.  Understanding the human dimensions of managing overabundant charismatic wildlife in australia.  _Biological Conservation_ 244, 108506 (2020).
  * [72] McKinnon, M. _et al._ Media coverage of lethal control: A case study of kangaroo culling in the australian capital territory.  _Human Dimensions of Wildlife_ 23, 90–99 (2018).
  * [73] Mehmet, M. & Simmons, P.  Kangaroo court? an analysis of social media justifications for attitudes to culling.  _Environmental Communication_ 12, 370–386 (2018).
  * [74] Fabres, G. & Brown Jr, W. L.  The recent introduction of the pest ant wasmannia auropunctata into new caledonia.  _Australian Journal of Entomology_ 17, 139–142 (1978).
  * [75] Nugent, G., Fraser, W. & Sweetapple, P.  Top down or bottom up? comparing the impacts of introduced arboreal possums and ‘terrestrial’ruminants on native forests in new zealand.  _Biological Conservation_ 99, 65–79 (2001).
  * [76] Shine, R., Ward-Fear, G. & Brown, G. P.  A famous failure: Why were cane toads an ineffective biocontrol in australia?  _Conservation Science and Practice_ 2, e296 (2020).


# Reinforcement Learning Trajectory Generation and Control for Aggressive
Perching on Vertical Walls with Quadrotors

Chen-Huan Pi, Kai-Chun Hu, Yu-Ting Huang, Stone Cheng

(October 2019)

###### Abstract

Micro aerial vehicles are widely being researched and employed due to their
relative low operation costs and high flexibility in various applications. We
study the under-actuated quadrotor perching problem, designing a trajectory
planner and controller which generates feasible trajectories and drives
quadrotors to desired state in state space. This paper proposes a trajectory
generating and tracking method for quadrotor perching that takes the
advantages of reinforcement learning controller and traditional controller.
The trained low-level reinforcement learning controller would manipulate
quadrotor toward the perching point in simulation environment. Once the
simulated quadrotor has successfully perched, the relative trajectory
information in simulation will be sent to tracking controller on real
quadrotor and start the actual perching task. Generating feasible trajectories
via the trained reinforcement learning controller requires less time, and the
traditional trajectory tracking controller could easily be modified to control
the quadrotor and mathematically analysis its stability and robustness. We
show that this approach permits the control structure of trajectories and
controllers enabling such aggressive maneuvers perching on vertical surfaces
with high precision.

###### Index Terms:

Aggressive flight, reinforcement learning, quadrotor

##  I Introduction

Micro Aircraft Vehicles (MAVs) are widely being researched and employed due to
their relative low operation costs and high flexibility in various
applications. The applications include search and rescue [2], wind turbine
inspection [13, 8], reconnaissance and mapping [3], target tracking,
environmental monitoring, etc.

Most of the quadrotor is powered by battery for motor driving and generate
lift, which gives a limited time for quadrotor tasks performing. Flying with a
power cable connected to power supply on ground [5, 1, 9] can solve the
problem but the quadrotor can only be operated in a specific and relative open
area. If the quadrotor can perch on arbitrary angle surface such as vertical
wall to reduce its battery consumption while keeping environment monitoring
using on-board camera or even for battery charging. The mission time can be
extended with a great amount of time.

The quadrotor is an under-actuated system, meaning it cannot control its
position while maintaining arbitrary attitude. To successfully perch on the
wall, quadrotors using gripper mechanism [16, 18, 7] to perch on vertical
walls had been proposed. However, the perching device would cost additional
energy consumption. Another approach for perching is using aggressive maneuver
to land without adding mechanism on quadrotor.

There are a few other researches discussing perching of traditional
quadrotors. To land a traditional quadrotor, the process can be separated into
two different parts, 1) Trajectory planning to the target landing point. 2)
Track the designed trajectory and land. To generate a feasible trajectory,
Quadratic Program problem were formulated under the dynamics constrain of
quadrotor [17, 11].

Figure 1: Aggressive perching maneuver using reinforcement learning trajectory
planner and tracking controller

Using Reinforcement Learning (RL) for low-level quadrotor control has been
proposed and shown its implementation result in real world quadrotor to
stabilize and hover or fly on a specific trajectory [6, 12]. The advantage of
using RL controller is due to its nonlinearity and fast computation. However,
there exist a gap showing on the tracking result while transferring the RL
controller on quadrotors from simulation to real world experiment. In reality
condition, the quadrotor hovering and tracking results didn’t perform as same
as the results in simulation. The reason is that the RL-based controller is
sensitive to the deviation of model physical parameters such as take-off
weight or motor characteristic. This situation can be harsher when considering
aggressive maneuver using RL controlling the motors directly.

To fix this problem, we propose a method that takes the advantages of RL
controller and traditional controller illustrated in Fig. 1. First, we train a
RL controller to learn how to perch the quadrotor to the set point by directly
command the thrust force of each motor. Second, we create a simulated
quadrotor which has the same initial states as the real quadrotor. The trained
low-level RL controller would manipulate quadrotor toward the landing point in
simulation environment. Third, we extract the trajectory information in
simulation and send it to traditional controller on real quadrotor. Generating
feasible trajectories via the trained RL controller would require less time,
and the traditional trajectory tracking controller could easily be modified to
control the quadrotor and mathematically analysis its stability and
robustness. Moreover, this combined structure of controller can be implement
to different physical parameters or different type of vehicle.

Figure 2: The quadrotor launch in random initial state. Stage I generate the
trajectory of perching. In Stage II, the quadrotor start to track along the
planned trajectory and switch to Stage III for attitude control while it
approach to landing position.

In Section II, the dynamics of quadrotor and reinforcement learning algorithm
were presented. Section III and Section IV describes the RL perching training
for trajectory generating and the traditional controller design to perform
perching task. Section V shows the successful application of the combined
structure to the quadrotor, and Section VI conclude the paper.

##  II Methodology

The first challenge of planning is that since the quadrotor is an under-
actuated system, the designed trajectory needs to be feasible with respect to
the system dynamic constrains. A possible technique is to solve a constrained
optimization, for example:

$\begin{split}\min_{a}&\int_{t_{i}}^{t_{f}}R(x,\dot{x},a)dt+T(x(t_{f}))\\\
&\dot{x}=f(x,a)\\\ &x(t_{i})=x_{0}\end{split}$ (1)

where $R$ is run time cost, $T$ terminal cost, $f$ system dynamic, $a$ input
action, and $x_{0}$ is the initial state. With this scheme, we can obtain a
suitable trajectory $x(t)$ with respect to time.

In this paper, we assuming dynamic trajectory planning problem can be
described by a Markov decision processes (MDP), with this property the best
decisions are independent to the history. Next, we utilize reinforcement
learning to search the feasible solutions in (1) while the searching results
were memorized by an artificial neural network. With this scheme, the
searching cost passes through the training phase.

###  II-A Reinforcement Learning

Reinforcement learning (RL) is one of basic machine learning paradigms, aim to
solve MDP through learning processes. In this section, we provide a short
introduction to reinforcement learning.

A Markov decision processes is described by a 4-tuple
$(\mathcal{S},\mathcal{A},\mathcal{P},r)$ where $\mathcal{S}$ is the state
space, $\mathcal{A}$ is action space,
$\mathcal{P}=\mathcal{P}(s_{t+1}|s_{t}=s\in\mathcal{S},a_{t}=a\in\mathcal{A})$
is transition probability and $r=r(s,a)\in\mathbb{R}$ is the reward function.
The goal of MDPs is to find a policy

$\pi=\pi(a\in\mathcal{A}|s)$ (2)

that maximize cumulative discounted reward

$\sum_{t=0}^{\infty}\gamma^{t}r(s_{t},a_{t}),$ (3)

where $\gamma\in(0,1)$ is discounting factor. To connect (1), we write MDP
statement as follow:

$\begin{split}\max_{\pi}\mathbb{E}&\left[\sum_{t=0}^{\infty}\gamma^{t}r_{t}\Big{|}a_{t}\sim\pi(a|s_{t})\right]\\\
&s_{t+1}\sim P(s|s_{t},a_{t})\\\ &s_{0}\sim P_{0}(s)\end{split}$ (4)

where $r_{t}=r(s_{t},a_{t})$, $P_{0}$ is the initial state distribution. For a
given policy $\pi$, the following recurrent relation hold

$V^{\pi}(s_{t})=\mathbb{E}\left[r_{t}+\gamma
V^{\pi}(s_{t+1})\Big{|}a_{t}\sim\pi(a|s_{t})\right],$ (5)

where

$V^{\pi}(s)=\mathbb{E}\left[\sum_{t=0}^{\infty}\gamma^{t}r_{t}\Big{|}a_{t}\sim\pi(a|s_{t}),s_{0}=s\right].$
(6)

The quantity

$A^{\pi}(x,a)=\mathbb{E}\left[r_{t}+\gamma
V^{\pi}(s_{t+1})|s_{t}=x,a_{t}=a\right]-V^{\pi}(s)$ (7)

can be used for determine that which action is good with respect to policy
$\pi$. In fact, we can guarantee [12] that a policy $\pi$ is better than
policy $\mu$ if

$\forall(s,a)\in\mathcal{S}\times\mathcal{A}$
$\left[\pi(a|s)-\mu(a|s)\right]A^{\mu}(s,a)>0$ (8a) or
$\left[\pi(a|s)-\mu(a|s)\right]A^{\pi}(s,a)>0.$ (8b)

Therefore, for a policy search iteration, Equation (8) provides an improvement
guideline to distinguish whatever an action should be memorized or not.

###  II-B Quadrotor Dynamics

The simplest way to modeling a multi-rotor is considered as a rigid body
dynamics

$\frac{d}{dt}\begin{pmatrix}\mathbf{R}\\\ \mathbf{\Omega}\\\ \mathbf{x}\\\
\mathbf{v}\\\
\mathbf{a}\end{pmatrix}=\begin{pmatrix}\mathbf{R}\hat{\mathbf{\Omega}}\\\
\mathbf{J}^{-1}\left(\tau-\mathbf{\Omega}\times\mathbf{J}\mathbf{\Omega}\right)\\\
\mathbf{v}\\\
\frac{1}{m}\mathbf{R}\left[0~{}0~{}T_{z}\right]^{T}+\mathbf{g}\\\
\mathbf{s}\end{pmatrix}$ (9)

where $\mathbf{R}$ is rotation matrix, $\mathbf{\Omega}$ is angular velocity
on body frame, and the hat map $\hat{}:\mathbb{R}^{3}\rightarrow SO(3)$ is
defined as the condition $\hat{x}y=x\times y,\forall x,y\in\mathbb{R}^{3}$
[10]. The $\mathbf{x}$, $\mathbf{v}$, and $\mathbf{a}$ are position, velocity,
acceleration with respect to inertia frame, $J$ is the momentum of inertia,
$m$ is the mass of rigid body. $\mathbf{g}$ is the gravity acceleration. The
force and moment model of each quadrotor motor can be simplified as

$F_{i}=k_{F}\omega^{2}_{i},M_{i}=k_{M}\omega^{2}_{i},$ (10)

where $\omega_{i}$ is the rotation speed of motors and $i=1,2,3,4$. $k_{f}$,
$k_{M}$ are the coefficient of lift and moment. The relation between moment
$\tau$, force $T_{z}$ of quadrotor and motor speed $\omega_{i}$ can be written
as

$\begin{pmatrix}\mathbf{\tau}\\\
T_{z}\end{pmatrix}=\begin{pmatrix}\frac{-l}{\sqrt{2}}k_{F}&\frac{l}{\sqrt{2}}k_{F}&\frac{l}{\sqrt{2}}k_{F}&\frac{-l}{\sqrt{2}}k_{F}\\\
\frac{l}{\sqrt{2}}k_{F}&\frac{-l}{\sqrt{2}}k_{F}&\frac{l}{\sqrt{2}}k_{F}&\frac{-l}{\sqrt{2}}k_{F}\\\
k_{M}&k_{M}&-k_{M}&-k_{M}\\\
-1&-1&-1&-1\end{pmatrix}\begin{pmatrix}\omega^{2}_{1}\\\ \omega^{2}_{2}\\\
\omega^{2}_{3}\\\ \omega^{2}_{4}\end{pmatrix}$ (11)

where $l$ is the physical dimension of quadrotor.

###  II-C Training and Inference

Integrating (9), (11) and Euler method, we can obtain the discrete time
transfer function $G$ of a quadrotor

$s_{t+1}=G(s_{t},a_{t}),$ (12)

where $s_{t}=(\mathbf{R},\mathbf{\Omega},\mathbf{x},\mathbf{v})_{t}$ and
$a_{t}=(T_{1},T_{2},T_{3},T_{4})_{t}$. We change our objective in terms of
deterministic transition and the main objective become

$\begin{split}&\max_{\pi}V^{\pi}(s_{0})\\\ &s_{t+1}=F(s_{t},a_{t})\\\
&s_{0}\sim P_{0}(s)\end{split}$ (13)

in practice.

Actor-Critic architecture is used in this paper. The value network (Critic),
was trained by a V-trace method [4], which is a modified temporal difference
(TD) learning [14]. The policy network (Actor), was trained according to (8)
through hinge loss. In addition, since (8) can be satisfied by adjusting the
likelihood per pair of state and action $(s,a)$, we replace the likelihoods by
its log of likelihood for simplifying the complexity in gradient calculation.
The summarized objective is listed as follows

$\begin{split}&\min_{\theta_{\pi}}\max\\{0,\xi-\hat{A}_{t}\log\frac{\pi(a_{t}|s_{t})}{\mu(a_{t}|s_{t})}\\},\\\
&\min_{\theta_{v}}\left|\hat{V}_{t}-V(s_{t};\theta_{v})\right|^{2},\\\
&\hat{A}_{t}=r_{t}+\gamma
V(s_{t+1})-V(s_{t})+\frac{\pi(a_{t+1}|s_{t+1})}{\mu(a_{t+1}|s_{t+1})}\hat{A}_{t+1},\\\
&\hat{V}_{t}=V(s_{t};\theta_{v})+\frac{\pi(a_{t}|s_{t})}{\mu(a_{t}|s_{t})}\hat{A}_{t},\\\
&\mathcal{T}=\\{(s_{i},a_{i},s_{i+1},\mu(a_{i}|s_{i}))|i=1\cdots
T\\}\in\mathcal{B},\end{split}$ (14)

where $\xi$ is the margin to make the samples that already satisfied (8) not
to provide losses for policy improvement [12].

Notice that since the simplified rigid body dynamic is independent to the
rotor number. The different physical configurations are only effect on $J$ and
$m$, the trained RL policy could be implemented directly on to quadrotors.
Unfortunately, training RL controller on simplified ideal environment and
transfer to real world still exists system uncertainty and system response
latency. These factors may lead to the quadrotor become unstable especially
under aggressive maneuver situation. To overcome this problem, we consider a
hybrid control strategy and introduce a tracking algorithm in next section.

##  III Simulation Environment and Trajectory Generating Result

In this section, we used the aforementioned dynamic equations and
reinforcement algorithms to train a quadrotor perching policy controlling the
four motors output directly which can successfully landing on the vertical
wall. The reward function for RL training is designed as follows

$\displaystyle\text{reward}=-\begin{bmatrix}w_{1}e^{-d_{x}}&w_{2}&w_{3}\end{bmatrix}\begin{bmatrix}\parallel\mathbf{R}_{e}\parallel&\parallel\mathbf{x}_{e}\parallel&\parallel
a\parallel\end{bmatrix}^{T},$ (15)

where $\mathbf{R}_{e}$, $\mathbf{x}_{e}$, and $a$ are angle error, position
error, and action output. $w_{1}$ to $w_{3}$ are the weights of the error, and
$d_{x}$ is the distance between quadrotor to the perching point on $x$-axis.
For RL value and policy neural network construction, 3-layer with 32 and 128
nodes were used. The policy network outputs four motors thrust command to
control the quadrotor perching task.

Once the quadrotor lands successfully in the simulator controlled by RL policy
neural network, we extract the trained controller for generating feasible
trajectory. The trajectory contains position, velocity, acceleration, jerk,
and yaw angle information and use designed traditional controller to track
along the planned path.

We design a sequence for controlling the quadrotor during operating the
mission and can be separated into following three stages and shown in Fig.3:

  * •

Stage I: Trajectories generation

  * •

Stage II: Path following

  * •

Stage III: Attitude control

Figure 3: Three stages of perching control sequence

In stage I, feasible trajectory would be generated from the current state of
quadrotor and can successfully perch to desired point eventually. During the
path following stage, the controller keeps tracking on the desired path with
the information of position, velocity, acceleration, jerk, and yaw angle until
the quadrotor approaches the launch point. As the distance $d_{x}$ is shorter
than designed threshold $\epsilon$, quadrotor switch to stage III and attitude
control would be triggered. During the attitude control stage, the controller
attempt to stable the quadrotor to final pose and perch on the set point of a
vertical plane. While the quadrotor reaching the original point with pitch
angle of $90^{\circ}$, it would be considered to be a successful cases to
landing onto the vertical surface.

##  IV Controller Design

In this section, we introduce our controller for trajectory tracking along
perching neural network behavior. The tracking and perching consist of two
controllers as follows,

1\. Trajectory Tracking Control: controlling the quadrotor center of mass to
follow the three-dimensional path including position, velocity, acceleration,
and jerk information from RL trajectory generator.

2\. Attitude Control: controlling the quadrotor to desired roll, pitch and yaw
angle from trajectory tracking controller command.

###  IV-A Trajectory Tracking Control

The tracking controller is based on proportional-derivative (PD) controller
and can be written as following expression,

$\mathbf{a}_{c}=\mathbf{K_{x}}(\mathbf{x}_{ref}-\mathbf{x}_{f})+\mathbf{K_{v}}(\mathbf{v}_{ref}-\mathbf{v}_{f})+\mathbf{a}_{ref}$
(16)

where $\mathbf{a_{c}}$ is the desired acceleration vector of quadrotor,
$\mathbf{x}_{ref},\mathbf{v}_{ref},\mathbf{a}_{ref}$ are the control
references from trajectory generator, $\mathbf{x}_{f}$, and $\mathbf{v}_{f}$
are the feedback, and $\mathbf{K_{x}},\mathbf{K_{v}}$ are the gains of each
component. The total thrust $T_{c}$ and desired attitude expressed in rotation
matrix form $(\mathbf{R}_{c})_{Traj}$ can be derived as

$\displaystyle T_{c}$
$\displaystyle=-m\left\|\mathbf{a}_{c}\right\|\cdot\mathbf{b_{z}}$ (17)
$\displaystyle(\mathbf{R}_{c})_{Traj}$
$\displaystyle=\left[\mathbf{b_{x}}_{c}~{}\mathbf{b_{y}}_{c}~{}\mathbf{b_{z}}_{c}\right],$

where

$\displaystyle\mathbf{e}_{y}$
$\displaystyle=[-\sin\psi_{ref},\cos\psi_{ref},0]^{T}$ (18)
$\displaystyle\mathbf{b_{z}}_{c}$
$\displaystyle=\frac{\mathbf{a}_{c}}{\|\mathbf{a}_{c}\|}$
$\displaystyle\mathbf{b_{x}}_{c}$
$\displaystyle=\frac{\mathbf{e}_{y}\times\mathbf{b_{z}}_{c}}{\|\mathbf{e}_{y}\times\mathbf{b_{z}}_{c}\|}$
$\displaystyle\mathbf{b_{y}}_{c}$
$\displaystyle=\frac{\mathbf{b_{z}}_{c}\times\mathbf{b_{x}}_{c}}{\|\mathbf{b_{z}}_{c}\times\mathbf{b_{x}}_{c}\|}$

$\psi_{ref}$ is the desired heading angle, and
$\|\mathbf{e}_{y}\times\mathbf{b_{z}}_{c}\|>0$ is always greater than 0
because the quadrotor would not operate upside down in trajectory tracking
control.

###  IV-B Attitude Control

The attitude controller not only receive the $(\mathbf{R}_{c})_{Traj}$ command
from trajectory tracking control but also switch to final perching attitude
$(\mathbf{R}_{c})_{P}$. The controller consists of dual-loop control
structure. The outer loop P control in (21) determines the desired angular
velocity of quadrotor, where $\mathbf{R}_{c}$ switches according to the
distance of quadrotor to perching point on $x$-axis $d_{x}$.

$\mathbf{R}_{c}=\left\\{\begin{matrix}(\mathbf{R}_{c})_{Traj}&,d_{x}\geq\epsilon\\\
(\mathbf{R}_{c})_{P}&,0<d_{x}<\epsilon\end{matrix}\right.$ (19)

The error of rotational angle between $\mathbf{R}_{c}$ and measurement
feedback $\mathbf{R}_{f}$ is given by

$\displaystyle\mathbf{e_{R}}$
$\displaystyle=-0.5(\mathbf{R}_{c}^{T}\mathbf{R}_{f}-\mathbf{R}_{f}^{T}\mathbf{R}_{c})^{\vee}$
(20)

where the vee map $\vee$ is the mapping
$\mathrm{SO}(3)\rightarrow\mathbb{R}^{3}$ [10], and the desired angular
velocity $\mathbf{\Omega}_{c}$ is derived as

$\displaystyle\mathbf{\Omega}_{c}$
$\displaystyle=\mathbf{K_{R}}\mathbf{e_{R}}+\mathbf{\Omega}_{ref}.$ (21)

The $\mathbf{\Omega}_{ref}$ is the feed-forward term given from trajectory
tracking controller, which use the jerk of the trajectory to compensate
angular velocity command of the quadrotor [15].

The inner loop uses a PID controller for desired angular acceleration shown in
(22).

$\displaystyle\mathbf{e_{\Omega}}$
$\displaystyle=-\mathbf{\Omega}_{f}+\mathbf{R}_{f}^{T}\mathbf{R}_{c}\mathbf{\Omega}_{c}$
(22) $\displaystyle\dot{\mathbf{\Omega}}_{c}$
$\displaystyle=\mathbf{K_{p}}\mathbf{e_{\Omega}}+\mathbf{K_{i}}\int\mathbf{e_{\Omega}}+\mathbf{K_{d}}\mathbf{\dot{e}_{\Omega}}$

where $\mathbf{\Omega}_{f}$ is the feedback of quadrotor angular velocity.

The command of rotor speed is determined according to (9) (11) with $T_{c}$
and $\dot{\mathbf{\Omega}}_{c}$.

##  V Simulation Results

To verify that our proposed perching trajectory generating and tracking
structure could successfully landing on the wall, we conduct 50 flight tests
starting from arbitrary initial position and velocity in a $1\mathrm{m}\times
1\mathrm{m}\times 1\mathrm{m}$ space in front of the perching point in the
simulation environment illustrates in Fig. 2. All trails perch to the set
position successfully. Fig. 4 shows four flight results start from different
initial condition. Table. I shows the mean and standard deviation (SD) of the
distance error and pitch angle while the The vehicles touch the wall. The
simulation experiment shows while the quadrotor reach to $x=0$, the mean
distance error to the landing point on $y$ and $z$ is $-0.47$ cm and $-1.74$
cm. The mean pitch angle is $-1.17^{\circ}$.

TABLE I: Simulation results of sampling 50 trails and the mean contact position and SD on $y,z$-axis and pitch angle $\theta$ while the quadrotor reach position $x=0$. | Mean | SD
---|---|---
$y$-axis |  $-0.47$ cm | 0.2 cm
$z$-axis |  $-1.74$ cm | 0.21 cm
pitch angle | $88.83^{\circ}$ | $0.62^{\circ}$

(a) Case 1

(b) Case 2

(c) Case 3

(d) Case 4

Figure 4: Four cases of flight result start from different initial condition.

##  VI Conclusion

In this paper, we present a quadrotor perching control structure which
combining RL trajectory planning and traditional tracking controller. The
complex nonlinear dynamic constrain optimization problem can be solved by
using trained RL trajectory planner to find a feasible path rapidly. A well
designed traditional controller can guarantee the stability of quadrotor. The
control structure can make up the uncertainty using RL control policy output
directly to motors when it has disturbance in operation or modeling
imperfection in RL training process. The perching simulation in Section V
demonstrates the quadrotor follows the trajectory generated by RL control
policy and perches to designed point with standard deviation $0.2$ cm on
$y,z$-axis and pitch angle $88.83^{\circ}$ in average on the wall.

Future work will focus on the real-world implementation, verify the overall
control structure perching performance, and observing the robustness while
facing disturbance and model uncertainties. However, the long term goal is to
perform tasks physical contact such as battery charging or environment
monitoring.

## References

  * [1] Michele Bolognini and L. Fagiano.  Lidar-based navigation of tethered drone formations in an unknown environment.  ArXiv, abs/2003.12981, 2020.
  * [2] K. Choutri, M. Lagha, and L. Dala.  A fully autonomous search and rescue system using quadrotor uav.  2020\.
  * [3] P. Durdevic, D. Ortiz-Arroyo, and Z. Yang.  Lidar assisted camera inspection of wind turbines: Experimental study.  2019 1st International Conference on Electrical, Control and Instrumentation Engineering (ICECIE), pages 1–7, 2019.
  * [4] Lasse Espeholt, Hubert Soyer, R. Munos, K. Simonyan, V. Mnih, Tom Ward, Yotam Doron, Vlad Firoiu, T. Harley, Iain Dunning, S. Legg, and K. Kavukcuoglu.  Impala: Scalable distributed deep-rl with importance weighted actor-learner architectures.  ArXiv, abs/1802.01561, 2018.
  * [5] L. Fagiano.  Systems of tethered multicopters: Modeling and control design.  IFAC-PapersOnLine, 50:4610–4615, 2017.
  * [6] J. Hwangbo, I. Sa, R. Siegwart, and M. Hutter.  Control of a quadrotor with reinforcement learning.  IEEE Robotics and Automation Letters, 2:2096–2103, 2017.
  * [7] A. Kalantari, Karan Mahajan, D. Ruffatto, and M. Spenko.  Autonomous perching and take-off on vertical walls for a quadrotor micro air vehicle.  2015 IEEE International Conference on Robotics and Automation (ICRA), pages 4669–4674, 2015.
  * [8] Evan Kaufman, Kuya Takami, Zhuming Ai, and Taeyoung Lee.  Autonomous quadrotor 3d mapping and exploration using exact occupancy probabilities.  2018 Second IEEE International Conference on Robotic Computing (IRC), pages 49–55, 2018.
  * [9] Taeyoung Lee.  Geometric controls for a tethered quadrotor uav.  2015 54th IEEE Conference on Decision and Control (CDC), pages 2749–2754, 2015.
  * [10] Taeyoung Lee, M. Leok, and N. McClamroch.  Geometric tracking control of a quadrotor uav on se(3).  49th IEEE Conference on Decision and Control (CDC), pages 5420–5425, 2010.
  * [11] Daniel Mellinger, Nathan Michael, and Vijay Kumar.  Trajectory generation and control for precise aggressive maneuvers with quadrotors.  The International Journal of Robotics Research, 31:664 – 674, 2012\.
  * [12] Chen-Huan Pi, Kai-Chun Hu, Stone Cheng, and I-Chen Wu.  Low-level autonomous control and tracking of quadrotor using reinforcement learning.  Control Engineering Practice, 95:104222, 2020.
  * [13] Bjorn E. Schafer, D. Picchi, T. Engelhardt, and D. Abel.  Multicopter unmanned aerial vehicle for automated inspection of wind turbines.  2016 24th Mediterranean Conference on Control and Automation (MED), pages 244–249, 2016.
  * [14] R. Sutton.  Learning to predict by the methods of temporal difference learning.  1988\.
  * [15] Ezra Tal and Sertac Karaman.  Accurate tracking of aggressive quadrotor trajectories using incremental nonlinear dynamic inversion and differential flatness.  2018 IEEE Conference on Decision and Control (CDC), pages 4282–4288, 2018.
  * [16] Justin Thomas.  Grasping, perching, and visual servoing for micro aerial vehicles.  2017\.
  * [17] Justin Thomas, M. Pope, Giuseppe Loianno, E. Hawkes, M. A. Estrada, Hao Jiang, M. Cutkosky, and V. Kumar.  Aggressive flight with quadrotors for perching on inclined surfaces.  Journal of Mechanisms and Robotics, 8:051007, 2016.
  * [18] Haijie Zhang, Jiefeng Sun, and Jianguo Zhao.  Compliant bistable gripper for aerial perching and grasping.  2019 International Conference on Robotics and Automation (ICRA), pages 1248–1253, 2019.