# Solving Two-Player General-Sum Game Between Swarms

Mukesh Ghimire1, Lei Zhang1, Wenlong Zhang2, Yi Ren1 and Zhe Xu1 *This work
was not supported by any organization1M. Ghimire, L. Zhang, Y. Ren, and Z. Xu
are with the School of Engineering, Matter, Transport, and Energy, Arizona
State University, Tempe, AZ 85287, USA. Email: {mghimire, lzhan300, yiren,
<EMAIL_ADDRESS>Zhang is with the School of Manufacturing Systems and
Networks, Arizona State University, Mesa, AZ 85212, USA. Email:
<EMAIL_ADDRESS>
###### Abstract

Hamilton-Jacobi-Isaacs (HJI) PDEs are the governing equations for the two-
player general-sum games. Unlike Reinforcement Learning (RL) methods, which
are data-intensive methods for learning value function, learning HJ PDEs
provide a guaranteed convergence to the Nash Equilibrium value of the game
when it exists. However, a caveat is that solving HJ PDEs becomes intractable
when the state dimension increases. To circumvent the curse of dimensionality
(CoD), physics-informed machine learning methods with supervision can be used
and have been shown to be effective in generating equilibrial policies in two-
player general-sum games. In this work, we extend the existing work on agent-
level two-player games to a two-player swarm-level game, where two sub-swarms
play a general-sum game. We consider the Kolmogorov forward equation as the
dynamic model for the evolution of the densities of the swarms. Results show
that policies generated from the physics-informed neural network (PINN) result
in a higher payoff than a Nash Double Deep Q-Network (Nash DDQN) agent and
have comparable performance with numerical solvers.

##  I INTRODUCTION

Swarms, or groups of robots have the ability to carry out complex tasks that
are difficult for a single agent. Swarms have been deployed to perform diverse
tasks, for example, surveillance and reconnaissance in a military capacity,
search and rescue, and even entertainment in the form of light shows [1].
Other applications include task allocations [2], complex formation [3], smart
farming [4], and decision making [5]. As new applications emerge every day, it
is imminent that there will be multiple groups of swarms that may be operating
in the same region with their own objective. One of the most trivial examples
would be a zero-sum game with two swarm groups with competing objectives. Most
of the work in swarms focuses on studying swarm behaviors in a more single-
agent setting, where interactions are limited within the swarm group. While
some work have proposed game theoretic frameworks to solve problems in
communication networks [6, 7], they do so in a distributed fashion. We are
interested in finding a high level control strategy for swarms with
arbitrarily large population.

Modeling of swarms can be broadly divided into macroscopic and microscopic
models. Macroscopic models are invariant to the number of agents whereas
individual-level microscopic models change with the number of agents [8], and
become intractable as the number of agents gets large [9]. As such macroscopic
models can be thought of as robust models that can be used to study large-
scale swarms. A building block for studying such swarms is the Kolmogorov
forward equation which describes the evolution of the density of a stochastic
process. We define density as the ratio of the population of swarms in a
region to its total population. Under certain conditions, the Kolmogorov
forward equation can be applied to encode the macroscopic behavior of a swarm.

Multi-agent interactions can often be represented as zero-sum or general-sum
differential games depending on the objective [10]. The Nash Equilibrial
values of these games are viscosity solutions to the Hamilton-Jacobi-Isaacs
(HJI) equations [11]. HJI PDEs are conventionally solved on a mesh and often
suffer from the curse of dimensionality as the state dimension increases [12].
We apply the differential game theory to the density control of multiple swarm
groups and study the interaction that emerges.

In this work, we formulate a general-sum game between swarm groups which are
modeled as a density that evolves in continuous time. We use the power of
neural network being the universal function approximator [13] to learn the
nash equilibrial values governed by the HJI equation. Our contribution to the
existing research on learning and control of the swarm systems is the
introduction of a game-theoretic approach to model interactions among multiple
swarm groups. We assess the effectiveness of our approach by comparing its
performance against both a widely employed reinforcement learning method for
games and a numerical solver.

##  II RELATED WORK

Swarm Control. Control of swarms is of great interest to the research
community. An existing challenge in this line of research is developing models
and control mechanisms for large-scale swarms [8]. [9] devised an optimal
control strategy for controlling a large-scale swarm to a target distribution.
[14] used a leader-follower framework for herding a robotic swarm to a desired
distribution. Another work that is in the spirit of our work is [15], which
applied Pontryagin’s Maximum Principle for control of a large-scale robotic
population in an optimal control setting. In addition, [6] and [7] proposed a
game-theoretic approach on a graph to solve the Coverage game. Another common
line of work includes using a Markov chain that models the evolution of the
density distribution. [16] used the Markov chain and provided a probabilistic
control algorithm for swarms of agents subject to some temporal logic
specifications. In this paper, we study the interaction between two large-
scale swarm groups in an environment where the swarm groups have their own
objective. We formulate a general-sum game between these swarm groups and
determine their optimal policy.

Solving HJI PDEs using Deep Learning. Recent works have considered using
autoregressive methods such as physics-informed machine learning to learn
values of zero-sum [17], and general-sum [10] differential games. [10]
extended [17] from zero-sum game with continuous value function to a general-
sum game with discontinuous values with respect to states and time. Nash
Equilibrium values of general-sum differential games satisfy HJI PDEs, which
makes the residual of HJI PDEs a good candidate for a loss function used in
physics-informed neural networks. By minimizing the PDE residuals alongside
the boundary conditions, a neural network, often a deep one, is trained to
predict the Nash equilibrial values associated with the game. Since Deep
Learning methods are well known for their scalability to high dimensional
problems (see [18], [19] for details), we leverage this method to solve the
problem of learning the value of general-sum games.

Multi-Agent Reinforcement Learning (MARL).  MARL, unlike single-agent
reinforcement learning, addresses the problem of decision-making involving
multiple agents that operate in a common environment. Standard Rollout
algorithm is used where the problem is reformulated as a single-agent by using
a joint action space [20]. Standard Rollout algorithm has been extended to
Multi-agent Rollout algorithm in order to reduce the complexity arising from
the joint action space [21]. While most algorithms are in the spirit of
optimal control, some works exist that have formulated multi-agent decision-
making problems as a game [22], [23]. [22] proposed a modified Q-learning
method where Q-values are defined on joint action space and are updated
following Nash Equilibrial strategies. [23] extended this idea from a tabular
method to a scalable one using the power of deep neural networks. In contrast
to the case studies that are discussed in these two related works, our case
studies have significantly larger action sets.

This paper is organized as follows. In Section III, we present the notations
that appear in the paper frequently. In Section IV, we discuss the fundamental
concepts, assumptions, and challenges that are crucial to the development of
the paper. We also motivate the problem formulation in this section. In
Section V, we discuss the algorithms that are central to the contributions
highlighted in the paper. In Section VI, we test our algorithm on different
case studies. Finally in Section VII, we conclude with some limitations of the
current work and the possible future directions.

##  III NOTATIONS

Borrowing notations from graph theory, we denote a directed graph by the tuple
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ containing a set of $M$ vertices,
$\mathcal{V}=\\{1,\dots,M\\}$, and a set of $N_{\mathcal{E}}$ edges,
$\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$, where
$e=(i,j)\in\mathcal{E}$ if there is an edge from vertex $i\in\mathcal{V}$ to
vertex $j\in\mathcal{V}$. A source map is defined as
$S:\mathcal{E}\rightarrow\mathcal{V}$ and a target map as
$T:\mathcal{E}\rightarrow\mathcal{V}$ for which $S(e)=i$ and $T(e)=j$,
whenever $e=(i,j)\in\mathcal{E}$. The graph $\mathcal{G}$ is said to be
bidirected if $(v,w)\in\mathcal{E}$ implies that $(w,v)\in\mathcal{E}$ for all
$v,w\in\mathcal{V}$.

We follow [10]’s notation system and the implementation therein for HJI
equations. Let $\mathcal{X}_{i}$ and $\mathcal{U}_{i}$ denote the state and
action space respectively for Player $i$. The time-invariant state dynamics of
Player $i$ is denoted by
$\dot{\mathbf{x}}_{i}(t)=f(\mathbf{x}_{i}(t),\mathbf{u}_{i}(t))$ for
$\mathbf{x}_{i}(t)\in\mathcal{X}_{i}$ and $\mathbf{u}_{i}\in\mathcal{U}_{i}$.
Given a time horizon $T$, the instantaneous and terminal losses of Player $i$
are denoted by $l$, and $c(\mathbf{x}_{i},\mathbf{x}_{-i})$ respectively. Note
that $l$ is a constant in this formulation, however, in general it is a
function of state $(l(\mathbf{x}_{i},\mathbf{x}_{-i}))$. For a complete-
information general-sum differential game between two players, the value
function for each player is
$\nu_{i}(\cdot,\cdot,\cdot):\mathcal{X}_{i}\times\mathcal{X}_{-i}\times[0,T]\rightarrow\mathbb{R}$.
We will adopt shorthands $f_{i}$, $l_{i}$, $c_{i}$, and $\nu_{i}$ respectively
to denote player-wise dynamics, losses, and the value. We use
$\mathbf{a}_{i}=(a_{i},a_{-i})$ to concatenate ego and the other player’s
variables. $\nabla_{x}$ denotes partial derivative with respect to $x$.

##  IV PRELIMINARIES

Swarm Modeling. We consider two homogeneous sub-swarms, each containing $N$
agents that occupy the regions or vertices $\mathcal{V}$. The swarms evolve in
continuous time over this region, which can be denoted by vertices on the
graph. We denote the graph by $\mathcal{G}=(\mathcal{V},\mathcal{E})$, with
vertices $\mathcal{V}$ denoting the regions where swarms reside, and
$\mathcal{E}$ representing the edges along which the swarms can transition.

Let, $X_{v}^{i}(t)$ be the number of agents of sub-swarm $i\in\\{1,2\\}$ in
the region $v\in\mathcal{V}$. The fraction (or empirical distribution) of the
sub-swarms $i$ at location $v\in\mathcal{V}$ at time $t$ is calculated as
$\frac{1}{N}X_{v}^{i}(t)$. Let
$\Omega=\\{\mathbf{y}\in\mathbb{R}^{n}:\sum_{i=1}^{n}y_{i}=1$}. As
$N\rightarrow\infty$, the empirical distribution converges to a deterministic
quantity $\mathbf{x}(t)\in\Omega$, which can be used as a state of the sub-
swarm $i$. We denote the state of the sub-swarm $i$ as a vector
$\boldsymbol{x}_{i}(t)$, whose each entry $x_{i,v}(t)$ denotes the fraction of
agents in the region $v\in\mathcal{V}$ at time $t$. We use the Kolmogorov
forward equation, also known as the mean-field model as in [14] to evolve the
state of the sub-swarm:

$\dot{\boldsymbol{x}}(t)=\sum_{e\in\mathcal{E}}u_{e}(t)\mathbf{B}_{e}\boldsymbol{x}(t),\quad\boldsymbol{x}(0)=\boldsymbol{x}^{0}\in\Omega$
(1)

where $\mathbf{B}_{e}$ are the control matrices with the following entries

$\displaystyle B_{e}^{ij}=\begin{cases}-1&\text{ if }i=j=S(e)\\\ 1&\text{ if
}i=T(e),j=S(e)\\\ 0&\text{ otherwise }\end{cases}$

Hamilton-Jacobi-Isaacs Equations. The Nash equilibrial values of a two-player
general-sum differential game, when exist, solve the following HJI equations
$(H)$ and satisfy the following boundary condition $(B)$ [24]:

$\displaystyle
H(\nu_{i},\nabla_{\mathbf{x}_{i}}\nu_{i},\mathbf{x}_{i},t):=\nabla_{t}\nu_{i}+\nabla_{\mathbf{x}_{i}}^{\top}\textbf{f}_{i}-l_{i}=0$
(2) $\displaystyle
B(\nu_{i},\mathbf{x}_{i}):=\nu_{i}(\mathbf{x}_{i},T)-c_{i}=0,\quad\text{ for
}i=1,2.$

The player’s policies for ego agent (and other) are derived by maximizing the
equilibrial Hamiltonian
$h_{i}(\mathbf{x}_{i},\nabla_{\mathbf{x}_{i}}\nu_{i},t)=\nabla_{\mathbf{x}_{i}}\nu_{i}^{\top}f_{i}-l_{i}$;
$u_{i}=\arg\max_{u\in\mathcal{U}_{i}}\\{h_{i}\\}$, and
$u_{-i}=\arg\max_{u\in\mathcal{U}_{-i}}\\{h_{-i}\\}$.

Pontryagin’s Maximum Principle: We can use the PMP equation to generate open-
loop policies, which is often more tractable than solving HJI equations. These
open-loop policies can be used as a basis for evaluating the learned closed-
loop policies. For a given initial state
$(\boldsymbol{x}_{1},\boldsymbol{x}_{2})\in\mathcal{X}_{1}\times\mathcal{X}_{2}$,
we obtain the open-loop policies by solving the following boundary value
problem (BVP) according to PMP:

$\displaystyle\dot{\mathbf{x}_{i}}$
$\displaystyle=\mathbf{f}_{i},\quad\mathbf{x}_{i}[0]=\bar{\mathbf{x}}_{i},$
(3) $\displaystyle\dot{\boldsymbol{\lambda}}_{i}$
$\displaystyle=-\nabla_{\mathbf{x}_{i}}\mathbf{h}_{i},\quad\boldsymbol{\lambda}_{i}[T]=-\nabla_{\mathbf{x}_{i}}c_{i},$
$\displaystyle u_{i}$
$\displaystyle=\arg\max_{u\in\mathcal{U}_{i}}\\{h_{i}\\},\quad i=1,2$

$\boldsymbol{\lambda}_{i}$ are the time-dependent co-states for both players
concatenated into one. The costates connect PMP and HJI via
$\boldsymbol{\lambda}_{i}=\nabla_{\mathbf{x}_{i}}\nu_{i}$. Note that the
solutions to Eq. (3) are unique to the initial states.

Learning Values of General-sum Differential Games. Directly learning values of
a zero-sum differential game in a self-supervised fashion was first explored
in [17]. In this approach, we directly fit a neural network to satisfy the
governing HJI PDE. Let
$\hat{\nu}_{i}(\cdot,\cdot,\cdot):\mathcal{X}_{i}\times\mathcal{X}_{-i}\times[0,T]\rightarrow\mathbb{R}$
be a neural network parameterized by $\theta$ that approximates $\nu_{i}$.
With
$\mathcal{D}=\left\\{\left(\boldsymbol{x}_{1}^{(k)},\boldsymbol{x}_{2}^{(k)},t^{(k)}\right)\right\\}_{k=1}^{K}$
representing the uniform samples in
$\mathcal{X}_{i}\times\mathcal{X}_{-i}\times[0,T]$, the loss function that
guides the learning of the general-sum value is:

$\displaystyle\min_{\theta}L_{1}(\hat{\nu_{i}};\theta):=\sum_{k=1}^{K}\left|H\left(\hat{\nu}_{i}^{(k)},\nabla_{\mathbf{x}_{i}}\hat{\nu}_{i}^{(k)},\mathbf{x}_{i}^{(k)},t^{(k)}\right)\right|$
(4)
$\displaystyle\quad+\;C_{1}\left|B\left(\hat{\nu}_{i}^{(k)},\mathbf{x}_{i}^{(k)}\right)\right|$

where, $\hat{\nu}_{i}^{(k)}$ is the output of the neural network,
$\hat{\nu}_{i}^{(k)}=\hat{\nu}_{i}\left(x_{i}^{(k)},x_{-i}^{(k)},t^{(k)}\right)$.
$C_{1}$ balances the HJI PDE loss and the boundary loss.

Note that at each training iteration, we compute the control policies for each
player by maximizing their equilibrial Hamiltonian. We refer the readers to
[10] for more details on challenges and different methods of learning values
using this approach. For the purpose of this paper, we only consider the self-
supervised learning method to learn the value function.

##  V METHODOLOGY

In this section, we present the algorithms for generating open-loop
trajectories via solving BVP and training the value network using self-
supervised learning method.

###  V-A BVP Solver

In this subsection, we discuss about generating open-loop trajectories. We use
scipy’s solve_BVP function to solve Eq. (3). To successfully solve the BVP, we
first compute the analytical expressions for the quantities in Eq. (3). These
are (1) the augmented dynamics which contains the state dynamics
$(\dot{\mathbf{x}})$ and the co-state dynamics $(\dot{\boldsymbol{\lambda}})$,
and (2) the boundary conditions. Note we also need a sub-routine to compute
the control policies for each of the agents. Furthermore, convergence requires
good guesses of the solution along the trajectory. To do so, given an initial
state, we solve the state and co-state trajectories by solving their
respective ODEs obtained from PMP. An alternative to this, albeit slower, is
implementing time-marching [25].

###  V-B Multi-Agent RL using Nash Q-Learning

We also formulate the 2-regions case as a Reinforcement Learning problem and
train the sub-swarms using Nash Q-Learning. Instead of using a tabular method,
we use a neural network, especially a Double Deep Q-Network that approximates
the Q values. We implement the vanilla Nash Q-Learning following [22]. The
Nash Q-function for a learning agent $i$ is defined on the state of the
system, and its action along with the actions of the other players, whereas in
the Q-learning algorithm, the Q-function is only a function of the agent.
Furthermore, in Q-learning, agents update their Q-values by accounting for the
future values as a result of maximizing their Q-values, whereas in Nash
Q-learning, agents update their Q-values by following a Nash equilibrium
strategy. More formally, the update is as follows:

$\displaystyle Q^{i}_{t+1}(s,a^{1},\dots,a^{n})$
$\displaystyle=(1-\eta_{t})Q^{i}_{t}(s,a^{1},\dots,a^{n})+$ (5)
$\displaystyle\eta[r^{i}_{t}+\beta NashQ^{i}_{t}(s^{\prime})]$

where, $\eta$ is the learning rate, and $\beta$ is the discount rate. $Nash$
is the operator of choice that computes the Nash equilibrium of the stage game
at state $s^{\prime}$. In this work, we use support enumeration as a nash
operator.

Unlike the vanilla Nash Q-learning, we use two neural networks that
approximate the Q-function of the sub-swarms in the 2-regions case. Each sub-
swarm has an action space of 4 – (0, 0), (0, 1), (1, 0), (1, 1). Note that
each agent must also know the Q-values of the other agent to compute the Nash
Equilibrium. However, this information is not available, so each agent makes a
conjecture about the other agent’s Q-values. We achieve this by defining the
Q-function on the joint action space of both sub-swarms. The result is a
$4\times 4$ matrix of Q-values for each sub-swarm. At each learning stage,
each sub-swarm solves a bi-matrix game as follows:

TABLE I: Bi-matrix stage game for Nash Q-Learning. |  |  | Agent $-i$ |
---|---|---|---|---
|  | (0, 0) | (0, 1) | (1, 0) | (1, 1)
Agent $i$ | (0, 0) | $a_{1},b_{1}$ | $a_{2},b_{2}$ | $a_{3},b_{3}$ | $a_{4},b_{4}$
(0, 1) | $a_{5},b_{5}$ | $a_{6},b_{6}$ | $a_{7},b_{7}$ | $a_{8},b_{8}$
| (1, 0) | $a_{9},b_{9}$ | $a_{10},b_{10}$ | $a_{11},b_{11}$ | $a_{12},b_{12}$
| (1, 1) | $a_{13},b_{13}$ | $a_{14},b_{14}$ | $a_{15},b_{15}$ | $a_{16},b_{16}$

where, $a$ is the q-value for agent $i$, and $b$ is the q-value for agent $-i$
as conjectured by agent $i$. Methods such as Lemke-Howson [26] or support
enumeration can be used to compute the Nash Equilibria for the bi-matrix game.
Note that we use NashPy’s implementation of support enumeration [27]. Note
that table I gets significantly larger when we consider higher dimensional
case studies. Solving the bi-matrix game becomes a challenge for these cases.
As a result we only compare the performance of Nash DDQN with PINN for the
case study with 2 regions. We describe the implementation details for Nash
DDQN briefly in Algorithm 1.

input : # Episodes $B>0$, Minibatch size $m>0$, Episode length $N$,
Exploration rate $\varepsilon$, Discount rate $\beta$, soft-update rate $\tau$

initialize : Replay Buffer $\mathcal{D}$, Q-networks
$(Q_{\theta_{1}},Q_{\theta_{2}})$, and respective target networks
$(Q_{\theta_{1^{\prime}}},Q_{\theta_{2^{\prime}}})$

1 for  _Episode b $\leftarrow$ 1 to B_ do

2 Reset Simulation, sample initial state $s$ randomly.

3 for  _Episode steps t $\leftarrow$ 1 to N_ do

4 Select actions $u_{1},u_{2}$ via $NashQ_{\theta_{i}}(s)$ with $\varepsilon$
random exploration;

5 Store $\mathcal{D}\leftarrow(s,u_{1},u_{2},r_{1},r_{2},s^{\prime})$;

6 Sample transitions $\mathcal{Y}=\\{y\\}_{i=1}^{m}$ from $\mathcal{D}$;

7 Set $Q_{i}\leftarrow Q_{\theta_{i}}(s,u_{1},u_{2}),i=\\{1,2\\}$;

8 Perform gradient descent steps on $(Q_{i}(s,u_{1},u_{2})-[r_{i}+\beta
NashQ_{\theta_{i}}(s^{\prime})])^{2}$, $i=\\{1,2\\}$;

9 Update target network parameters:
$\theta_{i^{\prime}}\leftarrow\tau\times\theta_{i}+(1-\tau)\times\theta_{i^{\prime}}$

10 end for

11

12 end for

Algorithm 1 Nash DDQN

###  V-C Self-supervised learning

In this subsection, we discuss the method of self-supervised learning
algorithm using PINN to learn the value function. We implement curriculum
learning by first learning the value at the final time and gradually
increasing the time horizon starting from the end time. We provide a
simplified algorithm for training the value network. In our experiment, we use
ADAM to optimize the neural network parameter $(\theta)$ with a decaying step-
size (learning rate) $\gamma\in[2e-5,1e-6]$. We use a neural network with 3
hidden layers containing 64 neurons each with tanh activation function, with
the final layer being linear. For the higher dimensional case study we
increase the hidden layer to 5 and the neurons to 128 at each layer. We use
NVIDIA A100 with 2 GPUs for all the training in this paper. The dimension of
input $\mathcal{X}$ to the neural network depends on the number of regions in
the graph. We further reduced the dimension of the system using the fact that
the densities of each sub-swarm sum to 1. Hence, for the case with 2 regions,
the value function is a 3-dimensional function instead of 5, including time.

input : $T$, $num\\_epoch$, $pretrain\\_iters$, $k$, $n$

output : $V_{\theta}$

1 initialize neural network $V_{\theta}$;

2 sample $k$ samples of $\mathcal{X}\in\Omega$;

3 pre-train $V_{\theta}$ at the boundary $(t=0)$ for $pretrain\\_iters$
iterations;

4 set $itr=0$;

5 while  _$itr\leq num\\_epoch$_ do

6 sample $k$ samples of $\mathcal{X}\in\Omega$;

7 sample $k$ samples of $t\sim\mathcal{U}(0,T(itr)/{num\\_epoch})$;

8 append $n$ samples at the boundary $(t=0)$ to $\mathcal{X}$ and to $t$;

9 compute total loss ($L_{1}$) using PDE loss $(H)$, and boundary loss $(B)$;

10 update $\theta\leftarrow\theta+\gamma\nabla_{\theta}L_{1}$

11 end while

Algorithm 2 Self-supervised Learning

##  VI CASE STUDIES

Objective. Now, let us define the general-sum game that we study in this paper
along with the case studies. Formally, the payoff and the value of the game
for a player $i$ is defined as:

$\displaystyle
J_{i}(\boldsymbol{x}_{i}(0),\boldsymbol{x}_{-i}(0),u_{i},u_{-i})$
$\displaystyle=g(\boldsymbol{x}_{i}(T),\boldsymbol{x}_{-i}(T))$
$\displaystyle\dot{\boldsymbol{x}}_{i}$
$\displaystyle=f(\boldsymbol{x}_{i},u_{i})$ (6) $\displaystyle V_{i}^{*}$
$\displaystyle=\max_{u_{i}}J_{i}(\boldsymbol{x}_{i}(0),\boldsymbol{x}_{-i}(0),u_{i},u_{-i})$

Note that we do not have an instantaneous loss in the definition of the value,
$(l_{i}=0,\forall i)$. The dynamics equation, $f$ is different for the case
studies, which depends on the underlying graph.

We test our method in three different graphs representing both low and high
dimensional systems. We consider two sub-swarms with the same control
capabilities interacting in an environment containing $M$ regions
$(|\mathcal{V}|=M)$. We first consider $M=2$ (as a toy case) and $M=4$ regions
with the former represented by a bidirected graph and the latter represented
by a directed graph as shown in Fig. 1. We then consider a high dimensional
case (21 D), with $M=10$ regions, to show that the proposed method is easily
scalable.

Figure 1: (a) Bidirected graph with 2 regions or vertices. (b) Directed graph
with 4 regions. (c) Directed graph with 10 regions. The maximum transition
rates (control bounds) for all the edges is a constant 1, and the minimum rate
is 0.

The state of each sub-swarm is denoted by $\mathbf{x}_{i}\in\mathbb{R}^{M}$,
with $M$ being the number of regions. The goal of each sub-swarm is to achieve
higher proportion in any one of the regions (vertices in the graph). The
terminal payoff to sub-swarm $i$ is $g(\mathbf{x}_{i},\mathbf{x}_{-i})$ is a
Boltzmann operator $(S_{\alpha})$ defined on the element-wise difference
between the states of the two sub-swarms. We use the Boltzmann operator
instead of the $\max$ operator to make the function differentiable and
continuous. The following represents the terminal payoff to the sub-swarm
group $i=1$. To compute the terminal payoff for sub-swarm group $i=2$, simply
reverse the order of difference.

$\displaystyle
g(\mathbf{x}_{1},\mathbf{x}_{2})=S_{\alpha}(\mathbf{x}_{1}-\mathbf{x}_{2})=\frac{\sum_{j=1}^{M}\left(\mathbf{x}_{1}-\mathbf{x}_{2}\right)_{j}e^{\alpha(\mathbf{x}_{1}-\mathbf{x}_{2})_{j}}}{\sum_{j=1}^{M}e^{\alpha(\mathbf{x}_{1}-\mathbf{x}_{2})_{j}}}$
(7)

where $\alpha$ is the temperature parameter and as $\alpha\rightarrow\infty$,
it behaves like the $\max$ operator. We set $\alpha=1$ for the experiments in
this paper (unless mentioned explicitly). The gradient of $S_{\alpha}$
resembles that of a softmax:

$\displaystyle\nabla_{{x}_{i}}S_{\alpha}(x_{1},\dots,x_{M})$
$\displaystyle=\frac{e^{\alpha}x_{i}}{\sum_{j=1}^{M}e^{\alpha
x_{j}}}\Big{[}1+$ (8)
$\displaystyle\quad\quad\alpha(x_{i}-S_{\alpha}(x_{1},\dots,x_{M}))\Big{]}$

An important thing to note is that the states in our formulation are not
independent as they must sum to 1. We highlight how this brings about some
crucial changes to the formulation using the 2-regions case study. After
reducing the state dimension from two-dimensional to one, denote the state of
sub-swarm 1 with $x$ and that of sub-swarm 2 with $y$, which represents their
density in region 1. Then, the final payoff becomes:

$\displaystyle c_{i=1}=S(x-y)$
$\displaystyle=\frac{(x-y)e^{(x-y)}+(y-x)e^{(y-x)}}{e^{(x-y)}+e^{(y-x)}}$
$\displaystyle=(x-y)\frac{e^{(x-y)}-e^{-(x-y)}}{e^{(x-y)}+e^{-(x-y)}}$
$\displaystyle=(x-y)\tanh{(x-y)}$

And, its gradient with respect to the state variables are:

$\displaystyle\frac{\partial c_{i=1}}{\partial x}$
$\displaystyle=\tanh{(x-y)}+(x-y)\operatorname{sech}^{2}(x-y)$
$\displaystyle\frac{\partial c_{i=1}}{\partial y}$
$\displaystyle=-\frac{\partial c_{i=1}}{\partial x}$

###  VI-A Sub-swarms in 2 regions

We first solve a simple case with only 2 regions and identical homogeneous
sub-swarms. We compare the open-loop policy obtained by solving the BVP using
Eq. (3) with the policy obtained from the neural network. First, we show that
the closed-loop policy from the neural network results in approximately the
same final payoff for the swarm groups. We take 1000 initial states and find
the absolute error between the values at the final time both using the BVP
solver and the neural network. Fig. 2 shows that the resulting payoffs are
close to each other.

Training. For 2 dimensional state space, we sample 65k data and pretrain at
the boundary for 10k iterations. We then implement curriculum learning by
slowly increasing the time horizon and train for another 110k iterations. We
also ensure that there are at least 10k data points at the boundary after the
pretrain stage.

Figure 2: Absolute Error of Values at final time starting from same initial
states between BVP solution and neural network. $(V_{1}^{BVP}$, $V_{2}^{BVP})$
represent values of sub-swarms 1 and 2 obtained from BVP solver, and
$(V_{1}^{NN},V_{2}^{NN})$ are obtained from value network.

Results. Using the learned value network, we generate the closed-loop
trajectories for a given initial condition representing the density of swarm
groups in region 1. The results show that each swarm group will occupy one
separate region at the end of the game. Fig. 4 shows the trajectories of the
swarm groups (top corresponds to Swarm Group 1), over the time horizon of 2.5
seconds starting from 4 different initial conditions. This solution maximizes
payoffs of both swarm groups and therefore is one of the Nash Equilibrium
solutions.

We also observe that for some initial states value network and the BVP solver
result in different equilibria. Since the problem is symmetric, it is possible
that the swarms can occupy any one of two regions to guarantee a maximum
payoff as long as they pick separate regions. We show one example of such case
in Fig. 3. To further test that these solutions are indeed a Nash Equilibrium,
we generate 100 different trajectories using the learned value network as a
closed-loop controller. These trajectories are then fed into the BVP solver as
an initial guess. The resulting solution from the BVP solver converged to the
solution obtained from the neural network for 97 out of the 100 total
trajectories.

Figure 3: A case demonstrating possible equilibria for the 2 regions case. (a)
is the initial condition of the sub-swarms. (b) and (c) are the final
distributions of the sub-swarms in the environment obtained from the value
network and the BVP solver respectively. Both distributions are equilibria for
the sub-swarms as they result in approximately equal payoffs. Multiple
equilibria exist due to the symmetric nature of the game. Figure 4:
Interaction in two regions. Density evolution of swarm groups starting from
four different initial conditions, $\textbf{x}_{0}^{(1)}$,
$\textbf{x}_{0}^{(2)}$, $\textbf{x}_{0}^{(3)}$, and $\textbf{x}_{0}^{(4)}$,
which represent the densities of the swarm groups in region 1. The top figure
corresponds to sub-swarm 1 and the bottom corresponds to sub-swarm 2. Each
sub-swarm aggregates in one region each, which maximizes their payoffs.

###  VI-B Sub-swarms in 4 regions

Next, we consider an environment with 4 regions as shown in Fig. 1(c). Number
of swarm groups is 2 with same transition rates across all edges.

Training. For 4 regions case, we increase the training iterations to 200k with
20k pretrain steps. In addition to curriculum learning, we also increase the
data points gradually starting from 30k in the beginning and adding 10k data
points after every 10k iterations after the pretrain stage.

Results. As in the previous case study, the sub-swarms eventually find one
region each where they dominate. Fig. 5 shows the evolution of densities at
the four regions for four different initial conditions. The densities at the
final time for these four initial conditions are shown in Fig. 6. All cases
lead to the scenario where sub-swarms maintain dominance in one region.

We use the initial states of the 150 closed-loop trajectories generated using
the value network as a controller and solve the game using the BVP solver. We
then compare the open-loop trajectories of the converged cases with the
initial states from the closed-loop trajectories and report the final payoff
differences. Note that some initial states when solved using the BVP solver
may not converge to a solution. While, for the 2 regions case, all 150 initial
states achieved convergence, for the 4 regions case, 139 out of the total 150
initial states reached convergence.

We observe that for some initial states, the final distributions are not
optimal (see Fig. 6(a) and (c)). Sub-swarm-2 in Fig. 6(a) could get a higher
payoff by aggregating in region 3 instead of two regions (3 and 4). This
behavior can be attributed to the temperature parameter in the payoff
function. Since hardening the temperature parameter led to difficulty in the
convergence of BVP, we tested the effect of the temperature parameter using
the value network. We trained a value network with $\alpha=20$ in the final
payoff function. The results show that the sub-swarms proceed to aggregate in
one region resulting in the desired behavior. Fig 7 shows the final
distributions of the sub-swarms for the same four initial conditions mentioned
before.

We compare the final payoffs across the learning algorithms for both 2 and 4
regions case studies. The results are reported in table II. A negative number
means that the BVP solution resulted in a higher final payoff than the
respective learning method. Final payoffs using PINN were close to the that
obtained from the BVP solution whereas Nash DDQN resulted in a significantly
lower payoff on average.

Figure 5: Interaction in four regions. Density evolution of swarm groups from four different initial conditions (a - d). Red plot corresponds to the density of sub-swarm 1 at that particular region, and the blue corresponds to that of sub-swarm 2. Sub-Swarms tend to find at least one region where they maintain dominance. Figure 6: (a-d) Represents the density distribution for two sub-swarms at the final time for four different initial conditions. The initial conditions are the same as in Fig. 5. Figure 7: Final distributions of the sub-swarms when the temperature parameter $\alpha=20$. For the same initial conditions as in the previous cases, the sub-swarms attempt to aggregate in one region, which is desirable. TABLE II: Comparison of final payoffs from solving the interaction using physics-informed neural network and Nash-DDQN with BVP solver. $\overline{V_{1}}$ denotes mean final payoff to sub-swarm 1. Environment | $\overline{V_{1}}$ | $\overline{V_{2}}$ | Algorithm | $\overline{\hat{V}_{1}-V_{1}}$ | $\overline{\hat{V}_{2}-V_{2}}$
---|---|---|---|---|---
Case-1 (2 regions) | 0.697 | 0.697 | PINN | -0.011 | -0.012
| Nash
---
DDQN
-0.162 | -0.281
| Case-2
---
(4 regions)
0.150 | 0.147 | PINN | 0.0001 | -0.002

###  VI-C High Dimensional Case - 10 regions

Finally, to end the paper, we present a higher dimensional case study as shown
in Fig. 1(c). There are 10 regions connected by a directed graph with 10
edges. Each sub-swarm group has to pick a 10-dimensional action at each time
step. The action space of each sub-swarm is $2^{10}$, which makes the problem
intractable to solve in the case of Nash Q-learning and BVP. Solving BVP
becomes intractable as it requires computing guesses for each of the actions.

Training. We slightly increase the neural network size to 5 hidden layers
containing 128 neurons each. We keep other hyper-parameters the same as those
in the 4 regions case.

Figure 8: (a-d) Shows the density distribution of two sub-swarms at the final
time in the environment with 10 regions. Similar to the 4 regions case, the
sub-swarms are able to find a region where they dominate. The regions where
the sub-swarms dominate are bounded with a box.

Results. The final distributions of the sub-swarms across the 10 regions for 4
different initial conditions are shown in Fig. 8. From the plot, we see that
there exists at least one region for both sub-swarms where their proportions
are higher compared to the other sub-swarm. We do not observe a pronounced
difference in densities of the sub-swarms as in the case of 4 and 2 regions. A
possible explanation for not achieving this behavior could be an insufficient
time horizon. It is possible that the time horizon of 2.5 seconds is not
enough for sub-swarms to be able to aggregate in one region.

##  VII CONCLUSIONS

In this work, we demonstrated the efficacy of physics-informed neural networks
in solving general sum games between swarms. Unlike conventional solvers that
suffer from the curse of dimensionality as the state dimension increases,
PINNs can be easily scaled with the state dimension. However, an open question
still exists regarding the verification of the policies as a result of the
black-box nature of the neural network, especially for the higher dimensional
cases where we lack an analytical solution. In this work, we used simple case
studies to demonstrate the concept of solving HJI PDEs in a swarm setting.
Future works will explore applications to real-life inspired problems such as
network communications. Another possible direction is the addition of an
uncertainty term in the system that introduces a new layer of complexity. In
the examples discussed in the paper, an implied assumption is made about the
information availability. Strategies can get complicated if we introduce
incomplete information in the system, which is often the case in reality. A
natural extension of the present work is solving incomplete information games
on swarms.

## References

  * [1] M. Schranz, M. Umlauft, M. Sende, and W. Elmenreich, “Swarm robotic behaviors and current applications,” Frontiers in Robotics and AI, p. 36, 2020.
  * [2] S. Berman, A. Halász, M. A. Hsieh, and V. Kumar, “Optimized stochastic policies for task allocation in swarms of robots,” IEEE transactions on robotics, vol. 25, no. 4, pp. 927–937, 2009.
  * [3] A. Kushleyev, D. Mellinger, C. Powers, and V. Kumar, “Towards a swarm of agile micro quadrotors,” Autonomous Robots, vol. 35, no. 4, pp. 287–300, 2013\.
  * [4] D. Albani, T. Manoni, A. Arik, D. Nardi, and V. Trianni, “Field coverage for weed mapping: toward experiments with a uav swarm,” in Bio-inspired Information and Communication Technologies: 11th EAI International Conference, BICT 2019, Pittsburgh, PA, USA, March 13–14, 2019, Proceedings 11, pp. 132–146, Springer, 2019.
  * [5] W. Liu, A. F. Winfield, J. Sa, J. Chen, and L. Dou, “Towards energy optimization: Emergent task allocation in a swarm of foraging robots,” Adaptive behavior, vol. 15, no. 3, pp. 289–305, 2007.
  * [6] X. Ai, V. Srinivasan, and C.-K. Tham, “Optimality and complexity of pure nash equilibria in the coverage game,” IEEE Journal on Selected Areas in Communications, vol. 26, no. 7, pp. 1170–1182, 2008.
  * [7] E. Paraskevas, D. Maity, and J. S. Baras, “Distributed energy-aware mobile sensor coverage: A game theoretic approach,” in 2016 American Control Conference (ACC), pp. 6259–6264, IEEE, 2016.
  * [8] K. Elamvazhuthi and S. Berman, “Mean-field models in swarm robotics: A survey,” Bioinspiration & Biomimetics, vol. 15, no. 1, p. 015001, 2019\.
  * [9] C. Sinigaglia, A. Manzoni, and F. Braghin, “Density control of large-scale particles swarm through pde-constrained optimization,” IEEE Transactions on Robotics, vol. 38, no. 6, pp. 3530–3549, 2022.
  * [10] L. Zhang, M. Ghimire, W. Zhang, Z. Xu, and Y. Ren, “Approximating discontinuous nash equilibrial values of two-player general-sum differential games,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 3022–3028, 2023.
  * [11] M. G. Crandall and P.-L. Lions, “Viscosity solutions of hamilton-jacobi equations,” Transactions of the American mathematical society, vol. 277, no. 1, pp. 1–42, 1983.
  * [12] I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamilton-jacobi projections,” journal of Scientific Computing, vol. 19, pp. 323–346, 2003.
  * [13] Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang, “The expressive power of neural networks: A view from the width,” Advances in neural information processing systems, vol. 30, 2017.
  * [14] K. Elamvazhuthi, Z. Kakish, A. Shirsat, and S. Berman, “Controllability and stabilization for herding a robotic swarm using a leader: A mean-field approach,” IEEE Transactions on Robotics, vol. 37, no. 2, pp. 418–432, 2021.
  * [15] D. Milutinović and P. Lima, “Modeling and optimal centralized control of a large-size robotic population,” IEEE Transactions on Robotics, vol. 22, no. 6, pp. 1280–1285, 2006.
  * [16] F. Djeumou, Z. Xu, M. Cubuktepe, and U. Topcu, “Probabilistic control of heterogeneous swarms subject to graph temporal logic specifications: A decentralized and scalable approach,” IEEE Transactions on Automatic Control, vol. 68, no. 4, pp. 2245–2260, 2022.
  * [17] S. Bansal and C. Tomlin, “DeepReach: A deep learning approach to high-dimensional reachability,” in IEEE International Conference on Robotics and Automation (ICRA), 2021.
  * [18] J. Lu, Z. Shen, H. Yang, and S. Zhang, “Deep network approximation for smooth functions,” SIAM Journal on Mathematical Analysis, vol. 53, pp. 5465–5506, jan 2021.
  * [19] Z. Shen, H. Yang, and S. Zhang, “Neural network approximation: Three hidden layers are enough,” Neural Networks, vol. 141, pp. 160–173, sep 2021.
  * [20] D. P. Bertsekas, “Rollout algorithms for discrete optimization: A survey,”
  * [21] D. Bertsekas, “Multiagent reinforcement learning: Rollout and policy iteration,” IEEE/CAA Journal of Automatica Sinica, vol. 8, no. 2, pp. 249–272, 2021.
  * [22] J. Hu and M. P. Wellman, “Nash q-learning for general-sum stochastic games,” Journal of machine learning research, vol. 4, no. Nov, pp. 1039–1069, 2003\.
  * [23] P. Casgrain, B. Ning, and S. Jaimungal, “Deep q-learning for nash equilibria: Nash-dqn,” Applied Mathematical Finance, vol. 29, no. 1, pp. 62–78, 2022\.
  * [24] A. W. Starr and Y.-C. Ho, “Nonzero-sum differential games,” Journal of optimization theory and applications, vol. 3, pp. 184–206, 1969.
  * [25] T. Nakamura-Zimmerer, Q. Gong, and W. Kang, “Adaptive deep learning for high-dimensional hamilton–jacobi–bellman equations,” SIAM Journal on Scientific Computing, vol. 43, no. 2, pp. A1221–A1247, 2021.
  * [26] C. E. Lemke and J. T. Howson, Jr, “Equilibrium points of bimatrix games,” Journal of the Society for industrial and Applied Mathematics, vol. 12, no. 2, pp. 413–423, 1964.
  * [27] V. Knight, katiemcgoldrick, M. Panayides, Y. Wang, A. S. Gaba, tokheim, I. F. Jr., O. Konovalov, J. Campbell, N. Glynatsi, P. Rivière, R. Baldevia, R. Szeto, arwheel, newaijj, and volume-on max, “drvinceknight/nashpy: v0.0.40,” Aug. 2023.