# Taming Waves: A Physically-Interpretable Machine Learning
Framework for Realizable Control of Wave Dynamics

Tristan Shah1, Feruza Amirkulova2, Stas Tiomkin1,∗
1Computer Engineering Department, 1Mechanical Engineering Department
Charles W. Davidson College Of Engineering, San Jose State University, CA, USA
{tristan.shah, feruza.amirkulova<EMAIL_ADDRESS>∗ Corresponding
Author

###### Abstract

Controlling systems governed by partial differential equations is an
inherently hard problem. Specifically, control of wave dynamics is challenging
due to additional physical constraints and intrinsic properties of wave
phenomena such as dissipation, attenuation, reflection, and scattering. In
this work, we introduce an environment designed for the study of the control
of acoustic waves by actuated metamaterial designs. We utilize this
environment for the development of a novel machine-learning method, based on
deep neural networks, for efficiently learning the dynamics of an acoustic PDE
from samples. Our model is fully interpretable and maps physical constraints
and intrinsic properties of the real acoustic environment into its latent
representation of information. Within our model we use a trainable perfectly
matched layer to explicitly learn the property of acoustic energy dissipation.
Our model can be used to predict and control scattered wave energy. The
capabilities of our model are demonstrated on an important problem in
acoustics, which is the minimization of total scattered energy. Furthermore,
we show that the prediction of scattered energy by our model generalizes in
time and can be extended to long time horizons. We make our code repository
publicly available.

###### Index Terms:

Data-driven control of PDE, deep neural networks, wave equation, forward
modelling, inverse design, perfectly matched layers, meta-materials,
representation learning, differentiable simulation.

##  I Introduction

Figure 1: Overall pipeline of our algorithm which shows the flattening of 2D
partial information into the initial conditions of a 1D wave simulation.
Information supplied to the wave encoder comes from the stack of ”wave images”
shown in blue at timesteps: $t_{\tau-2},t_{\tau-1},t_{\tau}$ is input to the
CNN. The source of acoustic energy is a point which appears in white.
Concentric rings emanate from this point and scatter upon the design shown in
grey. Information about the sequence of designs is input to the design encoder
and used to construct a speed of sound function $c^{z}(x,t)$ which influences
the latent solution. The relevant information is combined and used to
propagate the latent dynamics $\mathbf{z}(x,t_{\tau})$. Finally, scattered
energy is computed in latent space and used to approximate scattered energy in
the real environment.

Partial Differential Equations (PDE)s are useful mathematical tools for
describing how physical systems behave as a continuum. Long established PDEs
such as the classical wave equation, Navier-Stokes equation, heat equation,
and the Schrödinger equation govern systems which are the foundation of our
understanding of physics. When combined together and with other observations
of physical systems they can be used to describe complex phenomenon such as
wild fires [1], the global climate [2], or spread of disease [3].

Control over systems governed by PDEs usually requires the existence of an
analytical solution which is dependent on control parameters to influence the
system in a desirable direction [4, 5]. With the rise of powerful sample based
learning methods it is now possible to influence PDEs without deriving a
solution [6, 7], unlocking control of a much larger suite of problems in
physics than was previously accessible. However, employing algorithms from
machine learning to solve problems in PDE control is not always
computationally tractable. PDEs are difficult to model with neural networks
(NN) due to their infinite-dimensional tensor (function valued) states and
sensitive dynamics [8]. While machine learning (ML) methods have been shown to
be very successful at controlling simple scenarios, it is important to search
for ways to reduce the complexity of this type of problem for the development
of truly scalable methods.

Establishing control over the acoustic wave equation is important because it
provides pathways to developing new technologies. These technologies range
from acoustic cloaking for airborne and aquatic stealth, wave steering devices
for seismic wave manipulation, to ’super-focusing’ devices for high-resolution
ultrasound imaging and surgery. There exist some known analytical solutions to
this PDE [9, 5] which have been used to solve these problems. However, in our
work, we pursue an optimal control method that leverages ML.

In our research, we design an algorithm for prediction and control of acoustic
PDEs that could eventually be realized by a physical device. We allow our
model to observe partial information from its environment which could
realistically be collected by sensors. Furthermore, we consider a control
mechanism capable of being constructed and deployed. We also make advancements
in interpretability of dynamical models learnt from samples. We design our
model with a latent space structured by the wave equation, which allows us to
visualize and draw meaning from its state. We use the interpretable structure
of the wave equation to introduce a novel and trainable energy dissipation
layer derived from domain knowledge in acoustics. Overall our work represents
an advancement in the field of PDE control which will enable future works to
build more interpretable and robust artificial intelligence.

####  I-1 Contributions

The results of our work can be distilled into three major contributions:

  * •

Environment for studying the behavior of waves in free field domains in
response to a controller.

  * •

Algorithm for learning the dynamics of scattered energy in response to a
control signal.

  * •

A novel trainable dissipation layer which can learn to explicitly remove
energy from latent space.

We elaborate on our method in Section Method.

##  II Prior Work

In recent years, the control of systems governed by PDEs has become an
increasingly active topic in ML [5, 4, 7, 6, 10]. Research in this field falls
into two categories. One can be thought of as ”parameter optimization” of a
PDE solution. It requires a known solution to a PDE which can be used to form
an optimization objective based on some parameters. An algorithm such as
gradient-based optimization (GBO) [5] or Reinforcement Learning (RL) [4] is
used to optimize the design parameters. However, the requirement for a PDE
solution limits its applicability. The other category is dynamical control of
PDEs. This category does not require knowledge of a PDE solution. The method
takes a PDE, incorporates a control signal into its dynamics, and uses methods
like Reinforcement Learning (RL) [7] or Model Predictive Control (MPC) [6] to
select appropriate control actions. The limitation of this work is its use of
an RNN to learn PDE dynamics. RNNs lack interpretability [11] and tend to
overfit [12]. In the following text, we elaborate on these two categories and
overview the distinguishing features of the current work, which are explained
in detail in Section Method.

Parameter optimization has been used [5] to design arbitrarily shaped acoustic
cloaks. In their work, Cominelli et al. formulate an objective function using
a known solution to the wave equation. Their objective is dependent on
material properties in a region of a computational domain and predicts
scattered wave energy. Material properties are optimized using GBO such that
scattered energy is minimized, resulting in an acoustic cloak. Similarly, [4]
used the similar approach to design a cloaking device by optimizing the
position and radii of cylindrical scatterer configurations. Scattered wave
energy is minimized by RL agents which explore the design parameter space for
optimal configurations. Optimization of design parameters through known PDE
solutions is a robust and precise method but is limited by the existence of
such solutions. In many cases there are no analytical formulas to solve PDEs
making parameter optimization infeasible in many important scenarios.

Recently, a formulation that incorporates the dynamics of a PDE into an
optimal control problem has been explored [7]. The method in [7] has the
advantage of not requiring a known solution in order to be applicable;
widening the range of problems where it can be employed. The idea proposed in
[7] is the numerical integration of PDE dynamics while a control signal is
applied. The authors in [7] consider a particular PDE, the 1D Kuramoto-
Sivashinsky (KS) PDE, where control
$\phi:\Omega\times[0,T]\rightarrow\mathbb{R}^{m}$ is an additive force applied
to system variable $u:\Omega\times[0,T]\rightarrow\mathbb{R}^{n}$. In their
work $\phi$ is fixed over a time interval and is updated between subsequent
intervals. Equation (14) shows integration of the PDE over an interval where
$\phi$ is being applied concurrently with the dynamics resulting in a next
state $u_{\tau+1}$. This formulation is interesting because it unifies the
dynamics that we wish to control with the application of a control signal.
Methods from optimal control such as Reinforcement Learning (RL) or Model
Predictive Control (MPC) can be applied to manage the control signal for
optimization of any particular objective function given in Equation (13). An
important feature of their work is the efficient manner in which $\phi$ is
constructed. The authors in [7] consider $\phi$ as a superposition of four
Gaussian distributions with fixed variance and position. The effect of the
controller on $\phi$ is to modulate the intensity of each Gaussian through the
$\mathbf{a}_{\tau}[i]$ term. This procedure significantly reduces the
dimension of the control input from $\Omega$ to $\mathbb{R}^{4}$. Reduced
dimensional construction of a control signal is a highly relevant technique to
our work as we wish to control waves in a two spacial dimensions. Generating
control over a high resolution finite element grid may be computationally
intractable without a parameterized method. This work provides useful
formalism for framing a PDE as an optimal control problem but is limited by
the simplicity of the scenario it considers.

Other work [6] also considers dynamical control over the Navier-Stokes
equation. Their research focuses on designing a physically realizable control
algorithm. In their work, Bieker et al. trained a system of deep neural
networks to predict lift and drag coefficients on a configuration of rotating
cylinders over a time horizon. A recurrent neural network (RNN) was used to
approximate these quantities by taking in the state of the simulator and a
sequence of control inputs. Bieker et al. show that their RNN dynamics model
is able to make accurate predictions and performs optimal control with MPC.
There are two particular contributions made in this work which are relevant to
our research. Firstly, the model takes in information gathered from physically
realizable sensors rather than the full state of the fluid simulation. This
practice makes the model much more applicable in real world scenarios when
limited data is available. Secondly, the effect of the control input on the
PDE is locally sparse. In other words, the authors do not consider the ability
to influence the PDE at every point in the computational domain but only at
the interface between the cylinders and the ambient fluid. The effect of
sparse control is indirect compared to the alternative but more realistically
achievable with a physical device. This work is relevant to our research in
its proposal of parameterized construction of a control signal but is limited
to a relatively simple system which occurs in 1D. Additionally, the dynamics
model is not interpretable and does not incorporate domain knowledge as
helpful structure.

##  III Current Work

###  III-A Preliminaries

####  III-A1 Acoustic Metamaterials

Metamaterials are engineered materials that exert an influence over acoustic
waves in the extraordinary form of sound light cloaking, super-focusing,
negative refraction, negative mass density, and negative Poisson’s ratio to
name a few. The realization of novel acoustic materials with enhanced
capabilities provides the potential expanded functionality over a range of
applications including cloaking [13, 14], wave steering, and focusing. In this
work, we design metastructures focusing on sound direction control such as
cloaking by suppressing the scattered wave energy.

####  III-A2 Acoustic Multiple Scattering

Multiple scattering [14] from fluid media with dissimilar material properties
is of great importance to various applications in acoustics and elastic wave
propagation. For instance, biomedical ultrasound and non-distractive testing
make use of acoustic field scattering for localization, imaging, and object
detection. In this work, we consider a 2D planar motion and multiple
scattering by cylindrical structures in acoustic media. For demonstration
purposes, we consider a cluster of fluid cylinders located inside another
fluid medium. A realization of such composite fluid media can be achieved by
using effective medium theory [15]. The governing equation for sound waves is
a wave equation:

$\ddot{u}=c^{2}\nabla^{2}u$ (1)

where $\nabla^{2}$ is the Laplacian, $t$ is time, $u$ is the physical quantity
such as pressure or velocity potential, and $c$ is the sound speed in a fluid
medium. The boundary conditions at the interface around each scatterer are the
continuity of normal components of the particle velocity and normal stress
components.

###  III-B Free Field Scattered Wave Environment

The environment we introduce simulates interaction between acoustic waves and
metamaterial designs actuated by a controller. It serves as a platform for
studying control of wave dynamics in open space. A key feature of our
environment is its capability to simulate a ”free field” infinite spacial
domain. Many problems in acoustics rely on the assumption that the space in
which waves propagate is infinite [16]; no reflections occur along the
boundary of the domain. We make our environment conveniently accessible to the
ML community through an interface in the Julia Programming Language. Our goal
is to promote the use of problems in acoustics as a testing ground for PDE
control algorithms. Our environment and its implementation offer the technical
tools required to further explore this area and establish stronger connections
between acoustics and ML. The structure of our environment and important
features are described in the following text.

####  III-B1 Excitation Force

Our environment differs from environments in classical control because it has
a time-dependent excitation which continously adds energy to the system.
Moreover, the excitation is not directly controllable by the agent we consider
in this work. Sources in our environment can have an arbitrary shape defined
by a function $f(\mathbf{x})$ in Equations (5, 6) which oscillate with
frequency $\omega$. This flexible definition of the source function enables
our environment to simulate different scenarios such as plane waves or point
sources. It is important to note that if an excitation exists in a closed
environment the energy level will constantly increase. In order to study how
acoustic waves interact with a non-linearity in isolation from energy buildup,
dissipation must be incorporated [16]. Without such dissipation, energy will
reflect on the boundaries of the domain multiple times and interfere with
scattered energy in the area of interest.

####  III-B2 Perfectly Matched Layer

One of the most important distinguishing features of our environment is its
free field capability. We simulate infinite space in two dimensions using a
Perfectly Matched Layer (PML). A PML is an artificial material placed around
the boundary region of a computational domain that exponentially dampens waves
as they propagate [17]. Dampening waves causes their energy to approach zero
and never return to the domain. The main idea behind a PML is to apply complex
coordinate stretching in the frequency domain as waves enter the PML region.
The transformation shown in Equation 2 is applied to a wave equation in 1D and
scales spacial derivatives to become smaller farther into the PML. This
transformation can be applied to wave equations in any dimension as long as
the coordinate stretching procedure is applied to each spacial partial
derivative. The application of the PML transformation to the two-dimensional
wave equation results in fields $\psi_{x}$, $\psi_{y}$, and $\gamma$ which
dissipate energy along the $x$ and $y$ axes and corners respectively. The
evolution of these fields in time is determined by Equations (7-9). The
resulting dissipative fields influence the solution variables $u$, $v_{x}$,
and $v_{y}$ in Equations (4-6). Functions $\sigma_{x}$ and $\sigma_{y}$
determine the shape of the PML in their respective axes in 2D. These functions
have a value of 0.0 everywhere in the domain and slope up towards 1.0 inside
the boundary region:

$\frac{\partial}{\partial x}\rightarrow\frac{\partial}{\partial
x}\cdot\frac{1}{1+i\frac{\sigma_{x}}{\omega}}.$ (2)

####  III-B3 Design Representation

The particular type of metamaterial designs we implement in our environment
are fluid scatterers. Specifically, we consider cylindrical scatterers with
radial symmetry; though in principle other shapes are achievable in this
framework. The scatterers influence the behavior of an incident wave by
reflecting a portion of its energy. This phenomenon occurs through a
difference in material properties at the interface between the scatterer and
the ambient domain. A benefit of using such scatterers is their influence on
the wave environment is physically plausible and could be achieved by real
devices. This method of actuating the environment is spatially sparse because
the control cannot affect the wave at every point in the field. Sparse control
in our work is in contrast to other works which consider an agent with the
ability to affect the PDE solution at every point in space [7]. Not only is
non-sparse control usually not realistic, but it is also extremely
computationally complex to optimize a controller to output function valued
actions at the same scale as the computational domain. In our framework, a
design can be expressed by an indicator function defined by Equation (3) which
assigns the appropriate speed of sound to every point in the domain at time
$t$

$c_{\tau}(\mathbf{x},t)=\begin{cases}1032\;m/s&\mathbf{x}\in d_{\tau}(t),\\\
1531\;m/s&else.\end{cases}$ (3)

####  III-B4 Control by Agent

Actions in this framework represent an agent’s control over the metamaterial
design. For the particular design that we consider in our experiments, the
agent actuates its control by adjusting the radii of each scatterer. This
adjustment is selected at a time step and applied over the time interval at a
rate of $500m/s$. We use linear interpolation between the parameters (radii)
of designs to represent a smooth transition between the design before and
after the application of the action. The indicator function for material
properties can be evaluated at any time in this interpolation to obtain the
information necessary to propagate the wave simulator. In our case, we
consider the ambient domain as water with a speed of sound of $1531m/s$ and
the scatterers as a fluid with three times the speed of sound in the air of
$1032m/s$. These values were chosen for the universal significance of water,
however our simulator can support any nonzero value of sound speed. The
propagation of waves in our simulator involves numerical integration of the
PML-transformed Equation (4) with $c_{\tau}$ as the indicator function. Once
the time of the simulation has reached the end of the interpolation, a new
action is applied to the design which forms the next interpolation interval
$c_{\tau+1}$.

$\displaystyle\dot{u}$
$\displaystyle=c_{\tau}^{2}(\partial_{x}v_{x}+\partial_{y}v_{y})+\psi_{x}+\psi_{y}$
$\displaystyle\qquad\qquad\qquad\qquad-(\sigma_{x}+\sigma_{y})u-\gamma,$ (4)
$\displaystyle\dot{v}_{x}$ $\displaystyle=\partial_{x}(u+f\cdot\sin(2\pi\omega
t))-\sigma_{x}v_{x},$ (5) $\displaystyle\dot{v}_{y}$
$\displaystyle=\partial_{y}(u+f\cdot\sin(2\pi\omega t))-\sigma_{y}v_{y},$ (6)
$\displaystyle\dot{\psi}_{x}$
$\displaystyle=c_{\tau}^{2}\sigma_{x}\partial_{y}v_{y},$ (7)
$\displaystyle\dot{\psi}_{y}$
$\displaystyle=c_{\tau}^{2}\sigma_{y}\partial_{x}v_{x},$ (8)
$\displaystyle\dot{\gamma}$ $\displaystyle=\sigma_{x}\sigma_{y}u.$ (9)

The full derivation of these equations appears in supplementary materials.

####  III-B5 Real Dynamics

The dynamics of the system variables describing displacement
$u=u(\mathbf{x},t)$ and velocity $v_{x}=v_{x}(\mathbf{x},t)$ and
$v_{y}=v_{y}(\mathbf{x},t)$ are defined by Equations (4, 5, 6). The
displacement field is influenced by the divergence of the velocity field
weighted by the speed of sound function $c_{\tau}$. Here, the operators
$\partial_{x}$ and $\partial_{y}$ represent partial derivatives
$\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$. Dissipation
from the PML is incorporated through the inclusion of fields
$\psi_{x}=\psi_{x}(\mathbf{x},t)$, $\psi_{y}=\psi_{y}(\mathbf{x},t)$, and
$\gamma(\mathbf{x},t)$. Similarly, the velocity fields are also updated with
the appropriate partial derivatives of the displacement field. The energy
source is incorporated by the addition of a forcing term $f=f(\mathbf{x})$
multiplied by a sinusoidal function on time. Together the full state of this
PDE system can be represented by Equation (10). To ensure the stability of
numerical integration of this system a fourth order Runge-Kutta scheme is
used:

$\mathbf{w}=\begin{bmatrix}u&v_{x}&v_{y}&\psi_{x}&\psi_{y}&\gamma\end{bmatrix}^{\mathsf{T}}.$
(10)

####  III-B6 Partial State Observation

In our work we develop an environment which simulates realistic interaction
between a controller and acoustic waves. In such a setting full observability
of the environment is not achievable. We consider information which can be
obtained through sensors such as readings of the total pressure field at
points in the simulation domain. The environment state representations that we
use are in the form of reduced resolution images ($128\times 128$) from the
full discretized finite element grid ($700\times 700$) which are derived from
the displacement of the total wave $u(\mathbf{x},t)$ over the domain evaluated
at times $t_{\tau-2}$, $t_{\tau-1}$, and $t_{\tau}$. We do not include images
of the underlying velocity components $v_{x}$ and $v_{y}$ because they would
not be directly observable in a real environment. In our case we used the
reduced resolution displacement images as input to our algorithm but in
principle an arbitrary sensor reading can be simulated to provide more
realistic observation.

####  III-B7 Training Signal

Various quantities can be computed from the state of our simulator and used as
reward or cost signals for a control algorithm. One possibility is to control
the energy in the scattered acoustic field. Scattered waves are acoustic
fields which result from the interaction between a wave and an obstacle. This
field can be obtained by Equation (12). In this definition the total field $u$
is the full waveform which accounts for an incident wave plus the scattering
it produces through interaction with a boundary condition. The incident field
$u_{inc}$ is a solution which accounts for the behavior of waves emanating
from the source $f(\mathbf{x})$ without the presence of a design that produces
scattering. In other words $u_{inc}$ assumes the entire domain to have a
constant rate of sound propagation. These fields are represented by two
separate wave dynamics which share the same grid size, ambient speed of sound,
initial conditions, and source. To obtain the scattered energy field we
subtract the two solutions. The scattered field is important because this
quantity is directly dependent on the behavior of the metamaterial design and
by extension the controller. The particular objective we demonstrate in this
paper is the reduction in total scattered energy $\sigma_{sc}$ which is
defined by the integral:

$\sigma_{sc}(t)=\int_{\mathbf{x}\in\Omega}u_{sc}(\mathbf{x},t)^{2}d\Omega$
(11)

where

$u_{sc}=u-u_{inc}.$ (12)

###  III-C Algorithm

We have developed a model that predicts scattered wave energy $\sigma_{sc}(t)$
as agent’s actions are applied to a design in our simulator. The full
schematic of the model is shown in Fig 1 and details information pathways from
the wave and design sequence to the latent dynamics and decoder. Our model has
several features we will highlight in this section including its latent
dynamics, interpretability, and ability to make predictions with partial
information.

####  III-C1 Reduced-Dimensionality Latent Wave Dynamics

The latent dynamics (LD) component of our model shown in the integration of
$\mathbf{z}$ in Fig 1 is a compressed, lower dimensional analog to real
dynamics present in our environment. LD share physical correspondence with
real dynamics (RD) in the form of learnable initial conditions, excitation
force, material properties, and energy dissipation. The key distinction of LD
is its lower dimensional state which is 1D. To accommodate the compression of
the RD into a lower dimensional space while maintaining as much information as
possible we assume that the LD must occur in a physically large space.

####  III-C2 Wave Encoder

In our model partial state information of the wave simulation in 2D is
processed by the wave encoder (CNN) and transformed into full information
initial conditions and parameters of 1D LD. We simulate two separate LD (total
and incident) and subtract their displacement fields to obtain a scattered
energy field through Equation (12). To simulate both fields, the initial
conditions of the LD are output from the wave encoder
$(u^{z}_{tot},v^{z}_{tot},u^{z}_{inc},v^{z}_{inc})$ which are discretized
functions over 1D space. The wave encoder also produces a learnable excitation
force function $f^{z}$ which oscillates at a known frequency of the source in
RD. This force continuously adds energy to the LD. Lastly, the wave encoder
learns the shape function of the PML in latent space, allowing it to control
how much energy is dissipated.

####  III-C3 Design Encoder

Nonuniform distribution of material properties in space influences the RD and
results in scattering of acoustic waves by affecting the rate at which sound
propagates: $c$. The same effect can be achieved in LD through a 1D function
which defines material properties in the latent space. It is a natural choice
to use information from the state of the metamaterial design to construct
latent material properties $c^{z}(x,t)$. We use a design encoder (MLP) to map
the parameters (radii of each scatterer in the configuration) of the sequence
of design states to corresponding discretized functions in 1D. Linear
interpolation between the functions produces a continuous latent speed of
sound function $c^{z}(x,t)$ which is defined over the same time interval as
the actions are applied in RD. This interpolation between latent
representations of designs can be visualized in Fig 1. The latent speed of
sound $c^{z}$ function is applied only to the latent total field $u^{z}_{tot}$
while $u^{z}_{inc}$ has a constant rate of sound propagation which is assumed
to be the same as the ambient speed of sound in RD.

####  III-C4 Sinusoidal Embedding

We use a linear superposition of $\sin$ waves to construct continuous latent
fields with a fixed number of frequencies $\sum_{n=1}^{N}b_{n}\sin(\frac{n\pi
x}{L})$. Intuition from this design choice comes from the general solution to
the wave equation in 1D in which initial conditions can be described in this
manner. Sinusoidal embedding is applied to the outputs of both wave and design
encoders. Instead of producing the full latent fields as their raw outputs,
encoders learn a vector of frequency coefficients $b_{n}^{(i)}$ for each
field. The sinusoidal basis functions have fixed endpoints of zero at the
boundaries of the latent domain.

####  III-C5 Forward Pass

The initial conditions excitation force, and PML functions are obtained from
the wave encoder $\mathbf{z}(t_{\tau})$ and the interpolation of material
properties function $c^{z}(x,t)$ from the design encoder. Conditions are
propagated by the LD to obtain a full solution $\mathbf{z}(\cdot)$. LD can be
solved with and without a PML to test the necessity for dissipation of energy
when predicting $\sigma_{sc}$ in free field. To our knowledge our work is the
first to ever use a PML to dissipate energy in the latent dynamics of a neural
network. The LD solution is obtained by propagating wave dynamics over the
time interval that actions are applied and $\sigma_{sc}(t)$ is measured in RD.
Similarly to the RD a Runge-Kutta integrator is used to integrate LD. The
predicted scattered energy signal is computed from $\mathbf{z}(\cdot)$ by
obtaining $u^{z}_{sc}$ through Equation (12) and integrating the energy over
1D space $\hat{\sigma}(t)_{sc}=\int_{x}u^{z}_{sc}(x,t)^{2}dx$. Prediction
error between true and predicted signals is measured through mean squared
error loss $\int_{t}(\sigma_{sc}(t)-\hat{\sigma}_{sc}(t))^{2}dt$.
Additionally, $\int_{t}(\sigma_{tot}(t)-\hat{\sigma}_{tot}(t))^{2}dt$ and
$\int_{t}(\sigma_{inc}(t)-\hat{\sigma}_{inc}(t))^{2}dt$ are minimized as
auxiliary constraints.

####  III-C6 Backwards Pass

Obtaining gradients of the latent solution with respect to the prediction
error relies on auto-differentiation (AD). These gradients can be computed
using reverse mode differentiation in deep learning frameworks [18]. However,
in order to train the wave encoder and the design encoder, gradients w.r.t the
initial latent conditions, parameters $\mathbf{z}(x,t_{\tau})$, and latent
speed of sound $c^{z}(x,t)$ must be computed since these are the outputs of
the encoders. Standard AD methods are not best suited to take gradients
through the integration of a PDE. Instead we have used the adjoint sensitivity
method [19] which is faster and more memory-efficient. Adjoint sensitivity
works by treating backpropagation as the reverse integration of a differential
equation where the solution variable itself is a gradient. Reverse integration
of the gradient w.r.t. each time step of LD to the initial condition allows us
to obtain the information necessary for updating the parameters of the
encoders.

####  III-C7 Model Predictive Control

We use model predictive control (MPC) to demonstrate our model’s ability to
influence the quantity that it predicts. We optimize a sequence of actions
under the objective in Equation (13). In our case, the cost Equation (16) is
summed over the number of actions in the sequence $\tau=1,\dots,n$ where
$\beta$ is a Lagrange Multiplier that penalizes large actions. Random shooting
with $256$ shots and a prediction horizon of $10$ is used to generate and
select an action sequence that is optimal under the given objective. The first
action of the optimal sequence is applied to the RD and this process repeats.
We demonstrate the performance of our model in Section Experiments.

Equations (13)-(15) are adapted from the work of [7] and serve as a
formulation for framing PDE control as an optimal control problem. We minimize
the following functional:

$\min_{a}\;J(\mathbf{w},a)=\min_{a}\sum_{\tau=1}^{T}\ell(\mathbf{w}_{\tau},a_{\tau}),$
(13)

with

$\displaystyle\mathbf{w}_{\tau+1}=$
$\displaystyle\mathbf{w}_{\tau}+\int_{t_{\tau}}^{t_{\tau+1}}\mathcal{N}(\mathbf{w}(x,t),a(t))dt$
(14)

which is a discrete time version of the PDE in continuous time [7]:

$\displaystyle\frac{d\mathbf{w}}{dt}=\mathcal{N}(\mathbf{w}(x,t),a(t)),$ (15)

where $\mathbf{w}$ given by Equation (10) is a vector of fields describing the
state of the RD system. The sequence of actions is represented as the time
dependent function $a(t)$ that changes the state of the design $d(t)$ over a
time interval. In the above equations $\mathcal{N}$ represents a particular
PDE operator. In our case this operator is the PML-transformed first order
classical wave equation system given in Equations (4-9) which acts upon the
fields of $\mathbf{w}$. Minimization in Equation (13) occurs over control
trajectories, $a(t)$ which are obtained through interpreting a design
trajectory $d(t)$ as an interpolation of nonlinear speed of sound functions
$c(\mathbf{x},t)$ in RD. MPC involves optimizing a control trajectory under
objective Equation (13), where $\ell$ is defined Equation (16), then applying
the first action in the sequence $a(t_{\tau})$ to RD. The prediction of
$\hat{\sigma}_{sc}(t)$ occurs by representing $\mathbf{w}$ in a learned lower
dimensional state $\mathbf{z}$. This state is propagated with influence from
predicted latent speed of sound function $c^{z}(x,t)$ to obtain the full
solution $\mathbf{z}(\mathbf{x},t)$

$\ell(\mathbf{w}_{\tau},a_{\tau})=\int_{t_{\tau}}^{t_{\tau+1}}\hat{\sigma}_{sc}(t)dt+\beta\lVert
a(t_{\tau})\rVert_{2}^{2}.$ (16)

##  IV Experiments

In this section we demonstrate the capabilities of our model on the task of
prediction and control of total scattered energy. Several questions are
explored and answered:

  1. 1.

Can our model generalize to long planning horizons in test time when it is
trained only on a short horizon?

  2. 2.

Is it important to dissipate energy in latent space with a trainable PML
layer?

  3. 3.

Will a trained model of the wave dynamics be able to control $\sigma_{sc}(t)$?

####  IV-1 Prediction

The results of predicting $200$ actions ($0.20$s) with a model trained on a
horizon of $20$ actions are shown in Fig 2. We run our latent dynamics with
and without a PML in the LD integration process in Fig 1. It is clear that the
PML model is capable of predicting macro level trends in $\sigma_{sc}(t)$ over
the entire episode of $20000$ integration steps despite not being explicitly
trained to do so; thereby demonstrating qualities of temporal generalization.
Without a PML the model fails to successfully predict for longer time horizons
than it was exposed to in training. A buildup of energy from the source
$f^{z}$ occurs in latent space due to its inability to leave the closed
system. Dissipation of energy in LD allows for stable steady state behavior
which can be extended much longer than its training horizon. The SOTA baseline
we compare our method to is NeuralODE, a general purpose deep neural network
based ordinary differential equations solver. As shown in Fig 2 NeuralODE does
not immediately diverge from the ground truth signal but tracks it for
approximately $75$ actions ($0.075s$). After extending the prediction horizon
of NeuralODE beyond $75$ actions it begins to diverge. We hypothesize that
this divergence is due to NeuralODEs inability to explicitly dissipate latent
energy. Implicitly it must learn some degree of energy dissipation through
it’s trainable dynamics, however, not enough to make long term stable
predictions. Qualitatively, the predictions generated by NeuralODE are rather
simplistic. It generates a signal which generally follows the trend of ground
truth scattering for sub $75$ action horizons but it fails to capture high
frequency oscillations in $\sigma_{sc}$. The predictions generated by our
method appear as thicker lines. This is due to the underlying wave dynamics in
the latent space of our model which are able to simulate low and high
frequency waves.

Figure 2: Demonstration of our models capability to predict $\sigma_{sc}(t)$
over an entire episode from partial information and a sequence of $200$
actions taken in the environment compared against NeuralODE. These curves
represent $20000$ integration timesteps of the dynamics. Prediction
capabilities of our model with energy dissipation (blue), without (orange),
and NeuralODE (red) are compared against true scattered energy (blue).  Figure
3: Evaluation of our model’s prediction error as a function of prediction
horizon compared to NeuralODE. Mean squared error is evaluated at horizons
from $20$ to $200$ in intervals of $10$. Each evaluation is conducted on a
randomly sampled batch of $32$ data points unseen during training.

We investigate how the performance of our model scales with increased horizon
length compared to NeuralODE in Fig 3. Starting at the training horizon of
length $20$ ($0.02s$), the mean squared error loss distribution is visualized
using a line plot which shows the mean and one standard deviation. Prediction
horizon length is incremented by $10$ actions up to the entire episode length
of $200$ ($0.20s$). It can be observed that the mean and variance of loss
measured with the PML model is negligible. The error of NeuralODE’s
predictions is steady at first until approximately $75$ actions beyond which
it rapidly degrades. This comparison shows that our model is significantly
superior to NeuralODE in long term generalization.

####  IV-2 Effect of PML

Energy in latent space $u^{z}_{sc}(x,t)^{2}$ is visualized in Fig 4. Our model
with a trainable PML shown in Fig 4 (left) demonstrates that the scattered
energy field naturally segments into one major zone of activity and three
additional zones which appear very faintly. The PML is learned such that
energy dissipates on either side of each zone and that energy from neighboring
zones does not build up and interfere with one another. This enables steady
energy levels which result in the stable predictions in Fig 2 and Fig 3. Our
model without a PML (right) shows clear signs of energy buildup. The latent
space is unable to be divided into zones and constructive interference causes
rapid accumulation of energy. It is clear this model would not be capable of
making stable predictions without the ability to dissipate excess energy.

####  IV-3 Control

We used our model with a trainable PML to control $\sigma_{sc}(t)$ in the real
environment due to our observation of its generality in the previous result.
An MPC agent uses our model’s prediction of $\hat{\sigma}_{sc}(t)$ produced by
the LD in Fig 1 to inform the cost function in Equation (16). Optimal actions
selected by way of MPC and applied to the real environment. We show a
significant reduction (approximately $30\%$) in mean $\sigma_{sc}(t)$ over the
episodes generated with MPC compared to those with random control in Fig 5.
The results of this experiment show that our model’s predictions are accurate
enough to be used for control. With a more advanced control algorithm the
performance could be further improved and possibly result in an effective
acoustic cloaking device.

Figure 4: Comparison of the scattered energy fields of our model with and
without a latent PML. The field is visualized over $100$ actions ($0.1s$).

##  V Conclusion

We have developed an environment and algorithm which can serve as a basis for
ML research in continuous dynamical systems governed by wave equations. Our
environment is capable of simulating designs which can be manipulated in
continuous time and influence complex PDE systems with physically realizable
mechanisms. The novel algorithm we introduce operates on partial information
from the environment and incorporates domain knowledge from acoustics to
structure a latent representation of system dynamics. Our model is
interpretable through physical properties, therefore, a level of meaning can
be assigned to the evolution of its LD. The interpretable structure of the
latent space we are able to apply a novel and trainable PML layer which
actively dissipates energy. We discover that energy buildup in latent space is
an important attribute and must be dissipated to allow the latent dynamics to
maintain stability in the steady state. To our knowledge, a PML has not been
incorporated into neural networks before. Furthermore, our work is the first
to control of acoustic waves in the time domain with ML.

##  VI Acknowledgments

The authors thank the graduate students at the Computational Intelligence
Control and Information Laboratory at Charles W. Davidson College of
Engineering SJSU, and especially, Megha Jain and Ezgi Kaya for the discussions
and help with proofreading. This work was supported by the Research and
Innovation RSCA Fellowship grant through SJSU.

Figure 5: Comparison of $\sigma_{sc}$ produced in the environment when taking
random actions (blue) versus actions selected with random shooting MPC (green)
with a planning horizon of $10$ actions. Each curve is an average of 6
episodes run with their respective policies.

## References

  * [1] M. Á . Javaloyes, E. Pendás-Recondo, and M. Sánchez, “A general model for wildfire propagation with wind and slope,” _SIAM Journal on Applied Algebra and Geometry_ , vol. 7, no. 2, pp. 414–439, May 2023\. [Online]. Available: https://doi.org/10.1137%2F22m1477866
  * [2] S. H. Schneider and R. E. Dickinson, “Climate modeling,” _Reviews of Geophysics_ , vol. 12, no. 3, pp. 447–493, 1974. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/RG012i003p00447
  * [3] F. Majid, A. M. Deshpande, S. Ramakrishnan, S. Ehrlich, and M. Kumar, “Analysis of epidemic spread dynamics using a pde model and covid-19 data from hamilton county oh usa,” _IFAC-PapersOnLine_ , vol. 54, no. 20, pp. 322–327, 2021, modeling, Estimation and Control Conference MECC 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2405896321022382
  * [4] T. Shah, L. Zhuo, P. Lai, A. De La Rosa-Moreno, F. Amirkulova, and P. Gerstoft, “Reinforcement learning applied to metamaterial designa),” _The Journal of the Acoustical Society of America_ , vol. 150, no. 1, pp. 321–338, 07 2021. [Online]. Available: https://doi.org/10.1121/10.0005545
  * [5] S. Cominelli, D. E. Quadrelli, C. Sinigaglia, and F. Braghin, “Design of arbitrarily shaped acoustic cloaks through partial differential equation-constrained optimization satisfying sonic-metamaterial design requirements,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 478, no. 2257, Jan. 2022. [Online]. Available: https://doi.org/10.1098%2Frspa.2021.0750
  * [6] K. Bieker, S. Peitz, S. L. Brunton, J. N. Kutz, and M. Dellnitz, “Deep model predictive flow control with limited sensor data and online learning,” _Theoretical and Computational Fluid Dynamics_ , vol. 34, no. 4, pp. 577–591, Mar. 2020. [Online]. Available: https://doi.org/10.1007/s00162-020-00520-4
  * [7] S. Werner and S. Peitz, “Learning a model is paramount for sample efficiency in reinforcement learning control of pdes,” 2023.
  * [8] Y. Pan, A. massoud Farahmand, M. White, S. Nabi, P. Grover, and D. Nikovski, “Reinforcement learning with function-valued action spaces for partial differential equation control,” 2018.
  * [9] A. N. Norris, “Acoustic cloaking theory,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 464, no. 2097, pp. 2411–2434, Apr. 2008. [Online]. Available: https://doi.org/10.1098%2Frspa.2008.0076
  * [10] T. Wu, T. Maruyama, and J. Leskovec, “Learning to accelerate partial differential equations via latent global evolution,” 2022.
  * [11] A. Rahman, Y. Chang, and J. Rubin, “Interpretable additive recurrent neural networks for multivariate clinical time series,” 2021.
  * [12] W. Zaremba, I. Sutskever, and O. Vinyals, “Recurrent neural network regularization,” 2015.
  * [13] A. N.Norris, “Acoustic cloaking theory,” _prsa_ , vol. 464, pp. 2411–2434, 2008.
  * [14] F. A. Amirkulova and A. N. Norris, “The gradient of total multiple scattering cross-section and its application to acoustic cloaking,” _Journal of Theoretical and Computational Acoustics_ , p. 1950016, Mar. 2020.
  * [15] D. Torrent, F. Håkansson, A. Cervera, and J. Sánchez-Dehesa, “Homogenization of two-dimensional clusters of rigid rods in air,” _Phys. Rev. Lett._ , vol. 96, p. 204302, 2006. [Online]. Available: https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.96.204302
  * [16] S. G. Johnson, “Notes on perfectly matched layers (pmls),” 2021.
  * [17] D. K. Prigge, “Absorbing boundary conditions and numerical methods for the linearized water wave equation in 1 and 2 dimensions,” Ph.D. dissertation, University of Michigan, 2016.
  * [18] M. Innes, E. Saba, K. Fischer, D. Gandhi, M. C. Rudilosso, N. M. Joy, T. Karmali, A. Pal, and V. Shah, “Fashionable modelling with flux,” _CoRR_ , vol. abs/1811.01457, 2018. [Online]. Available: https://arxiv.org/abs/1811.01457
  * [19] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” 2019.

##  VII Appendix - PML Derivation

Step 1. Given the split field form of the acoustic wave equation with pressure
field $u:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$, velocity component in the
$x$ direction $v_{x}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$, and velocity
component in the $y$ direction
$v_{y}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$, and speed of sound
$c=\sqrt{ab}$.

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=b\frac{\partial
v_{x}}{\partial x}+b\frac{\partial v_{y}}{\partial y}$ (17)
$\displaystyle\frac{\partial v_{x}}{\partial t}$
$\displaystyle=a\frac{\partial u}{\partial x}$ (18)
$\displaystyle\frac{\partial v_{y}}{\partial t}$
$\displaystyle=a\frac{\partial u}{\partial y}$ (19)

Step 2. The transformations are described by complex coordinate stretching in
the frequency domain. The terms $\sigma_{x}$ and $\sigma_{y}$ are scalar
functions which describe where the PML region is active and $\omega$ is
frequency.

$\displaystyle\frac{\partial}{\partial
x}\rightarrow\frac{1}{1+i\frac{\sigma_{x}}{\omega}}\frac{\partial}{\partial
x}$ (20) $\displaystyle\frac{\partial}{\partial
y}\rightarrow\frac{1}{1+i\frac{\sigma_{y}}{\omega}}\frac{\partial}{\partial
y}$ (21)

Step 3. Convert (17) to frequency domain and apply (20) and then (21).

$-i\omega u=b\frac{\partial v_{x}}{\partial x}+b\frac{\partial v_{y}}{\partial
y}+b\frac{\partial v_{x}}{\partial
x}\frac{i\sigma_{y}}{\omega}+b\frac{\partial v_{y}}{\partial
y}\frac{i\sigma_{x}}{\omega}\\\
-u(\sigma_{x}+\sigma_{y})-i\frac{u\sigma_{x}\sigma_{y}}{\omega}$ (22)

Step 4. Every $\frac{i}{\omega}$ in (22) is converted to a field in the
frequency domain.

$\displaystyle-i\omega\Psi_{y}$ $\displaystyle=b\sigma_{y}\frac{\partial
v_{x}}{\partial x}$ $\displaystyle-i\omega\Psi_{x}$
$\displaystyle=b\sigma_{x}\frac{\partial v_{y}}{\partial y}$
$\displaystyle-i\omega\Omega$ $\displaystyle=u\sigma_{x}\sigma_{y}$
$\displaystyle-i\omega u$ $\displaystyle=b\frac{\partial v_{x}}{\partial
x}+b\frac{\partial v_{y}}{\partial
y}-u(\sigma_{x}+\sigma_{y})+\Psi_{x}+\Psi_{y}-\gamma$

Step 5. Convert the system of equations back into the time domain.

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=b\frac{\partial
v_{x}}{\partial x}+b\frac{\partial v_{y}}{\partial
y}-u(\sigma_{x}+\sigma_{y})+\Psi_{x}+\Psi_{y}-\gamma$ (23)
$\displaystyle\frac{\partial\Psi_{x}}{\partial t}$
$\displaystyle=b\sigma_{x}\frac{\partial v_{y}}{\partial y}$ (24)
$\displaystyle\frac{\partial\Psi_{y}}{\partial t}$
$\displaystyle=b\sigma_{y}\frac{\partial v_{x}}{\partial x}$ (25)
$\displaystyle\frac{\partial\gamma}{\partial t}$
$\displaystyle=u\sigma_{x}\sigma_{y}$ (26)

Step 6. Convert (18) to frequency domain and apply (20).

$-i\omega v_{x}+v_{x}\sigma_{x}=a\frac{\partial u}{\partial x}$ (27)

Step 7. Convert (19) to frequency domain and apply (20).

$-i\omega v_{y}+v_{y}\sigma_{y}=a\frac{\partial u}{\partial y}$ (28)

Step 8. Convert (27) and (28) to time domain.

$\displaystyle\frac{\partial v_{x}}{\partial t}$
$\displaystyle=a\frac{\partial u}{\partial x}-v_{x}\sigma_{x}$ (29)
$\displaystyle\frac{\partial v_{y}}{\partial t}$
$\displaystyle=a\frac{\partial u}{\partial y}-v_{y}\sigma_{y}$ (30)

The final system of equations used to simulate acoustic waves in 2D with a pml
is given by the equations in Step 5 and Step 8.