# Combining multitemporal optical and SAR data for LAI imputation with BiLSTM
network

W. Zhao, F. Yin, H. Ma, Q. Wu, J. Gomez-Dans, P. Lewis  W. Zhao, F. Yin, H.
Ma, Q. Wu, J. Gomez-Dans, P. Lewis
Dept. of Geography
University College London
Gower Street, London, WC1E 6BT UK
{weiying.zhao, feng.yin.15, hongyuan.ma, qingling.wu, j.gomez-dans,
<EMAIL_ADDRESS>

###### Abstract

Leaf area index (LAI) is an important biophysical parameter and plays a
significant role in winter wheat yield prediction. Persistent clouds
dramatically affect the acquisition of crop conditions with Sentinel-2 remote
sensing images, which may cause unreliable yield predictions. Synthetic
Aperture Radar (SAR) can provide all-weather imagery and the ratio between
cross- and co-polarized channels (C-band) has a high correlation with time
series LAI over winter wheat areas. Here, the time series Sentinel-1 VH/VV is
evaluated for imputing LAI, so as to increase its spatial-temporal density. We
propose using a bidirectional LSTM (BiLSTM) network to impute the time series
LAI. The half mean squared error of the predicted time series for each time
step is used as the loss function. During the training, the mean loss over the
observations in the mini-batch is calculated. Two test regions are selected,
one is in the south of German, and another is in the North China Plain. Only
the Sentinel-1 VH/VV and Sentinel-2 generated LAI which is acquired over the
growing season are used during the training. Sentinel-2 generated LAI, which
has been calibrated by in-situ measurements, is treated as the true value
during the training. Plenty of experimental results show that the LAI
imputation results provided by the BiLSTM method are much better than
traditional regression methods, such as exponential function and polynomial
function. It can capture the nonlinear dynamics between multiple time series.
BiLSTM is robust to impute the time series LAI in different winter wheat
fields with different growing conditions. It can even provide satisfactory
results when fewer Sentinel-2 images are available. Since LSTM only see the
information from the past, it cannot provide as good results as BiLSTM,
especially over the senescence period. In conclusion, the results indicate
that the BiLSTM network can be used to impute LAI with time-series Sentinel-1
VH/VV data and Sentinel-2 generated LAI data. This method can also be extended
to solve other time series missing value imputation problems.

_K_ eywords Bidirectional LSTM $\cdot$ imputation $\cdot$ time series $\cdot$
Sentinel-1 VH/VV $\cdot$ Sentinel-2 $\cdot$ LAI

##  1 Introduction

Leaf area index (LAI) can characterize crop canopies with dimensionless
quantity. Timely, continuous and accurate monitoring of crop LAI is critical
for winter wheat yield forecasting. Optical sensors with high spatial and
spectral resolutions have been widely used for generating LAI. However, the
passive acquisition model is easily affected by cloudy weather in the crop
growing period between seedling and flowering, which will lead to lots of
missing values [1]. Currently, more and more available multisensor,
multitemporal and multispectral data allow the estimation of missing values
and fill the temporal monitoring information gaps. SAR sensors are ideal to
fulfil this task, because of their all-time and all-weather capabilities, with
very good accuracy of the acquisition geometry and without effects of
atmospheric constituents for amplitude data. Thanks to the Copernicus
programme, more and more Sentinel-1 and Sentinel-2 data are acquired and free
access. In this paper, we mainly pay attention to the integration of
multitemporal Sentinel-1 Ground Range Detected (GRD) VH and VV ratio (VH/VV)
and limited available Sentinel-2 generated LAI to impute more LAI values.

There are three popularly used methods dealing with missing values: regression
[2], interpolation and matrix completion [3]. Interpolation methods can
provide a smooth result, but it needs the points which can describe the whole
temporal features. For time series Sentinel-2 generated LAI and time series
Sentinel-1, they have the capacity of describing the whole growth of winter
wheat, with fast increasing values in the spring or rainy season and a
decrease in the senescence period [2]. Multiple papers have proved the
temporal correlation between time series vegetation descriptors (LAI and NDVI)
and multitemporal Sentinel-1 data [1, 2, 4, 5].

Even though the winter wheat condition may vary in different fields,
multisensor time series have a high feature correlation in the same
geolocation. With the temporal correlated time series, we transfer the LAI
imputation problem to a multivariate time series regression task. We apply the
BiLSTM network to fulfil this task, it can take into account both the forward
and backward information.

The following section will briefly introduce BiLSTM theory. The preprocessing
of multisensor time series data will be introduced in section 3. In section 4,
we will evaluate the spatial and temporal imputation performance of the
trained BiLSTM network and compare it with the other 3 methods. Then, the
conclusion is drawn in the final section.

##  2 Methodology

BiLSTM is a popularly used non-parametric modelling technique. In the
following, we briefly introduce the general BiLSTM formulations and their
application to solving the time series LAI imputation problem.

---
Figure 1: Bidirectional LSTM architecture. $S$ represents the network state.
(The architecture is modified based on Christopher Olah’ blog1)

The architecture of BiLSTM network 111Neural Networks, Types, and Functional
Programming: http://colah.github.io/posts/2015-09-NN-Types-FP/ is shown as
Fig.1. A BiLSTM layer learns bidirectional long-term dependencies between time
steps of time series. These dependencies can be useful when learning the
complete time series at each time step. Relative insensitivity to gap length
is an advantage of LSTM over RNNs and other sequence learning methods in
numerous applications.

The large variance of the training data values may generate an unstable model,
which will cause a large generalization error. Thus, given a time series
vector $X$, we normalize the training predictors to have zero mean and unit
variance.

$x=\frac{X-\mathbf{E}(X)}{{Var(X)}}$ (1)

For each element in the input sequence, each layer computes the following
function:

$\begin{array}[]{ll}i_{t}=\sigma(W_{ii}x_{t}+b_{ii}+W_{hi}h_{(t-1)}+b_{hi})\\\
f_{t}=\sigma(W_{if}x_{t}+b_{if}+W_{hf}h_{(t-1)}+b_{hf})\\\
g_{t}=\tanh(W_{ig}x_{t}+b_{ig}+W_{hg}h_{(t-1)}+b_{hg})\\\
o_{t}=\sigma(W_{io}x_{t}+b_{io}+W_{ho}h_{(t-1)}+b_{ho})\\\ c_{t}=f_{t}\circ
c_{(t-1)}+i_{t}\circ g_{t}\\\ h_{t}=o_{t}\circ\tanh(c_{t})\\\ \end{array}$ (2)

where $i_{t}$, $f_{t}$, $g_{t}$, $o_{t}$ are the input gate, forget gate, cell
gate, and output gate, respectively. $c_{t}$ is the cell state at time $t$,
$h_{t}$ is the hidden state at time $t$, $x_{t}$ is the input at time $t$,
$\sigma$ is the sigmoid function, and $\circ$ represents the Hadamard product.

Then, the network follows by a fully connected layer
$\tilde{y}_{t}=\sigma(Wh_{t}+b)$, a dropout layer to avoid the overfitting
problem and a fully connected layer with the same length as the input
features, so as to capture all the information of the entire sequence at each
prediction.

The half-mean-squared-error of the predicted responses for each time step is
used as the imputation error:

$\mathcal{L}=\frac{1}{2M}\sum_{i=1}^{M}\sum_{j=1}^{R}(\hat{y}_{ij}-y_{ij})^{2}$
(3)

where $M$ is the sequence length and $R$ is the amount of responses. This
process acts as the regression layer and it will compute the loss over the
training data in the mini-batch. Finally, we can obtain the imputed $\hat{y}$
through minimizing $\mathcal{L}(\hat{y},y)$. The architecture of LSTM only
contains the information flow from the past to the forward (from $S_{0}$ to
$S_{i}$ in Fig.1).

Table 1: Part of multisensor time series with missing values and imputation results. Both the multisensor time series contain the missing values which are indexed by 0. All the valid Sentinel-2 LAI time series are larger than 0, while all the valid Sentinel-1 VH/VV values are smaller than 0. All of these values have the same geolocation. The valid values represent the growing condition of winter wheat. Imputation results | 0.32 | 0.36 | 0.43 | 0.53 | 0.76 | 0.96 | 1.38 | 1.8 | 2.16 | 2.57 | 2.73
---|---|---|---|---|---|---|---|---|---|---|---
Sentinel-2 LAI | 0.32 | 0 | 0.43 | 0.53 | 0 | 0 | 0 | 1.8 | 2.16 | 0 | 0
Sentinel-1 VH/VV | 0 | -8.28 | 0 | -7.85 | -7.56 | -7.21 | -6.84 | 0 | -6.12 | -5.83 | -5.65

During the application of the BiLSTM method, only the temporal points which
have at least one acquisition value are used (Tab.1). In the spatial domain,
the length of the time series may be different. After standardization, the
missing values will be treated as one kind of information. We regard missing
values as variables of the BiLSTM graph, which are involved in the
backpropagation process [6].

##  3 Study area and time series preprocessing

In this section, we introduce the preparation of the multisensor time series.
Two areas are selected to train and evaluate the methods. One is located in
the north of Munich, Germany. The other one is located in Hengshui area, Hebei
Province, China. All the Sentinel images are prepared using Google Earth
Engine. The preprocessing of Sentinel-1 and Sentinel-2 images is the same as
that used in the Sentinel toolbox. The future processing of the data will be
introduced as follows.

###  3.1 Preprocessing of Sentinel-2 LAI time series

Sentinel-2 is systematically acquired at high spatial resolution (10 m) with 5
revisiting days (with 2 satellites) in the same observation angle. Thanks to
Copernicus Programme, it is open access. After atmospheric correction of
Sentinel-2 TOA reflectance, we can retrieve LAI using inverse emulator 222SIAC
on GEE: https://github.com/MarcYin/SIAC_GEE, as shown in [7]. All the time
series LAI over the test areas are generated by Sentinel-2 images.

Only the areas in the Sentinel-2 images which are less affected by the clouds
are selected. 10 Sentinel-2 images are selected which were acquired over the
Hengshui area between March 2017 and June 2017. 12 Sentinel-2 images are used
which were acquired between March 2017 and August 2017 over the North of
Munich.

###  3.2 Preprocessing of Sentinel-1 time series

The Sentinel-1 GRD images which are obtained through interferometric wide
swath mode over the same area during a similar acquisition period as
Sentinel-2 images are collected. To reduce the speckle effect, we use the
multitemporal denoised images [8] to compute Sentinel-1 VH and VV ratio.
Specifically, given a time series of $M$ dB values
$[I_{1}(s),I_{2}(s),I_{3}(s),\cdot\cdot\cdot,I_{M}(s)]$ indexed by time $t$,
the denoised value at location $s$ can be calculated through:

$\hat{J}_{t}(s)=\frac{\mu^{{ML}}_{t}(s)}{M}\sum_{t^{\prime}=1}^{M}\frac{\emph{I}_{t^{\prime}}(s)}{\mu^{{ML}}_{t^{\prime}}(s)}$
(4)

where $\emph{I}_{t}(s)$ is the calibrated dB value, $J_{t}(s)$ the noise-free
value, $\mu^{{ML}}_{t}(s)$ is the initial estimation value using multi looking
(with a window size of 14$\times$14), $s$ is the position in one image and $t$
is the serial number of $M$ images. Then, the variance of the denoised value
can be computed with:

${Var}_{t}(s)=\hat{J}^{{ML}}_{t}(s)\frac{M+N-1}{MNL}$ (5)

Finally, we can compute the dual polarization ratio with is computed through
${\hat{J}^{{VH}}_{t}(s)}/{\hat{J}^{{VV}}_{t}(s)}$ [2].

To ensure all the time series SAR data $I_{t}(s)$ capture the same signal from
the object, we only use the images which are acquired through a similar
incidence angle. 214 images (orbit 117) are used for multitemporal SAR image
denoising over the Munich test area, while more than 220 images acquired
through ascending orbits over the Hengshui test area are used. All the SAR
images are well-registered with the Sentinel-2 images which are obtained at
the same place.

##  4 Experimental results and discussion

In this section, we demonstrate the spatial and temporal imputation
performance of the BiLSTM network by comparing it with some popularly used
methods.

In the test areas, most of the Sentinel-2 LAI and Sentinel-1 VH/VV values are
acquired at different times. To keep the spatial resolution of LAI, polynomial
and exponential regression parameters are computed once for each time series.
During the imputation, Sentinel-1 VH/VV time series are interpolated according
to Sentinel-2 acquisition time. Since non-LSTM methods cannot handle arbitrary
lengths of the time series, we only compare them with the test of real-time
series.

LSTM and BiLSTM have the same parameters before training. 60 hidden units are
used, followed by a fully connected layer of size 50, a dropout layer with a
0.5 probability of dropping out input elements and a fully connected layer
with the same length as the input features. The depth will be good for well-
leaning longer-term dependencies, that is partly the reason people stack LSTMs
in sequence-to-sequence models. The Sentinel acquisition frequency is very
high over Germany test areas, they can show all the temporal dynamic features
of winter wheat. We will use them to train the networks. During the whole
growing period, the length of the time series may be different in different
places, either caused by the weather or the sensor’s acquisition settings. To
ensure the trained network is robust to different situations, the training
time series has been randomly sampled with the shortest length equal to 8 and
the longest length equal to 69. Only the pixels which correspond to winter
wheat are used.

###  4.1 Temporal imputation

4 different time series (Fig.2) are selected to show the imputation difference
provided by the different methods. There are some gaps in the LAI time series
between April and May, which correspond to the rainy season. Sentinel-2 images
cannot capture the crop growing conditions because of the cloudy weather.
However, the LAI values acquired over this period are very important for crop
yield prediction.

For the fourth test area (fourth row of Fig.2), the saturation of SAR
backscattering values may be caused by the similar side-looking angle and
sowing direction. This similarity may cause SAR sensors to capture low-volume
backscattering values or by the large multi-looking window size during the
multitemporal single-channel denoising. This phenomenon popularly exists in
the North China Plain, especially when all the ascending orbits acquiring data
are used for the multitemporal denoising.

Figure 2: Different time series LAI imputation results in comparison. From
left to right: polynomial regression, exponential regression, LSTM, and
BiLSTM. Imputed results are shown using the red line. Four different places
(four rows) in the North China Plain are used. The imputation results and
index using the red line. Sentinel-1 VH/VV time series are smoothed using a
Gaussian kernel with a window size equal to 5. All the SAR data in the first
three rows are acquired through the same orbit, while the SAR data in the
fourth row are acquired through the ascending track with an orbit number equal
to 40 or 142.

When the multisensor time series have similar temporal features, polynomial
and exponential imputation methods can provide satisfying results. These two
methods are much simpler and faster than LSTM and BiLSTM. However, when the
two features have a complex relationship, they can not provide good results,
as shown in the fourth row of Fig.2. When the training samples contain too
much noise. The multisensor time series have a nonlinear relationship which is
too hard to describe using a low-degree regression model.

Table 2: LAI imputation performance using different methods. 20 points are used for the evaluation. The 20 points are distributed in the Hengshui area with two times field measurements. One is in March, another is in May.  Methods | RMSE (March) | RMSE (May)
---|---|---
BiLSTM | 0.68 | 1.7855
LSTM | 0.5912 | 2.2886
Polynomial | 0.9339 | 2.4088
Exponential | 0.9361 | 2.624

Generally, people always compute the regression coefficients between one
feature and another, and find the optimized parameter based on the least-
squares estimation. This will lead to the inverse results smaller than the
referenced target feature in the peak part.

LSTM and BiLSTM can take into account the two features and generate better
imputation results. As we can see from Fig.2, all the imputation results
provided by LSTM and BiLSTM are better than the former two methods. Even when
Sentinel-1 signals have complex changes, they still provide good results.
During the imputation, they can take the limited number of LAI points as the
weight to adjust the imputed time series.

###  4.2 Spatial performance of the imputation

All the imputation is performed in the temporal domain for each geolocated
point. In this section, we will analyze the imputation spatial performance by
comparing the results with Sentinel-2 provided LAI image and some in-situ
measurements. During the imputation, the parameters of polynomial and
exponential regression methods are separately computed for each time series.

As observed in Tab.2, both BiLSTM and LSTM obtain smaller root-mean-square
error (RMSE) results. With the reference of Fig 2, the big RMSE results
provided by LSTM may be caused by the larger imputed LAI values in May. The
big RMSE values provided by polynomial and exponential regression methods are
caused by the lower imputation results.

|  |
---|---|---
(1) Sentinel-2 RGB image | (2) Polynomial regression | (3) Exponential regression
|  |
(4) Sentinel-2 LAI | (5) LSTM network | (6) BiLSTM network
0  6
Figure 3: Different LAI imputation results in comparison. The acquisition time
of the Sentinel time series is shown at the bottom of the figure, with *
representing Sentinel-2 acquisition time and $\times$ representing Sentinel-1
acquisition time. Sentinel-2 image (1, 4) is acquired on 18/04/2017 (red *)
over North China Plain, and Sentinel-1 image is acquired on 19/04/2017 (blue
$\times$). LAI values range from 0 to 6.

From Fig.3 we can see that LSTM and BiLSTM can provide much better imputation
results in the spatial domain than polynomial and exponential regression
methods. Their imputation results have better spatial features, with smooth
values in the same field and clear line features in the boundary areas.
Compared with LSTM-provided results (Fig.3(5)), BiLSTM can also provide better
results in some narrow fields (Northeast of the test area in Fig.3). The
spatial distribution of the results provided by BiLSTM is much more similar to
the real LAI.

##  5 Conclusion

In this paper, we successfully applied the BiLSTM network to impute LAI with
multivariate time series data. Plenty of experiments show the effectiveness of
its imputation over the spatial and temporal domains. It is obvious that both
BiLSTM and LSTM imputation performance outperforms traditional regression
methods. It can well impute the missing values in the time series, especially
during the fast-growing season of winter wheat. Even when the temporal
Sentinel-1 VH/VV values have a complex change, the BiLSTM network still could
provide effective imputation results. In the future, we will take the spatial
correlation information into account for missing value imputation. The
spatiotemporal multisensor LAI imputation process has been separated into
several steps, we will find ways to merge them.

## References

  * [1] L. Pipia, J. Muñoz-Marí, E. Amin, S. Belda, G. Camps-Valls, and J. Verrelst.  Fusing optical and SAR time series for LAI gap fillingwith multioutput Gaussian processes.  Remote Sensing of Environment, 235:111452, 2019.
  * [2] A. Veloso, S. Mermoz, A. Bouvet, T. Le Toan, M. Planells, J.F. Dejoux, and E. Ceschia.  Understanding the temporal behavior of crops using sentinel-1 and sentinel-2-like data for agricultural applications.  Remote Sensing of Environment, 199:415–426, 2017.
  * [3] J. Yoon, W.R. Zame, and M. van der Schaar.  Estimating missing data in temporal data streams using multi-directional recurrent neural networks.  IEEE Transactions on Biomedical Engineering, 66(5):1477–1490, 2018\.
  * [4] Y. Kim, T. Jackson, R. Bindlish, H. Lee, and S. Hong.  Radar vegetation index for estimating the vegetation water content of rice and soybean.  IEEE Geoscience and Remote Sensing Letters, 9(4):564–568, 2011\.
  * [5] L. Stendardi, S.R. Karlsen, G. Niedrist, R. Gerdol, M. Zebisch, M. Rossi, and C. Notarnicola.  Exploiting time series of Sentinel-1 and Sentinel-2 imagery to detect meadow phenology in mountain regions.  Remote Sensing, 11(5):542, 2019.
  * [6] W. Cao, D. Wang, J. Li, H. Zhou, L. Li, and Y. Li.  BRITS: Bidirectional recurrent imputation for time series.  In Advances in Neural Information Processing Systems, pages 6775–6785, 2018.
  * [7] F. Yin, P. E. Lewis, J. Gomez-Dans, and Q. Wu.  A sensor-invariant atmospheric correction method: application to Sentinel-2/MSI and Landsat 8/OLI.  https://doi.org/10.31223/osf.io/ps957.
  * [8] S. Quegan and J.J. Yu.  Filtering of multichannel SAR images.  IEEE Transactions on Geoscience and Remote Sensing, 39(11):2373–2379, 2001.


# Learning a General Model of Single Phase Flow in Complex 3D Porous Media

Javier E. Santos
Los Alamos National Laboratory
Agnese Marcato
Los Alamos National Laboratory
Qinjun Kang
Los Alamos National Laboratory
Mohamed Mehana
Los Alamos National Laboratory
Daniel O’Malley
Los Alamos National Laboratory
Hari Viswanathan
Los Alamos National Laboratory
Nicholas Lubbers
Los Alamos National Laboratory

###### Abstract

Modeling effective transport properties of 3D porous media, such as
permeability, at multiple scales is challenging as a result of the combined
complexity of the pore structures and fluid physics - in particular,
confinement effects which vary across the nanoscale to the microscale. While
numerical simulation is possible, the computational cost is prohibitive for
realistic domains, which are large and complex. Although machine learning
models have been proposed to circumvent simulation, none so far has
simultaneously accounted for heterogeneous 3D structures, fluid confinement
effects, and multiple simulation resolutions. By utilizing numerous computer
science techniques to improve the scalability of training, we have for the
first time developed a general flow model that accounts for the pore-structure
and corresponding physical phenomena at scales from Angstrom to the
micrometer. Using synthetic computational domains for training, our machine
learning model exhibits strong performance (R2=0.9) when tested on extremely
diverse real domains at multiple scales.

##  1 Introduction

Estimating the transport properties of materials is challenging. On one hand,
laboratory experiments are costly and require specialized facilities and in-
hand processed samples that can not be re-used. On the other hand, numerical
simulations often require complex meshing, assumptions and approximations, and
are computationally intensive, especially for large complex domains.
Alternatively, machine learning (ML) models are demonstrated to learn complex
relations from large datasets with ease [1]. Most importantly, thanks to
advances in software and hardware, ML models are able to provide predictions
in new domains in milliseconds. This makes them attractive for forecasting
physical processes, especially in the geosciences, where real-time predictions
can be critical in field operations. However, training and testing machine
learning models with subsurface data has many challenges. First, ML models are
commonly optimized to 2D image benchmarks [2, 3]. On the other hand, imaging
technologies (CT, FIB-SEM, MRI) have had breakthroughs [4] that allow us to
obtain high definition samples that can be of size 2563 or more; 3D domains
are much more computationally challenging than most 2D images. Secondly,
subsurface formations exhibit a wide range of complex structures due to the
diverse rock genesis and diagenesis processes [5]. Finally, it is not
straightforward to embed the complexity of varying degrees of confinement
effects in the ML model. While much progress has been made in emulating
transport in porous media using ML [6, 7, 8, 9, 10, 11, 12, 13, 14], models in
the literature are not general enough to be used as a generic tool.

Fluid transport driven by an external force, for example, a pressure gradient
or the gravitational force, is affected by viscous, slip, and diffusive
effects depending on the characteristic length-scale of the available flow
paths. The primary way to quantify the impact of these effects on the ability
of a fluid to travel through a sample is the apparent permeability. This
quantity provides a volume-average scalar that summarizes all the micro- and
nano-scale effects into a single number, which can then be used in continuum
simulators or analytical approximations for forecasting field-scale processes.
A strict lower bound of the apparent permeability is the absolute
permeability, defined as the limit of the apparent permeability when
nanoconfinement effects vanish. Slip and diffusive effects can enhance
apparent permeability up to three orders of magnitude compared to absolute
permeability, depending on the geometry and pressure (as described in [15] and
shown in Figure 2).

The structure and connectivity of the pore structures have the leading impact
on the apparent permeability of a 3D porous medium. In micron-sized pores,
fluid transport is approximated by Darcy’s Law where viscous forces are
dominant with negligible contributions from inertial forces. During viscous
flow, the velocity profile of the fluid is non-uniform and affected by the
pore morphology and pressure gradient. This non-uniformity stems from the
viscous drag forces which are uneven across the flow paths. Fluids prefer the
least-drag force paths. For a tube with a circular cross-section, the viscous
forces are minimal at the center of the pore given rise to a parabolic-shaped
velocity profiles. In contrast, in smaller pores, transport is affected by
diffusion and slip at the solid-fluid boundary since the mean free path length
of the molecule is comparable to or larger than the pore size, increasing
molecular collisions [16] and enhancing permeability [17, 18, 19]. In nature,
permeable media hosts conduits and heterogeneities in a wide range of length-
scales. Subsurface formations have pores that can be a few Angstroms wide [20,
21] to fractures with apertures on the centimeter scale [22]. These wide range
of scales encompass flow regimes going from free-molecular to viscous
dominated [16]. These multiscale heterogeneities usually co-exist. For
example, in fractured shale a fluid can flow much faster through the fracture,
but the nanoporous matrix has a non-negligible contribution as shown by [23].

Laboratory experiments can be set-up to study deviations from Darcy’s law.
However, since the resolution of the instruments is on the same scale as the
phenomenon under investigation, running an experiment at subsurface conditions
is very challenging. In practice, samples are subject to extremely high
pressure gradients [24, 25, 26, 27, 28] but this does not reflect the
conditions that subsurface flow modellers are interested in. On the
computational side, analytical solutions considering fluid flow inside a tube
have been derived either from empirical observations or classical fluid
dynamics and can be used with pore-network models [29]. Although these pore-
network models try to capture the 3D connectivity of real samples, they
aggressively oversimplify the heterogeneity of the pore space by fitting it
with spheres and tubes. The high aspect ratio of the tubes will significantly
influence the solution of non-Darcy flows and as seen in SEM images of
mudrocks, pore structures vary from “spongey” kerogen patches, to long and
skinny microfractures, to stacks of platelets [30]. Finally, molecular
dynamics simulations can represent nanoconfinement effects accurately.
Nonetheless, the computational resources needed are vast, which limits the
calculations to very small volumes (boxes with side lengths of around 10 nm).

The lattice-Boltzmann method (LBM) can operate in larger domains, which yield
in representative sample sizes since it operates at a mesoscale. LBM
integrates physical insights from both the molecular- and micron-scales and
can predict larger-scale properties. In short, fluids are simulated as swarms
of particles, represented by particle distribution functions, which flow on a
discrete lattice. Landry [20] proposed the local effective viscosity model
(LEV-LBM) to simulate multiple scales encompassing multiple flow regimes by
utilizing a spatially-varying mean free path (sMFP). The sMFP accounts for the
reduction of the local kinetic viscosity caused by the confinement effect at
the nanoscale. One thing that sets apart the LBM simulation of that at the
nanoscale, is that flow fields do not rely exclusively on the pore-structure
but also on length scale, confinement pressure, temperature, and gas
properties. The LEV-LBM also accounts for slip and diffusive effects depending
on the scale. The main downside is that simulating a large complex sample with
confinement effects could take weeks running on 1,000 processors.

While the LBM offers detailed insights into fluid dynamics, there is a growing
interest in leveraging this data and the computational efficiency of neural
networks for permeability modeling. Neural networks are an attractive avenue
for modeling permeability since they can provide predictions in a split-second
enabling real-time field applications. Therefore, there have been many
attempts to train neural networks to predict permeability from images. In a
pioneering study, [31] employed a convolutional network (ConvNet) to learn a
mapping between a 2D image and its permeability, since then ConvNets have
remained the model of choice for this task. [32] provides a comprehensive
review of notable works to this date. Nevertheless, computational scaling with
domain size is still one of the main hurdles in training models with 3D images
of permeable media. Recently, we provided a proof of concept of a network that
allowed training and testing with large 3D images on a single graphics
processing unit, capturing representative elementary volumes within the
weights of the network [33].

In this letter, we present a general model for fluid transport through
permeable media building upon the work of MS-Net [33] by improving the model
design, scaling-up the parameter count (model size), and increasing the size
and diversity of the training set. For the first time, we used hundreds of
computationally large flow simulations in synthetic samples to train our
model, including a wide range of pore sizes, rock types, and boundary
conditions. The comprehensive dataset used for all the training and testing
data are publicly available [34]. To accommodate the larger model and dataset,
we likewise applied a number of computer science tools to ameliorate the
increased computational cost, which we detail. Our contributions are divided
into two main categories, 1) model architecture, where we the describe in
detail the model that allow for multi-scale processes to be accounted for and
2) training, which is essential for distilling information from such a diverse
set of samples properly and its vital for managing a model of this size and
complexity. This workflow is general and transferable to other properties and
domains. Finally, we show that the trained model is able to generalize the
problem of fluid flow through permeable media for more than one thousand
simulations, including many microCT images of real porous media.

Figure 1: Workflow figure.

##  2 Machine learning workflow

We describe the machine learning workflow in the following subsections.
Dataset describes the data used for both training and testing. Model
architecture describes the functional form of the machine learning model.
Training algorithms explains the various techniques required to scale the
training to large datasets and complex geometries, including the training
schedule which describes the actual train procedure. An schematic showing all
the stages of this process is shown in Figure 1.

###  2.1 Dataset

For a thorough review of the dataset, we refer the reader to [34]. The DRP-372
dataset includes more than 1200 flow simulations from a wide variety of
samples hosted in the Digital Rocks Portal [35]. We categorized the dataset
based on 6 lithology groups: spherepacks, processed-based geometries, shales,
carbonates, sandstones, and single fractures. The flow simulations were
carried-out in computational sizes of 2563 and 4803 where some of the latter
took up to 4 days running on 1,000 CPU cores to complete. The resolution of
these samples spanned from 0.5 nanometers to 5 micrometers. These simulations
were carried-out using a local-effective viscosity lattice-Boltzmann model
[20] which captures the impact of nanoconfinement effects, slip boundary
conditions, and diffuse boundary layers when applicable. The confinement
pressure was varied among individual simulations yielding a dataset that
captures a wide range of scales and flow regimes (Figure 2). The lower
pressure bound is set at 1 MegaPascal (MPa) where significant nanoconfinement
effects are present. The upper pressure was 20 MPa, where negligible
confinement effects were shown and the solution is the Navier-Stokes equation
solution. These simulations encompass flow regimes from the free molecular to
viscous flow. The data spans 9 orders of magnitude of apparent permeability.
The computational cost of these simulations required more than 5,000,000 core-
hours.

To train our model, we utilized 60 simple synthetic geometries described in
[36], including various types of sphere packs, both homogeneous and
heterogeneous, and fractures. These process-based samples are designed to
emulate various lithologies. These synthetic samples feature well-resolved
pores and interconnected structures in the flow direction. The main motivation
to use synthetic samples is to show that our model is able to abstract how
fluid flow behaves under different physical conditions and is able to be used
in a wide range of unseen real geometries.

Consequently, during testing, we evaluate the model’s performance
predominantly on real porous media, thereby highlighting its capacity to
transfer knowledge from synthetic to real-world scenarios. This dataset
includes 75 diverse samples from research groups worldwide [37, 38, 39, 40,
41, 42], captured using various imaging devices or simulated using different
process-based codes. By leveraging this extensive dataset, our study
rigorously evaluates the generalization capabilities of our model across a
wide range of real-world samples.

Figure 2: Flow regimes for simplified samples, in reality, these effects are
magnified due to the complexity of the 3D structures, for example, how
tortuous the paths are of how narrow some paths can get. Figure 3: The bottom
panel shows a cross-section of the nanoconfined simulation at five different
pressures. The normalized velocity (with respect to the maximum value)
contours are shown. It is visible that the sample at the lowest pressure (1
MPa) presents the higher degree of slip, shown as “fat“ velocity profile
(fewer contour lines), while the one at 20 MPa exhibits no slip. This can also
be seen in the top panels, where normalized mean free path (that could be
interpreted as the deviation from the Navier-Stokes solution) is shown.

###  2.2 Model architecture

Traditionally, analytical solutions have been proposed for permeability
estimation using average (effective) structural properties, the most prominent
example is the use of the porosity. Nevertheless, the average porosity of a
sample is a non-unique proxy (many samples have the same average porosity
value) and fails to capture the sample’s pore-space complexity and spatial
distribution. Another notable example is the fracture width used in the cubic
law to estimate single-fracture permeability. To get a more complete picture,
our proposed ML model uses the whole complex 3D structure, with its
thermodynamic conditions and scale. These are all included to aid the network
to compute the resultant apparent permeability under these conditions. We
explain in detail our network architecture in the following subsection.

We build our model using a larger version of the MS-Net [33] architecture,
which demonstrated a remarkable ability to learn from large porous media
samples. In short, the MS-Net is a system of convolutional networks that
process data at different resolutions. The system is trained end-to-end by
passing the outputs of the coarser networks sequentially as inputs to the
finer ones 111The details about the coarsening and refinement procedure can be
found in Section 2.4 of [33].. There are two main motivating factors for using
this architecture: 1) This is an elegant way to process large ($>$2563)
geometries without subsampling the domain, and 2) the coarsest model can see
the entire image at once and predict the bulk of the flow while the subsequent
ones refine it. In this work, we have implemented the following enhancements
to the previous model:

  * •

Input features: To capture simulations encompassing diverse thermodynamic
conditions, we introduce a spatially-varying 3D normalized mean free path in
our model. This feature, which arises from the Peng-Robinson equation,
accounts for the reduction of the local kinetic viscosity caused by the
confinement effects at the nanoscale (the detailed explanation can be found in
[34])222Initially, we tried training our network with only this input. The
trained model provided accurate predictions of permeability but the flow was
not spatially-smooth, since this map does not have enough variance, as shown
in Figure 3. We also include the 3D Euclidean distance of the pore-space. The
Euclidean distance labels each voxel with the distance to the closest solid
boundary.333The Euclidean distance outputs large values in the middle of the
pores and smaller values near the solid boundaries, is amenable to the
cylindrical profile of a fluid flowing on laminar flow. This feature aids to
get a spatially-smooth prediction. The computational cost of these two
features is negligible, even in very large domains.

  * •

Large number of scales: To significantly expand the model’s field of vision,
we augmented our system by incorporating six convolutional sub-networks. This
design decision enables the coarsest network (NN6 in Figure 1) to receive an
input that is $\frac{1}{32}$th the size of the original domain. Despite the
downsizing, the input retains crucial global structural features that affect
the fluid flow, which allows the kernels in the sixth network to see most of
the image at once.

  * •

Padding: To facilitate the model’s comprehension of flow continuity along the
flow direction, we incorporated periodic padding for the 3D convolution
operations. Meanwhile, since the faces perpendicular to this axis feature no-
flow boundaries in the simulation, we assigned zero-padding to these regions.
These design choices are particularly valuable in the networks operating at
coarser scales, as their field of vision is comparable to the input size (At
the coarsest scale, a 2563 domain corresponds has an input size of 83).

  * •

No flow in the solid-space: To explicitly ensure no flow within the solid
space, we introduced an additional constraint at the end of our model. The
final output of the model is multiplied by the input binary image, which
contains ones in the pore-space and zeros in the solid-space. This
multiplication effectively disregards the gradients of voxels outside the
pores. By significantly reducing the search-space, this approach leads to
notably faster training. Additionally, we conducted experiments where we
allowed gradients of the first solid voxel surrounding the grains, but
observed no discernible difference in the training process.

  * •

Fewer activation functions and normalization layers: To enhance accuracy, we
implemented a strategy recommended in [43] by reducing the number of
activation functions and normalization layers in our model. Specifically, we
removed the normalization layer preceding the last convolutional layer, as
well as the activation layer following it. This adjustment resulted in an
improvement in accuracy.

  * •

GeLU activation: To align with recent advancements in models, we transitioned
our remaining activations to Gaussian Error Linear Units (GeLU) [44]. GeLU
serves as a smoother alternative to ReLU and has gained popularity in the
field. By adopting GeLU, we achieved a slight increase in accuracy while
maintaining the same computational overhead.

  * •

Number of filters at each layer: Initially, we experimented with increasing
the depth of each network to enhance model capacity. However, this approach
proved challenging to train, likely due to increasingly difficult gradient
paths. Consequently, we opted for a substantial increase in the number of
parameters, significantly raising the count to 400M compared to the 10k of the
MS-Net. This substantial capacity boost resulted in a significant improvement
in accuracy. Notably, the augmented model’s training process required
adjustments to the training algorithm, as elaborated in Section 2.3, to
effectively handle the increased memory requirements.

In addition, we explored the inclusion of attention layers in the coarser
networks. While this design choice led to slight improvements in accuracy, it
also incurred significant computational costs. Furthermore, incorporating
attention mechanisms posed challenges in accommodating different image sizes,
making it less practical for our purposes. As a result, we opted for a fully
convolutional architecture, which better addressed the need for flexibility in
working with diverse image sizes.

###  2.3 Training algorithms

The development of an enhanced model, encompassing larger networks with
increased parameter counts, and training it on computationally larger domains
with more features, has yielded a notable improvement in accuracy for both
training and testing datasets. However, scaling up the model presented
significant challenges, particularly in addressing the memory constraints
imposed by the GPUs. In this section, we emphasize the algorithms that were
crucial in overcoming these obstacles and enabling the successful
implementation of the model.

  * •

Mixed-precision training: In order to accommodate the larger model and
facilitate training with bigger domains, we carried out the training process
using 16-bit precision. This approach offers several advantages, including
reduced memory consumption as the network parameters occupy half the space,
and accelerated data transfer due to operations requiring less bandwidth and
executing at a faster pace. However, it is important to note that training
stability can be compromised due to precision loss, which poses a significant
challenge in our specific case.

  * •

Activation checkpointing: To reduce GPU memory usage, we implemented
activation checkpointing in our model. This technique avoids the storage of
layer activations, recomputing them during the backward pass instead. By
sacrificing some computation time, we achieve a reduction in memory
requirements, allowing us to train larger models with the same memory
footprint.

  * •

Gradient accumulation: Even with the utilization of 16-bit precision, the size
of our model and samples restricts us from performing a backward pass with
more than one domain at a time. Since training with a single sample per step
(online training) would be extremely noisy, we use gradient accumulation to
“simulate” bigger batches. The idea behind gradient accumulation is simple:
after calculating the loss and gradients for each sample, instead of
immediately updating the model parameters, the gradients are accumulated
across consecutive samples. The model parameters are then updated based on the
cumulative gradient after a specified number of samples. This approach serves
the same purpose as having a multiple sample mini-batch, all while maintaining
the memory footprint of a single sample.

  * •

Loss function: Our dataset encompasses a wide range of scales, with
permeabilities spanning 8 orders of magnitude and individual pixel-wise
velocities spanning 24. To effectively convey the relative importance of each
sample to the model, we employed a loss function that calculates the pixel-
wise error relative to the variance of each velocity field. This approach
ensures that every sample carries the same level of "importance" during
training.

$\mathcal{L}=\sum_{i=1}^{g_{steps}}\sum_{s=1}^{6}\frac{{\rm
sum}((u_{i,s}-\hat{u}_{i,s})^{2})}{N_{vox,s}\cdot\sigma(u_{i})^{2}}.$ (1)

$g_{steps}$ are the number of gradient accumulation steps, $u_{i}$ represent
each of the 6 coarsened velocity fields, $N_{vox,s}$ the number of voxels per
scale and $\sigma(u_{i})^{2}$ the standard deviation of the velocity field.
Our loss function tries to reconstruct this. Solely computing the mean squared
error in the last layer would assign greater importance to samples with higher
velocities.

  * •

Input and output normalization: To enhance training speed and stability, we
normalized the Euclidean distance with a constant value of 50 and the velocity
with a value of 1e-9. These value choices accelerated the training process and
enhanced its stability. Additionally, we kept the normalized mean free path in
its original range of zero to one without further normalization.

  * •

Gradient clipping: To mitigate potential instabilities during training, we
employed gradient clipping, a technique that limits the magnitude of gradients
to a specified value. By rescaling the error derivative, the weight updates
are also appropriately adjusted, minimizing the risk of overflow or underflow.
This approach ensures that gradient descent exhibits stable behavior, even
when faced with a complex and irregular loss landscape. By implementing
gradient clipping, we successfully prevent training instabilities from
adversely affecting the model’s training path.

  * •

Optimizer: To enhance training stability and improve generalization to unseen
data, we utilized the AdamW optimizer during model training. One key benefit
of AdamW is the decoupling of weight decay from gradient updates, which
encourages the model to maintain smaller weights, promoting stability
throughout the training process and generalizing to unseen data better.

####  2.3.1 Training schedule

Given the significantly increased parameter count resulting from the
strategies detailed in the previous section, combined with the diversity of
the dataset, and the inherent instability of 16-bit training, network
initialization plays a critical role. However, when using the weight
initialization strategy proposed in [33], we faced challenges due to the wide
range of orders of magnitude in the individual pixel-wise velocity samples,
which lead to overflow or underflow issues shortly after a few training
iterations. To address this issue, we devised a solution that involved
"warming up" the model weights using the simplest geometry, spherepacks. To
ensure the model’s weight expressiveness in predicting a wide range of
velocity outputs, we trained the model using five 2563 spherepacks. These
geometries covered porosity ranges of 8 to 25% and pressures ranging from 1 to
20 Mpa (refer to Figure 3). Although these spherepacks had similar
permeability values, the velocities spanned a significantly broader range.
This diverse training data enabled the model’s weights to capture the
necessary complexity and variability required for training with more diverse
samples. This training process, conducted over 72 hours with a learning rate
of 1e-6 and a batch size of one, necessitated the incorporation of a
constraint. If the loss exceeded 1e-4, the corresponding sample was ignored,
and the network parameters were not updated, since very high losses triggered
network updates results in the training diverging. This constraint, enforced
through the PyTorch distributed data parallel engine, effectively identified
and froze the unused parameters. While this constraint was frequently
triggered during the early stages of training, it gradually diminished over
time, and disappearing completely after a few hours, allowing the training
process to effectively adjust the weights.

After achieving stabilization of the loss function, we implemented a multi-
stage training approach. Initially, we added one extra sample from each
spherepack, selected randomly at a specific pressure. During this stage, we
utilized an increased learning rate of 1e-5 and a batch size of two, which
continued for 24 hours. Subsequently, we included all the samples sized 2563
in the training process, further increasing the learning rate to 1e-4 and
expanded the batch size to 8 for an additional 24 hours. Finally, in the last
training stage, we included the remaining samples (the 4803 volumes) and
performed eight random augmentations (including a 180-degree flip in the flow
axis and 90, 180, 270-degree flips on the XY-axes). This resulted in
approximately 1000 samples of varying sizes to train on. We carried out this
process for approximately 168 hours. In the final stage, we fine-tune our
model by incorporating stochastic weight averaging (SWA) to further enhance
its performance. SWA offers notable advantages as it is easy to implement,
improves generalization, and imposes minimal computational overhead. This
aligns with the findings of [45], who also reported enhanced generalization
capabilities through the use of SWA. Throughout this last stage, we utilized a
10% validation split and saved the model with the best loss value observed in
this set. To prevent overfitting, early stopping was employed after 300 epochs
without improvement.

##  3 Results

Our primary focus is the apparent permeability of a 3D structure under varying
confinement pressures. While it is possible to directly train a model to
predict this quantity as a floating-point number, our preliminary experiments
revealed a tendency for such models to overfit quickly and exhibit poor
generalization to unseen data. Therefore, we trained our model using the
complete flow field, even though our main interest lies only in the apparent
permeability. In essence, the prediction of the flow field serves as an
auxiliary task to support the training process. We would like to highlight the
significant distinction in computational requirements between the creation of
the training data and the testing phase. Generating the training data was
accomplished with relative ease, requiring around 10,000 CPU hours due to the
simplicity of the geometries. In contrast, the testing phase demanded a
substantial investment of computational resources, totaling close to 5 million
hours of CPU time.

We evaluated the performance of our model on unseen samples by quantifying the
ratio between predicted and simulated apparent permeability (refer to Figure
4). The model’s performance was assessed across different lithologies and
pressures. Remarkably, the model demonstrated excellent performance on unseen
samples, despite being exclusively trained on synthetic data. We observed a
well-defined normal distribution for the process-based samples in the test
set, which closely resemble the training set. For the remaining lithologies,
we observed strong agreement between the predicted and simulated permeability
values. Notably, the single fracture lithology exhibited the lowest
performance, likely due to its distinct characteristics such as roughness,
tortuosity, and pinch points, which differ from the single-fracture
counterparts in the training set. Confinement effects are magnified as pore
size and pressure decrease, we noted that our model tends to exhibit stronger
agreement at lower confinement pressures, suggesting that the input feature of
mean free path captures more salient features under such conditions. In terms
of accuracy,most of the test samples demonstrate results within an order of
magnitude. This achievement is particularly remarkable considering that the
permeability of the test data spans eight orders of magnitude and local
phenomena encompasses regimes from free molecular to viscous dominated.

Figure 4: a) Performance of the model in the training and testing sets for
each lithology. Each black bar represents one sample. b) Performance of the
model for different confinement pressures c) Performance of the model in all
the samples of the test set. The coefficients of determination for the train
and test sets are shown in the lower right corner.

##  4 Conclusions

For decades, approximations such as the cubic law in fractures, variations of
the Kozeny-Carman equation for granular media, and empirical equations based
on laboratory measurements have been relied upon to estimate permeability.
However, these approaches are known to have significant limitations and
inherent flaws. In contrast, machine learning provides a promising avenue for
permeability estimation, offering enhanced accuracy and improved performance.
Our developed model represents a significant advancement in this field,
setting a new state-of-the-art benchmark for permeability estimation.

We have developed a data-driven model that utilizes a large and diverse
dataset to accurately predict the apparent permeability of complex pore
networks. This model can be easily fine-tuned with new data, enabling improved
accuracy for specific types of porous media, or can be directly utilized as
is. To improve performance and avoid overfitting, we carefully optimized the
model size and training methods. Our approach incorporates a comprehensive
dataset, encompassing a wide range of sample sizes, features, and boundary
conditions.

By employing a convolutional neural network, we trained our model on synthetic
data and tested it on realistic heterogeneous samples, spanning various
thermodynamic conditions, pore sizes, pore structures, and computational
sample sizes. The versatility of this model lends itself to numerous practical
and scientific applications. For instance, it can aid in establishing new
relationships for fluid flow in porous media, and quantifying the effects of
diagenetic processes over time. Since this model is extremely fast, The real-
time prediction aspect could be exploited to optimize image segmentation
parameters for binarization during x-ray scanning, estimate relative
permeability curves while imagining an experiment, optimize porous media
synthesis for specific properties, and permeability estimation from cuttings
to complement logging while drilling, among others. In a broader view, our
approach could be adopted for other phenomena and domains which will benefit
from real-time predictions of commonly challenging-to-estimate properties.

##  5 Open Research

All the data used to train and test our model is published as [34] and can be
found at https://www.digitalrocksportal.org/projects/372.

## References

  * [1] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, 5 2015.
  * [2] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  ImageNet: A large-scale hierarchical image database.  In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255. Ieee, IEEE, 6 2009.
  * [3] Alex Krizhevsky.  Learning Multiple Layers of Features from Tiny Images.  2009\.
  * [4] Dorthe Wildenschild and Adrian P. Sheppard.  X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems.  Advances in Water Resources, 51:217–246, 1 2013.
  * [5] FT Mackenzie.  Sediments, Diagenesis, and Sedimentary Rocks.  2005\.
  * [6] Agnese Marcato, Gianluca Boccardo, and Daniele Marchisio.  From Computational Fluid Dynamics to Structure Interpretation via Neural Networks: An Application to Flow and Transport in Porous Media.  Industrial & Engineering Chemistry Research, 3 2022.
  * [7] Ali Kashefi and Tapan Mukerji.  Point-cloud deep learning of porous media for permeability prediction.  Physics of Fluids, 33(9), 9 2021.
  * [8] Ying Da Wang, Traiwit Chung, Ryan T Armstrong, and Peyman Mostaghimi.  ML-LBM: Predicting and Accelerating Steady State Flow Simulation in Porous Media with Convolutional Neural Networks.  Transport in Porous Media, (0123456789), 4 2021.
  * [9] Allen K. Ting, Javier E. Santos, and Eric Guiltinan.  Using Machine Learning to Predict Multiphase Flow through Complex Fractures.  Energies 2022, Vol. 15, Page 8871, 15(23):8871, 11 2022.
  * [10] Taylr Cawte and Aimy Bazylak.  A 3D convolutional neural network accurately predicts the permeability of gas diffusion layer materials directly from image data.  Current Opinion in Electrochemistry, 35:101101, 10 2022.
  * [11] Baris Caglar, Guillaume Broggi, Muhammad A. Ali, Laurent Orgéas, and Véronique Michaud.  Deep learning accelerated prediction of the permeability of fibrous microstructures.  Composites Part A: Applied Science and Manufacturing, 158:106973, 7 2022.
  * [12] Hou Lin Zhang, Hao Yu, Xin Heng Yuan, Heng Yu Xu, Marembo Micheal, Jia Ning Zhang, Hong Lin Shu, Gao Cheng Wang, and Heng An Wu.  Permeability prediction of low-resolution porous media images using autoencoder-based convolutional neural network.  Journal of Petroleum Science and Engineering, 208:109589, 1 2022.
  * [13] Pengfei Tang, Dongxiao Zhang, and Heng Li.  Predicting permeability from 3D rock images based on CNN with physical information.  Journal of Hydrology, 606:127473, 3 2022.
  * [14] Bernard Chang, Javier Santos, and Masa Prodanovic.  ElRock-Net: Assessing the Utility of Machine Learning to Initialize 3D Electric Potential Simulations.  Society of Core Analysists, 2022.
  * [15] Ying Yin, Zhiguo Qu, Maša Prodanović, and Christopher J. Landry.  Identifying the dominant transport mechanism in single nanoscale pores and 3D nanoporous media.  Fundamental Research, 1 2022.
  * [16] F. Javadpour, D. Fisher, and M. Unsworth.  Nanoscale gas flow in shale gas sediments.  Journal of Canadian Petroleum Technology, 46(10):55–61, 2007.
  * [17] L. J. Klinkenberg.  The permeability of porous media to liquids and gases.  Drilling and Production Practice 1941, pages 200–213, 1941.
  * [18] F. Javadpour.  Nanopores and Apparent Permeability of Gas Flow in Mudrocks (Shales and Siltstone).  Journal of Canadian Petroleum Technology, 48(08):16–21, 8 2009.
  * [19] Wenhui Song, Jun Yao, Yang Li, Hai Sun, Lei Zhang, Yongfei Yang, Jianlin Zhao, and Hongguang Sui.  Apparent gas permeability in an organic-rich shale reservoir.  Fuel, 181:973–984, 10 2016.
  * [20] Christopher J. Landry, Maša Prodanović, and Peter Eichhubl.  Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability.  International Journal of Coal Geology, 159:120–134, 4 2016.
  * [21] Robert G. Loucks, Robert M. Reed, Stephen C. Ruppel, and Ursula Hammes.  Spectrum of pore types and networks in mudrocks and a descriptive classification for matrix-related mudrock pores.  AAPG Bulletin, 96(6):1071–1098, 6 2012.
  * [22] Randall Marrett, Orlando J. Ortega, and Celinda M. Kelsey.  Extent of power-law scaling for natural fractures in rock.  Geology, 27(9):799, 9 1999.
  * [23] Li Chen, Wenzhen Fang, Qinjun Kang, Jeffrey De’Haven Hyman, Hari S. Viswanathan, and Wen Quan Tao.  Generalized lattice Boltzmann model for flow through tight porous media with Klinkenberg’s effect.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 91(3):033004, 3 2015.
  * [24] SA Tison.  Experimental data and theoretical modeling of gas flows through metal capillary leaks.  Vacuum, 44(11-12):1171–1175, 11 1993.
  * [25] Luca Marino.  Experiments on rarefied gas flows through tubes.  Microfluidics and Nanofluidics, 6(1):109–119, 2009.
  * [26] Lawrence M. Anovitz and David R. Cole.  Characterization and Analysis of Porosity and Pore Structures.  Reviews in Mineralogy and Geochemistry, 80(1):61–164, 1 2015.
  * [27] Mark E. Curtis, Carl H. Sondergeld, Raymond J. Ambrose, and Chandra S. Rai.  Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging.  AAPG Bulletin, 96(4):665–677, 2012.
  * [28] C. H. Sondergeld, R. J. Ambrose, C. S. Rai, and J.. Moncrieff.  Micro-Structural Studies of Gas Shales.  In All Days. SPE, 2 2010.
  * [29] Ayaz Mehmani, Maša Prodanović, and Farzam Javadpour.  Multiscale, Multiphysics Network Modeling of Shale Matrix Gas Flows.  Transport in Porous Media, 99(2):377–390, 9 2013.
  * [30] Rui Xu, Masa Prodanovic, and Christopher Landry.  Pore scale study of water adsorption and subsequent methane transport in clay in the presence of wettability heterogeneity, 2020.
  * [31] Jinlong Wu, Xiaolong Yin, and Heng Xiao.  Seeing permeability from images: fast prediction with convolutional neural networks.  Science Bulletin, 63(18):1215–1222, 9 2018.
  * [32] Ying Da Wang, Martin J. Blunt, Ryan T. Armstrong, and Peyman Mostaghimi.  Deep learning in pore scale imaging and modeling.  Earth-Science Reviews, 215(September 2020):103555, 2021.
  * [33] Javier E. Santos, Ying Yin, Honggeun Jo, Wen Pan, Qinjun Kang, Hari S. Viswanathan, Maša Prodanović, Michael J. Pyrcz, and Nicholas Lubbers.  Computationally Efficient Multiscale Neural Networks Applied to Fluid Flow in Complex 3D Porous Media.  Transport in Porous Media, 140(1):241–272, 10 2021.
  * [34] Javier E. Santos, Bernard Chang, Alex Gigliotti, Ying Yin, Wenhui Song, Maša Prodanović, Qinjun Kang, Nicholas Lubbers, and Hari Viswanathan.  A Dataset of 3D Structural and Simulated Transport Properties of Complex Porous Media.  Scientific Data, 9(1):579, 10 2022.
  * [35] Masa Prodanovic, Maria Esteva, Matthew Hanlon, Gaurav Nanda, and Prateek Agarwal.  Digital Rocks Portal: a repository for porous media images, 2015.
  * [36] Javier E. Santos, Michael J. Pyrcz, and Maša Prodanović.  3D Dataset of binary images: A collection of synthetically created digital rock images of complex media.  Data in Brief, 40:107797, 2 2022.
  * [37] Tom Bultreys and Wesley De Boever.  Belgian Fieldstone, 2020.
  * [38] Tom Bultreys, Luc Van Hoorebeke, and Veerle Cnudde.  Estaillades Carbonate #2, 2018.
  * [39] Stefanie Van Offenwert, Tom Bultreys, and Veerle Cnudde.  Saturated solute transport micro-{CT} dataset in sintered glass and {Bentheimer} sandstone, 2019.
  * [40] Catherine Spurin, Samuel Krevor, Martin Blunt, and Tom Bultreys.  Decane and brine injected into {Estaillades} carbonate - steady-state experiments, 2021.
  * [41] Ryan Armstrong and Peyman Mostaghimi.  Moura Coal, 2017.
  * [42] Kamaljit Singh and Martin J Blunt.  High resolution {X}-ray micro-tomography datasets for in-situ effective contact angle analysis in carbonate rocks, 2018.
  * [43] Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell, Saining Xie, and Facebook AI Research.  A ConvNet for the 2020s.  1 2022.
  * [44] Dan Hendrycks and Kevin Gimpel.  Gaussian Error Linear Units (GELUs).  6 2016.
  * [45] Hwan-Jin Song, Soonyoung Roh, Juho Lee, Giung Nam, Eunggu Yun, Jongmin Yoon, and Park Sa Kim.  Benefits of Stochastic Weight Averaging in Developing Neural Network Radiation Scheme for Numerical Weather Prediction.  Journal of Advances in Modeling Earth Systems, 14(10):e2021MS002921, 10 2022.


# Real-Time Trajectory Optimization in Robot-Assisted Exercise and
Rehabilitation

Humberto De las Casas  Nicholas Chambers  Hanz Richter  Kenneth Sparks
Mechanical Engineering Department, Cleveland State University, Cleveland, OH
44115 USA.  Health and Human Performance Department, Cleveland State
University, Cleveland, OH 44115 USA.

###### Abstract

This work focuses on the optimization of the training trajectory orientation
using a robot as an advanced exercise machine (AEM) and muscle activations as
biofeedback. Muscle recruitment patterns depend on trajectory parameters of
the AEMs and correlate with the efficiency of exercise. Thus, improvements to
training efficiency may be achieved by optimizing these parameters. The
optimal regulation of these parameters is challenging because of the
complexity of the physiological dynamics from person to person as a result of
the unique physical features such as musculoskeletal distribution.
Furthermore, these effects can vary due to fatigue, body temperature, and
other physiological factors. In this paper, a model-free optimization method
using Extremum Seeking Control (ESC) as a real-time optimizer is proposed.
After selecting a muscle objective, this method seeks for the optimal
combination of parameters using the muscle activations as biofeedback. The
muscle objective can be selected by a therapist to emphasize or de-emphasize
certain muscle groups. The feasibility of this method has been proven for the
automatic regulation of an ellipsoidal curve orientation, suggesting the
existence of two local optimal orientations. This methodology can also be
applied to other parameter regulations using a different physiological effects
such as oxygen consumption and heart rate as biofeedback.

Keywords: Robotics; biomechanics; muscle activation; control systems; sport
physiology; rehabilitation, extremum seeking control; optimization.
Number of words in the document: 3352.

##  1 Introduction

The performance of human exercise depends on the training parameters of the
exercise machine. This performance is strictly important to increase
efficiency and reduce the probabilities of injuries. Only in 2019, 468,315
injuries were reported in the US [1]. It is estimated that more than 2 of 5
gym users have had at least one injury while working out [2]. Therefore, a
smart system for the automatic regulation of these parameters is worth
developing. In this work, we explore the feasibility of this approach.

Research in areas of human-machine interaction (HMI) such as human performance
and rehabilitation has received a lot of attention in recent years. For
instance, performance enhancements have been evidenced by using robots as
advanced exercise machines (AEMs). AEMs have proven to be highly valuable to
human performance and rehabilitation developments over the past years by
combining exercise physiology with technology [3]. AEMs have electric motors
and control systems able to produce controllable trajectories and impedances
even in microgravity. As a result, the applications of these devices can prove
to be beneficial for astronauts by allowing them to receive similar stimulus
as to that achieved under the Earth’s gravity forces [4, 5, 6, 7]. The ability
of the AEM to provide configurable trajectories and impedances makes it
possible to enhance training patterns with a higher versatility than
conventional training methods [8, 9, 10]. Previous researches using a powered
rowing machine [11] and 4 DOFs robot reported diversity on the
musculoskeletal, cardiovascular and cardiorespiratory systems through the
manual regulation of trajectory and impedance parameters [12, 13]. HMI
applications focusing on rehabilitation have also achieved performance
improvements using robotics to help people with disabilities to recover some
of the motor skills that they used to have [14, 15, 16].

The integration between information technology and robotics has been
successfully applied in several interactive environments to support the
improvement of human exercise and rehabilitation practices [9, 17]. The
objective of this technology is to enhance manual dexterity similarly to that
seen in virtual reality games [18]. One of the first robotic systems developed
for therapy was the MIT-Manus (MITM). This robot-aid therapy system has
assisted stroke patients by helping them to improve their motor skills for
over 30 years now. MITM works by fitting the arm into a brace attached to the
robot arm. Meanwhile, the subject is required to follow a target trajectory.
This activity effectively assists the person in creating new neural
connections which, over time, enable the re-learning of certain muscular
processes [19, 20]. Several extensions and new developments for MITM [8] have
been performed using manual regulation of the system parameters.

The research regarding robot-aid systems also goes into exercise training
capabilities using integrated smart capabilities. The integration of HMI
systems with smart controllers is known as Smart-HMI [21]. This novel
mechatronic discipline considers the human dynamics inside of the closed-loop
system. Thus an opportunity to unite human physiological understanding with
mechanical, electrical, and information technology is presented. The first
model-free training optimization has been reported using a single degree of
freedom (DOF) manipulator with variable speed as the AEM to optimize the
user’s power output [22, 23]. The optimization was performed using the
biomechanical configuration of the user which was real-time estimated with its
force-velocity relationship. The previous research was replicated using a
similar 1-DOF AEM and 2 Extremum Seeking Control (ESC) algorithms [24]. These
algorithms working independently were used to maximize the user’s power output
and estimate the user’s torque. Regarding other research into the area of
rehabilitation, biological model estimations have been successfully developed
and reported [25, 26]. The semi-automatic selection of machine parameters
based on biofeedback has made it possible to achieve better results when
compared to conventional rehabilitation practices.

Typical physiological measurements collected for human performance research
include electromyography (EMG), heart rate (HR), oxygen consumption (VO2), and
energy requirements. Despite the simplicity of their measurements,
physiological effects are complex as everyone has varying physiological
functioning. They depend on the unique features of each person (organismic
variables) such as force capacity, musculoskeletal distribution, body mass
index, and flexibility. Besides, they display time-varying dynamics due to
fatigue, body temperature, level of hydration, etc. As a result, modeling and
controlling muscle performance is challenging but necessary in some areas of
research [27]. Human dynamic models relating physical activity with muscle
activations have been reported [28]. However, these models do not take into
consideration the time-varying effects of training. They are only based on the
dynamics of the musculoskeletal distribution. Consequently, they are useful
and reliable for simulations, but not recommended for real-time experiments.

Figure 1: Exercise environment for the advanced training. $P$ and $P_{des}$
are the current and desired robot end-effector position (red and blue dot
respectively).

During the exercise protocol, the subject was required to follow the desired
position while receiving feedback from his or her current position. Thus, a
graphical user interface (GUI) is used to provide the required visual
feedback. The complete environment can be seen in Figure 1 where the blue line
represents the ellipsoidal trajectory of fixed axis lengths and variable
orientation that was designed based on the desired orientation using forward
dynamics. The user’s actual position $P$ (defined by the end-effector position
of the robot) is labeled with a red dot. The desired position $P_{des}$
(rotating periodically over the blue ellipsoidal trajectory) is labeled with a
blue dot. The user is required to track the blue dot as close as possible. EMG
was used to measure and provide biofeedback of the muscle activations to the
system. The model-free framework making use of the robot positions and muscle
activations regulates the orientation aiming to maximize muscle performance.
This muscle performance is defined as the moving average of the multiplication
of a muscle weight vector and the muscle activation derived as follows:

$\max_{\theta}\quad
J(t)=\frac{t_{s}}{t_{rev}}\sum^{t}_{(i=t-t_{rev}/t_{s})}\bigg{(}W_{m}M(t)\bigg{)},$
(1)

where $\theta$ is the ellipsoidal trajectory orientation, $t_{s}$ is the
sample time set at 5e-4 seconds, $t_{rev}$ is the frequency of revolution of
the reference cursor (blue dot), $t$ is the current time, $W_{m}$ is the
weight muscle vector, and $M$ is the vector of muscle activations. The weight
muscle vector receives the priority values for each muscle. This priority can
be selected by professionals or therapists to emphasize or de-emphasize
certain muscle groups. For this study, 3 possible priorities were selected to
be assigned arbitrarily. These priorities were low, medium, and high related
to priority gains of 1, 3, and 5 respectively. Therefore, muscles with higher
priority (higher gain) would receive more maximization focus.

##  2 Dynamics

The AEM used in this study was a WAM robot. The WAM is 4 DOFs articulated
linkage robot from Barrett Advanced Robotics. This lightweight cable-driven
manipulator is exceptionally dexterous, low-friction naturally backdrivable.
Therefore, the robot is able to provide an almost negligible interaction
resistance between the subject and the machine. The dynamics of this robot
[29] given in joint coordinates are derived as:

$D(q)\ddot{q}+C(q,\dot{q})\dot{q}+g(q)=\tau+J^{T}F_{ext},$ (2)

where $q=[q_{1},q_{2},q_{3},q_{4}]$ is a vector of joint displacements, $D(q)$
is the inertia matrix, $C(q,\dot{q})$ is the centripetal and Coriolis effects,
$g(q)$ is the gravity vector, $\tau$ is the control torque, $J$ is the
Jacobian, and $F_{ext}$ is the external force representing the interaction
force between subject and robot.

Muscle dynamics are mainly affected by time-invariant independent, organismic,
and intervening variables such as robotic parameters, musculoskeletal
distribution, performance status, level of hydration, and mood; and time-
variant variables like muscle temperature and fatigue. As a result, muscle
dynamics can be described as an unknown dynamic equation as follows:

$\dot{M}=f(M,u(t)),$ (3)

where $M=[M_{1},M_{2},M_{3},M_{4}]$ is a vector of muscle activations and
$u(t)$ is a nonlinear time-varying function including the ellipsoidal
orientation ($\theta$) and all the unknown and unmeasurable variables.

##  3 Control Architecture

Two controllers working simultaneously were involved in this project. The
first was a PD-gravity compensator controller working isolated on the AEM. The
second controller works on the complete HMI system relating the robot dynamics
with muscle performance.

###  3.1 Robot - PD-Gravity Compensation Control

In order to work on the planar ellipsoidal trajectory, 2 of the 4 joints of
the robot were fixed by applying a PD controller acting on the third and
fourth joints. Gravity compensation was applied to increase robot
maneuverability. The control law was defined as:

$\tau=P_{gain}e+D_{gain}\dot{e}+g(q),$ (4)

where $P_{gain}$ and $D_{gain}$ are diagonal PD gain matrices (zero diagonal
values for the first 2 DOFs), and $e$ is the error vector between the desired
and the current joint positions. The end-effector robot position ($P$) is sent
and showed in the GUI during the complete experiment.

###  3.2 Main Controller - Extremum Seeking Control (ESC)

The main controller is based on the use of the perturbation-based Extremum
Seeking Control (ESC) [30, 31]. The ESC scheme can be seen in Figure 2. ESC is
a branch of adaptive control developed for optimization. The aim of this
control is to enforce the output of a dynamical system to converge to an
unknown maximum or minimum operating points. Several ESC algorithms have been
developed and studied from the first appearance in the 1920’s [32, 31]. Some
of them are model-based methods. However, the popularity of ESC lies in its
model-free optimization algorithms and its ability to operate in real-time.

Figure 2: Scheme of ESC for SISO with deterministic perturbation, where
$\omega_{l}$ and $\omega_{h}$ are the low-pass and high-pass cutoff
frequencies respectively, $k$ is the gain, $a$ is the amplitude of the
perturbation, and $\omega_{1}$ is the frequency of the perturbation.

Comparable to any other optimization method, the basis of the ESC is provided
by the general optimal conditions [33]. For a general nonlinear system:

$\dot{x}=f(t,x,u),$ (5)

and the following performance output:

$y=J(t,\theta),$ (6)

where $t$ is the time, $x^{T}=[x_{1},x_{2},...,x_{n}]$ is the state vector,
$u$ is the input, and the functions $f:D\times\mathbb{R}\to\mathbb{R}^{n}$ and
$J:D\to\mathbb{R}$ are sufficiently smooth on $D\subseteq\mathbb{R}^{n}$.

Assuming stability on the system (or being stabilizable) at each of the
equilibrium points by a local feedback controller. Suppose the existence of
the following smooth control law:

$u=\alpha(x,\theta),$ (7)

parameterized by a scalar parameter $\theta$. The closed-loop system becomes:

$\dot{x}=f\bigg{(}t,x,\alpha(x,\theta)\bigg{)},$ (8)

therefore, the equilibrium points of the system will be parameterized by
$\theta$.

The existence of a local minimizer (or maximizer) is guaranteed under specific
assumptions about the closed-loop system. Assumptions and stability analysis
are derived in [32, 34].

ESC was selected to deal with the nonlinear time-varying muscle dynamics
including the unknown and unmeasurable parameters. This controller receives
processed muscle performance. Then, based on its variation, the controller
regulates the orientation of the ellipsoidal curve ($\theta$) to be
transformed using forward dynamics into the desired trajectory ($P_{des}$).

The ESC parameters used on this project (related to Figure 2) were 0.1, 1
$rad/s$, 0.1 $rad/s$ ,0.5 $rad/s$, and 1000 for the perturbation amplitude
($a$), perturbation frequency ($\omega_{1}$), low-pass filter cutoff frequency
($\omega_{l}$), high-pass filter cutoff frequency ($\omega_{h}$), and gain
($k$) respectively.

###  3.3 Closed-loop system scheme for the model-free optimization

Figure 3: Closed-loop system scheme for the model-free optimization approach.

The closed-loop system scheme can be seen in Figure 3. The process begins with
the GUI providing visual feedback to the user (see Figure 3-a). Next, the user
drives the end-effector of the robot following the desired ellipsoidal
trajectory (see Figures 6 and 3-b). As a result of the training, muscle
activations are produced and measured with EMGs as raw signals (see Figure
3-c). These signals are recorded at a frequency of 2 kHz and real-time
processed to obtain muscular activation (see Figure 3-d). The complete process
includes a bandpass second-order Butterworth filter between 30 and 950 Hz, a
full-wave rectification, and a low pass second-order Butterworth filter at 50
Hz. The muscular activations are then used to compute the overall muscle
performance which is then maximized (see Figure 3-e). Due to the high
sensitivity of the muscle signal, a second-order filter is applied at the
output of the performance function (see Figure 3-f). This filter introduces a
delay to the signal suppressing possible perturbations. The main controller
running the model-free optimization (see Figure 3-g) makes use of the
processed muscle performance to regulate the orientation of the visual
geometric path to be tracked by the subject (see Figure 3-h). While the
orientation is being adjusted to achieve maximal performance, the cycle
continues to repeat. The optimal solution is considered found after the
optimization variable remains between $\pm 10^{\circ}$ for at least 10 seconds
(convergence criterion).

##  4 Experimental Protocol

The complete HMI system used for the experiments can be seen in Figure 5. This
system includes a WAM robot as the AEM, a desktop PC with Xenomai Kernel, a
laptop PC with MATLAB and ControlDesk, a MicroLabBox real-time data
acquisition and control system (dSPACE GmbH, Paderborn, Germany), a set of
wireless EMG sensors (Trigno Wireless EMG, Delsys Inc.), and 2 DAQs (Labjack
U3-HV and Labjack T7). The PC with Xenomai Kernel was used to control the WAM
robot by CANbus communication and send the robot states (related to the human
arm position) to the laptop PC. The laptop PC was used to manage the setting
parameters of the MicroLabBox. The MicroLabBox was utilized to read the
variable states, perform the real-time processing, record the data, and run
the main control. The Delsys Trigno Wireless EMG system aided in collecting
the EMG activity of the four muscles selected for the study. The Labjack U3-HV
was used to transmit the robot states from one PC to the other. The Labjack T7
was used to transmit the desired trajectory from the main controller to a
visualization screen.

Based on the experimental protocol, a subject was required to follow an
ellipsoidal trajectory with variable orientation auto-regulated by using his
muscular activations as biofeedback. Based on the type of movement and the
range of motion, 4 muscles surrounding the glenohumeral joint [35] (see Figure
4) were selected in the following order:

  1. 1.

Lateral deltoid: for shoulder abduction and stabilization during flexion.

  2. 2.

Anterior deltoid: for shoulder flexion and inward rotation.

  3. 3.

Biceps brachii: for lower arm flexion, and assistance and stabilization during
shoulder flexion.

  4. 4.

Pectoralis major: for shoulder adduction in the frontal plane.

One healthy right-handed male subject of the age of 23 years, height of 180.4
cm, and weight of 93.2 kg performed 3 sets of experiments on 3 different days
(see Figure. 6). The experimental procedure began with warm-up and isometric
tests [36]. Isometric tests were used to assess muscle strength for the EMG
sensor calibration. To establish these contractions, the subject was required
to maintain different constant positions which muscles were capable of
producing maximal muscular contracting forces (highest activations) [37].
Then, the experiment proceeded with the trials. During each day, one set of 2
trials (one trial next to the other one) was performed by using the same
parameter configuration. The first set of experiments used the muscle weight
vector $W_{m}=[1,5,3,5]$. This set gave the highest priority to the anterior
deltoid and pectoralis major, medium priority to the biceps brachii, and the
lowest priority to the lateral deltoid. The second set used the muscle weight
vector $W_{m}=[3,5,1,1]$. This set gave the highest priority to the anterior
deltoid, medium priority to the lateral deltoid, and the lowest to biceps
brachii and pectoralis major. And the third set used the muscle weight vector
$W_{m}=[1,1,5,5]$). This latest set gave the highest priority to the biceps
brachii and pectoralis major and the lowest to the lateral and anterior
deltoid.

Figure 4: 4OptimX experiment configuration and EMG locations. Figure 5: Scheme
of the HMI system for training optimization. Figure 6: Experiment of human
training with AEM.

##  5 Results and Discussion

Figure 7: First set of experiments with weight factor $W_{m}$=[1, 5, 3, 5] Figure 8: Second set of experiments with weight factor $W_{m}$=[3, 5, 1, 1] Figure 9: Third set of experiments with weight factor $W_{m}$=[1, 1, 5, 5] Table 1: Experiments results - convergence of the solutions based on the selected weight factors (1, 3, and 5 representing the low, medium, and high priority gains respectively). Weight factor | Trial | Solution | Convergence
---|---|---|---
($W_{m}$) | (N∘) | (∘) | Time (s)
$[1,5,3,5]$ | 1 | 39.50 | 24.6
$[1,5,3,5]$ | 2 | 27.60 | 25.6
$[3,5,1,1]$ | 1 | 30.77 | 23.9
$[3,5,1,1]$ | 2 | -44.09 | 27.3
$[1,1,5,5]$ | 1 | -58.14 | 27.7
$[1,1,5,5]$ | 2 | 40.40 | 28.2

Upon completion of the 6 experimental trials, differences in the muscle
performance between trials from the same experimental set were revealed. A
higher performance value was exhibited in the second trials during all the
experimental sets. This result seems to be associated with the reduced effort
capacity of the muscles to perform activity when they start to get fatigued.
However, despite the difference in the muscle performance between trials, they
converged to similar solutions depending on the initial conditions and the
muscle weight vector. These results suggest that the formulated model-free
optimization method successfully enabled the user to exercise optimally.

It is important to consider that the optimal trajectory orientations don’t
only depend on the subject’s musculoskeletal system (lengths, weights,
distributions, etc), but also on his/her interaction with the machine.
Although the robot offers little resistance, the user required to overcome the
robot weight, inertia, and friction representing an additional resistance for
the user and affecting the results. Besides, muscle dynamics are permanently
changing. Therefore, ESC remains varying slightly even after reaching the
optimal solution. However, based on the established convergence criterion, the
solutions were computed and found in the neighborhood from $30^{\circ}$ to
$40^{\circ}$ and its symmetric respect to the vertical axis from $-60^{\circ}$
to $-45^{\circ}$ (see Figure 10). These results suggest that an inclined
ellipsoidal trajectory seems to provide a higher muscle demanding able to
maximize the activations for training performance.

Figure 10: Convergence solutions from the 3 sets of experiments.

It is important to note that the success of the model-free approach lies in
the accuracy of the user to track the desired trajectory and the selection of
the configuration parameters. Regarding the accuracy of the tracking, it is
recommended that the user is introduced and provided time to practice with the
robot prior to testing so as to familiarize with the mechanisms. During the
present study, tracking was acceptable and consistent between trials despite
the presence of minor tracking errors. The mean of the RMS tracking errors was
reported at 2.87 cm what proves the controller was robust enough to always
converge to the neighborhoods of the optimizers. Regarding the configuration
parameters, the framework might require some pre-tests to find a good
combination of parameters for each subject. For instance, higher gains or
frequencies might produce faster convergence, but at the cost of a higher
sensitivity which is not recommended on this approach because of undesired
performances. On the other hand, low gains or frequencies might never achieve
a convergence or not being able to deal with the time-varying dynamics due to
the fatigue and the thermogenic effect of the muscles. Nonetheless, the
parameters chosen for the 3 experimental were accurate enough to produce
robustness and a convenient convergence speed to deal with these variations.

##  6 Future Work

Some challenges encountered while conducting these experiments include the
high sensitivity of the EMGs. This proved problematic as the feedback resulted
in involuntary elbow rotations, producing alterations in the muscle dynamics.
The possibility of including human motion capture to feedback into the
controller could be promising as a means of providing more information on
concurrent motion dynamics. Another issue was associated with the electrical
noise produced by the other systems being used in the testing environment.
Aside from the current digital processes in this model, analog processing will
be included as future development. In addition, as an alternate solution to
muscle performance, estimation methods will be tested. Due to the high
computing cost of the high-frequency sampling, estimation methods such as
Kalman filtering show potential in improving the efficiency of this method.

Coincidentally in the 3 experimental sets, longer convergence times were
observed in the second trials (after a few minutes of working out). This
result suggests that there is a possible relationship between convergence time
and fatigue. It is known that muscles consist of many motor units that are not
fully active at the beginning of the workout, but they start to activate
together with the increasing of fatigue [38, 39]. The increase in the
activation of motor units produces an increase in the muscle activations
similarly to an increasing in the framework gain producing a higher
sensitivity, and thus undesired performances. These sudden changes might not
only delay the convergence but also they might even block it. Therefore, it
can be concluded that independently of the accuracy in the initial setting
parameter selection, recalibration might become needed after a few minutes of
training. The increase in the sensitivity previously observed might be solved
by decreasing the framework gain, thus future studies could include automatic
parameter calibrations to overcome this current limitation.

As a final conclusion, results suggest the feasibility of the formulated
model-free optimization for trajectory orientation by targeting muscular
performance. The feasibility of the approach will be also tested for impedance
regulation. As a long term objective, research on multi-variable optimization
of trajectory and impedance parameters using more physiological effects
related to the cardiovascular and cardiorespiratory systems (such as heart
rate, oxygen consumption, and energy expenditure) will be considered. This
approach will make it possible to optimize the training of an individual to
provide both an aerobic and anaerobic stimulus.

## Acknowledgements

We wish to confirm that there are no known conflicts of interest associated
with this publication and there has been no significant financial support for
this work that could have influenced its outcome. This research was supported
by the NSF (grant 1544702). This research has an IRB provided by Cleveland
State University (reference number 30305-RIC-HS) which covers for Ethical
Approval. Subjects on this research have read and signed informed consents.

## References

  * [1] Insurance Information Institute, Facts + statistics: Sports injuries, https://www.iii.org/fact-statistic/facts-statistics-sports-injuries, accessed: 2021-1-7.
  * [2] Protectivity Insurance, 41% of gym users have been injured whilst working out, https://www.protectivity.com/knowledge-centre/41-of-gym-users-have-been-injured-whilst-working-out/, accessed: 2021-1-7.
  * [3] D.-W. Jung, D.-S. Park, B.-S. Lee, M. Kim, Development of a motor driven rowing machine with automatic functional electrical stimulation controller for individuals with paraplegia: a preliminary study, in: Annals of Rehabilitation Medicine, 2012, pp. 379–385.
  * [4] M. A. A. Hargens, S. Parazynski, J. Friden, Muscle changes with eccentric exercise: Implications on earth and in space, Tech. rep., NASA (1989).
  * [5] M. Vogt, H. H. Hoppeler, Eccentric exercise: mechanism and effects when used as training regime or training adjunct, in: The American Physiological Society, 2014, pp. 1446–1454.
  * [6] L. Ploutz-Snyder, J. Ryder, K. English, F. Haddad, K. Baldwin, Risk of impaired performance due to reduced muscle mass, strength, and endurance, Tech. rep., NASA (March, 9 2015).
  * [7] S. Smith, S. Zwart, M. Heer, Human Adaptation to Spaceflight: The Role of Nutrition, NASA, 2014.
  * [8] H. Krebs, M. Ferraro, S. Buerger, M. Newbery, A. Makiyama, M. Sandmann, D. Lynch, B. Volpe, N. Hogan, Rehabilitation robotics: Pilot trial of a spatial extension for mit-manus, Journal of neuroengineering and rehabilitation 1 (2004) 5.  doi:10.1186/1743-0003-1-5.
  * [9] H. Krebs, J. Palazzolo, L. Dipietro, M. Ferraro, J. Krol, K. Rannekleiv, B. Volpe, N. Hogan, Rehabilitation robotics: Performance-based progressive robot-assisted therapy, Autonomous Robots 15 (2003) 7–20.  doi:10.1023/A:1024494031121.
  * [10] H. De las Casas, H. Richter, A. van den Bogert, Design and hybrid impedance control of a powered rowing machine, in: ASME 2017 Dynamic Systems and Control Conference, 2017.
  * [11] H. De las Casas, Desing and control of a powered rowing machine with programmable impedance, Master’s thesis, Cleveland State University (2017).
  * [12] H. De las Casas, K. Kleiss, H. Richter, K. Sparks, A. van den Bogert, Eccentric training with a powered rowing machine, Medicine in Novel Technology and Devices 2 (2019) 100008.
  * [13] H. De las Casas, S. Bianco, H. Richter, Targeted muscle effort distribution with exercise robots: Trajectory and resistance, Medical Engineering & Physics.
  * [14] C. Sicuri, G. Porcellini, G. Merolla, Robotics in shoulder rehabilitation, Muscles, ligaments and tendons journal 4 (2014) 207–13.
  * [15] W. H. Chang, Y.-H. Kim, Robot-assisted therapy in stroke rehabilitation, Journal of Stroke 15 (2013) 174–181.
  * [16] H. Schmidt, C. Werner, R. Bernhardt, S. Hesse, J. Krüger, Gait rehabilitation machines based on programmable footplates, Journal of NeuroEngineering and Rehabilitation.
  * [17] H. Krebs, N. Hogan, B. Volpe, M. Aisen, L. Edelstein, C. Diels, Overview of clinical trials with mit-manus: A robot-aided neuro- rehabilitation facility, Technology and health care : official journal of the European Society for Engineering and Medicine 7 (1999) 419–23.  doi:10.3233/THC-1999-7606.
  * [18] J. Hidler, D. Nichols, M. Pelliccio, K. Brady, Advances in the understanding and treatment of stroke impairment using robotic devices, Topics in stroke rehabilitation 12 (2005) 22–35.  doi:10.1310/RYT5-62N4-CTVX-8JTE.
  * [19] N. Hogan, H. I. Krebs, J. Charnnarong, P. Srikrishna, A. Sharon, MIT-MANUS: a workstation for manual therapy and training II, in: H. Das (Ed.), Telemanipulator Technology, Vol. 1833, International Society for Optics and Photonics, SPIE, 1993, pp. 28 – 34\.  doi:10.1117/12.142124.
URL https://doi.org/10.1117/12.142124

  * [20] N. Hogan, H. I. Krebs, J. Charnnarong, P. Srikrishna, A. Sharon, Mit-manus: a workstation for manual therapy and training. i, [1992] Proceedings IEEE International Workshop on Robot and Human Communication (1992) 161–165.
  * [21] F. Harashima, Human adaptive mechatronics, in: IEEE Workshop on Advanced Robotics and its Social Impacts, 2005., 2005, pp. 80–.  doi:10.1109/ARSO.2005.1511624.
  * [22] P. Y. Li, R. Horowitz, Control of smart exercise machines. i. problem formulation and nonadaptive control, IEEE/ASME Transactions on Mechatronics 2 (4) (1997) 237–247.  doi:10.1109/3516.653048.
  * [23] P. Y. Li, R. Horowitz, Control of smart exercise machines part ii: Self-optimizing control, IFAC Proceedings Volumes 29 (1) (1996) 5727 – 5732, 13th World Congress of IFAC, 1996, San Francisco USA, 30 June - 5 July.
  * [24] X. T. Zhang, D. M. Dawson, W. E. Dixon, B. Xian, Extremum-seeking nonlinear controllers for a human exercise machine, IEEE/ASME Transactions on Mechatronics 11 (2) (2006) 233–240.  doi:10.1109/TMECH.2006.871896.
  * [25] O. Barzilay, A. Wolf, Adaptive rehabilitation games, Journal of Electromyography and Kinesiology 23 (1) (2013) 182 – 189.
  * [26] Q. Wang, W. Chen, A. A. A. Timmermans, C. Karachristos, J. B. Martens, P. Markopoulos, Smart rehabilitation garment for posture monitoring, in: 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 2015, pp. 5736–5739.
  * [27] M. A. A. Rahman, M. Mizukawa, Model-based development and simulation for robotic systems with sysml, simulink and simscape profiles, International Journal of Advanced Robotic Systems 10 (2) (2013) 112.  doi:10.5772/55533.
  * [28] H. Warner, Simulation and control at the boundaries between humans and assistive robots, Ph.D. thesis, Cleveland State University (December 2019).
  * [29] S. H. M.W. Spong, M. Vidyasagar, Robot Modeling and Control, Wiley, 2005.
  * [30] M. Krstić, H.-H. Wang, Stability of extremum seeking feedback for general nonlinear dynamic systems, Automatica 36 (4) (2000) 595 – 601.
  * [31] M. Krstic, Extremum Seeking Control, Springer London, 2014.
  * [32] K. B. Ariyur, M. Krstic, Real Time Optimization by Extremum Seeking Control, John Wiley & Sons, Inc., New York, NY, USA, 2003.
  * [33] J. Nocedal, S. J. Wright, Numerical Optimization, Springer, New York, NY, 1999.
  * [34] C. Zhang, R. Ordonez, Extremum-Seeking Control and Applications A Numerical Optimization-Based Approach, Springer London Dordrecht Heidelberg New York, 2012\.  doi:10.1007/978-1-4471-2224-1.
  * [35] Y. Blache, M. Begon, B. Michaud, L. Desmoulins, P. Allard, F. Dal Maso, Muscle function in glenohumeral joint stability during lifting task, PLOS ONE 12 (12) (2017) 1–15.  doi:10.1371/journal.pone.0189406.
  * [36] N. Babault, M. Pousson, A. Michaut, J. Van Hoecke, Effect of quadriceps femoris muscle length on neural activation during isometric and concentric contractions, Journal of Applied Physiology 94 (3) (2003) 983–990.
  * [37] S. Garcia, J. Dueweke, C. Mendias, Optimal joint positions for manual isometric muscle testing, Journal of Sport Rehabilitationdoi:10.1123/jsr.2015-0118.
  * [38] J. Z. Liu, R. W. Brown, G. H. Yue, A dynamical model of muscle activation, fatigue, and recovery, Biophysical Journal 82 (5) (2002) 2344 – 2359.
URL http://www.sciencedirect.com/science/article/pii/S000634950275580X

  * [39] R. Enoka, Morphological features and activation patterns of motor units, Journal of clinical neurophysiology : official publication of the American Electroencephalographic Society 12 (6) (1995) 538—559.  doi:10.1097/00004691-199511000-00002.