# Label Assisted Autoencoder for Anomaly Detection in Power Generation Plants

Marcellin Atemkeng  Victor Osanyindoro  Rockefeller Rockefeller  Sisipho
Hamlomo  Jecinta Mulongo  Theophilus Ansah-Narh  Franklin Tchakounte  Arnaud
Nguembang Fadja  Department of Mathematics, Rhodes University, 6139 Makhanda,
South Africa  African Institute for Mathematical Sciences, Limbe, Cameroon
Department of Statistics, Rhodes University, 6139 Makhanda, South Africa
Ghana Space Science and Technology Institute, Ghana Atomic Energy Commission,
Accra, Ghana  Department of SFTI, School of Chemical Engineering and Mineral
Industries, University of Ngaoundéré, Cameroon  Dipartimento di Ingegneria,
University of Ferrara, Via Saragat 1, 44122 Ferrara, Italy

###### Abstract

One of the critical factors that drive the economic development of a country
and guarantee the sustainability of its industries is the constant
availability of electricity. This is usually provided by the national electric
grid. However, in developing countries where companies are emerging on a
constant basis including telecommunication industries, those are still
experiencing a non-stable electricity supply. Therefore, they have to rely on
generators to guarantee their full functionality. Those generators depend on
fuel to function and the rate of consumption gets usually high, if not
monitored properly. Monitoring operation is usually carried out by a (non-
expert) human. In some cases, this could be a tedious process, as some
companies have reported an exaggerated high consumption rate. This work
proposes a label assisted autoencoder for anomaly detection in the fuel
consumed by power generating plants. In addition to the autoencoder model, we
added a labelling assistance module that checks if an observation is labelled,
the label is used to check the veracity of the corresponding anomaly
classification given a threshold. A consensus is then reached on whether
training should stop or whether the threshold should be updated or the
training should continue with the search for hyper-parameters. Results show
that the proposed model is highly efficient for reading anomalies with a
detection accuracy of 97.20% which outperforms the existing model of 96.1%
accuracy trained on the same dataset. In addition, the proposed model is able
to classify the anomalies according to their degree of severity.

###### keywords:

Electric grid , Fuel consumption , Autoencoder, Label assisted autoencoder ,
Anomaly detection , Power generating plants

##  1 Introduction

About 3% of the world’s electrical energy is utilised by the information
communication technology companies [1] and the telecommunication industry is
one of the fastest growing industrial sectors among Agriculture, Banking,
Infrastructure, and Oil and Gas. The number of telecommunication industries
and the quest for expansion and growth has led to an increase in base stations
across targeted countries to boost their network coverage and enhance the
effective flow of communication. With the increase in the number of base
stations, the issue of base station management needs to be addressed. Grid
energy has been known to be the main source of power in developing countries
such as those in Africa, and it is expected that these base stations located
across different rural and urban areas will be powered by grid energy.
However, electricity is quite unstable in most parts of these developing
countries and this has forced base stations to look for a reliable alternative
source of energy. These alternatives include photovoltaic panels (PV), wind
turbines and diesel generators, but mostly generators due to a lack of space
for the installation of PV or wind turbines [2]. The high cost of fuel and its
transportation to supply stations located in rural areas has increased the
operational cost of these companies. These generators are being refilled
manually, thus creating room for irregular or unusual fuel consumption which
might be caused by several reasons such as fuel theft, fuel leakages or poor
maintenance of equipment. A study conducted in a base station in Cameroon has
shown that the design of the base station building, room cooling systems such
as air conditioners and careless handling of lights increase the rate of fuel
consumption in the base station [3]. Espadafor et al. [4] also attested that
generator performance can be affected by the age of the generator, the number
of loads powered by the generator and improper maintenance.

In Cameroon, TeleInfra LTD is one such company whose objective is to manage
base stations in various parts of the country. The services include
maintenance of base stations and refuelling of generators. Like any other
businesses relying on grid energy, the unstable power supply has resulted in
high operating costs as companies have to find alternative sources of power
supply to sustain the continuous operations of the business. The use of
alternative sources such as solar panels, hybrid energy and generators has
been implemented by TeleInfra for the sustainability of the business
performance. Data on fuel consumed such as working hours of the generator, the
quantity of fuel refuelled, the rate of consumption, generator maintenance,
and total fuel consumed, are collected from base stations [5].

Although many strategies have been proposed to improve energy saving such as
building a well-ventilated base station, the use of air conditioners as a
cooling system and heat pipes to remove hot air from the base station [6].
However, detecting such irregularities is challenging especially when there
are numerous base stations functioning at once. Anomaly detection in power
generating plants is aimed at detecting such irregularities in the behaviour
of the data provided. Although there are different algorithms for detecting
anomalies, machine learning algorithms are the most used and popular for
anomaly detection due to their ability for automation and their effectiveness
in the context of deep learning, especially when involving large datasets.
According to Goodfellow [7], deep learning is a variant of applied statistics
with an increasing emphasis on using computers in estimating complex functions
statistically, but a reduced emphasis on proving confidence intervals around
these functions. Machine learning algorithms come in several variants. In
supervised learning, models can make predictions on unlabelled data after they
have been trained on labelled data whereas, in unsupervised learning, models
can only make predictions on unlabelled data by learning similar features and
patterns embedded in the datasets. In reinforcement learning, a goal is given
and an agent undergoes training in an environment with the purpose of finding
an optimal solution to accomplishing the goal.

Anomalies can be an indicator of areas that require attention, and detecting
them has been quite popular among the research community. In the past, the
research community has conducted several anomaly detections ranging from
comprehensive to certain application domains with different machine learning
algorithms. Mulongo et al. [5] worked on a similar area, four different
supervised learning algorithms were used in their work for detecting an
anomaly in a power generation plant. However, in real-life scenarios, abnormal
behavioural patterns are very few compared to normal behaviour. For instance,
in the same work of Mulongo et al. [5], labels were created by the authors
based on certain criteria and only 35 per cent of the dataset turned out to be
anomalous and the authors had to duplicate the data to balance between the two
categories. These are key challenges in recognising anomalies with a
supervised learning algorithm since they greatly rely on labelling and balance
between normal and anomalous data patterns. The goal is to learn how
efficiently a generating plant should behave so that observations that do not
confine to the norm will be identified, and the required attention can be
given. Our work investigates an alternative approach for anomaly detection in
a power generation plant in an unsupervised manner based on the dataset used
in Mulongo et al. [5]. The aim is to compare performance with the results
obtained in Mulongo et al. [5]. The proposed unsupervised learning framework
is built from a modifying and fine-tuning autoencoder which is free of the
hassle of data labelling and balancing. An additional module is added to an
autoencoder, the new module uses some labelled data to check if each of the
observations is correctly classified and then one of the three steps below is
activated:

  * 1.

Update the threshold of the autoencoder which therefore increases the model
overall accuracy to an acceptable level

  * 2.

Update the interval of variation for numerical hyper-parameters, the best
values of the hyper-parameters are then explored in the new search space

  * 3.

Provide overall performance score

This paper is organised as follows: Section 2 explains anomaly detection and
investigates related works for anomaly detection in power grid plants. Section
3 discusses autoencoders which replicate the input data through a compressed
representation. This section also discusses the different evaluation metrics
used in this work. Section 4 proposes the label assisted autoencoder and
provides a detail discussion. The dataset and feature engineering are
discussed in Section 5. Section 6 discusses the results and limitations, and
Section 7 concludes the work.

##  2 Anomaly detection and related works

Anomalies also known as outliers often refer to instances or data samples that
are significantly distanced from the main body of an examined data [8]. These
distanced values often indicate a deviation from its established normal
pattern which can sometimes be a measurement error or an indication of a data
sample of a different population [9]. Outliers classification depends on the
type and domain of the given data as well as the data analyst. Since many
outliers are linked directly with abnormal behaviour, they are also referred
to as deviants, anomalies, or abnormalities in the literature of statistics
and data analysis [8]. According to [8], interpreting data is directly
associated with the detection of anomalous samples. Demestichas et al. [9]
suggested that it is essential to achieve the highest possible
interpretability level to properly select the best anomaly detection method
from different ranges of the relevant algorithm. There are two major
categories of anomalies depending on the given dataset; multivariate and
univariate [9]. Multivariate anomalies can be spotted in multi-dimensional
data while univariate anomalies are spotted in a single-dimensional data.
Besides the two categories of anomaly, there are other categories which depend
on the distribution of the given data. Data samples that are considered
anomalous when viewed against the entire dataset are point anomalies, while
data samples that are considered anomalous with respect to meta-information
related to the data sample are contextual anomalies [10]. In other words,
contextual anomalies are classified based on local neighborhoods, while point
anomalies are classified based on the overall dataset. Collective anomalies
denote anomalous data collection samples which together are considered an
anomalous pattern.

Fahim and Sillitti [11] provide two anomaly detection methods; statistical and
machine learning methods. The statistical method uses various algorithms such
as density-based, distance-based, parametric and statistical-based. However,
Trinh et al. [12] noted that one of the major challenges that are encountered
by this approach is the design of a suitable model that can accurately
separate normal data from unusual data points. On the other hand, machine
learning methods consist of both supervised and unsupervised learning
algorithms in which dataset can either be labelled for supervised learning or
unlabelled for unsupervised learning. Some advantages of this method are an
enhancement of detection speeds and its ability to handle complexity with the
less human intervention [13].

Many researchers have worked on different machine learning techniques for
anomaly detection, but most of the current works applied artificial neural
networks (ANNs) to classification tasks. The labelled data is used during the
training stage, and then the learned model is able to correctly classify
sample data never used during the training process. This technique is
generally classified under supervised machine learning techniques. Such an
example is trained in Mulongo et al [5] in which support vector machines (SVM
[14]), K-Nearest Neighbours (KNN [15]), Logistic Regression (LR [16]), and
MultiLayer Perceptron (MLP [17]) are used for anomaly detection associated
with the fuel consumed dataset from an energy company. However, the energy
sector is not the only place anomaly detection with supervised machine
learning has been applied, others include fraud detection in credit card [18]
attacks and anomaly detection in IoT sensors [19]. One of the main advantages
of supervised learning techniques is the ability to handle high-dimensional
datasets with high-performance [20]. However, there is a major problem with
this technique. When dealing with real life data, the majority of them
contains fewer anomalies data which is quite challenging and can cause the
problem of an unbalanced dataset. This is an issue for supervised learning
techniques since they greatly rely on labelled and balanced data. However,
unsupervised learning techniques can be used to address this problem, for
example, the autoencoder considers a specific kind of feed-forward neural
network that can be applied in outlier-based anomaly detection rather than
classification problems. Hawkins et al. [21] proposed an approach that
involved autoencoder for outlier detection, however many researchers have
investigated hybrid methods; e.g. [20] proposed an approach based on long
short term memory (LSTM) autoencoder and One-class SVM (OC-SVM). The approach
is used to detect anomalies based attacks in an unbalanced dataset. The idea
is to use the LSTM-autoencoder to train a model to learn the pattern in the
normal class (dataset without anomaly) so that-the model is able to replicate
the input data at the output layer with a small reconstruction error. When
there are anomalies in the data, the model fails to replicate the anomalous
samples. This arises when the reconstruct error is very high.

Another unsupervised learning technique is the k-means, Zhang et al. [22],
used the transformer model and the k-means clustering method for anomaly
detection. The k-means was also used in the work of Munz et al. [23], to
detect traffic anomalies, the main idea is to train data containing unlabelled
records and separate them into clusters of normal and anomalous data.

##  3 Autoencoder and performance metrics

###  3.1 Autoencoder

An autoencoder is a neural network that is trained to attempt to copy its
input to its output in an unsupervised manner. It consists of three parts: An
Encoder part, a code and a decoder part. The encoder compresses the input data
$\mathbf{X}$ defined in two dimensions as:

$\displaystyle\mathbf{X}$
$\displaystyle=\begin{pmatrix}\mathbf{X}_{1,1}&\mathbf{X}_{1,2}&\cdots&\mathbf{X}_{1,M}\\\
\mathbf{X}_{2,1}&\mathbf{X}_{2,2}&\cdots&\mathbf{X}_{2,M}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{X}_{N,1}&\mathbf{X}_{N,2}&\cdots&\mathbf{X}_{N,M}\end{pmatrix}\in\mathbb{R}^{N\times
M}$ (1)

where $\mathbf{X}_{:,i}\in\mathbb{R}^{N}$ is an observation with
$N\in\mathbb{R}$ entries among the $M\in\mathbb{R}$ observations. The encoder
compressed $\mathbf{X}$ to a lower-dimensional space
$\mathbf{H}_{:,i}\in\mathbb{R}^{Q}$. The decoder produces a predicted output
$\hat{\mathbf{X}}$,

$\displaystyle\hat{\mathbf{X}}$
$\displaystyle=\begin{pmatrix}\hat{\mathbf{X}}_{1,1}&\hat{\mathbf{X}}_{1,2}&\cdots&\hat{\mathbf{X}}_{1,M}\\\
\hat{\mathbf{X}}_{2,1}&\hat{\mathbf{X}}_{2,2}&\cdots&\hat{\mathbf{X}}_{2,M}\\\
\vdots&\vdots&\ddots&\vdots\\\
\hat{\mathbf{X}}_{N,1}&\hat{\mathbf{X}}_{N,2}&\cdots&\hat{\mathbf{X}}_{N,M}\end{pmatrix}\in\mathbb{R}^{N\times
M}$ (2)

by reconstructing the original data from the compressed representation. It
makes use of multiple layers and uses non-linear activation functions to learn
the non-linear relationship embedded in the data. The goal is to make the
reconstruction error as minimal as possible, which means finding the
parameters that make the reconstruction $\hat{\mathbf{X}}$ as close as
possible to the original input $\mathbf{X}$. Autoencoder can be applied to
various tasks such as anomaly detection [20, 24], generative model [25, 26],
clustering [27, 28], classification [29], recommendation systems [30] and
dimensionality reduction [31, 32]. Figure 1 depicts a simple architecture of
an autoencoder for which the detailed description is as follows:

Encoding: During the encoding process, assuming the input data
$\mathbf{X}_{:,i}$ is a high dimensional vector that is mapped to a low
dimensional vector $\mathbf{H}_{:,i}$ after filtering insignificant features.
This is expressed mathematically as:

$\mathbf{H}_{:,i}=f(\mathbf{X}_{:,i}),$ (3)

where $f$ is a neural network that is trained with sets of activation
functions weights and biases. Note that
$f\in\mathbb{R}^{N}\rightarrow\mathbb{R}^{Q}$, where $Q$ is the number of
compressed representations.

Decoding: During the decoding process, the compressed representation
$\mathbf{H}_{:,i}$ of $\mathbf{X}_{:,i}$ is used to generate the output
$\hat{\mathbf{X}}_{:,i}$ that maps back into the reconstruction of
$\mathbf{X}_{:,i}$:

$\displaystyle\hat{\mathbf{X}}_{:,i}$ $\displaystyle=g(\mathbf{H}_{:,i})$ (4)
$\displaystyle=g\big{(}f(\mathbf{X}_{:,i})\big{)},$ (5)

where $g$ is the decoding neural network with activation functions, weights,
and biases that could be completely independent of the corresponding
activation, weights, and biases of the encoding neural network. As discussed
in [33, 34], the autoencoder involves the search for $f$ and $g$ which
minimizes the average of the loss $\Delta_{\theta,\phi}$ between all the
samples $\mathbf{X}_{:,i}$, added for example, to a $l_{2}$ weighted
regularization term $a_{j}$:

$\displaystyle[\hat{f},\hat{g}]$
$\displaystyle=\operatorname*{arg\,min}_{f,g}\bigg{(}<\Delta_{\theta,\phi}\big{(}\mathbf{X}_{:,i},\hat{\mathbf{X}}_{:,i}\big{)}>+\lambda\sum_{j}a_{j}^{2}\bigg{)},$
(6)
$\displaystyle=\operatorname*{arg\,min}_{f,g}\bigg{(}<\Delta_{\theta,\phi}\big{(}\mathbf{X}_{:,i},g(f(\mathbf{X}_{:,i}))\big{)}>+\lambda\sum_{j}a_{j}^{2}\bigg{)},$
(7)

where $<\cdot>$ is the average operator and $\lambda$ is the hyper-parameter
that weights the regularization term, $a_{j}$. The higher the value of
$\lambda$ the greater the capacity penalty. Note that $j$ runs across the
hidden layers of the neural networks $f,g$. Finding the optimal neural
networks $f$ and $g$ involves updating their respective learning parameters
$\theta,\phi$ in each hidden layer of the neural networks so that the loss
$\Delta_{\theta,\phi}$ is smaller than a given limit. In this work, the mean
absolute loss:

$<\Delta_{\theta,\phi}\big{(}\mathbf{X}_{:,i},\hat{\mathbf{X}}_{:,i}\big{)}>=\frac{1}{M}\sum_{i=1}^{M}\bigg{|}\mathbf{X}_{:,i}-\hat{\mathbf{X}}_{:,i}\bigg{|}$
(8)

is used at the place of $<\cdot>$.

Figure 1: Simple architecture of an autoencoder (adapted from [35]).

###  3.2 Performance metrics

Performance metrics measure the ability of the overall quality of the model. A
single performance measure is not enough to validate an autoencoder,
therefore, different measures are performed and evaluated. The confusion
matrix as shown in Table 1 generates more meaningful measures to find the
detection accuracy, precision, recall, and F1 score. True Normal, $TN$
represents the number of observations in the normal class that are predicted
as normal by the model (i.e. below the threshold). True Anomaly, $TA$ is the
number of observations in the anomaly class that are predicted as an anomaly,
and are above the threshold. False Normal, $FN$ is the number of anomalous
observations that are below the threshold (i.e. predicted as normal classes).
False Anomaly, $FA$ is the number of normal observations that are above the
threshold (i.e. predicted as an anomaly).

Classification Class Distribution
---
|  Actual Normal |  Actual Anomaly
Predicted Normal |  True Normal |  False Normal
Predicted Anomaly |  False Anomaly |  True Anomaly
Table 1: Two class classification confusion matrix representation

The classification accuracy measures the general performance of the model by
producing the ratio of true prediction (true normal and true anomaly) out of
the total number of prediction:

$\displaystyle Accuracy=\frac{TN+TA}{TA+TN+FA+FN}.$ (9)

The Precision is the ratio of the true anomaly divided by the total number of
observations above the threshold (i.e. number or anomalies predicted):

$\displaystyle Precision=\frac{TA}{TA+FA}.$ (10)

The False positive rate, $FPR$ refers to the ratio of the normal samples above
the threshold to the actual number of anomalous samples:

$\displaystyle FPR=\frac{FA}{FA+TN}.$ (11)

True positive rate, $TPR$ also known as sensitivity or recall, is the ratio of
the number of anomalous samples above the threshold to the actual number of
samples in the anomaly class:

$\displaystyle TPR=\frac{TA}{TA+FN}.$ (12)

Specificity is a measure obtained from the outcome of the confusion matrix
which gives the ratio of true anomaly to the total negative class in the
sample:

$\displaystyle\text{Specificity}=\frac{TN}{TN+FA}.$ (13)

F-measure or $F_{1}$-score gives the harmonic mean between recall and
precision of the classifier. High F-measure indicates a better performance of
the classifier with no false alarm:

$\displaystyle
F\text{-measure}=2\left(\dfrac{\text{Precision}\times\text{Recall}}{\text{Precision}+\text{Recall}}\right).$
(14)

##  4 Proposed label assisted autoencoder

###  4.1 The architecture of the proposed model

The encoder takes as input a high-dimensional input data here a fixed-size
vector and reaches the latent space through mapping it to a a low-dimensional
representational vector. The decoder reconstructs the input data from the
reduced representation in the latent space. The final reconstruction error is
used to set a threshold to detect anomalies. An additional computation block
is added to check from a set of labelled data if the threshold is acceptable
to satisfy the required precision. This is a kind of validation block that
uses the threshold to decide if the threshold should be changed or if the
autoencoder should be trained to further minimize the reconstruction loss. The
entire architecture of the proposed assisted autoencoder is depicted in Figure
2 and a detailed description is provided below.

Figure 2: Overview of the proposed model. Each observation goes through the
autoencoder for training, the reconstruction error is measured and then a
threshold is used to make the decision whether the observation is an anomaly
or not. The labelling assistance module then takes over to check if each of
the observations is labelled, and then the label is used to check the veracity
of the corresponding anomaly classification. A consensus is then reached
whether training should stop or whether the threshold should be updated or the
training should continue with the search for hyper-parameters.

  1. 1.

The dataset is made up of observations and features which result in a
2-dimensional array of size $N\times M$, where $N$ and $M$ represent the
number of features and the number of observations respectively. All features
of each observation are collected and received as input in the algorithm
called nodes.

  2. 2.

The input data of size $N\times M$ is collected and reduced to a latent form
of size $N\times Q$ where $Q<M$, while the $N\times Q$ set is fed into the
decoder to produce a predicted output of size $N\times M$. Appendix Appendix A
shows the parameters and the architecture of the deep neural networks trained
in the encoder and decoder. The number of filters, size of filters and layers
are displayed in the encoding and decoding phases.

  3. 3.

After the decoding process, a reconstruction error is produced for each data
point. The reconstruction error refers to the measure of how much the
reconstructed input deviates from the original input. A threshold is set as a
decision point to decide the acceptable amount of deviation, and the features
of the observation that go beyond this threshold is classified as an anomaly.
The observation that is below this threshold is a normal data without anomaly.

  4. 4.

The label assisting module then takes over to verify if each of the
observations is labelled, then the label is used to verify the veracity of the
corresponding anomaly classification. The labelled observations are then
checked and an agreement is reached if the threshold is satisfied to obtain
the desired precision to detect anomalies, if not the threshold is updated or
the model is further trained to find the best hyper-parameters for the given
threshold.

To illustrate the above steps, an example of a scenario is described as
follows. Assume the input data is $\mathbf{X}_{:,i}$ with $N$ features and
their corresponding predictions given as:

$\displaystyle\mathbf{X}_{:,i}$
$\displaystyle=\mathbf{X}_{1,i},\mathbf{X}_{2,i},\mathbf{X}_{3,i},\cdots,\mathbf{X}_{N-1,i},\mathbf{X}_{N,i}$
(15) $\displaystyle=(0.1,0.2,0.3,\cdots,0.4,0.5)$ (16)
$\displaystyle\hat{\mathbf{X}}_{:,i}$
$\displaystyle=\hat{\mathbf{X}}_{1,i},\hat{\mathbf{X}}_{2,i},\hat{\mathbf{X}}_{3,i},\cdots,\hat{\mathbf{X}}_{N-1,i},\hat{\mathbf{X}}_{N,i}$
(17) $\displaystyle=(0.6,0.21,0.32,\cdots,0.61,0.53)$ (18)

The value of each of these features indicates the importance of the feature in
the anomaly classification; e.g. the running time of a generator is more
important than the generator capacity. The reconstruction loss for each
feature is then calculated as in Equation 7:

$\displaystyle L_{:,i}$
$\displaystyle=|\mathbf{X}_{:,i}-\hat{\mathbf{X}}_{:,i}|$ (19)
$\displaystyle=(0.5,0.01,0.02,0.21,0.03).$ (20)

Assume the maximum reconstruction loss is set to $0.2$, then the absolute
error classifies the sample under the normal category while if a priority is
given to the important feature, say $\mathbf{X}_{N-1,i}$ then the sample would
be labelled as an anomaly since the reconstruction loss of
$\mathbf{X}_{N-1,i}$ is $L_{N-1,i}=0.21$ which is beyond $0.2$. Figure 3
illustrates the reconstruction loss scenario adopted in this example.

###  4.2 Training and testing phases

Both the training and testing phases are equally important when it comes to
the generalisation of the model. The dataset is divided into training and
testing, depending on the proportion of choice, in our case, we used the ratio
of 3:1 for training and testing respectively. The training phase is divided
into two folds. Firstly, minimising the reconstruction error is the focus of
the training for the reconstructed outputs to converge to the input samples.
Secondly, the reconstruction error is calculated for each data point to find
an optimal threshold for anomaly detection.

The training phase consists of normalizing the dataset so that all the
features are reduced to a common scale without distorting the differences in
the range of the values. In this work, the mathematical measure used to
normalize the data is given by:

$\mathbf{X}_{:,\mathrm{iscaled}}=\frac{\mathbf{X}_{:,i}-\mathbf{X}_{:,\mathrm{min}}}{\mathbf{X}_{:,\mathrm{max}}-\mathbf{X}_{:,\mathrm{min}}},$
(21)

where $\mathbf{X}_{:,\mathrm{min}}$ and $\mathbf{X}_{:,\mathrm{max}}$
represent the data point with the minimum and maximum entries respectively.
This is followed by dividing and separating the anomalous samples from the
training set so that the algorithm learns to reconstruct the norms. Once the
training is completed, a reconstruction loss between the input and output is
calculated and a backpropagation strategy is applied to adjust the weights and
parameters of the model. The testing phase checks the performance of the model
on the unseen test data using the threshold obtained from the training phase.
The whole process is described in Figure 3.

Figure 3: Computing reconstruction loss for
$\mathbf{X}_{:,i}=(\mathbf{X}_{1,i},\mathbf{X}_{2,i},\mathbf{X}_{3,i},\cdots,\mathbf{X}_{N-1,i},\mathbf{X}_{N,i})$
using the predicted output
$\hat{\mathbf{X}}_{:,i}=(\hat{\mathbf{X}}_{1,i},\hat{\mathbf{X}}_{2,i},\hat{\mathbf{X}}_{3,i},\cdots,\hat{\mathbf{X}}_{N-1,i},\hat{\mathbf{X}}_{N,i})$.

##  5 Datasets

###  5.1 Data description

The dataset used in this paper is gathered from a Telecom base station
management company in Cameroon named TeleInfra and was subject to previous
studies for anomaly detection in [5]. The dataset is collected over the period
of one year, i.e. from September 2017 to September 2018. It consists of $6010$
observations from various base stations in 46 towns and villages (known as
clusters) across Cameroon. These stations mainly rely on generators as the
main supply of power. The dataset also consists of 17 variables which are
defined in both numerical and categorical forms. A detailed description of
each variable is shown in Table 2.

Feature Description
---
CONSUMPTION HIS |  The total fuel consumed between a specific period of time before the next refuelling is done.
CONSUMPTION_RATE |  The number of litres the generator consumes per hour.
Cluster |  The cities where the generator sites are located
CURRENT HOUR METER GE1 |  The hour meter reading of the generator.
Site Name |  Name of the site where each generator is located
EFFECTIVE_DATE_OF_VISIT |  The date of meter reading, refuelling and recording
PREVIOUS_DATE_OF_VISIT |  The previous date of visit
Months |  The month when the reading was taken
NUMBER_OF_DAYS |  The number of days before the next refuelling process.
GENERATOR_1_CAPACITY_(KVA) |  The capacity of the generator
POWER TYPE |  Type of power used in the power plant
PREVIOUS HOUR METER G1 |  The previous meter reading of the generator.
PREVIOUS_FUEL_QTE |  The total quantity of fuel left inside the generator tank on the previous date of the visit.
QTE_FUEL_FOUND |  The quantity of fuel found inside the generator tank before refuelling is done.
QTE_FUEL_ADDED |  The quantity of fuel added to the generator during refuelling process.
TOTALE_QTE_LEFT |  Quantity left in the generator after refuelling.
RUNNING_TIME |  The total number of hours the generator worked before the next refuelling is done
Table 2: Description of the different features in the dataset.

Anomalies are observed in different features of the dataset and the observed
anomalies are classified based on three indicators: (1) for a given time
period if the generator running time is zero and the quantity of fuel consumed
is not zero. (2) when the running time per day is more than 24 hours, and (3)
when the daily consumed quantity of fuel is more than the maximum consumption
a generator can consume. For a data sample to receive the anomaly tag 1, it
has to demonstrate at least one of the three anomaly indicators listed,
otherwise it is given the normal tag 0. A full workflow of the entire
labelling process is illustrated in Figure 4. During the labelling process,
output variables are assigned labels $0$ and $1$ representing the normal and
anomaly classes respectively. For a single generator, Figure 5 illustrates the
plot of the working hours per day, for example, all data samples above 24 h
threshold show anomalies in the running time of one of the generators and
these samples are assigned the label $1$ because it is known that one day only
has 24 hours. The 6010 observations are curated to remove missing samples
leaving 5902 observations with all information. Out of $5905$ observations, it
is observed that $3832$ samples are labelled as normal and $2073$ samples as
abnormal, resulting in 64.8% normal samples and 35.1% abnormalities in the
entire dataset. Figure 6 shows all the clusters and their respective total
fuel consumption including the degree of anomalies in the entire dataset.

###  5.2 Feature importance

Feature selection is performed by fitting the data using a random forest
classifier with 16 features. Note that one could use any other method to find
the most important feature. Since the relative importance of the most
important feature is too high (100% as shown in Figure 7) compared to other
features, any algorithm will predict the same feature as the most important.
Figure 7 shows that the feature “Running time per day” has the greatest
influence on the output and can be coined as the most important feature in the
dataset. Even though it is followed by the “Daily consumption within a period”
the huge difference between the “Running time per day” and the remaining
features shows that priority should be given to the feature “Running time per
day” when considering its reconstruction error for anomaly detection.

Figure 4: Flowchart showing how the labels are decided: For a data sample to
receive the anomaly tag 1, it has to demonstrate at least one of the three
anomaly indicators listed, otherwise it is given the normal tag 0.

###  5.3 Correlation

The correlation matrix is used to visualize the linear relationship between
two variables. The values produced by the covariance matrix range from -1 to
1, where 1 indicates a strong positive linear relationship, -1 indicates a
strong negative linear relationship and 0 indicates that there is no linear
relationship between the variables. Figure 8 shows that the key variable
“Running time per day” and the variable “Daily consumption within a period”
has a strong positive correlation, which is reasonable since the daily
quantity of fuel consumed by a generator is a function of the running time. A
strong correlation is also observed between the three pairs of features;
“Total quantity of fuel after refilling”, “Quantity of fuel found” and the
“Previous quantity of fuel recorded”. However, the feature “Previous hour
meter G1” has no significant correlation with any of the features.

Figure 5: Observed anomaly in the number of working hours in a day for a
single generator. For example, all data samples above 24 h threshold show
anomalies in the running time of the generator and these samples are assigned
the label $1$, because it is known that one day only has 24 hours.

##  6 Results

The results are observed from a number of different evaluation metrics.

###  6.1 Training

In order to assess the performance of our model during training, we set aside
10% of the training data as the validation set. The training and validation
loss are seen in Figure 9 as it trends down different epochs. We observe that
the validation loss is below the training loss. However, this can not
influence the predicted accuracy with further hyper-parameter tuning since the
difference between the two errors is negligible. The training loss accesses
the error in the model during training. Finding the model’s appropriate
threshold requires testing the model with the entire test dataset. The
confusion matrix is used to get the model accuracy and True positive rate
(TPR). They are plotted over a range of thresholds as seen in Figure 10(a) and
the level at which they both attain an average maximum point is the best
threshold for the model. As seen in Figure 10(a), the best threshold is 0.232
with an accuracy score of 0.962. However, if the sensitivity (TPR) or
predicted anomaly is a priority for the organization without minding the cost
they may secure while sorting out the False normal (predicted anomalies that
are actually), the threshold could then be reduced to 0.231. This is the point
where all anomalies are predicted (i.e. predicted anomalies are at 100%) with
a higher false normal (FN) which has an effect on the overall model accuracy.

Figure 6: Fuel consumed per cluster showing the degree of anomalies in the
dataset.

The threshold is also used to categorise these anomalies from mild to extreme
using their reconstruction error. Figure 10(b) shows the reconstruction error.
The samples with extreme threshold are prioritised when one seeks to find the
reasons for the presence of anomalies.

###  6.2 Model performance

The total number of test samples is 1,476 with 1,006 normal samples and the
remaining 470 are anomalous samples. Figure 11 illustrates the performance of
the proposed model at different thresholds based on the confusion matrix. The
proposed model is able to detect a total of 455 anomalous samples correctly
out of the 470 samples with the anomaly label, this accounts for 96.8% (TPR)
of the total anomaly samples. The model also detected a total number of 979
normal samples correctly out of 1006\. The model incorrectly classified 15
normal samples as anomalies (FN) and 27 anomalies as normal samples (FA).
These results show an accuracy of 97.15%, a precision of 94.40%, a recall of
96.81%, a specificity of 97.31% and a F1-score of 95.59%. Table 3 provides a
summary of the performance metrics.

Figure 7: Feature importance for the 16 variables fitted using Random Forest
Classifier. The feature “Running time per day” has the greatest influence on
the output and can be coined as the most important feature in the dataset.

###  6.3 Comparison with other models using the same datatset

Table 4 shows the performance of the proposed model compared to that discussed
in [5]. We used the latter paper to compare our work with because both are
implemented using the same Teneifera dataset for anomaly detection. The
proposed label assisted autoencoder has the best performance with an accuracy
of 97.2% and a recall of 96.8%. The multi-layer perception proposed in [5]
shows the most competitive performance with a higher F1-score, specificity and
precision. However, the label assisted autoencoder is flexible since if we
adjust the threshold, the recall will increase at the cost of specificity and
overall accuracy.

###  6.4 Anomaly classification

The reconstruction error for each data sample differs from one another (as
seen in Figure 10(b)) and provides an opportunity to classify these predicted
anomalies according to their reconstruction error. In this work, 4 classes
$A,B,C$, and $D$ are considered. Class $A$ represents anomalies that are
slightly above the threshold, class $B$ represents anomalies that are above
twice the threshold, class $C$ represents anomalies that are above four times
the threshold, and class D represents anomalies that are above eight times the
threshold. This implies that each class has twice a threshold compared to its
predecessor. Table 5 shows the classes and their corresponding thresholds
showing that 28.25% of the test dataset belongs to the anomaly category of
class $A$, 2.03% belongs to class $B$, 0.20% belongs to class $C$ and 0.34%
belongs to $D$.

###  6.5 Discussion

Feature importance has played an important role in this model, the feature
“Running time per day” has the most significant importance of 100, and “Daily
consumption within a period” coming in the second place with an important
measure of 16, from this we gave priority to the reconstruction error of the
“Running time per day”. To further justify our choice, the correlation matrix
in Figure 8 shows a strong positive correlation of 0.74 between the two
features. Using the reconstruction error of the key variable, we were able to
train and compare the proposed model with the work in [5]. A recall score of
96.8% outperformed all the models proposed in [5] in detecting anomalies.
Also, as shown in Table 3, our model shows a recall score of approximately
100% with a decrease in the overall accuracy of 85% when the threshold is
decreased to 0.231.

Figure 8: Correlation matrix of all numerical features. The “Running time per
day” and “Daily consumption within a period” have a strong positive
correlation, which is reasonable since the daily quantity of fuel consumed by
a generator is a function of the running time.

##  7 Conclusion

Telecommunication industry is one of the dominant information communication
technology industries that rely on a huge amount of electric power supply for
their operations, and thus it is indispensable in their daily dealings.
However, its availability in underdeveloped countries, particularly in Africa,
has been a constant source of contention. Despite the industry’s rise through
the creation of base stations, they have had to turn to alternative energy
sources such as the use of gasoline or diesel with generators, and the use of
solar power, to name a few.

Figure 9: Training and validation loss as a function of the number of
iterations. The mean absolute error is used to measure the loss.

TeleInfra telecommunication company established in Cameroon is one of such
companies hooked on these challenges due to the state of power supply in the
country. The telecommunication equipment that is fixed in different parts of
the rural and urban areas in Cameroon requires an uninterrupted supply of
electricity to achieve the goal of establishing strong and seamless
communication channels in the country, however, the country’s electrical
generation is mostly based on hydropower (73%), with perpetual power
interruptions, particularly during the dry seasons when water levels are low
[36]. The consequence of the diversification to alternative sources of power,
particularly the usage of generators posed another challenge of irregularities
or anomalies in fuel consumption at the base stations due to the observed high
consumption rate in the power generation plants. TeleInfra telecommunication
company is faced with the challenge of unaccounted high fuel consumption for
their operations at the base stations. Since they solely depend on generating
plants as their major source of power supply, they necessarily have to
continually refill these generators and these are done manually. Such
activities are known to have emanated in possible cases of pilferage of fuel
due to the observed anomalies in fuel consumption. As a result, it is
essential to investigate the likely factors contributing to the anomalies by
collecting data on fuel consumption at each of the base stations for the
purpose of minimizing the costs of operation.

We have proposed a label assisted autoencoder-based deep-learning technique
for detecting anomalies in the fuel consumption datasets of the base station
management company namely TeleInfra. In our proposed model, an autoencoder is
used to generate an encoded representation of the input features and construct
the output from the encoded representation to look like the input features of
the series of decoders. The maximum reconstruction error from the trained
model is obtained from the training set and it is set as a threshold for
detecting anomalies on the test dataset. The anomaly detector identifies each
data sample from the testing set as anomaly when they exceed the threshold
assigned. Results showed that our proposed model is highly efficient for
reading anomalies with a detection accuracy of 97.20% and outperforms existing
supervised learning models. The proposed model is flexible and the threshold
is adjustable according to the need of the user which can classify anomalies
from severe to mild.

This work opens future research possibilities which could involve using
different variations of autoencoders such as long short term memory
autoencoders and memory-augmented autoencoders combined with our proposed
label assisted unit. The latter does not require the feature importance
analysis for selecting the best reconstruction error.

(a) Threshold detection

(b) Reconstruction Error

Figure 10: Threshold detection (10(a)) using the reconstruction error (10(b)).
The threshold is used to categorise anomalies from mild to extreme using their
reconstruction error.

(a) Threshold = 0.232

(b) Threshold = 0.231

Figure 11: Detection of results based on confusion matrix. Model performance
---
Threshold |  Accuracy |  F1-Measure |  Recall |  Precision |  Specificity
0.232 |  0.972 |  0.956 |  0.968 |  0.944 |  0.973
0.231 |  0.850 |  0.810 |  1.00 |  0.680 |  0.780
Table 3: Evaluation performance. Model performance
---
Paper |  Techniques |  Accuracy |  F1-Score |  Recall |  Precision |  Specificity
Mulongo et al. [5] |  LR |  0.708 |  0.811 |  0.709 |  0.943 |  0.699
|  SVM |  0.949 |  0.962 |  0.962 |  0962 |  0.925
|  KNN |  0.851 |  0.888 |  0.887 |  0.890 |  0.783
|  MLP |  0.961 |  0.971 |  0.954 |  0.988 |  0.976
Our Model |  AE |  0.972 |  0.956 |  0.968 |  0.944 |  0.973
Table 4: Comparison to similar models using the same dataset. Categorizing
Anomalies
---
Class |  Threshold |  Predicted Number of Samples |  Percentage of Test Data
A |  0.232 |  417 |  28.25%
B |  0.464 |  30 |  2.03%
C |  0.928 |  3 |  0.20%
D |  1.856 |  5 |  0.34%
Table 5: Four categories of anomalies.

## Disclosure Statement

We hereby state that no known competing financial interests or personal ties
could have influenced the research presented in this study.

## References

  * [1] I. Humar, X. Ge, L. Xiang, M. Jo, M. Chen, J. Zhang, Rethinking energy efficiency models of cellular networks with embodied energy, IEEE network 25 (2) (2011) 40–49.
  * [2] J. Lorincz, I. Bule, Renewable energy sources for power supply of base station sites, International Journal of Business Data Communications and Networking (IJBDCN) 9 (3) (2013) 53–74.
  * [3] A. Ayang, P.-S. Ngohe-Ekam, B. Videme, J. Temga, Power consumption: base stations of telecommunication in sahel zone of cameroon: typology based on the power consumption model and energy savings, Journal of Energy 2016.
  * [4] F. J. Espadafor, J. B. Villanueva, M. T. García, Analysis of a diesel generator crankshaft failure, Engineering Failure Analysis 16 (7) (2009) 2333–2341.
  * [5] J. Mulongo, M. Atemkeng, T. Ansah-Narh, R. Rockefeller, G. M. Nguegnang, M. A. Garuti, Anomaly detection in power generation plants using machine learning and neural networks, Applied Artificial Intelligence 34 (1) (2020) 64–79.
  * [6] W. Jinggang, K. Ligai, D. Meixia, Z. Jin, G. Xiaoxia, Feasibility analysis using natural source cooling the idc plant, in: 2009 Chinese Control and Decision Conference, IEEE, 2009, pp. 2579–2584.
  * [7] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * [8] C. C. Aggarwal, An introduction to outlier analysis, in: Outlier analysis, Springer, 2017, pp. 1–34.
  * [9] K. Demestichas, T. Alexakis, N. Peppes, E. Adamopoulou, Comparative analysis of machine learning-based approaches for anomaly detection in vehicular data, Vehicles 3 (2) (2021) 171–186.
  * [10] M. A. Hayes, M. A. Capretz, Contextual anomaly detection framework for big sensor data, Journal of Big Data 2 (1) (2015) 1–22.
  * [11] M. Fahim, A. Sillitti, Anomaly detection, analysis and prediction techniques in iot environment: A systematic literature review, IEEE Access 7 (2019) 81664–81681.
  * [12] H. D. Trinh, E. Zeydan, L. Giupponi, P. Dini, Detecting mobile traffic anomalies through physical control channel fingerprinting: A deep semi-supervised approach, IEEE Access 7 (2019) 152187–152201.
  * [13] S. Omar, A. Ngadi, H. H. Jebur, Machine learning techniques for anomaly detection: an overview, International Journal of Computer Applications 79 (2).
  * [14] S. Vishwanathan, M. N. Murty, Ssvm: a simple svm algorithm, in: Proceedings of the 2002 International Joint Conference on Neural Networks. IJCNN’02 (Cat. No. 02CH37290), Vol. 3, IEEE, 2002, pp. 2393–2398.
  * [15] L. E. Peterson, K-nearest neighbor, Scholarpedia 4 (2) (2009) 1883.
  * [16] T. G. Nick, K. M. Campbell, Logistic regression, Topics in biostatistics (2007) 273–301.
  * [17] D. W. Ruck, S. K. Rogers, M. Kabrisky, Feature selection using a multilayer perceptron, Journal of Neural Network Computing 2 (2) (1990) 40–48.
  * [18] S. Bhattacharyya, S. Jha, K. Tharakunnel, J. C. Westland, Data mining for credit card fraud: A comparative study, Decision Support Systems 50 (2011) 602–613.
  * [19] M. Hasan, M. M. Islam, M. I. I. Zarif, M. Hashem, Attack and anomaly detection in iot sensors in iot sites using machine learning approaches, Internet of Things 7 (2019) 100059.
  * [20] M. Said Elsayed, N.-A. Le-Khac, S. Dev, A. D. Jurcut, Network anomaly detection using lstm based autoencoder, in: Proceedings of the 16th ACM Symposium on QoS and Security for Wireless and Mobile Networks, 2020, pp. 37–45.
  * [21] S. Hawkins, H. He, G. Williams, R. Baxter, Outlier detection using replicator neural networks, in: International Conference on Data Warehousing and Knowledge Discovery, Springer, 2002, pp. 170–180.
  * [22] J. Zhang, H. Zhang, S. Ding, X. Zhang, Power consumption predicting and anomaly detection based on transformer and k-means, Frontiers in Energy Research (2021) 681.
  * [23] G. Münz, S. Li, G. Carle, Traffic anomaly detection using k-means clustering, in: GI/ITG Workshop MMBnet, Vol. 7, 2007, p. 9.
  * [24] C. Zhou, R. C. Paffenroth, Anomaly detection with robust deep autoencoders, in: Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, 2017, pp. 665–674.
  * [25] A. Mishra, S. Krishna Reddy, A. Mittal, H. A. Murthy, A generative model for zero shot learning using conditional variational autoencoders, in: Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2018, pp. 2188–2196.
  * [26] G. R. Johnson, R. M. Donovan-Maiye, M. M. Maleckar, Generative modeling with conditional autoencoders: Building an integrated cell, arXiv preprint arXiv:1705.00092.
  * [27] P. Ge, C.-X. Ren, D.-Q. Dai, J. Feng, S. Yan, Dual adversarial autoencoders for clustering, IEEE transactions on neural networks and learning systems 31 (4) (2019) 1417–1424.
  * [28] X. Guo, X. Liu, E. Zhu, J. Yin, Deep clustering with convolutional autoencoders, in: International conference on neural information processing, Springer, 2017, pp. 373–382.
  * [29] W. Luo, J. Li, J. Yang, W. Xu, J. Zhang, Convolutional sparse autoencoders for image classification, IEEE transactions on neural networks and learning systems 29 (7) (2017) 3289–3294.
  * [30] D. Ferreira, S. Silva, A. Abelha, J. Machado, Recommendation system using autoencoders, Applied Sciences 10 (16) (2020) 5510.
  * [31] M. Sakurada, T. Yairi, Anomaly detection using autoencoders with nonlinear dimensionality reduction, in: Proceedings of the MLSDA 2014 2nd workshop on machine learning for sensory data analysis, 2014, pp. 4–11.
  * [32] W. Wang, Y. Huang, Y. Wang, L. Wang, Generalized autoencoder: A neural network framework for dimensionality reduction, in: Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2014, pp. 490–497.
  * [33] U. Michelucci, An introduction to autoencoders, arXiv preprint arXiv:2201.03898.
  * [34] Y. Wei, J. Jang-Jaccard, W. Xu, F. Sabrina, S. Camtepe, M. Boulic, Lstm-autoencoder based anomaly detection for indoor air quality time series data, IEEE Sensors Journal.
  * [35] D. Bank, N. Koenigstein, R. Giryes, Autoencoders, arXiv preprint arXiv:2003.05991.
  * [36] E. Muh, S. Amara, F. Tabet, Sustainable energy policies in cameroon: A holistic overview, Renewable and Sustainable Energy Reviews 82 (2018) 3420–3429.  doi:https://doi.org/10.1016/j.rser.2017.10.049.

## Appendix A

Architecture of deep neural networks trained in the encoder and decoder. The
number of filters, size of filters and layers are displayed in the encoding
and decoding phases.

Figure 12: Architecture of the deep neural networks trained in the encoder and
decoder.


# Direct and inverse modeling of soft robots by learning a condensed FEM model

Etienne Ménager∗, Tanguy Navez∗, Olivier Goury∗ and Christian Duriez∗ Etienne
Ménager and Tanguy Navez contributed equally to this work. ∗Univ. Lille,
Inria, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Finite Element Method (FEM) is a powerful modeling tool for predicting the
behavior of soft robots. However, its use for control can be difficult for
non-specialists of numerical computation: it requires an optimization of the
computation to make it real-time. In this paper, we propose a learning-based
approach to obtain a compact but sufficiently rich mechanical representation.
Our choice is based on nonlinear compliance data in the actuator/effector
space provided by a condensation of the FEM model. We demonstrate that this
compact model can be learned with a reasonable amount of data and, at the same
time, be very efficient in terms of modeling, since we can deduce the direct
and inverse kinematics of the robot. We also show how to couple some models
learned individually in particular on an example of a gripper composed of two
soft fingers. Other results are shown by comparing the inverse model derived
from the full FEM model and the one from the compact learned version. This
work opens new perspectives, namely for the embedded control of soft robots,
but also for their design. These perspectives are also discussed in the paper.

##  I INTRODUCTION

This last decade, large advances in the modeling and control of soft robots
have been achieved. In the early 2010s, when the soft robotics community
started to experience a very strong growth, it was considered that modeling
and controlling soft robots was a challenge for which no general solution
existed [1][2][3]. Since then, solutions have been proposed that can be
classified into three categories: geometric methods, mechanical models and
learning approaches. Geometric methods are based on strong hypotheses
concerning the behavior of the robot, in particular on hypotheses of constant
or piecewise constant curvatures [4]. They are therefore efficient for some
types of robots, especially continuum robots, and for some applications,
especially when the actuation is well distributed along the main axis of the
robot. The mechanical models are based on the mechanics of continuous media.
Depending on the assumptions applicable to the robot, the work is based on the
theory of beams (non-deformable section) of Euler- Bernoulli [5] and Cosserat
[6], or on the more general theories of 3D deformable solids, with in
particular the use of FEM model [7]. These approaches are more general but
require expertise in numerical modeling, especially to make the computation
fast enough for robotic applications. Machine learning approaches have been
applied to soft robots [12] using the numerous methods and software libraries
that have been developed in recent years. From input-output data, an algorithm
learns and then predicts the robot behavior. Learning approaches sometimes use
mechanical models to generate the needed data, but in this case, the model
simulation is often used as a black box. So far in the literature, mechanical
and learning approaches are relatively distinct. In this paper, we propose, on
the contrary, a ”white box” approach where we take advantage of the structure
of the mechanical model to perform targeted learning.

Contribution: This paper introduces a compact mechanical representation (MR)
of a soft robot based on a FEM model. The compact representation is based on
the data of a condensed compliance matrix of the robot in the actuator /
effector spaces. This matrix varies in a non-linear manner in the workspace of
the robot, but we show that we can learn these data with a standard method,
the MultiLayer Perceptron (MLP) [8]. The method is evaluated on two FEM
simulations of soft robots by testing the obtained inverse models in quasi-
static configuration. We also show how the individually learned models can be
combined through the case of a gripper composed by 2 soft fingers which model
has been learned individually. This method opens interesting perspectives
because it enables to obtain a mechanical model usable for the control and
which is differentiable from data obtained on a finite element model of a soft
robot.

##  II BACKGROUND AND RELATED WORK

###  II-A Soft robot quasi-statics behavior with Finite Element Method

The mechanical behavior of soft robot deformations can be described through
continuum mechanics for which there are no analytical solutions in the general
case. Non linear FEM is one of the numerical methods employed to obtain a
converging approximate solution. We focus on quasi-static formulation without
contacts. Based on solid mechanics, the quasi-static equilibrium is written:

$\displaystyle K(x)\Delta x=f_{ext}-f_{int}(x)+H_{a}^{T}\lambda_{a}$ (1)

where $\Delta x$ are small displacement of nodes around their current position
$x$, $K(x)$ the tangent stiffness matrix depending on the current positions of
the FEM nodes, $f_{ext}$ are external forces such as gravity forces, $f_{int}$
are non linear internal forces of the deformable structure computed from the
material laws chosen for the robot. Actuation constraints
$H_{a}^{T}\lambda_{a}$ are formulated as Lagrangian multipliers. This is
convenient both for forward and inverse modeling.

In the inverse problem, the constraints are not known in advance and we face
quadratic optimization problems solved for the actuation with dedicated
solvers [9] [10]. A soft robot motion is then controlled by specific points
called effectors. The most common use is their control in position where the
squared norm of the shift $\delta_{e}=x_{effector}-x_{goal}$ between effectors
and goals is minimized. The resolution is done in two steps. First, a free
configuration of the robot is solved without any constraint:

$\displaystyle K(x)\Delta x^{free}_{i}=f_{ext}-f_{int}(x_{i-1})$ (2)

After solving this system for $\Delta x^{free}$, it is then possible to
evaluate $\delta_{e}^{free}$ which corresponds to the violation of the
objective if there is no actuation. This violation is updated by a linear
approximation of the nodes displacement using FEM interpolation. Using the
Schur complements gives:

$\displaystyle\delta_{e}=W_{ea}\lambda_{a}+\delta_{e}^{free}$ (3)
$\displaystyle\delta_{a}=W_{aa}\lambda_{a}+\delta_{a}^{free}$

where the $W_{ij}=H_{i}K^{-1}(x)H_{j}^{T}\text{ for }i,j\in{e,a}$ are the
Schur complements [11] i.e. the compliance matrices projected in constraint
space. They gather the mechanical coupling between effectors $e$ and actuators
$a$. The $\delta_{i}\text{ for }i\in{e,a}$ are homogeneous to distances and
computed from the same distance vector $\delta_{i}^{free}$ evaluated at the
previous equilibrium.

Solving the inverse problem finally amounts in solving the following
optimisation problem:

$\displaystyle arg_{\lambda_{a}}$ $\displaystyle min\quad||\delta_{e}||^{2}$
(4) $\displaystyle s.t.$
$\displaystyle\delta_{min}\leq\delta_{a}\leq\delta_{max}\text{ (Actuators
course constraint)}$
$\displaystyle\lambda_{min}\leq\lambda_{a}\leq\lambda_{max}\text{ (Actuation
effort constraint)}$

In the direct problem, when the actuator forces intensity $\Delta\lambda_{a}$
is varying then the corresponding motion of the effectors
$\Delta\delta_{e}=W_{ea}\Delta\lambda_{a}$ can be computed. As actuators are
also controlled in position, varying the position of the actuators
$\Delta\delta_{a}$ lead to obtain the motion created on the effectors through
$\Delta\delta_{e}=W_{ea}W_{aa}^{-1}\Delta\delta_{a}$. This equation provides
the Jacobian $J=W_{ea}W_{aa}^{-1}$ of the soft robot direct kinematics.

This background shows that both direct and inverse models of soft robots can
be derived as long as an approximation of equations 3 given the effective
state of the actuation $\delta_{a}$ is available. In the case of a cable
actuator, this value corresponds to the effective length of the cable.

###  II-B State Representation Learning and robotics

State Representation Learning (SRL) algorithms [15] are applied for learning
small dimension features called states for fully describing an agent. These
methods are widely used in robotics for obtaining a compact representation of
robots when many descriptive sensors are involved. Small dimensional states
are often used for solving control tasks using methods such as Reinforcement
Learning. By keeping only the useful features, better time and computation
efficiency are reached. This is a critical point for this kind of methods,
since testing actions in the environment or retrieving state is expensive.

Mathematically speaking, an agent is considered within an environment with
which it can interacts using actions $a_{t}$ belonging to a given action
space. Each action makes the agent go from a state $s_{t}$ to a state
$s_{t+1}$, unknown but assumed to exist. The agent is not aware of his own
state, but only of observations $o_{t}$ belonging to a given observation
space. The idea behind SRL is to learn a small dimension representation
$\widetilde{s}_{t}$ of the state $s_{t}$ from previous observations. This
amounts to find a mapping $\psi$ such that $\widetilde{s}_{t}=\psi(o_{1:t})$.

In machine learning and more particularly in the robotics and control fields
many methods have been implemented [15]: observation reconstruction such as
Principal Component Analysis [16], learning a direct or inverse model for
example using direct model and encoder [17], using prior knowledge in the cost
function to constrain the state space like robotic prior [18], or integrating
relational biases in data [20] [21]. These methods introduce a compact reduced
model of a robot and tools for learning or reconstructing it. Compared to
these different end-to-end methods, our method enables to learn a mechanical
model that is not specific to a single task.

##  III MATHEMATICAL MODELING

###  III-A Mechanical representation of soft robot

In this work, a high dimensional state of the robot $s_{t}$ is available
through simulation. When using FEM methods, the state of the robot is the
position of all the vertices of its mesh. For reducing the dimension of the
robot state, the global compliance matrix $W$ projected in the constraint
space (see section II-A) is a good candidate because it is a small dimensional
matrix, independent of the size of the FEM mesh and it expresses a direct
relationship between actuators and effectors. In combination with the free
displacement vector $\delta^{free}$, it is relevant enough for modeling and
controlling a soft robot. However, the computation of $W$ is expensive because
it involves the multiplication and inversion of several large matrices and
this needs to be computed at every simulation step. This is an incentive for
learning it rather than computing it from mechanical simulation. Learning
these mechanical quantities rather than non-supervised one makes it possible
to obtain interpretive robot representation useful for simulating and
controlling the robot using mechanical equations.

The learning problem is then resumed as the following: given observations
$o_{t}=(\delta_{a,t},W_{t},\delta^{free}_{t})$ where $\delta_{a,t}$, $W_{t}$
and $\delta^{free}_{t}$ are respectively the effective state reached by the
actuation constraint when applying a specific actuation and corresponding
projected compliance matrix and free displacement of the robot, the goal is to
reconstruct the sub-part of the observation $(W_{t},\delta^{free}_{t})$ from
initial observation $o_{0}$ and current constraint state $\delta_{a,t}$. This
leads to learn a function $F$ such that

$\widetilde{W}_{t},\widetilde{\delta}^{free}_{t}=F(\delta_{a,t},W_{0},\delta^{free}_{0})$
(5)

where $W_{0}$ is the compliance matrix at time $t=0$. This is a good
compliance initialization since it already captures the mechanical coupling
between constraints and the compliance matrix varies smoothly from that
initial position during simulation. This mathematical formulation exclude
cases when the robot may have several configurations for a same actuation i.e.
cases where the simulation falls into a local minimum.

Figure 1: Schematic of the different learning steps involved in the MR
learning framework on the example of a robot with 3 actuator constraints. Step
1: Extraction of the robot constraints and creation of the data. Step 2:
Prediction of $(\widetilde{W}^{tri},\widetilde{\delta}_{t}^{free})$ and
comparison between the prediction and the true value of $W^{tri}$ and
$\delta^{free}$.

###  III-B Learning models

As described in section II-A, both actuation and effector constraints are
considered in this work. The learned compliance matrix $W$ is a symmetric
matrix that represents the relations between the different constraints. It is
formed by the Schur complements introduced in Eq. 3:

$W_{t}=\begin{bmatrix}W_{ee,t}&W_{ea,t}\\\ W_{ae,t}&W_{aa,t}\\\ \end{bmatrix}$
(6)

The magnitude of both the values of $W_{t}$ and $d^{free}_{t}$ depends of the
type of constraints and the measurement units chosen for the simulation. For
avoiding vanishing gradients during the learning process, both the input and
output data are element-wise standardized. This standardization is computed
once using all the training set. The values of $W_{t}$ and $d^{free}_{t}$
related to the actuator space and those related to the effector space can be
of different order of magnitude. In order not to overload the notations, we
consider standardized inputs for learning and not standardized inputs for
control afterward.

As the matrix is symmetric, the networks are trained in a supervised manner
through minimizing the reconstruction loss

$L=\sum_{t}[||\widetilde{W}_{t}^{tri}-W_{t}^{tri}||^{2}+||\widetilde{\delta}_{a,t}^{free}-\delta_{a,t}^{free}||^{2}]$
(7)

where $W_{t}^{tri}$ is the upper triangular part of the matrix $W_{t}$,
$\widetilde{W}_{t}^{tri}$ and $\widetilde{\delta}_{a,t}^{free}$ are the
predicted values of $W_{t}^{tri}$ and $\delta_{a,t}^{free}$ at time $t$ given
by the networks.

Although we also have considered a more complex architecture for the MR
learning described in Section V, a single MLP is sufficient for learning a MR
for control tasks without considering contacts. A MLP is a type of neural
network organized in several layers in which information flows from the input
to the output layer only. In each layer, a learnable affine transformation of
the input is followed by the applications of a non linear function according
to

$X_{n+1}=\sigma_{n}(\omega_{n}X_{n}+b_{n})$ (8)

where $X_{n+1}$ is the output of the layer $n$, $\sigma$ a non linear
activation function and $(\omega_{n},b_{n})$ the weight and bias of the layer
$n$. In our implementation, $X_{0}$ is the raw data vectorized and
concatenated to form a single vector, $\sigma$ is the $ReLu$ function and
$(\omega_{n},b_{n})$ are the learnable weights of the network. The output
$X_{N}$ where $N$ is the number of layers of the MLP is the vectorized couple
$(W_{t},\delta_{t}^{free})$.

The MR learning framework is summarized in Fig. 1.

###  III-C Robot control pipeline from mechanical representation

The soft robot control problem is expressed as a quadratic optimization
problem provided in Eq. 4. The matrices $\widetilde{W}_{ea}$ and
$\widetilde{W}_{aa}$ as well as the vector $\widetilde{\delta}_{a}^{free}$ and
$\widetilde{\delta}_{e}^{free}$ are learned variables. The position of the
goal does not depend on the configuration of the robot in the formulation of
$\delta_{e}^{free}$. Even if $\widetilde{\delta}_{e}^{free}$ is learned for a
specific goal position, this value can then be easily updated for other goals.
For controlling a soft robot using MR, the complete pipeline therefore
consists of an off-line learning phase followed by an online control algorithm
applied on learned variables.

####  III-C1 Off-line learning

: As the framework is a supervised learning framework, pre-computation steps
where data are collected and the learning model is trained are first
performed. The data $(W^{tri}_{t},\delta_{t}^{free},\delta_{a,t})$ are
obtained through simulation and stored in a database. Each data point is
associated to a different robot configuration after applying an action
$a_{t}$. The range of possible values of $a_{t}$ is bounded in a continuous
interval. Sampling the actuation space is done using homogeneous grid search
for building the training dataset and random search for the test dataset. The
number of points of the test set is chosen to amount for $25\%$ of the number
of points of the training set. Once the data is acquired, learning consists in
training the network to predict the quantities
$(W^{tri}_{t},\delta_{t}^{free})$ from the quantities
$(W^{tri}_{0},\delta_{0}^{free},\delta_{a,t})$ using Eq. 8 and 7.

####  III-C2 Online control

: The trained network is used for controlling the soft robot according to the
following scheme:

  1. 1.

$\delta_{a,t}$, $W^{tri}_{0}$ and $\delta_{0}^{free}$ are recovered from FEM
simulation. During the first iteration of a soft robot control experiment,
$\delta_{a,t}$ is obtained for a zero actuation effort applied on the
actuators.

  2. 2.

Prediction of
$(\widetilde{W}^{tri}_{t},\widetilde{\delta}_{a,t}^{free})=F(W^{tri}_{0},\delta_{0}^{free},\delta_{a,t})$
using the learned network. $\widetilde{W}_{t}$ is reconstructed from
$\widetilde{W}_{t}^{tri}$. In order to entirely reconstruct
$\delta_{t}^{free}$ , $\delta_{e,t}^{free}$ is build from the distances
between goals and effectors at time $t$.

  3. 3.

The inverse optimization problem from Eq. 4 is solved using a quadratic
optimization solver. The result is the actuation effort vector $\lambda_{a,t}$
for controlling the robot. In our implementation, the solver from the C++
library qpoases [11] is used, but others can be used.

  4. 4.

The actuation effort $\lambda_{a,t}$ is applied in the simulation and a FEM
time step is computed for updating the actuation state $\delta_{a,t}$.

  5. 5.

The control loop starts again from step 1. until convergence, i.e. until an
accuracy or stability criterion is satisfied.

TABLE I: Data acquisition parameters for each robot. $n_{l}$ sample points resulting in $n_{s}$ samples per cables. Cables displacement are sampled in $[s_{a,min},s_{a,max}]$. Robot | $n_{l}$ | Cables
---|---|---
|  | $s_{a,min}$ | $s_{a,max}$ | $n_{s}$
Diamond | 6561 | 0 | 30 | 11
Finger | 15000 | 0 | 20 | 81

##  IV VALIDATION AND DEMONSTRATION

###  IV-A Robots

The method is illustrated and validated using two FEM models of very different
soft robots driven by cables. Although the method can be applied to other
types of actuation, cable actuators have been chosen for showing that the
unilateral aspect of this type of actuation can be managed. Indeed a cable can
be pulled but not pushed. The two used robots are introduced Fig. 2 and
described below:

  * •

The Diamond robot is a soft parallel silicone robot, driven by 4 cables. The
FEM model of this robot has been validated in [7].

  * •

The Gripper robot uses two identical silicone Fingers, validated in [22]. They
are made up of three segments connected by accordion-shaped joints and
actuated by 2 cables. As this robot is made of two identical Fingers, the
mechanical matrices of the system can be built from the matrices computed on a
single Finger.

These two robots are simulated in SOFA [13] using components from SOFA
plugins. The main interest of these robots is that it is possible to
physically build them and control them from the simulation using an inverse
controller. Simulation results act here as ground truth, allowing both to
easily retrieve training data and to apply control schemes. The transition
from FEM simulation to real robots has already been addressed in above
mentioned works.

For each robot, the data are retrieved from the simulation as described in
III-C. Data acquisition parameters for each robot are summarized in Table I.
In these examples, the actuators considered are cables. During data
acquisition, the space of displacements i.e. the space of relative lengths
$s_{a,t}$ imposed on the cables is discretized and the actuation state
$\delta_{a,t}$ is observed as the reached length of the cables at equilibrium.
Note that many more samples are considered for the Finger than for the
Diamond. Indeed, since more extreme positions and then more deformations are
considered when the Finger closes on itself, more data are required for
capturing larger evolutions of the learned mechanical features.

Figure 2: SOFA simulation of the two robots considered in this study. 3D
simulation meshes and effectors are shown respectively in blue and as black
dots. a) The Diamond robot controlled in position using an effector located on
its tip. b) The Gripper robot controlled using effectors located at the tip of
each finger.

Figure 3: Control results for the Diamond. a) Trajectory of the end effector
in the $(x,y)$ plane for 30 different goals (orange dots) when using learned
mechanical matrices (blue cross) or computed mechanical matrices (green
cross). b) Illustration of the trajectory to be performed by the Diamond end
effector. The end effector of the Diamond is first put in its initial position
represented by a black dot. It must then successively reach different targets
distributed in a circle around this initial position as shown by the orange
dots.

###  IV-B Learning results

A small MLP of 3 layers with 400 weights to learn is sufficient for learning
the MR of each robot. For evaluating the learning process, the average loss is
computed on the test set of random robots configurations. Learning results are
given in Table II.

TABLE II: Initial and final test loss values for Diamond and Finger robots. The best trained network encountered is kept at best epoch. Robot | Initial Loss | Final Loss | Best Epoch
---|---|---|---
Diamond | 8.1e-1 | 1.2e-4 | 6600
Finger | 1.0 | 4.4e-6 | 3730

Although there is more constraints in the Diamond problem, it takes much less
samples data for training the MLP for the Diamond than for the Finger. This is
because the projected compliance matrix of the Finger robot varies much more
than that of the Diamond robot. For the Finger, the influence of the actuators
on the end effector changes a lot depending on the configuration of the robot,
i.e. depending on how bent the robot is. On the contrary, the effector always
has similar behavior in terms of constraints for the Diamond.

###  IV-C Inverse modeling results

In order to evaluate the use of the learned MRs for inverse modeling, a
trajectory is considered in the form of one or several goal positions to be
reached by the effector for each robot.As the soft robot dynamics is neglected
in this study, the results are only valid for trajectories time-scales in the
quasi-static regime.

####  IV-C1 Diamond robot control

A circular trajectory of the effector in a horizontal plane is considered for
the Diamond. This trajectory is chosen for being convenient for two dimensions
visualization although the robot end effector can be controlled in all three
directions of space. Then the positions obtained for each goal using our MR
learning based control pipeline and the classical inverse model are compared.
Control results are shown in Fig. 3.

In most configurations the positions reached by the tips when the compact
learned model is used is very similar to the one when the full FEM model is
considered. This shows that the compact model is a good representation of our
system. On some positions, the error between the learned model and the full
model are a bit larger: this comes from the learning errors specific to the
method. This error can be further decreased by increasing the discretization
of the actuation space during data collection and by judicious choices of
networks hyperparameters. This example validate the hypothesis that we can use
learned quantities $(W_{t},\delta^{free}_{t})$ for performing control tasks.

####  IV-C2 Gripper robot control

Figure 4: Position obtained using mechanical coupling between Fingers $1$ and
$2$ and learned Finger MR. The position of the effector of Finger $1$
$P_{e_{1}}^{t}$ is given by a goal point while the position of the effector of
Finger $2$ $P_{e_{2}}^{t}$ is imposed by the coupling. This coupling is
represented in the picture by an orange object to grasp.

The training done for a single Finger robot, as shown in Fig. 2.b., can be
transferred to the Gripper robot. The learned MR of the two Fingers are
mechanically coupled by adding a force between the effectors. This force
corresponds to a configuration in which the two Fingers catch an object by
enabling one Finger to reach a position while forcing the other Finger to
follow it. To take this coupling into account, the optimization problem from
Eq. 4 is reformulated according to:

$\displaystyle arg_{\lambda_{a,1}^{t},\lambda_{a,2}^{t},\lambda_{e}^{t}}$
$\displaystyle min\quad||P_{e_{1}}^{t}-P_{goal}||^{2}$ (9) $\displaystyle
s.t.$ $\displaystyle P_{e_{1}}^{t}+\beta=P_{e_{2}}^{t}$

where $P_{e_{1}}^{t}$ and $P_{e_{2}}^{t}$ are the respective positions of the
effectors of Fingers $1$ and $2$ at time $t$, $P_{goal}$ the target position
of $P_{e_{1}}^{t}$, $\beta$ a distance describing the size of the object,
$\lambda_{a,1}^{t}$ and $\lambda_{a,2}^{t}$ the forces applied on the cable
actuators of respectively Fingers $1$ and $2$ at time $t$ and
$\lambda_{e}^{t}$ the coupling force shared on the two Fingers. Both the
quantities $P_{e_{1}}^{t}$ and $P_{e_{2}}^{t}$ can be written with the
compliance matrices projected in the constraint space as

$\displaystyle
P_{e_{1}}^{t}=P_{e_{1}}^{t-1}+\widetilde{W}_{ea,1}^{t}(\lambda_{a,1}^{t}-\lambda_{a,1}^{t-1})+\widetilde{W}_{ee,1}^{t}(\lambda_{e}^{t}-\lambda_{e,1}^{t-1})$
(10) $\displaystyle
P_{e_{2}}^{t}=P_{e_{2}}^{t-1}+\widetilde{W}_{ea,2}^{t}(\lambda_{a,2}^{t}-\lambda_{a,2}^{t-1})-\widetilde{W}_{ee,2}^{t}(\lambda_{e}^{t}-\lambda_{e}^{t-1})$

where the $\widetilde{W}_{ea}$ and $\widetilde{W}_{ee}$ values are the
predicted projected compliance values. The quantities at time $t-1$ are known
quantities that are not involved in the optimization.

This optimization problem enables to mechanically couple the two Fingers using
the values learned for a single Finger. Indeed it is possible to build
$\widetilde{W}_{2}$ from the same network used to predict $\widetilde{W}_{1}$
by paying attention to the orientation of the Fingers. Solving this
optimization problem applied to the Gripper robot enables to catch an object.
Starting from the position given in the Fig. 2.b, the resolution of Eq. 9
leads to the result given in the Fig. 4. The movement of the gripped object
between the two Fingers can be controlled. However, there is a limitation when
the force applied to the effectors of the Fingers becomes important for
extreme actuation positions. To solve this problem, it would be necessary to
integrate the possibility of putting forces on the effector during the
learning stage.

##  V DISCUSSION, PERSPECTIVES AND CONCLUSION

In this paper we introduced a method for representing a soft robot in a
compact way. This representation takes the form of MR learning representing
the quantities $W_{t}$ (compliance matrix projected in the constraint space)
and $\delta_{t}^{free}$ (free value of the constraints) given the constraint
state $\delta_{a,t}$ and the initial value $W_{0}$ and $\delta^{free}_{0}$.
The particularity of this method is that these quantities are learned from a
finite set of robot configurations. We showed that this method enables for
controlling robot with quadratic optimisation. Moreover, since prediction is
made using a low-dimensional network, it is very suitable for embedding the
control algorithm on a microprocessor where it is not possible to fully
simulate the FEM model. Although the method has been illustrated on soft cable
robots, nothing prevents it to be used with a wider variety of actuators such
as pneumatic chambers or for rigid robots. In addition, even though the
framework has been introduced with the quasi-static assumption, it can be
extended to soft robot dynamics by taking into account mass and damping
matrices in the expression of $W$.

The model has currently some limitations. First, we do not take into account
the contacts between the environment and the robot. Each contact point
represents a new constraint that must be integrated into the input data. MLP
networks can not handle dynamic resizing of the input and output data. We will
test Graph Neural Network for learning a model independent of the number of
constraints and directly integrate relational biases in its data structure. As
described in [19], message passing could be used to incrementally update the
attributes of the graph on which we learn the matrix. In the MR learning case,
the value $W_{ij}$ in $W$ are the influence of the constraint $i$ on the
constraint $j$. These values therefore represent relationships between
quantities, and it seems natural to us to represent these data in the form of
graphs where the nodes $v_{i}$ represent the constraints defined by its value
$(\delta_{a},\delta^{free})_{i}$ and the edges $e_{ij}=W_{ij}$ are the
interaction between constraint $i$ and constraint $j$. This graph could be
initialized with $W_{0}$ and $\delta^{free}_{0}$, and the update with message
passing to predict $W_{t}$ and $\delta^{free}_{t}$. We believe that this model
will enable handling the contacts in future works.

Secondly, although some research works are interested in creating a reduced
model of a soft robot [16], actuation space scaling is usually the blocking
factor. Indeed collecting data by discretizing the actuation space is
exponential in the number of actuators and becomes very quickly impractical.
For managing cases with more actuators, a discretization of the effector space
could be considered when the number of actuators increase. We also showed with
the example of the Gripper robot that individually learned model could be
combined for building the model of a larger robot. Thus another solution when
we can obtain very distinct modules constituting a robot is to learn
separately their MR.

Finally, this approach has the particularity of providing a differentiable
function $F$ linking the output of the network $W_{t}$ and $\delta_{t}^{free}$
with its input, i.e. with the actuation state $\delta_{a,t}$ of the robot or
the initial compliance matrix projected in the constraint space $W_{0}$. The
creation of a design loss from $W_{t}$ and/or $\delta_{t}^{free}$ would allow
to directly optimize the robot design without the need to simulate it. This is
a critical point in design optimization loop when many designs have to be
evaluated. The creation of a control loss would also allow to find suitable
actuation efforts for a task to be performed. Finally, since the network
represents the whole robot in a compact way and since it enables to predict
the evolution of this representation according to the robot’s actuation state,
it is possible to use it in a model based Reinforcement Learning approach for
learning high level control tasks.

##  VI ACKNOWLEDGMENT

Tanguy Navez would like to acknowledge the support of both the Region Hauts de
France and the ANR through the program AI_PhD@Lille co-funded by the
University of Lille and INRIA Lille - Nord Europe.

## References

  * [1] D. Trivedi, C. Rahn, W. Kier, I. Walker, ”Soft robotics: Biological inspiration, state of the art, and future research” in Applied bionics and biomechanics, vol5 99-117, IOS Press, 2008.
  * [2] H. Lipson, ”Challenges and opportunities for design, simulation, and fabrication of soft robots” in Soft Robotics Journal vol1, p 21-27, Mary Ann Liebert, Inc., 2014.
  * [3] M. Carmel, ”Soft robotics: a perspective—current trends and prospects for the future” in Soft Robotics Journal vol1, p 5-11, Mary Ann Liebert, Inc., 2014.
  * [4] R. Webster III, B. Jones, ”Design and kinematic modeling of constant curvature continuum robots: A review”, in The International Journal of Robotics Research, vol29, p1661-1683, SAGE Publications, 2010.
  * [5] Y.Shapiro, K.Gabor, A. Wolf, ”Modeling a hyperflexible planar bending actuator as an inextensible Euler–Bernoulli beam for use in flexible robots” in Soft Robotics Journal, vol2., p 71-79, Mary Ann Liebert, Inc., 2015.
  * [6] F. Renda, F.Boyer, J. Dias, L.Seneviratne, ”Discrete cosserat approach for multisection soft manipulator dynamics” in IEEE Transactions on Robotics, vol34., p 1518-1533, IEEE, 2018.
  * [7] C.Duriez, ”Control of elastic soft robots based on real-time finite element method”, IEEE international conference on robotics and automation, 2013.
  * [8] S. Hyakin, “Neural Networks: A Comprehensive Foundation.”, 1994.
  * [9] E. Coevoet, A. Escande, and C. Duriez, “Softrobots locomotion and manipulation control using fem simulation and quadratic programming,” in 2019 2nd IEEE International Conference on Soft Robotics (RoboSoft), pp. 739–745,IEEE, 2019.
  * [10] E. Coevoet, A. Escande, and C. Duriez,“Optimization-based inverse model of softrobots with contact handling,” in IEEE Robotics and Automation Letters, vol. PP, pp. 1–1, 022017.
  * [11] H.J. Ferreau, C. Kirches, A. Potschka, H.G. Bock and M. Diehl, ”qpOASES: A parametric active-set algorithm for quadratic programming”, in Mathematical Programming Computation, vol. 6, pp. 327-363, 2014.
  * [12] , D.Kim, S-H. Kim, T. Kim, B. Kang, M. Lee, W. Park, S.Ku, D. Kim, J.Kwon, H. Lee, and others ”Review of machine learning methods in soft robotics”, Plos one, vol 16, 2021
  * [13] F. Faure, C. Duriez, H. Delingette, J. Allard,B. Gilles, S. Marchesseau, H. Talbot, H. Courtecuisse, G. Bousquet, I. Peterlik,et al., “Sofa: A multi-model framework for interactive physical simulation,” in Soft tissue biomechanical modeling for computer assisted surgery, pp. 283–321, Springer, 2012.
  * [14] E. Coevoet, T. Morales-Bieze, F. Largilliere,Z. Zhang, M. Thieffry, M. Sanz-Lopez, B. Carrez, D. Marchal, O. Goury, J. Dequidt,et al.,“Software toolkit for modeling, simulation, and control of soft robots,” in Advanced Robotics,vol. 31, no. 22, pp. 1208–1224, 2017.
  * [15] T. Lesort, N. Diaz-Rodriguez, J.-F. Goudou and D. Filliat, ”State representation learning for control: An overview”, in Neural Networks, 108(2018):370-392, 2018.
  * [16] O. Goury and C. Duriez, ”Fast, generic, and reliable control and simulation of soft robots using model order reduction”, in IEEE Transactions on Robotics, vol. 34, pp. 1565-1576, 2018.
  * [17] J. M. Bern; Y. Schnider; P. Banzet; N. Kumar; S. Coros,“Soft Robot Control With a Learned Differentiable Model”, in IEEE International Conference on Soft Robotics, pp. 417-423, 2020.
  * [18] T. Lesort, M. Suerin, X. Li, N. Diaz-Rodriguez and D. Filliat, ”Unsupervised state representation learning with robotic priors: a robustness benchmark”, in arXiv, 2017.
  * [19] P. Battaglia, J. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, C. Gulcehre, F. Song, A. Ballard, J. Gilmer, G. Dahl, A. Vaswani, K. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. Botvinick, O. Vinyals, Y. Li, R. and Pascanu, ”Relational inductive biases, deep learning, and graph networks”, in arXiv, 2018.
  * [20] J. Kim, J. Park, S. Choi, and S. Ha, ”Learning Robot Structure and Motion Embeddings using Graph Neural Networks”, in arXiv, 2021.
  * [21] J. Almeida, P. Schydlo, A. Dehban, and J. Santos-Victor, ”SENSORIMOTOR GRAPH: Action-Conditioned Graph Neural Network for Learning Robotic Soft Hand Dynamics” in arXiv, 2021.
  * [22] S. Escaida Navarro, S. Nagels, H. Alagi,L.-M. Faller, O. Goury, T. Morales-Bieze, H. Zangl, B. Hein,R. Ramakers, W. Deferme, G. Zheng ans C. Duriez,“ A model-based sensor fusion approach for force and shape estimation in soft robotics”, in IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 5621-5628, 2020.