# Token-Modification Adversarial Attacks for Natural Language Processing: A
Survey

Tom Roth1,2 Yansong Gao2 Alsharif Abuadbba2 Surya Nepal2 &Wei Liu1 1University
of Technology Sydney, Australia
2CSIRO’s Data61, Sydney, Australia {thomas.roth<EMAIL_ADDRESS>{garrison.gao, sharif.abuadbba<EMAIL_ADDRESS>
###### Abstract

There are now many adversarial attacks for natural language processing
systems. Of these, a vast majority achieve success by modifying individual
document tokens, which we call here a token-modification attack. Each token-
modification attack is defined by a specific combination of fundamental
components, such as a constraint on the adversary or a particular search
algorithm. Motivated by this observation, we survey existing token-
modification attacks and extract the components of each. We use an attack-
independent framework to structure our survey which results in an effective
categorisation of the field and an easy comparison of components. We hope this
survey will guide new researchers to this field and spark further research
into the individual attack components.

##  1 Introduction

Adversarial machine learning investigates how a model can be exploited by an
adversary through a carefully crafted adversarial attack Kurakin et al.
(2017). The majority of adversarial machine learning research has been for
image processing, but there has been increasing interest in exploring the
natural language processing (NLP) domain. Motivated by this recent growth, we
survey the recent work in this area and recommend several research directions.

An adversarial example is an example created by adding a perturbation to an
original example so that a model processes the original as expected but fails
on the adversarial in some prespecified manner. As our running example,
Wallace et al. (2020) wish to perturb a document that reads “Save me it’s over
100°F” in order to introduce a specific token 22°C into the output of an
English-to-German translation model. Their attack algorithm creates the
adversarial example by replacing “100°F” with “102°F”. While the model
translates the original example correctly as “Rette mich, es ist über 100°F”,
it now translates the adversarial example as “Rette mich, es ist über 22°C”,
inadvertently introducing the specific token.

Figure 1: Example of a token-modification adversarial attack on a machine
translation model.

Usually text adversarial examples must also meet a set of constraints. For
example, one constraint might specify that the semantic meaning of the
original example is maintained. In this case, “Save me it’s above 100°F” would
be valid, but not the proposed “Save me it’s over 102°F”. A second constraint
might require that the adversarial example appear non-suspicious to a human.
This rules out “Save fffffff it’s qqqqq over 102°F”, but our example “Save me
it’s over 102°F” would pass.

The most studied type of text adversarial attacks are token-modification
attacks: those modify individual document tokens of some granularity (e.g.
words) to achieve attack success. Our running example is a token-modification
attack. In these attacks, the adversary does not change the rough syntactic
structure of the document. Rather, they treat a document as an ordered set of
tokens, and repeatedly apply transformations to the individual tokens, such as
replacing a word with its synonym, deleting a character, or inserting a phrase
into the document. The attack then stops once the adversary reaches its goal,
or alternately, if some stopping criteria is reached.

Previous surveys—the two most comprehensive being Zhang et al. (2020) and
Alshemali and Kalita (2020)—categorise text attacks by token granularity,
adversary goal, model access, NLP task attacked, and so on. This approach
over-emphasises the differences between attacks. Many attack algorithms only
differ from each other in their use of one component (e.g. the attacks
proposed by Zhang et al. (2019)), which means algorithms can change their
categorisation with only trivial modifications. We solve this problem by
categorising individual attack components rather than the attacks themselves.
We are the first survey to take this approach.

##  2 Attack Framework

Recently Morris et al. (2020a) proposed a framework that defines four
components to a textual adversarial attack: (a) a goal function that specifies
if the attack was successful; (b) a set of transformations that define the
precise actions available to the adversary; (c) a search method that specifies
the order transformations are applied; and (d) a set of constraints to define
a valid adversarial example. Each text attack becomes a constrained search
problems defined by its choice of components.

Morris et al. (2020a) use their framework to write a software library to
quickly implement text adversarial attacks. Motivated by this work, we extend
this framework and systematically survey the range of attack components. Due
to space requirements, we restrict our survey to token-modification
attacks—since they are the predominant attack type—but our framework is
broadly applicable to other attack classes. Our taxonomy and covered attacks
are summarized in Table 1.

Notation. Let $Ex$ be the original example with label $y_{true}$. All examples
are documents, and each document is made up of a set of tokens
$\\{t_{1},\cdots,t_{n}\\}$, for which a token $t_{i}$ has the numeric
representation ${\bm{e}}_{t_{i}}$ (e.g. word embeddings). The document $Ex$ is
transformed to a numeric representation ${\bm{x}}$ and then acted on by a
model $f$. Likewise, the adversarial example is $Ex^{\prime}$ (also with label
$y_{true}$), its numeric representation is ${\bm{x}}^{\prime}$, and we denote
its output $f({\bm{x}}^{\prime})$ by $y^{\prime}$. We use $l$ to represent a
generic loss function that acts on a single example.

|  | Category
---
| Papers
---
Adv. Goals |  | Classification
---
| Ren et al. (2019); Jin et al. (2019); Li et al. (2020b, 2019)
---
Gao et al. (2018); Eger et al. (2019); Zang et al. (2020); Alzantot et al.
(2020)
Zhang et al. (2019); Liang et al. (2018); Hsieh et al. (2019); Xu and Du
(2020)
Wang et al. (2019); Samanta and Mehta (2017); Garg and Ramakrishnan (2020)
| Seq2Seq
---
| Hsieh et al. (2019); Michel et al. (2019); Camburu et al. (2020)
---
Jia and Liang (2017); Ebrahimi et al. (2018a); Belinkov and Bisk (2018)
Cheng et al. (2020); Wallace et al. (2020); Papernot et al. (2016); Zhao et
al. (2018)
Transformations |  | Any
---
replacement
| Ebrahimi et al. (2018b, a); Gao et al. (2018); Cheng et al. (2020)
---
Hsieh et al. (2019); Papernot et al. (2016); Wallace et al. (2020)
| Human errors
---
| Belinkov and Bisk (2018); Li et al. (2019); Hsieh et al. (2019)
---
| Visually similar
---
| Eger et al. (2019); Li et al. (2019)
---
| Synonym list
---
| Ren et al. (2019); Zang et al. (2020); Tan et al. (2020); Samanta and Mehta
(2017)
---
| Embedding space
---
| Li et al. (2019); Xu and Du (2020); Alzantot et al. (2020)
---
Jin et al. (2019); Wang et al. (2019); Li et al. (2020b)
Zhang et al. (2019); Hsieh et al. (2019); Michel et al. (2019)
| Language model
---
predictions
| Alzantot et al. (2020); Zhang et al. (2019); Li et al. (2020b)
---
Garg and Ramakrishnan (2020); Zhang et al. (2019)
Search Methods |  | All at once
---
| Tan et al. (2020); Ebrahimi et al. (2018a); Belinkov and Bisk (2018); Eger
et al. (2019)
---
| Document order
---
| Zhang et al. (2019)
---
| Random order
---
| Yoo et al. (2020); Hsieh et al. (2019)
---
| Importance
---
ranking order
| Li et al. (2019); Xu and Du (2020); Li et al. (2020b); Samanta and Mehta
(2017)
---
Jin et al. (2019); Gao et al. (2018); Garg and Ramakrishnan (2020); Ren et al.
(2019)
Hsieh et al. (2019); Ebrahimi et al. (2018b, a)
| Beam/greedy search
---
| Ebrahimi et al. (2018b); Wallace et al. (2020); Ebrahimi et al. (2018a);
Michel et al. (2019)
---
| Genetic algorithms
---
| Alzantot et al. (2020); Wang et al. (2019)
---
| Particle swarm opt.
---
| Zang et al. (2020)
---
Constraints |  | Readability
---
| Belinkov and Bisk (2018)
---
| Semantic
---
| Zang et al. (2020); Cheng et al. (2020); Ebrahimi et al. (2018b); Wang et
al. (2019)
---
Jin et al. (2019); Xu and Du (2020); Tan et al. (2020); Li et al. (2020b)
Garg and Ramakrishnan (2020); Zhang et al. (2019); Alzantot et al. (2020)
| Distance
---
| Xu and Du (2020); Ebrahimi et al. (2018b, a)
---
Gao et al. (2018); Li et al. (2019); Hsieh et al. (2019)
| Performance
---
| Xu and Du (2020); Ebrahimi et al. (2018b); Ren et al. (2019); Samanta and
Mehta (2017)
---

Table 1: Our proposed categorization of fundamental components of token-
modification attacks.

##  3 Adversary Goals

The adversary goal describes what the adversary intends to achieve. They are
specific to a task so we categorise by this below.

###  3.1 Classification

Attacks on text classifiers are usually either targeted or untargeted attacks.
Targeted attacks attempt to induce the model to classify $Ex^{\prime}$ as a
specific class $y_{target}$ different to $y_{true}$. They do this by
minimising $l(f({\bm{x}}),y_{target})$. Untargeted attacks only require that
$Ex$ is misclassified, and their algorithms typically aim to maximise
$l(f({\bm{x}}),y_{true})$. The two are equivalent for binary classification,
and for mulitclass classification, converting between the two is usually
straightforward 111For example, instead of perturbing ${\bm{x}}$ towards a
specific decision boundary, a untargeted attack perturbs ${\bm{x}}$ towards
the nearest decision boundary..

###  3.2 Sequence-to-Sequence

There are a wide range of sequence-to-sequence (seq2seq) tasks in NLP,
including machine translation, dialog systems, text summarisation, paraphrase
generation, speech recognition, speech synthesis and image generation from
captions. Unlike classification, there are no standardised adversary goals for
seq2seq attacks, so we borrow terminology and broadly classify these as either
untargeted or targeted.

#### Untargeted attacks.

These attacks simply attempt to degrade the output sequence quality by as much
as possible, as quantified by some metric. This is the most common goal for
seq2seq attacks Hsieh et al. (2019); Michel et al. (2019); Camburu et al.
(2020); Jia and Liang (2017); Ebrahimi et al. (2018a).

#### Targeted attacks.

These attacks attempt to induce specific changes in the output sequence. They
include: (a) changing a specific output token to any other token Papernot et
al. (2016); Ebrahimi et al. (2018a); Zhao et al. (2018), (b) changing every
output token from its original value Cheng et al. (2020), (c) changing a set
of specific output tokens to another set of specific tokens Wallace et al.
(2020); Ebrahimi et al. (2018a), (d) introducing a set of specific output
tokens anywhere in the output sequence Hsieh et al. (2019); Cheng et al.
(2020); Zhao et al. (2018), or (e) finding an incorrect input sequence that
generates a specific output sequence Wallace et al. (2020). These attacks all
involve carefully defining an attack-dependent loss function and subsequently
maximising or minimising it.

##  4 Transformations

A transformation is a function that transforms one ordered set of tokens to
another. We can place these into broad categories, which we call here
transformation types. There are four transformation types commonly used by
token-modification attacks: replacing one token with another, inserting a
token, deleting a token, or changing token ordering. Replacing a token is by
far the most popular transformation type, so we give it special focus below.

When replacing a token $t$ with a replacement $t_{r}$, the attack algorithm
will select $t_{r}$ from a set of replacements. We denote this set by
${\mathbb{S}}(t)$, parameterised to indicate its dependence on $t$. Specifying
a method to construct ${\mathbb{S}}(t)$ is a core component of many attack
algorithms but identifying the best method remains an open problem. The
standard approach is to first generate an initial set of candidates and then
to apply constraints to filter the set. For example, an algorithm could use a
dictionary to generate a set of synonyms and then filter out any using a
different part-of-speech (POS) tag to $t$.

Below we describe different methods that have been used to generate the
initial set of candidates, while deferring filtering constraints to Section 6.
Attacks are free to use a combination of these methods.

###  4.1 Any Replacement

The adversary replaces $t$ with any token they like of the same granularity.
This is both the simplest approach and the most effective at achieving the
adversary’s goal. Yet $Ex^{\prime}$ is very likely to be ungrammatical and to
lose its semantic meaning, thus this approach is usually not followed. It is
most often seen with character-level attacks that allow any alphanumeric
character (upper or lower case) or punctuation as a valid replacement Ebrahimi
et al. (2018b, a); Gao et al. (2018). These attacks make no attempt to
maintain semantics and focus instead on making the perturbation as small as
possible. It is also sometimes seen when adversarial attacks for images are
applied naively to the text embeddings (see Liang et al. (2018) for an
example).

###  4.2 Human Errors

The adversary replaces $t$ with a token that resembles an innocent human-made
mistake, such as a typo. This is intended to disguise $Ex^{\prime}$ as a
poorly-written document so that it is not considered suspicious by a human.
One method is to replace letters with adjacent characters on a QWERTY
keyboard, as to mimic a keyboard typo Belinkov and Bisk (2018); Li et al.
(2019). Another is to replace a word with a commonly-misspelled variant (e.g.
achieve with acheive). Common human misspellings are easy to find: there are
resources (e.g. Wikipedia) that list them, or alternatively, Belinkov and Bisk
(2018) extract their list from corpora tracking the edit history of documents.

###  4.3 Visually Similar

The adversary replaces $t$ with a token that looks similar (e.g. replacing I
with l) with the intention that a human overlooks the error. Again, usually
here $t$ is a character. This can also include many-to-one or one-to-many
character replacements (e.g. replacing cl with d). These attacks tend to be
either almost indistinguishable or very obvious, depending on the font used
and the choice of replacement.

The simplest method is to create a mapping between each token in the
vocabulary and a manually-chosen replacement for each Eger et al. (2019); Li
et al. (2019). For example, Eger et al. (2019) replace characters by a version
of themselves with a diacritic (e.g. c to č). This work also proposes a
different method that uses Unicode descriptions of characters to find
replacements of the same character and case222For example, a has Unicode
description LATIN SMALL LETTER A, and a potential replacement à has Unicode
description LATIN SMALL LETTER A WITH GRAVE. An alternative is to use visual
embeddings that treat each character as a grid of pixels. Each character is
then represented numerically, such as with a raw pixel vector Eger et al.
(2019), or as a dense representation, like from the intermediate layers of a
CNN Li et al. (2020a). At this point ${\mathbb{S}}(t)$ can be constructed by
taking the closest points to the character in that latent space.

###  4.4 Synonym List

The adversary replaces $t$ with a replacement $t_{r}$ taken from a predefined
list of tokens that have the same meaning. It is most frequently used with
word-level attacks. Tan et al. (2020) use a set of inflectional perturbations
of words which all have the same base meaning. More common is using is a set
of synonyms from a thesaurus, lexical database, knowledge database or similar.
Examples include Ren et al. (2019) who use WordNet Fellbaum (1998) as their
lexical database, and Zang et al. (2020) which use the knowledge database
HowNet Dong et al. (2010). This method is often unable to distinguish the
sense in which a word is used, so replaced words often read strangely or do
not fit well in the context of $Ex^{\prime}$. Sometimes an attack will add
constraints to try to mitigate this problem (e.g. Zang et al. (2020)) but in
practice the constraints do not solve the problem very well.

###  4.5 Embedding Space

The adversary replaces token $t$ with a token $t_{r}$ corresponding to a
nearby point in a token embedding space. The adversary assumes that since
tokens co-occurring in similar contexts cluster together in embedding space,
these tokens also have similar meaning, and hence the nearby points around $t$
form a natural set of synonyms. This is most commonly used with word
embeddings. Li et al. (2019); Xu and Du (2020) both do this with GloVe
Pennington et al. (2014) embeddings.

Yet in practice, close points in a word embedding space (e.g. GloVe Pennington
et al. (2014)) include similar words (like word synonyms), but also related
words, such as antonyms (e.g. east and west, or fast and slow) or words with
similar concepts but different meanings (e.g. English, American, and
Australian). Hence this approach usually leads to transformations that change
the semantic content of $Ex$. This is mitigated by instead using counter-
fitted embeddings Mrkšić et al. (2016). Counter-fitting is a fast post-
processing step for word embeddings that pushes together similar words while
also pushing apart related words. The result is a word embedding where only
synonyms are close to each other. Hence better results are obtained if the
adversary uses a counter-fitted word embedding space. Alzantot et al. (2020);
Jin et al. (2019); Wang et al. (2019); Li et al. (2020b) all use counter-
fitted GloVe embeddings.

What makes a nearby token? Common methods include taking the $k$-nearest
neighbours (KNN) to the original token (e.g. Li et al. (2019); Xu and Du
(2020), or to consider tokens within a distance333Usually this is Euclidean
distance or cosine distance/similarity. $\delta$ to $t$ (e.g. Zhang et al.
(2019)), or to do both and take the intersection of the sets (e.g. Alzantot et
al. (2020); Jin et al. (2019)). Hsieh et al. (2019) does this too, but with
document embeddings instead of word embeddings. Wang et al. (2019) create
synonym clusters in word-embedding space (using their own clustering
algorithm) intending that any word in the cluster can replace a “central”
word, although their examples indicate problems with this approach.

A problem with this class of approaches is that token embeddings can only use
a single point per token, but many tokens are used in a range of contexts.
This means replacements are often for the wrong sense of a token.

###  4.6 Language Model Predictions

The adversary uses a language model $LM$ to find replacements for $t$. The
usual way is through masked language modelling: replace $t$ in $Ex$ with a
mask token, get predictions from $LM$ for the mask token, and use the top $k$
predictions as ${\mathbb{S}}(t)$. Alzantot et al. (2020); Zhang et al. (2019);
Li et al. (2020b); Garg and Ramakrishnan (2020) variations on this approach.
Alzantot et al. (2020); Zhang et al. (2019) use both language model
predictions and word embeddings to rank replacements.

This approach has two main advantages. Firstly, generated replacements are
higher quality since the language model accounts for token context. Secondly,
rare tokens are seldom suggested, which is advantageous since they tend to
appear suspicious Ippolito et al. (2020). There are disadvantages too: the
attack requires a complex and memory-intensive language model, the language
model is slower to generate replacements, and the attack cannot cache synonyms
between documents (since they depend on document context). Nonetheless, the
quality increase in replacements is usually high enough that the advantages
outweigh the disadvantages.

##  5 Search Methods

The search method attempts to find a successful sequence of transformations
that solves the constrained search problem. Search methods rank
transformations with a scoring function, a function estimating the effect of
the transformation on the loss function.

The simplest scoring function directly calculates the loss increase of the
transformation. This measures the impact exactly, however the approach is
slow—particularly for large models. Hence a popular alternative is to use a
gradient-based scoring function that uses some variant of $\partial
l/\partial{\bm{x}}$: a first-order approximation of the effect on $l$ when
changing elements of ${\bm{x}}$.444Some attacks (e.g. Papernot et al. (2016))
instead use the Jacobian matrix $\partial f(x)/\partial{\bm{x}}$ and directly
change model predictions, but the loss-function approach is predominant..
Since the gradient approximates the direction of steepest loss increase,
transformations can be ranked by their dot product $\partial
l/\partial{\bm{x}}\cdot(g({\bm{x}})-{\bm{x}})$, where the vector
$g({\bm{x}})-{\bm{x}}$ measures the relative change in input-space of a
transformation $g$. Both the direction and the magnitude of this vector are
important, so dot product is a better metric to use than cosine similarity.

The gradient approximation only needs a single backwards pass through the
model, but requires the adversary can access model parameters (white-box
attack). It also has two implicit assumptions: that $l$ is close to linear
around ${\bm{x}}$, and that all other features are fixed when evaluating the
importance of a single feature Singla et al. (2019). In most cases probably
neither assumption holds, but regardless the search method usually can find a
successful adversarial example.

Wallace et al. (2020) use a hybrid approach which has advantages of both
approaches: the speed of the gradient approximation, and the accuracy of
direct-checking. They first use the gradient approximation to create a
shortlist of the top $k$ transformations across all tokens in $Ex$. They then
find the best transformation by directly checking the loss change of each.

Below we present common search algorithms. Search methods either stop when the
adversary goal is reached or when the constraints prevent them from searching
further.

###  5.1 All at Once

Rather than apply transformations one-by-one, this method creates
$Ex^{\prime}$ by simultaneously applying a transformation to all applicable
tokens in $Ex$. This is done in Tan et al. (2020); Ebrahimi et al. (2018a),
who both replace each applicable word $w\in Ex$. Belinkov and Bisk (2018)
allow for multiple runs of this procedure, with each run selecting a different
transformation. Eger et al. (2019) implement a variation where they replace
each token with a probability $p$.

###  5.2 Document Order

This approach transforms tokens in document order. Attacks choose the best
transformation for the first applicable token $t_{1}$ in $Ex$, then for
$t_{2}$, and so on until the attack is successful. The algorithm can cycle
through the document until it is successful. Zhang et al. (2019) take this
approach, but rather than select the best transformation, they have predefined
probabilities for several different transformation types.555For each word $w$,
they either replace it with a word from ${\mathbb{S}}(w)$, delete it, or
insert a word in front of it, with predefined probabilities 0.5, 0.25, and
0.25 respectively.

###  5.3 Random Order

This search method randomly selects and replaces a token from the document.
This method is used as a baseline by Yoo et al. (2020); Hsieh et al. (2019),
where tokens are either replaced by the best replacement or by a random
replacement.

###  5.4 Importance Ranking Order

This method ranks the influence of each token on the loss function and then
applies transformations in that order. This keeps the perturbation small since
more influential tokens are transformed first. Yoo et al. (2020) call this an
importance ranking.

There are several common methods to calculate an importance ranking. The first
uses gradient information. Intuitively, elements of ${\bm{x}}$ with a high
gradient magnitude $\partial l/\partial{\bm{x}}$ affect $l$ more, with the
direction depending on the sign of the gradient. The quantity $\partial
l/\partial{\bm{x}}$ cannot be used directly to rank tokens since the length of
${\bm{x}}$ is a multiple of the number of tokens. The common workaround is to
calculate importance of a token $t$ by the magnitude of its loss gradient
vector $||\partial l/\partial{\bm{e}}_{t_{i}}||_{2}$. The drawback is that
using only gradient magnitude and not direction means the transformation
impact cannot be predicted precisely. Nevertheless, this remains a sound
heuristic and many attacks use this method Li et al. (2019); Xu and Du (2020);
Samanta and Mehta (2017).

Another method is to systematically perturb each token in the document and
measure the change in loss. There are a number of different perturbation
methods. One is to remove the token entirely Jin et al. (2019); Garg and
Ramakrishnan (2020); Li et al. (2019). Another is to replace it with a
“neutral” token that theoretically does not affect model predictions. The
neutral token has been the out-of-vocab/unknown token Ren et al. (2019), the
zero vector Ebrahimi et al. (2018b), and the mask token Li et al. (2020b).
Another two methods were proposed by Gao et al. (2018). They first split the
document into a prefix, which contains the document up to (but not including)
the token, and a suffix, which contains the document from after the token
until the end of the document. The two metrics are then the temporal head
score, given by $J(f(\text{prefix})-J(f(\text{prefix + token}))$ and the
temporal tail score, given by $J(f(\text{token +
suffix})-J(f(\text{suffix}))$.

A third method calculates the importance ranking with attention scores. This
method only works if the victim model has an attention mechanism. Hsieh et al.
(2019) use the attention scores of the first encoder-decoder stack, or if the
model has multiple heads, the average of the attention layer over each head.
The best importance ranking ranked scores from highest to lowest.

How do these methods compare? Yoo et al. (2020) reported similar performance
testing these methods666They did not test the attention method., though the
gradient approach was more computationally efficient777Systematic perturbation
methods need many forward passes through the model, but gradient methods only
need single forwards and backward pass..

###  5.5 Beam Search and Greedy Search

Beam search is a search algorithm that keeps track of multiple promising paths
at once. The beam width parameter controls the number of paths. Larger beam
widths find better solutions, but take more time to evaluate. Setting a beam
width of one gives greedy search, where the best transformation is selected at
each step.

Beam search has been used as the search method in Ebrahimi et al. (2018b, a),
and greedy search in Michel et al. (2019); Wallace et al. (2020). Beam search
has been reported to have a high success rate while also being relatively
query-efficient for shorter documents Yoo et al. (2020).

###  5.6 Genetic Algorithms

Genetic algorithms solve optimisation problems by evolving a population of
candidate solutions over time, using a fitness function to perform “natural
selection” until one is successful. Variants have been applied to text attacks
by Alzantot et al. (2020); Wang et al. (2019). The algorithm starts by
applying $S$ different transformations to $Ex$ to create a candidate
population of size $S$. Then at each step a new candidate set is created by:
(a) applying transformations to the original $S$ candidates; (b) calculating
the loss for each, which is called the “fitness” score; (c) sampling $S$ pairs
of these documents with probability proportional to their fitness score; (d)
creating $S$ new candidates through applying a crossover function to the two
parents. The crossover function in Alzantot et al. (2020) simply selects each
word randomly from one of the two parent documents, while Wang et al. (2019)
instead split each parent document at a point and take one fragment from each.

###  5.7 Particle Swarm Optimisation

Particle swarm optimisation is a search method with both continuous and
discrete variants. It starts by initialising a number of “particles”. Each
particle has a position and a velocity, and keeps track of both its best
position and the best position achieved across all particles. Particles move
with a random component and a probability to move towards either the global
best or local best, depending on the algorithm variant.

Zang et al. (2020) use particle swarm optimisation as their search method. Yoo
et al. (2020) reported that this search method reliably leads to a good
solution, but is query inefficient.

##  6 Constraints

Attacks can consider a set of constraints. Weaker constraints increase attack
success, but also find an inferior adversarial example. Constraints can have
many benefits on $Ex^{\prime}$: reducing semantic drift from $Ex$, increasing
readability or naturalness, keeping $Ex$ and $Ex^{\prime}$ “close” together,
ensuring $Ex^{\prime}$ has correct grammar, and so on.

We group constraints into rough categories based on their intended effect.

###  6.1 Readability

These constraints intend to increase human readability of $Ex^{\prime}$. Some
are only applicable for character-level attacks. One of these is to apply a
spellchecker to $Ex^{\prime}$ and see if the result still achieves the
adversary’s goal. Belinkov and Bisk (2018) use other constraints. They only
perform character transformations on words of a certain length, and avoid
changing the first and last letters of a word.

More generally, adversarial attacks might require $Ex^{\prime}$ have correct
grammar, as checked by grammar software. Sometimes this constraint is not
appropriate: for example, the attack by Tan et al. (2020) deliberately mimics
non-native English speakers, so the constraint would not be used then. Another
constraint requires adversarial examples maintain some number of
linguistically acceptable sentences Morris et al. (2020a). This can be
measured with a model trained on a linguistic acceptability dataset.

###  6.2 Semantic

These constraints aim to maintain the semantic content of the original
document. A widely used constraint only allows one transformation on each
token (e.g. Cheng et al. (2020); Ebrahimi et al. (2018b)), or more generally,
up to $k$ transformations per token Wang et al. (2019)). Avoiding multiple
transformations on a token eliminates a likely source of semantic drift.
Another type of constraint specifies the tokens a transformation can act on.
Here the adversary wishes to avoid changing non-influential tokens that are
also vital for document structure. Sometimes these tokens are approximated by
a list of stop words, and hence one version of this constraint explicitly
forbids replacing stop words Jin et al. (2019); Xu and Du (2020). This is
sometimes reformulated as only allowing replacements for certain word types.
For example, Tan et al. (2020) only allows nouns, verbs and adjectives to be
replaced, while Zang et al. (2020) also allows replacing adverbs. A variant
used by Li et al. (2020b) instead disallows stop words in the set of
replacements.

Other constraints filter the set of replacements to improve their quality,
where quality is judged as similarity to the original token, or the extent to
which the document reads naturally. One method requires replacements to have
the same part of speech tag as the original Jin et al. (2019); Zang et al.
(2020); Xu and Du (2020). This helps to eliminate potential synonyms that use
an incorrect word sense. Another measures entropy of a language model for each
replacement when it is put into context, and only considers replacements above
a threshold Alzantot et al. (2020).

Sometimes the distance between original and adversarial document
encodings888Document encodings are numeric representations of documents,
similar to how word embeddings are numeric representations of words. is used
to estimate semantic similarity. At each attack step, the adversary finds
document encodings of $Ex$ and $Ex^{\prime}$, and then computes the cosine
similarity of the two. Then if the cosine similarity is below a threshold, the
adversary avoids transforming the token Jin et al. (2019), or alternatively,
they can rank different transformations by cosine similarityGarg and
Ramakrishnan (2020). A commonly used document encoder for this purpose is the
Universal Sentence Encoder Cer et al. (2018).

An adversary might also apply task-specific constraints to $Ex^{\prime}$ to
maintain semantic similarity. For example, for sentiment analysis Zhang et al.
(2019) do not replace any words that they determine express sentiment. They
also avoids replacing negation words.

###  6.3 Distance

These constraints constrain the distance between $Ex$ and $Ex^{\prime}$, where
the distance metric is usually Hamming or Levenshtein edit distance999Other
metrics include n-gram based metrics (e.g. BLEU Papineni et al. (2002)) and
Word Mover Distance (as proposed by Zhang et al. (2020)) but neither one has
been widely used.. Hamming distance measures the minimum number of token
replacements required to go from $Ex$ to $Ex^{\prime}$. Levenshtein distance
is similar but also allows inserting and deleting tokens. Several attacks use
an edit distance constraint Ebrahimi et al. (2018a); Gao et al. (2018); Li et
al. (2019).

###  6.4 Performance

These constraints attempt to improve attack performance They include: adding
the padding and unknown tokens to the synonym set Xu and Du (2020),
lemmatisation/delemmatisation of candidate synonyms to create more potential
replacements Zang et al. (2020), or only allowing character transformations if
the transformed word does not exist in the vocabulary101010This is not
intended to maintain readability, but rather to ensure word-level models
represent the replacement tokens with out-of-vocab tokens, rather than with
another word. Ebrahimi et al. (2018b).

Some attacks prioritise transformations that introduce tokens over-represented
in other classes to $y_{true}$. The reasoning is that these transformations
are likely to be effective while also appearing natural to humans. One example
is Ren et al. (2019) who replace named entities with the most frequent same-
type named entities from a different class. Samanta and Mehta (2017) first
assign each example a category (e.g. movie genre for a dataset of movie
reviews), and then prioritise word replacements that are both used frequently
in that category and also in different-class examples.

##  7 Conclusions and Future Directions

Token-modification attacks are simple permutations of our presented attack
components. We encourage further research into these components, such as
developing new scoring functions, finding better sets of replacements, or
identifying better search algorithms. We also make two further suggestions.

#### Improve semantic similarity constraints.

As it stands semantic similarity is best evaluated by a human Morris et al.
(2020b) since current automated semantic similarity metrics are not yet close
to human judgement. Improving these automated metrics will in turn increase
the text quality of attacks substantially.

#### Extend to multimodal systems.

The majority of adversarial attacks attack systems of a single modality, such
as images or text. Yet real-world data is often multimodal, consisting of
images, videos, and text (e.g., web pages, emails, social media). Extending
text attacks to use multimodal information is another promising future
direction.

## References

  * Alshemali and Kalita (2020) B. Alshemali and J. Kalita.  Improving the Reliability of Deep Neural Networks in NLP: A Review.  Knowledge-Based Systems, 191:105210, 2020.
  * Alzantot et al. (2020) M. Alzantot, Y. Sharma, A. Elgohary, B. J. Ho, M. B. Srivastava, and K. W. Chang.  Generating Natural Language Adversarial Examples.  Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing, EMNLP 2018, pages 2890–2896, 2020.
  * Belinkov and Bisk (2018) Y. Belinkov and Y. Bisk.  Synthetic and Natural Noise Both Break Neural Machine Translation.  In 6th International Conference on Learning Representations, {ICLR} 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.
  * Camburu et al. (2020) O.-M. Camburu, B. Shillingford, P. Minervini, T. Lukasiewicz, and P. Blunsom.  Make Up Your Mind! Adversarial Generation of Inconsistent Natural Language Explanations.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 4157–4165, Online, 7 2020. Association for Computational Linguistics.
  * Cer et al. (2018) D. Cer, Y. Yang, S. y. Kong, N. Hua, N. Limtiaco, R. St. John, N. Constant, M. Guajardo-Céspedes, S. Yuan, C. Tar, Y. H. Sung, B. Strope, and R. Kurzweil.  Universal sentence encoder for English.  EMNLP 2018 - Conference on Empirical Methods in Natural Language Processing: System Demonstrations, Proceedings, abs/1803.1:169–174, 2018.
  * Cheng et al. (2020) M. Cheng, J. Yi, P.-Y. Chen, H. Zhang, and C.-J. Hsieh.  Seq2Sick: Evaluating the Robustness of Sequence-to-Sequence Models with Adversarial Examples.  Proceedings of the AAAI Conference on Artificial Intelligence, 34(04):3601–3608, 2020.
  * Dong et al. (2010) Z. Dong, Q. Dong, and C. Hao.  HowNet and Its Computation of Meaning.  In Coling 2010: Demonstrations, pages 53–56, Beijing, China, 8 2010\. Coling 2010 Organizing Committee.
  * Ebrahimi et al. (2018a) J. Ebrahimi, D. Lowd, and D. Dou.  On Adversarial Examples for Character-Level Neural Machine Translation.  In Proceedings of the 27th International Conference on Computational Linguistics, pages 653–663, Santa Fe, New Mexico, USA, 8 2018\. Association for Computational Linguistics.
  * Ebrahimi et al. (2018b) J. Ebrahimi, A. Rao, D. Lowd, and D. Dou.  HotFlip: White-Box Adversarial Examples for Text Classification.  In I. Gurevych and Y. Miyao, editors, ACL 2018, Melbourne, Australia, July 15-20, 2018, Volume 2: Short Papers, pages 31–36. Association for Computational Linguistics, 2018.
  * Eger et al. (2019) S. Eger, G. G. İşgüder, A. Rücklé, J.-U. Lee, C. Schulz, M. Mesgar, K. Swarnkar, E. Simpson, and I. Gurevych.  Text Processing Like Humans Do: Visually Attacking and Shielding NLP Systems.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 1634–1647, Minneapolis, Minnesota, 6 2019. Association for Computational Linguistics.
  * Fellbaum (1998) C. Fellbaum.  WordNet: An Electronic Lexical Database.  Bradford Books, 1998.
  * Gao et al. (2018) J. Gao, J. Lanchantin, M. L. Soffa, and Y. Qi.  Black-Box Generation of Adversarial Text Sequences to Evade Deep Learning Classifiers.  Proceedings - 2018 IEEE Symposium on Security and Privacy Workshops, SPW 2018, pages 50–56, 2018.
  * Garg and Ramakrishnan (2020) S. Garg and G. Ramakrishnan.  BAE: BERT-based Adversarial Examples for Text Classification.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 6174–6181, Online, 11 2020. Association for Computational Linguistics.
  * Hsieh et al. (2019) Y.-L. Hsieh, M. Cheng, D.-C. Juan, W. Wei, W.-L. Hsu, and C.-J. Hsieh.  On the Robustness of Self-Attentive Models.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 1520–1529, Florence, Italy, 7 2019. Association for Computational Linguistics.
  * Ippolito et al. (2020) D. Ippolito, D. Duckworth, C. Callison-Burch, and D. Eck.  Automatic Detection of Generated Text is Easiest when Humans are Fooled.  In ACL, pages 1808–1822, 2020.
  * Jia and Liang (2017) R. Jia and P. Liang.  Adversarial Examples for Evaluating Reading Comprehension Systems.  In Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, pages 2021–2031, Copenhagen, Denmark, 9 2017. Association for Computational Linguistics.
  * Jin et al. (2019) D. Jin, Z. Jin, J. T. Zhou, and P. Szolovits.  Is BERT Really Robust? A Strong Baseline for Natural Language Attack on Text Classification and Entailment.  The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020, pages 8018–8025, 2019.
  * Kurakin et al. (2017) A. Kurakin, I. J. Goodfellow, and S. Bengio.  Adversarial Machine Learning at Scale.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net, 2017.
  * Li et al. (2019) J. Li, S. Ji, T. Du, B. Li, and T. Wang.  TextBugger: Generating Adversarial Text Against Real-world Applications.  In 26th Annual Network and Distributed System Security Symposium, NDSS 2019, San Diego, California, USA, February 24-27, 2019, number February. The Internet Society, 2019.
  * Li et al. (2020a) J. Li, T. Du, S. Ji, R. Zhang, Q. Lu, M. Yang, and T. Wang.  TextShield: Robust Text Classification Based on Multimodal Embedding and Neural Machine Translation.  In 29th USENIX Security Symposium (USENIX Security 20), pages 1381–1398. USENIX Association, 8 2020.
  * Li et al. (2020b) L. Li, R. Ma, Q. Guo, X. Xue, and X. Qiu.  BERT-ATTACK: Adversarial Attack Against BERT Using BERT.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 6193–6202, Online, 11 2020. Association for Computational Linguistics.
  * Liang et al. (2018) B. Liang, H. Li, M. Su, P. Bian, X. Li, and W. Shi.  Deep Text Classification Can be Fooled.  In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, {IJCAI-18}, pages 4208–4215. International Joint Conferences on Artificial Intelligence Organization, 2018\.
  * Michel et al. (2019) P. Michel, X. Li, G. Neubig, and J. Pino.  On Evaluation of Adversarial Perturbations for Sequence-to-Sequence Models.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 3103–3114, Minneapolis, Minnesota, 6 2019. Association for Computational Linguistics.
  * Morris et al. (2020a) J. Morris, E. Lifland, J. Y. Yoo, J. Grigsby, D. Jin, and Y. Qi.  TextAttack: A Framework for Adversarial Attacks, Data Augmentation, and Adversarial Training in NLP.  pages 119–126, 2020.
  * Morris et al. (2020b) J. X. Morris, E. Lifland, J. Lanchantin, Y. Ji, and Y. Qi.  Reevaluating Adversarial Examples in Natural Language.  2020\.
  * Mrkšić et al. (2016) N. Mrkšić, D. Ó Séaghdha, B. Thomson, M. Gašić, L. Rojas-Barahona, P.-H. Su, D. Vandyke, T.-H. Wen, and S. Young.  Counter-fitting Word Vectors to Linguistic Constraints.  In Proceedings of HLT-NAACL, 2016.
  * Papernot et al. (2016) N. Papernot, P. D. McDaniel, A. Swami, and R. E. Harang.  Crafting Adversarial Input Sequences for Recurrent Neural Networks.  In J. Brand, M. C. Valenti, A. Akinpelu, B. T. Doshi, and B. L. Gorsic, editors, 2016 IEEE Military Communications Conference, MILCOM 2016, Baltimore, MD, USA, November 1-3, 2016, pages 49–54. IEEE, 2016.
  * Papineni et al. (2002) K. Papineni, S. Roukos, T. Ward, and W.-J. Zhu.  BLEU: A Method for Automatic Evaluation of Machine Translation.  In Proceedings of the 40th Annual Meeting on Association for Computational Linguistics, ACL ’02, page 311–318, USA, 2002. Association for Computational Linguistics.
  * Pennington et al. (2014) J. Pennington, R. Socher, and C. D. Manning.  GloVe: Global Vectors for Word Representation.  In In EMNLP, 2014.
  * Ren et al. (2019) S. Ren, Y. Deng, K. He, and W. Che.  Generating Natural Language Adversarial Examples through Probability Weighted Word Saliency.  Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 1085–1097, 7 2019.
  * Samanta and Mehta (2017) S. Samanta and S. Mehta.  Towards Crafting Text Adversarial Samples.  arXiv, abs/1707.0, 2017.
  * Singla et al. (2019) S. Singla, E. Wallace, S. Feng, and S. Feizi.  Understanding Impacts of High-Order Loss Approximations and Features in Deep Learning Interpretation.  In K. Chaudhuri and R. Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5848–5856. PMLR, 2019.
  * Tan et al. (2020) S. Tan, S. Joty, M.-Y. Kan, and R. Socher.  It’s Morphin’ Time! Combating Linguistic Discrimination with Inflectional Perturbations.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 2920–2935, Online, 7 2020. Association for Computational Linguistics.
  * Wallace et al. (2020) E. Wallace, M. Stern, and D. Song.  Imitation Attacks and Defenses for Black-box Machine Translation Systems.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 5531–5546, Online, 11 2020. Association for Computational Linguistics.
  * Wang et al. (2019) X. Wang, H. Jin, and K. He.  Natural Language Adversarial Attacks and Defenses in Word Level.  2019\.
  * Xu and Du (2020) J. Xu and Q. Du.  TextTricker: Loss-Based and Gradient-Based Adversarial Attacks on Text Classification Models.  Eng. Appl. Artif. Intell., 92:103641, 2020.
  * Yoo et al. (2020) J. Y. Yoo, J. X. Morris, E. Lifland, and Y. Qi.  Searching for a Search Method: Benchmarking Search Algorithms for Generating NLP Adversarial Examples, 2020.
  * Zang et al. (2020) Y. Zang, F. Qi, C. Yang, Z. Liu, M. Zhang, Q. Liu, and M. Sun.  Word-level Textual Adversarial Attacking as Combinatorial Optimization.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 6066–6080, Online, 7 2020. Association for Computational Linguistics.
  * Zhang et al. (2019) H. Zhang, H. Zhou, N. Miao, and L. Li.  Generating Fluent Adversarial Examples for Natural Languages.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 5564–5569, Florence, Italy, 7 2019. Association for Computational Linguistics.
  * Zhang et al. (2020) W. E. Zhang, Q. Z. Sheng, A. Alhazmi, and C. Li.  Adversarial Attacks on Deep-learning Models in Natural Language Processing.  ACM Transactions on Intelligent Systems and Technology, 11(3):1–41, 2020.
  * Zhao et al. (2018) Z. Zhao, D. Dua, and S. Singh.  Generating Natural Adversarial Examples.  In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.


#  Data-Driven Reduced-Order Aeroelastic Modeling of Highly Flexible Aircraft
by Parametric Dynamic Mode Decomposition

Tianyi He, Weihua Su He is with Department of Mechanical and Aerospace
Engineering, Utah State University, Logan, Utah 84322-4130, USA. Email:
<EMAIL_ADDRESS>He is the corresponding author.Su is with Department of
Aerospace Engineering and Mechanics, The University of Alabama, Tuscaloosa,
Alabama 35487-0280, USA. Email<EMAIL_ADDRESS>
###### Abstract

This paper presents a method of data-driven parametric Dynamic Mode
Decomposition (p-DMD) to derive a linear parameter-varying reduced-order model
(LPV-ROM) for the nonlinear aeroelasticity of highly flexible aircraft. It
directly uses the data snapshots obtained at varying flight conditions, and
encodes the physical understanding of the nonlinear model’s polynomial
dependency on flight conditions to produce a polynomial-dependent LPV-ROM.
Therefore, this method can handle not only the equilibrium flight conditions
but also the cases of continuously-varying flight conditions. In the numerical
studies, a highly flexible cantilever wing and a slender vehicle built based
on it are first studied with fixed angles of attack as the scheduling
parameter. The comparisons between traditional linearization-based parametric
modeling and the data-driven p-DMD modeling are performed to verify the
modeling accuracy. The results demonstrate that the current p-DMD modeling
method can capture the aeroelastic and flight dynamic responses of highly
flexible aircraft in both time and frequency domains. In addition, the
proposed p-DMD method is applied to the highly flexible aircraft in a
perturbed longitudinal flight with varying angles of attack as the scheduling
parameter. The nonlinear aeroelastic and flight dynamic data are compared with
the simulation results of the data-driven p-DMD model. The comparison results
demonstrate that it can accurately capture the non-equilibrium (or transient)
aeroelastic and flight dynamic behaviors of such slender vehicles.

###### Index Terms:

Data-Driven Modeling; Linear Parameter-Varying System; Dynamic Mode
Decomposition; Highly Flexible Aircraft

##  I Introduction

High-aspect-ratio slender wings have been demanded to design next-generation
high-altitude, long-endurance aircraft due to the benefits of improved vehicle
performance and efficiency. However, the geometrically nonlinear effect
pertinent to the highly flexible wings must be adequately addressed in
modeling the open-loop dynamic behaviors, analyzing their stability
characteristics, and the control system development for such slender vehicles.
Accurately predicting the coupled nonlinear aeroelasticity and flight dynamics
of highly flexible aircraft is essential for their design, optimization, and
control. However, the complex fluid-structure interactions (FSI) induced by
the geometric nonlinearity, as well as the coupled nonlinear aeroelastic and
flight dynamic behavior, pose significant challenges to the nonlinear
aeroelastic modeling and active control of such vehicles.

One of the significant challenges is the nonlinear elastic deformation of
highly flexible wings, which makes the FSI more complicated, even in the
scenario of a simple wing geometry. Linearization is the most commonly used
method to predict aircraft stability and other open-loop performance, and
obtain control-oriented models [1]. However, linearization is conducted
assuming small vibrations about an equilibrium flight condition. Therefore, it
cannot well address the vehicle’s dynamic behavior with arbitrary large wing
deformations. Linear parameter varying (LPV) modeling and control have been
demonstrated to be a powerful tool that addresses the nonlinear aeroelasticity
of flexible wing [2, 3]. An LPV model can be established by interpolating
multiple linearized models. After that, an LPV controller can be designed to
be parameter-dependent to smoothly schedule controller gains at varying flight
conditions or vibration magnitude [4].

Another challenge is the development of reduced-order models (ROMs) to enable
fast computations of optimization and control, usually at the expense of
fidelity. The high-fidelity structural model can be analyzed and reduced to a
lower-order model, while the accuracy is almost maintained by keeping the most
significant modes. An eigenvalue problem is usually solved, and a set of
orthogonal modes are selected as states of the reduced-order model. Model
reductions are repeated at gridded flight conditions, and interpolations are
conducted to capture parametric dependency to get a parametric model over a
range of flight conditions. It’s quite often to encounter the cases at varying
flight conditions: 1) System stability varies from stable to unstable, or
varies from unstable to stable [3]; 2) A complex conjugate pair merge to two
real-valued poles or two real-valued poles split to complex conjugate pair
[5]. Therefore, the eigenvalues and eigenvectors need to be tracked and
aligned at gridded flight conditions [6] to ensure the consistency of states
among the fixed models.

Recent years have seen a growing interest in data-driven modeling techniques
for the aeroelastic analysis of flexible wings. There are many data-driven
methods of aeroelastic modeling using neural networks [7], Kriging [8]. These
methods leverage a large amount of data sets obtained from experiments or
numerical simulations, train a model to learn the input-output relationships
without a physical model, and achieve accurate predictions of aeroelastic
behaviors. These approaches are particularly well-suited for complex,
nonlinear systems where analytical models may be insufficiently precise or
difficult to establish. However, the computational complexity is usually not
time-efficient. Besides, these models lack an in-depth understanding and
interpretation of the physical system. Therefore, these models are unreliable
for flight control design, which is safety critical.

State-space model by the data-driven method is of particular interest from the
control point of view. Identification methods or mode decomposition can be
used to find the optimal reduced-order model that best approximates the
system. Eigensystem realization algorithm (ERA) [9] was developed for modal
parameter identification and model reduction from testing data. Autoregressive
with exogenous input (ARX) model was derived to identify the input-output
relationships in discrete-time domain [10, 11]. Note that neither ARX nor ERA
encodes the information of the physical states. If the physical state
information is available, a physically-interpretable state-space model can be
obtained via several methods. Proper orthogonal decomposition (POD) aims at
finding a proper low-dimensional subspace of spatially orthogonal modes that
best approximate the data in the $l_{2}$-norm manner [12, 13]. Sparse
identification of nonlinear dynamics (SINDy) can produce a physically
interpretable nonlinear model based on a library of nonlinear basis functions
[14, 15].

Dynamic Mode Decomposition (DMD) is an equation-free method that can reveal
the evolution of the spatial modes and temporal dynamics from the data
snapshots. The well-established DMD algorithm uses data at perturbed
equilibrium points, analyzes the mode, and renders a Linear Time-Invariant
(LTI) model [16, 5]. In order to obtain ROMs, the projection-based model
reduction method was often employed to keep the most significant modes [17]. A
dynamic mode decomposition with control (DMDc) was applied on a flexible wing
at different operating points and interpolated into a parametric model over a
range of operating regimes [5]. A balanced mode decomposition (BMD) was
derived using an oblique projection, where the projection basis space was
shared to achieve state consistency, and the testing projection was chosen as
parameter-dependent to obtain parametric models [16]. However, these works
were based on the assumption of small wing deformations, and the persistently
excited inputs were used to perturb the equilibrium points to obtain data
snapshots. The transient system dynamics caused by the varying scheduling
parameter were missing in the parametric model from interpolation.

This paper presents an innovative method that directly uses data snapshots at
varying flight conditions to derive a parametric LPV-ROM. The modeling
accuracy is verified by comparing it with the traditional linearization-
interpolation approach and nonlinear full-order models on the test bed of a
slender wing and a flexible aircraft in our previous work [3]. The results
reveal that, at fixed flight conditions of equilibrium, where a steady state
is achieved, these two methods render consistent responses in the frequency
domain and time domain. Furthermore, this paper demonstrates that the proposed
method is better than traditional DMD because it can also address the scenario
of varying flight conditions that are not equilibrium, where the aircraft
exhibit transient dynamic behaviors.

The novelties of the proposed method and the contributions of this paper are
three-fold:

  1. 1.

The proposed p-DMD is applicable not only at the equilibrium points, where a
steady state is achieved, to derive an LTI-ROM, but also at the non-
equilibrium scenarios to derive a polynomial-dependent LPV-ROM. The assumption
of small deformation is relaxed because no small perturbation about the
equilibrium is needed. Therefore, the proposed method can address the coupled
aeroelasticity and flight dynamics for the highly flexible wing with large
deformations and unsteady dynamics. Besides, p-DMD can capture the transient
system dynamics at varying flight conditions.

  2. 2.

It’s well known that traditional parametric modeling by modal alignment, order
reduction, and interpolations is challenging for large-dimension systems,
usually resulted from high-fidelity finite-element modeling. However, the
proposed method avoids these steps and derives parametric models directly from
data snapshots, which makes p-DMD promising to extract ROMs from high-fidelity
data.

  3. 3.

The physical understanding of parameter dependency is encoded into the
parametric Hankel matrices of the data-driven LPV representation (see Eqns.
(10) and (15)). The parametric ROM is directly obtained over a range of
scheduling parameters without interpolations of grid models. Last, the state
consistency is inherent under varying flight conditions because the same
projection eigenvectors are used (see Eqn. (19)).

In the rest of this paper, the traditional method of linearization-
interpolation approach of aeroelastic modeling is briefly reviewed in Section
II. The data-driven reduced-order LPV modeling is then derived and described
in Section III. After that, a numerical study and results will be provided to
demonstrate the effectiveness of the proposed method in Sections IV and V. At
last, conclusions are made in Section VI.

##  II Nonlinear and Linearized Equations of Motion

For control system development of highly flexible aircraft, one must follow a
systematic approach based on a nonlinear aeroelastic formulation. On top of
that, the authors have developed reduced-order LPV models for highly flexible
aircraft. This section briefly reviews the fundamentals of nonlinear
aeroelastic formulation and LPV modeling by linearization and interpolation at
equilibrium points.

There have been numerous endeavors that examine the nonlinear aeroelasticity
of highly flexible aircraft. While more refined models based on shell
formulations may capture more details of deformations in morphing and
geometrically complicated wings, a geometrically nonlinear beam theory is
usually sufficient to establish an aeroelastic formulation for highly flexible
aircraft. Palacios et al. [18] comprehensively discussed three types of beam
formulations for the structural, aeroelastic, and flight dynamic analysis of
highly flexible aircraft. By coupling a strain-based geometrically nonlinear
beam formulation with the finite-state inflow theory [19, 20, 21], Cesnik and
his colleagues established the nonlinear aeroelastic and flight dynamic
formulation, now known as UM/NAST [22, 23, 24].

By following the strain-based approach [23, 24], the coupled nonlinear
aeroelasticity and flight dynamics of highly flexible aircraft (see Fig. 1) is
governed by

$\begin{gathered}\mathbf{M}_{FF}(\bm{\varepsilon})\ddot{\bm{\varepsilon}}+\mathbf{M}_{FB}(\bm{\varepsilon})\dot{\bm{\beta}}+\mathbf{K}_{FF}\bm{\varepsilon}=\mathbf{R}_{F}\\\
\mathbf{M}_{BF}(\bm{\varepsilon})\ddot{\bm{\varepsilon}}+\mathbf{M}_{BB}(\bm{\varepsilon})\dot{\bm{\beta}}=\mathbf{R}_{B}\\\
\dot{\bm{\zeta}}=-\frac{1}{2}\bm{\Omega}_{\bm{\zeta}}(\bm{\beta})\bm{\zeta}\\\
\dot{\mathbf{p}}_{B}^{G}=\begin{bmatrix}\mathbf{C}^{GB}(\bm{\zeta})&\mathbf{0}\end{bmatrix}\bm{\beta}\\\
\dot{\bm{\mathbf{\lambda}}}=\mathbf{F}_{1}(\dot{\bm{\varepsilon}},\bm{\beta})\bm{\lambda}+\mathbf{F}_{2}(\dot{\bm{\varepsilon}},\bm{\beta})\begin{Bmatrix}\dot{\bm{\varepsilon}}\\\
\bm{\beta}\\\
\end{Bmatrix}+\mathbf{F}_{3}\begin{Bmatrix}\ddot{\bm{\varepsilon}}\\\
\dot{\bm{\beta}}\\\ \end{Bmatrix}\\\ \end{gathered}$ (1)

where the independent variables are the elastic strain along the flexible
members $\bm{\varepsilon}$ and the aircraft rigid-body velocity $\bm{\beta}$.
The rigid body’s propagation is further described by the quaternion
$\bm{\zeta}$ and inertial position $\mathbf{p}_{B}^{G}$. $\bm{\lambda}$ is the
inflow state [19, 20, 21] from the unsteady aerodynamics. Finally,
$\bm{\delta}$ consists of the deflections of trailing-edge control surfaces:
elevator, aileron, and rudder, respectively. The load vector consists of the
contributions of initial strain, damping, gravity, aerodynamics, and
propulsion (thrust force), i.e.,

$\begin{array}[]{c}\begin{Bmatrix}\mathbf{R}_{F}\\\ \mathbf{R}_{B}\\\
\end{Bmatrix}=\begin{Bmatrix}\mathbf{K}_{FF}\bm{\varepsilon}^{0}\\\
\mathbf{0}\\\
\end{Bmatrix}+\begin{Bmatrix}\mathbf{R}_{F}^{\textrm{damp}}(\bm{\varepsilon},\dot{\bm{\varepsilon}},\bm{\beta})\\\
\mathbf{R}_{B}^{\textrm{damp}}(\bm{\varepsilon},\dot{\bm{\varepsilon}},\bm{\beta})\\\
\end{Bmatrix}+\begin{Bmatrix}\mathbf{R}_{F}^{\textrm{grav}}(\bm{\zeta})\\\
\mathbf{R}_{B}^{\textrm{grav}}(\bm{\zeta})\\\ \end{Bmatrix}\\\
+\begin{Bmatrix}\mathbf{R}_{F}^{\textrm{aero}}(\bm{\varepsilon},\dot{\bm{\varepsilon}},\ddot{\bm{\varepsilon}},\bm{\beta},\dot{\bm{\beta}},\bm{\lambda},\bm{\delta})\\\
\mathbf{R}_{B}^{\textrm{aero}}(\bm{\varepsilon},\dot{\bm{\varepsilon}},\ddot{\bm{\varepsilon}},\bm{\beta},\dot{\bm{\beta}},\bm{\lambda},\bm{\delta})\\\
\end{Bmatrix}+\begin{Bmatrix}\mathbf{R}_{F}^{\textrm{prop}}(T)\\\
\mathbf{R}_{B}^{\textrm{prop}}(T)\\\ \end{Bmatrix}\,.\end{array}$ (2)

The detailed expressions are elaborated in [23, 24, 3].

Figure 1: Flexible wing with large, nonlinear deformation.

The traditional linearization-interpolation approach to derive an LPV model is
summarized here. First, the nonlinear dynamical equation is linearized at a
fixed equilibrium flight condition [23, 3]. The second step is to perform the
modal coordinate transformation on the high-dimensional full-order models
(FOMs) at grid scheduling parameters. As a result, a series of LTI-FOMs in
modal coordinates are obtained. The third step is to align the vibration modes
to depict the modal evolution so that model reduction can be performed by
keeping the most significant modes in terms of $H_{2}$ or $H_{\infty}$ norm.
In the final step, the evolution of modal LTI-ROMs is interpolated into a
polynomial parameter-dependent LPV-ROM [6]. An LPV controller and smooth
switching LPV controller can be further developed against gust disturbance and
model uncertainty [4].

The parametric state-space equation reads as

$\displaystyle\dot{\mathbf{x}}$
$\displaystyle=\mathbf{A}(\theta)\mathbf{x}+\mathbf{B}(\theta)\mathbf{u}\hfill$
(3) $\displaystyle\mathbf{y}$
$\displaystyle=\mathbf{C}(\theta)\mathbf{x}+\mathbf{D}(\theta)\mathbf{u}\hfill$

where the state variable and control input are

$\displaystyle\mathbf{x}^{T}$
$\displaystyle=\begin{Bmatrix}\Delta\bm{\varepsilon}^{T}&\Delta\dot{\bm{\varepsilon}}^{T}&\Delta\bm{\beta}^{T}&\Delta\bm{\zeta}^{T}&\left(\Delta\mathbf{p}_{B}^{G}\right)^{T}&\Delta\bm{\lambda}^{T}\\\
\end{Bmatrix}\hfill$ (4) $\displaystyle\mathbf{u}^{T}$
$\displaystyle=\begin{Bmatrix}\Delta\bm{\delta}^{T}&\Delta T\\\
\end{Bmatrix}\,.\hfill$

The scheduling parameter $\theta$ is usually a measurable quantity relevant to
the varying flight condition, and output matrices $\mathbf{C},\mathbf{D}$ are
selected according to the user’s need. The discrete-time model can be further
converted from continuous-time model (3).

This approach realizes the modeling of interactions among rigid-body dynamics,
structural dynamics, and aerodynamics. The coupling affects the equivalent
stiffness and inertial matrices under varying flight conditions. The model-
based approach has several characteristics. First, this approach produces
physically explainable models where the states of strains, rigid-body
dynamics, and aerodynamics are aligned and analyzed. The state consistency
across flight regimes is guaranteed. Second, the local stability of the
coupled system can be easily analyzed by the linearized system matrices.

However, there are several limitations in the linearization-based approach: 1)
The linearization and model reduction process usually need heavy computations
if refined meshing is conducted at structures and aerodynamics nodes. Modal
alignment to enable state consistency is difficult for high-dimensional
systems, which is especially critical in high-fidelity modeling and unsteady
aerodynamics. 2) The partial derivative of aerodynamic loading
($\mathbf{R}_{F}^{\textrm{aero}}$) relative to strains, strain rates need to
be known. In the low-fidelity case, a complex analytical function can lead to
close-form partial derivatives. However, the derivatives are hard to derive in
high-fidelity and complex-geometry cases. 3) The biggest drawback of the
linearization-based approach is that the deformation is assumed to be small
around equilibrium. With geometrically nonlinear deformation or at transient
flight conditions, the linearization approach cannot be applied. The LPV model
from interpolation neglects the variation rate of the scheduling parameter and
may cause pitfalls of the gain-scheduling controller.

Therefore, in this study, we propose an innovative data-driven modeling
approach for highly flexible aircraft that is capable of describing the
geometrically nonlinear deformation in a dynamic flight condition. The input
of the proposed algorithm is the data snapshots of states, scheduling
parameters (flight conditions), and control inputs of persistent excitation.
The outcome of the algorithm is a polynomial-dependent LPV-ROM. Users can
select the dimension of the LPV-ROM and the polynomial order to describe the
nonlinear aeroelasticity based on the need for accuracy levels. The data
snapshots are not required to be collected at equilibrium points, but can also
be collected at varying flight conditions.

##  III Reduced-Order LPV Modeling by p-DMD

###  III-A Problem Statement

Based on Eq. (1), one can form a generic discrete-time nonlinear parameter
varying model Eq. (5) to describe the flight dynamic and aeroelastic behavior
of a highly flexible aircraft.

$\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\mathbf{f}\left(\mathbf{x}_{k},\mathbf{u}_{k},\mathbf{\theta}_{k}\right)$
(5) $\displaystyle\mathbf{y}_{k}$ $\displaystyle=\mathbf{C}\mathbf{x}_{k}$

where $\mathbf{x}_{k}\in\mathbf{R}^{n_{x}}$,
$\mathbf{u}_{k}\in\mathbf{R}^{n_{u}}$, $\mathbf{y}_{k}\in\mathbf{R}^{n_{y}}$
are the states, control inputs, and outputs, with $k$ being the time index.
$\mathbf{\theta}_{k}$ is a vector of time-varying scheduling parameters of
flight conditions, which are assumed to be known in real-time. The output
matrix $\mathbf{C}$ can be properly selected to obtain some of the states as
outputs, for example, strains or rigid-body velocities, which renders a
constant $\mathbf{C}$ and $\mathbf{D}=0$.

The aircraft’s aeroelastic behavior may vary with flight conditions, such as
flight speed and angle of attack. This study considers only a single
parameter, but it can be easily extended to a multiple-dimensional case. From
the nonlinear model (1), all states (strains, strain rates, quaternions,
rigid-body motions, aerodynamic states) are included in the parametric FOM.
Over a range of flight conditions $\theta_{k}\in\mathcal{P}\in\mathbf{R}$, the
full-order nonlinear model (5) is approximated by a parametric FOM (6), and
the matrices $\mathbf{A}(\mathbf{\theta}_{k})\in\mathbf{R}^{n_{x}\times
n_{x}},\mathbf{B}(\mathbf{\theta}_{k})\in\mathbf{R}^{n_{x}\times n_{u}}$ are
assumed as $n_{p}$-order polynomials dependent on $\theta_{k}$ in (7).

$\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\mathbf{A}\left(\mathbf{\theta}_{k}\right)\mathbf{x}_{k}+\mathbf{B}\left(\mathbf{\theta}_{k}\right)\mathbf{u}_{k}$
(6) $\displaystyle\mathbf{y}_{k}$ $\displaystyle=\mathbf{C}\mathbf{x}_{k}$

$\mathbf{A}\left(\mathbf{\theta}_{k}\right)=\mathbf{A}_{0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{A}_{i},\quad\mathbf{B}\left(\mathbf{\theta}_{k}\right)=\mathbf{B}_{0}+\sum_{i=1}^{n_{\mathrm{p}}}(\mathbf{\theta}_{k})^{i}\mathbf{B}_{i}.$
(7)

The technical problem is to find a reduced-order LPV state-space model (8) to
approximate the nonlinear aeroelastic model, where matrices
$\mathbf{A}_{r}(\theta_{k})\in\mathbf{R}^{n_{z}\times
n_{z}},\mathbf{B}_{r}(\theta_{k})\in\mathbf{R}^{n_{z}\times n_{u}}$, $n_{z}\ll
n_{x}$, and $\mathbf{z}_{k}$ is the state vector of LPV-ROM.

$\displaystyle\mathbf{z}_{k+1}$
$\displaystyle=\mathbf{A}_{r}\left(\theta_{k}\right)\mathbf{z}_{k}+\mathbf{B}_{r}\left(\theta_{k}\right)\mathbf{u}_{k}$
(8) $\displaystyle\mathbf{y}_{k}$ $\displaystyle=\mathbf{C}_{r}\mathbf{z}_{k}$

$\mathbf{A}_{r}\left(\theta_{k}\right)=\mathbf{A}_{r0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{A}_{ri},\quad\mathbf{B}_{r}\left(\theta_{k}\right)=\mathbf{B}_{r0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{B}_{ri}.$
(9)

###  III-B Parametric Data-Driven Representations

Suppose the data snapshots are obtained as
$\left\\{\mathbf{u}_{k},\mathbf{x}_{k},\theta_{k}\right\\}_{N_{d}}$ during
time steps $k=1$ to $N_{d}$. The parametric LPV model (6), (7) can be
equivalently expressed by (10) with parametric Hankel matrices.

$\begin{array}[]{rl}\mathbf{x}_{k+1}=&(\mathbf{A}_{0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{A}_{i})x_{k}+(\mathbf{B}_{0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{B}_{i})\mathbf{u}_{k}\\\
=&\underbrace{\left[\begin{array}[]{cccc}\mathbf{A}_{0}&\mathbf{A}_{1}&\cdots&\mathbf{A}_{n_{\mathrm{p}}}\end{array}\right]}_{\mathcal{A}}\left[\begin{array}[]{c}\mathbf{x}_{k}\\\
(\theta_{k})^{1}\otimes\mathbf{x}_{k}\\\ \vdots\\\
(\theta_{k})^{n_{\mathrm{p}}}\otimes\mathbf{x}_{k}\end{array}\right]\\\
+&\underbrace{\left[\begin{array}[]{cccc}\mathbf{B}_{0}&\mathbf{B}_{1}&\cdots&\mathbf{B}_{n_{\mathrm{p}}}\end{array}\right]}_{\mathcal{B}}\left[\begin{array}[]{c}\mathbf{u}_{k}\\\
(\theta_{k})^{1}\otimes\mathbf{u}_{k}\\\ \vdots\\\
(\theta_{k})^{n_{\mathrm{p}}}\otimes\mathbf{u}_{k}\end{array}\right]\\\
=&\mathcal{A}\underbrace{\left[\begin{array}[]{c}\mathbf{x}_{k}\\\
\Theta_{k}\otimes\mathbf{x}_{k}\end{array}\right]}_{\Omega_{\mathbf{x}}(\mathbf{\theta}_{k})}+\mathcal{B}\underbrace{\left[\begin{array}[]{c}\mathbf{u}_{k}\\\
\Theta_{k}\otimes\mathbf{u}_{k}\end{array}\right]}_{\Omega_{\mathbf{u}}(\mathbf{\theta}_{k})}\\\
\end{array}$ (10)

where,
$\Theta_{k}=\left[(\theta_{k})^{1},\cdots,(\theta_{k})^{n_{\mathrm{p}}}\right]^{T}$.

######  Remark 1 (parametric Hankel matrix).

The parametric data-driven representation of a nonlinear parameter-varying
system has additional polynomial coefficient matrices
$\mathbf{A}_{1},\cdots,\mathbf{A}_{n_{\mathrm{p}}},\mathbf{B}_{1},\cdots,\mathbf{B}_{n_{\mathrm{p}}}$
other than $\mathbf{A}_{0},\mathbf{B}_{0}$. The data snapshot formulates
parametric Hankel matrices $\Omega_{\mathbf{x}}(\mathbf{\theta}_{k})$,
$\Omega_{\mathbf{u}}(\mathbf{\theta}_{k})$, which encodes the physical
understanding of the nonlinear aeroelasticity by the Kronecker product kernel.

Note that the dimension of $\Theta_{k}$ is $n_{p}$, and the dimension of
$\mathbf{x}_{k}$ is $n_{x}$. The following matrices can be constructed by
stacking data snapshots of the sampling window from 1 to $N_{d}$.

$\displaystyle\mathbf{U}=\left[\begin{array}[]{lll}\mathbf{u}_{1},&\cdots,&\mathbf{u}_{N_{\mathrm{d}}}\end{array}\right],\mathbf{U}^{\mathrm{kr}}=\left[\Theta_{1}\otimes\mathbf{u}_{1},\cdots,\Theta_{N_{\mathrm{d}}}\otimes\mathbf{u}_{N_{\mathrm{d}}}\right],$
(11)
$\displaystyle\mathbf{X}=\left[\begin{array}[]{lll}\mathbf{x}_{1},&\cdots,&\mathbf{x}_{N_{\mathrm{d}}}\end{array}\right],\mathbf{X}^{\mathrm{kr}}=\left[\Theta_{1}\otimes\mathbf{x}_{1},\cdots,\Theta_{N_{\mathrm{d}}}\otimes\mathbf{x}_{N_{\mathrm{d}}}\right],$
$\displaystyle\mathbf{X}^{+}=\left[\begin{array}[]{lll}\mathbf{x}_{2}&\cdots&\mathbf{x}_{N_{\mathrm{d}}+1}\end{array}\right].$

$\begin{array}[]{rl}\mathbf{X}^{+}&=\left[\begin{array}[]{ccc}\mathbf{x}_{2}&\cdots&\mathbf{x}_{N_{\mathrm{d}}+1}\end{array}\right]\\\
&=\mathcal{A}\left[\begin{array}[]{ccccc}\begin{array}[]{c}\mathbf{x}_{1}\\\
\Theta_{1}\otimes\mathbf{x}_{1}\end{array}&\cdots\begin{array}[]{c}\mathbf{x}_{N_{d}}\\\
\Theta_{N_{d}}\otimes\mathbf{x}_{N_{d}}\end{array}\end{array}\right]+\mathcal{B}\left[\begin{array}[]{ccc}\begin{array}[]{c}\mathbf{u}_{1}\\\
\Theta_{1}\otimes\mathbf{u}_{1}\end{array}&\cdots\begin{array}[]{c}\mathbf{u}_{N_{d}}\\\
\Theta_{N_{d}}\otimes\mathbf{u}_{N_{d}}\end{array}\end{array}\right].\end{array}$
(12)

Therefore, we have the full-order state evolution represented by the finite-
horizon data snapshot,

$\mathbf{X}^{+}=\mathcal{A}\left[\begin{array}[]{c}\mathbf{X}\\\
\mathbf{X}^{\mathrm{kr}}\end{array}\right]+\mathcal{B}\left[\begin{array}[]{c}\mathbf{U}\\\
\mathbf{U}^{\mathrm{kr}}\end{array}\right]=\left[\mathcal{A}\quad\mathcal{B}\right]\left[\begin{array}[]{c}\mathbf{X}\\\
\mathbf{X}^{\mathrm{kr}}\\\ \mathbf{U}\\\
\mathbf{U}^{\mathrm{kr}}\end{array}\right].$ (13)

Based on this, the full-order parametric matrices of the polynomial-dependent
LPV model can be derived as

$\left[\mathcal{A}\quad\mathcal{B}\right]=\mathbf{X}^{+}\left[\begin{array}[]{c}\mathbf{X}\\\
\mathbf{X}^{\mathrm{kr}}\\\ \mathbf{U}\\\
\mathbf{U}^{\mathrm{kr}}\end{array}\right]^{\dagger}$ (14)

where, $\dagger$ denotes pseudo-inverse. Therefore, the full-order model (6)
and LPV full-order model can be denoted by a data-driven parametric
representation in (15).

$\mathbf{x}_{k+1}=\mathbf{X}^{+}\left[\begin{array}[]{c}\mathbf{X}\\\
\mathbf{X}^{\mathrm{kr}}\\\ \mathbf{U}\\\
\mathbf{U}^{\mathrm{kr}}\end{array}\right]^{\dagger}\left[\begin{array}[]{c}\mathbf{x}_{k}\\\
\Theta_{k}\otimes\mathbf{x}_{k}\\\ \mathbf{u}_{k}\\\
\Theta_{k}\otimes\mathbf{u}_{k}\end{array}\right]$ (15)

###  III-C Model Order Reduction by p-DMD

Taking the singular value decomposition (SVD) on the data snapshots, one can
obtain the model approximation by keeping the $n_{z}$ numbers of most
significant singular values in $\Sigma_{r}$,

$\left[\mathbf{X};\mathbf{X}^{\mathrm{kr}};\mathbf{U};\mathbf{U}^{\mathrm{kr}}\right]=\mathbf{U}\Sigma\mathbf{V}^{T}\approx\mathbf{U}_{r}\Sigma_{r}\mathbf{V}_{r}^{T}.$
(16)

Therefore, the full-order parametric matrices can be represented by

$\left[\mathcal{A},\mathcal{B}\right]=\mathbf{X}^{+}\mathbf{V}_{r}\Sigma_{r}^{-1}\mathbf{U}_{r}^{T}$
(17)

From the SVD of one-step shifted states sequence, the projection basis can be
obtained by extracting the dominant singular values and the corresponding
orthonormal vectors $\hat{U}_{z}$.

$\mathbf{X}^{+}=\hat{\mathbf{U}}\hat{\Sigma}\hat{\mathbf{V}}^{T}\approx\hat{\mathbf{U}}_{z}\hat{\Sigma}_{z}\hat{\mathbf{V}}_{z}^{T}$
(18)

Therefore, the reduced-order matrices are obtained from the projection basis
from $\hat{U}_{z}$ is

$\left[\mathcal{A}_{r},\mathcal{B}_{r},\mathbf{C}_{r}\right]=\left[\hat{\mathbf{U}}_{z}^{T}\mathcal{A}\hat{\mathbf{U}}_{z},\hat{\mathbf{U}}_{z}^{T}\mathcal{B},\mathbf{C}\hat{\mathbf{U}}_{z}\right].$
(19)

where,
$\mathcal{A}_{r}=\left[\mathbf{A}_{r0},\mathbf{A}_{r1},\cdots,\mathbf{A}_{rn_{p}}\right],\mathcal{B}_{r}=\left[\mathbf{B}_{r0},\mathbf{B}_{r1},\cdots,\mathbf{B}_{rn_{p}}\right]$.
With the obtained matrices $\mathcal{A}_{r},\mathcal{B}_{r}$, the reduced-
order parameter-dependent system matrices can be extracted as

$\mathbf{A}_{r}\left(\theta_{k}\right)=\mathbf{A}_{r0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{A}_{ri},\quad\mathbf{B}_{r}\left(\theta_{k}\right)=\mathbf{B}_{r0}+\sum_{i=1}^{n_{\mathrm{p}}}(\theta_{k})^{i}\mathbf{B}_{ri}.$
(20)

######  Remark 2.

The data snapshots
$\left(\mathbf{U},\mathbf{U}^{kr},\mathbf{X},\mathbf{X}^{kr},\mathbf{X}^{+}\right)$
are not required to be around the steady-state (equilibrium) flight condition
and small deformation. Therefore, p-DMD has improved capability than
traditional DMD [16, 5] in terms of varying flight conditions (scheduling
parameters) and large structure deformation. Besides, the transient dynamics
caused by varying scheduling parameters can be captured by p-DMD as long as
the data snapshot is obtained at dynamic flight conditions, which is studied
and verified in Section V.

######  Remark 3 (state consistency).

A constant projection vector $\hat{\mathbf{U}}_{z}$ is used for order
reduction projection from full-order states, so state consistency is
inherently guaranteed. The proposed method avoids performing DMD at multiple
operating points and interpolating, so eigen-tracking in interpolation to
obtain state consistency is not needed. On the contrary, the traditional DMD
cannot always render system matrices with state consistency across multiple
fixed scheduling parameters unless the same projection basis can be
successfully found and implemented [5, 16]. Besides, the projection basis is
formulated from snapshot of states, so the LPV-ROM’s states are physically
interpretative as the dominant behaviors of states.

##  IV Numerical Study Case 1: Highly Flexible Cantilever Wing

The highly flexible rectangular wing provided in [25] is studied for the
current study. The physical and geometrical properties of the wing are listed
in Table I. It is divided into ten elements with a clamped root. A trailing-
edge control surface on the wing is defined from its $60\text{\,}\mathrm{\char
37\relax}$ to $90\text{\,}\mathrm{\char 37\relax}$ span, occupying
$20\text{\,}\mathrm{\char 37\relax}$ of the chord.

TABLE I: Properties of a highly flexible wing Quantity | Value
---|---
Length |  16
Chord |  1
Spanwise ref. axis location (from L.E.) |  50 of chord
Center of gravity (from L.E.) |  50 of chord
Torsional rigidity |  1e4
Flat bending rigidity |  2e4
Chord bending rigidity |  4e6
Mass per unit length |  0.75
Rotational Inertia per unit length |  0.1

###  IV-A Fixed Scheduling Parameter (Angle of Attack)

The flexible wing is placed in a freestream of $U_{\infty}=\qty{20}{\per}$ at
the altitude of 20000. The angle of attack measured at the wing root is
selected as the scheduling parameter for LPV models. The models span at grid
angles of attack ranging from $-10^{\circ}$ to $10^{\circ}$ with an increment
of $0.5^{\circ}$, resulting in 41 grid models. The states of the full-order
model include four strains, four strain rates, and six aerodynamic states for
each beam element. Therefore, the full-order state-space model has a dimension
of 140. The outputs are selected as the out-of-plane bending strains of all
ten elements. Since it is cantilevered at the root, the rigid-body dynamics
are not included in the study of the slender wing.

The setup of simulations is as follows to obtain the data snapshots at an
equilibrium point of fixed angle of attack. A chirp signal in Fig. 2 is
applied to the main-wing flap at one fixed angle of attack to excite the
aeroelastic modes, and the data snapshot of states and inputs are collected to
conduct the data-driven LPV-ROM modeling algorithm. The chirp signal consists
of frequency components spanning from $0.1$ Hz to $10$ Hz. The full-order
linearized aeroelastic system is simulated at the sampling rate of
$T_{s}=\qty{0.001}{}$.

Figure 2: Excitation input of a chirp signal with sweep frequency from 0.1-10
Hz.

####  IV-A1 Model order versus modeling accuracy

The p-DMD method is repeatedly applied at grid angles of attack, and the order
of ROMs versus modeling accuracy is explored. The first 20 most significant
singular values $\Sigma$ are plotted in Fig. 3. It can be seen the first 12
singular values can retain over $95\%$ of the summation of all singular
values, and is selected as the order of ROMs.

The $v$-gap metrics between the ROMs and FOMs at fixed parameters are plotted
in Fig. 4 to evaluate the model accuracy of the ROMs. The gap metric indicates
the distances between two systems and ranges between 0 and 1. A value close to
0 means that the FOM and ROM are close, and the stabilizing controller of ROM
can stabilize FOM. Fig. 4 shows the maximum gap metric occurs around $0.1$ Hz,
and has a very small value of $10^{-5}$ level. Therefore, the LPV-ROM model of
order $12$ can accurately capture the full-order model in the selected
frequency range.

At the angles of attack $\alpha$ equals $None$ and $None$, respectively, the
time- and frequency-domain responses of linearized FOM, 10-order LPV-ROM, and
12-order LPV-ROM are compared in Figs. 5-8. Note that, the plotted variables
are the perturbed magnitudes from the trim points. It can be easily observed
that the 10-order LPV-ROM has an obvious modeling error in the frequency range
of concern, while the 12-order LPV-ROM model can capture the FOM in the
frequency range of 040rad/s ($\approx 6.3$ ). Therefore, the 12-order LPV-ROM
is the “minimum order” model that can accurately describe the system dynamics.
The polynomial order is selected as 4 to formulate LPV-ROM.

Figure 3: The 20 largest singular values of grid models at fixed scheduling
parameters. Figure 4: $v$-gap metric surface of frequency responses at fixed
models. Gap metric close to 0 meaning that two systems are close.

Figure 5: Frequency responses of $u$ to $y_{1}$ (top), $y_{5}$ (middle) and
$y_{10}$ (bottom), with 10th, 12th-order ROMs at angle of attack
$\alpha=-10^{\circ}$.

Figure 6: Time-domain responses of $y_{1}$ (top), $y_{5}$ (middle) and
$y_{10}$ (bottom), with 10th, 12th-order ROMs at angle of attack
$\alpha=-10^{\circ}$.

Figure 7: Frequency responses of $u$ to $y_{1}$ (top), $y_{5}$ (middle) and
$y_{10}$ (bottom), with 10th, 12th-order ROMs at angle of attack
$\alpha=0^{\circ}$.

Figure 8: Time-domain responses of $y_{1}$ (top), $y_{5}$ (middle) and
$y_{10}$ (bottom), with 10th, 12th-order ROMs at angle of attack
$\alpha=0^{\circ}$.

####  IV-A2 Verification between p-DMD modeling and linearization-based
modeling

In this section, the results from the authors’ previous work of linearization-
based modeling are compared with the data-driven p-DMD modeling. The nonlinear
model (1) is linearized and reduced order at grid angles of attack to obtain a
series of LTI-FOM. The root loci of the ROM are shown in Fig. 9.

For reasonable comparisons, the equilibrium conditions at multiple angles of
attacks are focused. The data snapshot are obtained by simulating the
aeroelastic model at small perturbations of steady state. The comparisons are
plotted in Figs. 10 and 11. The data-driven LPV-ROM is found to match the
linearization-based LPV-ROM, verifying that the data-driven approach can
capture the physical characteristics of the aeroelasticity.

It is worth mentioning the scenarios of $\alpha=-10^{\circ},+10^{\circ}$, the
eigenvalues have a positive real part (see root loci in Fig. 9), indicating
that the aeroelastic system is unstable at these angles of attack. In the
linearized-based approach, the $\sigma$-shifted $H_{2}$ norm is calculated by
shifting the eigenvalue by $\sigma$ to obtain a stable system [6]. However,
the proposed p-DMD directly finds the dominant modes by analyzing the singular
values of parametric Hankel matrices. The p-DMD approach can address the
unstable dynamics.

To conclude this section, the p-DMD method can accurately capture the local
dynamics of the highly flexible wing at a fixed angle of attack, and the model
accuracy is verified with the well-established linearization-based approach in
frequency- and time-domain simulations.

Figure 9: Root loci of ROMs of slender wing at grid points of angle of attack.

Figure 10: Frequency and time-domain responses comparison of data-driven ROMs
and linearization-based ROMs at $\alpha=-10^{\circ}$.

Figure 11: Frequency and time-domain responses comparison of data-driven ROMs
and linearization-based ROMs at $\alpha=0^{\circ}$.

##  V Numerical Study Case 2: Highly Flexible Aircraft

In this section, we study and show the results of data-driven modeling for a
highly flexible aircraft with coupled elastic deformation and rigid-body
motions. The aircraft is built based on the slender wing in Section IV, as
shown shown in Fig. 12. All members of the aircraft are modeled as slender
beams. However, the body and tail members are 100 times stiffer than the main
wing, although they all have the same inertial properties. The chord length of
the horizontal and vertical tails is 0.5, while their beam reference lines are
defined at a quarter of the chord from the leading edge. There is no
aerodynamics coupled on the body. Similar to the main wing, trailing-edge
control surfaces are also defined on the horizontal tail, from its 25 to the
full span. The flaps on the horizontal tail are used as elevators to trim the
vehicle for level flight and excite the longitudinal motion.

Figure 12: Platform of flexible vehicle.

###  V-A Fixed Scheduling Parameter (Angle of Attack)

In this simulation, the highly flexible aircraft is set cruising in a
freestream of $U_{\infty}=\qty{22}{\per}$ at the altitude of 20000. For LPV
modeling, the body angle of attack at the origin of the body frame is used as
the scheduling parameter. The LPV models span at grid angles of attack ranging
from $-7.5^{\circ}$ to $7.5^{\circ}$ with an increment of $0.5^{\circ}$. The
total number of grid models is 31. The flexible vehicle model has 34 beam
elements and four control inputs (two control surfaces of the main wing, one
elevator, and one thrust). The symmetric wing vibrations and longitudinal
rigid-body motion are considered. Instead of being used as ailerons, the
surfaces on the main wing deflect symmetrically when excited. The FOM includes
four strains, four strain rates, and six aerodynamic states for each beam
element, in addition to six rigid-body velocities, four quaternions, and three
inertial positions of the vehicle. Therefore, at a fixed angle of attack, the
full-order $A$ matrix is of dimension $471\times 471$, the $B$ matrix is of
dimension $471\times 3$, and the 26 outputs are selected as the out-of-plane
bending strains on all main-wing elements (20 in total) plus longitudinal and
vertical body velocities, pitch angle and rate, longitudinal and vertical body
positions.

At a fixed angle of attack, the chirp signal in Fig. 2 is fed as input
perturbation of the main-wing flaps to excite the aeroelastic modes, and the
data snapshots of states and inputs are collected to conduct the data-driven
LPV-ROM modeling algorithm. The chirp signal consists of frequency components
ranging from 0.110. The full-order linearized aeroelastic system is also
simulated at the sampling rate of $T_{s}=\qty{0.001}{}$.

The 20 largest singular values of the state-Hankel matrix are plotted in Fig.
13, among which the 12 largest singular values are determined to be able to
capture over $95\%$ accuracy of FOMs. The root loci of flexible vehicles at
varying angles of attack is Fig. 14, which is obtained from model reduction
from FOMs. The v-gap metric comparison between the LTI-ROMs and LPV-ROMs by
applying p-DMD on flexible vehicles is plotted in Fig. 15. These gap metrics
are close to 0 at most of the regime, indicating the data-driven LPV-ROM can
capture the local dynamics with excellent accuracy.

The time-domain responses of data-driven LPV-ROM and full-order linearized
model are compared in Fig. 16. The strain of wing root, longitudinal velocity,
and pitch rate all match with an error less than $2\%$. Note that at
$\alpha=-7.5,7.5^{\circ}$, the aeroelastic model becomes unstable (see Fig.
14, one pair of eigenvalues have positive real parts), but the p-DMD can still
render an accurate data-driven model. The frequency-domain response
comparisons are omitted for the sake of space limit.

Figure 13: The 20 largest singular values of flexible aircraft at fixed
scheduling parameters. Figure 14: Root loci of ROMs of flexible aircraft at
grid angles of attack. Figure 15: $v$-gap metric surface of frequency
responses at fixed models. Gap metric close to 0 meaning that two systems are
similar.

Figure 16: Time-domain responses of LPV-ROM and FOM at angle of attack
$\alpha=-7.5,+7.5^{\circ}$. (a)-(b): Bending strains at wing root; (c)-(d):
Longitudinal velocity; (e)-(f): Pitch rate.

###  V-B Varying Scheduling Parameter (Angle of Attack)

In this section, the capability to model aeroelasticity at varying angles of
attack and non-equilibrium conditions is demonstrated.

The nonlinear equation of motion Eq. (1) is numerically integrated in the time
domain to obtain the vehicle’s nonlinear transient behavior. The nominal
flight speed is still 22 at the altitude of 20000. The simulation starts with
a static operating condition, where the body pitch angle is $\theta_{B}=$$$,
elevator deflection is $\delta_{e}=$$$, and thrust is $T=\qty{51.07}{}$. A
chirp signal with a magnitude of $5^{\circ}$, sweep frequency $0.01$-$2$ Hz is
applied to the elevator to excite the rigid-body pitch motions and aeroelastic
dynamics. The simulation is conducted with a sampling time of $0.005$s, and
the full-order nonlinear model is simulated within a $10$-second period. All
the state and input data are collected. The scheduling parameter (angle of
attack) is calculated in real-time by $\alpha=-\arcsin(\beta_{3}/U_{\infty})$.

The collected data snapshot is applied with the p-DMD to obtain the data-
driven LPV-ROM. The 20 largest singular values of the parametric Hankel matrix
are plotted in Fig. 17. The data-driven LPV-ROM is chosen to have the
polynomial order of 4, and dimension of 12.

Figure 17: Square roots of 20 largest singular values of the parametric Hankel
matrix.

After that, the model accuracy of LPV-ROM is verified using other control
inputs. Another chirp signal with a magnitude of $5^{\circ}$ and sweep
frequency $0.1$-$2$ Hz is fed into the full-order nonlinear model and LPV-ROM,
and their responses are compared. Fig. 18 shows the vehicle’s time-domain
elastic and rigid-body responses. Two models yield reasonably similar
responses with a model error of $2.7\%$. From the results, it is worth
emphasizing that the bending strains and rigid-body velocities start from a
non-zero initial condition and remain at a non-equilibrium flight condition.
This demonstrates that the proposed p-DMD can address the non-equilibrium
conditions. In addition, the angle of attack calculated from the two models’
responses matches accurately in Fig. 18, which continuously varies between
$0.5^{\circ}$-$3.5^{\circ}$. This study demonstrates that the proposed p-DMD
can address the case of varying scheduling parameters. As stated in the
novelties of this work, since no interpolation and linearization are used in
the data-driven approach, the proposed method can derive parametric ROMs for
flexible aircraft with large wing deformations and unsteady dynamics.

Figure 18: Time-domain responses of LPV-ROM and nonlinear FOM at non-
equilibrium flight condition.

##  VI Conclusion

This paper presented a data-driven approach for generating accurate linear
parameter-varying reduced-order models (LPV-ROM) regarding the nonlinear
aeroelasticity and flight dynamics of highly flexible aircraft. This method
encoded the physical understanding of the nonlinear dynamics of highly
flexible aircraft into a parametric Hankel matrix and performed parametric
dynamic mode decomposition (p-DMD) to render LPV-ROMs. Both the fixed and
varying scheduling parameter scenarios could be addressed by the p-DMD. The
method was applied to and verified with a highly flexible cantilever wing and
a complete free-flight aircraft. The local modeling accuracy with fixed
scheduling parameters was verified by cross-validation with the linearized
model. The global modeling accuracy was verified by comparing it with the
nonlinear full-order model regarding the time-domain responses in transient
flight conditions. This method is promising to use data from high-fidelity
analysis to directly develop accurate models for flexible aircraft.

## References

  * [1] M. C. van Shoor, S. H. Zerweckh, and A. H. von Flotow, “Aeroelastic stability and control of a highly flexible aircraft,” in Proceedings of the 30th AIAA Structures, Structural Dynamics and Materials Conference, (Mobile, AL), Apr. 1989.  AIAA-1989-1187.
  * [2] J. M. Barker and G. J. Balas, “Comparing linear parameter-varying gain-scheduled control techniques for active flutter suppression,” Journal of Guidance, Control, and Dynamics, vol. 23, no. 5, pp. 948–955, 2000\.
  * [3] T. He and W. Su, “Gust alleviation of highly flexible aircraft with model predictive control,” in AIAA SCITECH 2023 Forum, p. 0586, 2023.
  * [4] T. He, G. Zhu, S. S.-M. Swei, and W. Su, “Smooth-switching lpv control for vibration suppression of a flexible airplane wing,” Aerospace Science and Technology, vol. 84, pp. 895–903, 2019.
  * [5] N. Fonzi, S. L. Brunton, and U. Fasel, “Data-driven nonlinear aeroelastic models of morphing wings for control,” Proceedings of the Royal Society A, vol. 476, no. 2239, p. 20200079, 2020.
  * [6] A. K. Al-Jiboory, G. Zhu, S. S.-M. Swei, W. Su, and N. T. Nguyen, “Lpv modeling of a flexible wing aircraft using modal alignment and adaptive gridding methods,” Aerospace science and technology, vol. 66, pp. 92–102, 2017.
  * [7] K. Li, J. Kou, and W. Zhang, “Deep neural network for unsteady aerodynamic and aeroelastic modeling across multiple mach numbers,” Nonlinear Dynamics, vol. 96, pp. 2157–2177, 2019.
  * [8] S. Timme, S. Marques, and K. Badcock, “Transonic aeroelastic stability analysis using a kriging-based schur complement formulation,” AIAA journal, vol. 49, no. 6, pp. 1202–1213, 2011.
  * [9] J.-N. Juang and R. S. Pappa, “An eigensystem realization algorithm for modal parameter identification and model reduction,” Journal of guidance, control, and dynamics, vol. 8, no. 5, pp. 620–627, 1985.
  * [10] W. Zhang and Z. Ye, “Reduced-order-model-based flutter analysis at high angle of attack,” Journal of Aircraft, vol. 44, no. 6, pp. 2086–2089, 2007.
  * [11] W. Zhang, Z. Ye, and C. Zhang, “Aeroservoelastic analysis for transonic missile based on computational fluid dynamics,” Journal of Aircraft, vol. 46, no. 6, pp. 2178–2183, 2009.
  * [12] J. L. Lumley, “The structure of inhomogeneous turbulent flows,” Atmospheric turbulence and radio wave propagation, pp. 166–178, 1967.
  * [13] G. Berkooz, P. Holmes, and J. L. Lumley, “The proper orthogonal decomposition in the analysis of turbulent flows,” Annual review of fluid mechanics, vol. 25, no. 1, pp. 539–575, 1993.
  * [14] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the national academy of sciences, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [15] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” Proceedings of the Royal Society A, vol. 474, no. 2219, p. 20180335, 2018.
  * [16] A. Iannelli, U. Fasel, and R. S. Smith, “The balanced mode decomposition algorithm for data-driven lpv low-order models of aeroservoelastic systems,” Aerospace Science and Technology, vol. 115, p. 106821, 2021.
  * [17] P. Benner, S. Gugercin, and K. Willcox, “A survey of projection-based model reduction methods for parametric dynamical systems,” SIAM review, vol. 57, no. 4, pp. 483–531, 2015.
  * [18] R. Palacios, J. Murua, and R. Cook, “Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft,” AIAA Journal, vol. 48, no. 11, pp. 2648–2659, 2010.
  * [19] D. A. Peters and M. J. Johnson, “Finite-state airloads for deformable airfoils on fixed and rotating wings,” in Symposium on Aeroelasticity and Fluid/Structure Interaction (P. P. Friedmann and J. C. I. Chang, eds.), vol. AD-44 of ASME Winter Annual Meeting, (New York, NY), pp. 1–28, American Society of Mechanical Engineers, 1994.
  * [20] D. A. Peters, S. Karunamoorthy, and W.-M. Cao, “Finite state induced flow models part i: Two-dimensional thin airfoil,” Journal of Aircraft, vol. 32, no. 2, pp. 313–322, 1995.
  * [21] D. A. Peters and C.-J. He, “Finite state induced flow models part ii: Three-dimensional rotor disk,” Journal of Aircraft, vol. 32, no. 2, pp. 323–333, 1995.
  * [22] C. M. Shearer and C. E. S. Cesnik, “Nonlinear flight dynamics of very flexible aircraft,” Journal of Aircraft, vol. 44, no. 5, pp. 1528–1545, 2007.
  * [23] W. Su and C. E. S. Cesnik, “Nonlinear aeroelasticity of a very flexible blended-wing-body aircraft,” Journal of Aircraft, vol. 47, no. 5, pp. 1539–1553, 2010.
  * [24] W. Su and C. E. S. Cesnik, “Dynamic response of highly flexible flying wings,” AIAA Journal, vol. 49, no. 2, pp. 324–339, 2011.
  * [25] M. J. Patil, D. H. Hodges, and C. E. S. Cesnik, “Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft,” Journal of Aircraft, vol. 38, no. 1, pp. 88–94, 2001.