# Energy networks for state estimation with random sensors using sparse labels

Yash Kumar
Department of Mechanical Engineering
Delhi Technological University
Shahbad Daulatpur, Main Bawana Road, Delhi-110042, India
<EMAIL_ADDRESS>
&Souvik Chakraborty
Department of Applied Mechanics
Indian Institute of Technology Delhi
Hauz Khas - 110042, New Delhi, India
<EMAIL_ADDRESS>

###### Abstract

State estimation is required whenever we deal with high-dimensional dynamical
systems, as the complete measurement is often unavailable. It is key to
gaining insight, performing control or optimizing design tasks. Most deep
learning-based approaches require high-resolution labels and work with fixed
sensor locations, thus being restrictive in their scope. Also, doing Proper
orthogonal decomposition (POD) on sparse data is nontrivial. To tackle these
problems, we propose a technique with an implicit optimization layer and a
physics-based loss function that can learn from sparse labels. It works by
minimizing the energy of the neural network prediction, enabling it to work
with a varying number of sensors at different locations. Based on this
technique we present two models for discrete and continuous prediction in
space. We demonstrate the performance using two high-dimensional fluid
problems of Burgers’ equation and Flow Past Cylinder for discrete model and
using Allen–Cahn equation and Convection-diffusion equations for continuous
model. We show the models are also robust to noise in measurements.

_K_ eywords State estimation $\cdot$ Differentiable implicit layers $\cdot$
Dynamical systems

##  1 Introduction

State estimation is the ability to recover flow based on a few measurements.
It is an inverse problem and arises in many engineering applications such as
remote sensing, medical imaging, ocean dynamics, reservoir modeling, and blood
flow modeling. Uses of fluid estimation include flow control [1, 2], cardiac
blood flow modeling [3, 4, 5], ship wake identification [6], climate
prediction [7], optimizing machine design for low-drag vehicles, efficient
turbo-machines, etc. Few challenges faced in the processes are limited
sensors, sparse label data, moving sensors, ill-posed problems, noisy
measurements, etc. This work focuses on learning from moving sparse label data
and sensor measurements with a deep learning-based model using an implicit
optimization layer for network training. Sparse fluid is encountered in
various situations. One reason is that storing high-resolution data generated
during direct numerical simulation is challenging due to limited storage
space. It makes analysis, sharing, and visualization difficult. Another
important reason is that real data is hard to measure on a full scale, like
cardiovascular blood flow data obtained from flow magnetic resonance imaging
(MRI) [8, 9, 10].

For high dimensional state estimation problems like fluids, popular approaches
include library-based approaches observer dynamical system stochastic
approaches. Library-based methods use offline data, and the library consists
of generic modes such as Fourier, wavelet, discrete cosine transform basis, or
data specific Proper orthogonal decomposition (POD) or Dynamic mode
decomposition (DMD) modes, or training data. Library-based approaches using
sparse representation assume state can be expressed as the combination of
library elements. In an observer dynamical system, we assume the system’s
dynamics to produce a full state and update it based on new measurements to
reduce estimation error forming a closed feedback loop. The estimate is
maintained by Kalman filtering [11, 12, 13]. Tu et al. [14] applied dynamic
mode decomposition [15, 16] as a reduced-order model to Kalman smoother
estimate to identify coherent structures. Buffoni et al. [17] used a nonlinear
observer-based on Galerkin projection of Navier-Stokes equation to estimate
POD coefficients. Stochastic estimation was proposed by Adrian [18] for a
turbulence study where the conditional mean was approximated using a power
series. Extension to this method [19, 20, 21] can be found in the literature.
Bonnet et al. [22] extended stochastic approach to estimate POD coefficients.
A linear mapping between sensors and coefficients was assumed. These
approaches allow more flexibility in sensor placements and have been applied
for flow control over airfoil [23] and analyzing isotropic turbulence [24,
25].

We consider a problem with sparse label data whose position may vary with
time, for which POD can not be performed. Thus it restricts us from using
traditional POD-based approaches. Other deep learning-based approaches like
[26, 27, 28] give accurate predictions with fewer sensors but require high-
resolution labels for training. Gao et al. [29] uses physics-based loss for
super-resolution using spare data and thus assumes a fixed number of sensors
and their positions. In this work, we propose an Energy network for state
estimation with random sensors (ENSERS), a technique to learn a model from
spare training labels capable of predicting full states given a varied number
of sensors at random locations. We present two models trained using this
technique. The first one produces discrete high-dimensional predictions in
space. Second, produce continuous prediction utilizing the information of
coordinates. We demonstrate the results corresponding to four high complexity
problems: 2-dimensional (2D) coupled Burgers’ equation, transient flow,
Allen–Cahn equation, and Convection-diffusion equation.

The remainder of the paper is organized as follows. In Section 2, details on
the problem statement is provided. Details on the proposed approach are
provided in Section 3.1. Section 3 and 4 give details on discrete and
continuous formulation with two numerical examples in each to illustrate the
performance of the proposed approach. Finally, Section 5 provides the
concluding remarks.

##  2 Problem statement

Consider a dynamical system obtained by partial discretization of the
$d$-dimensional governing differential equations:

$J_{t}(x,t)=F(x,J(x,t)),\quad J^{n}=J(x,t_{n}),\quad x\in\Omega$ (1)

The simulation time domain is discretized by $L$ steps and the space domain is
discretized by $\omega$ segments resulting in

$Z=\\{J^{l}_{m}\in\mathbb{R}^{\omega}\ |\ l=0,...,L-1,\ m=0,...,M-1\\}$ (2)

where $Z\in\mathbb{R}^{L\times M\times\omega}$, $l$ is time step index, $M=$
number of system’s state variables. e.g. $m=0$ represents x-velocity and $m=1$
represents y-velocity in $\S$ 3.2 of 2D coupled Burgers’ equation. We consider
sensor location and data location, represented by integers in discrete domain
respectively as

$S=\\{\lambda^{l}_{m}\in\mathbb{Z}^{p}\cap[0,\omega-1]^{p}\ |\ l=0,...,L-1,\
m=0,...,M-1\\}$ (3) $T=\\{\pi^{l}_{m}\in\mathbb{Z}^{h}\cap[0,\omega-1]^{h}\ |\
l=0,...,L-1,\ m=0,...,M-1\\}$ (4)

where $p=$ number of sensors, $h=$ number of data nodes,
$S\in\mathbb{Z}^{L\times M\times p}\cap[0,\omega-1]^{L\times M\times p}$,
$T\in\mathbb{Z}^{L\times M\times h}\cap[0,\omega-1]^{L\times M\times h}$. The
corresponding sensor values and data values are

$X=\\{\sigma^{l}_{m}\in\mathbb{R}^{p}\ |\ l=0,...,L-1,\ m=0,...,M-1\\}$ (5)
$\Phi=\\{\psi^{l}_{m}\in\mathbb{R}^{h}\ |\ l=0,...,L-1,\ m=0,...,M-1\\}$ (6)

$\sigma^{l}_{m}=\Lambda^{l}_{m}J^{l}_{m}$ (7)
$\psi^{l}_{m}=\Pi^{l}_{m}J^{l}_{m}$ (8)

where $X\in\mathbb{R}^{L\times M\times p}$, $\Phi\in\mathbb{R}^{L\times
M\times h}$. $\Lambda^{l}_{m}\in\mathbb{R}^{p\times\omega}$ and
$\Pi^{l}_{m}\in\mathbb{R}^{h\times\omega}$ are measurement matrices composed
of one-hot row vectors. $\Lambda^{l}_{m}$ and $\Pi^{l}_{m}$ are defined as

$\Lambda^{l}_{mij}=\left\\{\begin{array}[]{ c l }1&\quad\textrm{if
}\lambda^{l}_{mi}=j\\\ 0&\quad\textrm{otherwise}\end{array}\right.$ (9)
$\Pi^{l}_{mij}=\left\\{\begin{array}[]{ c l }1&\quad\textrm{if
}\pi^{l}_{mi}=j\\\ 0&\quad\textrm{otherwise}\end{array}\right.$ (10)

Note that integer senor locations $\lambda^{l}_{m}$ are selected randomly and
are kept fixed during training. Similarly, a different set of sensor locations
$S$ is selected for testing the network. In this work, we aim to use sensor
data from a set of $\gamma$ system states and produce $\gamma$ high
dimensional states. Thus we divided the sensor values and data into chunks,
each with $\gamma$ states.

$\chi^{k}=\\{\sigma^{k*z+i}_{m}\in\mathbb{R}^{p}\ |\ i=0,...,\gamma-1,\
m=0,...,M-1\\}$ (11)

$\phi^{k}=\\{\psi^{k*z+i}_{m}\in\mathbb{R}^{h}\ |\ i=0,...,\gamma-1,\
m=0,...,M-1\\}$ (12)

where $\chi^{k}\in\mathbb{R}^{\gamma\times M\times p}$ are sensor values for
$\gamma$ time steps, $\phi^{k}\in\mathbb{R}^{\gamma\times M\times h}$ are data
values for $\gamma$ time steps, $k=0,...,N-1$, $N$ is number of train samples,
$z$ is time steps between first state $\sigma^{k*z}_{m}$ of each tensor
$\chi^{k}$.

##  3 Discrete space models

###  3.1 Proposed approach

Figure 1: Network architecture of proposed Ensers model.

In this section, we propose a novel deep learning-based framework for state
estimation. Prediction of $\gamma$ high dimensional states is done via feed-
forward neural network (FNN) $\Gamma$ using optimized reduced state vector
$\xi^{k}\in\mathbb{R}^{\varsigma}$:

$\Gamma(\xi^{k})=D^{k}=\\{D^{ki}_{m}\in\mathbb{R}^{\omega}\ |\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$ (13)

where $D^{k}\in\mathbb{R}^{\gamma\times M\times\omega}$ is a third-order
tensor of predicted states. The output vector of FNN has a dimension of
$\gamma*M*\omega$ which is reshaped into a tensor of shape $\gamma\times
M\times\omega$. The reduced state is obtained by solving the following
minimization problem using sensor data $\chi^{k}$ and predicted states by the
neural network.

$\xi^{k}=\arg\min_{\tilde{\xi^{k}}}\|vec(\chi^{k}-\rho(\tilde{\xi^{k}}))\|_{2}^{2}$
(14a)

$\rho(\tilde{\xi^{k}})=\\{\tilde{Q}^{ki}_{m}\in\mathbb{R}^{p}\ |\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$ (14b)

$\tilde{Q}^{ki}_{m}=\Lambda^{kz+i}_{m}\tilde{D}^{ki}_{m}$ (14c)

$\Gamma(\tilde{\xi^{k}})=\\{\tilde{D}^{ki}_{m}\in\mathbb{R}^{\omega}\ |\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$ (14d)

where $\tilde{Q}^{ki}_{m}\in\mathbb{R}^{p}$ are values of predicted states at
sensor locations $\lambda^{kz+i}_{m}$,
$\Lambda^{kz+i}_{m}\in\mathbb{R}^{p\times\omega}$ is measurement matrix. Note
that in practice $\xi^{k}$ is obtained by a few steps of gradient descent
instead of global minimization. This inner optimization loop is implemented
using the library ‘higher’ [30] in PyTorch. Also, dimension $\varsigma$ of the
reduced state vector is quite small (e.g. 8 in first experiment $\S$ 3.2)
therefore the time required for gradient descent steps is negligible. The
network $\Gamma$ is trained by minimizing data loss and physics-based loss
$P(\chi^{k})$ across training samples $N$.

$\theta^{*}=\arg\min_{\theta}\sum_{k=0}^{N-1}\|vec(\phi^{k}-\Upsilon(\chi^{k}))\|_{2}^{2}+P(\chi^{k})$
(15a)

$\Upsilon(\chi^{k})=\\{Y^{ki}_{m}\in\mathbb{R}^{h}\ |\ i=0,...,\gamma-1,\
m=0,...,M-1\\}$ (15b)

$Y^{ki}_{m}=\Pi^{kz+i}_{m}D^{ki}_{m}$ (15c)

where $D^{ki}_{m}\in\mathbb{R}^{\omega}$ are predicted states by network
$\Gamma$ from Eq. (13), $Y^{ki}_{m}\in\mathbb{R}^{h}$ are values of predicted
states at data locations $\pi^{kz+i}_{m}$, $\Upsilon(\chi^{k})$ are predicted
states composed of $\gamma$ time steps. Fig. 1 shows the network architecture
during training. Equations (13) and (14) together forms the implicit
optimization layer shown in the Fig. 1. Training and testing procedure is
shown in Algorithm 1 and 2 respectively. Note that for demonstration purpose,
in all Algorithms a batch size of $1$ is considered but in practice batch size
is selected based on problem as mentioned in each experiment section. Also we
use Huber loss function during testing as it is more robust to noise than mean
squared error (MSE) loss function.

Methods based on implicit optimization layers come under the category of
Optimization-based Modeling architectures [31] and are well-studied for
generic classification, and structured prediction tasks [32, 33, 34, 35]. The
most common way of training such models is through Unrolled Differentiation
[36, 37, 38, 33]. It is done by introducing an optimization procedure such as
gradient descent into the inference procedure. Other ways of training include
Implicit argmin differentiation using implicit function theorem but need
argmin operations to be convex. This method can be found in works of [39, 40].
In this work, we use Unrolled Differentiation technique for training.

1

2

31:Inputs: $Z,S,T$. {Eq. (2), Eq. (3), Eq. (4)}

42:Set Hyper-parameters: $\eta_{o}$: outer learning rate, $\eta_{i0}$: inner
learning rate at epoch=0, $\hat{\eta_{i}}$: inner learning rate rate, $\beta$:
batch size, $I_{o}$: outer iterations, $I_{i}$: inner iterations,
$\hat{\zeta}$: physics penalty rate, $\zeta_{0}$: physics penalty at epoch=0,
$N$: number of train samples, $z$, $\gamma$.

53:Calculate data-set: $\chi$, $\phi$, $\Lambda$, $\Pi$

64:Initialize: Neural network model: $\Gamma(\cdot;\theta)$

75:for $\iota_{o}=0$ to $I_{o}-1$ do {Outer optimization Loop}

6: $\eta_{i}=\eta_{i0}+\iota_{o}\hat{\eta_{i}}$ {Schedule inner learning rate}

7: $\zeta=\zeta_{0}+\iota_{o}\hat{\zeta}$ {Schedule physics penalty}

88: for $k=0$ to $N-1$ do

9: Initialize: $\tilde{\xi^{k}}$

910: for $\iota_{i}=0$ to $I_{i}-1$ do {Inner optimization Loop}

11: $\tilde{D^{k}}=\Gamma(\tilde{\xi^{k}})$ {Eq. (14d)}

12: $\tilde{Q^{k}}\leftarrow$ DrawValuesAtSensorLocations($\tilde{D^{k}}$)
{Eq. (14c)}

13: $\mathfrak{L}=MSE(\chi^{k},\tilde{Q^{k}})$

14: $\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}\leftarrow$
Backprop($\mathfrak{L}$)

15:
$\tilde{\xi^{k}}=\tilde{\xi^{k}}-\eta_{i}\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}$

16: end for

17: $\xi^{k}=\tilde{\xi^{k}}$

18: $D^{k}=\Gamma(\xi^{k})$ {Eq. (13)}

1019: $Y^{k}\leftarrow$ DrawValuesAtDataLocations($D^{k}$) {Eq. (15c)}

20: $\mathcal{L}=MSE(\phi^{k},Y^{k})+\zeta P(D^{k})${Calculate loss}

21: $\frac{\partial\mathcal{L}}{\partial\theta}\leftarrow$
Backprop($\mathcal{L}$)

22:
$\theta\leftarrow\theta-\eta_{o}\frac{\partial\mathcal{L}}{\partial\theta}$
{Update weights}

23: end for

1124:end for

1225:Output: Trained network $\Gamma(\cdot;\theta^{*})$.

Algorithm 1 Training ENSERS

1

21:Inputs: Trained network $\Gamma(\cdot;\theta^{*})$, $S$. {Eq. (3)}

32:Set Hyper-parameters: $\eta_{i}$: inner learning rate, $I_{i}$: inner
iterations, $\hat{N}$: number of test samples, $z$, $\gamma$.

43:Calculate data-set: $\chi$, $\Lambda$

54:for $k=0$ to $\hat{N}-1$ do

5: Initialize: $\tilde{\xi^{k}}$

66: for $\iota_{i}=0$ to $I_{i}-1$ do {Inner optimization Loop}

7: $\tilde{D^{k}}=\Gamma(\tilde{\xi^{k}})$ {Eq. (14d)}

8: $\tilde{Q^{k}}\leftarrow$ DrawValuesAtSensorLocations($\tilde{D^{k}}$) {Eq.
(14c)}

9: $\mathfrak{L}=MSE(\chi^{k},\tilde{Q^{k}})$

10: $\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}\leftarrow$
Backprop($\mathfrak{L}$)

11:
$\tilde{\xi^{k}}=\tilde{\xi^{k}}-\eta_{i}\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}$

712: end for

13: $\xi^{k}=\tilde{\xi^{k}}$

14: $D^{k}=\Gamma(\xi^{k})$ {Eq. (13)}

815:end for

916:Output: Predicted states $\\{D^{ki}_{m}\ |\ k=0,...,N-1,\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$.

Algorithm 2 Testing ENSERS

####  3.1.1 Physics-based loss function

The physics-based loss function is used in the approach because of spare
training labels. Training any network with just spare labels will produce
garbage values on nodes without labels. Physics-based loss functions have been
used to train neural networks for solving PDEs. A popular class of methods is
PINNs [41]. The basic idea here is to place a neural network prior to the
state variable and then estimate the neural network parameters by using a
physics-informed loss function. Several improvements to the originally
proposed PINN can also be found in the literature. For example, Zhu et al.
[42] developed convolutional PINN for time-independent systems. Geneva and
Zabaras [43] used physics constrained auto-regressive model for surrogate
modeling of dynamical systems. We use Runge-Kutta methods with $q$ stages for
defining loss between $\gamma$ state predictions. Let

$V^{n}=D^{k0}$ (16a)

$V^{n+c_{i}}=D^{ki},\quad i=1,...,q$ (16b)

$V^{n+1}=D^{k(\gamma-1)}$ (16c)

where $D^{ki}$ are network prediction. General form of Runge-Kutta methods
with q stages applied to Eq. (1):

$V^{n+c_{i}}=V^{n}-\Delta t\sum_{j=1}^{q}a_{ij}F(V^{n+c_{j}}),\quad i=1,..,q$
(17a)

$V^{n+1}=V^{n}-\Delta t\sum_{j=1}^{q}b_{j}F(V^{n+c_{j}})$ (17b)

We use the Implicit Runge-Kutta methods with q stages and thus parameters
$\\{a_{ij},b_{j},c_{j}\\}$ are chosen accordingly. Now, shifting second term
on right hand side (RHS) in Eq. (17) to left hand side (LHS) and replacing
exact operator $F$ with numerical gradient based operator $\hat{F}$

$\hat{W_{i}}=V^{n+c_{i}}+\Delta
t\sum_{j=1}^{q}a_{ij}\hat{F}(V^{n+c_{j}}),\quad i=1,..,q$ (18a)

$\hat{W_{q+1}}=V^{n+1}+\Delta t\sum_{j=1}^{q}b_{j}\hat{F}(V^{n+c_{j}})$ (18b)

$P(\chi^{k})=\sum_{i=1}^{q+1}\|vec(\hat{W_{i}}-V^{n})\|_{2}^{2}$ (18c)

where $\hat{W_{i}}$ are different estimates of $V^{n}$, $P(\chi^{k})$ is the
physics based loss function. For calculating loss the spatial gradients are
approximated using Sobel filter 2D convolutions [44]. See $\S$ A for
additional details. Note that physics-based loss function is only used during
training of network.

In next section, we present two examples to show model proposed is able to
learn from sparse moving data labels. We illustrate the performance of the
proposed approach with plots of prediction. We show that model is robust
against noisy sensor measurements by showing error corresponding to various
noise level. Error used as a quantitative metric in plots in defined as

$\epsilon^{k}=\frac{\left\|D^{k\gamma^{*}}_{m}-J^{k*z+\gamma^{*}}_{m}\right\|_{2}}{\left\|J^{k*z+\gamma^{*}}_{m}\right\|_{2}}\quad
k=0,...,\hat{N}$ (19)

where $\epsilon^{k}\in\mathbb{R}$ represents the error,
$\hat{N}numberoftestsamples$, $J^{k*z+\gamma^{*}}_{m}\in\mathbb{R}^{\omega}$
are the true state and $D^{k\gamma^{*}}_{m}\in\mathbb{R}^{\omega}$ are the
predicted state using the proposed approach. $\left\|\cdot\right\|_{2}$
represents the L2 norm.

###  3.2 Experiment: 2D coupled Burgers’ equation

As the first example, we consider the the 2D coupled Burgers’ system. It has
the same convective and diffusion form as the in-compressible Navier-Stokes
equations. It is an important model for understanding of various physical
flows and problems, such as hydrodynamic turbulence, shock wave theory, wave
processes in thermo-elastic medium, vorticity transport, dispersion in porous
medium. The governing equations for Burgers’ equation takes the following
form:

$\bm{u}_{t}+\bm{u}\cdot\nabla\bm{u}-\nu\Delta\bm{u}=0,$ (20)

with periodic boundary condition

$\begin{split}\bm{u}\left(x=0,y,t\right)&=\bm{u}\left(x=L,y,t\right),\\\
\bm{u}\left(x,y=0,t\right)&=\bm{u}\left(x,y=L,t\right).\end{split}$ (21)

Eq. (20) can be written in expanded form as

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+v\frac{\partial u}{\partial y}-\nu(\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=0$ (22)
$\displaystyle\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial
x}+v\frac{\partial v}{\partial y}-\nu(\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial y^{2}})=0,$

where $\nu$ is viscosity, $u$ and $v$ are the $x$ and $y$ components of
velocity. We consider $\\{x,y\\}\in[0,1]$. The initial condition is defined
using truncated Fourier series with random coefficients:

$\bm{u}(x,y,t=0)=\frac{2\bm{w}(x,y)}{\max_{\\{x,y\\}}|\bm{w}(x,y)|}+\bm{c},$
(23)

where

$\bm{w}(x,y)=\sum_{i=-L}^{L}\sum_{j=-L}^{L}\bm{a}_{ij}\sin(2\pi(ix+jy))+\bm{b}_{ij}\cos(2\pi(ix+jy)),$
(24)

where $\bm{a}_{ij},\bm{b}_{ij}\sim\bm{N}(0,\mathbf{I}_{2})$, $L=4$ and
$\bm{c}\sim\bm{\mathcal{U}}(-1,1)\in\mathbb{R}^{2}$.

####  3.2.1 Data-set and Model Parameters

We use FeNICS [45] computing platform to solve the partial differential
equations (22) for generating data-set. We discretize the spatial domain with
$64\times 64$ grid and use a time-step of $0.005$. Parameters related to data-
set considered are displayed in Table 1. We use noisy measurements to test the
approach. Noise level is measured in signal to noise ratio (SNRdB) in
decibels(dB) and is represented by $SNR_{dB}$. Signal after adding noise $r$
is formulated as:

$r=s+(\frac{P}{2*10^{\frac{SNR_{dB}}{10}}})^{0.5}*\mathcal{N}$ (25a)
$P=\frac{\sum s^{2}}{n}$ (25b)

where, $s$ is noise-free signal, $\mathcal{N}$ is random variable with
standard normal distribution. Plots depicting different noise levels used for
testing is shown in Fig. 2.

Figure 2: Plots of 2D Burgers velocity for visualizing different noise levels. SNRdB represents signal to noise ratio in decibels. Table 1: Data-set parameters for 2D Burgers problem. $\gamma$ |  $M$ |  $L$ |  $z$ |  $\omega$ |  domain
---|---|---|---|---|---
5 |  2 |  50 |  2 |  3969 |  $63\times 63$
Table 2: Network architecture of proposed Ensers for 2D coupled Burgers’ problem. Layer |  Input |  Output |  Activation
---|---|---|---
FC |  8 |  64 |  Softplus
FC |  64 |  64 |  Softplus
FC |  64 |  $\omega*M*\gamma$ |  linear
Table 3: Training hyper parameters of Ensers for 2D Burgers problem. $\eta_{o}$ |  $\eta_{i0}$ |  $\hat{\eta_{i}}$ |  $\beta$ |  $I_{o}$ |  $I_{i}$ |  $\zeta_{0}$ |  $\hat{\zeta}$ |  $p$ |  $h$ |  $N$ |  $\varsigma$
---|---|---|---|---|---|---|---|---|---|---|---
0.0002 |  0.1 |  0.006 |  11 |  2001 |  4 |  0.005 |  0.0001 |  32 |  800 |  22 |  8
Table 4: Testing hyper-parameters of Ensers for 2D Burgers problem. $\eta_{i}$: inner learning rate, $I_{i}$: inner iterations $\eta_{i}$ |  $I_{i}$ |  $p$ |  $\hat{N}$
---|---|---|---
5 |  100 |  [4, 16] |  12

FNN used in the model is a shallow network with three fully connected (FC)
layers. We use Softplus [46] activation function which has smooth first
derivatives. Thus it helps avoid discontinuities because unrolling the
inference procedure involves computing
$\nabla_{\theta}\nabla_{\xi}(\mathfrak{L})$. Network architecture is displayed
in Table 2. We use 800 data nodes for training from each system variable data
which is $20\%$ of high-resolution data. Training hyper-parameters and
training data are summarized in table 3. For testing, we use sensor locations
different from the ones used during training.

####  3.2.2 Results and Discussions

Figure 3: Ensers predictions of a 2D coupled Burgers’ test case. (Top to
bottom) x-velocity FEM target solution, x-velocity Ensers prediction,
x-velocity L1error, y-velocity FEM target solution, y-velocity Ensers
prediction and y-velocity L1error.

Figure 4: Violin plot representing error $\epsilon$ distribution in prediction
vector of x and y velocity for different noise levels in sensor measurements
for 2D coupled Burgers’ test case.

Fig. 3 shows the prediction of Ensers for x and y velocity components at
various times in simulation with 16 sensors. We see that it produces results
that accurately captures the current state of the system at most points in the
domain. If a network is trained without a physics-based loss function then it
will produce garbage values at points where data is absent. Fig 4 shows violin
plot of error vector $\epsilon$ defined in Eq. (19) between target and
prediction with various noise level in sensor measurement for test case. In
Eq. (19), $m$ for x velocity and y velocity are $0$ and $1$ respectively and
$\gamma^{*}=2$. The plot represents the distribution of error $\epsilon$ in
the prediction vector. Greater spread corresponding to a point on the y-axis
corresponds to more values present in the error vector $\epsilon$ around that
value of the point on the y-axis. We see the model is robust to noise as the
mean error in the plot are close for cases of no noise($SNR_{dB}=None$) and
high noise($SNR_{dB}=10$).

###  3.3 Experiment: Flow Past Cylinder

As the second example, we consider the flow past a cylinder problem. It is a
well known canonical problem and is characterized by periodic laminar flow
vortex shedding. System is governed by incompressible, laminar, Newtonian
fluid equations:

$\frac{\partial(u)}{\partial x}+\frac{\partial(v)}{\partial y}=0$ (26a)

$\frac{\partial(u)}{\partial t}+u\frac{\partial(u)}{\partial
x}+v\frac{\partial(u)}{\partial y}=-\frac{1}{\rho}\frac{\partial(p)}{\partial
x}+\nu\left(\frac{\partial^{2}(u)}{\partial
x^{2}}+\frac{\partial^{2}(u)}{\partial y^{2}}\right)$ (26b)

$\frac{\partial(v)}{\partial t}+u\frac{\partial(v)}{\partial
x}+v\frac{\partial(v)}{\partial y}=-\frac{1}{\rho}\frac{\partial(p)}{\partial
y}+\nu\left(\frac{\partial^{2}(v)}{\partial
x^{2}}+\frac{\partial^{2}(v)}{\partial y^{2}}\right)$ (26c)

####  3.3.1 Data-set and Model Parameters

A schematic representation of the computational domain is shown in Fig. 5(a).
The circular cylinder is considered to have a diameter of $1$ unit. The center
of the cylinder is located at a distance of $8$ units from the inlet. The
outlet is located at a distance of $25$ units from the center of the cylinder.
The sidewalls are at $4$ units distance from the center of the cylinder. At
the inlet boundary, a uniform velocity of $1$ unit along the $X$-direction is
applied. Pressure boundary condition with $P=0$ is considered at the outlet. A
no-slip boundary at the cylinder surface is considered. Coordinate of the
snapshot cutout stretches from $[1.5,-2]\times[5.5,2]$ which is discretized
into $64\times 64$ points in $x$ and $y$ directions (see Fig. 5(b)).

Figure 5: (a) Schematic representation of the computational domain with boundary conditions at the inlet and the outlet. The cylinder has a diameter of 1 unit. A no-slip boundary is considered at the cylinder wall. Zero pressure gradient at the inlet and zero velocity gradient at the outlet are considered. (b) Schematic of the problem domain with snapshot cutout of $a\times a$. For flow past cylinder problems, $a=4$ units. The schematics are not to scale. Figure 6: Plots of velocity of Flow Past Cylinder for visualizing different noise levels. SNRdB represents signal to noise ratio in decibels. Table 5: Data-set parameters for Flow Past Cylinder problem. $\gamma$ |  $M$ |  $L$ |  $z$ |  $\omega$ |  $\nu$ |  Re |  domain
---|---|---|---|---|---|---|---
5 |  3 |  25 |  1 |  4096 |  0.005 |  200 |  $64\times 64$

The data-set is generated by using Unsteady Reynolds-averaged Navier Stokes
(URANS) simulation in OpenFoam [47]. The overall problem domain is discretized
into 63420 elements with finer mesh near the cylinder. Time step $\delta
t=0.02$ units is considered. Code for OpenFOAM simulation can be found at
[48]. Parameters related to the data-set considered are displayed in Table 5.

We use a shallow network in the model with two fully connected(FC) layers.
Network architecture is displayed in Table 6. We use 500 data nodes for
training from each system variable data which is $12.2\%$ of high-resolution
data. Training and testing hyper-parameters are shown in table 7 and 8
respectively. Sensor locations for testing are different from the ones used
during training. For this problem physics-based loss function is extended to
include the continuity equation:

$P(\chi^{k})=\sum_{i=1}^{q+1}\|vec(W_{i}-V^{n})\|_{2}^{2}+\sum_{i=0}^{\gamma-1}\|vec(\frac{\partial(V^{n}_{0})}{\partial
x}+\frac{\partial(V^{n}_{1})}{\partial y})\|_{2}^{2}$ (27)

Table 6: Network architecture of proposed Ensers for Flow Past Cylinder problem. Layer |  Input |  Output |  Activation
---|---|---|---
FC |  8 |  64 |  Softplus
FC |  64 |  $\omega*M*\gamma$ |  linear
Table 7: Training hyper parameters of Ensers for Flow Past Cylinder. $\eta_{o}$ |  $\eta_{i0}$ |  $\hat{\eta_{i}}$ |  $\beta$ |  $I_{o}$ |  $I_{i}$ |  $\zeta_{0}$ |  $\hat{\zeta}$ |  $p$ |  $h$ |  $N$ |  $\varsigma$
---|---|---|---|---|---|---|---|---|---|---|---
0.0003 |  0.1 |  0.002 |  9 |  3001 |  5 |  0.01 |  0.0006 |  16 |  500 |  18 |  8
Table 8: Testing hyper parameters of Ensers for Flow Past Cylinder problem. $\eta_{i}$: inner learning rate, $I_{i}$: inner iterations $\eta_{i}$ |  $I_{i}$ |  $p$ |  $\hat{N}$
---|---|---|---
5 |  100 |  [4, 16] |  12

####  3.3.2 Results and Discussions

Figure 7: Ensers predictions of a Flow Past Cylinder test case. (Top to
bottom) x-velocity target solution, x-velocity Ensers prediction, x-velocity
L1error, y-velocity target solution, y-velocity Ensers prediction, y-velocity
L1error, pressure target solution, pressure Ensers prediction and pressure
L1error.

Figure 8: Violin plot representing error $\epsilon$ distribution in prediction
vector of x,y velocity and pressure for different noise levels in sensor
measurements for Flow Past Cylinder test case.

Fig. 7 shows prediction of Ensers for pressure, x and y velocity components
with $16$ sensors. We see that it produces results that accurately capture the
current state of the system at most points in the domain. Fig. 8 shows violin
plot of error defined in Eq. (19) between target and prediction with various
noise levels in sensor measurement for the test case. In Eq. (19), $m$ for x
velocity, y velocity and pressure are $0$, $1$ and $2$ respectively and
$\gamma^{*}=2$. The plot represents the distribution of error in prediction
vectors. We see model is robust to noise as mean error in Fig 8 are close for
cases of no noise($SNR_{dB}=None$) and high noise($SNR_{dB}=10$).

##  4 Continuous space models

###  4.1 Proposed approach

In this section, we propose a novel deep learning-based continuous framework
for state estimation. Prediction of $\gamma$ high dimensional states is done
via multiple passes from a feed-forward neural network (FNN) $\Gamma$ for
every collocation point with coordinate vector $X_{r}$ using optimized reduced
state vectors $\xi^{k}\in\mathbb{R}^{\varsigma}$:

$\begin{split}\Gamma(\Xi^{k}_{r})=\\{D^{ki}_{mr}\in\mathbb{R}\ |\
i=0,...,\gamma-1,\ m=0,...,M-1\\},\quad\ r=0,...,\omega\end{split}$ (28)

$\Xi^{k}_{r}=\\{\xi^{k},X_{r}\\},\quad\ r=0,...,\omega$ (29)

where output vector of FNN $\Gamma$ has dimension $\gamma*M$, $D^{ki}_{mr}$ is
predicted states value at $r^{th}$ collocation point, $X_{r}$ are coordinates
of $r^{th}$ collocation point. Output of $\omega$ FNNs are combined to form
third-order tensor $D^{k}\in\mathbb{R}^{\gamma\times M\times\omega}$. Reduced
state is obtained by solving the following minimization problem using sensor
data $\chi^{k}$ and predicted states by the neural network.

$\xi^{k}=\arg\min_{\tilde{\xi^{k}}}\|vec(\chi^{k}-\rho(\tilde{\Xi^{k}}))\|_{2}^{2}$
(30a)

$\rho(\tilde{\Xi^{k}})=\\{\tilde{Q}^{ki}_{m}\in\mathbb{R}^{p}\ |\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$ (30b)

$\tilde{Q}^{ki}_{m}=\Lambda^{kz+i}_{m}\tilde{D}^{ki}_{m}$ (30c)

$\begin{split}\Gamma(\tilde{\Xi^{k}_{r}})=\\{\tilde{D}^{ki}_{mr}\in\mathbb{R}\
|\ i=0,...,\gamma-1,\ m=0,...,M-1\\},\\\ \quad r=0,...,\omega\end{split}$
(30d)

$\tilde{\Xi^{k}_{r}}=\\{\tilde{\xi^{k}},X_{r}\\},\quad r=0,...,\omega$ (30e)

where $\tilde{Q}^{ki}_{m}\in\mathbb{R}^{p}$ are values of predicted states at
sensor locations $\lambda^{kz+i}_{m}$. Note that in practice $\xi^{k}$ is
obtained by a few steps of gradient descent instead of global minimization.
The network is trained by minimizing data loss and physics-based loss
$P(\chi^{k})$ across training samples $N$.

$\theta^{*}=\arg\min_{\theta}\sum_{k=0}^{N-1}\|vec(\phi^{k}-\Upsilon(\chi^{k}))\|_{2}^{2}+P(\chi^{k})$
(31a)

$\Upsilon(\chi^{k})=\\{Y^{ki}_{m}\in\mathbb{R}^{h}\ |\ i=0,...,\gamma-1,\
m=0,...,M-1\\}$ (31b)

$Y^{ki}_{m}=\Pi^{kz+i}_{m}D^{ki}_{m}$ (31c)

where $D^{ki}_{m}\in\mathbb{R}^{\omega}$ are predicted states from Eq. (28),
$Y^{ki}_{m}\in\mathbb{R}^{h}$ are values of predicted states at data locations
$\pi^{kz+i}_{m}$, $\Upsilon(\chi^{k})\in\mathbb{R}^{\gamma\times M\times p}$
are predicted states composed of $\gamma$ time steps, $\theta$ are network
parameters. Fig. 1 shows the network architecture during training. Equations
(28) and (30) together forms the implicit optimization layer shown in the Fig.
1 for the continuous model. Training and testing procedure is shown in
Algorithm 3 and 4 respectively.

1

2

31:Inputs: $Z,S,T$. {Eq. (2), Eq. (3), Eq. (4)}

42:Set Hyper-parameters: $\eta_{o}$: outer learning rate, $\eta_{i0}$: inner
learning rate at epoch=0, $\hat{\eta_{i}}$: inner learning rate rate, $\beta$:
batch size, $I_{o}$: outer iterations, $I_{i}$: inner iterations,
$\hat{\zeta}$: physics penalty rate, $\zeta_{0}$: physics penalty at epoch=0,
$N$: number of train samples, $z$, $\gamma$.

53:Calculate data-set: $\chi$, $\phi$, $\Lambda$, $\Pi$

64:Initialize: Neural network model: $\Gamma(\cdot;\theta)$

75:for $\iota_{o}=0$ to $I_{o}-1$ do {Outer optimization Loop}

6: $\eta_{i}=\eta_{i0}+\iota_{o}\hat{\eta_{i}}$ {Schedule inner learning rate}

7: $\zeta=\zeta_{0}+\iota_{o}\hat{\zeta}$ {Schedule physics penalty}

88: for $k=0$ to $N-1$ do

9: Initialize: $\tilde{\xi^{k}}$

910: for $\iota_{i}=0$ to $I_{i}-1$ do {Inner optimization Loop}

1011: $\tilde{\Xi^{k}_{r}}=\\{\tilde{\xi^{k}},X_{r}\\},\quad r=0,...,\omega$

12: $\tilde{D^{k}_{r}}=\Gamma(\tilde{\Xi^{k}_{r}}),\quad r=0,...,\omega$ {Eq.
(30d)}

13: $\tilde{Q^{k}}\leftarrow$ DrawValuesAtSensorLocations($\tilde{D^{k}}$)
{Eq. (30c)}

14: $\mathfrak{L}=MSE(\chi^{k},\tilde{Q^{k}})$

15: $\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}\leftarrow$
Backprop($\mathfrak{L}$)

16:
$\tilde{\xi^{k}}=\tilde{\xi^{k}}-\eta_{i}\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}$

1117: end for

18: $\hat{X_{r}}=X_{r}$

19: $\Xi^{k}_{r}=\\{\tilde{\xi^{k}},\hat{X_{r}}\\},\quad r=0,...,\omega$

20: $D^{k}_{r}=\Gamma(\Xi^{k}_{r}),\quad r=0,...,\omega$ {Eq. (28)}

21: $Y^{k}\leftarrow$ DrawValuesAtDataLocations($D^{k}$) {Eq. (31c)}

22: $\mathcal{L}=MSE(\phi^{k},Y^{k})+\zeta P(D^{k})${Calculate loss}

23: $\frac{\partial\mathcal{L}}{\partial\theta}\leftarrow$
Backprop($\mathcal{L}$)

24:
$\theta\leftarrow\theta-\eta_{o}\frac{\partial\mathcal{L}}{\partial\theta}$
{Update weights}

25: end for

1226:end for

1327:Output: Trained network $\Gamma(\cdot;\theta^{*})$.

Algorithm 3 Training ENSERS

1

21:Inputs: Trained network $\Gamma(\cdot;\theta^{*})$, $S$. {Eq. (3)}

32:Set Hyper-parameters: $\eta_{i}$: inner learning rate, $I_{i}$: inner
iterations, $\hat{N}$: number of test samples, $z$, $\gamma$.

43:Calculate data-set: $\chi$, $\Lambda$

54:for $k=0$ to $\hat{N}-1$ do

5: Initialize: $\tilde{\xi^{k}}$

66: for $\iota_{i}=0$ to $I_{i}-1$ do {Inner optimization Loop}

77: $\tilde{\Xi^{k}}=\\{\tilde{\xi^{k}},X_{r}\\},\quad r=0,...,\omega$

8: $\tilde{D^{k}_{r}}=\Gamma(\tilde{\Xi^{k}_{r}}),\quad r=0,...,\omega$ {Eq.
(30d)}

9: $\tilde{Q^{k}}\leftarrow$ DrawValuesAtSensorLocations($\tilde{D^{k}}$) {Eq.
(30c)}

10: $\mathfrak{L}=MSE(\chi^{k},\tilde{Q^{k}})$

11: $\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}\leftarrow$
Backprop($\mathfrak{L}$)

12:
$\tilde{\xi^{k}}=\tilde{\xi^{k}}-\eta_{i}\frac{\partial\mathfrak{L}}{\partial\tilde{\xi^{k}}}$

813: end for

14: $\Xi^{k}_{r}=\\{\tilde{\xi^{k}},X_{r}\\},\quad r=0,...,\omega$

15: $D^{k}_{r}=\Gamma(\Xi^{k}_{r}),\quad r=0,...,\omega$ {Eq. (28)}

916:end for

1017:Output: Predicted states $\\{D^{ki}_{m}\ |\ k=0,...,N-1,\
i=0,...,\gamma-1,\ m=0,...,M-1\\}$.

Algorithm 4 Testing ENSERS

####  4.1.1 Physics-based loss function

The physics-based loss function for continuous space models differs from
discrete formulation due to how RHS of Eq. (1) i.e. $F$ is evaluated. In this
case, we use automatic differentiation for calculating gradients w.r.t.
coordinates used in $F$. For example in first case of Allen–Cahn equation
$\hat{F}$ is evaluated as:

$\hat{F_{r}}=0.0001\frac{\partial^{2}u}{\partial\hat{X_{r}}^{2}}+5u^{3}-5u,\quad
r=0,...,\omega$ (32)

where $\hat{X_{r}}$ is coordinate vector concatenated with optimized reduced
state vector at end of inner optimization loop, see line $18$ in algorithm
(3). In second case of Convection-diffusion equation $\hat{F}$ is evaluated
as:

$\hat{F_{r}}=a(x,y)\frac{\partial
u}{\partial\hat{X_{r}^{0}}}+b(x,y)\frac{\partial
u}{\partial\hat{X_{r}^{1}}}+c\frac{\partial^{2}u}{\partial\hat{X_{r}^{0}}^{2}}+d\frac{\partial^{2}u}{\partial\hat{X_{r}^{1}}^{2}},\quad
r=0,...,\omega.$ (33)

where $\hat{X_{r}^{0}},\hat{X_{r}^{1}}$ are x and y coordinate respectively,
$a,b,c,d$ are defined in Eq. 36.

###  4.2 Experiment: Allen–Cahn equation

We consider the Allen–Cahn equation along with periodic boundary conditions.
The Allen–Cahn equation is a well-known equation from the area of reaction-
diffusion systems. It describes the process of phase separation in
multicomponent alloy systems, including order-disorder transitions.

$\displaystyle u_{t}-0.0001u_{xx}+5u^{3}-5u=0,\quad x\in[-1,1],t\in[0,1],$
(34) $\displaystyle u(0,x)=x^{2}cos(\pi x),$ $\displaystyle u(t,-1)=u(t,1),$
$\displaystyle u_{x}(t,-1)=u_{x}(t,1).$

####  4.2.1 Data-set and Model Parameters

Data-set is generated by simulating the Allen–Cahn equation (34) using
conventional spectral methods. Starting from an initial condition
$u(0,x)=x^{2}cos(\pi x)$ and assuming periodic boundary conditions
$u(t,-1)=u(t,1)$ and $u_{x}(t,-1)=u_{x}(t,1)$, we integrated Eq. (34) up to a
final time $t=1.0$ using the Chebfun package [49] with a spectral Fourier
discretization with 512 modes and a fourth-order explicit Runge–Kutta temporal
integrator with time-step $\Delta t=10^{-5}$. For more details on the data-set
see [41]. Plots depicting different noise levels used for sensor measurement
during testing are shown in Fig. 9.

Figure 9: Data plots of Allen–Cahn equation for visualizing different noise levels. SNRdB represents signal to noise ratio in decibels. Table 9: Data-set parameters for Allen–Cahn equation. $\gamma$ |  $M$ |  $L$ |  $z$ |  $\omega$ |  domain
---|---|---|---|---|---
5 |  1 |  50 |  1 |  128 |  $128$

The network considered in the continuous formulation has significantly fewer
weights compared to the discrete one. This is because the output of the
network is predicted at a single point. Network architecture is shown in table
10. The input size of the network is equal to the sum of the size of the
reduced state $\varsigma$ and the dimension of coordinates. For training we
use 10 data nodes $p$ which is $7.81\%$ of $128$ total nodes $\omega$. Similar
to previous cases inner loop learning rate $\eta_{i}$ and physics loss penalty
$\zeta$ are increased linearly with each epoch. Other training and testing
hyperparameters are shown in table 11 and 12 respectively.

Table 10: Network architecture of proposed Ensers for Allen–Cahn equation. Layer |  Input |  Output |  Activation
---|---|---|---
FC |  7 |  128 |  Tanh
FC |  128 |  128 |  Tanh
FC |  128 |  128 |  Tanh
FC |  128 |  128 |  Tanh
FC |  128 |  $M*\gamma$ |  linear
Table 11: Training hyper parameters of Ensers for Allen–Cahn equation. $\eta_{o}$ |  $\eta_{i0}$ |  $\hat{\eta_{i}}$ |  $\beta$ |  $I_{o}$ |  $I_{i}$ |  $\zeta_{0}$ |  $\hat{\zeta}$ |  $p$ |  $h$ |  $N$ |  $\varsigma$
---|---|---|---|---|---|---|---|---|---|---|---
0.001 |  0.005 |  $6\mathrm{e}{-5}$ |  10 |  1201 |  8 |  0.01 |  $2\mathrm{e}{-5}$ |  10 |  10 |  40 |  6
Table 12: Testing hyper parameters of Ensers for Allen–Cahn equation. $\eta_{i}$: inner learning rate, $I_{i}$: inner iterations $\eta_{i}$ |  $I_{i}$ |  $p$ |  $\hat{N}$
---|---|---|---
0.1 |  50 |  [6, 16] |  12

####  4.2.2 Results and Discussions

Fig. 10 shows the prediction of Ensers with 16 sensors. A noticeable benefit
of the continuous formulation is that prediction is smooth compared to
discrete cases. Fig. 11 shows violin plot of L2 error defined in Eq. (19)
between target and prediction with various noise levels in sensor measurement
for the test case. In Eq. (19), $m=0$ and $\gamma^{*}=2$. Similar to discrete
cases, the model is robust to noise in measurements. Fig. 11 shows a marginal
increase in error with noise in the case of 16 sensors. Also, error bars with
6 sensors are similar to 16 sensors with low noise and increase slowly with
noise.

Figure 10: Ensers predictions of a Allen–Cahn equation test case with 16
sensors and no noise.

Figure 11: Violin plot representing error $\epsilon$ distribution in
prediction vector of data, for different noise levels in sensor measurements
for Allen–Cahn equation test case.

###  4.3 Experiment: Convection-diffusion equations

Next we consider a 2-dimensional linear variable-coefficient convection-
diffusion equation on $\omega=[0,2\pi]\times[0,2\pi]$,

$u_{t}=a(x,y)u_{x}+b(x,y)u_{y}+cu_{xx}+du_{yy}\ (t,x,y)\in[0,0.2]\times\omega$
(35)

$\displaystyle a(x,y)=0.5(cos(y)+x(2\pi-x)sin(x))+0.6$ (36) $\displaystyle
b(x,y)=2(cos(y)+sin(x))+0.8$ $\displaystyle c=0.2,d=0.3$

Convection-diffusion equations are classical PDEs that are used to describe
physical phenomena where particles, energy, or other physical quantities are
transferred inside a physical system due to two processes namely diffusion and
convection. These equations are widely applied in many scientific areas and
industrial fields, such as pollutants dispersion in rivers or atmosphere,
solute transferring in a porous medium, and oil reservoir simulation. We
consider variables convection and diffusion coefficients Eq. (36) in this
experiment.

####  4.3.1 Data-set and Model Parameters

Data is generated by solving the problem (35) using a high precision numerical
scheme with a pseudo-spectral method for spatial discretization and 4th order
Runge-Kutta for temporal discretization (with time step size $\delta t=0.01$).
We assume periodic boundary conditions and the initial value

$u(x,y,0)=\sum_{|k|,|l|\leq N}\lambda_{k,l}cos(kx+ly)+\gamma_{k,l}sin(kx+ly)$
(37)

where $N=9,\lambda_{k,l},\gamma_{k,l}\sim\mathcal{N}(0,0.02)$, and k and l are
chosen randomly. For more details on the data-set see [50]. Noisy data used
for sensor measurement during testing is shown in Fig. 12 with various noise
levels at a particular time. We use a simulation of length 40 as shown in
table 13 along with other data parameters.

Figure 12: Data plots of Convection-diffusion equation for visualizing different noise levels. SNRdB represents signal to noise ratio in decibels. Table 13: Data-set parameters for Convection-diffusion equation. $\gamma$ |  $M$ |  $L$ |  $z$ |  $\omega$ |  domain
---|---|---|---|---|---
5 |  1 |  40 |  1 |  4096 |  $64\times 64$
Table 14: Network architecture of proposed Ensers for Convection-diffusion equation. Layer |  Input |  Output |  Activation
---|---|---|---
FC |  8 |  128 |  Tanh
FC |  128 |  256 |  Tanh
FC |  256 |  256 |  Tanh
FC |  256 |  128 |  Tanh
FC |  128 |  $M*\gamma$ |  linear
Table 15: Training hyper parameters of Ensers for Convection-diffusion equation. $\eta_{o}$ |  $\eta_{i0}$ |  $\hat{\eta_{i}}$ |  $\beta$ |  $I_{o}$ |  $I_{i}$ |  $\zeta_{0}$ |  $\hat{\zeta}$ |  $p$ |  $h$ |  $N$ |  $\varsigma$
---|---|---|---|---|---|---|---|---|---|---|---
0.0005 |  0.005 |  $5\mathrm{e}{-5}$ |  11 |  1001 |  10 |  0.005 |  $1\mathrm{e}{-5}$ |  32 |  1024 |  33 |  6
Table 16: Testing hyper parameters of Ensers for Convection-diffusion equation. $\eta_{i}$: inner learning rate, $I_{i}$: inner iterations $\eta_{i}$ |  $I_{i}$ |  $p$ |  $\hat{N}$
---|---|---|---
0.2 |  50 |  [4, 16] |  12

Network architecture is shown in table 14. Similar to the previous case a deep
network with $5$ layers and a small output size is used. The activation
function is Tanh, which also has a continuous first derivative. For training,
we use $p=1024$ training nodes i.e. $25\%$ of total nodes. Training and
testing hyperparameters are shown in tables k and ll respectively. Inner loop
learning rate $\eta_{i}$ and physics loss penalty $\zeta$ are increased
linearly with each epoch.

####  4.3.2 Results and Discussions

Figure 13: Ensers predictions of a Convection-diffusion equation test case.
(Top to bottom) Target solution, Ensers prediction, L1error.

Figure 14: Violin plot representing error $\epsilon$ distribution in
prediction vector of x,y velocity and pressure for different noise levels in
sensor measurements for Convection-diffusion equation test case.

Fig. 13 shows the prediction of the Ensers with 16 sensors. The model is able
to predict the state accurately and with little distortion. Fig. 14 shows
violin plot of L2 error defined in Eq. (19) between target and prediction with
various noise levels in sensor measurement for the test case. In Eq. (19),
$m=0$ and $\gamma^{*}=2$ i.e. we use middle prediction out $\gamma$ time
steps. We see the model is robust to noise as the mean error in the plots are
close for cases of no noise($SNR_{dB}=None$) and high noise($SNR_{dB}=10$) for
16 sensors case and increases slightly noise for 4 sensor cases. Training
continuous model for two dimensions takes more time than discrete because of
the increased number of collocation points and multiple automatic
differentiation required for calculating loss function.

##  5 Conclusions

In this work, we develop a novel technique to learn a deep learning model from
spare moving training labels for which model reduction is nontrivial. The
method uses an implicit optimization layer for minimizing the energy of the
solution implemented through the technique of unrolled differentiation. We
proposed two formulations based on this technique for discrete and continuous
prediction in space. For learning from spare training labels we included a
physics-based loss function calculated via convolutional filters in the
discrete formulation and via automatic differentiation in the continuous
formulation. Where most deep learning-based methods assume fixed sensors
Ensers is capable of predicting full states given a varied number of sensors
at random locations.

We demonstrate the model performance using two-fluid problems of 2D coupled
Burgers’ equation and Flow Past Cylinder for the discrete case. For the
continuous case, we used two problems namely Allen–Cahn equation and the
Convection-diffusion equation. Model is shown to be robust against noisy
sensor measurements. Future work can be aimed at quantifying uncertainty in
such networks. Another direction can be to train networks for future states
only using initial and boundary conditions.

## Acknowledgements

SC acknowledges the financial support received from Science and Engineering
Research Board (SERB) via project no. SRG/2021/000467 and Indian Institute Of
Technology–Delhi in form of seed grant.

## Reproducibility

The codes associated with the paper will be released on acceptance.

##  Appendix A Convolution operators for gradient and laplacian terms

Sobel Filter used to estimate 1st-order gradient is:

$E=\left[\begin{array}[]{ccccc}1&-8&0&8&-1\\\ 2&-16&0&16&-2\\\
3&-24&0&24&-3\\\ 2&-16&0&16&-2\\\ 1&-8&0&8&-1\\\ \end{array}\right]$ (38a)

$\frac{\partial}{\partial x}=E\times\frac{1}{9*12\delta x}$ (38b)

Filter used to estimate laplacian is:

$\frac{\partial}{\partial y}=E^{T}\times\frac{1}{9*12\delta y}$ (38c)

$\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}filter=\left[\begin{array}[]{ccccc}0&0&-1&0&0\\\ 0&0&16&0&0\\\
-1&16&-60&16&-1\\\ 0&0&16&0&0\\\ 0&0&-1&0&0\\\
\end{array}\right]\times\frac{1}{12\delta x\delta y}$ (39)

## References

  * Cordier [2011] L. Cordier, Springer Vienna, 2011\. doi:10.1007/978-3-7091-0758-4_1.
  * Semaan et al. [2019] R. Semaan, M. Y. El Sayed, R. Radespiel, Springer International Publishing, Cham, 2019.
  * Sankaran et al. [2012] S. Sankaran, M. Esmaily Moghadam, A. M. Kahn, E. E. Tseng, J. M. Guccione, A. L. Marsden,  Patient-specific multiscale modeling of blood flow for coronary artery bypass graft surgery,  Ann Biomed Eng 40 (2012) 2228–2242.
  * Yakhot et al. [2007] A. Yakhot, T. Anor, G. E. Karniadakis,  A reconstruction method for gappy and noisy arterial flow data,  IEEE Trans Med Imaging 26 (2007) 1681–1697.
  * Kissas et al. [2020] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, P. Perdikaris,  Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4d flow mri data using physics-informed neural networks,  Computer Methods in Applied Mechanics and Engineering 358 (2020) 112623\. URL: https://www.sciencedirect.com/science/article/pii/S0045782519305055. doi:https://doi.org/10.1016/j.cma.2019.112623.
  * Graziano et al. [2016] M. D. Graziano, M. D’Errico, G. Rufino,  Ship heading and velocity analysis by wake detection in sar images,  Acta Astronautica 128 (2016) 72–82. URL: https://www.sciencedirect.com/science/article/pii/S0094576516300303. doi:https://doi.org/10.1016/j.actaastro.2016.07.001.
  * Kalnay [2002] E. Kalnay, Atmospheric Modeling, Data Assimilation and Predictability, Cambridge University Press, 2002. doi:10.1017/CBO9780511802270.
  * Ong et al. [2015] F. Ong, M. Uecker, U. Tariq, A. Hsiao, M. T. Alley, S. S. Vasanawala, M. Lustig,  Robust 4d flow denoising using divergence-free wavelet transform,  Magnetic Resonance in Medicine 73 (2015) 828–842. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.25176. doi:https://doi.org/10.1002/mrm.25176. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.25176.
  * Callaghan and Grieve [2017] F. M. Callaghan, S. M. Grieve,  Spatial resolution and velocity field improvement of 4D-flow MRI,  Magn Reson Med 78 (2017) 1959–1968.
  * Fathi et al. [2018] M. F. Fathi, A. Bakhshinejad, A. Baghaie, D. Saloner, R. H. Sacho, V. L. Rayz, R. M. D’Souza,  Denoising and spatial resolution enhancement of 4D flow MRI using proper orthogonal decomposition and lasso regularization,  Comput Med Imaging Graph 70 (2018) 165–172.
  * Bishop [2006] C. M. Bishop, Pattern recognition and machine learning, springer, 2006.
  * Reif et al. [1999] K. Reif, S. Gunther, E. Yaz, R. Unbehauen,  Stochastic stability of the discrete-time extended kalman filter,  IEEE Transactions on Automatic control 44 (1999) 714–728.
  * Wan et al. [2001] E. A. Wan, R. Van Der Merwe, S. Haykin,  The unscented kalman filter,  Kalman filtering and neural networks 5 (2001) 221–280.
  * Tu et al. [2012] J. Tu, J. Griffin, A. Hart, C. Rowley, L. Cattafesta, L. Ukeiley,  Integration of non-time-resolved piv and time-resolved velocity point sensors for dynamic estimation of velocity fields,  Experiments in Fluids 54 (2012). doi:10.1007/s00348-012-1429-7.
  * Schmid and Sesterhenn [2008] P. Schmid, J. Sesterhenn,  Dynamic Mode Decomposition of numerical and experimental data,  in: APS Division of Fluid Dynamics Meeting Abstracts, volume 61 of APS Meeting Abstracts, 2008, p. MR.007.
  * ROWLEY et al. [2009] C. W. ROWLEY, I. MEZIĆ, S. BAGHERI, P. SCHLATTER, D. S. HENNINGSON,  Spectral analysis of nonlinear flows,  Journal of Fluid Mechanics 641 (2009) 115–127. doi:10.1017/S0022112009992059.
  * Buffoni et al. [2008] M. Buffoni, S. Camarri, A. Iollo, E. Lombardi, M. Salvetti,  A non-linear observer for unsteady three-dimensional flows,  Journal of Computational Physics 227 (2008) 2626 – 2643. URL: http://www.sciencedirect.com/science/article/pii/S0021999107004858. doi:https://doi.org/10.1016/j.jcp.2007.11.005.
  * Adrian [1975] R. Adrian, On the role of conditional averages in turbulence theory., 1975.
  * Guezennec [1989] Y. G. Guezennec,  Stochastic estimation of coherent structures in turbulent boundary layers,  Physics of Fluids A: Fluid Dynamics 1 (1989) 1054–1060. URL: https://doi.org/10.1063/1.857396. doi:10.1063/1.857396. arXiv:https://doi.org/10.1063/1.857396.
  * Ewing and Citriniti [1999] D. Ewing, J. H. Citriniti,  Examination of a lse/pod complementary technique using single and multi-time information in the axisymmetric shear layer,  in: J. N. Sørensen, E. J. Hopfinger, N. Aubry (Eds.), IUTAM Symposium on Simulation and Identification of Organized Structures in Flows, Springer Netherlands, Dordrecht, 1999, pp. 375–384.
  * Naguib et al. [2001] A. M. Naguib, C. E. Wark, O. Juckenhöfel,  Stochastic estimation and flow sources associated with surface pressure events in a turbulent boundary layer,  Physics of Fluids 13 (2001) 2611–2626. URL: https://doi.org/10.1063/1.1389284. doi:10.1063/1.1389284. arXiv:https://doi.org/10.1063/1.1389284.
  * Bonnet et al. [1994] J. P. Bonnet, D. R. Cole, J. Delville, M. N. Glauser, L. S. Ukeiley,  Stochastic estimation and proper orthogonal decomposition: Complementary techniques for identifying structure,  Experiments in Fluids 17 (1994) 307–314. URL: https://doi.org/10.1007/BF01874409. doi:10.1007/BF01874409.
  * Pinier et al. [2007] J. Pinier, J. Ausseur, M. Glauser, H. Higuchi,  Proportional closed-loop feedback control of flow separation,  Aiaa Journal - AIAA J 45 (2007) 181–190. doi:10.2514/1.23465.
  * Adrian [1979] R. Adrian,  Conditional eddies in isotropic turbulence,  Physics of Fluids 22 (1979). doi:10.1063/1.862515.
  * Tung and Adrian [1980] T. C. Tung, R. J. Adrian,  Higher-order estimates of conditional eddies in isotropic turbulence,  The Physics of Fluids 23 (1980) 1469–1470. URL: https://aip.scitation.org/doi/abs/10.1063/1.863130. doi:10.1063/1.863130. arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.863130.
  * Kumar et al. [2021] Y. Kumar, P. Bahl, S. Chakraborty, State estimation with limited sensors – a deep learning based approach, 2021. arXiv:2101.11513.
  * Erichson et al. [2019] N. B. Erichson, L. Mathelin, Z. Yao, S. L. Brunton, M. W. Mahoney, J. N. Kutz, Shallow learning for fluid flow reconstruction with limited sensors and limited data, 2019. arXiv:1902.07358.
  * Nair and Goza [2020] N. J. Nair, A. Goza,  Leveraging reduced-order models for state estimation using deep learning,  Journal of Fluid Mechanics 897 (2020) R1. doi:10.1017/jfm.2020.409.
  * Gao et al. [2021] H. Gao, L. Sun, J.-X. Wang,  Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels,  Physics of Fluids 33 (2021) 073603. URL: https://doi.org/10.1063/5.0054312. doi:10.1063/5.0054312. arXiv:https://doi.org/10.1063/5.0054312.
  * Grefenstette et al. [2019] E. Grefenstette, B. Amos, D. Yarats, P. M. Htut, A. Molchanov, F. Meier, D. Kiela, K. Cho, S. Chintala,  Generalized inner loop meta-learning,  arXiv preprint arXiv:1910.01727 (2019).
  * Amos [2019] B. Amos, Differentiable Optimization-Based Modeling for Machine Learning, Ph.D. thesis, Carnegie Mellon University, 2019\.
  * Goodfellow et al. [2013] I. Goodfellow, M. Mirza, A. Courville, Y. Bengio,  Multi-prediction deep boltzmann machines,  Advances in Neural Information Processing Systems (2013).
  * Stoyanov et al. [2011] V. Stoyanov, A. Ropson, J. Eisner,  Empirical risk minimization of graphical model parameters given approximate inference, decoding, and model structure,  in: G. Gordon, D. Dunson, M. Dudík (Eds.), Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, volume 15 of Proceedings of Machine Learning Research, PMLR, Fort Lauderdale, FL, USA, 2011, pp. 725–733. URL: https://proceedings.mlr.press/v15/stoyanov11a.html.
  * Brakel et al. [2013] P. Brakel, D. Stroob, t, B. Schrauwen,  Training energy-based models for time-series imputation,  Journal of Machine Learning Research 14 (2013) 2771–2797. URL: http://jmlr.org/papers/v14/brakel13a.html.
  * Lecun et al. [2006] Y. Lecun, S. Chopra, R. Hadsell, A tutorial on energy-based learning, 2006.
  * Utama et al. [2018] D. Utama, A. N., M. Iqbal,  An optimal generic model for multi-parameters and big data optimizing: a laboratory experimental study,  Journal of Physics: Conference Series 978 (2018) 012045. doi:10.1088/1742-6596/978/1/012045.
  * Belanger et al. [2017] D. Belanger, B. Yang, A. McCallum, End-to-end learning for structured prediction energy networks, 2017. arXiv:1703.05667.
  * Metz et al. [2016] L. Metz, B. Poole, D. Pfau, J. Sohl-Dickstein,  Unrolled generative adversarial networks (2016).
  * Johnson et al. [2016] M. J. Johnson, D. K. Duvenaud, A. Wiltschko, R. P. Adams, S. R. Datta,  Composing graphical models with neural networks for structured representations and fast inference,  in: D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, R. Garnett (Eds.), Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016. URL: https://proceedings.neurips.cc/paper/2016/file/7d6044e95a16761171b130dcb476a43e-Paper.pdf.
  * Jordan-Squire [2015] C. Jordan-Squire,  Convex optimization over probability measures,  2015\.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational Physics 378 (2019) 686–707. doi:https://doi.org/10.1016/j.jcp.2018.10.045.
  * Zhu et al. [2019] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, P. Perdikaris,  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data,  Journal of Computational Physics 394 (2019) 56–81. URL: https://www.sciencedirect.com/science/article/pii/S0021999119303559. doi:https://doi.org/10.1016/j.jcp.2019.05.024.
  * Geneva and Zabaras [2020] N. Geneva, N. Zabaras,  Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks,  Journal of Computational Physics 403 (2020) 109056. URL: https://www.sciencedirect.com/science/article/pii/S0021999119307612. doi:https://doi.org/10.1016/j.jcp.2019.109056.
  * Sobel and Feldman [1973] I. Sobel, G. Feldman,  A 3×3 isotropic gradient operator for image processing,  Pattern Classification and Scene Analysis (1973) 271–272.
  * Logg et al. [2012] A. Logg, K.-A. Mardal, G. N. Wells, et al., Automated Solution of Differential Equations by the Finite Element Method, Springer, 2012. doi:10.1007/978-3-642-23099-8.
  * Glorot et al. [2011] X. Glorot, A. Bordes, Y. Bengio,  Deep sparse rectifier neural networks,  in: G. Gordon, D. Dunson, M. Dudík (Eds.), Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, volume 15 of Proceedings of Machine Learning Research, PMLR, Fort Lauderdale, FL, USA, 2011, pp. 315–323. URL: https://proceedings.mlr.press/v15/glorot11a.html.
  * Jasak [2009] H. Jasak,  Openfoam: open source cfd in research and industry,  International Journal of Naval Architecture and Ocean Engineering 1 (2009) 89–94.
  * Nithin Adidela and Sudhakar [????] R. S. Nithin Adidela, Y. Sudhakar, Laminar flow over a circular cylinder simulation with openfoam v7, https://github.com/nithinadidela/circular-cylinder, ????
  * Driscoll et al. [2014] T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun guide, 2014\.
  * Long et al. [2018] Z. Long, Y. Lu, X. Ma, B. Dong,  PDE-net: Learning PDEs from data,  in: J. Dy, A. Krause (Eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, PMLR, 2018, pp. 3208–3216. URL: https://proceedings.mlr.press/v80/long18a.html.