date$-$ date$-$ 11institutetext: Department of Ocean Engineering and Naval
Architecture
Indian Institute of Technology Kharagpur, India

# Evaluation of machine learning algorithms for predictive Reynolds stress
transport modeling

J P PANDA  H V WARRIOR<EMAIL_ADDRESS>
(_; _)

###### Abstract

The application machine learning (ML) algorithms to turbulence modeling has
shown promise over the last few years, but their application has been
restricted to eddy viscosity based closure approaches. In this article we
discuss rationale for the application of machine learning with high-fidelity
turbulence data to develop models at the level of Reynolds stress transport
modeling. Based on these rationale we compare different machine learning
algorithms to determine their efficacy and robustness at modeling the
different transport processes in the Reynolds Stress Transport Equations.
Those data driven algorithms include Random forests, gradient boosted trees
and neural networks. The direct numerical simulation (DNS) data for flow in
channels is used both as training and testing of the ML models. The optimal
hyper-parameters of the ML algorithms are determined using Bayesian
optimization. The efficacy of above mentioned algorithms is assessed in the
modeling and prediction of the terms in the Reynolds stress transport
equations. It was observed that all the three algorithms predict the
turbulence parameters with acceptable level of accuracy. These ML models are
then applied for prediction of the pressure strain correlation of flow cases
that are different from the flows used for training, to assess their
robustness and generalizability. This explores the assertion that ML based
data driven turbulence models can overcome the modeling limitations associated
with the traditional turbulence models and ML models trained with large
amounts data with different classes of flows can predict flow field with
reasonable accuracy for unknown flows with similar flow physics. In addition
to this verification we carry out validation for the final ML models by
assessing the importance of different input features for prediction.

###### keywords:

Machine learning, Pressure strain correlation,Turbulence modeling, CFD

$R_{ij}$ | Reynolds stress tensor
---|---
$b_{ij}$ | Reynolds stress anisotropy
$\delta_{ij}$ | Kronecker delta
$k$ | turbulent kinetic energy
$\rho$ | density
$\nu$ | Kinetic viscosity
$U_{i}$ | Mean velocity component
$u_{i}$ | Fluctuating velocity component
$P_{ij}$ | Production of turbulence
$T_{ijk}$ | Diffusive transport
$\epsilon_{ij}$ | Dissipation rate tensor
$\phi_{ij}$ | Pressure strain correlation

##  1 Introduction

Computational fluid dynamics(CFD) utilizes numerical techniques to model and
solve problems that involve flow of fluids. CFD based methods[Versteeg and
Malalasekera (2007)] have became popular with the advancement in computational
facilities and formulation of models for turbulence[Pope (2000)]. The CFD
based techniques are widely used in the early years of 21st century. The
success of CFD techniques is largely dependent on the development of the
techniques to model the turbulence. Mainly there are four broad classes of CFD
approaches in the turbulence modeling framework, those are eddy viscosity
models[Pope (1975), Pope (2000), Panda et al. (2020)], Reynolds stress
transport models[Mishra and Girimaji (2017), Panda and Warrior (2018)], Large
eddy simulation[Pope (2000)] and Direct numerical simulation[Lee and Moser
(2015), Lee and Moser (2018)]. In industrial applications eddy viscosity based
are widely used, since those have least level of complexity and a few number
of equations are solved for the modeling of the flow field. In recent years
with booming in computational facilities researchers have started using LES
for flow prediction in complex and larger domains. But LES can not be used for
modeling and prediction complex flows of Engineering interest, since the cost
of simulations of LES is very high. The LES and DNS data are mainly used for
development and calibration of new turbulence models. In contrast to all the
above mentioned approaches, the Reynolds stress transport models[Panda and
Warrior (2018), Mishra and Girimaji (2017)] are superior to eddy viscosity
models, since in these models equations are solved for all the components of
Reynolds stresses from which the Reynolds stress field is approximated. The
building block of such models are the Reynolds stress transport
equations[Panda (2019)], in which, the two important terms that need to be
modeled are the pressure strain correlation[Panda (2020)] and the dissipation
term. The models for all these terms are developed by calibrating with few
cases of experiments such as grid turbulence, grid turbulence with mean strain
and shear flows. So, when these models are applied to other realistic cases of
flow predictions may provide inaccurate results.

To overcome these difficulties recent emphasis of turbulence modelers have
been shifted towards development of data-driven turbulence models[Duraisamy et
al. (2019)]. Machine learning is a process of function approximation, the
function can be approximated by correlating the input parameters with the
output. There are various machine learning algorithms. Few important of those
are neural networks, Random forests[Wang et al. (2017), Kaandorp and Dwight
(2020), Luan and Dwight (2020)], gradient boosted trees[Wu et al. (2020)],
Gene expression programming[Zhao et al. (2020)] and sparse symbolic
regression[Schmelzer et al. (2020)]. In turbulence modeling, the problem is
posed as one of supervised learning, where the model attempts to minimize the
prediction error over a training data-set. This data required for training,
validation and testing the ML models can be obtained from experimental
investigations, high fidelity DNS and LES data-sets of turbulent flows. In the
context of supervised learning, data driven approaches of turbulence modelling
can be classified into three main categories a) quantification of
uncertainties in the RANS models (with this approach the uncertainties in
Reynolds stress tensor can be quantified), b) finding the discrepancies in the
model coefficients and magnitude of the terms in the governing equation and c)
direct modeling of turbulence parameters such as Reynolds stress(EVM),
pressure strain correlation (RSTM) and the subgrid scale stress(LES) in terms
of the mean flow parameters. Among the three mentioned approaches of data-
driven turbulence modeling, the models developed through the third approach
can be trained with huge database of turbulence and can be applied for
prediction of unseen flow cases of similar flow physics.

[Singh et al. (2017)] developed model augmentations for Spalart-Allmaras(SA)
turbulence model using adjoint based full field inference using experimentally
measured lift coefficient data. These models forms are reconstructed using
neural networks, and applied in CFD solver, to predict flow in different
operating conditions. [Tracey et al. (2015)] used a shallow neural network
(one hidden layer) to model the source terms of the SA turbulence model.
[Parish and Duraisamy (2016)] learned a turbulence production term using
machine learning and applied that to the $k$-equation of the $k-\omega$
turbulence model. [Maulik et al. (2018)] used neural networks to model eddy
viscosities for RANS simulations. [Ling et al. (2016)] used a special kind of
neural network called tensor basis neural network to model the Reynolds stress
anisotropy and also compared the model prediction against a simple multi layer
perception. The optimal hyper-parameters of the model were recommended using
Bayesian optimization. [Zhu et al. (2019)] used neural networks to construct a
mapping function between the turbulent eddy viscosity and the mean flow
variables and ML model completely replaces the partial differential equation
model. [Wu et al. (2018)] proposed a physics based implicit treatment to model
Reynolds stresses using random forests. The optimal eddy viscosity and non-
linear part of the Reynolds stresses were both predicted. They used seven
distinct physics based features to model the Reynolds stresses, those are
strain rate tensor, rotation rate tensor, pressure gradient, turbulence
kinetic energy gradient, wall distance based Reynolds number, turbulence
intensity and ratio of turbulent time-scale to mean strain time-scale. They
used different approaches to achieve in-variance in machine learning.
[Weatheritt and Sandberg (2016)] used Gene expression programming to formulate
non-linear consecutive stress-strain relationship. The mathematical model was
created by using high fidelity and uncertainty measures. The learning method
has the capability to produce a constraint free model. [Weatheritt and
Sandberg (2017)] used symbolic regression to model the algebraic form of the
Reynolds stress anisotropy tensor. The equations were trained using hybrid
RANS/LES data and the new model was employed in RANS closure to test the
prediction of model in 3D geometries. [Schmelzer et al. (2020)] discovered
algebraic Reynolds stress models using sparse regression and they have used
high fidelity LES/DNS data for training and cross validation of the model.
There case of separated flows were considered, those are periodic hills,
converging-diverging channel and curved backward facing step. The prediction
of the machine learnt model was better than the $k-omega$ SST model. [Huijing
et al. (2021)] used the model developed by [Schmelzer et al. (2020)] to
predict the fully three dimensional high Reynolds number flows, e.g. wall
mounted cubes and cuboids. [Fang et al. (2020)] used neural networks to model
the Reynolds stress anisotropy using neural networks. They proposed different
modification of the neural network structure to accommodate effect of Reynolds
number, non-locality and wall effects into the modeling basis. With such
distinct feature injection, significant improvement of model prediction was
observed. [Beck et al. (2019)] proposed a novel data-driven strategy for
turbulence modeling for LES using artificial neural networks.

The methodology used in many of the studies discussed till now uses large high
fidelity datasets obtained from DNS or LES simulations[Parashar et al.
(2020)]. This data is used to train and validate machine learning based
models, that can vary from deep learning based models or ensembled meta-
models. The training of this model involves the inference of optimal
coefficients for the closure of the turbulence model. In these studies, the
turbulence model form corresponds to 2-equation eddy viscosity based models
(EVM) or Algebraic Reynolds Stress Models (ARSM). While these studies have had
considerable success, this outlined methodology may be hamstrung by the
discordance between the high fidelity of datasets and the fidelity of the
baseline model form. As an illustration DNS data is able to replicate the high
degree of anisotropy in turbulent flows. However eddy viscosity based models
are not capable of replicating a high degree of turbulence anisotropy due to
the linear eddy viscosity hypothesis[Mishra et al. (2019)]. The eddy viscosity
hypothesis assumes that turbulence anisotropy is a linear function of the
instantaneous mean strain rate. Thus the anisotropy predicted by any eddy
viscosity based model is forced to lie on the plane strain segment on the
barycentric triangle [Edeling et al. (2018)]. As a result the information
regarding the turbulence anisotropy in the high fidelity data is ignored by
the machine learning algorithm, due to the model form of the baseline eddy
viscosity closure. Additionally the high fidelity data reflects the
complicated dependence of turbulent statistics on streamline curvature and the
mean rate of rotation. This is important information to model the mean and
frame rotation effects on turbulence evolution. However eddy viscosity based
closures model the Reynolds stresses as functions of the mean rate of strain
only[Pope (2001)]. Thus the information regarding the effect of streamline
curvature or frame rotation is ineffectual because of the form of the baseline
model [Mishra and Girimaji (2019)]. Algebraic Reynolds Stress Models assume
that turbulent diffusive and convective fluxes are relatively small
(explicitly that the flow is source dominated [Gatski and Speziale (1993)]).
This is a restrictive assumption and is not valid for most turbulent flows. In
light of these illustrations the adoption of a different baseline model
formulation, flexible enough to take advantage of the high fidelity
information in the data may be advisable. An alternative that can meet these
requirements is the Reynolds Stress Modeling approach. Instead of assuming the
form of relationship between the Reynolds stresses and the mean gradient, the
Reynolds stress models use the Reynolds Stress Transport Equations as
evolution equations for components of turbulent anisotropy. Explicit
computation of the evolution of components of turbulent anisotropy results in
improved representation of anisotropy. Reynolds stress models can account for
directional effects of the Reynolds stresses. Reynolds stress models can
represent the limiting states of turbulent flow, for example the return to
isotropy of turbulence in decaying turbulent flows, and the Rapid Distortion
Limit where turbulence behaves like an elastic medium (instead of a viscous
medium). Because of the separate modeling of the turbulent transport
processes, Reynolds stress models account for effects of flow stratification,
buoyancy, streamline curvature, etc. Thus in a machine learning framework,
utilizing Reynolds Stress Models as the baseline model can enable the
algorithm to exploit a substantially higher proportion of physics and
information in the high fidelity data.

However there has been little research to extend the potential of Reynolds
Stress Modeling approach by utilizing machine learning algorithms. This is a
central novelty of this investigation.

In the Reynolds Stress Modeling approach separate models are formulated for
the terms in the Reynolds Stress Transport Equation, where each such term
represents a different turbulence transport process. These transport processes
include turbulent diffusion, rotational effects, rate of dissipation and the
pressure strain correlation. While high fidelity models for all these terms
are important, accurate and robust modeling of the pressure strain correlation
term is a long standing challenge in turbulence modeling. The pressure strain
correlation term represents physics responsible for the transfer of energy
between different components of the Reynolds stress tensor [Mishra and
Girimaji (2014)]. It is responsible for the non-local interactions in
turbulent flows, the initiation of instabilities in rotation dominated flows,
the return to isotropy observed in decaying flows, etc [Mishra and Girimaji
(2013)]. While classical models have been developed for the pressure strain
correlation term, such physics driven models have many limitations in their
ability to account for streamline curvature effects, realizability
requirements, their performance in complex engineering flows [Mishra and
Girimaji (2017)]. In this work we have modeled the pressure strain correlation
of turbulence using three different machine learning approaches, those are
artificial neural network, random forest and gradient boosted decision trees.
The ML models developed with above mentioned algorithms are trained and tested
for DNS data of turbulent channel flow at different Reynolds numbers. We have
grouped the data-sets in 4 different combination to perform 4 distinct
training and testing of the ML models. The Bayesian hyper-parameter
optimization was also used to find the best hyper-parameters of the ML models
that perform well for unknown prediction problems. One main advantage of
hyper-parameter optimization is that it reduces the chance of over-fitting in
tree based algorithms and also enhance the generalizability of the model.

The remainder of the article is summarized as follows: Section 2 provides a
background on the turbulence modeling framework. In Section 3 the limitations
of the turbulence models in predicting complex flow fields is discussed.
Section 4 describes the different machine learning algorithms used in modeling
the turbulence parameters. In section 5 the detailed methods of optimization
of the hyper-parameters of the ML models are discussed. In the subsequent
sections the training-testing of the models and feature importance are
discussed followed by conclusion and future scope.

##  2 Reynolds stress transport modeling framework

Different fidelities of turbulence simulation approaches can be differentiated
based on the scales of motion that are directly simulated based on their
governing equations and the scales of motion that are modeled. For example, in
Direct Numerical Simulation, all scales of motion are computed. In Reynolds
Averaged Navier Stokes (RANS) approaches, all turbulent scales of motion are
modeled. RANS closures can carry out this modeling of turbulent scales of
motion in different manner. One equation models may use the concept of mixing
length to model turbulent scales of motion. Two-equation models like the
$k-\epsilon$ and the $k-\omega$ models use the eddy viscosity hypothesis to
model turbulent motions. Reynolds Stress Models are based off the Reynolds
Stress Transport Equations. Here, each turbulent transport process is present
as a separate term in the evolution equations. In Reynolds stress transport
modeling, such transport terms such as pressure strain correlation and
dissipation tensors can be modeled using data driven approaches. In this
section we discuss the traditional method of modeling of the pressure strain
correlation starting from the basic governing equations.

The Reynolds stress transport model[Mishra and Girimaji (2010), Panda and
Warrior (2018), Panda et al. (2017)] does not rely upon any ad hoc definitions
of turbulent stresses in terms of strain field, rather the stress field is
directly computed from the modeled transport equation of Reynolds stress
components. The Reynolds stress transport equation has four terms in its right
side, those are production, transport, dissipation and the pressure strain
correlation term. The pressure strain correlation term mainly accounts for the
complex flow features resulting from the stream line curvature and flow
separation. The Reynolds stress model much more reliable and accurate than its
two equation counter parts. The Reynolds stress transport equation has the
form[Panda (2019)]:

$\displaystyle\partial_{t}\overline{u_{i}u_{j}}+U_{k}\frac{\partial\overline{u_{i}u_{j}}}{\partial
x_{k}}=P_{ij}-\frac{\partial T_{ijk}}{\partial
x_{k}}-\epsilon_{ij}+\phi_{ij},$ $\displaystyle\mbox{where},$ $\displaystyle
P_{ij}=-\overline{u_{k}u_{j}}\frac{\partial U_{i}}{\partial
x_{k}}-\overline{u_{i}u_{k}}\frac{\partial U_{j}}{\partial x_{k}},$
$\displaystyle
T_{kij}=\overline{u_{i}u_{j}u_{k}}-\nu\frac{\partial\overline{u_{i}u_{j}}}{\partial{x_{k}}}+\delta_{jk}\overline{u_{i}\frac{p}{\rho}}+\delta_{ik}\overline{u_{j}\frac{p}{\rho}}$
$\displaystyle,\epsilon_{ij}=-2\nu\overline{\frac{\partial u_{i}}{\partial
x_{k}}\frac{\partial u_{j}}{\partial x_{k}}},$
$\displaystyle\phi_{ij}=\overline{\frac{p}{\rho}(\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})}$

$P_{ij}$ denotes the production of turbulence, $T_{ijk}$ is the diffusive
transport, $\epsilon_{ij}$ is the dissipation rate tensor and $\phi_{ij}$ is
the pressure strain correlation. The pressure fluctuations are governed by a
Poisson equation:

$\frac{1}{\rho}{\nabla}^{2}(p)=-2\frac{\partial{U}_{j}}{\partial{x}_{i}}\frac{\partial{u}_{i}}{\partial{x}_{j}}-\frac{\partial^{2}u_{i}u_{j}}{\partial
x_{i}\partial x_{j}}$ (1)

The fluctuating pressure term is split into a slow and rapid pressure term
$p=p^{S}+p^{R}$. Slow and rapid pressure fluctuations satisfy the following
equations

$\frac{1}{\rho}{\nabla}^{2}(p^{S})=-\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}{(u_{i}u_{j}-\overline{u_{i}u_{j}})}$ (2)
$\frac{1}{\rho}{\nabla}^{2}(p^{R})=-2\frac{\partial{U}_{j}}{\partial{x}_{i}}\frac{\partial{u}_{i}}{\partial{x}_{j}}$
(3)

It can be seen that the slow pressure term accounts for the non-linear
interactions in the fluctuating velocity field and the rapid pressure term
accounts for the linear interactions[Mishra and Girimaji (2015)]. The pressure
strain correlation is modeled using rational mechanics approach. The rapid
term can be modeled as[Pope (2000)]

$\phi_{ij}^{R}=4k\frac{\partial{U}_{l}}{\partial{x_{k}}}(M_{kjil}+M_{ikjl})$
(4)

where,

$M_{ijpq}=\frac{-1}{8\pi
k}\int\frac{1}{r}\frac{\partial^{2}R_{ij}(r)}{\partial r_{p}\partial r_{p}}dr$
(5)

where, $R_{ij}(r)=\langle u_{i}(x)u_{j}(x+r)\rangle$ For homogeneous
turbulence the complete pressure strain correlation can be written as

$\phi_{ij}=\epsilon
A_{ij}(b)+kM_{ijkl}(b)\frac{\partial\overline{v}_{k}}{\partial{x_{l}}}$ (6)

The most general form of slow pressure strain correlation is given by

$\phi^{S}_{ij}=\beta_{1}b_{ij}+\beta_{2}(b_{ik}b_{kj}-\frac{1}{3}II_{b}\delta_{ij})$
(7)

Established slow pressure strain correlation models including the models of
[Rotta (1951)] and [Sarkar and Speziale (1990), Panda et al. (2017), Warrior
et al. (2014)] use this general expression. Considering the rapid pressure
strain correlation, the linear form of the model expression is

$\displaystyle\frac{\phi^{R}_{ij}}{k}=C_{2}S_{ij}+C_{3}(b_{ik}S_{jk}+b_{jk}S_{ik}-$
(8)
$\displaystyle\frac{2}{3}b_{mn}S_{mn}\delta_{ij})+C_{4}(b_{ik}W_{jk}+b_{jk}W_{ik})$

Here $b_{ij}=\frac{\overline{u_{i}u_{j}}}{2k}-\frac{\delta_{ij}}{3}$ is the
Reynolds stress anisotropy tensor, $S_{ij}$ is the mean rate of strain and
$W_{ij}$ is the mean rate of rotation. Rapid pressure strain correlation
models like the models of [Mishra and Girimaji (2017), Speziale et al. (1991),
Panda and Warrior (2018)] use this general expression. For this work the
Reynolds stress model of [Speziale et al. (1991)] is used which has the form:

$\displaystyle\phi_{ij}^{(R)}=(C_{1}-C_{1}^{*}II^{0.5})KS_{ij}+$ (9)
$\displaystyle C_{2}K(b_{ik}S_{jk}+b_{jk}S_{ik}-2/3b_{mn}S_{mn}\delta_{ij})$
$\displaystyle+C_{3}K(b_{ik}W_{jk}+b_{jk}W_{ik})$

(a)

(b)

(c)

Figure 1: Evolution of turbulence kinetic energy for rotating shear flows, a)
$W/S=0$, b) $W/S=0.25$ and c) $W/S=0.5$, Symbols LES data, dashed lines
predictions of LRR model, dotted lines predictions of SSG model.

The closure coefficients are taken as $C_{1}=0.8$, $C_{1}^{*}=1.3$,
$C_{2}=1.25$ and $C_{3}=0.4$.

##  3 Limitation of the existing pressure strain correlation model in the
prediction of complex turbulent flows

Most of the Reynolds stress transport models available in literature have some
tuned coefficients, those are obtained from the calibration of set of
equations against few standard experimental/DNS data-sets like flow in
channels, grid turbulence and grid turbulence, so when those models are tested
for complex flow configurations produce unrealistic results. The Reynolds
stress transport models are also associated with another type of uncertainty,
i.e. those can not accurately capture the non-local nature of flow(The flow at
one point may be affected by flow physics of upstream points), since the model
form does not have any terms for accommodation of the non-local nature of
flow. In order to overcome these limitations of the Reynolds stress transport
models data driven models of the pressure strain correlation model can be
developed using machine learning approaches.

We have performed numerical simulations to check the predictive capability of
different pressure strain correlation models in rotation dominated flow
fields. It is noticed that, there is a strong disparity between the model
predictions and the LES data particularly at higher values of rotation to
strain ratio. The results of those comparisons are presented in fig.1. fig.1
a, b and c correspond to rotation to strain ratio of $0$, $0.25$ and $0.5$
respectively. The symbols in the fig.1 represent the LES results of[Bardino et
al. (1983)].

##  4 Description of the data set

The dataset used in the modeling of the pressure strain correlation of
turbulence was obtained from the Oden Institute Turbulence File Server[Lee and
Moser (2015)], in which DNS simulations were performed for flow in channels
with friction Reynolds number ranging from 550 to 5200. The DNS simulations
were performed with a B-spline collocation method in the wall-normal direction
and for the stream-wise and span-wise directions Fourier-Galerkin method was
used. The detailed information on the on the methods employed in the DNS of
channel flow is available in [Lee and Moser (2015)]. We have used four
distinct data-sets for training and testing of the ML models. The four
training and testing data-sets are presented in table1.

Case | Training Set | Testing Set
---|---|---
1 | $Re_{\lambda}=550,1000,2000$ | $Re_{\lambda}=5200$
2 | $Re_{\lambda}=550,1000,5200$ | $Re_{\lambda}=2000$
3 | $Re_{\lambda}=550,2000,5200$ | $Re_{\lambda}=1000$
4 | $Re_{\lambda}=1000,2000,5200$ | $Re_{\lambda}=550$
Table 1: Four training and test cases for the turbulent channel flow

##  5 Physics based input features for ML models

The turbulence parameter e.g. pressure strain correlation can be accurately
modeled using suitable input features. The most important features for the
modeling of the pressure strain correlation can be the Reynolds stress
anisotropy, dissipation, velocity gradient and the turbulence kinetic energy.
The input features must be wisely chosen so that there is no over-fitting.
Secondly this ensures that physics based constraints are met in the final
model. For instance due to Galilean in-variance we should ensure that the
features in the modeling basis also obey this requirement. The functional
mapping for the ML models can be written as:

$\displaystyle\phi_{uv}=f_{1}(b_{uv},\epsilon,\frac{du}{dy},k)$ (10)

We have used the formula:
$\alpha^{*}=\frac{\alpha-\alpha_{min}}{\alpha_{max}-\alpha_{min}}$ for
normalizing the inputs to the ML models, so that those will be in the range 0
and 1. This avoids clustering of training in one direction and enhance
convergence in the training.

##  6 Data driven turbulence modeling with machine learning

In turbulence modeling frameworks such as eddy viscosity and Reynolds stress
transport models the Reynolds stress and pressure strain correlation terms
plays major role and with accurate modeling of those two parameters, the model
form uncertainties can be eradicated. In this work, we have modeled above
mentioned terms using three different machine learning algorithms, those are
deep neural networks, Random forests and gradient boosted trees. In turbulence
modeling and fluid dynamics applications deep neural networks are widely used.
Next, we have used random forests in the models, those have the capability to
check the feature importance, i.e. which of the input features has dominant
correlation with the output. We also have considered gradient boosted trees in
the modeling of the turbulence parameters. The basic difference between RF and
GB trees is that RF creates random training samples from full training set
based on Bootstrap Aggregating(Bagging) and a weak learner is trained in
parallel and in GB trees ensemble methods consist in fitting several weak
learners sequentially as shown in figure. 2.

(a)

(b)

Figure 2: a) Bagging vs b) Boosting

###  6.1 Artificial Neural Networks

Artificial Neural Networks(ANN)(fig.3) are the machine learning systems
inspired from the biological neural networks, these are also known as multi-
layer perception(MLP). The biological neural networks(BNN) are the circuits
that carry out a specific task when activated. These are population of neurons
interconnected by synapses. Similar to BNNs, ANN/MLP have artificial neurons.
A very basic unit of a MLP is a perceptron as shown in fig.3a. The input data
from a (l-1)th layer are multiplied by a weight $W$, which are linearly
combined and allowed to pass through a non-linear activation function $\eta$:

$\displaystyle q_{i}^{l}=\eta(\sum_{j}W_{ij}^{l}q_{j}^{l-1}).$ (11)

A MLP can be constructed by combining a number of perceptrons. A typical MLP
is shown in fig.3b. The weights in a MLP are optimized to minimize a cost
function $E$ by back propagation. There are several optimization techniques by
which the weights can be calculated. Those are gradient descent, Quasi-
Newton, Stochastic Gradient Descent or the Adaptive moment Estimation etc. The
activation functions, which can be used in MLP are sigmoid
$\eta(\beta)=1/(1+e^{-\beta})$, hyperbolic tangent(tanh)
$\eta(\beta)=(e^{-\beta}-e^{-\beta}).(e^{-\beta}+e^{-\beta})$ and RELU
$\eta(\beta)=max[0,\beta]$.

(a)

(b)

Figure 3: a) A perceptron b) An artificial neural network(Multi layer
perception)

###  6.2 Random forests(RF)

The random forest algorithm was proposed by [Breiman (2001)]. The basic
building block of a random forest is a decision tree. The structure of the
decision tree is shown in fig.4a. The boxes in the decision tree represent
group of features and data. The decision tree seeks if/then/else rules for
obtaining the desired output. Random forest(RF)(fig.4b) regression combines
performance of multiple decision trees for predicting a output variable. It is
a assemble learning technique, works on the concept of bagging method. In RF
regression trees are constructed using a subset of random sample drawn with
replacement from training data. In RF to depict the growth of the tress random
vectors are generated and the trees are not allowed to prun. To perform
splitting of the dataset, a random combination of features is selected at each
and every node.

Figure 4: Architecture of of a random forest, a)a single decision tree b)a
random forest

###  6.3 Gradient boosted decision trees(GBDT)

Similar to random forest algorithm, the building block of GBDT are the
decision trees. Random forests utilizes method of bagging to combine many
decision trees to create an ensemble. Bagging simply means combining the
decision trees in parallel. However, boosting means combining decision trees
in series to achieve strong learner. The decision trees are the weak learners.
The boosting algorithms learn slowly, since trees are added sequentially. Each
tree in the boosting trees focus on errors from previous one, making the
boosting algorithm an efficient and accurate model. Boosting process is slow,
since the trees in the GBDT are added sequentially.

##  7 Hyper-parameter optimisation

The hyper-parameter optimization is the process of finding the best hyper-
parameters for a machine learning algorithm that returns best performance
while measured on a validation set. The hyper-parameters are turned by
application engineers during model developments. There are mainly two methods
of hyper-parameter tuning, those are either performed manually or done by some
optimization algorithm. The automated algorithms are grid search, random
search and Bayesian optimization. The grid search and random search not
informed by past evaluations. The Bayesian optimization keep track of past
evaluations, hence we have considered the Bayesian optimization for automated
hyper-parameter tuning.

###  7.1 Manual search

The hyper-parameters of different machine learning algorithms, which need to
finely tuned for improving the predictive capability of any machine learning
model. For neural networks the important hyper-parameters are number of layers
and number of neurons in each layer. Similarly for random forest and gradient
boosted decision trees the hyper-parameters are number of decision trees and
maximum depth of the trees. Few other parameters are also important in tree
algorithms, those are minimum sample leaf and learning rate for gradient
boosting decision trees. However fine tuning of all such parameters using
manual search method is tedious task and also there is also the chance of over
fitting. In manual search, we have considered $R^{2}$ as the criteria for
finding the hyper-parameters. By manual search the optimal hyper-parameters
are tuned as follows: for neural network: 5 layers with 10 neurons in each
layer and for both RF and GBDT the number of decision trees and max depth are
taken as 5 and 10 respectively. By considering above mentioned hyper-
parameters the $R^{2}$ and MSE values for test data are calculated as follows
for case 4: ANN: $0.9897$, $9.97*10^{-6}$, RF: $0.9841$, $1.66*10^{-5}$, GBDT:
$0.981$, $1.99*10^{-5}$.

###  7.2 Bayesian optimization

In Bayesian optimization, the parameter space is sampled and a gradient
boosted tree was constructed and trained. The performance of the model is
evaluated using validation data. The optimal hyper-parameters are those, that
yield the lowest error on the validation data. The model validation error was
treated as a sample from Gaussian process at each hyper-parameter setting. In
contrast to other optimization techniques, the Bayesian optimization requires
relatively lesser number of model evaluations, since at each step the
information from all previous states is utilized to inform the GP model of the
validation error. This method is highly useful, when computational cost of
evaluation of objective function is very high. The detailed steps in Bayesian
optimization are as follows:

1\. Build a surrogate probability model of the objective function
2\. Finding the hyper-parameters that performs best on the surrogate model
3\. Application of these hyper-parameters to the true objective function
4\. Updating of the surrogate model incorporating the new results
5\. The steps 2 and 4 must be repeated until max iterations or time is reached

(a)  Figure 5: Histogram of the data variability for training data in case 1, There are very few training samples for high values of pressure strain correlation. ML models | Parameter I | parameter II
---|---|---
MLP | N. of neurons | Layers
RF | N. of decision trees | Depth
GBDT | N. of decision trees | Depth
Table 2: The hyper-parameters of the machine learning models. (N. stands for
number).

(a)

(b)

(c)

(d)

Figure 6: Loss convergence in MLP training. a) case 4 b) case 3 c) case 2 d) case 1. Solid line and dashed line correspond to training and validation losses respectively. | Case 1 | Case 2 | Case 3 | Case 4
---|---|---|---|---
Models | Train $R^{2}$ | Test $R^{2}$ | Train $R^{2}$ | Test $R^{2}$ | Train $R^{2}$ | Test $R^{2}$ | Train $R^{2}$ | Test $R^{2}$
MLP | 0.967 | 0.978 | 0.978 | 0.984 | 0.921 | 0.891 | 0.992 | 0.892
RF | 0.971 | 0.978 | 0.949 | 0.966 | 0.962 | 0.988 | 0.984 | 0.885
GBDT | 0.982 | 0.871 | 0.968 | 0.958 | 0.951 | 0.955 | 0.982 | 0.871
Table 3: $R^{2}$ of $\phi_{uv}^{*}$ predictions by training-prediction cases

The hyper-parameters mentioned in table 2 are only a limited set of important
parameters, through which the model can be developed. The models with those
limited tuned parameters may not perform well for unknown prediction data-sets
on which the model is trained. To enhance the predictive capability of any
machine learning model at least few other hypo-parameters must be considered
in the model development stage. However, with increase in number of hyper-
parameters, manual tuning may be a tedious process and advanced optimization
techniques like Bayesian optimisation must be used for fine turning of the
hyper-parameters. For brevity, we have only considered the hyper-parameter
optimization of gradient boosted regression tree using Bayesian optimization.
In addition to depth and number of decision trees, three other hypo-parameters
are considered, those are max features, minimum simple leaf and minimum sample
split. From Bayesian optimization, the optimal hyper-parameters(maximum depth,
Number of estimators, minimum sample split, minimum sample leaf, maximum
features) are (66, 63, 9, 5, sqrt). The maximum depth correspond to highest
depth of the individual regression tree. The best value of this parameter
depends on the interaction of the input variables. Number of estimators are
the number of boosting stages used in the model development. The larger is the
number of the estimators, better is the performance of the the model. The
minimum sample leaf is the minimum number of minimum number of samples
required to split an internal node. The other terms has their usual meaning. A
detailed discussion of all these hyper-parameters of GBDT is available in
[Pedregosa et al. (2011)]. The $R^{2}$ values of the optimized GBDT was found
to be 0.9915 and 0.9076 for the training and testing respectively for case 4.
The predictions of the GBDT with optimized hyper-parameters for unknown flow
cases are summarized in section 9.

##  8 Training of the ML models

All the three ML models are trained for the DNS data of turbulent channel flow
for four different test cases. The four different cases are shown in table 1.
Since the turbulence statistics are available for four different friction
Reynolds numbers, in each training and testing phase, we have considered data
for three friction Reynolds numbers as training and the other was considered
for testing. We have considered two different approaches for finding the
hyper-parameters of the ML models as discussed in section 6. Using the hyper-
parameters obtained from manual search, the $R^{2}$ scores obtained from
different ML models are presented in table 3. For case 1, the ML models were
trained with data for friction Reynolds numbers 550, 1000, 2000 and tested for
5200. The maximum value of train $R^{2}$ is obtained from the results of case
4 with MLP. The loss convergence history of MLP for the four different cases
are shown in figure.6. In most of the training cases, the $R^{2}$ magnitude is
greater than $0.95$, which signifies a good correlation between the actual and
predicted value of pressure strain term ($\phi$) for training samples (A
$R^{2}$ value of 0 and 1 correspond to no-correlation and perfect correlation
between the actual and predicted values of the pressure strain term). In
figure.7, the results of DNS simulation samples versus data driven model
predictions of the pressure strain correlation are presented. There is very
good match between the DNS results and DNS model predictions. There are very
few points in the data-set with high values of the pressure strain
correlation, this signifies the fact that, the relative inaccuracy of the
model predictions for high values of the pressure strain correlation is due to
very few such training data samples(figure5).

(a)

(b)

(c)

Figure 7: DNS simulation samples versus data driven model predictions of the
pressure strain correlation for the ML models training for case 4, a) Random
forest b) Gradient boosted trees c) Neural network.

(a)

(b)

(c)

(d)

Figure 8: Prediction of the pressure strain correlation using different
machine learning algorithms.Solid line, dashed line and dotted line represent
predictions of neural network, gradient boosted trees and random forests
respectively. a) case 4, b) case 3, c) case 2, d) case 1.

##  9 Testing of the trained ML models

The trained models were tested against testing data for all the four cases as
discussed in table.1. The results of testing of the ML models are presented in
8. From all four figures of figure 8 it is clear that, the GBDT predictions
are better in comparison to the predictions of MLP and RF. The neural networks
failed to predict the pressure strain correlation for $y+$ values greater than
25 for case 4. The testing results are presented in table 3. In most of the
testing cases the $R^{2}$ value was found to be greater than $0.95$, but in
few cases like case 4, the $R^{2}$ value falls below $0.9$ for the testing
case. This is because of over-fitting of the models with the training data.
The over-fitting problem is common in random forests. The advantage of GBDT
over MLP is that the former needs less number of weights in developing the
correlation between the inputs and the output. Since the ultimate aim of
turbulence modeler is to apply the ML based based into CFD solver, models with
less number of weights are always preferred. Since model with large number of
weights, may end with divergence in CFD solver. So we have tried to enhance
the predictive capability of the GBDT with Bayesian optimization of hyper-
parameters as discussed in section 7. The results of predicted values $\phi$
with the optimized GBDT is presented in 9. The solid line in the figure
correspond to the predictions of the optimized GBDT and the dashed line
correspond the prediction of the GBDT with the hyper-parameters obtained from
manual tuning. from the figure it is clear that, the predictions of the
optimized GBDT matches well with the DNS results. Finally, we have tested the
GBDT model against a fully unknown flow e.g. Couette flow in channels. The
model predictions are presented in figure.10. It is noticed that the GBDT
predictions are matching well the DNS results of [Lee and Moser (2018)].

(a)  Figure 9: Prediction of the pressure strain correlation for case 4, using
the optimized GBDT. The hyper-parameters of the GBDT were tuned using Bayesian
optimization. (a)  Figure 10: Prediction of the pressure strain correlation
for turbulent Couette flow at $Re_{\lambda}=550$.

##  10 Conclusions

In this article, we discuss rationale for the application of machine learning
with high-fidelity turbulence data to develop models at the level of Reynolds
stress transport modeling. Then, we ascertain the efficacy of different
machine learning algorithms at creating surrogate models for the pressure
strain correlation. We have modeled the pressure strain correlation of
turbulence using three different machine learning approaches, those are Random
forests, gradient boosted tress and artificial neural networks. The input
features to the ML models were chosen from the traditional modeling basis of
the pressure strain correlation, those are mean strain, turbulence kinetic
energy, Reynolds stress anisotropy and turbulence dissipation. The ML models
were trained and tested for DNS data of turbulent channel flow at different
friction Reynolds numbers. The optimal values of ML model hyper-parameters
were optimized using manual search and Bayesian optimization approaches. The
feature importances of different input features of the random forest were
obtained using mean decrease in impurity method. It was noticed that the mean
strain has larger correlation with the pressure strain term. The ML models
developed by using data-driven approaches can be utilized in computational
fluid dynamics solvers for improved flow predictions in turbulent channel
flows. in future course of work, for achieving universality, large amount of
DNS or experimental data can be fed to the ML models during training and
development.

## References

  * Bardino et al. (1983) Bardino, J., J. H. Ferziger, and W. C. Reynolds (1983).  Improved turbulence models based on large eddy simulation of homogeneous, incompressible turbulent flows.  Report Stanford Univ..
  * Beck et al. (2019) Beck, A., D. Flad, and C.-D. Munz (2019).  Deep neural networks for data-driven les closure models.  Journal of Computational Physics 398, 108910.
  * Breiman (2001) Breiman, L. (2001).  Random forests.  Machine learning 45(1), 5–32.
  * Duraisamy et al. (2019) Duraisamy, K., G. Iaccarino, and H. Xiao (2019).  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics 51, 357–377.
  * Edeling et al. (2018) Edeling, W. N., G. Iaccarino, and P. Cinnella (2018).  Data-free and data-driven rans predictions with quantified uncertainty.  Flow, Turbulence and Combustion 100(3), 593–616.
  * Fang et al. (2020) Fang, R., D. Sondak, P. Protopapas, and S. Succi (2020).  Neural network models for the anisotropic reynolds stress tensor in turbulent channel flow.  Journal of Turbulence 21(9-10), 525–543.
  * Gatski and Speziale (1993) Gatski, T. B. and C. G. Speziale (1993).  On explicit algebraic stress models for complex turbulent flows.  Journal of fluid Mechanics 254, 59–78.
  * Huijing et al. (2021) Huijing, J. P., R. P. Dwight, and M. Schmelzer (2021).  Data-driven rans closures for three-dimensional flows around bluff bodies.  Computers & Fluids, 104997.
  * Kaandorp and Dwight (2020) Kaandorp, M. L. and R. P. Dwight (2020).  Data-driven modelling of the reynolds stress tensor using random forests with invariance.  Computers & Fluids, 104497.
  * Lee and Moser (2015) Lee, M. and R. D. Moser (2015).  Direct numerical simulation of turbulent channel flow up to $re_{t}au=5200$.  Journal of Fluid Mechanics 774, 395–415.
  * Lee and Moser (2018) Lee, M. and R. D. Moser (2018).  Extreme-scale motions in turbulent plane couette flows.  Journal of Fluid Mechanics 842, 128–145.
  * Ling et al. (2016) Ling, J., A. Kurzawski, and J. Templeton (2016).  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 155–166.
  * Luan and Dwight (2020) Luan, Y. and R. P. Dwight (2020).  Influence of turbulence anisotropy on rans predictions of wind-turbine wakes.  In Journal of Physics: Conference Series, Volume 1618, pp. 062059. IOP Publishing.
  * Maulik et al. (2018) Maulik, R., O. San, A. Rasheed, and P. Vedula (2018).  Data-driven deconvolution for large eddy simulations of kraichnan turbulence.  Physics of Fluids 30(12), 125109.
  * Mishra et al. (2019) Mishra, A. A., K. Duraisamy, and G. Iaccarino (2019).  Estimating uncertainty in homogeneous turbulence evolution due to coarse-graining.  Physics of Fluids 31(2), 025106.
  * Mishra and Girimaji (2019) Mishra, A. A. and S. Girimaji (2019).  Linear analysis of non-local physics in homogeneous turbulent flows.  Physics of Fluids 31(3), 035102.
  * Mishra and Girimaji (2010) Mishra, A. A. and S. S. Girimaji (2010).  Pressure–strain correlation modeling: towards achieving consistency with rapid distortion theory.  Flow, turbulence and combustion 85(3-4), 593–619.
  * Mishra and Girimaji (2013) Mishra, A. A. and S. S. Girimaji (2013).  Intercomponent energy transfer in incompressible homogeneous turbulence: multi-point physics and amenability to one-point closures.  Journal of Fluid Mechanics 731, 639–681.
  * Mishra and Girimaji (2014) Mishra, A. A. and S. S. Girimaji (2014).  On the realizability of pressure-strain closures.  Journal of fluid mechanics 755, 535.
  * Mishra and Girimaji (2015) Mishra, A. A. and S. S. Girimaji (2015).  Hydrodynamic stability of three-dimensional homogeneous flow topologies.  Physical Review E 92(5), 053001.
  * Mishra and Girimaji (2017) Mishra, A. A. and S. S. Girimaji (2017).  Toward approximating non-local dynamics in single-point pressure–strain correlation closures.  Journal of Fluid Mechanics 811, 168–188.
  * Panda (2019) Panda, J. (2019).  A review of pressure strain correlation modeling for reynolds stress models.  Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 0954406219893397.
  * Panda (2020) Panda, J. (2020).  A reliable pressure strain correlation model for complex turbulent flows.  J Appl Fluid Mech 13, 1167–1178.
  * Panda et al. (2020) Panda, J., K. Sasmal, S. Maity, and H. Warrior (2020).  A simple nonlinear eddy viscosity model for geophysical turbulent flows.  J. Appl. Fluid Mech. 13, 1167–1178.
  * Panda and Warrior (2018) Panda, J. and H. Warrior (2018).  A representation theory-based model for the rapid pressure strain correlation of turbulence.  Journal of Fluids Engineering 140(8).
  * Panda et al. (2017) Panda, J., H. Warrior, S. Maity, A. Mitra, and K. Sasmal (2017).  An improved model including length scale anisotropy for the pressure strain correlation of turbulence.  ASME Journal of Fluids Engineering 139(4), 044503\.
  * Parashar et al. (2020) Parashar, N., B. Srinivasan, and S. S. Sinha (2020).  Modeling the pressure-hessian tensor using deep neural networks.  Physical Review Fluids 5(11), 114604.
  * Parish and Duraisamy (2016) Parish, E. J. and K. Duraisamy (2016).  A paradigm for data-driven predictive modeling using field inversion and machine learning.  Journal of Computational Physics 305, 758–774.
  * Pedregosa et al. (2011) Pedregosa, F., G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay (2011).  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Pope (1975) Pope, S. (1975).  A more general effective-viscosity hypothesis.  Journal of Fluid Mechanics 72(2), 331–340.
  * Pope (2000) Pope, S. (2000).  Turbulent Flows.  New York: Cambridge University Press.
  * Pope (2001) Pope, S. B. (2001).  Turbulent flows.
  * Rotta (1951) Rotta, J. (1951).  Statistische theorie nichthomogener turbulenz.  Z. Phys. 129, 547–572.
  * Sarkar and Speziale (1990) Sarkar, S. and C. G. Speziale (1990).  A simple nonlinear model for the return to isotropy in turbulence.  Physics of Fluids A: Fluid Dynamics 2(1), 84–93.
  * Schmelzer et al. (2020) Schmelzer, M., R. P. Dwight, and P. Cinnella (2020).  Discovery of algebraic reynolds-stress models using sparse symbolic regression.  Flow, Turbulence and Combustion 104(2), 579–603.
  * Singh et al. (2017) Singh, A. P., S. Medida, and K. Duraisamy (2017).  Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils.  AIAA journal 55(7), 2215–2227.
  * Speziale et al. (1991) Speziale, C. G., S. Sarkar, and T. B. Gatski (1991).  Modelling the pressure–strain correlation of turbulence: an invariant dynamical systems approach.  Journal of fluid mechanics 227, 245–272.
  * Tracey et al. (2015) Tracey, B. D., K. Duraisamy, and J. J. Alonso (2015).  A machine learning strategy to assist turbulence model development.  In 53rd AIAA aerospace sciences meeting, pp. 1287.
  * Versteeg and Malalasekera (2007) Versteeg, H. K. and W. Malalasekera (2007).  An introduction to computational fluid dynamics: the finite volume method.  Pearson education.
  * Wang et al. (2017) Wang, J.-X., J.-L. Wu, and H. Xiao (2017).  Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data.  Physical Review Fluids 2(3), 034603.
  * Warrior et al. (2014) Warrior, H., S. Mathews, S. Maity, and K. Sasmal (2014).  An improved model for the return to isotropy of homogeneous turbulence.  ASME Journal of Fluids Engineering 136(3), 034501\.
  * Weatheritt and Sandberg (2016) Weatheritt, J. and R. Sandberg (2016).  A novel evolutionary algorithm applied to algebraic modifications of the rans stress–strain relationship.  Journal of Computational Physics 325, 22–37.
  * Weatheritt and Sandberg (2017) Weatheritt, J. and R. Sandberg (2017).  The development of algebraic stress models using a novel evolutionary algorithm.  International Journal of Heat and Fluid Flow 68, 298–318.
  * Wu et al. (2020) Wu, J., J. Li, X. Qiu, and Y. Liu (2020).  A comparative analysis of multi-machine learning algorithms for data-driven rans turbulence modelling.  In Journal of Physics: Conference Series, Volume 1684, pp. 012043. IOP Publishing.
  * Wu et al. (2018) Wu, J.-L., H. Xiao, and E. Paterson (2018).  Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework.  Physical Review Fluids 3(7), 074602.
  * Zhao et al. (2020) Zhao, Y., H. D. Akolekar, J. Weatheritt, V. Michelassi, and R. D. Sandberg (2020).  Rans turbulence model development using cfd-driven machine learning.  Journal of Computational Physics 411, 109413.
  * Zhu et al. (2019) Zhu, L., W. Zhang, J. Kou, and Y. Liu (2019).  Machine learning methods for turbulence modeling in subsonic flows around airfoils.  Physics of Fluids 31(1), 015105.


# Gas Source Localization Using Physics-Guided Neural Networks

Victor Prieto Ruiz, Patrick Hinsen,
Thomas Wiedemann, and Dmitriy Shutin
_Institute of Communications and Navigation,_
_German Aerospace Center (DLR),_ Oberpfaffenhofen, Germany  Constantin
Christof
This work was supported by the EU Project TEMA. The TEMA project has received
funding from the European Commission under HORIZON EUROPE (HORIZON Research
and Innovation Actions) under grant agreement 101093003 (HORIZON-
CL4-2022-DATA-01-01). Views and opinions expressed are however those of the
author(s) only and do not necessarily reflect those of the European Union or
Commission. Neither the European Commission nor the European Union can be held
responsible for them. _School of Computation, Information and Technology_ ,
_Technical University of Munich,_
Garching, Germany

###### Abstract

This work discusses a novel method for estimating the location of a gas source
based on spatially distributed concentration measurements taken, e.g., by a
mobile robot or flying platform that follows a predefined trajectory to
collect samples. The proposed approach uses a Physics-Guided Neural Network to
approximate the gas dispersion with the source location as an additional
network input. After an initial offline training phase, the neural network can
be used to efficiently solve the inverse problem of localizing the gas source
based on measurements. The proposed approach allows avoiding rather costly
numerical simulations of gas physics needed for solving inverse problems. Our
experiments show that the method localizes the source well, even when dealing
with measurements affected by noise.

###### Index Terms:

gas source localization, robotic olfaction, physics-guided neural network,
inverse problem.

##  I Introduction

The problem of Gas Source Localization (GSL) (e.g., of a gas leak) by means of
measurements taken by mobile agents has received increasing attention in the
last decades due to rapid improvements in robotics and gas sensor technology.
While biologically inspired localization approaches, such as chemotaxis or
anemotaxis, have initially served as simple guiding principles, these
strategies have been found to be not robust enough for efficient source
localization [1, Sec. 4.1]. More recent probabilistic model-based approaches,
such as infotaxis, have proved more reliable in patchy and turbulent regimes.
However, they rely on assumptions on the source distribution, making them
highly dependent on prior knowledge [2, Sec. 6]. Further, the incorporated
models often come with high computational costs. For a detailed review of
different GSL-approaches, we refer to [1] and the references therein.

In the present work, we propose a strategy which approximates the complex gas
dispersion physics through the use of an easy-to-evaluate neural network. More
precisely, a Physics-Guided Neural Network (PGNN) [3, Secs. 1,2] is used as a
surrogate model for the advection-diffusion Partial Differential Equation
(PDE). This network is trained to learn and efficiently emulate the
complicated functional dependency between the spatial gas concentration and
the source location.

Our surrogate model approach offers several advantages: First, it is flexible
and can be applied to various types of PDEs with different boundary conditions
and parameters. Second, our surrogate model can be differentiated easily with
respect to the source location, thanks to automatic differentiation libraries.
This makes the model extremely useful in the inverse problem setting that
requires a gradient-based numerical optimization method. Third, the
computationally expensive network training can be done once offline.
Afterwards, the network can be used for solving GSL-problems with different
observations, without the need for retraining or a costly PDE-solver. This
opens up the possibility for real-time applications, even those using
hardware-constrained robots. Lastly, our method is also robust to noisy
measurements and, thus, able to handle inaccuracies caused by sensor
limitations.

##  II Methodology

###  II-A Gas Model

We assume that the gas occupies the spatial domain
$\Omega=($0\text{\,}\mathrm{km}$,$1\text{\,}\mathrm{km}$)\times($0\text{\,}\mathrm{km}$,$1\text{\,}\mathrm{km}$)$.
On the set $\Omega$, we model the gas spread by the linear advection-diffusion
PDE. For the sake of simplicity, we consider a time-invariant (equilibrium)
problem with homogeneous Dirichlet BCs. This is given by the PDE

$\displaystyle\partial_{t}u(\boldsymbol{x})=\kappa\Delta
u(\boldsymbol{x})-\boldsymbol{v}\cdot\nabla
u(\boldsymbol{x})+s(\boldsymbol{x})$ $\displaystyle=0$ $\displaystyle\text{in
}\Omega,$ (1) $\displaystyle u(\boldsymbol{x})$ $\displaystyle=0$
$\displaystyle\text{on }\partial\Omega.$

Here, $\partial\Omega$ is the boundary of $\Omega$,
$u\colon\Omega\to\mathbb{R}$ is the gas concentration, $s$ is a source term
(either a real-valued function or a measure on $\Omega$),
$\kappa=$1\text{\,}{\mathrm{km}}^{2}\text{\,}{\mathrm{h}}^{-1}$$ is the
diffusion coefficient, and
$\boldsymbol{v}=[$3\text{\,}\mathrm{km}\text{\,}{\mathrm{h}}^{-1}$,$3\text{\,}\mathrm{km}\text{\,}{\mathrm{h}}^{-1}$]^{\top}$
is the advection velocity. We assume a uniform (known) wind flow.

Figure 1:  Problem setup. The gas concentration field is shown as an example.
(A ball of radius $20\text{\,}\mathrm{m}$ around $\boldsymbol{p}$ is excluded
to avoid the singularity.) The dotted black circle represents
$D_{\mathrm{O}}$, the trajectory of the observations. The source location
$\boldsymbol{p}$ is marked by a black dot and the set $\mathcal{P}$ of
possible source positions by the blue hatched region.

Since we are interested in GSL, we assume that there is a location parameter
set
$\mathcal{P}=[$0.35\text{\,}\mathrm{km}$,$0.65\text{\,}\mathrm{km}$]\times[$0.35\text{\,}\mathrm{km}$,$0.65\text{\,}\mathrm{km}$]\subset\Omega$
which fully represents all possible points $\boldsymbol{p}\in\mathbb{R}^{2}$
where the source could be located. We shall use $s_{\boldsymbol{p}}$ and
$u_{\boldsymbol{p}}$ to denote the source term and the gas concentration field
associated with the source position $\boldsymbol{p}$, respectively. That is,
$u_{\boldsymbol{p}}$ is the solution of (1) with source term
$s_{\boldsymbol{p}}$.

For the purposes of this paper, we restrict our attention to ideal point
sources, i.e., Dirac delta distributions. Keeping with previous notation, this
means that we consider source terms of the form
$s_{\boldsymbol{p}}:=\delta_{\boldsymbol{p}}$, $\boldsymbol{p}\in\mathcal{P}$.
Note that, for such $s_{\boldsymbol{p}}$, the solutions $u_{\boldsymbol{p}}$
of (1) possess a singularity at $\boldsymbol{p}$. This means in particular
that $u_{\boldsymbol{p}}$ cannot be evaluated at $\boldsymbol{p}$.

###  II-B Source Localization Problem

In a source localization problem governed by (1), the goal is to identify the
location $\boldsymbol{p}\in\mathcal{P}$ associated with the source term
$s_{\boldsymbol{p}}$ from measurements of $u_{\boldsymbol{p}}$ taken by a
mobile robot on an observation domain $D_{\mathrm{O}}\subset\Omega$. The set
$D_{\mathrm{O}}$ can be considered as the robot trajectory in the region of
interest; cf. Fig. 1.

We suppose that the gas concentration measurements are taken at
$N_{\mathrm{O}}$ observation points
$\boldsymbol{x}_{1},...,\boldsymbol{x}_{N_{\mathrm{O}}}\in D_{\mathrm{O}}$. A
measurement $y_{j}$ at $\boldsymbol{x}_{j}$ is the gas concentration
$u_{\boldsymbol{p}}(\boldsymbol{x}_{j})$ perturbed by zero-mean additive
Gaussian noise with standard deviation $\sigma$. To avoid evaluating the
solution of (1) at its singularity, we henceforth assume that $\mathcal{P}\cap
D_{\mathrm{O}}=\emptyset$. We remark that this restriction can be avoided by
replacing the point sources $\delta_{\boldsymbol{p}}$ with a suitable
mollification.

To determine the source location based on measurements
$y_{1},...,y_{N_{\mathrm{O}}}$, we consider the minimization problem

$\min_{\boldsymbol{p}^{\prime}\in\mathcal{P}}\leavevmode\nobreak\
\mathcal{J}(\boldsymbol{p}^{\prime}):=\sum_{j=1}^{N_{\mathrm{O}}}(y_{j}-u_{\boldsymbol{p}^{\prime}}(\boldsymbol{x}_{j}))^{2}.$
(2)

Note that standard optimization algorithms for (2) require evaluations of the
objective $\mathcal{J}$ and its derivatives at numerous points
$\boldsymbol{p}^{\prime}$. Due to the definition of
$u_{\boldsymbol{p}^{\prime}}$, for each of these evaluations, one has to solve
the PDE (1) with parameter $\boldsymbol{p}^{\prime}$ and, in the case that
derivatives are required, an additional adjoint PDE. This is often
computationally prohibitive, in particular if the aim is to solve (2) with
limited hardware, e.g., directly on-board a robot. The main idea of our
approach is to resolve this problem by setting up a reusable neural network
model for the functional dependency $(\boldsymbol{x},\boldsymbol{p})\mapsto
u_{\boldsymbol{p}}(\boldsymbol{x})$ which is trained once offline, i.e., we
use a surrogate for $u_{\boldsymbol{p}^{\prime}}$ in (2).

###  II-C Neural Network (NN) Surrogate

####  II-C1 Network Architecture

As our surrogate model for the functional dependency
$(\boldsymbol{x},\boldsymbol{p})\mapsto u_{\boldsymbol{p}}(\boldsymbol{x})$,
we consider a Multilayer Perceptron (MLP), a fully-connected feed-forward
neural network. A popular architecture for an $L$-layer MLP
$\mathcal{N}\colon\mathbb{R}^{n_{\text{in}}}\to\mathbb{R}^{n_{\text{out}}}$
with input variable $\boldsymbol{z}$ is the following:

$\displaystyle\boldsymbol{l}_{0}:=\boldsymbol{z},\quad\boldsymbol{l}_{k}$
$\displaystyle:=\rho(\boldsymbol{W}_{k}\boldsymbol{l}_{k-1}+\boldsymbol{b}_{k}),\,\,k=1,...,L-1,$
(3) $\displaystyle\mathcal{N}(\boldsymbol{z})$
$\displaystyle:=\boldsymbol{W}_{L}\boldsymbol{l}_{L-1}+\boldsymbol{b}_{L}.$

Here, $\boldsymbol{W}_{k}\in\mathbb{R}^{n_{k}\times
n_{k-1}},\boldsymbol{b}_{k}\in\mathbb{R}^{n_{k}}$ are the $k$-th layer’s
weight matrix and bias vector, respectively, with
$n_{0},...,n_{L}\in\mathbb{N}$ the number of neurons at each layer
($n_{0}=n_{\text{in}}$, $n_{L}=n_{\text{out}}$), and
$\rho\colon\mathbb{R}\to\mathbb{R}$ is the nonlinear activation function,
evaluated componentwise when applied to a vector. We note that the definition
(3) follows common practice and has the last network layer as a scaling layer
without an activation function.

Figure 2: The surrogate model’s performance over its training period. Blue:
$H_{1}$-loss on training set. Orange: Mean Square Error (MSE) of function
values on test set (withheld from NN-training).

In our application, the network concatenates both the spatial coordinate
$\boldsymbol{x}\in\Omega$ and the source location
$\boldsymbol{p}\in\mathcal{P}$ as a single input
$\boldsymbol{z}:=[\boldsymbol{x};\boldsymbol{p}]$, and it outputs an
approximation of $u_{\boldsymbol{p}}(\boldsymbol{x})$. This means
$n_{\text{in}}=n_{0}=2+2=4,n_{\text{out}}=n_{L}=1$.

For the numerical experiments in this paper, we used a network with 5 layers,
30 neurons in each hidden layer ($n_{1},...,n_{4}=30$), and the softplus
activation function.

####  II-C2 Enforcing Boundary Conditions

We use the hard-enforcement approach suggested in [4] to realize the BCs. Our
approximator multiplies the MLP output $\mathcal{N}$ from (3) by a cutoff
function $\psi$ which vanishes on the spatial boundary $\partial\Omega$. This
ensures homogeneous Dirichlet BCs. Thus, we define

$\tilde{u}_{\boldsymbol{p}}(\boldsymbol{x}):=\psi(\boldsymbol{x})\mathcal{N}([\boldsymbol{x};\boldsymbol{p}]).$
(4)

As the function $\psi$, we use the so-called “$R$-$m$ approximate distance
function” to the boundary $\partial\Omega$ [4].

(a)

(b)

Figure 3: NN training results. Fig. (a) shows the difference between the NN-
surrogate $\tilde{u}_{\boldsymbol{p}}$ and the reference solution
$u_{\boldsymbol{p}}$ in $\Omega$, normalized by the maximum value of
$u_{\boldsymbol{p}}$ excluding a ball of radius $20\text{\,}\mathrm{m}$ around
the source, in the situation of Fig. 1. Fig. (b) shows the (Euclidean)
position error of the location $\boldsymbol{p}$ calculated by solving (2) with
our NN-surrogate approach with increasing amounts of White Gaussian Noise.

####  II-C3 Network Training

To train the NN-surrogate $\tilde{u}_{\boldsymbol{p}}$, we use a physics-
guided approach. This supervised learning method compares the network output
and its spatial gradients with a reference solution and attempts to make them
as similar as possible. We thus define the _Physics-Guided $H_{1}$-Loss_ by

$L_{\textbf{PG}}:=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\left[(\tilde{u}_{\boldsymbol{p}_{i}}-u_{\boldsymbol{p}_{i}})(\boldsymbol{q_{i}})\right]^{2}+\lVert(\nabla\tilde{u}_{\boldsymbol{p}_{i}}-\nabla
u_{\boldsymbol{p}_{i}})(\boldsymbol{q_{i}})\rVert_{2}^{2}.\vspace{-0.1cm}$

Here, $\tilde{u}_{\boldsymbol{p}}$ is the model in (4), $u_{\boldsymbol{p}}$
denotes a reference solution computed, e.g., by the Finite Element Method
(FEM), $\nabla$ is the spatial gradient,
$\\{(\boldsymbol{p}_{i},\boldsymbol{q}_{i})\\}\subset\mathcal{P}\times\Omega$
is a set of $N_{t}$ evaluation points, and $L_{\textbf{PG}}$ is understood as
a function of the weights and biases in $\mathcal{N}$. Note that the
$H_{1}$-loss, as opposed to the naive Mean Square Error (MSE), includes
information about the gradient of the reference solution.

In our numerical experiments, we used $N_{t}$ = 80000 points
$(\boldsymbol{p}_{i},\boldsymbol{q}_{i})$ in $L_{\textbf{PG}}$, sampled
randomly from $\mathcal{P}\times\Omega$. Since $u_{\boldsymbol{p}}$ has a
singularity at $\boldsymbol{p}$, we excluded points $\boldsymbol{q}_{i}$ with
a distance smaller than $0.02\text{\,}\mathrm{km}$ from the source location
$\boldsymbol{p}_{i}$. For the generation of the reference data, we used a
discretization based on the FEniCS library [5] on a Friedrichs–Keller
triangulation on a $241\times 241$ grid and linear Lagrangian elements. The
NN-surrogate was trained and implemented in Python using the Pytorch library
[6]. We used the ADAM Optimizer with learning rate $l=10^{-3}$ to train
$\tilde{u}_{\boldsymbol{p}}$ over 1500 epochs of the dataset. We additionally
computed the MSE on another set of 20000 test points randomly sampled from
$\mathcal{P}\times\Omega$, excluding singularities, to check the surrogate’s
generalization performance. The training results can be seen in Fig. 2.

###  II-D Optimization Algorithm

Once the NN-surrogate $\tilde{u}_{\boldsymbol{p}}$ has been trained, the
inverse problem (2), with $u_{\boldsymbol{p}}$ replaced by
$\tilde{u}_{\boldsymbol{p}}$, can be solved with standard optimization
algorithms. In our experiments, we used the LBFGS-B method [7] for this
purpose. The method had information about the gradients
$\nabla_{\boldsymbol{p}}\tilde{u}_{\boldsymbol{p}}$ via the automatic
differentiation functionalities of the Pytorch library [6]. Note that the
optimization step only requires gradients of the surrogate model
$\tilde{u}_{\boldsymbol{p}}$. Measurements of the concentration gradients do
not have to be taken in our GSL-approach and we do not require that the actual
concentration field is smooth. (It should, of course, be such that it is
sufficiently well described by (1).)

##  III Results

To test our method’s performance, we first checked that the surrogate model
provides a sensible approximation. In addition to the loss metric in the
training curve, we used a comparison with a reference solution for this
purpose. In Fig. 3(a), we see the difference between modeled values and ground
truth for a fixed source position, normalized by the maximum concentration
value of the reference function outside of its singular region. We can see
that the difference is within 2% of the maximum source value in all of
$\Omega$.

Next, we tested the performance of our PGNN in the GSL-problem (2). To do
this, we chose a set $\\{\boldsymbol{x}_{j}\\}$ of $N_{\mathrm{O}}=100$ evenly
spaced points on the circle of radius $0.25\text{\,}\mathrm{km}$ around
$($0.5\text{\,}\mathrm{km}$,$0.5\text{\,}\mathrm{km}$)$ as the observation
domain $D_{\mathrm{O}}$. Note that this choice of $D_{\mathrm{O}}$ is purely
for demonstration purposes. More complicated $D_{\mathrm{O}}$ could be
considered here without problems. We then sampled 324 source positions
$\boldsymbol{p}$ inside the set $\mathcal{P}$ on a grid. For each of these
source positions, first, we used the FEM to compute the ground truth values
$y_{j}$ at the observation points $\boldsymbol{x}_{j}$ without noise, and then
solved the inverse problem (2), localizing the gas source to get the estimate
$\hat{\boldsymbol{p}}$. The source positioning error is measured as
$\lVert\boldsymbol{p}-\hat{\boldsymbol{p}}\rVert_{2}$. We additionally
repeated these tests, adding white Gaussian noise of standard deviation
$\sigma\in\\{0.00625,0.0125,0.025\\}$ to the observations. Each inverse
problem was solved in the span of 1 second.

The results of this survey are depicted in Fig. 3(b) in a violin plot of the
error. We can see that, without noise, the network can identify the correct
source location up to an error of $6\text{\,}\mathrm{m}$ in the worst case and
up to an error of $4\text{\,}\mathrm{m}$ on average. Given that the resolution
of the FEM-approximation is $4.14\text{\,}\mathrm{m}$, this performance is
close to optimal. Even with additive noise corresponding to $\sim 6\%$ of the
highest values encountered on the circle, the solver is still able to give
source locations within $28\text{\,}\mathrm{m}$ of the real value. We also
notice that the median error (middle line) increases much slower than the
maximum error.

##  IV Conclusion

In this paper, we have demonstrated how PGNNs can be used to solve GSL-
problems in an accurate, computationally cheap, and robust manner. Our
approach is flexible and offers various possibilities for extensions, e.g., to
three dimensions and unknown wind/diffusion parameters. In future work, we
plan to include dynamic path optimization techniques and realistic sensor
models into our framework. Further open points are the validation of our
method in real world experiments and the detailed comparison with existing
GSL-techniques—two issues that are beyond the scope of the present paper.

## References

  * [1] Adam Francis, Shuai Li, Christian Griffiths and Johann Sienz  “Gas source localization and mapping with mobile robots: A review”  In _J. of Field Rob._ 39.8, 2022, pp. 1341–1373  DOI: https://doi.org/10.1002/rob.22109
  * [2] Tao Jing, Qing-Hao Meng and Hiroshi Ishida  “Recent Progress and Trend of Robot Odor Source Localization”  In _IEEJ Transactions on Electrical and Electronic Engineering_ 16.7, 2021, pp. 938–953  DOI: https://doi.org/10.1002/tee.23364
  * [3] Salah A. Faroughi et al.  “Physics-Guided, Physics-Informed, and Physics-Encoded Neural Networks and Operators in Scientific Computing: Fluid and Solid Mechanics”  In _J. Comput. Inf. Sci. Eng._ , 2024, pp. 1–45  DOI: 10.1115/1.4064449
  * [4] N. Sukumar and Ankit Srivastava  “Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks”  In _Comput. Methods Appl. Mech. Eng._ 389, 2022, pp. 114333  DOI: https://doi.org/10.1016/j.cma.2021.114333
  * [5] Martin S. Alnaes et al.  “The FEniCS Project Version 1.5”  In _Arch. Num. Softw._ 3, 2015  DOI: 10.11588/ans.2015.100.20553
  * [6] Adam Paszke et al.  “Pytorch: An imperative style, high-performance deep learning library”  In _Adv. Neural Inf. Process. Syst._ 32, 2019
  * [7] Richard H. Byrd, Peihuang Lu, Jorge Nocedal and Ciyou Zhu  “A limited memory algorithm for bound constrained optimization”  In _SIAM J. Sci. Comp._ 16.5 Society for IndustrialApplied Mathematics Publications, 1995, pp. 1190–1208  DOI: 10.1137/0916069

  *[GSL]: Gas Source Localization
  *[PGNN]: Physics-Guided Neural Network
  *[PDE]: Partial Differential Equation
  *[PDEs]: Partial Differential Equation
  *[BCs]: boundary condition
  *[NN]: Neural Network
  *[MLP]: Multilayer Perceptron
  *[FEM]: Finite Element Method
  *[MSE]: Mean Square Error