[orcid=0000-0002-6444-2102] [1]

[orcid=0000-0001-7125-1710]

[orcid=0000-0002-4534-3483]

[orcid=0000-0001-5532-6604]

[cor1]Corresponding author. Email address<EMAIL_ADDRESS>(P.M.S.
Sánchez)

# Robust Federated Learning for execution time-based device model
identification under label-flipping attack

Pedro Miguel Sánchez Sánchez  Alberto Huertas Celdrán  José Rafael Buendía
Rubio  Gérôme Bovet  Gregorio Martínez Pérez  Department of Information and
Communications Engineering, University of Murcia, Murcia 30100, Spain
Communication Systems Group (CSG), Department of Informatics (IfI), University
of Zurich UZH, 8050 Zürich, Switzerland  Cyber-Defence Campus within
armasuisse Science & Technology, CH—3602 Thun, Switzerland

###### Abstract

The computing device deployment explosion experienced in recent years,
motivated by the advances of technologies such as Internet-of-Things (IoT) and
5G, has led to a global scenario with increasing cybersecurity risks and
threats. Among them, device spoofing and impersonation cyberattacks stand out
due to their impact and, usually, low complexity required to be launched. To
solve this issue, several solutions have emerged to identify device models and
types based on the combination of behavioral fingerprinting and Machine/Deep
Learning (ML/DL) techniques. However, these solutions are not appropriated for
scenarios where data privacy and protection is a must, as they require data
centralization for processing. In this context, newer approaches such as
Federated Learning (FL) have not been fully explored yet, especially when
malicious clients are present in the scenario setup. The present work analyzes
and compares the device model identification performance of a centralized DL
model with an FL one while using execution time-based events. For experimental
purposes, a dataset containing execution-time features of 55 Raspberry Pis
belonging to four different models has been collected and published. Using
this dataset, the proposed solution achieved 0.9999 accuracy in both setups,
centralized and federated, showing no performance decrease while preserving
data privacy. Later, the impact of a label-flipping attack during the
federated model training is evaluated, using several aggregation mechanisms as
countermeasure. Zeno and coordinate-wise median aggregation show the best
performance, although their performance greatly degrades when the percentage
of fully malicious clients (all training samples poisoned) grows over 50%.

###### keywords:

Model Identification Federated Learning Adversarial Attacks Secure Aggregation
Insider Threat

##  1 Introduction

Currently, there exist a vast number of devices deployed all over the world,
from smart cars, traffic lights, security systems, to smart homes and
industries. The IoT market has grown to a total of 31 billion connected
devices by 2020, with a forecast of $\approx$30 billion devices connected to
each other by 2023, according to Cisco [7]. One of the main reasons of this
growth is the fourth industrial revolution or Industry 4.0, with the explosion
of a set of technologies and paradigms such as 5G, machine and deep learning
(ML/DL), robotics, and cloud computing.

The emergence of such technologies poses new challenges to be solved in order
to ensure a safe and efficient environment [17]. In this sense, there are
billions of connected devices, many of them performing critical tasks where
failures can be fatal, such as autonomous car driving or industrial
operations. In addition, the growing popularity of these technologies makes
them a desirable target for cybercriminals. Between the possible security
threats affecting resource constrained devices, device impersonation is one of
the most serious problems of large organizations with proprietary hardware
where one device model could be impersonated for malicious purposes, such as
industrial espionage. In addition, there are a multitude of counterfeit
devices on the market, some of which are difficult to differentiate from the
original [15].

To solve these issues, device model and type identification based on
performance fingerprinting arises as a solution [20]. The main benefit of
device model identification is to prevent third-party attacks such as
spoofing, as well as to identify malicious or counterfeit devices. Although
there are numerous works in the literature exploring the identification of
models from different performance characteristics, such as execution-time,
network connections or system logs, and leveraging ML/DL for data processing,
these solutions mostly require data centralization, making them not suitable
for scenarios where data leakage protection and privacy is critical. In this
sense, Federated Learning (FL) based techniques have recently gained enormous
prominence [24]. In FL approaches, the training data of the ML/DL models
remain private and while the locally trained models are shared. Later, these
models are aggregated (usually by a central party) into a joint model that
goes back to the clients for further training, repeating the process in a
cyclic fashion. This approach improves both the privacy of the data, as it
does not leave the client but also the communication overhead, as sharing only
model parameters is usually less resource consuming that sharing the complete
data used for training.

In addition, there are few datasets modeling the performance of IoT devices
for identification [20], and any of them is focused on execution time
performance or FL-based scenarios. Moreover, most of the current solutions in
the literature do not explore the impact of possible adversarial attacks
targeting the ML/DL models during their generation and deployment [16]. These
attacks may happen when one of the clients participating in the federation
acts in a malicious way sending corrupted model updates. These problems have
additional importance in FL setups, where the control of the clients is no
longer under the entity generating the joint ML/DL model.

Therefore, this work explores the following three main areas to improve the
completeness of the literature: (i) the identification of device models using
centralized Machine Learning (ML) algorithms and execution time data, (ii) the
decentralization of this training using the Federated Learning (FL)
techniques, and (iii) the use of the Adversarial Machine Learning (AML)
techniques to evaluate and improve the robustness of the generated models. In
this sense, its main contributions are:

  * •

An execution time-based performance dataset collected in 55 different
Raspberry Pi (RPi) devices from 4 different models, and intended for model
identification. This dataset is generated using physical devices under normal
functioning, reflecting a real scenario where many devices are operating.

  * •

The comparison between a centralized and a federated Multi-Layer Perceptron
(MLP) model with identical configuration, only changing its training approach.
It is showed how the federated setup maintains an almost identical model
identification accuracy of 0.9999, without losing performance and improving
data leakage protection and privacy.

  * •

The comparison of different aggregation methods as countermeasure for the
federated model under a label-flipping attack. Federated averaging,
coordinate-wise median, Krum and Zeno aggregation methods are compared,
showing median and Zeno the best results regarding attack resilience.

The remainder of this paper is structured as follows. Section 2 describes the
closest works in the literature, motivating this research. Section 3 explains
the procedure followed to extract the model identification data. Later,
Section 4 compares the performance of a DL-based classifier when it is trained
from a centralized and from a federated approach. Section 5 explains the
adversarial setup followed to test the solution resilience against attacks.
Finally, Section 6 draws the conclusions extracted from the present research
and future lines to explore.

##  2 Related Work

This section will review how the problem of device identification has been
addressed to date from different approaches and techniques. Likewise, some
works in the literature on Federated Learning and Adversarial Machine Learning
will be analyzed.

Device type and model identification has been widely explored in the
literature, with varied data sources and ML/DL-based processing techniques
[20]. As one of the closest works to the present one, the authors of [4]
proposed a novel challenge-response fingerprinting framework called STOP-AND-
FRISK (S&F) to identify classes of Cyber-Physical Systems (CPS) devices and
complement traditional CPS security mechanisms based on hardware and
OS/kernel. It is exposed that unauthorized and spoofed devices may include
manipulated pieces of software or hardware components that may adversely
affect CPS operations or collect vital CPS metrics from the network. Another
interesting paper showing a fingerprinting technique using hardware
performance is [19]. Such a technique is based on the execution times of
instruction sequences available in API functions. Due to its simplicity, this
method can also be performed remotely. Additionally, network is the main data
source employed in the literature for device model and type identification
[13], as it can be collected from an external gateway.

Regarding the application of FL in device identification, the authors of [8]
leveraged FL for device type identification using network-based features.
Here, authors experienced a slightly reduced performance compared to a
centralized setup, 0.851 F1-score in the centralized and 0.849 in the
federated, but the training process was faster and safer. Additionally, in
[14], the authors performed application type classification based on network
traffic using FL to build the models. Although the authors of [21] proposed a
distributed solution for network-based model identification, data is shared
with an aggregator that performs clustering for model inference. Therefore, no
privacy in preserved in this solution.

Moreover, datasets available in the literature for device type or model
identification are focused in dimensions such as network connection [3] or
radio frequency fingerprinting [1]. However, there are not execution time-
based datasets modeling device performance for identification, just some
benchmark datasets focused in other tasks [22].

Concerning adversarial ML in FL, the authors of [10] exposed the impossibility
of the central server to control the clients of the federated network. A
malicious client could send poisoned model updates to the server in order to
worsen learning performance. A new framework for federated learning is
proposed in which the central server learns to detect and remove malicious
model updates using a detection model. Finally, the authors of [18] considered
the presence of adversaries in their solution for FL-based network attack
detection. However, no model identification experiments were carried out.

In conclusion, although each research topic, namely hardware based model
identification, federated learning, and adversarial ML, has been separately
explored. To the best of our knowledge, and as Table 1 shows, there is no work
in the literature analyzing device model identification from a federated
learning perspective. Besides, there is not a dataset focused on model
identification based on execution time-based features. Furthermore, there is
no solution evaluating the impact of adversarial attacks when some clients are
malicious, together with the main aggregation-based attack mitigation
techniques.

Table 1: Comparison of the most relevant model identification literature works. Work |  Model Identification |  FL |  Adv. attack |  Conclusions
---|---|---|---|---
[4] |  ✓ (Hardware-based) |  ✗ |  ✗ |  0.9873 average accuracy using correlation-based algorithms to recognize 11 device classes.
[19] |  ✓ (Hardware-based) |  ✗ |  ✗ |  +200 computers individually identified based on execution-time statistical comparison.
[8] |  ✓ (Network-based) |  ✓ |  ✗ |  0.882 accuracy using a federated LSTM network to identify 10 IoT device types.
[14] |  – (App identification) |  ✓ |  ✗ |  0.92 accuracy using a federated CNN to identify user-level applications.
[21] |  ✓ (Network-based) |  ✗ (Distributed) |  ✗ |  $\approx$0.97 accuracy for clustering-based IoT device type classification.
This work |  ✓ (Hardware-based) |  ✓ |  ✓ (Label-flipping) |  0.9999 accuracy identifying RPi models and adversarial impact analysis.

##  3 Scenario and Dataset Creation

This section describes the scenario and the procedure followed to generate the
execution time dataset used in the present work. Besides, it provides some
insights on the data distribution that can be useful to understand the model
identification performance.

###  3.1 Scenario description

In total, a setup of 55 Raspberry Pis from different models but identical
software images are employed for data collection, running using Raspbian 10
(buster) 32 bits as OS and Linux kernel 5.4.83. The generated dataset is
composed of 2.750.000 vectors (55 devices * 50000 vectors per device). Each
vector has two labels associated, one regaring the individual device that
generated it, and another regarding the model of this device. Data collection
was performed under normal device functioning and default frequency and power
configuration, where the CPU frequency is automatically adjusted according to
the workload. The list of devices contained in the dataset is showed in Table
2.

Table 2: Devices employed in data collection. N devices |  Model |  N samples
---|---|---
12 |  RPi 4 Model B |  660000
22 |  RPi 3 Model B+ |  1210000
5 |  RPi 2 Model B |  275000
16 |  RPi 1 Model B |  880000

###  3.2 Dataset creation

The generated dataset has been made publicly available [6] for download and
research of other authors. The published data includes both identifiers for
RPi model and for individual devices, so new research could be done regarding
individual device identification.

For the device performance dataset generation, the CPU performance of the
device was leveraged as data source. In this sense, the time to execute a
software-based random number generation function was measured in microseconds.

To minimize the impact of noise and other processes running in the device, the
monitored function was executed in groups of 1000 runs a total number of 50000
times per group. Then, for each 1000-run group, a set of statistical features
was calculated, generating a performance fingerprint composed of 50000 vectors
per device. In total, 13 statistical features are calculated: maximum,
minimum, mean, median, standard deviation, mode sum, minimum decrease, maximum
decrease, decrease summation, minimum increase, maximum increase and increase
summation. Decrease and increase values are calculated as the negative or
positive difference between two consecutive values in each 1000-run group.
Besides, the device model is added as label. Table 3 shows an example of a
vector in the dataset belonging to a Raspberry Pi 4 device.

Table 3: Dataset vector example for a RPi4 device. Values represent the time required to execute a function expressed in microseconds. Min |  Max |  Mean |  Median |  Std Dev |  Mode |  Sum |  Min Decr. |  Max Decr. |  Decr. sum |  Min Incr. |  Max Incr. |  Incr. sum |  Model
---|---|---|---|---|---|---|---|---|---|---|---|---|---
2.1 |  12.7 |  4.2 |  4.5 |  1.3 |  3.5 |  4221.8 |  -113.74 |  -8.7 |  -0.001 |  117.8 |  0.001 |  7.81 |  RPi4

###  3.3 Data exploration

Figure 1 shows the data distribution for min, max, mean and median features.
It can be observed how the values vary according to the model that generated
the vector, resulting in a presumably good model identification performance.

Figure 1: Min, Max, Mean and Median feature distributions

##  4 Centralized vs Federated Model Identification Performance

This sections seeks to evaluate firstly the performance of the generated
dataset when identifying the different device models in a DL-based centralized
setup, and secondly the performance variation when the model is generated in a
distributed manner, following a FL-based approach.

###  4.1 Centralized setup

For the centralized experiment, the dataset described in Section 3 is divided
in 80% for training/validation and 20% for testing, without data suffling.
Min-max normalization is applied then using the training data to set the
boundaries.

$x^{\prime}=\frac{x-min(x)}{max(x)-min(x)}$ (1)

To measure the centralized classification performance, a (MLP) classifier is
implemented. After several iterations testing different number of layers and
neurons per layer, the chosen MLP architecture is composed of 13 neurons in
the input layer (one per feature), two hidden layers with 100 neurons each one
using relu (Rectified Linear Unit) as activation function [2], and 4 neurons
in the output softmax layer (one per model class). Adam [9] was used as
optimizer with a 0.001 learning rate, and 0.9 and 0.999 as first and second-
order moments. Table 4 shows the details of the model.

Table 4: MLP architecture for model identification. Layer |  Neurons |  Activation
---|---|---
1 |  13 |  -
2 |  100 |  relu
3 |  100 |  relu
4 |  4 |  softmax
Optimizer |  Adam |

With this setup, the MLP is trained for 100 epochs using early stopping if no
validation accuracy improvement occurs in 20 epochs.

Figure 2 shows the confusion matrix resultant of the evaluation of the test
datastet. As it can be seen, almost a perfect identification is achieved, with
only 15 samples being misclassified out of $\approx$550000 (0.999972
accuracy). These results are aligned with the expectations, as having
different CPUs in each RPi model makes the execution time of the same
functions different between them. However, model identification performance is
not the main focus of the present work, where the priority is to prove the
effectiveness of a federated setup and the impact of adversarial attacks and
countermeasures.

Figure 2: Centralized evaluation confusion matrix.

###  4.2 Federated scenario and results

Once the centralized model has been obtained, the decentralized model is
implemented using FL to compare the performance of both approaches. The FL
approach is based in horizontal FL, where the clients have datasets with the
same features but from different data samples.

For implementation, the IBM Federated Learning library [11] is used, which
incorporates the necessary tools to perform the training in a decentralized
manner.

####  4.2.1 Scenario

For the decentralization of the training phase, a scenario has been created in
which there are 5 independent organizations in which the available data are
distributed. Each of them has a certain number of devices belonging to
different models, but not all of them have information on all models, i.e.
there are organizations that only have devices of type 4 model, others that
only have devices of type 2 and 3 models, etc. Figure 3 provides the details
of the device distribution in each organization. Therefore, the 5
organizations intend to generate a global model capable of identifying all the
existing device models among all of them. This setup leads to an scenario of
Non-IID (Non-Independent and Identically Distributed) data, harder to solve
with FL as model aggregation will be negatively influenced in the aggregated
models are very different to each other.

Figure 3: Data division in the federated learning scenario.

####  4.2.2 Federated architecture design

In order to test the performance of a FL-based setup, first it is necessary to
define the architecture to be implemented. In this sense, Figure 4 shows the
organization of the different client which will hold the data and upload their
local models to the aggregator in order to cyclically build a common model
capable of making predictions based on the local data of all clients.

Figure 4: Designed federated architecture for experimentation.

####  4.2.3 Performance Evaluation

In order to fairly compare the models, the MLP architecture to be trained will
be the same as the one used in the centralized model (see Table 4, i.e. the
layers will have 13, 100, 100, 100, and 4 neurons, from input to output. As
aggregation method, Federated Averaging is applied as proposed in [12]. As
initialization step, the aggregation server performs two tasks: (1) to
initialize the weights of the model that the clients will start to train, so
all clients start from the same setup; (2) to retrieve for each client its
min-max values of each feature for common dataset normalization, having a min-
max normalization for each dataset x in organization o $\in x$ defined as:

$x^{\prime}_{o}=\frac{x_{o}-min(x_{i\in[n]})}{max(x_{i\in[n]})-min(x_{i\in[n]})}$
(2)

Algorithm 1 defines the iterative training process for the model generation,
assuming previous dataset normalization. Each client performs local updates of
the model and returns them to the server for aggregation, repeating then the
process for the desired number of rounds.

0: Server executes:

0: initialize $w_{0}$

0: for each round t = 1,2… do

0: $m\leftarrow max(C\cdot K,1)$

0: $S_{t}\leftarrow$ (random set of $m$ clients)

0: for each client $k\in S_{t}$ in parallel do

0: $w^{k}_{t+1}\leftarrow$ ClientUpdate($k,w_{t}$)

0: $w_{t+1}\leftarrow\sum_{k=1}^{K}\frac{n_{k}}{n}w_{t+1}^{k}$ //Aggregation

0: ClientUpdate($k,w$): //Run on client $k$

0: $B\leftarrow$ (split $P_{k}$ into batches of size $B$)

0: for each local epoch i from 1 to W do

0: for batch b $\in B$ do

0: $w\leftarrow w-\eta\bigtriangledown l(w;b)$ //Local update

0: return $w$ to server

Algorithm 1 FederatedAveraging. The K clients are indexed by k; B is the local
minibatch size, E is the number of local epochs, and $\eta$ is the learning
rate; $w$ are the model weights; $P_{k}$ is the local dataset of client $k$.
[12]

The training process was executed for 90 federated rounds, with one epoch per
round. Figure 5 shows the evolution of the local validation accuracy for each
one of the clients during the training process. It can be appreciated how the
maximum performance is reached around epoch 50, and then the accuracy scores
for each client keep oscillating between 0.95 and 1 until round 90.

Figure 5: Federated training validation accuracy evolution.

Regarding performance, Figure 6 shows the results of the test dataset
evaluation, the same dataset than in the centralized setup. Here, the results
are almost identical, with only 17 errors in $\approx$550000 test samples and
an accuracy of 0.999969.

Figure 6: Federated evaluation confusion matrix.

From the previous results, a main conclusion can be extracted: no performance
loss has been introduced in the resultant model due to the application of a
FL-based approach. Besides, as no data has left each organization in the
process, the privacy of the information has been kept private successfully.

##  5 Adversarial Attack and Robust Aggregation

After testing the effectiveness of Federated Learning, its robustness will be
tested using adversarial attacks, specifically the label-flipping technique,
using different aggregation algorithms in order to see which one best fits the
proposed scenario in the presence of attacks.

###  5.1 Label flipping attack

The label-flipping adversarial technique is applied during the training
process, using the same scenario described above with the difference that this
time part of the data will be poisoned.

In this sense, the federated training is carried out by poisoning 25, 50, 75
and 100% of the data of 1, 2 and 3 different organizations, representing 20%,
40% and 60% malicious clients, respectively. These configurations are used
because potential malicious clients may not poison all their data and just one
portion, in order to go undetected and make their activity more difficult to
identify. So, a total of 12 adversarial scenarios have been created (4 poising
percentages * 3 possible malicious organizations). This setup is generated by
modifying the labels of the training data, changing the value of each label to
a random value between 1 and 4 that is not the value of the original label.
The poisoned organizations are ORG1, ORG2 and ORG4 (in that order for 1, 2, 3
malicious clients).

Figure 7 shows the results when FedAvg is applied as aggregation algorithm in
the 12 previous adversarial scenarios (as well as when no label-flipping
attack is applied).

Figure 7: FedAvg performance against label-flipping attack. (X axis depicts
poisoning percentage, Y axis depicts accuracy).

As can be seen, aggregation by averaging offers good performance up to a 50%
poisoning, maintaining the accuracy over 0.9. However, accuracy drops rapidly
to hit rates close to 0% when the poisoning is 75% or higher. Therefore,
FedAvg cannot be considered a robust aggregation method in the presence of the
label-flipping attack. Next, 3 different aggregation methods will be analyzed
in the following in order to check which one offers better performance.

###  5.2 Robust aggregation methods

Next, several aggregation methods focused on improving the model resilience to
malicious clients will be evaluated and compared to the default FedAvg
algorithm.

####  5.2.1 Coordinate-wise median aggregation

Coordinate-wise median [25] follows the scheme of the aggregation by average
with the difference that the combination of the weights is done by calculating
the median of each weight of the local models. In short, following Algorithm
1, the averaging aggregation step is substituted by a median operation.

Next, Figure 8 depicts the accuracy results when the different attack setups
are applied when using median aggregation. Coordinate-wise median follows a
similar pattern to FedAvg aggregation, dropping from 50% poisoning rate.
However, it has performed better especially when there is only one poisoned
organization (20% malicious clients). While FedAvg dropped to 0.20-0.40, the
median has remained around $\approx$0.9.

Figure 8: Coordinate-wise median aggregation performance against label-
flipping attack. (X axis depicts poisoning percentage, Y axis depicts
accuracy).

####  5.2.2 Krum aggregation

The idea behind Krum [5] is to select one of the m local models that is most
similar to the rest as the global model. The idea is that even if the selected
model is a poisoned model the impact would not be so great since it would be
similar to other models that are probably not poisoned. The aggregator
calculates the sum of the distances between each model and its closest local
models. Krum selects the local model with the smallest sum of distances as the
global model. Figure 9 shows the results when Krum is applied as aggregation
algorithm.

Figure 9: Krum performance against label-flipping attack. (X axis depicts
poisoning percentage, Y axis depicts accuracy).

As can be seen, Krum has remained constant for all configurations with an
accuracy of 0.6896 . This is because this aggregation method chooses a single
local model as the global model and discards the information from the rest of
the local models. Therefore, what is happening is that it always chooses the
same local model, and this one belongs to an organization that has not been
poisoned, so the hit rate remains constant. Figure 10 shows that the resulting
global model only recognizes device models of types 0 and 2.

On the other hand, this organization has not been poisoned, which explains
that the performance remains constant since the resulting global model is
identical regardless of the percentage of poisoning. Therefore, it can be
concluded that Krum is selecting the resulting local model of organization 3
in all scenarios, loosing the information regarding the classes not seen in
this organization (see Figure 3).

Figure 10: Krum confusion matrix during label-flipping attack.

####  5.2.3 Zeno aggregation

Zeno [23] is suspicious of potentially malicious organizations and uses a
ranking-based preference mechanism. The number of malicious organizations can
be arbitrarily large, and only the assumption that ’clean’ organizations exist
(at least one) is used. Each organization is ranked based on the estimated
descent of the loss function. The algorithm then aggregates the organizations
with the highest scores. The score roughly indicates the reliability of each
organization. In this sense, it could be seen as a combination of Krum and
averaging aggregation mechanisms. Figure 11 shows the results when Zeno is
applied as aggregation algorithm.

Figure 11: Zeno performance against label-flipping attack. (X axis depicts
poisoning percentage, Y axis depicts accuracy).

In this case, Zeno has outperformed the aggregation by median and Krum when
only one client is malicious (20% of the total), achieving 0.9072 accuracy.
When there is only one poisoned organization Zeno remains constant without
being altered by this attack. Figure 12 shows the confusion matrix of Zeno
when one client is malicious. It can be appreciated how the performance
decrease comes from the impossibility of classifying the second class, the one
under represented in the scenario as there are only 5 RPi2 in the dataset.

Figure 12: Zeno confusion matrix with one malicious client.

When there are 2 or 3 poisoned organizations, Zeno performance drops once the
poisoning rate reaches 75% and 100%, but it still manages to maintain an
acceptable performance above 0.50, considering the degree of the attack.
Figure 13 compares the performance evolution of Zeno and coordinate-wise
median with different number of poisoned organizations. As it can be
appreciated, median performance is higher in all scenarios until the poisoning
percentage goes above 50%. After that, Zeno shows a better or equal
performance in all cases, being the greater difference when three
organizations are completely malicious (60% malicious clients).

(a) One organization poisoned.

(b) Two organizations poisoned.

(c) Three organizations poisoned.

Figure 13: Zeno and coordinate-wise median aggregation with different number
of poisoned clients (X axis depicts poisoning percentage, Y axis depicts
accuracy).

##  6 Conclusions and Future Work

In the present work it has been demonstrated that it is possible to identify
device models using only statistical data concerning the CPU execution time of
the device. An MLP model has been obtained capable of identifying four RPi
device models with a 99.99% accuracy rate. Besides, the effectiveness of
Federated Learning technique has been tested against centralized learning. For
this setup, a scenario has been proposed where a total of 5 organizations aim
to create a model capable of identifying the device models without sharing the
actual data with each other. The resulting model has obtained identical
performance in both cases, centralized and distributed. Thus taking advantage
of the benefits offered by Federated Learning, training a data privacy and
data security preserving model, while maintaining the efficiency of the model
obtained through a traditional approach. On the other hand, different
aggregation algorithms have been tested in order to check which one best fits
the proposed scenario facing a label-flipping attack. Zeno has turned out to
be the best performing aggregation method in the presence of attacks due to
combining the Krum and mean aggregation methods. By selecting the m best
models and aggregating them using mean aggregation, less information is lost
than with Krum by ignoring certain organizations that are considered
malicious. Finally, the data collected for the previous experimentation has
been made publicly available due to the lack of performance fingerprinting
datasets focused on device identification, and prepared for FL-based setups.

As future work, the efforts will be focused on experimentation with more types
of device models with more complex scenarios such as making each device a
single client instead of being grouped into organizations. On device
identification, it is planned to focus on identifying individual devices with
a high hit rate and not just identifying device models, as well as testing
other modes of identification by collecting data from other hardware elements
than the CPU. It would also be interesting to poison the local model weights
instead of the local data (model poisoning) or experiment with other
adversarial attack techniques such as Evasion attacks, where the goal is to
trick the model once it is trained and not to poison the training process.

## Data Availability

The datasets generated during and/or analysed during the current study are
available in the RPi model device identification Mendeley Data repository,
https://doi.org/10.17632/vr9wztmfxg.2

## Acknowledgment

This work has been partially supported by (a) the Swiss Federal Office for
Defense Procurement (armasuisse) with the TREASURE (R-3210/047-31) and
CyberSpec (CYD-C-2020003) projects and (b) the University of Zürich UZH.

## References

  * Abbas et al. [2021] Abbas, S., Nasir, Q., Nouichi, D., Abdelsalam, M., Talib, M.A., Waraga, O.A., et al., 2021\.  Improving security of the internet of things via rf fingerprinting based device identification system.  Neural Computing and Applications , 1–17.
  * Agarap [2018] Agarap, A.F., 2018.  Deep learning using rectified linear units (relu).  arXiv preprint arXiv:1803.08375 .
  * Aksoy and Gunes [2019] Aksoy, A., Gunes, M.H., 2019\.  Automated iot device identification using network traffic, in: ICC 2019-2019 IEEE International Conference on Communications (ICC), IEEE. pp. 1–7.
  * Babun et al. [2021] Babun, L., Aksu, H., Uluagac, A.S., 2021.  Cps device-class identification via behavioral fingerprinting: From theory to practice.  IEEE Transactions on Information Forensics and Security 16, 2413–2428.
  * Blanchard et al. [2017] Blanchard, P., El Mhamdi, E.M., Guerraoui, R., Stainer, J., 2017\.  Machine learning with adversaries: Byzantine tolerant gradient descent, in: Proceedings of the 31st International Conference on Neural Information Processing Systems, pp. 118–128.
  * Bovet and Sánchez Sánchez [2021] Bovet, G., Sánchez Sánchez, P.M., 2021\.  Rpi model device identification.  URL: http://dx.doi.org/10.17632/vr9wztmfxg.2, doi:10.17632/vr9wztmfxg.2.
  * Cisco [2020] Cisco, 2020.  Cisco annual internet report (2018–2023) white paper.  https://www.cisco.com/c/en/us/solutions/collateral/executive-perspectives/annual-internet-report/white-paper-c11-741490.html.  [Online; accessed 22-November-2021].
  * He et al. [2021] He, Z., Yin, J., Wang, Y., Gui, G., Adebisi, B., Ohtsuki, T., Gacanin, H., Sari, H., 2021.  Edge device identification based on federated learning and network traffic feature engineering.  IEEE Transactions on Cognitive Communications and Networking .
  * Kingma and Ba [2014] Kingma, D.P., Ba, J., 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * Li et al. [2020] Li, S., Cheng, Y., Wang, W., Liu, Y., Chen, T., 2020\.  Learning to detect malicious clients for robust federated learning.  arXiv preprint arXiv:2002.00211 .
  * Ludwig et al. [2020] Ludwig, H., Baracaldo, N., Thomas, G., Zhou, Y., Anwar, A., Rajamoni, S., Ong, Y., Radhakrishnan, J., Verma, A., Sinn, M., et al., 2020\.  Ibm federated learning: an enterprise framework white paper v0. 1.  arXiv preprint arXiv:2007.10987 .
  * McMahan et al. [2017] McMahan, H.B., Moore, E., Ramage, D., Hampson, S., y Arcas, B.A., 2017.  Communication-efficient learning of deep networks from decentralized data.  arXiv:1602.05629.
  * Meidan et al. [2017] Meidan, Y., Bohadana, M., Shabtai, A., Guarnizo, J.D., Ochoa, M., Tippenhauer, N.O., Elovici, Y., 2017.  Profiliot: A machine learning approach for iot device identification based on network traffic analysis, in: Proceedings of the symposium on applied computing, pp. 506--509.
  * Mun and Lee [2021] Mun, H., Lee, Y., 2021.  Internet traffic classification with federated learning.  Electronics 10, 27\.
  * Negka et al. [2019] Negka, L., Gketsios, G., Anagnostopoulos, N.A., Spathoulas, G., Kakarountas, A., Katzenbeisser, S., 2019\.  Employing blockchain and physical unclonable functions for counterfeit iot devices detection, in: Proceedings of the International Conference on Omni-Layer Intelligent Systems, pp. 172--178.
  * Nguyen et al. [2019] Nguyen, T.D., Marchal, S., Miettinen, M., Fereidooni, H., Asokan, N., Sadeghi, A., 2019\.  DÏot: A federated self-learning anomaly detection system for IoT, in: 39th IEEE International Conference on Distributed Computing Systems, pp. 756--767.  doi:10.1109/ICDCS.2019.00080.
  * Nižetić et al. [2020] Nižetić, S., Šolić, P., González-de, D.L.d.I., Patrono, L., et al., 2020.  Internet of things (iot): Opportunities, issues and challenges towards a smart and sustainable future.  Journal of Cleaner Production 274, 122877.
  * Rey et al. [2021] Rey, V., Sánchez, P.M.S., Celdrán, A.H., Bovet, G., Jaggi, M., 2021.  Federated learning for malware detection in iot devices.  arXiv preprint arXiv:2104.09994 .
  * Sanchez-Rola et al. [2018] Sanchez-Rola, I., Santos, I., Balzarotti, D., 2018.  Clock around the clock: Time-based device fingerprinting, in: 2018 ACM SIGSAC Conference on Computer and Communications Security, pp. 1502--1514.
  * Sánchez Sánchez et al. [2021] Sánchez Sánchez, P.M., Jorquera Valero, J.M., Huertas Celdrán, A., Bovet, G., Gil Pérez, M., Martínez Pérez, G., 2021\.  A survey on device behavior fingerprinting: Data sources, techniques, application scenarios, and datasets.  IEEE Communications Surveys Tutorials 23, 1048--1077.  doi:10.1109/COMST.2021.3064259.
  * Thangavelu et al. [2019] Thangavelu, V., Divakaran, D.M., Sairam, R., Bhunia, S.S., Gurusamy, M., 2019.  DEFT: A distributed IoT fingerprinting technique.  IEEE Internet of Things Journal 6, 940--952.
  * Varghese et al. [2021] Varghese, B., Wang, N., Bermbach, D., Hong, C.H., Lara, E.D., Shi, W., Stewart, C., 2021.  A survey on edge performance benchmarking.  ACM Comput. Surv. 54.  URL: https://doi.org/10.1145/3444692, doi:10.1145/3444692.
  * Xie et al. [2019] Xie, C., Koyejo, S., Gupta, I., 2019.  Zeno: Distributed stochastic gradient descent with suspicion-based fault-tolerance, in: International Conference on Machine Learning, PMLR. pp. 6893--6901.
  * Yang et al. [2019] Yang, Q., Liu, Y., Chen, T., Tong, Y., 2019.  Federated machine learning: Concept and applications.  ACM Transactions on Intelligent Systems and Technology 10, 1--19.
  * Yin et al. [2018] Yin, D., Chen, Y., Kannan, R., Bartlett, P., 2018.  Byzantine-robust distributed learning: Towards optimal statistical rates, in: International Conference on Machine Learning, PMLR. pp. 5650--5659.


# Data-driven low-dimensional model of a sedimenting flexible fiber

Andrew J. Fox  Michael D. Graham<EMAIL_ADDRESS>Department of Chemical
and Biological Engineering,
University of Wisconsin-Madison,
Madison, WI 53706, USA

###### Abstract

The dynamics of flexible filaments entrained in flow, important for
understanding many biological and industrial processes, are computationally
expensive to model with full-physics simulations. This work describes a data-
driven technique to create high-fidelity low-dimensional models of flexible
fiber dynamics using machine learning; the technique is applied to
sedimentation in a quiescent, viscous Newtonian fluid, using results from
detailed simulations as the data set. The approach combines an autoencoder
neural network architecture to learn a low-dimensional latent representation
of the filament shape, with a neural ODE that learns the evolution of the
particle in the latent state. The model was designed to model filaments of
varying flexibility, characterized by an elasto-gravitational number
$\mathcal{B}$, and was trained on a data set containing the evolution of
fibers beginning at set angles of inclination. For the range of $\mathcal{B}$
considered here (100-10000), the filament shape dynamics can be represented
with high accuracy with only four degrees of freedom, in contrast to the 93
present in the original bead-spring model used to generate the dynamic
trajectories. We predict the evolution of fibers set at arbitrary angles and
demonstrate that our data-driven model can accurately forecast the evolution
of a fiber at both trained and untrained elasto-gravitational numbers.

††preprint: APS/123-QED

##  I Introduction

The dynamics of flexible filaments in flow are an active area of interest,
with applications in many biological and industrial systems. [1, 2]
Biopolymers, such as DNA, can behave like elastic fibers suspended in flow,
and understanding their dynamics is important for designing and understanding
various medical procedures. [3] In paper manufacturing, the development and
processing of fibers from wood pulp is vital for creating product at given
specifications and for optimizing production to reduce waste. [4] In
particular, the relationship and competition between hydrodynamic forces
acting on the particle and elastic forces acting along the fiber can produce a
range of interesting behavior. [5, 6, 7]

In recent years, modeling the sedimentation dynamics of flexible filaments has
been explored through experiments and numerical analyses. Xu and Nadim [8]
investigated the sedimentation of a flexible slender particle with small
deformation in quiescent viscous fluid. They performed a theoretical analysis
of a slender body undergoing small amplitudes of deformation using the Euler-
Bernoulli beam theory, showing a ”U-shape” conformation at long settling
times. Schlagberger and Netz [9] furthered this analysis by performing
computational simulations of the system. Using numerical computations of a
settling filament of a bead-spring-model, they expanded this analysis to
understand the dynamics at higher amplitudes of deformation. Lagomarsino  _et
al._ [10] continued the use of bead-spring models to investigate the
sedimentation of highly flexible filaments. They analyzed very flexible
filaments with a bead-spring model, revealing a stable ”W-shape” confirmation
for sufficiently high degrees of flexibility. Li  _et al._ [11] used
analytical and numerical analysis to understand the settling dynamics of an
elastic fiber when oriented along its major axis. They applied a slender body
theory to highly flexible filaments, identifying a buckling instability when
sedimenting in this configuration. Shojaei and Dehghani [12] investigated the
rotational dynamics of a flexible filament settling in quiescent flow. Using
numerical simulations of a slender beam model, they identified different
rotational behaviors of filaments at different degrees of elasticity. Delmotte
_et al._ [13] developed a generalized formulation of bead spring models with
which to investigate the dynamics of flexible fibers outside of Stokes flow.
Their approach allowed them to explore the dynamics of elastic filaments at
finite Reynolds numbers. Marchetti  _et al._ [14] performed both experiments
and numerical computations to analyze the settling behavior of flexible
filaments. Using a bead-spring model, they successfully demonstrated a
quantitative agreement between these studies, identifying distinct regimes of
particle deformation. Cunha  _et al._ [15] compared the dynamics of non-
Brownian and Brownian flexible filaments settling in quiescent fluid. Using
experiments and numerical simulations, they demonstrated good agreement
between the observed and computed dynamics for both filament regimes.

As the interplay between the hydrodynamic and elastic forces in these studies
can be complex and the entire particle structure must be adequately captured,
full physics simulations of fluid dynamics can be computationally expensive
[16, 17]. One potential solution to model fluid and particle mechanics without
long computational times is to develop low-dimensional, data-driven models[18,
19, 20]. These models, created though machine learning, can capture the
evolution of a dynamical system with significantly fewer dimensions than full
physics simulations, with minimal loss of prediction accuracy[21, 22, 23].
These models use latent representations of the full-field data to learn a
parameterization of the manifold on which the dynamics of a system exist,
primarily through the use of undercomplete autoencoder neural networks[24, 25,
26]. By combining this with a time-integrating neural network, either in
discrete or continuous time, the full dynamics of the system can be captured
with much less computational expense[27, 28].

Here, we will develop data-driven models to forecast the dynamics of a
flexible fiber settling in a quiescent Newtonian fluid in the Stokes flow
regime. To our knowledge, such modeling approaches have not been applied to
deformable-particle dynamics problems. Specifically, we will model the
trajectory and shape evolution of a filament, beginning lineally-aligned at
rest at an arbitrary angle of initial orientation, and evolving until the
fiber reaches a terminal configuration. We will train autoencoder neural
networks to learn a low-dimensional representation of the shape of a fiber,
then use a time-integrating neural network to learn both the evolution of the
filament shape and the position; these will be combined to produce a single
low-dimensional model of the fiber dynamics. We will demonstrate that the
fiber structure and dynamics can be captured at significantly lower dimension
than the full-field state. We will demonstrate that our model can forecast the
evolution of fibers at arbitrary angles of initial inclination and can predict
the dynamics of fiber shape and position at a range of elasto-gravitational
numbers with minimal error.

##  II Formulation

###  II.1 Physics-based model

Consider a flexible filament of length $L$ entrained in a quiescent, viscous
Newtonian fluid of viscosity $\mu$, driven by a gravitational force,
$\mathbf{F}_{g}$, and experiencing a viscous drag, $\mathbf{F}_{drag}$,
dependent on the instantaneous fiber shape. The fiber begins aligned in the
$xy$-plane at rest, set at an arbitrary initial angle of inclination,
$\theta_{0}$, relative to the $x$-axis. The filament is allowed to settle
under an external force (typically considered to be gravity), whereupon the
elasticity of the fiber induces a deformation of its shape, trending towards a
“U-shape” bending[8]. At long times, the shape of the fiber will approach a
terminal “U-shape”, with the particular shape dependent on the balance between
gravitational and elastic forces, given by an elasto-gravitational number
$\mathcal{B}$[14]. A diagram of the problem is shown in Figure 1.

Figure 1: A flexible filament of length $L$ settling under an external force
$\mathbf{F}_{g}$ in a quiescent Newtonian fluid. The filament begins at an
arbitrary initial inclination, $\theta_{0}$, relative to the $x$-axis, and
evolves until it reaches a terminal shape. The fiber is modeled as a series of
$N$ beads of radius $a$ connected by springs, with the center bead, which has
position $\mathbf{c}(t)$, shown in red.

To model this system, we apply a bead-spring model, wherein the fiber is
discretized into a series of $N=31$ spherical beads of radius $a$, shown by
Marchetti  _et al._ [14] to adequately capture the range of rotational
dynamics of interest in this work. The spheres are connected by massless
springs of equilibrium length, $l_{0}=2a$ with an equilibrium angle formed by
three beads of $\theta=\pi$. The fluid is modeled as a viscous Newtonian fluid
of viscosity $\mu$ governed by the Stokes equations, with fluid and particle
inertia neglected. The evolution of the filament shape and position is
governed by a force balance, with each bead experiencing gravitational and
viscous drag forces, and the beads interacting through multi-body hydrodynamic
interactions and elastic forces. The velocity of bead $\alpha$, given by
$\dot{\mathbf{X}}^{\alpha}$, at position $\mathbf{X}^{\alpha}$ and interacting
with beads $\beta$ at position $\mathbf{X}^{\beta}$ is given by [14]

$\dot{X}_{i}^{\alpha}=\sum_{\beta}\mathcal{M}_{ij}^{\alpha\beta}\left(F^{\beta}_{j}-\frac{\partial\mathcal{U}}{\partial
X_{j}^{\beta}}\right).$ (1)

Here $\mathcal{M}_{ij}^{\alpha\beta}$ is the position-dependent $N$-body
mobility tensor, which accounts for hydrodynamic interactions between spheres
$\alpha$ and $\beta$ [29], and $\mathbf{F}^{\beta}$ is the external force
exerted on each bead by gravity, which in a viscous fluid is balanced by
Stokes drag and is given by $\mathbf{F}^{\beta}=6\pi\mu a\mathbf{U}_{s}$,
where $\mathbf{U}_{s}$ is the Stokes settling velocity. The elastic potential,
$\mathcal{U}$, stems from a discrete form of the wormlike chain model, written
as [9]

$\mathcal{U}=\sum_{\gamma}\left[aS\left(\frac{\mathbf{X}^{\gamma,\gamma+1}}{2a}-1\right)^{2}+\frac{B}{2a}\left(1-\cos\theta^{\gamma,\gamma+1}\right)\right].$
(2)

For an isotropic elastic cylinder, the stretching modulus, $S$, and bending
modulus $B$, depend solely on the cylinder radius and Youngs modulus, $E$,
such that $S=E\pi a^{2}$ and $B=E\pi a^{4}/4$, respectively; this formulation
allows the stretching and bending of the filaments to be governed by a single
physical parameter and has been shown to accurately capture the dynamics of
flexible fibers[14]. The geometric dependence of the elastic potential is
given by the distance between neighboring beads $\gamma$ and $\gamma+1$,
$\mathbf{X}^{\gamma,\gamma+1}=\mathbf{X}^{\gamma+1}-\mathbf{X}^{\gamma}$, and
the angle between neighboring springs $\mathbf{X}^{\gamma,\gamma+1}$ and
$\mathbf{X}^{\gamma,\gamma-1}$, $\theta^{\gamma,\gamma+1}$. The inter-particle
mobility tensor is selected as the Rotne-Prager-Yamakawa tensor, given for
$\alpha\neq\beta$ by [14]

$\mathcal{M}_{ij}^{\alpha\beta}=\frac{1}{6\pi\mu
a}\left\\{\frac{3}{4}\left[\frac{\delta_{ij}}{\frac{X^{\alpha,\beta}}{a}}+\frac{\frac{X_{i}^{\alpha,\beta}X_{j}^{\alpha,\beta}}{a^{2}}}{\left(\frac{X^{\alpha,\beta}}{a}\right)^{3}}\right]\right.+\left.\frac{3}{2}\left[\frac{\delta_{ij}}{3\left(\frac{X^{\alpha,\beta}}{a}\right)^{3}}-\frac{\frac{X_{i}^{\alpha,\beta}X_{j}^{\alpha,\beta}}{a^{2}}}{\left(\frac{X^{\alpha,\beta}}{a}\right)^{5}}\right]\right\\}.$
(3)

The chosen tensor accounts for hydrodynamic interactions between particles up
to order $O\left(\frac{a}{X^{\alpha,b}}\right)^{3}$, where
$X^{\alpha,\beta}=|\mathbf{X}^{\alpha,\beta}|=|\mathbf{X}^{\beta}-\mathbf{X}^{\alpha}|$
is the distance between beads $\alpha$ and $\beta$. A full Rotne-Prager-
Yamakawa tensor, which provides a separate formulation for
$X^{\alpha,\beta}<2a$ and is positive definite for all particle
configurations, has been observed by Marchetti  _et al._ [14] to not
significantly change the filament dynamics in the parameter range of interest
in this study, as the spheres in the bead-spring model do not significantly
overlap. The self-mobility $\mathcal{M}_{ij}^{\alpha\alpha}$ is given by the
Stokes mobility[14]. The bead velocity equation can be nondimensionalized by
the length scale $a$ and the force scale $6\pi\mu aU_{s}$, resulting in

$\dot{\hat{X}}_{i}^{\alpha}=\sum_{\beta}\hat{\mathcal{M}}_{ij}^{\alpha\beta}\left(\hat{F}_{j}^{\beta}-\mathcal{E}\frac{\partial\hat{\mathcal{U}}}{\partial
X_{j}^{\beta}}\right),$ (4)

where $\hat{\mathbf{X}}(t)$ is the dimensionless bead positions and
$\hat{F}_{j}^{\beta}$ is a downward unit vector. Now a dimensionless modulus
$\mathcal{E}=\frac{E\pi a^{2}}{6\pi\mu aU_{S}}$ appears as the only parameter
in the formulation. This can be related to the elasto-gravitational number by
$\mathcal{B}=32\Delta\rho g\ell^{3}/Ea^{2}\equiv 24(a/l)^{-3}/\mathcal{E}$;
here, the volume of the filament is calculated via a bent cylinder
encompassing the bead-spring system, the use of which has been shown to
successfully match the behaviors of flexible fibers in experiments[14].

In this paper, we investigated the settling dynamics of settling fibers with
elasto-gravitational numbers in the range $100\leq\mathcal{B}\leq 10000$. At
each $\mathcal{B}$, we calculated the trajectory of filaments beginning at 23
initial orientations, evenly spaced in the range $0\leq\theta_{0}\leq
11\pi/24$. Integration of the positions of the bead-spring system was
performed using the “dopri” integrator of the ”ode” solver in Python, an
explicit Runge-Kutta method of order (4)5, in a method adapted from Marchetti
_et al._ [14]. The shape and position of each filament was allowed to evolve
until the fiber reached a terminal shape, dependent on $\mathcal{B}$, taking a
total settling time $T=T(\mathcal{B},\theta_{0})$.

The evolution of the shape of filaments of $\mathcal{B}=1000$ at twelve
initial orientations is shown in Figure 2. In each trajectory, the filament
bends into a “U-shape” configuration and rotates into a position such that the
shape is symmetric in the $x$-direction. The terminal shape is identical for
all initial orientations for filaments of a given $\mathcal{B}$.

Figure 2: The evolution of the shape of a filament settling in a quiescent
Newtonian fluid from an initial orientation to a common terminal shape at
$\mathcal{B}=1000$ for all initial angles of orientation within the training
data set.

Shown in Figure 3 is the trajectory of the center bead of a modeled filament,
beginning at twelve different initial angles of inclination. The initial
positions of the filaments are arbitrary, and are shown here set at $y=0$ and
offset in the $x$-direction to differentiate between trajectories. The total
distance traveled by each filament depends on the initial orientation, as the
time required to reach the terminal shape depends on the bending and rotation
required. The settling time and positional translation increases as the
initial angle of orientation increases, as further time is necessary to
reorient into a symmetric shape.

Figure 3: The trajectories $\mathbf{c}(t)$ of the center bead of a filament
settling in a quiescent Newtonian fluid at $\mathcal{B}=1000$ for all initial
angles of inclination within the training data set. The initial positions are
denoted by symbol “$\blacktriangledown$” and the terminal positions are
denoted by the symbol “$\times$”.

In the range of $\mathcal{B}$ investigated in this study, the flexible
filament undergoes three distinct regimes of rotational dynamics. Examples of
these dynamics are shown in Figure 4, with shape evolution shown for
$\theta_{0}=\pi/4$ at seven elasto-gravitational numbers within the scope of
this investigation. At low $\mathcal{B}$, the shape dynamics are primarily
rotational, with the stiff fiber bending only slightly in the terminal shape.
At moderate $\mathcal{B}$, the rotational dynamics exist in a bending regime,
where the left and right sides of the filament bend toward each other to reach
the terminal shape and only slight rotation occurs. For high $\mathcal{B}$,
the filament undergoes a “snaking” behavior, where the lower side of the fiber
bends over itself and slides along the higher side to reach its final shape.
The terminal shape, as shown, strongly and solely depends on the given elasto-
gravitational number, with higher $\mathcal{B}$ producing more pronounced
”U-shapes” which are identical, to an arbitrary degree of precision, within
each $\mathcal{B}$.

Figure 4: The evolution of the shape of a filament settling in a quiescent
Newtonian fluid from a common initial angle of orientation of $\pi/4$ to a
terminal shape for all $\mathcal{B}$ within the training data set.

###  II.2 Data-driven model

Using the data from the bead-spring model, data-driven models were trained to
forecast the evolution of the shape and position of a filament beginning at an
arbitrary initial angle of orientation. A total of 227 trajectories were
generated with $100\leq\mathcal{B}\leq 10000$ and $0\leq\theta\leq 22\pi/48$.
The calculated trajectories were separated into three data sets: the training
data, test data set A, and test data set B. A breakdown of the distribution of
the trajectories into each data set is shown in Table II.2. The training data
(green on Table II.2) consisted of trajectories of filaments with elasto-
gravitational numbers $\mathcal{B}=[100,200,400,1000,2000,4000,100000]$ at
twelve initial angles of inclination, evenly space in the range
$0\leq\theta_{0}\leq 22\pi/48$. Test data set A (yellow) consisted of
trajectories of fibers at the same elasto-gravitational numbers as the
training data set at eleven untrained initial angles of inclination, evenly
spaced in the range $\pi/48\leq\theta_{0}\leq 21\pi/48$. Test data set A
allows for evaluating the ability of the data-driven models to interpolate the
dynamics of the settling filaments to different initial angles within trained
$\mathcal{B}$. Test data set B (red) consisted of filament trajectories at
elasto-gravitational numbers $\mathcal{B}=[150,300,700,1500,3000,7000]$. Test
data set B allows for evaluating the ability of the data-driven models to
interpolate the dynamics of the settling filaments to different initial angles
and elasto-gravitational numbers.

—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—c—  $\mathcal{B}$
$48\theta_{0}/\pi$ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
100 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
150 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
200 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
300 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
400 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
700 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
1000 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
1500 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
2000 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
3000 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
4000 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]
7000 [fill=red] [fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
[fill=red] [fill=red] [fill=red] [fill=red] [fill=red]
10000 [fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green] [fill=yellow] [fill=green]
[fill=yellow] [fill=green] [fill=yellow] [fill=green] [fill=yellow]
[fill=green] [fill=yellow] [fill=green]

Table 1: Distribution of settling filament evolutions into training and
testing data sets by elasto-gravitational number and initial angle of
inclination. The training data set is denoted by the green cells, testing data
set A by the yellow cells, and testing data set B by the red cells.

Figure 5: Block diagram for data-driven model combining the autoencoder and
temporal-evolution scheme. The temporal-evolution neural network, expanded in
red, can be separated into two distinct neural networks forecasting the
evolution of latent representation of the shape, $\mathbf{h}(t+\tau)$, and the
shape-dependent change in position, $\mathbf{c}$; in practice, these can be
forecasted by a single neural network.

A single model was developed to forecast the dynamics of the fiber across all
$\mathcal{B}$ within the scope of this investigation. A block diagram of the
data-driven low-dimensional model architecture is shown in Figure 5. In the
model, the position of the beads in a filament $\hat{\mathbf{X}}(t)$
[connection (a) in the block diagram], normalized by the half filament length
$L/2$, is separated into the shape of the filament $\bar{\mathbf{X}}(t)$
[(b)], where the position of the center bead of the fiber subtracted from each
bead position, and the position of center bead $\mathbf{c}(t)$ [(c)], such
that $\bar{\mathbf{X}}(t)=\hat{\mathbf{X}}(t)-\mathbf{c}(t)$. The terminal
shape of the filament $\bar{\mathbf{X}}_{term}[(d)](\mathbf{B})$ (i.e.
$\bar{\mathbf{X}}(t)$ as $t\rightarrow\infty$) is then subtracted from the
filament shape to give the residual shape of the filament $\check{\mathbf{X}}$
[(e)]. An autoencoder neural network architecture, composed of an encoder and
a decoder neural networks, learns a nonlinear coordinate transformation to
reduce the dimensionality of the fiber shape to a low-dimension latent state
and reconstruct the full shape from the latent state. The encoder
$\mathcal{E}(\check{\mathbf{X}}(t))$ first compresses the shape to a low-
dimensional latent state, $\mathbf{h}(t)$ [(g)], denoted as

$\mathbf{h}(t)=\mathcal{E}(\check{\mathbf{X}}(t);\mathcal{B};\theta_{E}),$ (5)

where the latent dimension $d_{h}$ is much less than the ambient dimension
$d_{\check{\mathbf{X}}}=3N=93$. The decoder neural network,
$\mathcal{D}(\mathbf{h}(t))$ then uses the low-dimensional latent state to
reconstruct the full state of the fiber shape,
$\tilde{\check{\mathbf{X}}}(t)$, such that

$\tilde{\check{\mathbf{X}}}(t)=\mathcal{D}(\mathbf{h}(t);\mathcal{B};\theta_{D})$
(6)

Both the encoder and decoder can receive the elasto-gravitational number as an
additional input [(f) & (k)], the necessity of which will be explored below.

The encoder and decoder networks are trained simultaneously, learning the
neural network weights $\theta_{E}$ and $\theta_{D}$ using gradient descent to
minimize the reconstruction error, given by

$\mathcal{L}\left(\check{\mathbf{X}};\theta_{E},\theta_{D}\right)=\left\langle\left\|\check{\mathbf{X}}-\mathcal{D}\left(\mathcal{E}\left(\check{\mathbf{X}};\mathcal{B};\theta_{E}\right);\mathcal{B};\theta_{D}\right)\right\|_{2}^{2}\right\rangle,$
(7)

where $\langle\cdot\rangle$ is the average over the training data.

A time-integrating neural network is then trained to use the current latent
state and the parameter value $\mathcal{B}$ [(h)] to forecast the evolution of
the latent space one time step $\tau$ into the future, $\mathbf{h}(t+\tau)$
[(i)], while simultaneously forecasting the new position, $\mathbf{c}(t+\tau)$
[(j)].

In this study, we perform time-integration though the use of a neural ordinary
differential equation (ODE), a neural network trained to forecast the right-
hand side of the temporal evolution of the latent state, given as

$\frac{d\mathbf{h}}{dt}=\mathcal{F}_{h}(\mathbf{h}(t);\mathcal{B};\theta_{h})-A\mathbf{h}.$
(8)

The final term on the right-hand side provides an explicit damping, which has
been shown to stabilize the latent state dynamics for neural ODEs[19]. The
evolution of the center bead position can be similarly modeled by a neural
ordinary differential equation, wherein the neural network reconstructs the
right-hand side of

$\frac{d\mathbf{c}}{dt}=\mathcal{F}_{c}(\mathbf{h}(t);\mathcal{B};\theta_{c}),$
(9)

where $\mathbf{c}(0)=\mathbf{0}$. It should be noted that while both the
change in latent shape and position depend on the current latent shape of the
fiber, both are independent of the current position, allowing the model to
forecast the evolution of the filament at any arbitrary position.

The latent state and center bead positional change are then integrated forward
in time to forecast the temporal evolution of the latent state. The latent
state one time step forward, $\mathbf{h}(t+\tau)$, can then be calculated as

$\tilde{\mathbf{h}}(t+\tau)=\mathbf{h}(t)+\int_{t}^{t+\tau}\left(\mathcal{F}_{h}\left(\mathbf{h}\left(t^{\prime}\right);\mathcal{B};\theta_{h}\right)-A\mathbf{h}\right)dt^{\prime},$
(10)

and similarly the position one time step forward, $\mathbf{c}(t+\tau)$, as

$\tilde{\mathbf{c}}(t+\tau)=\mathbf{c}(t)+\int_{t}^{t+\tau}\mathcal{F}_{c}\left(\mathbf{h}\left(t^{\prime}\right);\mathcal{B};\theta_{F}\right)dt^{\prime}.$
(11)

In practice, the evolution of the latent state and positional change can be
forecasted by a single neural ODE,
$\mathcal{F}(\bar{\mathbf{h}}(t);\mathcal{B};\theta_{F})$, where the latent
state and the position vector are concatenated as
$\bar{\mathbf{h}}(t)=[\mathbf{h}(t),\mathbf{c}(t)]$. The time-integrating
neural network is trained to learn the neural network weights $\theta_{F}$
using a standard stochastic gradient descent to minimize the forecasting
error, with a loss given by

$\mathcal{L}\left(\mathbf{h};\theta_{F}\right)=\left\langle\left\|\mathbf{h}(t+\tau)-\tilde{\mathbf{h}}(t+\tau)\right\|_{2}^{2}\right\rangle+\left\langle\left\|\mathbf{c}(t+\tau)-\tilde{\mathbf{c}}(t+\tau)\right\|_{2}^{2}\right\rangle.$
(12)

The two components of the loss are given equal weighting, as the relative
magnitudes of each term are sufficiently comparable that neither dominates in
training. The resulting forecast of the latent state can then be decoded to
produce a predicted latent shape $\tilde{\check{\mathbf{X}}}(t+\tau)$ [(l)]
and combined with the terminal shape filament shape [(m)] to recover the
forecast of the filament shape $\tilde{\bar{\mathbf{X}}}(t+\tau)$ [(n)]. The
forecasts of the shape and position can then be recombined to produce a full
forecast of the evolution of the settling fiber
$\tilde{\hat{\mathbf{X}}}(t+\tau)$ [(o)]. By time-integrating the neural ODE
from a given initial condition, the dynamics of the full filament can be
modeled in the low-dimensional state until the terminal position is reached.

##  III Results and discussion

We first sought to determine the necessary autoencoder architecture with which
to model the evolution of a filament at multiple elasto-gravitational numbers.
We evaluated the performances of four distinct architectures. The first
architecture, denoted as “No $\mathcal{B}$” did not receive the elasto-
gravitational number in either the encoder or the decoder, such that
$\tilde{\check{\mathbf{X}}}=D\left(E\left(\check{\mathbf{X}};\theta_{E}\right);\theta_{D}\right)$
and connections (f) and (k) are removed. The second, referred to as “Encoder
$\mathcal{B}$”, received the elasto-gravitational number as an input solely to
the encoder, such that
$\tilde{\check{\mathbf{X}}}=D\left(E\left(\check{\mathbf{X}};\mathcal{B};\theta_{E}\right);\theta_{D}\right)$
and only connection (f) is removed. Similarly, the third architecture,
referred to as “Decoder $\mathcal{B}$”, received the elasto-gravitational
number as an input solely to the decoder such that
$\tilde{\check{\mathbf{X}}}=D\left(E\left(\check{\mathbf{X}};\theta_{E}\right);\mathcal{B};\theta_{D}\right)$
and only connection (f) is removed. The final architecture, denoted as “Double
$\mathcal{B}$” received the elasto-gravitational number as an input to both
the encoder and the decoder, such that
$\tilde{\check{\mathbf{X}}}=D\left(E\left(\check{\mathbf{X}};\mathcal{B};\theta_{E}\right);\mathcal{B};\theta_{D}\right)$
and both connections (f) and (k) are included. Block diagrams of these
architectures are shown in Figure 6.b-e. Each architecture was used to train
autoencoders at latent dimensions $1\leq d_{h}\leq 6$ with the training data
set. The reconstruction loss across all of the testing data was calculated for
the models at each latent dimension, and the accuracy of the coordinate
transformations learned by the encoders and decoders evaluated.

Figure 6: a) Loss over the testing data for each autoencoder architecture as a
function of latent dimension. Block diagrams for the autoencoder neural
network architectures (b) “No $\mathcal{B}$”, (c) “Encoder $\mathcal{B}$”, (d)
“Decoder $\mathcal{B}$”, and (e) “Double $\mathcal{B}$”.

Shown in Figure 6 is the reconstruction loss for test data sets A and B at
each latent dimension. The performance of each autoencoder architecture is
similar, with the loss dropping precipitously at a latent dimension between
two and four, with higher latent dimensions showing no significant improvement
in the reconstruction loss. A latent dimension of four was selected as the
standard size for the remainder of this study, to ensure maximum accuracy
while still achieving substantial dimension reduction as $d_{h}=4\ll
d_{\check{\mathbf{X}}}=93$, with an additional three degrees of freedom needed
to account for filament position. Of note, none of the autoencoder
architectures resulted in significantly lower reconstruction losses,
indicating that any of the latent state parametrizations can be used to
develop a dynamical model. By heuristic observations, the “Decoder
$\mathcal{B}$” was observed to produce the most consistent results, and as
such was selected for continued use in this study, although it should be noted
that all architectures can produce accurate models.

Once the appropriate autoencoder architecture and required number of degrees
of freedom were determined, the data-driven model utilizing neural ODEs was
applied to forecast the evolution of all filament evolutions in test data sets
A and B. as described in table II.2. The initial bead positions of the
filament in each trajectory were input to the model and the shape and
positional change forecasted using Eqs. 8-9 over the settling period, $T$,
which is dependent on the initial orientation and elasto-gravitational number.
The time-dependent error $E(t)$ for each forecast was calculated as

$E(t)=\frac{\left\langle\left\|\bar{\mathbf{X}}(t)-\tilde{\bar{\mathbf{X}}}(t)\right\|_{2}^{2}\right\rangle}{D_{shape}}+\frac{\left\langle\left\|\bar{\mathbf{c}}(t)-\tilde{\bar{\mathbf{c}}}(t)\right\|_{2}^{2}\right\rangle}{D_{pos}}.$
(13)

The first term in the error formula calculates the error in the forecasted
shape of the filament and is normalized by the error between two filaments
offset by one filament length,
$D_{shape}=\|\bar{\mathbf{X}}-(\bar{\mathbf{X}}+L)\|_{2}^{2}$. The second term
calculates the error in the forecasted position of the center bead, normalized
by the total distance covered by the true data,
$D_{pos}=\left\|\bar{\mathbf{c}}(t)-\bar{\mathbf{c}}(0)\right\|_{2}^{2}$. The
best and worst forecasts at each $\mathcal{B}$ were identified by selecting
the forecasts with the lowest and highest maximum instantaneous error,
respectively. Here, we will show representations of the forecasts of these
selected trajectories, and the accuracy of the remaining trajectories within
the test data sets can be inferred to fall between the bounding forecasts
shown. Videos of all forecasted trajectories in the test data sets can be
found in the Supplemental Material for further elucidation[30].

We first compared the predicted evolution of the shape of the filament over
the settling period to the true trajectory. Shown in Figure 7 are snapshots of
the shape of the filaments from the best and worst forecasts for reconstructed
trajectories in test data set A over the settling period, overlaying the true
shape of the filament at that time. For all trajectories, the predicted shape
of the filament closely matches the true shape across the settling period,
showing our model has successfully interpolated the shape evolution to new
initial angles of orientation. The best forecasts across test data set A show
little deviation from the true evolutions, with the underlying true shapes
nearly indistinguishable from the predictions. The terminal shapes of the
forecasts all trend towards the true terminal shape, demonstrating that the
model has learned the appropriate long time behavior of the residual shape.
The worst forecast of the shape over this test data set occurs at
$\mathcal{B}=10000$, where deviations between the true and reconstructed shape
are apparent; this error is caused by a slight lag in the forecasted
reorientation of the filament, resulting in a time-delay between the
instantaneous shapes. It should be noted that the initial orientations of the
best and worst forecasts are reasonably evenly distributed across the initial
angles of inclinations, showing that our data-driven model performs equally
well across the initial angle distribution and does not favor any obvious
subset of inclinations.

Figure 7: Evolution of the shape of a filament settling in quiescent Newtonian
fluid given an initial orientation from the best and worst forecasts (red) and
the true evolution (black) at a given $\mathcal{B}$. Here, $\mathcal{B}$ is
within the training data set, and the initial angles of inclination are not.

The predicted shapes of the filaments over the settling period in test data
set B were then compared to the true shapes from the trajectories, shown in 8
with the reconstructed shapes overlaid on the true from the best and worst
forecast at each $\mathcal{B}$. As with the forecasts from test data set A,
the forecasted shapes here closely match the true shapes, with only minor
deviations observed across the test data set, demonstrating that our model has
successfully interpolated the shape evolution to new elasto-gravitational
numbers. The reconstructed shapes from the best forecasts almost completely
match the true shapes, with only small errors between the shapes apparent.
Again, the predicted terminal shapes are all indistinguishable from the true,
indicating that our model has correctly identified the long time collapse of
the residual shape at untrained $\mathcal{B}$-dependent terminal shapes. The
worst forecasts in test data set B, here at $\mathcal{B}=$ 150 & 7000, show
noticeable deviations between the true and reconstructed shapes at
intermediate times, due again to a lag in the start of the reorientation
behavior between the two trajectories. Once more, the initial angle of
orientations for the best and worst forecasts show no grouping towards common
$\theta_{0}$, meaning that, at these interpolated $\mathcal{B}$, the model
still does not show preference towards forecasts of any discernible subset of
inclination angles.

Figure 8: Evolution of the shape of a filament settling in quiescent Newtonian
fluid given an initial orientation from the best and worst forecasts (red) and
the true evolution (black) at a given $\mathcal{B}$. Here, neither
$\mathcal{B}$ nor the initial angles of inclination are within the training
data set.

We next investigated the accuracy of the forecasts of the center bead
trajectory by our data-driven model. Shown in Figure 9 is a comparison between
the true and forecasted trajectories of the center bead $\mathbf{c}(t)$ in the
filament evolution from test data set A in 9.a and from test data set B in
9.b, with the trajectories from best and worst forecasts at each $\mathcal{B}$
on the left and right, respectively. As the initial position of the filament
is arbitrary, we have set all trajectories to begin on the $x$-axis, offset in
the $x$-direction for visual clarity. As can be seen, the forecasted
positional evolution of the center bead closely tracks the true positions for
all predicted trajectories. The best forecasts of the trajectories at each
$\mathcal{B}$ in both data sets are nearly identical to the true trajectories,
with only slight positional errors between them; many of the trajectories from
the worst forecasts across the test data sets show a similar degree of
accuracy to the true positions. In some of the worst forecasts, the
reconstructed trajectories are appreciably different from the true; however,
in contrast to the total distance traversed by the center bead, the error
between the true and forecasted positions is comparatively small across the
calculated trajectories.

Figure 9: Trajectory of the center bead $c(t)$ a filament settling in
quiescent Newtonian fluid given an initial orientation from the best (solid
lines) and worst (dashed lines) forecasts and the true evolution at a given
$\mathcal{B}$. The initial positions are denoted by symbol
“$\blacktriangledown$” and the terminal positions are denoted by the symbol
“$\times$”, with the true and predicted positions in black and red,
respectively. In (a), the $\mathcal{B}$ are within the training data set, but
at initial angles of inclination not within the training data; in (b), both
$\mathcal{B}$ and the initial angles of inclination are not within the
training data.

Finally, we investigated the error in the forecast across the investigated
elasto-gravitational numbers in both test data sets. In Figure 10, the
ensemble-averaged error over the settling period is shown at each
$\mathcal{B}$ for test data set A in 10.a and for test data set B in 10.b, as
well as the ensemble-averaged error across the $\mathcal{B}$ in each data set.
At short times, the instantaneous error shows a sudden increase associated
with the positional error, due to small total distance covered compared to the
initial error between the change in position for the true and reconstructed
trajectories; as time increases, the error associated with the positional
error decreases precipitously as the total distance covered grows much larger
than the slight difference in position. At intermediate times, lags between
the true and forecasted reorientation behaviors cause instantaneous error due
to difference between the shapes, causing an increase in error in a limited
number of trajectories. At long times, the error trends toward a small but
finite value, as slight differences between forecasted terminal shape and
terminal velocity cause an unavoidable deviation between the trajectories. At
all $\mathcal{B}$ in both test data sets, the error remains small for the
entire settling period, with each ensemble-averaged error remaining below 0.1
in all cases.

Figure 10: Ensemble-average error vs. time for forecasts of the evolution of a
filament settling in quiescent Newtonian fluid at each $\mathcal{B}$ (gray)
and averaged over all $\mathcal{B}$ (red). In (a), the $\mathcal{B}$ are
within the training data set, but at initial angles of inclination not within
the training data; in (b), both $\mathcal{B}$ and the initial angles of
inclination are not within the training data.

##  IV Conclusions

In this work, we have developed a data-driven modeling technique for learning
a low-dimensional model for the dynamics of an elastic fiber at arbitrary
initial angle of inclination settling under gravity in a quiescent, viscous
Newtonian fluid. We have demonstrated that our proposed neural network
architecture can capture that filament shape in a latent state with
significantly fewer dimensions than the full state space, and we have shown
that the architecture needed to adequately reduce the dimension is aided by,
but does not require, the elasto-gravitational number of the filament as an
additional input. We have shown that the model can forecast the evolution of
filaments at untrained angles of inclination, and can interpolate the
evolution behavior to trajectories of fibers at unseen $\mathcal{B}$. The
forecasts of the shape of the settling fiber closely match the true shapes
across the settling period for all investigated $\mathcal{B}$, with only small
deviations due to a lag in reorientation time. The predicted trajectories of
the center bead position closely match the true trajectories, with the
positional error remaining small in comparison to the total distance covered
by the filament. We have shown that, across all trajectories, in both test
data sets, the instantaneous ensemble-averaged error remains low for the
entire settling period.

In this study, we have developed a technique to model the dynamics of a
flexible particle, and we have demonstrated that our proposed model is able to
accurately forecast the evolution of a simple elastic structure. It remains an
open question as to the optimal architecture with which to model a settling
filament. In this paper, we developed our model at a latent dimension that,
while significantly lower than the ambient dimension of the system, is not
known to be the minimal dimension of this system. Future study can reveal if,
given improved modeling techniques, the necessary latent dimension for this
system can decrease [31]. Additionally, changing the type of neural network
used in the model architecture may lead to more accurate, lower dimensional
models. Here, we only used fully-connected, feed-forward dense neural networks
in our data-driven models. Given the localized nature of the data of the
structure of an elastic particle, using an architecture that utilizes more
localized information processing may produce better models, especially for
objects with more geometric complexity. Recent studies have shown that graph
neural networks have been shown to improve models of multiphase flows[32, 33].
Given the inherent graph structure of a discretized flexible particle, graph
neural networks could leverage this structure to improve model forecasts.

Now that we have demonstrated that our technique can model a simple elastic
particle in a simple flow, future investigation should determine the ability
of such models to forecast the dynamics of more complicated particles in
differing flow fields. For the flexible filament, this can include modeling
the dynamics of such particles in flow fields such as shear flow and
extensional flow. Previous studies have examined the behavior of filaments in
these flow fields using numerical analysis and physics-based simulations[34,
35, 36, 37]. This can then be extended to more complex particles, such as
elastic sheets and deformable spheres, as well as investigating particles with
varying shapes at rest. This can include developing low-dimensional models for
biological systems, where cells suspended in flow behave as elastic particles.
In the long term, such models could be expanded to forecast the evolution of
multiple elastic particles entrained in flow, capturing the inter-particle
hydrodynamic interactions with a single, low-dimensional model.

## References

  * Du Roure  _et al._ [2019] O. Du Roure, A. Lindner, E. N. Nazockdast, and M. J. Shelley, Dynamics of flexible fibers in viscous flows and fluids, Annual Review of Fluid Mechanics 51, 539 (2019).
  * Yu and Graham [2023] Y. Yu and M. Graham, Free-space and near-wall dynamics of a flexible sheet sedimenting in stokes flow, arXiv , 2310.08722 (2023).
  * Shelley [2016] M. J. Shelley, The dynamics of microtubule/motor-protein assemblies in biology and physics, Annual Review of Fluid Mechanics 48, 487 (2016).
  * Lundell  _et al._ [2011] F. Lundell, L. D. Soderberg, and P. H. Alfredsson, Fluid mechanics of papermaking, Annual Review of Fluid Mechanics 43, 195 (2011).
  * Graham [2011] M. D. Graham, Fluid dynamics of dissolved polymer molecules in confined geometries, Annual Review of Fluid Mechanics 43, 273 (2011).
  * Yu and Graham [2021] Y. Yu and M. Graham, Coil-stretch-like transition of elastic sheets in extensional flows, Soft Matter 17, 543 (2021).
  * Yu and Graham [2022] Y. Yu and M. Graham, Wrinkling and multiplicity in the dynamics of deformable sheets in uniaxial extensional flow, Physical Review Fluids 7, 023601 (2022).
  * Xu and Nadim [1994] X. Xu and A. Nadim, Deformation and orientation of an elastic slender body sedimenting in a viscous liquid, Physics of Fluids 6, 2889 (1994).
  * Schlagberger and Netz [2005] X. Schlagberger and R. R. Netz, Orientation of elastic rods in homogeneous stokes flow, Europhysics Letters 70, 129 (2005).
  * Lagomarsino  _et al._ [2005] M. C. Lagomarsino, I. Pagonabarraga, and C. P. Lowe, Hydrodynamic induced deformation and orientation of a microscopic elastic filament, Physical Review Fluids 94, 148104 (2005).
  * Li  _et al._ [2013] L. Li, F. Manikantan, D. Saintillan, and S. E. Spagnolie, The sedimentation of flexible filaments, Journal of Fluid Mechanics 735, 705 (2013).
  * Shojaei and Dehghani [2015] B. Shojaei and H. Dehghani, The sedimentation of slim flexible particles in stokes flow, Indian Journal of Science and Technology 8, 1 (2015).
  * Delmotte  _et al._ [2015] B. Delmotte, E. Climent, and F. Plouraboué, A general formulation of bead models applied to flexible fibers and active filaments at low reynolds number, Hournal of Computational Physics 286, 14 (2015).
  * Marchetti  _et al._ [2018] B. Marchetti, V. Respa, A. Lindner, O. du Roure, E. Bergougnoux, L. Guazzelli, and C. Duprat, Deformation of a flexible fiber settling in a quiescent viscous fluid, Physical Review Fluids 3, 104102 (2018).
  * Cunha  _et al._ [2022] L. H. P. Cunha, J. Zhao, F. C. MacKintosh, and S. L. Biswal, Settling dynamics of Brownian chains in viscous fluids, Physical Review Fluids 7, 034303 (2022).
  * Cheng  _et al._ [2023] X. Cheng, C. Caruso, W. A. Lam, and M. D. Graham, Marginated aberrant red blood cells induce pathologic vascular stress fluctuations in a computational model of hematologic disorders, Science Advances 9, eadj6423 (2023).
  * Ebrahimi and Bagchi [2023] S. Ebrahimi and P. Bagchi, A computational study of red blood cell deformability effect on hemodynamic alteration in capillary vessel networks, Scientific Reports 12, 4304 (2023).
  * Linot and Graham [2022] A. J. Linot and M. D. Graham, Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations, Chaos 32, 073110 (2022).
  * Linot  _et al._ [2023] A. J. Linot, J. W. Burby, Q. Tang, P. Balaprakash, M. D. Graham, and R. Maulik, Stabilized neural ordinary differential equations for long-time forecasting of dynamical systems, Journal of Compuational Physics 474, 111838 (2023).
  * Vinuesa and Brunton [2022] R. Vinuesa and S. L. Brunton, Enhancing computational fluid dynamics with machine learning, Nature Computational Science 2, 358 (2022).
  * Fox  _et al._ [2023] A. J. Fox, C. R. Constante-Amores, and M. D. Graham, Predicting extreme events in a data-driven model of turbulent shear flow using an atlas of charts, Physical Review Fluids 8, 094401 (2023).
  * Linot and Graham [2020] A. J. Linot and M. D. Graham, Deep learning to discover and predict dynamics on an inertial manifold, Physical Review E 101, 062209 (2020).
  * Omata and Shirayama [2019] N. Omata and S. Shirayama, A novel method of low-dimensional representation for temporal behavior of flow fields using deep autoencoder, AIP Advances 9, 015006 (2019).
  * Zeng  _et al._ [2023a] K. Zeng, C. E. De Jesus, Perez De Jesus, A. J. Fox, and M. D. Graham, Autoencoders for discovering manifold dimension and coordinates in data from complex dynamical systems, arXiv , 2305.01090 (2023a).
  * Perez De Jesus and Graham [2022] C. E. Perez De Jesus and M. D. Graham, Data-driven low-dimensional dynamic model of Kolmogorov flow, arXiv , 2210.16708 (2022).
  * Brunton  _et al._ [2020] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Machine learning for fluid mechanics, Annual Review of Fluid Mechanics 52, 477 (2020).
  * Zeng  _et al._ [2023b] K. Zeng, A. J. Linot, and M. D. Graham, Data-driven control of spatiotemporal chaos with reduced-order neural ode-based models and reinforcement learning, Proceedings of the Royal Society A 478, 20220297 (2023b).
  * Srinivasan  _et al._ [2019] P. A. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, and R. Vinuesa, Predictions of turbulent shear flows using deep neural networks, Physical Review Fluids 4, 054603 (2019).
  * Graham [2018] M. D. Graham,  _Microhydrodynamics, Brownian Motion, and Complex Fluids_, 1st ed., Cambridge University Press (Cambridge University Press, 2018).
  * [30] See supplemental material at [url will be inserted by publisher] for videos of the forecasted in the test data sets. each video shows the full evolution of the true and forecasted trajectories, with the evolution of the shape and position of the filaments shown alongside for viewing convenience. each video is titled by the elasto-gravitational number, $\mathcal{B}$ and the normalized initial angle of inclination, $48\theta_{0}/\pi$.
  * Zeng  _et al._ [2024] K. Zeng, C. E. Perez De Jesus, A. J. Fox, and M. D. Graham, Autoencoders for discovering manifold dimension and coordinates in data from complex dynamical system, arXiv , 2305.01090v3 (2024).
  * Ma and Pan [2024] Z. Ma and W. Pan, Shape deformation, disintegration, and coalescence of suspension drops: Efficient simulation enabled by graph neural networks, International Journal of Multiphase Flow 178, 104845 (2024).
  * Ma  _et al._ [2022] Z. Ma, Z. Ye, and W. Pan, Fast simulation of particulate suspensions enabled by graph neural network, Computer Methods in Applied Mechanics and Engineering 400, 115496 (2022).
  * Słowicka  _et al._ [2012] A. M. Słowicka, M. L. Ekiel-Jeżewska, K. Sadlej, and E. Wajnryb, Dynamics of fibers in a wide microchannel, Journal of Compuational Physics 136, 044904 (2012).
  * Słowicka  _et al._ [2013] A. M. Słowicka, E. Wajnryb, and M. L. Ekiel-Jeżewska, Lateral migration of flexible fibers in poiseuille flow between two parallel planar solid walls, The European Physical Journal E 36, 31 (2013).
  * Słowicka  _et al._ [2015] A. M. Słowicka, E. Wajnryb, and M. L. Ekiel-Jeżewska, Dynamics of flexible fibers in shear flow, The Journal of Chemical Physics 143, 124904 (2015).
  * Farutin  _et al._ [2016] A. Farutin, T. Piasecki, A. M. Słowicka, C. Misbah, E. Wajnryb, and M. L. Ekiel-Jeżewska, Dynamics of fibers in a wide microchannel, Soft Matter 12, 7307 (2016).