# AdverSAR: Adversarial Search and Rescue via Multi-Agent Reinforcement
Learning

Aowabin Rahman1, Arnab Bhattacharya1, Thiagarajan Ramachandran1, Sayak
Mukherjee1, Himanshu Sharma1,
Ted Fujimoto2, Samrat Chatterjee3
Optimization and Control Group1, Data Analytics Group2, Data Sciences and
Machine Intelligence Group3
Pacific Northwest National Laboratory
Richland, USA

###### Abstract

Search and Rescue (SAR) missions in remote environments often employ
autonomous multi-robot systems that learn, plan, and execute a combination of
local single-robot control actions, group primitives, and global mission-
oriented coordination and collaboration. Often, SAR coordination strategies
are manually designed by human experts who can remotely control the multi-
robot system and enable semi-autonomous operations. However, in remote
environments where connectivity is limited and human intervention is often not
possible, decentralized collaboration strategies are needed for fully-
autonomous operations. Nevertheless, decentralized coordination may be
ineffective in adversarial environments due to sensor noise, actuation faults,
or manipulation of inter-agent communication data. In this paper, we propose
an algorithmic approach based on adversarial multi-agent reinforcement
learning (MARL) that allows robots to efficiently coordinate their strategies
in the presence of adversarial inter-agent communications. In our setup, the
objective of the multi-robot team is to discover targets strategically in an
obstacle-strewn geographical area by minimizing the average time needed to
find the targets. It is assumed that the robots have no prior knowledge of the
target locations, and they can interact with only a subset of neighboring
robots at any time. Based on the centralized training with decentralized
execution (CTDE) paradigm in MARL, we utilize a hierarchical meta-learning
framework to learn dynamic team-coordination modalities and discover emergent
team behavior under complex cooperative-competitive scenarios. The
effectiveness of our approach is demonstrated on a collection of prototype
grid-world environments with different specifications of benign and
adversarial agents, target locations, and agent rewards.

###### Index Terms:

Search and Rescue, Multi-agent Reinforcement Learning, Adversarial
Reinforcement Learning, Critical Infrastructure Security.

##  I Introduction

Multi-agent autonomous teams are becoming increasingly ubiquitous in a variety
of Search and Rescue (SAR) missions that involve locating and rescuing targets
in emergency situations. For example, military and first responders often
deploy multi-robot teams to locate missing targets or survivors during
disaster scenarios. In remote locations where communication with human
operators may not be possible, a SAR team needs to search large areas of
secluded geographical terrain in a fully autonomous and decentralized manner
[1]. This can lead to a significant increase in time, effort, and operational
costs when missing-target locations are not known a priori with any degree of
certainty. In such settings, inter-agent communications play a key role in
enabling collaborative and coordinated decision-making [2]. This is because
information shared by other robots (e.g. location of identified targets,
terrain mapping) can help a robot reason about the current state of the world
and narrow down it’s search area rather than exploring independently using
local information. Team performance in cooperative settings can significantly
improve if agents systematically coordinate their exploration.

However, multi-agent teams not only need to learn on how to cooperate, but
also, on how to communicate effectively for cooperation. Even minor
manipulations of messages shared within a team can significantly reduce the
benefits of coordination [3, 4]. This destabilizing effect of perturbing
communication data is further amplified in problems with sparse rewards [5].
For example, in a SAR mission, the agents receive rewards only when some or
all of the missing targets have been discovered and rescued (i.e., when a
mission is completed successfully). Stealthy adversaries can strategically
broadcast falsified sensor readings that prevents a SAR team from reaching or
locating all of the targets despite active coordination between the agents
[1].

In this paper, we propose an adversarial multi-agent reinforcement learning
(MARL) approach to learn efficient coordination and collaboration strategies
to accomplish SAR mission objectives in the presence of adversarial
communications. Our training approach is based on the well-known centralized
training with decentralized execution (CTDE) paradigm [6] in MARL that uses
full system information of all agents to generate rich features for learning
effective decentralized agent-level policies. We have two major contributions
in this paper. First, we developed a a general reward structure for MARL that
allows effective adversarial modeling and inference in sparse-reward settings.
Second, we developed an adversarial multi-agent reinforcement learning
algorithm that allows both cooperative and adversarial agents to be trained
simultaneously. The effectiveness of our algorithmic approach is demonstrated
for different environment and team specifications by reducing the overall
target-exploration time of cooperative agents in adversarial presence.

##  II Problem Description

Our problem setup comprises of a team of cooperative agents tasked with
covering a search area efficiently to locate missing assets (or targets). The
locations of the missing assets are not known a priori. Therefore, agents need
to coordinate their respective search areas to avoid visiting the same
locations multiple times. A mission is successful when all of the missing
targets have been located within a pre-specified mission completion time. We
make three assumptions in our model. First, it is assumed that all agents have
homogeneous sensing capabilities. Second, agents cannot explicitly communicate
with each other during mission execution. Finally, it is assumed that target
locations do not change over the search time of the SAR team. The problem is
formulated as a Decentralized Partially Observable Markov Decision Process
(Dec-POMDP) model [7] defined using a tuple $M=(S,A,O,P,R,Z,\gamma)$; here,
$S$ denotes the finite set of states for all agents, $A$ is the joint set of
actions, $O$ is the bounded set of observations,
$P(s^{\prime}|s,a)=P(s^{\prime}=s_{t+1}|s=s_{t},a=a_{t})$ is the transition
function, $R$ is the reward function, $Z(o_{t}|s_{t},a_{t})$ is the
observation function and $\gamma$ is the discount factor. The goal of a
general Dec-POMDP is to solve for a set of decentralized policies,
$\pi=\\{\pi^{1},\pi^{2},...\pi^{N}\\}$. The value function for a given policy
$\pi$ starting from an initial state $s_{0}$ is expressed as [8]:

$V^{\pi}=\mathbb{E}_{\pi}\big{(}\sum_{t=0}^{\infty}R(s_{t},a_{t})|s_{0}\big{)}$
(1)

The solution to the Dec-POMDP is a policy $\pi^{*}$ that maximizes the value
function $\pi^{*}=\operatorname*{arg\,max}_{\pi}V^{\pi}$. Next, we propose an
approach to model adversarial inference in a SAR context.

##  III State-of-the-art Model for SAR Exploration

###  III-A Summary of Current Approach

For the naive case, i.e. SAR without the presence of adversaries, we leverage
a two-tier MARL algorithm proposed in [9]. We would like to note that our key
contributions in this paper concern modeling of adversarial presence and
modifications to the algorithm by Iqbal and Sha [9] to account for adversarial
presence – which are detailed in section IV. Nonetheless, we provide a brief
overview of the the two-tier MARL approach in [9] as follows:

  * •

An upper-tier policy selection model: A meta-policy model selects a single
policy from a set of policies that all agents will adopt. Each policy from the
set corresponds to a different collaboration strategy for exploration between
the cooperative agents. We will discuss later on how each collaboration
strategy is mathematically modeled as an intrinsic reward. The parameters of
the upper-tier model, $\phi$, are learned via policy-gradient updates [9].

  * •

A lower-tier soft actor-critic (SAC) model: The lower-tier model optimizes
individual agent policies (corresponding to each collaborative exploration
strategy/novelty function). The model uses a central critic-network shared
across all agents, and a policy model specific to each agent. For improved
sample efficiency, the model is trained using off-policy learning methods,
i.e. all policies and value functions are learned using all data collected,
regardless of which policy and exploration strategy were used to generate a
given data sample.

###  III-B Intrinsic Rewards

As mentioned in Section III-A, each modality of team collaboration can be
modeled as an intrinsic reward function. Since the SAR mission is a sparse-
reward problem, i.e. the number of missing targets is small relative to the
size of the domain, the intrinsic rewards serve a mechanism for training the
agents to explore the domain efficiently. First, we define a novelty function
$f_{i}$, which denotes how novel a given agent $i$ considers its own
observations at a given location, based on past visits to that location. In a
discrete environment domain (such as Gridworld), this can simply be the
inverse of the state visits [9]. An intrinsic reward for agent $i$, denoted by
$g_{i}$, is defined as a function that quantifies how novel all other agents
consider the observations of agent $i$, and is expressed as

$g:[f_{1}(\mathbf{o_{i}}),f_{2}(\mathbf{o_{i}}),...f_{n}(\mathbf{o_{i}})]\mapsto
r_{intr}.$ (2)

Here, $r_{intr}\in\mathbb{R}$ is a scalar value. For this paper, we consider
three different exploration modalities, defined using three different
intrinsic reward function $g_{i}$:

  * •

Minimum strategy rewards an agent when it explores a location no other agent
has explored. The intrinsic reward for this case can be expressed as

$g_{i}=\min_{j\in\\{1...N\\}}f_{j}(o_{i}).$

  * •

Covering strategy rewards an agent when it explores a location that is more
novel than the average agent and is defined as

$g_{i}=f_{i}(o_{i})\mathbbm{1}[f_{i}(o_{i})>\mu(o_{i})],$

where $\mu_{i}=\sum_{j}f_{j}(o_{i})/N$.

  * •

Burrowing strategy is the reverse of the covering strategy, as it rewards
agents for exploring a location less novel than the average agent. This
strategy leads to agents identifying the dead-ends quickly, which could lead
to a reduced overall exploration time, and is expressed as

$g_{i}=f_{i}(o_{i})\mathbbm{1}[f_{i}(o_{i})<\mu(o_{i})],$

where $\mathbbm{1}(\cdot)$ denotes an indicator function.

Additional details on how to construct intrinsic reward functions is provided
in [9].

###  III-C Centralized Training with Decentralized Execution

Centralized Training with Decentralized Execution (CTDE) is an MARL training
paradigm in which the agents share information during training (e.g.
observations of other agents), but act on their own local observations during
execution/evaluation [10, 6]. CTDE is useful in avoiding the non-stationarity
issues that often arise in training multi-agent systems. Since a centralized
critic in an actor-critic algorithm has access is used to observations of all
agents, the agents are less likely to encounter dynamic changes in the
environment – which can lead to improved stability during training [11]. In
the context of our problem, a global state and an intrinsic reward function
are used to train the policies of each agent, both of which require each agent
having access to information of other agents. However, during execution, each
agent acts based on its own policy, which only requires access to local
observation specific to that agent.

##  IV Proposed Algorithmic Approach

###  IV-A Adversarial Modeling

Here, we describe the SAR problem in the adversarial context. We introduce
$N_{adv}$ adversarial agents (where $N_{adv}\in N$) in our domain, in addition
to the cooperative agents already in the domain ($N_{c}\in N$). The overall
goal of the adversary is to increase the time taken for the cooperative agents
to find the missing assets. We make the following considerations in defining
the adversary’s reward structure.

  * •

The adversarial team has full knowledge of the location of the missing assets,
and at each time-step, $t$, incurs a reward that discourages the cooperative
agents to complete their mission. In our setup, we propose the following
reward function for each adversarial agent:

$r_{adv}^{ext}=\alpha\sum_{i=1}^{N_{c}}\sum_{m=1}^{M_{nf}}|x_{i}-x_{m}|+|y_{i}-y_{m}|\vspace{-0.5mm}$
(3)

In Equation 3, $(x_{i},y_{i})$ denotes the position of agent $i$ at a given
time $t$, $(x_{m},y_{m})$ denotes the position of target $m$, $M_{nf}$ denotes
the number of missing assets not discovered by the cooperative agents at time-
step $t$, and $\alpha$ is a normalizing factor defined as
$\alpha=\frac{K}{N_{c}(L+W)}$. Note that $N_{c}$ is the number of cooperative
agents, $L$ and $W$ are the dimensions of a discrete environment, and $K$ is a
hyper-parameter such that $0<K\leq 1$. In this paper, we set $K=1$.

  * •

The adversarial agents are rewarded for delaying the cooperative agents from
completing the mission. Unlike the cooperative agents, the adversarial agents
are not rewarded for locating the missing assets.

  * •

If, at a given location in the search environment, an adversarial agent
encounters a missing asset, it will spoof and falsify the location of the
missing asset.

Table I summarizes the reward structure for the cooperative and the
adversarial agents, which will be referred to as the “baseline” reward
structure. Next, in section IV-B, we describe a modified reward structure to
enhance coordination under adversarial influence.

TABLE I: Reward structure of cooperative agents vs. adversarial agents | Cooperative agents | Adversarial Agents
---|---|---
External Rewards: |  |
• Time penalty at each timestep: | $-0.1$ | $+0.1$
• Locating asset (not located already by other assets): | $+10$ | $0$
• Completing task (i.e. locating all assets) | $+10$ | $0$
• Not completing task | -10 | 0
• Adversarial reward at each timestep, $t$ | N/A | Eq. 3

###  IV-B Reward Structure for Optimal Coverage

We aim to extend the MARL algorithm for a general optimal coverage problem,
and propose a “modified” reward structure for both the cooperative and the
adversarial teams. One key limitation of the existing model is that the
primary extrinsic reward that the agents receive are dependent on target
locations. This makes the trained two-tier MARL model difficult to generalize
during evaluation for a new environment map and target locations. As such, we
define a reward structure to encourage cooperative agents to explore the
environment, irrespective of the location of the targets. We propose that at
each time-step $t$, the cooperative agents are rewarded for exploring states
that are novel, i.e. states that have not been visited by the cooperative
team. Therefore, during training, we assign a secondary reward $r_{sec}$
signal defined as:

$r_{sec,coop}=\sum_{i=1}^{N_{c}}\mathbbm{1}[v_{N_{C}}(x_{i},y_{i})=1].$ (4)

In Equation 4, $v_{N_{C}}$ denotes the total number of visits by the
cooperative agents at a given location ($x_{i}$, $y_{i}$) in the grid, i.e.,
$v_{N_{C}}=\sum_{i}^{N_{c}}v_{i}$ . Thus, Equation 4 indicates that the
cooperative team is rewarded when a given cooperative agent moves to a state
not visited by the team previously. For the adversarial agents, the secondary
reward signal is defined as

$r_{sec,adv}=\sum_{i=1}^{N_{c}}\mathbbm{1}[v_{N_{C}}(x_{i},y_{i})>v_{thresh})].$
(5)

Here, $v_{thresh}$ denotes a threshold value for the number of visits, beyond
which a state visit is considered redundant (we set $v_{thresh}=1$). The
intuition behind the adversarial reward signal is that adversarial agents
should encourage cooperative agents to visit redundant states that have been
previously visited. The cumulative reward function is a weighted sum of the
cooperative and adversarial team signals defined as

$r=r_{sec}+\beta r_{intr},$ (6)

where $\beta$ is time-varying parameter defined as

$\beta$ = {$0.1$, $t\leq 0.4T_{max}$,$0.1\exp(-kt)$, otherwise.

Here, $T_{max}$ denotes the maximum number of time-steps per training episode.
Note that the proposed reward structure for a general optimal-coverage problem
does not require the terminal state to be dependent on location of targets.
Thus, during training, we train the cooperative and the adversarial agents for
a given map with no targets. During inference, however, we revert to the SAR
setting requiring cooperative agents to locate their targets.

TABLE II: Reward structure for optimal coverage (as detailed in section IV-B.  | Cooperative agents | Adversarial Agents
---|---|---
External Rewards: |  |
• Time penalty at each timestep: | $0$ | $0$
• Locating asset (not located already by other assets): | N/A | N/A
• Cooperative secondary reward at each timestep | Eq. 4 | N/A
• Adversarial reward at each timestep, $t$ | N/A | Eq. 5

###  IV-C Training Algorithm

We leverage an adversarial MARL training procedure to train the cooperative
agents in presence of adversarial agents. The adversarial actions inject
continuous perturbations to the performance of the cooperative agents and,
therefore, once the cooperative agents are trained in such scenarios, we can
guarantee sufficient robustness to malicious influences. We simultaneously
train the adversary’s action to minimize the reward obtained by the
cooperative agents, thereby, making the cooperative agents aware of an
intelligently trained adversary. To implement this idea, we follow [12] to
iteratively train the cooperative and adversarial agents. That is, we first
train the cooperative agents keeping the adversary’s policy fixed, followed by
the cooperative agent’s policy being fixed and adversarial policy being
updated (as preesented in algorithm 1). We hypothesize that training the
adversarial agents alternatively with the cooperative agents, as opposed to
training both cooperative and adversarial agents together, will improve
adversarial inference and effectiveness of the cooperative team’s policies.

Algorithm 1 Adversarial training for coordinated multi-agent exploration.

Initialize environment with $N$ agents: $N_{c}$ cooperative agents and
$N_{adv}$ adversarial agents

Initialize replay buffer $D$ and $t_{update}\leftarrow 0$

nitialize choice for random asset locations

for t = 1…total number of timesteps do

if episode is completed then

for j = 1…number of policies do

UpdateSelector(R, h)

end for

Sample policy head: $h\sim\prod$

Set $t_{ep}\leftarrow 0$ and reward $R\leftarrow 0$

end if

Select actions for each agent $n$, $a_{i}\sim\pi_{n}^{h}(\cdot|{o_{n}})$,
$n\in N$

Send actions to environment to fetch
$s^{\prime},\mathbf{o^{\prime}},r_{coop},r_{adv}$

$R\leftarrow R+\gamma^{t_{ep}}r$

Store tuple
($\mathbf{s},\mathbf{o},\mathbf{a},\mathbf{s^{\prime}},\mathbf{o^{\prime}},r_{coop}$)
in $D_{1}$

Store tuple
($\mathbf{s},\mathbf{o},\mathbf{a},\mathbf{s^{\prime}},\mathbf{o^{\prime}},r_{adv}$)
in $D_{2}$

Increment $t_{update}+=1$ and $t_{ep}+=1$

if $t_{update}$ == steps per update then

for j = 1..$N_{iter,c}$ do

Sample minibatch, $B\sim D_{1}$

Update critics using all coop agents, $n\in N_{c}$

Update policies for all coop agents, $n\in N_{c}$

Update targets for coop agents: $\bar{\phi_{c}},\bar{\psi_{c}}$

end for

for j = 1..$N_{iter,adv}$ do

Sample minibatch, $B\sim D_{2}$

Update critics using all adv. agents, $n\in N_{adv}$

Update policies for all adv. agents, $n\in N_{adv}$

Update targets for adv. agents: $\bar{\phi_{adv}},\bar{\psi_{adv}}$

end for

Reinitialize $t_{update}=0$

end if

end for

##  V Numerical Results

###  V-A Experimental Setup

We consider a 2D Grid-world domain (see Figure 1) of size $20\times 20$ as our
simulation environment. Each cell within the domain represents the local
position ($x_{i},y_{i}$) of an agent $i$. At each time $t$, an agent can
choose an action $a_{i}$ from a maximum of four actions; each action
corresponding to taking a single step among four directions (going left,
right, up or down).

Figure 1: Schematic of Grid-world environment with two cooperative agents

In our experiments, we aim to investigate: (a) how adversarial interference
can derail performance of cooperative agents in a SAR mission; and (b) to what
extent can the proposed training algorithm mitigate adversary’s actions. We
perform four case studies, as detailed in Table III. We compare cases I vs. II
to investigate objective (a) and cases III vs. IV to investigate objective
(b). We evaluated the performance of the cooperative agents using the flow-
time metric, which is defined as the number of time-steps needed for the
cooperative team to locate the targets. For each case, we compute the mean
flow-time by taking the average across 12 parallel instantiations of the
environment. To investigate how well our models generalize, we evaluate the
trained models on a different Grid-world map different to the one used during
training, for two different target locations (see Figure 2).

Figure 2: Grid-world maps used for training and inference: fig. (a) shows the map used for training (case III), whereas fig (b) and (c) show the inference maps with different sets of locations for the targets. Note that for case I, II and IV, we use a map identical to (a) for training, except that the target locations are not included. TABLE III: Case studies performed  Case Study | Description |  | Coop. agents
---
| Adversarial agents
---
| Reward
---
Structure
I |  | Algorithm with reward structure for optimal coverage.
---
SAC trained with two cooperative agents.
2 during training + inference | None | Modified
II |  | Training identical to case I. During inference, the
---
third cooperative agent is swapped for an adversarial agent.
| 3 during training
---
2 during inference
1 | Modified
III |  | Algorithm with reward structure specific to agents locating
---
missing assets. SAC trained with two cooperative and one
adversarial agent
2 during training + inference | 1 | Baseline
IV |  | Algorithm with reward structure for optimal coverage.
---
SAC trained with two cooperative and one adversarial agent.
2 during training + inference | 1 | Modifed

###  V-B Adversarial Inference in a Cooperative SAR Setting

To illustrate the impact of adversarial interference, we consider the case
where we first train the SAC model with three cooperative agents, but during
inference, we replace one of the cooperative agents with an adversarial agent
(case II). We compare this case against one where we train and infer with two
cooperative agents (i.e. the number of cooperative agents remain identical
across both cases). As reported in Table IV, the mean flow-times for case I is
8500, while for case II, the cooperative agent was unable to locate the target
within 18,000 time-steps (which we set as the maximum number of time-steps
during inference). Figure 3 illustrates the comparison between Cases I and II
during inference, for one of the 12 parallel instantiations. For both cases I
and II, the agents adopt a “burrowing” strategy. While counter-intuitive, this
strategy allows the agent to quickly identify dead ends in the map and then
diverge to cover the domain. In Case I, the cooperative agents (in yellow and
green) progressively cover more novel states to locate both of the missing
assets, as seen from figure 3(a)-(c). However, the cooperative agents are
thwarted, as they do not cover the sub-region where the missing assets are
located. This is primarily because the right side of the domain was covered by
a third cooperative agent during training, which, during inference, is swapped
with an adversarial agent. Another reason for this could be that the adversary
(blue agent) in Case II adopts a policy where it follows a trajectory close to
one of the cooperative agents (in yellow). Since the observation of each agent
includes whether another agent is present within three grid spaces, the
adversarial presence in case II can “spoof” the cooperative agent and prevent
it from covering novel states. Thus, we observe that the presence of an
adversary can significantly hinder a cooperative MARL task such as SAR.
Therefore, it is critical to incorporate the presence of an adversary while
training a multi-agent setup.

Figure 3: Illustration of agent trajectories during inference for a sample instantiation: (a), (b) and (c) correspond to preliminary, intermediate and terminal stages for cases I and II. For case II, the RL algorithm had been trained with three cooperative agents; but during inference, one of the coooperative agents was swapped for an adversarial agent (in blue). We observe that the adversarial agent prevents the cooperative agents from reaching the target locations.  TABLE IV: Mean flow-time during inference for each of the four cases Case | Mean flow-time (timesteps)
---|---
Map A | Map B
I | 1625 | 1302
II | >18000 | >18000
III | >18000 | >18000
IV | 290 | 789

###  V-C Adversarial Threat Mitigation using Proposed Algorithm

We now investigate how effective the reward structures for optimal coverage
(for the cooperative and adversarial agents) and the algorithm proposed in
section III are in mitigating adversarial interference (case III), which we
compare against a case in which we use the reward structure for cooperative
structure as defined in [9], and the adversarial reward as defined in Section
3. We observe that in Case III, the cooperative agents head directly for the
locations where the missing assets were located during training, since the
episode rewards in this case were primarily dependent on finding the missing
assets. As the location of the missing assets have changed from training to
inference, this lack of generalizability leads to the cooperative agents
failing to explore the domain efficiently and locate one of the missing assets
within the allotted 18000 time-steps. In Case IV, where the rewards are not
goal-conditioned (i.e. dependent on the location of the targets) and the
impact of the adversary has been incorporated, the mean flow-time is 289.5.
This supports our hypothesis that the proposed reward structure allows agents
to generalize novel states better.

Figure 4: Illustration of agent trajectories during inference for a sample
instantiation: (a), (b) and (c) correspond to preliminary, intermediate and
terminal stages for cases III and IV. In case III, the reward structure and
the terminal state during training are dependent on agents reaching the
targets. During inference when the map and the target locations are altered,
the cooperative agents failed to reach the targets within 18,000 time-steps.
In case IV, with the modified reward structure, the agents explored the space
more efficiently during inference and able to locate the targets.

##  VI Conclusions and Future Work

In this paper, we have shown that adversarial interference can significantly
hinder the performance of a cooperative team in a SAR mission. We developed a
reward structure that focuses on optimal coverage and not conditioned on
target location, and an adversaril MARL training algorithm to learn
cooperative team behavior that mitigates adversarial impact and improve
generalizability in different SAR scenarios. As next steps, we seek to
integrate our algorithm and demonstrate results in a higher fidelity 3D
environment such as Microsoft AirSim. Moreover, we had assumed that the agents
execute their decentralized policies based on their local observations, i.e.
no communication between the agents were considered during execution. One
possible research direction would be to investigate how adversarial
interference could affect a cooperative team, where each agent is able to
communicate with other agents within a given communication radius. The code is
open-source and hosted in github.com/aowabinr/SAR_hst.

## Acknowledgment

The research described here is part of the Mathematics for Artificial
Reasoning in Science Initiative at Pacific Northwest National Laboratory
(PNNL). It was conducted under the Laboratory Directed Research and
Development Program at PNNL, a multiprogram national laboratory operated by
Battelle for the U.S. Department of Energy under contract DE-AC05-76RL01830.

## References

  * [1] J. P. Queralta, J. Taipalmaa, B. C. Pullinen, V. K. Sarker, T. N. Gia, H. Tenhunen, M. Gabbouj, J. Raitoharju, and T. Westerlund, “Collaborative multi-robot search and rescue: Planning, coordination, perception, and active vision,” _IEEE Access_ , vol. 8, pp. 191 617–191 643, 2020.
  * [2] L. Zhou and P. Tokekar, “Multi-robot coordination and planning in uncertain and adversarial environments,” _Current Robotics Reports_ , vol. 2, no. 2, pp. 147–157, 2021.
  * [3] J. Tu, T. Wang, J. Wang, S. Manivasagam, M. Ren, and R. Urtasun, “Adversarial attacks on multi-agent communication,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 7768–7777.
  * [4] A. Behjat, H. Manjunatha, P. K. Kumar, A. Jani, L. Collins, P. Ghassemi, J. Distefano, D. Doermann, K. Dantu, E. Esfahani _et al._ , “Learning robot swarm tactics over complex adversarial environments,” in _2021 International Symposium on Multi-Robot and Multi-Agent Systems (MRS)_. IEEE, 2021, pp. 83–91.
  * [5] D. Wang, B. Ding, and D. Feng, “Meta reinforcement learning with generative adversarial reward from expert knowledge,” in _2020 IEEE 3rd International Conference on Information Systems and Computer Aided Education (ICISCAE)_. IEEE, 2020, pp. 1–7.
  * [6] R. Lowe, Y. I. Wu, A. Tamar, J. Harb, O. Pieter Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [7] D. Silver, “Lecture 2: Markov decision processes,” _UCL. Retrieved from https://www.davidsilver.uk/wp-content/uploads/2020/03/MDP.pdf_ , 2015.
  * [8] S. Omidshafiei, A.-A. Agha-Mohammadi, C. Amato, and J. P. How, “Decentralized control of partially observable markov decision processes using belief space macro-actions,” in _2015 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2015, pp. 5962–5969.
  * [9] S. Iqbal and F. Sha, “Coordinated exploration via intrinsic rewards for multi-agent reinforcement learning,” _arXiv preprint arXiv:1905.12127_ , 2019\.
  * [10] X. Lyu, Y. Xiao, B. Daley, and C. Amato, “Contrasting centralized and decentralized critics in multi-agent reinforcement learning,” _arXiv preprint arXiv:2102.04402_ , 2021.
  * [11] G. Papoudakis, F. Christianos, A. Rahman, and S. V. Albrecht, “Dealing with non-stationarity in multi-agent deep reinforcement learning,” _arXiv preprint arXiv:1906.04737_ , 2019.
  * [12] L. Pinto, J. Davidson, R. Sukthankar, and A. Gupta, “Robust adversarial reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 2817–2826.


# Extending turbulence model uncertainty quantification using machine learning

Marcel Matha and Christian Morsbach
Institute of Propulsion Technology, German Aerospace Center (DLR)
Linder Höhe, 51147 Cologne, Germany
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>

###### Abstract

In order to achieve a more virtual design and certification process of jet
engines in aviation industry, the uncertainty bounds for computational fluid
dynamics have to be known. This work shows the application of a machine
learning methodology to quantify the epistemic uncertainties of turbulence
models. The underlying method in order to estimate the uncertainty bounds is
based on an eigenspace perturbation of the Reynolds stress tensor in
combination with random forests.

##  1 Introduction

As a compromise between computational time and accuracy, Reynolds-averaged
Navier-Stokes (RANS) simulation is the workhorse in the industrial design
process of turbomachinery. The derivation of the RANS equations reveals an
unclosed term, called the Reynolds stress tensor. This tensor has to be
approximated in Computational-Fluid-Dynamics (CFD) simulations by applying
turbulence models. The prediction accuracy of the simulation is highly
dependent on this kind of models, which mostly accomplished transport
equations, mimicking main turbulence flow physics. Although RANS-based models,
such as linear eddy viscosity models, are widely used for complex engineering
flows, they suffer from the inability to replicate fundamental turbulent
processes. Turbulence models are one of the main limitations in striving for
reliable, environmental-friendly designs, due to general simplifying
assumptions during formulation of closure models. These simplifications are
the result of data observation and physical intuition, leading to a
significant degree of epistemic uncertainty.

In recent years, the interest in quantifying these uncertainties, leading to
more reliable simulative results, has grown. The group of Iaccarino proposed
an eigenspace perturbation framework, which is based on the inability of
common linear eddy viscosity models to deal with Reynolds stress tensor
anisotropy [4] [6]. The emergence of machine learning strategies guided the
path towards data driven approaches also for the turbulence modelling
community [3]. Heyse et al. recently enhanced the uncertainty estimation based
on the eigenspace perturbation approach by adding a data-driven method [5].

In order to obtain more sophisticated and trustworthy simulative results, we
investigate data-driven enhancements of the Reynolds stress perturbation
approach with DLR’s CFD solver suite TRACE. TRACE is being developed by the
Institute of Propulsion Technology with focus on turbomachinery and offers a
parallelized, multi-block CFD solver for compressible RANS equations.111TRACE
User Guide, trace-portal.de/userguide

##  2 Eigenspace perturbation framework

###  2.1 Data-free approach

RANS turbulence models are utilized in order to determine the Reynolds stress
tensor $\tau_{ij}=\overline{u_{i}^{\prime}u_{j}^{\prime}}$ in terms of mean
flow quantities, whereas $u^{\prime}$ is the fluctuating part of the velocity
($u=\overline{u}+u^{\prime}$). The overbar indicates the time averaging of the
flow quantities. The symmetric Reynolds stress tensor can be expressed by
applying an eigenspace decomposition as

$\tau_{ij}=k\left(a_{ij}+\frac{2}{3}\delta_{ij}\right)=k\left(v_{in}\Lambda_{nl}v_{jl}+\frac{2}{3}\delta_{ij}\right)\textrm{
.}$ (1)

Equation 1 includes the split into the anisotropy tensor $a_{ij}$ and the
isotopic part of $\tau_{ij}$, while
$k=\frac{1}{2}\overline{u_{i}^{\prime}u_{i}^{\prime}}$ being the turbulent
kinetic energy. The eigenspace decomposition provides the eigenvector matrix
$v$ and the diagonal eigenvalue matrix $\Lambda$. Iaccarino and co-workers
proposed a strategy to perturb the eigenvalues and eigenvectors in Equation 1,
resulting in a perturbed state of the Reynolds stress tensor

$\tau_{ij}^{*}=k\left(v_{in}^{*}\Lambda_{nl}^{*}v_{jl}^{*}+\frac{2}{3}\delta_{ij}\right)\textrm{
.}$ (2)

Every physical, realizable state of the Reynolds stress tensor, meaning to be
positive semidefinite, can be mapped onto barycentric coordinates

$\mathbf{x}=\mathbf{x}_{1c}\frac{1}{2}\left(\lambda_{1}-\lambda_{2}\right)+\mathbf{x}_{2c}\left(\lambda_{2}-\lambda_{3}\right)+\mathbf{x}_{3c}\frac{1}{2}\left(3\lambda_{3}+1\right)\quad\mathrm{with}\quad\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}\textrm{
,}$ (3)

which is essentially a linear transform according to
$\mathbf{x}=\mathbf{B}\lambda$ within the triangle defined by the corners
$\mathbf{x}_{1c}$, $\mathbf{x}_{2c}$ and $\mathbf{x}_{3c}$ (see 1) [1] .
The eigenvalue perturbation makes use of the fact, that the Reynolds stress
tensor is limited between the one, two and three component corner of the
barycentric triangle, where the limiting states reflect the number of non-zero
eigenvalues of the tensor. As shown in Figure 1(a), the eigenvalue
perturbation is defined as a shift in barycentric coordinates towards each of
the limiting states to location $\mathbf{x}^{*}$, according to

$\mathbf{x}^{*}=\mathbf{x}+\Delta_{B}\left(\mathbf{x}_{(t)}-\mathbf{x}\right)\textrm{
.}$ (4)

The relative distance $\Delta_{B}$ controls the magnitude of eigenvalue
perturbation towards the corner state
$\mathbf{x}_{(t)}\in{\mathbf{x}_{\mathrm{1C}},\mathbf{x}_{\mathrm{2C}},\mathbf{x}_{\mathrm{3C}}}\\}$.
The perturbed eigenvalues $\lambda_{i}^{*}$ can be remapped afterwards.

(a) Relation between unperturbed and perturbed state

(b) Usage of data-driven perturbation strength

(c) Componentwise perturbation correction

Figure 1: Schematic representation of the eigenvalue perturbation approach

The creation of the perturbed eigenvector matrix $\Lambda^{*}$ is purely
motivated on the manipulation of the turbulent production term
$P_{k}=-\tau_{ij}\frac{\partial u_{i}}{\partial x_{j}}$ and is not discussed
further during this extended abstract. Taking only the eigenvalue perturbation
into account, a designer needs to run three distinct RANS simulations in order
to get the full estimate of the turbulence model uncertainty bounds, since the
methodology seeks to sample from the three extremal states in the barycentric
map.

###  2.2 Data-driven approach

The data-free approach is a purely physics-based methodology, aiming for
understandable uncertainty bounds for the turbulence modelling community. One
of the major drawbacks of the proposed method is the fact, that a user has to
choose the degree of uncertainty by selecting $\Delta_{B}$ before each
investigation. The perturbation amplitude, closely related to the degree of
uncertainty, has only to become significant in flow regime, which contravene
the assumptions during formulation of the turbulence model. Thus, enabling a
spatially varying perturbation of the Reynolds stress tensor seems to be
worthwhile.
The prediction of the correct anisotropy was investigated in an a priori study
using random forests by Wang et al. [11]. Heyse et al. propose a strategy to
combine a random forest model with the forward eigenvalue perturbation
approach [5]. Physical flow features are extracted to train a random forest in
order to predict the local perturbation strength
$p=|\mathbf{x}_{\mathrm{Data}}-\mathbf{x}_{\mathrm{RANS}}|=|\mathbf{x}^{*}_{\mathrm{RANS}}-\mathbf{x}_{\mathrm{RANS}}|$
(see Figure 1(b)). Forward propagating CFD simulations follow training the
model, where the predicted perturbation strength is used to modify the
Reynolds stress towards the same three limiting states as in the data-free
approach.
For the sake of interpretability and usability, an ensemble learning random
regression forest is chosen in the presented work as well [2]. The python
library scikit-learn is used to train the random forest and evaluate its
prediction. We will investigate the proposed combination of a machine learning
model and the perturbation approach within TRACE. The aim is to go one step
further and merge the work of Wang et al. [11] and Heyse et al. [5].
As a first step a machine learning model should predict the componentwise
perturbation strength $\mathbf{p}_{\mathrm{corr}}=(p_{x_{\mathrm{corr}}}$,
$p_{y_{\mathrm{corr}}})$ with respect to data, as illustrated in Figure 1(c).
The effect of applying corrected eigenvalues of the Reynolds stresses can be
observed by propagating them through the solver.
Furthermore, the orientation of the Reynolds stress tensor in terms of
mutually orthogonal eigenvectors, can be described as a rigid body rotation
with three degrees of freedom. RANS results show a discrepancy to the
orientation of the eigenvectors obtained from high-fidelity data, e.g. scale
resolving simulations. Consequently, intrinsic Tait-Bryan angles
$\alpha,\beta,\gamma$ in the $z-y^{\prime}-x^{\prime\prime}$ convention
according to Rothmayr and Hodges [10] will be used to describe the eigenvector
rotation from RANS results towards the data. The combination of corrected
eigenvalues and eigenvectors completes a predictive framework to model
anisotropy discrepancy and propagate the predicted Reynolds stresses to mean
flow quantities. In terms of selected flow features, we extended physical
quantities to the proposed exhaustive invariant feature list in the work of
Ling et al. [8] and Wang et al. [11].

##  3 Channel flow results

The fully developed turbulent channel flow offers a possibility to demonstrate
the coupling of the uncertainty quantification (UQ) perturbation approach with
a machine learning model. The discrepancies between RANS and direct numerical
simulations (DNS) of Lee and Moser [7] at different Reynolds numbers
$Re_{\tau}$ serve as training and testing data for the random forests.

Table 1: Selected hyperparameters for the random forest regressors | Target quantities
---|---
Hyperparameter | $p$ | $\mathbf{p}_{\mathrm{corr}}$ |  $\mathbf{p}_{\mathrm{corr}}$ \+ $\alpha,\beta,\gamma$
max tree depth | 6 | 9 | 9
min sample count | 6 | 4 | 4
max active features | 3 | 3 | 3
number of trees | 30 | 15 | 30

Before the actual training of the model was conducted, the impact of four
different hyperparameters on the accuracy of the random regression forest
model is evaluated: the maximum tree depth, the minimum sample count, the
maximum number of active features and the number of trees. The selected
parameters for the respective target quantities are based on hyperparameter
studies (see Table 1). The random forest regression models are trained on
$Re_{\tau}\in\\{180,550,2000,5200\\}$ using the mean squared error to
determine the quality of each split, whereas the models are evaluated in
forward CFD-simulations at $Re_{\tau}=1000$. The discrepancies with respect to
barycentric coordinates are moderate in the channel center and start to
increase close to the wall due to strong anisotropy of turbulence. The two-
equation, linear eddy viscosity Menter SST $k$-$\omega$ turbulence model is
considered as the baseline model in the current study [9].

(a) Streamwise velocity profile

(b) Barycentric coordinates

Figure 2: Performance of random forest models in combination with the
perturbation uncertainty quantification approach at channel flow simulations
of $Re_{\tau}=1000$

The advantage of the data-driven UQ approach, determining a relative
perturbation strength $p$ can be seen in Figure 4. The predicted uncertainty
bounds in Figure 2(a) are significantly less conservative compared to the
data-free method with a chosen $\Delta_{B}=1.0$, while the traces of the
barycentric coordinates look plausible for each of the perturbed solutions
(see Figure 2(b)). As expected, the baseline model aligns itself with the
plane strain line, whereas the perturbed solutions tend to the limiting states
of turbulence in the corners. Additionally, the random forest model,
predicting the perturbation direction vector $\mathbf{p}_{\mathrm{corr}}$, was
evaluated in forward CFD solver mode. Although it struggles to predict the
$x-$component of the perturbation, the resulting eigenvalues of the Reynolds
stress tensor in the simulation run are satisfying, as illustrated in Figure
2(b). Unfortunately, guaranteeing eigenvalues close to the truth does not seem
to be sufficient in order to obtain the correct quantity of interest, in this
case the velocity profile.

(a) rotation around $z$-axis ($\alpha$)

(b) rotation around $y^{\prime}$-axis ($\beta$)

(c) rotation around $x^{\prime\prime}$-axis ($\gamma$)

(d) model accuracy for $\alpha$

(e) model accuracy for $\beta$

(f) model accuracy for $\gamma$

Figure 3: Overview of the eigenvector rotation (RANS $\rightarrow$ DNS) for
channel flow

Moreover, we start to explore the discrepancy field regarding the eigenvector
orientation of the RANS and the DNS data in Figure 3(a), 3(b) and 3(c). Unlike
the differences in barycentric coordinates, the eigenvector misalignment seems
to be significant in the log-law region. A random forest regressor was trained
to predict the Tait-Bryan angles, as discussed above. Except for some
outliers, the model is able to predict the correct rotation angles within an
acceptable range (see Figure 3(d), 3(e) and 3(f)). A combination of
propagating true eigenvalues and eigenvectors
($\mathbf{p}_{\mathrm{corr}}+\alpha,\beta,\gamma$) in the forward CFD
application led to an improved prediction of the Reynolds stresses, as
illustrated in Figure 4(a), although actually only the anisotropy tensor was
corrected without modifying the turbulent kinetic energy explicitly (see
Equation 2). Due to the continual response of the equations to the changed
Reynolds stresses, the turbulence model produces more accurate turbulent
kinetic energy compared with the DNS data (see Figure 4(c)).

(a) Main components of the Reynolds stress tensor

(b) Shear stress component of the Reynolds stress tensor

(c) Turbulent kinetic energy

Figure 4: Evaluation of model predicted turbulent quantities at channel flow
of $Re_{\tau}=1000$

It might be possible, that the non-smoothness of the model prediction causes
an issue on the derivatives of the Reynolds stresses, entering the RANS
equations, leading to a descrepancy of the flow field [11]. In order to answer
the research question, whether a propagation of corrected Reynolds stresses is
able to improve the accuracy of quantities of interests, we propagated
interpolated Reynolds stresses of the DNS through the RANS equations and
observed a perfect match for the velocity profile at $Re_{\tau}=1000$.
Additionally, adding random noise to these interpolated Reynolds stresses and
not noticing significant change in the velocity profile (not shown here),
refute the non-smoothness theory. Therefore, the different behavior of the
shear stress component between $10\leq y^{+}\leq 40$ in Figure 4(b) is
accountable for the misalignment of the velocity profile for the simulations
containing a machine learned correction of the eigenvalues
($\mathbf{p}_{\mathrm{corr}}$) and of the entire anisotropy tensor
($\mathbf{p}_{\mathrm{corr}}+\alpha,\beta,\gamma$). Further work needs to deal
with an additional correction for the turbulent kinetic energy based on a
machine learning model, which potentially lead to more accurate quantities of
interest.

##  4 Conclusion

This work presents the possibility to estimate data-driven uncertainty bounds
for turbulence model with state-of-the-art methods in DLR’s CFD solver TRACE.
Moreover, recent approaches in literature acted as thought-provoking impulse
to train machine learning models on parameters of the initial perturbation
approach. The existing concerns were confirmed, that propagating the Reynolds
stress corrections is challenging and might not easily lead to correct mean
flow quantities of interest. Nevertheless, getting to know the turbulence
model uncertainties and especially localize flow regions, where those are
present, will help making design exploration and certification more reliable.

#### Impact statement

Complex engineering designs always feature turbulent flows. Although
computational power is increasing massively during the past years, scale
resolving simulations for design optimization studies seems to be still a long
way off. RANS simulation, whose accuracy is heavily dependent on turbulence
models, will remain state-of the art in the upcoming years. The present work
aims to consolidate the arisen methods in the field of turbulence model
uncertainty quantification in combination with machine learning methods.
Additionally, designing and implementing a framework to easily conduct
uncertainty estimation for turbulence models was a major goal of this work.
The authors do not see any ethical concern regarding the present work.

## References

  * Banerjee et al. [2007] S. Banerjee, R. Krahl, F. Durst, and C. Zenger.  Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches.  _Journal of Turbulence_ , 8:N32, 2007.  doi: 10.1080/14685240701506896.  URL https://doi.org/10.1080/14685240701506896.
  * Breiman [2001] L. Breiman.  Random forests.  _Machine Learning_ , 45:5–32, 10 2001.  doi: 10.1023/A:1010950718922.
  * Duraisamy et al. [2019] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence modeling in the age of data.  _Annual Review of Fluid Mechanics_ , 51(1):357–377, Jan 2019.  ISSN 1545-4479.  doi: 10.1146/annurev-fluid-010518-040547.  URL http://dx.doi.org/10.1146/annurev-fluid-010518-040547.
  * Emory et al. [2013] M. Emory, J. Larsson, and G. Iaccarino.  Modeling of structural uncertainties in reynolds-averaged navier-stokes closures.  _Physics of Fluids_ , 25(11):110822, 2013.  doi: 10.1063/1.4824659.  URL https://doi.org/10.1063/1.4824659.
  * Heyse et al. [2021] J. Heyse, A. Mishra, and G. Iaccarino.  Estimating rans model uncertainty using machine learning.  _Journal of the Global Power and Propulsion Society_ , pages 1–14, 05 2021.  doi: 10.33737/jgpps/134643.
  * Iaccarino et al. [2017] G. Iaccarino, A. Mishra, and S. Ghili.  Eigenspace perturbations for uncertainty estimation of single-point turbulence closures.  _Physical Review Fluids_ , 2, 02 2017.  doi: 10.1103/PhysRevFluids.2.024605.
  * Lee and Moser [2015] M. Lee and R. D. Moser.  Direct numerical simulation of turbulent channel flow up to $\mathit{Re}_{{\it\tau}}\approx 5200$.  _Journal of Fluid Mechanics_ , 774:395–415, 2015.  doi: 10.1017/jfm.2015.268.
  * Ling et al. [2016] J. Ling, R. E. Jones, and J. Templeton.  Machine learning strategies for systems with invariance properties.  _J. Comput. Phys._ , 318:22–35, 2016.
  * Menter et al. [2003] F. Menter, M. Kuntz, and R. Langtry.  Ten years of industrial experience with the sst turbulence model.  _Heat and Mass Transfer_ , 4, 01 2003.
  * Roithmayr and Hodges [2016] C. Roithmayr and D. Hodges.  Dynamics: Theory and application of kane’s method.  2016\.
  * Wang et al. [2018] J.-X. Wang, J.-L. Wu, H. Xiao, and E. Paterson.  Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework.  _Physical Review Fluids_ , 3(7), Jul 2018.  ISSN 2469-990X.  doi: 10.1103/physrevfluids.3.074602.  URL http://dx.doi.org/10.1103/PhysRevFluids.3.074602.