# AMORE-UPF at SemEval-2018 Task 4: BiLSTM with Entity Library

Laura Aina Carina Silberer11footnotemark: 1 Ionut-Teodor
Sorodoc11footnotemark: 1
Matthijs Westera11footnotemark: 1 Gemma Boleda
Universitat Pompeu Fabra
Barcelona, Spain
<EMAIL_ADDRESS>denotes equal contribution.

###### Abstract

This paper describes our winning contribution to SemEval 2018 Task 4:
_Character Identification on Multiparty Dialogues_. It is a simple, standard
model with one key innovation, an _entity library_. Our results show that this
innovation greatly facilitates the identification of infrequent characters.
Because of the generic nature of our model, this finding is potentially
relevant to any task that requires effective learning from sparse or
unbalanced data.

##  1 Introduction

SemEval 2018 Task 4 is an entity linking task on multiparty
dialogue.111https://competitions.codalab.org/competitions/17310 It consists in
predicting the referents of nominals that refer to a person, such as she, mom,
Judy – henceforth _mentions_. The set of possible referents is given
beforehand, as well as the set of mentions to resolve. The dataset used in
this task is based on Chen and Choi (2016) and Chen et al. (2017), and
consists of dialogue from the TV show _Friends_ in textual form.

Our main interest is whether deep learning models for tasks like entity
linking can benefit from having an explicit entity library, i.e., a component
of the neural network that stores entity representations learned during
training. To that end, we add such a component to an otherwise relatively
basic model – a bidirectional LSTM (long short-term memory; Hochreiter and
Schmidhuber 1997), the standard neural network model for sequential data like
language. Training and evaluating this model on the task shows that the entity
library is beneficial in the case of infrequent entities.222Source code for
our model and for the training procedure is published on
https://github.com/amore-upf/semeval2018-task4.

##  2 Related Work

Previous entity linking tasks concentrate on linking mentions to Wikipedia
pages (Bunescu and Paşca 2006; Mihalcea and Csomai 2007 and much subsequent
work; for a recent approach see Francis-Landau et al. 2016). By contrast, in
the present task (based on Chen and Choi 2016; Chen et al. 2017) only a list
of entities is given, without any associated encyclopedic entries. This makes
the task more similar to the way in which a human audience might watch the TV
show, in that they are initially unfamiliar with the characters. What also
sets the present task apart from most previous tasks is its focus on
multiparty dialogue (as opposed to, typically, newswire articles).

A task that is closely related to entity linking is _coreference resolution_ ,
i.e., the task of clustering mentions that refer to the same entity (e.g., the
CoNLL shared task of Pradhan et al. 2011). Since mention clusters essentially
correspond to entities (an insight central to the approaches to coreference in
Haghighi and Klein 2010; Clark and Manning 2016), the present task can be
regarded as a type of coreference resolution, but one where the set of
referents to choose from is given beforehand.

Since our main aim is to test the benefits of having an entity library, in
other respects our model is kept more basic than existing work both on entity
linking and on coreference resolution (e.g., the aforementioned approaches, as
well as Wiseman et al. 2016; Lee et al. 2017, Francis-Landau et al. 2016). For
instance, we avoid feature engineering, focusing instead on the model’s
ability to learn meaningful entity representations from the dialogue itself.
Moreover, we deviate from the common strategy to entity linking of
incorporating a specialized coreference resolution module (e.g., Chen et al.
2017).

##  3 Model Description

We approach the task of character identification as one of multi-class
classification. Our model is depicted in Figure 1, with inputs in the top left
and outputs at the bottom. In a nutshell, our model is a bidirectional LSTM
(long short-term memory, Hochreiter and Schmidhuber 1997) that processes the
dialogue text and resolves mentions, through a comparison between the LSTM’s
hidden state (for each mention) to vectors in a learned entity library.

The model is given chunks of dialogue, which it processes token by token. The
$i$th token $\mathbf{t}_{i}$ and its speakers $S_{i}$ (typically a singleton
set) are represented as one-hot vectors, embedded via two distinct embedding
matrices ($\mathbf{W}_{t}$ and $\mathbf{W}_{s}$, respectively) and finally
concatenated to form a vector $\mathbf{x}_{i}$ (Eq. 1; see also Figure 1). In
case $\mathbf{S}_{i}$ contains multiple speakers, their embeddings are summed.

$\mathbf{x}_{i}=\mathbf{W}_{t}\ \mathbf{t}_{i}\
\|\displaystyle\sum_{\mathbf{s}\in S_{i}}\mathbf{W}_{s}\ \mathbf{s}$ (1)

We apply an activation function $f$ $(=\text{tanh})$. The hidden state
$\overrightarrow{\mathbf{h}}_{\\!i}$ of a _uni_ directional LSTM for the $i$th
input is recursively defined as a combination of that input with the LSTM’s
previous hidden state $\overrightarrow{\mathbf{h}}_{\\!i-1}$. For a _bi_
directional LSTM, the hidden state $\mathbf{h}_{i}$ is a concatenation of the
hidden states $\overrightarrow{\mathbf{h}}_{\\!i}$ and
$\overleftarrow{\mathbf{h}}_{\\!i}$ of two unidirectional LSTMs which process
the data in opposite directions (Eq. 2; see also Figure 1). In principle, this
enables a bidirectional LSTM to represent the entire dialogue with a focus on
the current input, including for instance its relevant dependencies on the
context.

$\mathbf{h}_{i}=\text{BiLSTM}(f(\mathbf{x}_{i}),\overrightarrow{\mathbf{h}}_{\\!i-1},\overleftarrow{\mathbf{h}}_{\\!i+1})$
(2)

In the model, learned representations of each entity are stored in the entity
library $\mathbf{E}\in\mathbb{R}^{N\times k}$ (see Figure 1): $\mathbf{E}$ is
a matrix which represents each of $N$ entities through a $k$-dimensional
vector, and whose values are updated (only) during training. For every token
$\mathbf{t}_{i}$ that is tagged as a mention,333For multi-word mentions this
is done only for the last token in the mention. we map the corresponding
hidden state $\mathbf{h}_{i}$ to a vector
$\mathbf{e}_{i}\in\mathbb{R}^{1\times k}$. This extracted representation is
used to retrieve the (candidate) referent of the mention from the entity
library: The similarity of $\mathbf{e}_{i}$ to each entity representation
stored in $\mathbf{E}$ is computed using cosine, and softmax is then applied
to the resulting similarity profile to obtain a probability distribution
$\mathbf{o}_{i}\in[0,1]^{1\times N}$ over entities (‘class scores’ in Figure
1):

$\mathbf{o}_{i}=\text{softmax}(\text{cosine}(\mathbf{E},\underbrace{(\mathbf{W}_{o}\
\mathbf{h}_{i}+\mathbf{b})}_{\mathbf{e}_{i}})$ (3)

Figure 1: The AMORE-UPF model (bias not depicted).

At testing time, the model’s prediction $\hat{c}_{i}$ for the $i$th token is
the entity with highest probability:

$\hat{c}_{i}=\text{argmax}(\mathbf{o}_{i})$ (4)

We train the model with backpropagation, using negative log-likelihood as loss
function. Besides the BiLSTM parameters, we optimize $\mathbf{W}_{t}$,
$\mathbf{W}_{s}$, $\mathbf{W}_{o}$, $\mathbf{E}$ and $\mathbf{b}$. We refer to
this model as AMORE-UPF, our team name in the SemEval competition. Note that,
in order for this architecture to be successful, $\mathbf{e}_{i}$ needs to be
as similar as possible to the entity vector of the entity to which mention
$t_{i}$ refers. Indeed, the mapping $\mathbf{W}_{o}$ should effectively
specialize in “extracting” entity representations from the hidden state
because of the way its output $\mathbf{e}_{i}$ is used in the model—to do
entity retrieval. Our entity retrieval mechanism is inspired by the attention
mechanism of Bahdanau et al. (2016), that has been used in previous work to
interact with an external memory Sukhbaatar et al. (2015); Boleda et al.
(2017). To assess the contribution of the entity library, we compare our model
to a similar architecture which does not include it (NoEntLib). This model
directly applies softmax to a linear mapping of the hidden state (Eq. 5,
replacing Eq. 3 above).

$\mathbf{o}_{i}=\text{softmax}(\mathbf{W}_{o}\ \mathbf{h}_{i}+\mathbf{b})$ (5)

##  4 Experimental Setup

### Data

We use the training and test data provided for SemEval 2018 Task 4, which span
the first two seasons of the TV show Friends, divided into scenes (train:
$374$ scenes from $47$ episodes; test: $74$ scenes from $40$ episodes).444The
organizers also provided data divided by episodes rather than scenes, which we
didn’t use. In total, the training and test data contain 13,280 and 2,429
nominal mentions (e.g., Ross, I; Figure 2), respectively, which are annotated
with the ID of the entity to which they refer (e.g., 335, 183). The utterances
are further annotated with the name of the speaker (e.g., Joey Tribbiani).
Overall there are 372 entities in the training data (test data: 106).

Joey Tribbiani (183): ”…see Ross, because I think you love her .” 335 183 335
306

Figure 2: Example of the data provided for the SemEval 2018 Task 4. It shows
the speaker (first line) of the utterance (second line) and the ids of the
entities to which the target mentions (underlined) refer (last line).

Our models do not use any of the provided automatic linguistic annotations,
such as PoS or named entity tags.

We additionally used the publicly available 300-dimensional word vectors that
were pre-trained on a Google News corpus with the word2vec Skip-gram model
(Mikolov et al., 2013).555The word vectors are available at
https://code.google.com/archive/p/word2vec/.

### Parameter settings

Using 5-fold cross-validation on the training data, we performed a random
search over the hyperparameters and chose those which yielded the best mean
F1-score. Specifically, our submitted model is trained in batch mode using the
Adam optimizer Kingma and Ba (2014) with a learning rate of $0.0005$. Each
batch covers $24$ scenes, which are given to the model in chunks of $757$
tokens. The token embeddings ($\mathbf{W}_{t}$) are initialized with the
word2vec vectors. Dropout rates of $0.008$ and $0.0013$ are applied on the
input $\mathbf{x}_{i}$ and hidden layer $\mathbf{h}_{i}$ of the LSTM,
respectively. The size of $\mathbf{h}_{i}$ is set to $459$ units, the
embeddings of the entity library $\mathbf{E}$ and speakers $\mathbf{W}_{s}$
are set to $k~{}\mbox{= 134}$ dimensions.

Other configurations, including randomly initialized token embeddings, weight
sharing between $\mathbf{E}$ and $\mathbf{W}_{s}$, self-attention Bahdanau et
al. (2016) on the input layer, a uni-directional LSTM, and rectifier or linear
activation function $f$ on the input embeddings did not improve performance.

For the final submission of the answers for the test data, we created an
ensemble model by averaging the output (Eq. 3) of the five models trained on
the different folds.

##  5 Results

| all entities | main entities
---|---|---
Models | F1 | Acc | F1 | Acc
AMORE-UPF | $\mathbf{41.1^{**}}$ | $\mathbf{74.7^{**}}$ | $79.4$ | $77.2$
NoEntLib | $26.4$ | $71.6$ | $\mathbf{79.5}$ | $\mathbf{77.5}$
Table 1: Results obtained for the submitted AMORE-UPF model and a variant of
it that does not use an entity library (NoEntLib). Best results are in
boldface. Differences with respect to the 2nd row marked by ‘**’ are
significant at the $0.001$ probability level (see text).

Two evaluation conditions were defined by the organizers – _all entities_ and
_main entities_ – with macro-average F1-score and label accuracy as the
official metrics, and macro-average F1-score in the _all entities_ condition
applied to the leaderboard. The _all entities_ evaluation has 67 classes: 66
for entities that are mentioned at least 3 times in the test set and one
grouping all others. The _main entities_ evaluation has 7 classes, 6 for the
main characters and one for all the others. Among all four participating
systems in this SemEval task our model achieved the highest score on the _all
entities_ evaluation, and second-highest on the _main entities_ evaluation.

Figure 3: Distribution of all $2,429$ target mentions in the test data in
terms of their part-of-speech.

Table 1 gives our results in the two evaluations, comparing the models
described in Section 4. While both models perform on a par on main entities,
AMORE-UPF outperforms NoEntLib by a substantial margin when all characters are
to be predicted (+15 points in F1-score, +3 points in accuracy; Table
1).666The mean difference between the single models (trained on a single fold)
and the ensemble AMORE-UPF is between -1.3 points (accuracy main entities, std
$=1.3$) and -4.2 points (F1-score all entities, std $=1.3$). The difference
between the models with/without an entity library are statistically
significant based on approximate randomization tests Noreen (1989), with the
significance level $p<0.001$. This shows that the use of an entity library can
be beneficial for the linking of rarely mentioned characters.

Figure 3 shows that most of the target mentions in the test data fall into one
of five grammatical categories. The dataset contains mostly pronouns (83%),
with a very high percentage of first person pronouns (44%). Figures 4 and 5
present the accuracy and F1-score which the two models described above obtain
on all entities for different categories of mentions. The entity library is
beneficial when the mention is a first person pronoun or a proper noun (with
an increase of 30 points in F1-score for both categories; Figure 4), and
closer inspection revealed that this effect was larger for rare entities.

##  6 Discussion

The AMORE-UPF model consists of a bidirectional LSTM linked to an entity
library. Compared to an LSTM without entity library, NoEntLib, the AMORE-UPF
model performs particularly well on rare entities, which explains its top
score in the _all entities_ condition of SemEval 2018 Task 4. This finding is
encouraging, since rare entities are especially challenging for the usual
approaches in NLP, due to the scarcity of information about them.

Figure 4: F1-score of the models on all entities depending on the part-of-
speech of the target mentions. Figure 5: Aaccuracy of the models on all
entities depending on the part-of-speech of the target mentions.

We offer the following explanation for this beneficial effect of the entity
library, as a hypothesis for future work. Having an entity library requires
the LSTM of our model to output some representation of the mentioned entity,
as opposed to outputting class scores more or less directly as in the variant
NoEntLib. Outputting a meaningful entity representation is particularly easy
in the case of first person pronouns and nominal mentions (where the effect of
the entity library appears to reside; Figure 4): the LSTM can learn to simply
forward the speaker embedding unchanged in the case of pronoun I, and the
token embedding in the case of nominal mentions. This strategy does not
discriminate between frequent and rare entities; it works for both alike. We
leave further analyses required to test this potential explanation for future
work.

Future work may also reveal to what extent the induced entity representations
may be useful in others, to what extent they encode entities’ attributes and
relations (cf. Gupta et al. 2015), and to what extent a module like our entity
library can be employed elsewhere, in natural language processing and beyond.

## Acknowledgments

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 715154), and from the Spanish Ramón y Cajal programme
(grant RYC-2015-18907). We are grateful to the NVIDIA Corporation for the
donation of GPUs used for this research. We are also very grateful to the
Pytorch developers. This paper reflects the authors’ view only, and the EU is
not responsible for any use that may be made of the information it contains.

## References

  * Bahdanau et al. (2016) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. 2016.  Neural Machine Translation by Jointly Learning to Align and Translate.  _Computing Research Repository_ , abs/1409.0473v7.
  * Boleda et al. (2017) Gemma Boleda, Sebastian Padó, Nghia The Pham, and Marco Baroni. 2017.  Living a Discrete Life in a Continuous World: Reference in Cross-modal Entity Tracking.  In _IWCS 2017 — 12th International Conference on Computational Semantics — Short papers_.
  * Bunescu and Paşca (2006) Razvan Bunescu and Marius Paşca. 2006.  Using Encyclopedic Knowledge for Named Entity Disambiguation.  In _11th conference of the European Chapter of the Association for Computational Linguistics_.
  * Chen et al. (2017) Henry Y. Chen, Ethan Zhou, and Jinho D. Choi. 2017.  Robust Coreference Resolution and Entity Linking on Dialogues: Character Identification on TV Show Transcripts.  In _Proceedings of the 21st Conference on Computational Natural Language Learning (CoNLL 2017)_ , pages 216–225.
  * Chen and Choi (2016) Yu-Hsin Chen and Jinho D Choi. 2016.  Character Identification on Multiparty Conversation: Identifying Mentions of Characters in TV Shows.  In _Proceedings of the 17th Annual Meeting of the Special Interest Group on Discourse and Dialogue_ , pages 90–100.
  * Clark and Manning (2016) Kevin Clark and Christopher D. Manning. 2016.  Improving Coreference Resolution by Learning Entity-Level Distributed Representations.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 643–653.
  * Francis-Landau et al. (2016) Matthew Francis-Landau, Greg Durrett, and Dan Klein. 2016.  Capturing Semantic Similarity for Entity Linking with Convolutional Neural Networks.  _Computing Research Repository_ , abs/1604.00734.
  * Gupta et al. (2015) Abhijeet Gupta, Gemma Boleda, Marco Baroni, and Sebastian Padó. 2015.  Distributional Vectors Encode Referential Attributes.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 12–21.
  * Haghighi and Klein (2010) Aria Haghighi and Dan Klein. 2010.  Coreference Resolution in a Modular, Entity-Centered Model.  In _Human Language Technologies: The 2010 Annual Conference of the North American Chapter of the Association for Computational Linguistics_ , pages 385–393. Association for Computational Linguistics.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780.
  * Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba. 2014.  Adam: A Method for Stochastic Optimization.  _Computing Research Repository_ , abs/1412.6980.
  * Lee et al. (2017) Kenton Lee, Luheng He, Mike Lewis, and Luke Zettlemoyer. 2017.  End-to-End Neural Coreference Resolution.  _Computing Research Repository_ , abs/1707.07045.
  * Mihalcea and Csomai (2007) Rada Mihalcea and Andras Csomai. 2007.  Wikify!: Linking Documents to Encyclopedic Knowledge.  In _Proceedings of the sixteenth ACM Conference on Information and Knowledge Management_ , pages 233–242.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Distributed Representations of Words and Phrases and Their Compositionality.  In _Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 2_ , pages 3111–3119.
  * Noreen (1989) E.W. Noreen. 1989.  _Computer-Intensive Methods for Testing Hypotheses: An Introduction_.  Wiley.
  * Pradhan et al. (2011) Sameer Pradhan, Lance Ramshaw, Mitchell Marcus, Martha Palmer, Ralph Weischedel, and Nianwen Xue. 2011.  CoNLL-2011 Shared Task: Modeling Unrestricted Coreference in OntoNotes.  In _Proceedings of the Fifteenth Conference on Computational Natural Language Learning: Shared Task_ , pages 1–27. Association for Computational Linguistics.
  * Sukhbaatar et al. (2015) Sainbayar Sukhbaatar, Arthur Szlam, Jason Weston, and Rob Fergus. 2015.  End-To-End Memory Networks.  _arXiv_ , pages 1–11.
  * Wiseman et al. (2016) Sam Wiseman, Alexander M Rush, and Stuart M Shieber. 2016.  Learning Global Features for Coreference Resolution.  _Computing Research Repository_ , abs/1604.03035.


# An Integrated Mechanical Intelligence and Control Approach Towards Flight
Control of Aerobat

Eric Sihite, Atefe Darabi, Pravin Dangol, Andrew Lessieur, and Alireza
Ramezani1 1SiliconSynapse Laboratory, ECE Department, Northeastern University,
Boston, MA, USA. emails: {e.sihite, darabi.a, dangol.p, lessieur.a,
a.ramezani} @northeastern.edu

###### Abstract

Our goal in this work is to expand the theory and practice of robot locomotion
by addressing critical challenges associated with the robotic biomimicry of
bat aerial locomotion. Bats are known for their pronounced, fast wing
articulations, e.g., bats can mobilize as many as forty joints during a single
wingbeat, with some joints reaching to over one thousand degrees per second in
angular speed. Copying bats flight is a significant ordeal, however, very
rewarding. Aerial drones with morphing bodies similar to bats can be safer,
agile and energy efficient owing to their articulated and soft wings. Current
design paradigms have failed to copy bat flight because they assume only
closed-loop feedback roles and ignore computational roles carried out by
morphology. To respond to the urgency, a design framework called Morphing via
Integrated Mechanical Intelligence and Control (MIMIC) is proposed. In this
paper, using the dynamic model of Northeastern University’s Aerobat, which is
designed to test the effectiveness of the MIMIC framework, it will be shown
that computational structures and closed-loop feedback can be successfully
used to mimic bats stable flight apparatus.

##  I Introduction

In recent years, there has been increased focus on making our residential
spaces smarter, safer, more efficient and closer to the materialization of the
concept of smart cities [1]. As a result, safety and security robots are
gaining ever growing importance [2] and drive a lucrative market. Smart cities
market was valued at USD 624.81 billion in 2019 and is expected to reach USD
1712.83 billion by 2025 [3].

Among these security robots, ground robots (wheeled and legged) are used the
most despite known limitations such as not being able to reach a vantage point
for surveillance, having limited operation time, having the potential to
collide with and harm humans in crowded spaces (e.g., sidewalks, airports,
etc.), not being able to negotiate rough terrain (e.g., street curbs, stairs,
bumps, etc.) or large obstacles (e.g., a wall or bush), and possessing slow
mobility [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Despite their superior
mobility, which makes them extremely suitable for civic surveillance and
monitoring applications, the contribution of Micro Aerial Vehicles (MAVs) to
this market has remained very limited. Current state-of-the-art MAVs with fast
rotating propellers and rigid structures pose extreme dangers to humans, e.g.,
they can cause penetrating injury, laceration resulting in blood loss and
massive destruction of the human body [15]. Furthermore, they cannot survive
unavoidable crashes in unstructured environments of cities or operate for more
than thirty to forty minutes. That said, the application of safety features in
these systems have not solved the problems and their operations in residential
spaces have remained limited by strict rules from the Federal Aviation
Administration (FAA) [15]. To be able to use MAVs in civic applications, there
is a need to transform their safety and efficiency.

Figure 1: Illustration of Northeastern University’s Aerobat.

The research goal in this work is to expand the theory and practice of robot
locomotion by addressing critical challenges associated with the robotic
biomimicry of bat aerial locomotion. The resulting MAVs with morphing bodies
will be safe, agile and energy efficient owing to their articulated, soft
wings and will autonomously operate with long operation lifespans. Bat
membranous wings possess unique functions [16] that make them a good example
to take inspiration from and transform safety, agility and efficiency of
current aerial drones. In contrast with other flying vertebrates, bats have an
extremely articulated musculoskeletal system (Fig. 2-B) which is key to their
body impact survivability and their impressively adaptive and multimodal
locomotion behavior [17]. Bats exclusively use this capability with their
structural flexibility to generate the controlled force distribution on each
membrane wing. Wing flexibility, complex wing kinematics, and fast muscle
actuation allow these creatures to change their body configuration within a
few tens of milliseconds. These characteristics are crucial to their unrivaled
agility and energetic efficiency [18].

Complex locomotion styles achieved through such synchronous movements of many
joints are showcased by several other species. There is an urgency for new
paradigms providing insight into how these animals mobilize and regulate so
many joints with small brains (limited computation power) and small muscles
(limited actuation power). That said, widely used paradigms have failed to
copy bat flight [19, 20] because they assume only closed-loop feedback roles
and ignore computational roles carried out by morphology. To respond to the
urgency, a design framework called Morphing via Integrated Mechanical
Intelligence and Control (MIMIC) is proposed in this work. In this work, using
simulation results, it will be shown that the MIMIC framework can be
successfully used to mimic bats flight apparatus known for their pronounced,
fast wing articulations, e.g., bats can mobilize as many as forty joints
during a single wingbeat, with some joints reaching to over one thousand
degrees per second in angular speed.

We expand upon our previous work with the bat robots [21, 22, 23, 24, 25]
where we will apply the MIMIC framework to our most recent morphing wing
design in [26]. This design captures the elbow flexion and extension that
allow the wing to fold during upstroke. The MIMIC framework will be the
continuation of our past attempts at developing control framework for bio-
inspired flapping wing drones in [27, 23, 28, 29, 30, 31, 32].

Figure 2: The kinetic sculpture (KS) design and the motivation to mimic a
bat’s natural flapping gait. (A) Illustrates the KS made of monolithically
fabricated rigid and flexible materials. (B) Depicts bat flapping gait with
the wing folding/expansion within one flapping cycle. (C) Illustrates
simulated Aerobat wingbeat cycle. (D) Shows the sensitivity analysis resutls.
Feedback-Driven Components (FDCs) used in the KS can change the behavior of
the structure which can be leveraged for flight design purpose.

This paper is outlined as follows: a brief overview of our bat robot and
testing platform, the Aerobat as shown in Fig. 1, followed by the dynamic
modeling, control, optimization, simulation, and concluding remarks.

##  II Overview of Aerobat, a Platform to Test the MIMIC Framework

Our robotic bat arm-wing, shown in Fig. 2-A and Fig. 3, is a mechanical
structure with computational roles designed to mimic the flying maneuvers of
biological bats [26]. We also refer to this structure with the term: Kinetic
Sculpture. The computational roles of this structure are actively modulated
using closed-loop feedback as this will be discussed later in this paper. This
arm-wing is composed of rigid links and flexible joints monolithically
fabricated by using PolyJet 3D printing technology. The mechanism for driving
the arm-wing which produces the flapping motion is composed of several gears,
cranks, and four-bar mechanisms which are actuated by a single motor. The
resulting flapping motion follows some of the biologically meaningful degrees-
of-freedom (DoF) present in the natural bat’s flapping flight, which in our
mechanism are the wing plunging motion and elbow extension/retraction, as
shown in Fig. 2-C. The arm-wing extends during the downstroke and retracts
during the upstroke which serves to minimize the negative lift and increase
the flapping gait efficiency. The dynamic modeling and simulation of Aerobat
flapping under four meaningful DoFs (plunging, elbow flexion/extension,
mediolateral movement, and feathering) was investigated in [32] where the wing
is modeled as a rigid plate and using optimization to find a steady flapping
gait and upside-down perching maneuver.

A more sophisticated structure is proposed to add control and morphological
freedom into this design. We propose to use Feedback-Driven Components (FDC)
to adjust the length of some links as shown in Fig. 3. This serves as a way to
change the arm-wing morphology resulting in a change in the end effector
trajectories, as shown in Fig. 2-D. The trajectories of interest are the
humerus and radius links end effectors, which are shown as $L_{5}$ and
$L_{12}$ in Fig. 3 respectively. The joint angles of these two links represent
the shoulder and elbow angles, which correspond to the plunging motion and the
elbow extension respectively. This framework is relevant to the concept of
morphological computation, where a simple control action can influence and
achieve complex manipulation in the body morphology which is commonly seen in
nature [33].

In order to simplify the system and facilitate control, the mechanical
structure is modeled as a network of massed and massless links connected by
compliant joints. Such mechanical systems that include elements with
comparatively small or zero mass and inertial distribution are considered as
singular mechanical systems. In this system, the massless network is used as a
guide for the massed system through linear/torsional spring and damper to
simulate the flexible joints of the kinetic sculpture, as shown in Fig. 4.
Here, only the humerus and radius links, in addition to the robot’s body, are
modeled as a massed system which significantly simplifies the dynamic modeling
and simulation.

Figure 3: Shows the wing structure which is composed of a network of
mechanical linkages. The system consists of 12 links
($\\{L_{1},\dots,L_{12}\\}$) and 17 joints/hinges ($\\{j_{1},...,j_{17}\\}$).
Joints 1 and 9 are the gears’ center of rotation where they will be driven up
to 10 Hz in angular speed. Four FDCs adjust the length of the linkage to
change the resulting end-effector trajectories and facilitate control (see the
sensitivity analysis results shown above.

##  III Dynamic Modeling of Mechanical Computational Structures in Aerobat

This section outlines the derivations of the equations of motion for the
massless and massed subsystems as described in Section II. There are two
primary frames of reference considered in this paper: the inertial and body
frames. The body frame is defined such that the $x$, $y$, and $z$ axes point
towards the robot’s front, left and top sides, respectively. The rotation from
body frame to inertial frame is defined as $\bm{x}=R_{B}\,\bm{x}^{B}$, where
the superscript $B$ represents the vectors defined in body frame.

The wing structure is modeled as a planar linkage mechanism and the rotation
matrix $R(\theta)$ is used to rotate from link’s local frame to the body
frame, where $\theta$ is the absolute angle of the link with respect to the
$y$-axis of body frame. The wing conformation parameters names and values used
to derive the equations in this paper are the same as the wing design in [26],
where the length $l_{i}$ is used to describe the length parameters of link $i$
($L_{i}$) shown in Fig. 3. For links that are not simply straight, the
dimensions involved are represented with additional alphabet subscripts (e.g.
$l_{3a}$, $l_{3b}$, etc).

###  III-A Massless Subsystem Kinematics Formulation

The equations of motion of the massless subsystem can be modeled from its
kinematics. This wing flaps in a planar motion, so for the sake of simplicity,
the kinematic equations are derived in 2D coordinates in the body’s $y$-$z$
plane. The superscript $B$ for the body frame representation is omitted to
simplify the derivations. Also, because of the symmetry of the left and right
wings, their kinematics follow the same formulation, therefore a general form
of equations are presented in this section.

As shown in Fig. 3, the massless links includes three closed kinematic chains.
Additionally, the crank gears are coupled to each other to have the same
angular velocity. Let $\bm{p}_{i}$ and $\theta_{i}$ represent the positions
and absolute angles of joint $i\in\\{1,\dots,17\\}$ about the body frame as
labeled in Fig. 3. The massless linkage kinematics can be derived from the
closed kinematic loop of the four-bar linkages. Derive the positions of the
joints 3, 11, and 15 as a function of $\theta_{i}$ following the kinematic
chain from the two closest stationary joints (i.e. joints 1, 4, 9, 12, and
14), represented using the subscripts $A$ and $B$, (e.g. $\bm{p}_{3A}$). For
example, the position of joint 3, $\bm{p}_{3}$, can be derived as follows

$\begin{gathered}\bm{p}_{3A}=\bm{p}_{1}+R(\theta_{1})[l_{1},0]^{\top}+R(\theta_{2})[l_{2},0]^{\top}\\\
\bm{p}_{3B}=\bm{p}_{4}+R(\theta_{4})[-l_{3a}-l_{3b},0]^{\top}\end{gathered}$
(1)

where $R(\theta_{i})$ is a 2D rotation matrix and $l_{j}$ are the length
variables representing the conformation of link $j$. Then the linkage
kinematic constraints can be formed as

$\bm{q}_{1}=[\theta_{1},\theta_{2},\theta_{4},\theta_{9},\theta_{10},\theta_{12},\theta_{13},\theta_{14},l_{3b},l_{3c},l_{8b},l_{10b}]^{\top}$
$\begin{gathered}\bm{C}(\bm{q}_{1})=\begin{bmatrix}\bm{p}_{3A}-\bm{p}_{3B}\\\
\bm{p}_{11A}-\bm{p}_{11B}\\\ \bm{p}_{15A}-\bm{p}_{15B}\\\
\theta_{1}-\theta_{9}-\Delta\phi\end{bmatrix}=\bm{0}_{7\times
1},\end{gathered}$ (2)

where $\bm{q}_{1}$ is formed from the non-constant variables and $\Delta\phi$
is the phase difference between the two crank gears ($\theta_{1}$ and
$\theta_{9}$). The linkage dimensions are constant except for $l_{3b}$,
$l_{3c}$, $l_{8b}$, and $l_{10b}$ which represent the length of the FDCs.

The equation of motion for the massless subsystem can be derived by taking the
second time derivative of $\bm{C}(\bm{q}_{1})$, which can be rearranged into
the following form

$\bm{\ddot{C}}=M_{A}(\bm{q}_{1})\ddot{\bm{q}}_{1}+\bm{h}_{A}(\bm{q}_{1},\dot{\bm{q}}_{1})=\bm{0}_{7\times
1}.$ (3)

The input to this subsystem is the acceleration of the crank gear driven by
the motor ($\theta_{1}$) and the FDCs. These accelerations can be formed,
using abuse of notation, as follows

$\displaystyle M_{B}$ $\displaystyle=\begin{bmatrix}1&\bm{0}_{1\times 11}\\\
\bm{0}_{4\times 7}&I_{4\times 4}\end{bmatrix},$ $\displaystyle
M_{B}\ddot{\bm{q}}_{1}$ $\displaystyle=\bm{u}_{1},$ (4)

where $\bm{u}_{1}=[u_{g},u_{3b},u_{3c},u_{8b},u_{10b}]^{\top}$ represents
accelerations, $u_{g}$ is the acceleration of the crank gear $\theta_{1}$.
Finally, combining (3) and (4) results in the following equation of motion

$\begin{gathered}M_{1}\bm{\ddot{q}}_{1}+\bm{h}_{1}=B_{1}\bm{u}_{1}\\\
M_{1}=\begin{bmatrix}M_{A}\\\
M_{B}\end{bmatrix},\,\,\bm{h}_{1}=\begin{bmatrix}\bm{h}_{A}\\\ \bm{0}_{5\times
1}\end{bmatrix},\,\,B_{1}=\begin{bmatrix}\bm{0}_{7\times 5}\\\ I_{5\times
5}\end{bmatrix}\end{gathered}$ (5)

###  III-B Massed Subsystem Dynamics Formulation

Figure 4: Shows the wing massed subsystem which is composed of the humerus and
radius links. The external forces and torques imposed by the massless
kinematic chain and flexible hinges act on this subsystem.

The equation of motion of the massed subsystem can be derived using the Euler-
Lagrangian dynamic formulation. Let $\bm{x}_{i}$ represent the Center of Mass
(CoM) position of the massed subsystem,
$i\in\mathcal{L_{M}}=\\{B,H_{L},H_{R},R_{L},R_{R}\\}$, where $B$, $H$, and $R$
represent the body, humerus and radius while the subscripts $L$ and $R$
represent the left and right wing components, respectively. Similar to Section
III-A, only the general form of the formulations are presented by omitting the
$L$ and $R$ subscripts due to the wing symmetry. The positions of the humerus
and radius links are defined as follows

$\displaystyle\bm{x}_{H}$
$\displaystyle=\bm{x}_{B}+R_{B}\,(\bm{p}_{7}^{B}+R_{x}(\theta_{7})\,\bm{l}_{H}/2)$
(6) $\displaystyle\bm{x}_{R}$
$\displaystyle=\bm{x}_{H}+R_{B}\,(R_{x}(\theta_{7})\,\bm{l}_{H}/2+R_{x}(\theta_{8})\,\bm{l}_{R}/2),$

where $R_{x}$ is the Euler rotation about $x$ axis, $\bm{x}_{B}$ is the body
inertial position,
$\bm{l}_{H}=[0,l_{h}\cos(\alpha),l_{5b}+l_{h}\sin(\alpha)]^{\top}$, and
$\bm{l}_{R}=[0,l_{r},0]$. Additionally, $\theta_{7}=\theta_{s}-\alpha$ and
$\theta_{8}=\theta_{e}+\theta_{s}+\alpha$, in which $\theta_{s}$ and
$\theta_{e}$ are the shoulder and elbow angles of the left wing respectively.
$\alpha$, $l_{h}$, $l_{r}$, and $l_{5b}$ are the constant physical parameters
of the wing structure. Here, $l_{h}$ and $l_{r}$ are the length of the humerus
and radius links respectively while $\alpha$ corresponds to the initial angle
of the shoulder joint at $\theta_{7}=0$.

The generalized coordinates vector is composed of both linear and rotational
states, therefore both Hamiltonian and Euler-Lagrangian principles are
required to derive the equations of motion of the massed links. In this case,
in order to avoid gimbal lock, the body angular velocity is derived by using
the modified Euler-Lagrangian equation of motion in SO(3). The kinetic and
potential energies of the massed parts and their corresponding Lagrangian,
$L=T-U$, can be derived as follows

$\displaystyle T$
$\displaystyle=\textstyle{\frac{1}{2}\sum_{i\in\mathcal{L_{M}}}m_{i}{\dot{\bm{x}}_{i}}^{\top}{\dot{\bm{x}}_{i}}+\bm{\omega}_{i}^{\top}I_{i}\bm{\omega}_{i}}$
(7) $\displaystyle U$
$\displaystyle=\textstyle\sum_{i\in\mathcal{L_{M}}}m_{i}[0,0,g]\bm{x}_{i}$

where $m_{i}$ is the mass, $I_{i}$ the inertia matrix, and $\bm{\omega}_{i}$
the angular velocity in body frame. The angular velocity of every massed link
could be computed as the sum of angular velocity of the link and body,
$\bm{\omega}_{i}=[\dot{\theta}_{i},0,0]^{\top}+\bm{\omega}_{B}$ for $i\neq B$.

The equation of motion of the body rotation in SO(3) is defined as follows

$\displaystyle
R_{B}^{\top}=[\bm{r}_{B1},\bm{r}_{B2},\bm{r}_{B3}],\qquad\tfrac{d}{dt}R_{B}=R_{B}\,[\bm{\omega}_{B}]_{\times}$
(8) $\displaystyle\textstyle{\frac{d}{dt}\left(\frac{\partial
L}{\partial\bm{\omega}_{B}}\right)+\bm{\omega}_{B}\times\frac{\partial
L}{\partial\bm{\omega}_{B}}+\sum_{j=1}^{3}\bm{r}_{Bj}\times\frac{\partial
L}{\partial\bm{r}_{Bj}}=\bm{u}_{2B}}$

where $[\,\cdot\,]_{\times}$ represents the skew symmetric operator and
$\bm{u}_{2B}$ is the generalized force component about the body rotation.
Finally, the Euler-Lagrange equations of motion are derived for the rest of
the generalized coordinates
$\bm{q}_{2}=[\theta_{s_{L}},\theta_{e_{L}},\theta_{s_{R}},\theta_{e_{R}},\bm{x}_{B}^{\top}]^{\top}$,
as follows

$\textstyle{\derivative{t}\frac{\partial
L}{\partial\dot{\bm{q}}_{2}}-{\frac{\partial
L}{\partial\bm{q}_{2}}}=\bm{u}_{2A}}.$ (9)

Then combining (8) and (9), the equation of motion of the massed system can be
formulated as

$\begin{gathered}M_{2}[\ddot{\bm{q}}_{2}^{\top},\dot{\bm{\omega}}_{B}^{\top}]^{\top}+\bm{h}_{2}=B_{i}\bm{u}_{i}+B_{s}\bm{u}_{s}+\bm{u}_{a}\\\
\dot{R}_{B}=R_{B}\,[\bm{\omega}_{B}]_{\times},\end{gathered}$ (10)

where $M_{2}$ is the matrix associated with the mass and inertia, and
$\bm{h}_{2}$ is the vector of gravity and coriolis forces. This system is
subjected to several force components, where $\bm{u}_{i}$ is the internal
spring and damping forces of the flexible joints (joints 7 and 8),
$\bm{u}_{s}$ is the actuation forces caused by the massless subsystem, and
$\bm{u}_{a}$ is the generalized aerodynamic force. Finally, combining (5) and
(10) forms the full dynamic model for the numerical simulation.

###  III-C Forces on the Massed System

Figure 5: Illustrates the application of aerodynamic strip theory. The overall
aerodynamic force is calculated both on spanwise and chordwise strip elements.
Then, they are integrated to obtain the resultant aerodynamic force.

The flexible joints 7 and 8 are modeled as the combination of a torsional
spring and damping which can be derived as
$\bm{u}_{i}=-(k_{i}(\bm{\theta}_{j}-\bm{\theta}_{j,0})+b_{i}\dot{\bm{\theta}_{j}})$,
where $k_{i}$ and $b_{i}$ are the stiffness and damping coefficients
respectively,
$\bm{\theta}_{j}=[\theta_{sL},\theta_{eL},\theta_{sR},\theta_{eR}]^{\top}$ is
the vector containing all massed joint angles, and $\bm{\theta}_{j,0}$ is the
nominal angle of the torsional springs. The $B_{i}=[I_{4\times 4},0_{4\times
10}]^{\top}$ matrix is setup to actuate these joint angles directly. The
flexible joints and linkages of the massless subsystem actuates the mass
system by the spring and damper model shown in Fig. 4. Let $\bm{u}_{s,6}$ and
$\bm{u}_{s,17}$ represent the forces acting on joints 6 and 17, defined as
follows

$\displaystyle\bm{u}_{s,6}$
$\displaystyle=(k_{g}\,(|\bm{p}_{5}-\bm{p}_{6}|-l_{4})+b_{g}(\dot{\bm{p}}_{5}-\dot{\bm{p}}_{6})^{\top}\bm{e}_{6})\bm{e}_{6}$
(11) $\displaystyle\bm{u}_{s,17}$
$\displaystyle=(k_{g}\,(|\bm{p}_{16}-\bm{p}_{17}|-l_{11})+b_{g}(\dot{\bm{p}}_{16}-\dot{\bm{p}}_{17})^{\top}\bm{e}_{17})\bm{e}_{17},$

where $\bm{e}_{6}=\frac{\bm{p}_{5}-\bm{p}_{6}}{|\bm{p}_{5}-\bm{p}_{6}|}$ and
$\bm{e}_{17}=\frac{\bm{p}_{16}-\bm{p}_{17}}{|\bm{p}_{16}-\bm{p}_{17}|}$ are
unit vectors, while $k_{g}$ and $b_{g}$ are spring and damping coefficients,
respectively. The generalized coordinates can be derived using the virtual
displacement and velocity formulation as follows

$B_{s,6}=\left(\partial\dot{\bm{p}}_{6}/\partial\dot{\bm{q}}_{2}^{\prime}\right)^{\top},\qquad
B_{s,17}=\left(\partial\dot{\bm{p}}_{17}/\partial\dot{\bm{q}}_{2}^{\prime}\right)^{\top}$
(12)

where
$\dot{\bm{q}}_{2}^{\prime}=[\dot{\bm{q}}_{2}^{\top},\bm{\omega}_{B}^{\top}]^{\top}$.
Then $B_{g}$ and $\bm{u}_{g}$ in (10) can be formulated as follows

$\displaystyle B_{g}$
$\displaystyle=[B_{g,6_{L}},B_{g,17_{L}},B_{g,6_{R}},B_{g,17_{R}}]$ (13)
$\displaystyle\bm{u}_{g}$
$\displaystyle=[\bm{u}_{g,6_{L}}^{\top},\bm{u}_{g,17_{L}}^{\top},\bm{u}_{g,6_{R}}^{\top},\bm{u}_{g,17_{R}}^{\top}]^{\top},$

where the subscripts $L$ and $R$ represent the left and right wing
respectively.

In (10), $\bm{u}_{a}$ represents the generalized aerodynamic force vector
acting on the wing composed of lift and drag forces which are generated by the
interaction of the wing surface with the airflow. The flapping induces a non-
uniform distribution of velocity across the wing surface, therefore the
aerodynamic forces must be calculated and then integrated across the wing
surface, as shown in Fig. 5. Due to the nature of the flapping motion, the
wingtip can also be a leading edge so we also integrate the aerodynamic forces
about the wingspan as well. The force is assumed to be concentrated at the
aerodynamic center located a quarter chord or wingspan away from the leading
edge.

The lift $d\bm{f}_{L}$ and drag $d\bm{f}_{D}$ forces for every segment can be
found by the following formulas

$\displaystyle d\bm{f}^{B}_{L}$ $\displaystyle=\tfrac{1}{2}\rho
v_{r}^{2}\,C_{L}(\beta)\,\bm{e}^{B}_{L}\,dS$ (14) $\displaystyle
d\bm{f}^{B}_{D}$ $\displaystyle=\tfrac{1}{2}\rho
v_{r}^{2}\,C_{D}(\beta)\,\bm{e}^{B}_{D}\,dS,$

where $\rho$ is the air density, $\bm{e}_{L}$ and $\bm{e}_{D}$ are the
directions of the lift and drag forces, $v_{r}$ is the relative speed of the
segment with respect to the wind speed about the axis $\bm{e}_{L}$ and
$\bm{e}_{D}$, $\beta$ is the angle of attack, $C_{L}$ and $C_{D}$ are the lift
and drag coefficients, and $dS$ the projected wing segment area as shown in
Fig. 5. The chordwise integration is defined about the segment surface area
$dS=c\,s\,d\hat{r}$ where $c$ is the airfoil chord length, $s$ the wing span
length ($s_{H}$ and $s_{R}$ for humerus and radius, respectively), and
$d\hat{r}$ the spanwise segment length. On the other hand, the spanwise
integration is defined with $dS=c\,s\,d\hat{c}$ where $d\hat{c}$ represents
the chordwise segment length.

The position where the aerodynamic forces are applied, $\bm{p}_{a}$, can be
found by using the following formulas

$\displaystyle\bm{p}_{a,H}^{B}$
$\displaystyle=\bm{p}_{7}^{B}+R_{x}(\theta_{s})\,\bm{l}_{a,H}$ (15)
$\displaystyle\bm{p}_{a,R}^{B}$
$\displaystyle=\bm{p}_{7}^{B}+R_{x}(\theta_{s})\,\left(\bm{l}_{H}^{B}+R_{x}(\theta_{e})\,\bm{l}_{a,R}\right),$

where $\bm{l}_{a,H}$ and $\bm{l}_{a,R}$ are the length vectors from their
respective joints to the aerodynamic forces, as illustrated in Fig. 5. The
relative speed between wing segment and the airflow is computed by

$v_{r}^{2}=(\bm{v}_{w}^{\top}R_{B}\bm{e}^{B}_{L})^{2}+(\bm{v}_{w}^{\top}R_{B}\bm{e}^{B}_{D})^{2},$
(16)

in which $\bm{v}_{w}=\bm{v}_{a}-\bm{v}_{\infty}$ is the relative wind speed,
where $\bm{v}_{\infty}$ is the true wind speed, and
$\bm{v}_{a}=\dot{\bm{p}}_{a}$ is the inertial velocity of the segment at the
location where the aerodynamic force is applied.

The lift and drag coefficients are defined following the results from [34]
which is experimentally tested for a MAV mimicking a fruit fly. These
coefficients are

$\displaystyle C_{L}(\beta)$
$\displaystyle=0.225+1.58\,\sin(2.13\,\beta-7.2^{\circ})$ (17) $\displaystyle
C_{D}(\beta)$ $\displaystyle=1.92-1.55\,\cos(2.04\,\beta-9.82^{\circ}),$

where $\beta$ is defined in degrees. The angle of attack $\beta$ can be
determined by using the following formula

$\beta=-\mathrm{atan2}(\bm{v}_{w}^{\top}R_{B}\bm{e}^{B}_{L},\bm{v}_{w}^{\top}R_{B}\bm{e}^{B}_{D}).$
(18)

Then the generalized aerodynamic force acting on each wing segment can be
derived as

$d\bm{u}_{a,j}(\hat{x})=\left(\frac{\partial\dot{\bm{p}}_{a,j}(\hat{x})}{\partial{\dot{\bm{q}}_{2}^{\prime}}}\right)^{\top}R_{B}\,\left(d\bm{f}^{B}_{L,j}(\hat{x})+d\bm{f}^{B}_{D,j}(\hat{x})\right),$
(19)

where $\hat{x}$ can either be $\hat{r}$ or $\hat{c}$, and the set
$j\in\mathcal{W}=\\{H_{L},H_{R},R_{L},R_{R}\\}$ is the set of all wing
segments (left/right and humerus/radius). Finally, the resulting generalized
aerodynamic forces acting on the system can be solved through the following
integration

$\displaystyle\bm{u}_{a}=\sum_{j\in\mathcal{W}}\left(\int_{0}^{1}d\bm{u}_{a,j}(\hat{r})\right)+\sum_{j\in\mathcal{W_{R}}}\left(\int_{0}^{1}d\bm{u}_{a,j}(\hat{c})\right),$
(20)

where the set $\mathcal{W_{R}}=\\{R_{L},R_{R}\\}$ is the set of the radius
wing segments.

##  IV Controller Design, Optimization, and Numerical Simulation

The wing mechanism is controlled by adjusting the flapping rate and the length
of the FDCs, which are the components of the massless subsystem, through the
input $\bm{u}_{1}$ in (5). In this view, kinetic sculptures deliver dynamic
morphing capabilities which are the key features in bats flight apparatus
while FDCs take supervisory roles to stabilize the flight dynamics. Let
$\dot{\theta}_{1}$ be the speed of the motor driven crank gear while
$\bm{l}=[l_{3b},l_{3c},l_{8b},l_{10b}]^{\top}$ be the vector containing the
length of FDCs. The flapping frequency and the FDC lengths can be adjusted by
a simple PD controller as shown below

$\displaystyle u_{g}$ $\displaystyle=K_{d1}(\omega_{ref}-\dot{\theta}_{1})$
(21) $\displaystyle\bm{u}_{p}$
$\displaystyle=K_{p2}\left(\bm{l}_{ref}-\bm{l}\right)-K_{d2}\,\dot{\bm{l}},$

where $K_{pi}$ and $K_{di}$ are the control gains, $\omega_{ref}$ is the
desired flapping frequency,
$\bm{u}_{p}=[u_{3b},u_{3c},u_{8b},u_{10b}]^{\top}$, and $\bm{l}_{ref}$ is the
desired FDC length vector. The flapping rate is set to be a constant value of
10 Hz which is the approximate flapping frequency of the Egyptian fruit bat
(rousettus aegyptiacus) which is the basis of our robot’s design. On the other
hand, the $\bm{l}_{ref}$ is found using optimization framework which will be
outlined in the following sections.

###  IV-A Open Loop Gait Optimization

Figure 6: Illustrates the optimized gait. This gait yields stable and periodic
trajectories. The gait produces adequate thrust and lift forces. However, as
it is seen above, the pitch angle is slowly drifting up.

This optimization is setup to find a gait which has a stable limit cycle by
optimizing the body initial pitch and FDC lengths. The FDC lengths are set to
be constant in this optimization, which results in a periodic flapping gait
with no feedback stabilization. The feedback into the system will be added
after this stable limit cycle is found which will be outlined in the next
section.

The cost function for the optimization is determined from a numerical
simulation done over a set period of time. Let the time evolution of the
simulated states be
$\bm{x}_{k+1}=\bm{f}(\bm{x}_{k},\bm{l}_{ref,k},\bm{u}_{k})$, where the
simulation states $\bm{x}_{k}$ contains all of the relevant states from the
equation of motions derived in Section III. The angular momentum is used as a
metric to stabilize the flapping gait which is derived using the following
formulations

$\bm{\Pi}=\textstyle\sum_{i\in\mathcal{L_{M}}}R_{B}\,I_{i}\,\omega_{i}-m_{i}(\bm{x}_{i}-\bm{x}_{CoM})\times\dot{\bm{x}}_{i},$
(22)

where $\bm{x}_{CoM}$ represents the averaged CoM position of the entire robot.
Then the stable limit cycle is found by solving the following optimization
problem:

$\displaystyle\min_{\bm{l}_{ref},\theta_{y}}$ $\displaystyle\quad
J=\sum^{N}_{k=1}(w_{1}\,\bm{\Pi}_{k}^{\top}\,\bm{\Pi}_{k}+w_{2}\,\dot{\bm{x}}_{B,k}^{\top}\,\dot{\bm{x}}_{B,k})\,\Delta
t$ (23) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\bm{l}_{min}\leq\bm{l}_{ref}\leq\bm{l}_{max},$

where the cost $J$ is setup to minimize the angular momentum $\bm{\Pi}$ and
body velocity $\dot{\bm{x}}_{B}$, $\theta_{y}$ is the initial pitch angle,
$\bm{l}_{min}$ and $\bm{l}_{max}$ are the bounds on $\bm{l}_{ref}$, $\Delta t$
is the simulation time step, $N$ is the simulation final step, and $w_{i}$ is
the weight factor in the cost function. The bounds are set at
$\bm{l}_{min}=0.8\,\bm{l}_{0}$ and $\bm{l}_{max}=1.2\,\bm{l}_{0}$, where
$\bm{l}_{0}$ is the initial lengths of FDCs in the wing structure. The robot
is initialized at rest and subjected to a wind speed of
$\bm{v}_{\infty}=[-2,0,0]^{\top}$ m/s. The wing mechanism is set at an offset
of -10 cm from the body CoM about the body $x$-axis and the wing chord length
is set at 20 cm. The remaining parameters are the same as the kinetic
sculpture shown in Fig. 2-A and our previous work in [26].

The optimization in (23) is run using an RK4 simulation with a simulation end
time of 1 second. The gait optimization result can be seen in Fig. 6, where a
relatively stable gait is found with the optimal parameters
$\bm{l}_{ref}=[7.8,10.5,6.2,7.2]$ mm and $\theta_{y}=33^{\circ}$. The
simulation result shows that this gait produces positive lift and thrust while
having a periodic and approximately constant velocity of $[-1.52,0,-0.96]$
m/s. A longer simulation time of 4 seconds is also shown in Fig. 6 where the
pitch angle is shown to drifts slowly upwards and the robot is also losing
forward speed over time. The drift in the pitch angle will be stabilized by
using the FDCs which will be outlined in the next section.

###  IV-B Pitch Stabilization Optimization

Figure 7: Shows the closed-loop simulation results based on the MIMIC
framework, i.e., integrating mechanical intelligence and control. The flight
controller has successfully stabilized the pitch angle after a transient
period of approximately 2 seconds.

A pitch stabilization controller can then be implemented to regulate the slow
pitch drift that exist in the optimized gait. We consider the following
controller

$\displaystyle\bm{l}_{ref}=\bm{l}_{ref,zp}+K_{c}\,(\theta_{y,ref}-\theta_{y}),$
(24)

where the $\bm{l}_{ref,zp}$ is the constant zero-path flight reference found
in the gait optimization (23), $\theta_{y,ref}$ is the pitch angle reference,
and $K_{c}$ is the controller gain matrix. The gain for the controller in (24)
can be found using the following optimization

$\displaystyle\min_{K_{c}}$
$\displaystyle\quad\textstyle{J=\sum^{N}_{k=1}(w_{1}\,\bm{\Pi}_{k}^{\top}\,\bm{\Pi}_{k}+w_{2}\,\dot{\bm{x}}_{B,k}^{\top}\,\dot{\bm{x}}_{B,k}}$
(25)
$\displaystyle\qquad\qquad\qquad+w_{3}\,(\theta_{y,ref}-\theta_{y_{k}})^{2})\,\Delta
t$ $\displaystyle\mathrm{s.t.}$ $\displaystyle\quad K_{c,min}\leq K_{c}\leq
K_{c,max}.$

This optimization has an additional cost $w_{3}$ for weighting the robot’s
pitch versus a constant pitch reference. The FDC lengths are also constrained
using the saturation $\bm{l}_{min}\leq\bm{l}_{ref}\leq\bm{l}_{max}$ to prevent
the lengths from going unbounded.

The optimization found an optimal controller gain of
$K_{c}=[0.42,-0.26,-0.38,-0.097]^{\top}$ and the simulation result is shown in
Fig. 7 which has a transient state up to $t=2$s before reaching a steady limit
cycle. The pitch is stable near the target pitch angle of $33^{\circ}$
throughout the simulation. This shows that the controller has successfully
achieved pitch stabilization by utilizing the FDCs.

##  V Conclusions and Future Work

Copying bat high-dimensional flight apparatus is nearly impossible by
considering classical feedback design paradigms based on sensing, computing
and actuation. Bat musculoskeletal system possesses many active and passive
joints, e.g., bat wings possess over 40 joints. In this work, we offered a
solution towards the robotic biomimicry of bat aerial locomotion. We proposed
a design framework called Morphing via Integrated Mechanical Intelligence and
Control (MIMIC). We leveraged computational structures called Kinetic
Sculptures (KS) to subsume part of the responsibility of closed-loop feedback
under mechanical intelligence. Then, we extended our previous works by
considering Feedback-Driven Components (FDC) in the design of KS which possess
supervisory roles and are used for flight stabilization. We used the dynamical
model of Northeastern University’s Aerobat to show the successful
stabilization of the Aerobat’s longitudinal dynamics. This robot possesses an
articulated wing structure with many active and passive joints. In addition,
the robot is tail-less which means its longitudinal dynamics are open-loop
unstable.

## References

  * [1] J.-N. Lee and K.-C. Kwak, “A trends analysis of image processing in unmanned aerial vehicle,” _International Journal of Computer, Information Science and Engineering_ , vol. 8, no. 2, pp. 2–5, 2014.
  * [2] I. Pavlidis, V. Morellas, P. Tsiamyrtzis, and S. Harp, “Urban surveillance systems: from the laboratory to the commercial world,” _Proceedings of the IEEE_ , vol. 89, no. 10, pp. 1478–1497, 2001.
  * [3] “Smart Cities Market | Growth, Trends, Forecast (2020-2025).” [Online]. Available: https://www.mordorintelligence.com/industry-reports/smart-cities-market
  * [4] I. Tiddi, E. Bastianelli, E. Daga, M. d’Aquin, and E. Motta, “Robot–city interaction: Mapping the research landscape—a survey of the interactions between robots and modern cities,” _International Journal of Social Robotics_ , vol. 12, no. 2, pp. 299–324, 2020.
  * [5] T. Koolen, S. Bertrand, G. Thomas, T. de Boer, T. Wu, J. Smith, J. Englsberger, and J. Pratt, “Design of a Momentum-Based Control Framework and Application to the Humanoid Robot Atlas,” _International Journal of Humanoid Robotics_ , vol. 13, no. 01, p. 1650007, Mar. 2016.
  * [6] S. Fahmi, C. Mastalli, M. Focchi, and C. Semini, “Passive Whole-Body Control for Quadruped Robots: Experimental Validation Over Challenging Terrain,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 3, pp. 2553–2560, Jul. 2019.
  * [7] C. Dario Bellicoso, C. Gehring, J. Hwangbo, P. Fankhauser, and M. Hutter, “Perception-less terrain adaptation through whole body control and hierarchical optimization,” in _2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids)_ , Nov. 2016, pp. 558–564, iSSN: 2164-0580.
  * [8] T. Bretl and S. Lall, “Testing Static Equilibrium for Legged Robots,” _IEEE Transactions on Robotics_ , vol. 24, no. 4, pp. 794–807, Aug. 2008\.
  * [9] D. Pardo, L. Möller, M. Neunert, A. W. Winkler, and J. Buchli, “Evaluating Direct Transcription and Nonlinear Optimization Methods for Robot Motion Planning,” _IEEE Robotics and Automation Letters_ , vol. 1, no. 2, pp. 946–953, Jul. 2016.
  * [10] C. Mastalli, M. Focchi, I. Havoutis, A. Radulescu, S. Calinon, J. Buchli, D. G. Caldwell, and C. Semini, “Trajectory and foothold optimization using low-dimensional models for rough terrain locomotion,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2017, pp. 1096–1103.
  * [11] C. Mastalli, I. Havoutis, A. W. Winkler, D. G. Caldwell, and C. Semini, “On-line and on-board planning and perception for quadrupedal locomotion,” in _2015 IEEE International Conference on Technologies for Practical Robot Applications (TePRA)_ , May 2015, pp. 1–7.
  * [12] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and Trajectory Optimization for Legged Systems Through Phase-Based End-Effector Parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, Jul. 2018.
  * [13] B. Aceituno-Cabezas, C. Mastalli, H. Dai, M. Focchi, A. Radulescu, D. G. Caldwell, J. Cappelletto, J. C. Grieco, G. Fernández-López, and C. Semini, “Simultaneous Contact, Gait, and Motion Planning for Robust Multilegged Locomotion via Mixed-Integer Convex Optimization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2531–2538, Jul. 2018.
  * [14] H. Dai and R. Tedrake, “Planning robust walking motion on uneven terrain via convex optimization,” in _2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids)_ , Nov. 2016, pp. 579–586, iSSN: 2164-0580.
  * [15] D. Arterburn, M. Ewing, R. Prabhu, F. Zhu, D. Francis _et al._ , “Faa uas center of excellence task a4: Uas ground collision severity evaluation.” 2017\.
  * [16] H. Tanaka, H. Okada, Y. Shimasue, and H. Liu, “Flexible flapping wings with self-organized microwrinkles,” _Bioinspiration & Biomimetics_, vol. 10, no. 4, p. 046005, 2015.
  * [17] D. K. Riskin, D. J. Willis, J. Iriarte-Díaz, T. L. Hedrick, M. Kostandov, J. Chen, D. H. Laidlaw, K. S. Breuer, and S. M. Swartz, “Quantifying the complexity of bat wing kinematics,” _Journal of Theoretical Biology_ , vol. 254, no. 3, pp. 604–615, 2008.
  * [18] A. Azuma, _The biokinetics of flying and swimming_. Springer Science & Business Media, 2012.
  * [19] J. W. Bahlman, S. M. Swartz, and K. S. Breuer, “Design and characterization of a multi-articulated robotic bat wing,” _Bioinspiration & Biomimetics_, vol. 8, no. 1, p. 016009, 2013.
  * [20] J. Colorado, A. Barrientos, C. Rossi, and K. S. Breuer, “Biomechanics of smart wings in a bat robot: morphing wings using sma actuators,” _Bioinspiration & biomimetics_, vol. 7, no. 3, p. 036006, 2012.
  * [21] A. Ramezani, S.-J. Chung, and S. Hutchinson, “A biomimetic robotic platform to study flight specializations of bats,” _Science Robotics_ , vol. 2, no. 3, 2017.
  * [22] J. Hoff, A. Ramezani, S.-J. Chung, and S. Hutchinson, “Synergistic design of a bio-inspired micro aerial vehicle with articulated wings.” in _Robotics: science and systems_ , 2016.
  * [23] J. Hoff, A. Ramezani, S.-J. Chung, and S. Hutchinson, “Optimizing the structure and movement of a robotic bat with biological kinematic synergies,” _The International Journal of Robotics Research_ , vol. 37, no. 10, pp. 1233–1252, 2018.
  * [24] J. Hoff, A. Ramezani, S.-J. Chung, and S. Hutchinson, “Reducing versatile bat wing conformations to a 1-dof machine,” in _Conference on Biomimetic and Biohybrid Systems_. Springer, 2017, pp. 181–192.
  * [25] A. Ramezani, “Towards biomimicry of a bat-style perching maneuver on structures: the manipulation of inertial dynamics,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 7015–7021.
  * [26] E. Sihite, P. Kelly, and A. Ramezani, “Computational structure design of a bio-inspired armwing mechanism,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5929–5936, 2020.
  * [27] A. Ramezani, X. Shi, S.-J. Chung, and S. Hutchinson, “Lagrangian modeling and flight control of articulated-winged bat robot,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 2867–2874.
  * [28] J. Hoff, U. Syed, A. Ramezani, and S. Hutchinson, “Trajectory planning for a bat-like flapping wing robot,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 6800–6805.
  * [29] A. Ramezani, X. Shi, S.-J. Chung, and S. Hutchinson, “Nonlinear flight controller synthesis of a bat-inspired micro aerial vehicle,” in _AIAA Guidance, Navigation, and Control Conference_ , 2016, p. 1376.
  * [30] A. Ramezani, S. U. Ahmed, J. Hoff, S.-J. Chung, and S. Hutchinson, “Describing robotic bat flight with stable periodic orbits,” in _Conference on Biomimetic and Biohybrid Systems_. Springer, 2017, pp. 394–405.
  * [31] U. A. Syed, A. Ramezani, S.-J. Chung, and S. Hutchinson, “From rousettus aegyptiacus (bat) landing to robotic landing: Regulation of cg-cp distance using a nonlinear closed-loop feedback,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3560–3567.
  * [32] E. Sihite and A. Ramezani, “Enforcing nonholonomic constraints in aerobat, a roosting flapping wing model,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 5321–5327.
  * [33] H. Hauser, A. J. Ijspeert, R. M. Füchslin, R. Pfeifer, and W. Maass, “The role of feedback in morphological computation with compliant bodies,” _Biological cybernetics_ , vol. 106, no. 10, pp. 595–613, 2012.
  * [34] S. P. Sane and M. H. Dickinson, “The control of flight force by a flapping wing: lift and drag production,” _Journal of experimental biology_ , vol. 204, no. 15, pp. 2607–2626, 2001.


# Model-independent partial wave analysis using a massively-parallel fitting
framework

L Sun1 R Aoude2 A C dos Reis2 M Sokoloff3 1 School of Physics and Technology,
Wuhan University, Wuhan, Hubei Province, 430072 China  2 Centro Brasileiro de
Pesquisas Físicas (CBPF), Rio de Janeiro, 22290-180 Brazil  3 University of
Cincinnati, Physics Department, ML0011, Cincinnati OH 45221-0011, USA
<EMAIL_ADDRESS>
###### Abstract

The functionality of GooFit, a GPU-friendly framework for doing maximum-
likelihood fits, has been extended to extract model-independent $\cal S$-wave
amplitudes in three-body decays such as $D^{+}\rightarrow h^{+}h^{+}h^{-}$. A
full amplitude analysis is done where the magnitudes and phases of the $\cal
S$-wave amplitudes are anchored at a finite number of $m^{2}(h^{+}h^{-})$
control points, and a cubic spline is used to interpolate between these
points. The amplitudes for $\cal P$-wave and $\cal D$-wave intermediate states
are modeled as spin-dependent Breit-Wigner resonances. GooFit uses the Thrust
library, with a CUDA backend for NVIDIA GPUs and an OpenMP backend for threads
with conventional CPUs. Performance on a variety of platforms is compared.
Executing on systems with GPUs is typically a few hundred times faster than
executing the same algorithm on a single CPU.

##  1 Introduction

The physics of scalar mesons below 2 $\mathrm{\,Ge\kern-1.00006ptV}$ is a
long-standing puzzle in light meson spectroscopy [1], this is partly due to
the fact that these mesons often have large widths and overlap with each
other, which make them hard to model. Contributions from Light scalars can be
extracted by using Dalitz plot analyses [2] in three-body hadronic decays of
charm mesons. First developed by the E791 Collaboration [3], the Model
Independent Partial Wave Analysis (MIPWA) technique extracts $\cal S$-wave
amplitudes in the Dalitz plot analysis with no assumption about the nature of
the $\cal S$-wave. To pin down its large numbers of free parameters, the MIPWA
technique requires the large samples of three-body decay events that have
become increasingly available at the $B$-factories and the LHC. Notably, the
LHCb experiment has recorded charm decays with sample sizes exceeding any
previous experiment by more than an order of magnitude, offering a unique
opportunity to study $\cal S$-wave structures with unprecedented levels of
precision. Analysing very large statistics samples requires disproportionately
more computing power. Running all the calculations in a single Central
Processing unit (CPU) thread would take a prohibitively long time.

Originally, a Graphics Processing Unit (GPU) was a specialised electronic
circuit designed to rapidly manipulate and alter memory to accelerate the
creation of images for output to a video display. The highly parallel
structure makes GPUs more effective than CPUs for algorithms where large
blocks of data are processed in parallel. Today, the functionality of GPUs has
been extended to enable highly parallel computing for scientific and other
more general applications. An open-source framework called GooFit [4] has been
developed to exploit the processing power of these GPUs for parallel function
evaluation, particularly in the context of maximum likelihood fits. This paper
reports an extension of the GooFit framework to support MIPWA of a three-body
decay with vastly reduced processing time.

##  2 MIPWA method

Essentially all studies of three-body hadronic $D_{(s)}$ decays employ the
same technique: the unbinned maximum likelihood fit of the Dalitz plot, in
which the quantum mechanical matrix element governing the decay is represented
by a coherent sum of amplitudes [1]. These amplitudes correspond to the
possible intermediate states in the decay chain $D_{(s)}\rightarrow
Rh_{3},R\rightarrow h_{1}h_{2}$ ($h=K,\ \pi$). The amplitudes are grouped
according to the orbital angular momentum $L$ of $R$ and $h_{3}$ in the rest
frame of the $D_{(s)}$, and the total amplitude is

$A(s_{12},s_{13})=\sum_{L}\sum_{k}c_{k}^{L}A_{k}^{L}(s_{12},s_{13})\,,$ (1)

where $s_{12(13)}\equiv m^{2}(h_{1}h_{2(3)})$. The amplitudes $A_{k}^{L}$ are
weighted by constant complex coefficients $c_{k}^{L}$, and the series is
truncated at $L=2$. For a resonance with non-zero spin, the amplitude $A_{k}$
is most often described using a relativistic Breit-Wigner function multiplied
by a real spin-dependent angular factor [5]. For the $\cal S$-wave, two
qualitatively different approaches exist. In the Isobar model, the $\cal
S$-wave is treated as a sum of a constant non-resonant term and Breit-Wigner
functions for the scalar resonances. While the Isobar model provides
reasonably good descriptions for narrow resonances, it fails to describe the
overlap of broad resonances. Additionally, the physical interpretation of the
constant non-resonant term is problematic.

To address these issues, the MIPWA describes the ensemble of scalar components
using a purely phenomenological set of parameters derived from the data. The
$s_{12(13)}$ mass spectrum is divided into $N-1$ slices with $N$ boundary
points separating the slices and at the two ends of the spectrum. The $\cal
S$-wave is represented by a generic complex function
$A_{0}(s)=a_{0}(s)e^{i\phi_{0}(s)}$. At each of the $k$ boundary points,
$A_{0}(s=s_{k})=a_{k}e^{i\phi_{k}}$ where $a_{k}$ and $\phi_{k}$ are real
parameters. Between the $N$ boundary points, the $\cal S$-wave is parametrised
by a cubic spline [6] in the complex plane. The set of {$a_{k},\phi_{k}$} are
free parameters, along with the coefficients $c_{l}^{L}$ of the higher spin
terms, are determined in the MIPWA fit. The large data sets studied, along
with the large number of parameters ($2N$) required to describe the $\cal
S$-wave, make maximising the likelihood a computationally intensive problem.
Interpolating between the boundary points leads to correlations between the
parameters, and hence to non-linear behavior.

##  3 MIPWA method with GooFit

GooFit provides an interface to allow probability density functions (PDFs) to
be evaluated in parallel, using either GPUs or multicore CPUs as back-ends.
While the original intention of GooFit was to utilise the massive
computational power of NVIDIA GPUs based on the proprietary Compute Unified
Device Architecture language (CUDA) [7], the Thrust parallel algorithms
library [8] also supports an OpenMP backend for conventional CPUs.

GooFit has been used to perform a number of amplitude analyses, for example
that of Ref. [9]. For a time-integrated Dalitz plot analysis, GooFit provides
the DalitzPlotPdf class to model the Dalitz plot PDF. The DalitzPlotPdf object
contains a list of ResonancePdf objects to describe resonant amplitudes as
well as a constant nonresonant amplitude. The original GooFit package provided
only the Isobar model to parameterise the $\cal S$-wave.

The ResonancePdf class has been extended to add support for the MIPWA method.
The entire $\cal S$-wave is treated as a single resonance (ResonancePdf
object), with a set of free parameters {$a^{k}_{0},\phi^{k}_{0}$} to be
determined in the fit.

##  4 MIPWA in $D^{+}\rightarrow K^{+}K^{+}K^{-}$ decay 111Charge conjugation
is implied throughout.

The $D^{+}\rightarrow K^{+}K^{+}K^{-}$ decay offers a good opportunity to
study the $K^{+}K^{-}$ $\cal S$-wave amplitude directly. LHCb collected a
large $D^{+}\rightarrow K^{+}K^{+}K^{-}$ data sample from
$2.1~{}\mbox{\,fb}^{-1}$ of $\sqrt{s}=7$ $\mathrm{\,Te\kern-1.00006ptV}$ $pp$
collisions recorded by the experiment during 2012. About 100K candidates were
selected with a signal purity of 90%. The resonant structure of the
$K^{+}K^{-}$ $\cal S$-wave amplitude was studied for the first time using an
Isobar model, and the results were first presented during the CHARM 2016
workshop [10]. This analysis indicates that the $\cal S$-wave component
accounts for about 90% of the decay rate, while the $\cal P$-wave
($\phi(1020)$ resonance) makes up the rest.

A specific goal of the work presented here is to provide a tool to extract the
$K^{+}K^{-}$ $\cal S$-wave by performing an MIPWA on the same LHCb data
sample. In this case, the Dalitz plot is symmetrised along the two axes of
$m^{2}(K^{+}K^{-})$. The $K^{+}K^{-}$ mass squared range is divided into 39
slices using 40 boundary points. Because the final state contains two
indistinguishable $K^{+}$ mesons, the $K^{+}K^{-}$ amplitudes interfere with
themselves, and this sort of interference allows the $\cal S$-wave phase to be
determined. To test the GooFit extension to support MIPWA fits, samples of
100K signal events are generated from a Dalitz plot PDF (“Toy MC”). The mass
projections from the fit to a test sample can be seen in Fig. 1. First, the
fit quality of the MIPWA method is tested by comparing the fitted
{$a_{k},\phi_{k}$} values with the inputs. In each fit iteration, the Dalitz
plot normalisation method (“DalitzPlotPdf::normalise()”) is called so that the
integral of the total PDF is equal to one. The normalisation integral is
calculated numerically based on evenly distributed grid points in the Dalitz
plot plane. Figure 2 shows the improvement in the fit quality as the
normalisation grid spacing is reduced. Although the finer granularity
increases numbers of calculations, the high speed of the GPUs makes this
problem tractable.

Figure 1:  The projections of $m^{2}(K^{+}K^{-})$ (left) and $m^{2}(K^{+}K^{+})$ (right) from $D^{+}\rightarrow K^{+}K^{+}K^{-}$. In each plot, the toy MC signal events (points with error bars ) are shown together with the total fit (blue line), $\phi$ resonance (red line), and $\cal S$-wave (magenta line) determined from the MIPWA. Figure 2:  Differences between the fitted values and input ones for each boundary point. Shown are set of points with different grid spacings for the normalisation: 0.01 ${\mathrm{\,Ge\kern-1.00006ptV}}^{2}$ (purple), 0.004 ${\mathrm{\,Ge\kern-1.00006ptV}}^{2}$ (red), and 0.001 ${\mathrm{\,Ge\kern-1.00006ptV}}^{2}$ (blue). Table 1: GPU performance of the MIPWA fit. Platform | GPU Model | Chip | CUDA cores | Run time (sec.)
---|---|---|---|---
tWorkstation | Tesla K40c | GK110BGL | 2880 | 76
Desktop PC | GeForce GTX 980 | 2nd gen. Maxwell (GM204) | 2048 | 67
Laptop ASUS N56V | GeForce GT 650M | GK107 | 384 | 179

As shown in Table 1, the GooFit performance on three different GPU platforms
has been measured by running the MIPWA fit over the same toy MC sample of 100
K events (see Fig. 1 for the fit projections). For comparison, the MIPWA fit
using an older code with the same functionality that runs on one CPU core
takes about eight hours to complete. Perhaps surprisingly, an older generation
mobile GPU (a GeForce GT 650M with 384 cores) provides excellent performance;
a newer HPC GPU board (a Tesla K40c with 2880 cores) provides better
performance, but not in proportion to the number of cores; a high-end gamer
board (a GeForce GTX 980) provides the best performance, albeit by a small
margin. Based on other studies with GooFit, significantly better performance
on the new P100 boards is anticipated which utilise NVIDIA’s new Pascal GPU
architecture.

Table 2: Specifications of the testing CPU platforms. Asterisks next to the number of cores indicate hyperthreading - two virtual processors per physical core. Name | Chip type | # of Cores | Clock [GHz] | RAM [GiB]
---|---|---|---|---
Cerberus | Intel Xeon E5520 | 8* | 2.27 | 24
Goofy | Intel Xeon CPU E5-2680 v3 | 24* | 2.50 | 120
Figure 3: Timing results and speedups for the test MIPWA fit. Tested on two
different CPU platforms as shown in Table 2: Goofy (Intel Xeon CPU E5-2680 v3
x 2) and Cerberus (Intel Xeon CPU E5520 x 2).

In addition to timing GooFit’s new MIPWA performance on GPUs using Thrust’s
CUDA backend, its performance on two different CPU platforms is timed using
Thrust’s OpenMP backend, as shown in Fig. 3. With two Intel Xeon E5-2680 v3
CPUs, the fit uses 791 seconds for one OpenMP thread, 50 seconds for 24
threads. The speedup is almost linear with the increase of the number of
OpenMP threads, up until the number of threads equals the numbers of physical
cores.

##  5 Summary

This paper describes an extension of GooFit to support MIPWA fits for three-
body decays, and have achieved speedups of a few hundred by using GPUs. The
main branch of GooFit’s source code is in a GooFit repository at
https://github.com/GooFit/GooFit, while the updated code with MIPWA support is
in a personal GooFit branch at https://github.com/liang-sun/GooFit.

###  5.1 Acknowledgments

This work was performed with support from NSF Award PHY-1414736. NVidia
provided K40 GPUs for our use through its University Partnership program. The
Ohio Supercomputer Center made their “Oakley” computer farm available for
development, for testing, and for GooFit outreach workshops.

## References

## References

  * [1] C. Patrignani et al. (Particle Data Group), Chin. Phys. C 40, 100001 (2016).
  * [2] Dalitz R H, 1953 Phil. Mag. 44 1068\.
  * [3] E. M. Aitala et al. (E791 Collaboration), Phys. Rev. D 73, 032004 (2006).
  * [4] R. Andreassen et al., IEEE Access 2, 160 (2014).
  * [5] S. Kopp et al. (CLEO Collaboration), Phys. Rev. D 63, 092001 (2001).
  * [6] W. H. Press et al., Numerical Recipes: The Art of Scientific Computing, 3rd Edition, Cambridge University Press, 2007.
  * [7] NVIDIA Corporation, 2015, NVIDIA CUDA C programming guide, Version 7.
  * [8] Nathan Bell and Jared Hoberock, Thrust: A Productivity-Oriented Library for CUDA, available from the Documentation link at https://thrust.github.io/, accessed 29 January, 2017.
  * [9] J. P. Lees et al. (BABAR Collaboration), Phys. Rev. D 93, 112014 (2016).
  * [10] LHCb Collaboration, LHCb-CONF-2016-008 (2016).


# A lower–bound for the number of conjugacy classes of $A_{n}$

Xandru Mifsud  Department of Mathematics, University of Malta, Msida, Malta
<EMAIL_ADDRESS>
###### Abstract

We establish a sharp lower–bound for the number of conjugacy classes
$k(A_{n})$ in the alternating group $A_{n}$, for $n\geq 3$. Namely, we show
that
$k\left(A_{n}\right)\geq\frac{k\left(A_{7}\right)}{\log_{2}\left|A_{7}\right|}\cdot\log_{2}\left|A_{n}\right|$
with equality if, and only if, $n=7$. The observations leading towards this
result were obtained through a genetic algorithm developed to search for
groups having certain properties.

###### keywords:

finite group, conjugacy class, alternating group

###### MSC:

[2020] 20E45

margin=1in

##  1 Introduction

One of the problems in group theory is the classification of finite groups by
the number of conjugacy classes. Given a group $G$, we denote its number of
conjugacy classes by $k(G)$ and its order by $|G|$. Unless otherwise stated,
all groups considered throughout this work are finite. All groups with
$k(G)\leq 14$ have been classified, with the cases for $k(G)=13$ and $k(G)=14$
being resolved recently in [1]. Consequently, bounding $k(G)$ is an important
problem in group theory.

A number of bounds for $k(G)$ have been established in the literature.
Recently, Maróti [2] showed that if $G$ is a finite group whose order is
divisible by a prime $p$, then $k(G)\geq 2\sqrt{p-1}$. Of notable interest are
bounds establishing explicit relations between $k(G)$ and $|G|$. Erdős and
Turán had established that $k(G)\geq\log_{2}\log_{2}|G|$. One of the main
problems in this area is resolving the following conjecture on the relation
between $k(G)$ and $\log_{2}|G|$.

######  Conjecture 1.1.

There exists a constant $C>0$ such that for any finite group $G$, $k(G)\geq
C\log_{2}|G|$.

One of the sharpest bounds obtained in this regard is the Pyber–Keller bound,
first established by Pyber [3] and then further sharpened by Keller [4].

######  Theorem 1.2 ([3, 4]).

There exists a (explicitly computable) constant $C>0$ such that every finite
group $G$ with $|G|\geq 4$ satisfies

$k(G)\geq\dfrac{C}{\left(\log_{2}\log_{2}|G|\right)^{7}}\cdot\log_{2}|G|.$

In light of Conjecture 1.1, we prove the following sharp–bound of this type,
for alternating groups.

######  Theorem 1.3.

Let $n\geq 3$. Then,

$k\left(A_{n}\right)\geq\dfrac{k\left(A_{7}\right)}{\log_{2}\left|A_{7}\right|}\cdot\log_{2}\left|A_{n}\right|$

with equality if, and only if, $n=7$.

###  1.1 Overview

We begin by introducing some preliminaries in Section 2, primarily some
textbook results on the relations between the number of conjugacy classes of
$S_{n}$ and the number of partitions of a positive integer $n$, denoted by
$p(n)$. This is followed by a proof of Theorem 1.3 in Section 3, and a short
discussion in Section 4 on the result and some open problems.

##  2 Preliminaries: Relations between $k(S_{n})$, $k(A_{n})$ and $p(n)$

A classical result in the study of $S_{n}$ is that the conjugacy classes of
$S_{n}$ correspond to the permutation cycle structures, and in turn to the
partitions of $n$. Hence $k(S_{n})=p(n)$. While a generating–function for the
partition function $p(n)$ is known, no closed–form expression is known.
Consequently, the study of (asymptotic) bounds and recurrence relations for
$p(n)$ represents a highly active area of research in number theory. Maróti
[5] established a number of lower bounds for $p(n)$, namely the following,
which will be useful in our proof for Theorem 1.3.

######  Theorem 2.1 ([5]).

For every positive integer $n$,

$\dfrac{e^{2\sqrt{n}}}{14}<p(n).$

Recall that $A_{n}$ is the group of even permutations in $S_{n}$. The
splitting criterion for conjugacy classes in $A_{n}$ states that given a
conjugacy class of even permutations in $S_{n}$: the conjugacy class is either
a conjugacy class in $A_{n}$ or else it ‘splits’ into two conjugacy classes in
$A_{n}$, with splitting occurring if and only if the cycle structure
corresponding to the conjugacy class consists of cycles having distinct odd
lengths. Consequently, the number of conjugacy classes in $S_{n}$ that split
in $A_{n}$ is equal to the number of partitions of $n$ into distinct odd
parts, which we term as $\alpha(n)$.

Let $\beta(n)$ be the number of conjugacy classes of even permutations in
$S_{n}$ that do not split in $A_{n}$. A classical result relating $k(A_{n})$
to $k(S_{n})$, via $\alpha(n)$ and $\beta(n)$, is the following.

######  Theorem 2.2.

Let $n\geq 3$. Then, $k(S_{n})=\alpha(n)+2\beta(n)$ and
$k(A_{n})=2\alpha(n)+\beta(n)$.

Consequently, we get the following corollary.

######  Corollary 2.3.

Let $n\geq 3$. Then, $\dfrac{k(S_{n})}{2}\leq k(A_{n})$.

###### Proof.

By Theorem 2.2, we have $\dfrac{k(S_{n})}{2}=\dfrac{\alpha(n)}{2}+\beta(n)\leq
2\alpha(n)+\beta(n)=k(A_{n})$, since $\alpha(n)\geq 0$. ∎

Further aspects of these results are given by Dénes, Erdős and Turán [6].

##  3 Proof of Theorem 1.3

We are now in a position to prove our main result. We will show that

$\dfrac{k(A_{7})}{\log_{2}|A_{7}|}=\min_{n\geq
3}\dfrac{k(A_{n})}{\log_{2}|A_{n}|}.$ (3.1)

For $3\leq n\leq 11$, by explicit computation111Using GAP [7]: gap> for n in
[3..11] do Print(n, ” : ”, NrConjugacyClasses(AlternatingGroup(n)) /
Log2(Float(Order(AlternatingGroup(n)))), ”n”); od; (see Table 1), we have that
$\frac{k(A_{7})}{\log_{2}|A_{7}|}\leq\frac{k(A_{n})}{\log_{2}|A_{n}|}$.

$n$ | $\frac{k(A_{n})}{\log_{2}|A_{n}|}$ | $n$ | $\frac{k(A_{n})}{\log_{2}|A_{n}|}$ | $n$ | $\frac{k(A_{n})}{\log_{2}|A_{n}|}$
---|---|---|---|---|---
3 | $1.892\dots$ | 6 | $0.824\dots$ | 9 | $1.030\dots$
4 | $1.115\dots$ | 7 | $0.796\dots$ | 10 | $1.154\dots$
5 | $0.846\dots$ | 8 | $0.979\dots$ | 11 | $1.278\dots$
Table 1: Values of $\frac{k(A_{n})}{\log_{2}|A_{n}|}$ for $3\leq n\leq 11$,
computed using GAP [7].

We prove the case when $n\geq 12$. Firstly note that
$\log_{2}|A_{n}|\leq(n-2)\log_{2}n$ since,

$\log_{2}\frac{n!}{2}=-1+\sum\limits_{r=1}^{n}\log_{2}r=\sum\limits_{r=3}^{n}\log_{2}r\leq(n-2)\log_{2}n$

noting that $\log_{2}(x)$ is monotonically increasing, $\log_{2}1=0$ and
$\log_{2}2=1$.

Using this inequality and Corollary 2.3, we obtain

$\dfrac{k(A_{n})}{\log_{2}|A_{n}|}\geq\dfrac{k(S_{n})}{2(n-2)\log_{2}n}$

and by Theorem 2.1 and the fact that $p(n)=k(S_{n})$, we have that:

$\dfrac{k(A_{n})}{\log_{2}|A_{n}|}>t(n),\ \ \
t(n)=\dfrac{e^{2\sqrt{n}}}{28(n-2)\log_{2}n}.$ (3.2)

Since, for large enough $n$, $28(n-2)\log_{2}n<<e^{2\sqrt{n}}$, then in
particular there exists a positive integer $N$ such that for all $n\geq N$,
$1\leq t(N)\leq t(n)$. In particular for all $n\geq 12$,

$1<1.016\ldots=t(12)\leq t(n)$

and consequently by (3.2) and the case for $A_{7}$ in Table 1,

$\dfrac{k(A_{n})}{\log_{2}|A_{n}|}>1>\dfrac{k(A_{7})}{\log_{2}|A_{7}|}$

and therefore (3.1) follows.

By the strict inequality in (3.2) and Table 1, it follows that equality holds
if, and only if, $n=7$. ∎

##  4 Concluding remarks

  1. (1)

This lower–bound was initially suggested through empirical analysis, using a
genetic algorithm developed using the Python API of SageMath [8]. Davies et
al. [9] recently highlighted the effectiveness of guiding mathematical
intuition in this manner, through the use of artificial intelligence.

  2. (2)

Since the number of characters of an Abelian group $G$ is $|G|$ and
$\frac{k(A_{7})}{\log_{2}|A_{7}|}<1$, then when $G$ is non–trivial we have
that $\frac{k(G)}{\log_{2}|G|}\geq 1>\frac{k(A_{7})}{\log_{2}|A_{7}|}$.

  3. (3)

Based on the result of Theorem 1.3 and (2) above, as well as further empirical
analysis carried out using our algorithm, we hypothesise that the constant $C$
in Conjecture 1.1 is ‘close’ to $\frac{k(A_{7})}{\log_{2}|A_{7}|}$, such that
it holds with high probability. Hence we pose the following problem.

######  Problem 4.1.

Which non–Abelian finite groups $G$ satisfy
$\frac{k(G)}{\log_{2}|G|}<\frac{k(A_{7})}{\log_{2}|A_{7}|}$?

Case–by–case analysis yields that $M_{11}$ (the Mathieu group of degree 11) is
the smallest such group satisfying Problem 4.1.

If $G$ is a group as in Problem 4.1 which embeds into $A_{n}$ for some
$n\in\mathbb{N}$, then $k(G)$ is a lower–bound for $k(A_{n})$, as follows.

######  Corollary 4.2.

Let $G$ be a finite non–Abelian subgroup of $A_{n}$ for some $n\in\mathbb{N}$.
If $\frac{k(G)}{\log_{2}|G|}<\frac{k(A_{7})}{\log_{2}|A_{7}|}$, then
$k(G)<k(A_{n})$.

###### Proof.

Follows immediately from the fact that $|G|\leq|A_{n}|$ and (3.1), namely
since:

$k(G)<\log_{2}|G|\cdot\dfrac{k(A_{7})}{\log_{2}|A_{7}|}\leq\log_{2}|G|\cdot\dfrac{k(A_{n})}{\log_{2}|A_{n}|}<k(A_{n}).$

∎

## Acknowledgements

The author would like to thank Prof. Josef Lauri for introducing the author to
Conjecture 1.1 and Adriana Baldacchino for their helpful suggestions in the
development of the genetic algorithm used to deduce the bound.

## References

  * [1] A. Vera–López, J. Sangroniz, The finite groups with thirteen and fourteen conjugacy classes., Math. Nachr. 280 (2007) 676–694.
  * [2] A. Maróti, A lower bound for the number of conjugacy classes of a finite group, Advances in Mathematics 290 (2016) 1062–1078.
  * [3] L. Pyber, Finite groups have many conjugacy classes, J. London Math. Soc. 46 (1992) 239–249.
  * [4] T. M. Keller, Finite groups have even more conjugacy classes, Isr. J. Math. 181 (2011) 433–444.
  * [5] A. Maróti, On elementary lower bounds for the partition function, Integers 3 (2003) Paper A10, 9 p., electronic only.
  * [6] J. Dénes, P. Erdős, P. Turán, On some statistical properties of the alternating group of degree $n$, Enseignement Math. 15 (1969) 89–99.
  * [7] The GAP Group, GAP – Groups, Algorithms, and Programming, Version 4.11.1 (2021).
URL https://www.gap-system.org

  * [8] The Sage Developers, SageMath, the Sage Mathematics Software System (Version 9.6) (2022).
URL https://www.sagemath.org

  * [9] A. Davies, P. Veličković, L. Buesing, S. Blackwell, D. Zheng, N. Tomašev, R. Tanburn, P. Battaglia, C. Blundell, A. Juhász, M. Lackenby, G. Williamson, D. Hassabis, P. Kohli, Advancing mathematics by guiding human intuition with ai, Nature 600 (7887) (2021) 70–74.

##  Appendix A Searching for groups through genetic algorithms

We briefly outline in this section the rationale behind our algorithmic work
which helped in the deduction of the main result of this short–paper. First
and foremost, we designed a score function associated with the parameters
which we are interested in. In particular, let $\mathcal{G}$ be a collection
of finite permutation groups. We defined a score–function
$s\colon\mathcal{G}\to\mathbb{R}^{+}$ which maximises
$\frac{\log_{2}|G|}{k(G)}$ for $G\in\mathcal{G}$, ie. the highest scoring
groups are those which minimise $\frac{k(G)}{\log_{2}|G|}$, since we are after
Conjecture 1.1.

We then generate a random population of permutation groups, each one
constructed from a random set of permutations and obtaining the permutation
group they generate using the `PermutationGroup` functionality of `SageMath`.
Optionally, one could also enrich the population with ‘known’ groups such as
cyclic groups, alternating groups, etc.

This population is then evolved through our genetic algorithm over a number of
pre–determined iterations, keeping the highest scoring groups from one
iteration to the next.

To evolve the population, besides scoring and selecting the groups in the
population, we include mechanisms to introduce new groups from the highest
scoring ones into the mix, by a number of mutation and cross–over operations.

###  A.1 Cross–over operations

We cross–over between two groups in one of two methods:

  1. i.

either by introducing the direct–product of the two permutation groups into
the population,

  2. ii.

or by considering a generating set of each of the two permutation groups,
randomly choosing a subset of permutations from each, and adding to the
population the group generated by the union of these resulting subsets.

In this manner we aim at being able to search for groups of larger order than
in the current population (through the direct–product) and which inherit
traits from two different high–scoring groups (by mixing their generators).

###  A.2 Mutation operations

The mutation operations considered here take a group and return a new one, by
mutating the generating set of the existing group and then returning whatever
group is generated by the mutated generating set. We outline here a number of
mutations which we explored:

  1. i.

a mutation could be the random removal of a generator from the existing
generating set,

  2. ii.

or the addition of a new randomly generated permutation to the existing
generating set,

  3. iii.

or the addition of a new generator of the form $ghg^{-1}h^{n}$ for some
existing generators $g$ and $h$,

  4. iv.

or the addition of a new generator of the form $(gh)^{n}$ for some existing
generators $g$ and $h$.

In this manner we aim at enriching the population which both ‘random’ as well
as well–understood structure.

###  A.3 Conclusion

In this manner, starting from different randomly generated populations, we
consistently converged towards the group $A_{7}$ as being the highest scoring
one, with other alternating groups scoring similarly (but slightly below),
thereby suggesting the result proved in this paper once the link was
established.

Further searches with a larger population and increased iterations eventually
yielded the Mathieu group of degree 11, $M_{11}$, as a higher scoring group
than $A_{7}$, as noted earlier.