# A Continuation Method for Large-Scale Modeling and Control: from ODEs to
PDE, a Round Trip

Denis Nikitin, Carlos Canudas-de-Wit and Paolo Frasca D. Nikitin, C. Canudas-
de-Wit and P. Frasca are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP,
GIPSA-lab, 38000 Grenoble, France.

###### Abstract

In this paper we present a continuation method which transforms spatially
distributed ODE systems into continuous PDE. We show that this continuation
can be performed both for linear and nonlinear systems, including
multidimensional, space- and time-varying systems. When applied to a large-
scale network, the continuation provides a PDE describing evolution of
continuous state approximation that respects the spatial structure of the
original ODE. Our method is illustrated by multiple examples including
transport equations, Kuramoto equations and heat diffusion equations. As a
main example, we perform the continuation of a Newtonian system of interacting
particles and obtain the Euler equations for compressible fluids, thereby
providing an original alternative solution to Hilbert’s 6th problem. Finally,
we leverage our derivation of the Euler equations to control multiagent
systems, designing a nonlinear control algorithm for robot formation based on
its continuous approximation.

###### Index Terms:

Control of Large-Scale Networks, Multiagent Systems, PDE

##  I Introduction

Most of the systems we encounter in real life consist of such a large number
of particles that the direct analysis of their interaction is impossible. In
such cases, simplified models are used that aggregate the behavior of a set of
particles and replace them with a continuous representation.

The first model describing a system of moving and interacting particles was
the Euler equation for liquids and gases which worked with a continuous field
[1]. This equation was one of the first Partial Differential Equations (PDEs).
Since then, many PDEs have been created that are models for describing
different physical processes, such as the Maxwell equation for electromagnetic
field [2].

With the advent of computers more and more attention has been paid to the
discretization, that is, the process of transforming PDE into a system of
ordinary differential equations in order to be able to numerically solve them
[3]. A correct discretization satisfies the condition that when the size of
the resulting ODE system increases, its solution will converge to the original
PDE. The choice of an apparent discretization method depends heavily on the
PDE type such that the solution would preserve the properties of PDE [4].

Although discretization is a widely known and widely used method, the inverse
problem of transforming an ODE system into PDE is more rarely studied. Our
work focuses on this particular problem, with the aim of filling this gap and
providing a counterpart to the discretization procedure. This can be useful
since PDEs provide a much more compact way of describing the system, which in
many cases is easier to analyze analytically than the corresponding ODE
system.

The idea of replacing the system with its compact and simplified
representation is widely used, especially for the ODE systems describing
large-scale networks. Probably the most known approach of this type is a model
reduction technique which transforms a network into a smaller one while
conserving the properties and the dynamics, see [5, 6]. Also the reduction
towards the average state was studied in [7, 8].

Apart from various model reduction techniques large-scale networks are studied
by mean field methods in case of the all-to-all interaction topology. In this
situation the effect of the network on each node is the same, therefore it is
enough to use an equation for a single agent together with parameters of a
state of the whole network, see [9] for a review with application to Kuramoto
networks.

The idea of mean field can be further extended to track not only a single
agent’s state, but the whole probability distribution over all agents’ states
in the network. This method is called population density approach [10] and is
mostly used to model large biological neural networks. In other fields of
physics it is described as a probability density of a physical system’s state
and is constructed by projection to probablistic space, see [11].

Large-scale systems can be also simplified by studying the approximations to
their probability densities, represented by moments. E.g., in [12, 13] a
moment-based approach was taken to control crowds dynamics. Different
applications and issues of the method of moments are covered in [14]. Also
shape parameters can be used to simplify the model and describe a shape of the
solution as in [15].

Mean field and population density approaches are suitable in the case when the
interaction topology between nodes is all-to-all. In other cases, the
continuous representation of a system requires more sophisticated tools. A
recently emerged theory of graphons studies graph limits, i.e. structural
properties that the graph possesses if the number of nodes tends to infinity
while preserving interaction topology, see [16]. Using graphons it is possible
to describe any dense graph as a linear operator in continuum space [17]. This
method was further used to control large-scale linear networks [18] and to
study sensitivity of epidemic networks [19].

However, we are interested in systems that are spatially distributed and which
have a position-dependent interaction, such as traffic in the city, power
networks, robot formations, etc. By applying population density method or
graphon theory to such system we would end up with a continuous model which
looses the spatial structure of the problem.

Our idea is to replace the original spatially distributed ODE system by a
continuous PDE whose state and space variables preserve the state and space
variables of the original system. We develop a method for linear spatially
invariant ODEs which transform them into PDEs with the help of finite
differences. We name this method as a continuation, since it is exactly
opposite to the discretization procedure. Further we show how the continuation
converges to the original system in sense of spectrum. Using computational
graph formalism [20] we extend the method to nonlinear systems and further to
space-dependent systems and systems with boundaries.

The continuation method allows to recover a PDE which describes the same
physical system as the original ODE network. Such a description can be very
helpful both for analysis and control design purposes. Indeed, one can use an
obtained PDE to design a continuous control which, being discretized back,
results in a control law for the original ODE system: this design framework is
illustrated in Fig. 1.

ODEcontinuationPDEcontroldesigncontrolled PDEdiscretizationcontrolled ODE
Figure 1: Proposed framework for control design based on the continuation
method and a continuous representation of the system.

Moving to the examples of the application of the continuation method, we
tackle the Hilbert’s 6th problem, questioning how one can rigorously transform
a system of interacting particles into the Euler PDE. We provide our treatment
of this problem using continuation, deriving the Euler PDE from a system of
Newton’s laws for the case of long-range interaction forces. Further, we show
how the method can be applied to the multiagent control, providing a simple
control algorithm to stabilize a robotic formation along the desired
trajectory, performing a maneuver of passing through a window. The control is
derived on a level of a PDE representation and then it is discretized to be
implemented on every agent in accordance with the scheme in Fig. 1.

We start in Section II by defining a continuation for linear ODEs, discussing
questions of accuracy, convergence and choice of the particular model. Section
III continues to nonlinear models, utilizing the computational graph
formalism. In Section IV the method is extended to much broader class of
systems, including multidimensional or space- and time-varying systems and
also discussing boundary conditions. Section V is devoted to a derivation of
the Euler PDE from a system of interacting particles. Finally, Section VI
applies the method to control a robotic formation.

##  II Method for linear systems

The simplest class of systems for which the transformation of ODE into PDE can
be performed is given by linear ODE systems corresponding to the dynamics of
states of nodes, which are aligned on the 1D line in space and depend only on
some fixed set of their neighbours.

$\rho_{i}$$\rho_{i+s_{2}}$$\rho_{i+s_{1}}$$\rho_{i+s_{j}}$$a_{2}$$a_{1}$$a_{j}$$x$$x_{i-1}$$x_{i}$$x_{i+1}$$x_{i+s_{j}}$$\Delta
x$$\Delta x$ Figure 2: System of nodes aligned in 1D line with dynamics given
by (1) with $s_{1}=-1$ and $s_{2}=1$.

Let the node $i$ have a state $\rho_{i}\in\mathbb{R}$ and a geographical
position $x_{i}\in\mathbb{R}$ such that for every $i$ the distance between two
consecutive nodes in space is constant, $x_{i+1}-x_{i}=\Delta x$ (the
assumption of $\Delta x$ being constant will be relaxed later on). Then the
systems of our interest take the form

$\dot{\rho}_{i}=\sum\limits_{j=1}^{N}a_{j}\rho_{i+s_{j}},$ (1)

that is $\dot{\rho}_{i}$ linearly depends only on $N$ neighbouring nodes
shifted by $s_{j}\in\mathbb{Z}$ for $j\in\\{1..N\\}$, and $a_{j}\in\mathbb{R}$
are the system gains, see Fig. 2. This type of systems belongs to the class of
linear spatially invariant systems [21], which is a natural class for
distributed control.

###  II-A Motivating Example

We start by considering the most simple ODE system of class (1) which has
spatial dependence:

$\dot{\rho}_{i}=\frac{1}{\Delta x}\left(\rho_{i+1}-\rho_{i}\right).$ (2)

Comparing with (1), here $N=2$, $a_{1}=1/\Delta x$, $a_{2}=-1/\Delta x$,
$s_{1}=1$ and $s_{2}=0$. This equation describes a transport of some quantity
along the line, and is usually referred as a Transport ODE. Equation (2) is
sometimes studied on its own, but often it comes as a result of a
discretization process applied to another equation,

$\frac{\partial\rho}{\partial t}=\frac{\partial\rho}{\partial x}.$ (3)

This equation belongs to a class of PDEs, which is usually thought to be more
difficult class of equations to study than ODEs. However, equation (3)
describes a perfect transport of information with finite propagation speed
along the line, which can be studied much more easily in PDE form than in ODE,
as it perfectly conserves the form of a solution, performing only a shift
along the line as time increases. We will refer to this equation as a
Transport PDE.

Equation (2) can be obtained from (3) by the discretization process, which has
been a well-established mathematical tool for centuries. Nevertheless, up to
now there was no strict procedure describing a general process which could
render equation (3) from (2). In the next subsections we explore more how the
general discretization procedure is defined for linear systems and how it
should be inverted to obtain a continuation process.

###  II-B Discretization

The discretization of PDEs is usually performed by a finite difference method,
where the partial derivatives are approximated by finite differences. For
example, in the case of Transport ODE,

$\frac{\partial\rho}{\partial x}\approx\frac{1}{\Delta
x}\left(\rho_{i+1}-\rho_{i}\right).$

This approximation is valid when $\Delta x$ is small. Indeed, assuming that
the solution to PDE is given by a smooth function $\rho(x)$ and using Taylor
series, we can write

$\rho_{i+1}=\rho(x_{i+1})=\rho(x_{i})+\frac{\partial\rho}{\partial x}\Delta
x+\frac{\partial^{2}\rho}{\partial x^{2}}\frac{\Delta x^{2}}{2}+...,$ (4)

where all partial derivatives are calculated in $x_{i}$. Thus, subtracting
$\rho_{i}$ and dividing by $\Delta x$, we get

$\frac{\partial\rho}{\partial x}=\left[\frac{1}{\Delta
x}\left(\rho_{i+1}-\rho_{i}\right)\right]-\frac{\partial^{2}\rho}{\partial
x^{2}}\frac{\Delta x}{2}-...,$ (5)

which means that the residual belongs to the class $O(\Delta x)$ of all
functions which go to zero at least as fast as $\Delta x$. Thus, taking
$\Delta x$ sufficiently small, one can ensure the arbitrary accuracy of the
approximation, provided all the partial derivatives are bounded.

Accuracy can be further increased by taking different points where the
function is sampled, called stencil points. For example, writing

$\rho_{i-1}=\rho(x_{i-1})=\rho(x_{i})-\frac{\partial\rho}{\partial x}\Delta
x+\frac{\partial^{2}\rho}{\partial x^{2}}\frac{\Delta x^{2}}{2}-...,$ (6)

subtracting (6) from (4) and dividing by $2\Delta x$, we get

$\frac{\partial\rho}{\partial x}=\left[\frac{1}{2\Delta
x}\left(\rho_{i+1}-\rho_{i-1}\right)\right]-\frac{\partial^{3}\rho}{\partial
x^{3}}\frac{\Delta x^{2}}{6}+....$ (7)

Thus, using stencil points $\\{i-1,i+1\\}$ to approximate the first-order
derivative in the point $i$ the obtained residual belongs to the class
$O(\Delta x^{2})$, which means that this discretization of the Transport PDE
is accurate to the second order.

In general, if one wants to approximate the derivative of order $m$ in point
$i$ using $N$ stencil points $\\{i+s_{1},i+s_{2},...,i+s_{N}\\}$ with $m<N$ in
form

$\frac{\partial^{m}\rho}{\partial
x^{m}}\approx\sum\limits_{j=1}^{N}\hat{a}_{j}\rho_{i+s_{j}}$ (8)

where coefficients $\hat{a}_{j}$ are unknown, one can define
$S_{N,N}\in\mathbb{R}^{N\times N}$, $\hat{a}\in\mathbb{R}^{N}$ and
$c\in\mathbb{R}^{N}$ by

$S_{N,N}=\begin{pmatrix}1&\cdots&1\\\ s_{1}&\cdots&s_{N}\\\
\vdots&\ddots&\vdots\\\
s_{1}^{N-1}&\cdots&s_{N}^{N-1}\end{pmatrix},\;\hat{a}=\begin{pmatrix}\hat{a}_{1}\\\
\hat{a}_{2}\\\ \vdots\\\ \hat{a}_{N}\end{pmatrix},\;c=\frac{m!}{\Delta
x^{m}}\begin{pmatrix}0\\\ \vdots\\\ 1\\\ \vdots\\\ 0\end{pmatrix},$

where $c$ has $1$ on the position $m+1$, and solve a linear system

$\hat{a}=S_{N,N}^{-1}\>c.$ (9)

The system (9) can be trivially obtained by writing Taylor series for all
points $\rho_{i+s_{1}}...\rho_{i+s_{N}}$ and summing them in a linear
combination as in (8). The obtained order of accuracy is at least $O(\Delta
x^{(N-m)})$, and sometimes can be higher if some of the higher derivatives are
also eliminated (as in case of (7)).

###  II-C Continuation

Essentially the same process can be applied to the equation (1) to get the PDE
version. For every term in the summation in (1) we can write

$\rho_{i+s_{j}}=\rho(x_{i+s_{j}})=\rho(x_{i})+\frac{\partial\rho}{\partial
x}\Delta xs_{j}+\frac{\partial^{2}\rho}{\partial x^{2}}\frac{\Delta
x^{2}s_{j}^{2}}{2}+...$ (10)

Thus, assume we state the problem of finding the PDE approximation of (1) in
form

$\sum\limits_{j=1}^{N}a_{j}\rho_{i+s_{j}}\approx\sum\limits_{k=0}^{d}c_{k}\frac{\Delta
x^{k}}{k!}\frac{\partial^{k}\rho}{\partial x^{k}},$ (11)

where $d$ is the highest order of derivative (order of continuation) we want
to use. Note that zero is also included in the right summation, since the
function itself can be used in the resulting PDE. Then, introducing
$S_{d+1,N}\in\mathbb{R}^{(d+1)\times N}$, $a\in\mathbb{R}^{N}$ and
$c\in\mathbb{R}^{d+1}$ by

$S_{d+1,N}=\begin{pmatrix}1&\cdots&1\\\ s_{1}&\cdots&s_{N}\\\
\vdots&\ddots&\vdots\\\
s_{1}^{d}&\cdots&s_{N}^{d}\end{pmatrix},\;\;\;a=\begin{pmatrix}a_{1}\\\
a_{2}\\\ \vdots\\\ a_{N}\end{pmatrix},\;\;\;c=\begin{pmatrix}c_{0}\\\ c_{1}\\\
\vdots\\\ c_{d}\end{pmatrix},$ (12)

the vector of unknown coefficients $c$ can be found by direct multiplication,

$c=S_{d+1,N}\;a,\quad\text{or}\quad
c_{k}=\sum\limits_{j=1}^{N}a_{j}s_{j}^{k}\quad\forall k\in\\{0,..,d\\}.$ (13)

Once (13) is solved, we write the PDE approximation to (1):

$\frac{\partial\rho}{\partial t}=\sum\limits_{k=0}^{d}c_{k}\frac{\Delta
x^{k}}{k!}\frac{\partial^{k}\rho}{\partial x^{k}}.$ (14)

As an example, applying (13) to the Transport ODE (2) renders coefficients
$c_{0}=0$ and $c_{k}=1/\Delta x$ for all $k>0$, and choosing $d=1$ we obtain
the Transport PDE (3).

###  II-D Accuracy of continuation

Procedures (9) and (13) look very similar from the algebraic point of view,
however they are qualitatively different in the way how the problem is
formulated and how we should interpret their results.

The discretization procedure tries to find the best approximation to a
continuous and smooth function $\rho$ and its derivatives. What is most
important, the discretization step $\Delta x$ is usually an adjustable
parameter which can be set by a system engineer arbitrarily small to satisfy
the desired performance. Thus the notion of accuracy of a discretization is
used to describe how fast the solution of the discretized equation tends to
the solution of the original equation when $\Delta x$ tends to zero. In some
sense this means quality of the discretization, since the higher order of
accuracy means that the engineer can take a larger $\Delta x$ to achieve the
same error and thus use the smaller number of states in the discretized
system.

Instead, when the original system is given by the ODE, the nodes have fixed
locations, thus $\Delta x$ is a true constant representing properties of an
underlying physical system and it cannot be changed by an engineer. In turn
this means that the accuracy defined as a class $O(\Delta x^{(N-m)})$ cannot
measure quality of the approximation as $\Delta x$ does not behave as an
arbitrarily small value. Moreover, in the ODE case the system state $\rho_{i}$
is known only on a given set of points $i$, thus in general the approximation
$\rho(x)$ can be non-smooth or discontinuous. Even if we assume the smoothness
at the initial moment of time, the dynamics can render its derivatives
unbounded. As a result the series in (10) can be non-convergent.

We know however that the systems (1) and (14) are connected by finite
difference methods (9) and (13). We will use this fact as a definition of a
PDE approximation to an ODE and of an accuracy of this approximation.

######  Definition 1.

Discretization of PDE to ODE is called valid if it is performed according to
the finite difference method (9).

######  Definition 2.

Continuation of ODE to PDE is called valid if there exists a valid
discretization of the obtained PDE to the original ODE.

Definitions 1 and 2 say that if we use continuation on some ODE and then
perform discretization at the same stencil points, we should arrive at the
same ODE. This procedure sets a constraint on the minimum order of the PDE:

######  Theorem 1.

The PDE (14) that is obtained from the ODE (1) with $N$ stencil points is
valid if and only if

$d+1\geqslant N.$ (15)

Proof of necessity:  Indeed, assume $d+1<N$. This means that the PDE
approximation (14) has the highest derivative at most of the order $d$. Thus
we can augment the vector of coefficients $c$ by $N-d-1$ zeros corresponding
to higher-order derivatives, which is equivalent to the augmentation of
$S_{d+1,N}$ with $N-d-1$ zero rows since $c$ was defined by (13). Now,
applying the discretization process (9) to $c$ we should arrive at the same
vector $a$ of the parameters of the ODE system. Since this should be true for
any $a$, we substitute $S_{N,N}^{-1}$ and augmented $S_{d+1,N}$ and obtain a
condition

$\footnotesize{\begin{pmatrix}1&1&\cdots&1\\\
s_{1}&s_{2}&\cdots&s_{N}\\\\[-2.84526pt] \vdots&\vdots&\ddots&\vdots\\\
s_{1}^{N-1}&s_{2}^{N-1}&\cdots&s_{N}^{N-1}\end{pmatrix}^{-1}\begin{pmatrix}1&1&\cdots&1\\\
s_{1}&s_{2}&\cdots&s_{N}\\\\[-5.69054pt]
\vdots&\vdots&\ddots&\vdots\\\\[-1.42262pt]
s_{1}^{d}&s_{2}^{d}&\cdots&s_{N}^{d}\\\ 0&0&\cdots&0\\\\[-5.69054pt]
\vdots&\vdots&\ddots&\vdots\end{pmatrix}=I,}$

which is impossible to satisfy since the second matrix is singular. Therefore
there is no valid discretization process for the PDE obtained by continuation
with order $d$ such that $d+1<N$, which by definition means that such
continuation is not valid.

Proof of sufficiency:  Case $d+1=N$ is trivial, since the equations (9) and
(13) are equivalent in this situation. This obviously provides a validity of
the continuation procedure.

Now assume $d+1>N$. Then the obtained vector of coefficients $c$ is of higher
dimension than $a$. The discretization (9) cannot be applied directly, since
there is not enough stencil points to express the finite differences for the
derivatives of order higher than $N-1$. However we can increase the set of
stencil points. Let us choose additional $d+1-N$ stencil points
$\bar{s}_{N+1},\bar{s}_{N+2},...,\bar{s}_{d+1}$. Applying continuation (13) to
the original ODE (1) and then (9) to the obtained PDE using the augmented set
of stencil points we get a new ODE gains $\bar{a}$ which are expressed as

$\footnotesize{\bar{a}=\begin{pmatrix}1&\cdots&1&1&\cdots&1\\\
s_{1}&\cdots&s_{N}&\bar{s}_{N+1}&\cdots&\bar{s}_{d+1}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
s_{1}^{d}&\cdots&s_{N}^{d}&\bar{s}_{N+1}^{d}&\cdots&\bar{s}_{d+1}^{d}\end{pmatrix}^{-1}\begin{pmatrix}1&\cdots&1\\\
s_{1}&\cdots&s_{N}\\\ \vdots&\ddots&\vdots\\\
s_{1}^{d}&\cdots&s_{N}^{d}\end{pmatrix}a.}$

It is now clear that the first $N$ elements of $\bar{a}$ are exactly $a$ and
the rest is zero, irrespective of the chosen additional points $\bar{s}_{j}$.
Thus the artificially introduced stencil points do not appear in the
discretized PDE, rendering the same ODE as the original one. ∎

The latter part of the proof of Theorem 1 means that the PDE obtained by the
process (13) with $d+1>N$ has more information that one with $d+1=N$, since it
provides exact Taylor approximations not only on the given set of points, but
in the additional $d+1-N$ points which can be chosen arbitrary. This property
can be used to define the order of accuracy as $d+1-N$.

######  Definition 3.

Order of accuracy of a valid continuation process of ODE to PDE is defined as
the number of additional points in which the corresponding discretization
process can be made exact simultaneously, i.e. $d+1-N$.

For example, a continuation

$\rho_{i+1}-\rho_{i}\quad\rightarrow\quad\Delta x\frac{\partial\rho}{\partial
x}$

is of order of accuracy 0, since trying to discretize the PDE on any larger
set of stencil points except from $\\{i,i+1\\}$ will give different ODE. At
the same time a continuation

$\rho_{i+1}-\rho_{i-1}\quad\rightarrow\quad 2\Delta
x\frac{\partial\rho}{\partial x}+\left(0\cdot\frac{\partial^{2}\rho}{\partial
x^{2}}\right)$

has order of accuracy 1 because the second derivative vanishes (thus $d=2$),
and it is possible to discretize the PDE on a set of stencil points of size 3
(with one additional point), for example $\\{i-1,i,i+1\\}$. From now on we
will consider only valid continuations, and most often in practice we will use
continuations of order of accuracy $0$.

###  II-E Convergence of continuation

It is clear that the higher order of continuation is taken, the better the
original ODE operator (1) is approximated by the PDE (14). It is possible to
study the convergence properties by shifting the problem to the frequency
domain using the Fourier transform.

Let us define the function $a(x)$ as

$a(x)=\sum\limits_{j=1}^{N}a_{j}\delta\left(x+s_{j}\Delta x\right),$

where $\delta(x)$ is the Dirac delta function. Therefore, for any integrable
function $\rho(x)$ the equation (1) is equivalent to the following system with
convolution

$\dot{\rho}(x)=(a\star\rho)(x).$ (16)

Use now the Fourier transform, defined as

$\hat{f}(\omega)=\int_{-\infty}^{\infty}f(x)e^{-ix\omega}dx$ (17)

for any integrable function $f(x)$ and for any frequency
$\omega\in\mathbb{R}$. It is known that the Fourier image of a convolution is
a multiplication. Therefore the system (16) is just

$\dot{\hat{\rho}}(\omega)=\hat{a}(\omega)\hat{\rho}(\omega),\qquad\hat{a}(\omega)=\sum\limits_{j=1}^{N}a_{j}e^{is_{j}\Delta
x\omega},$ (18)

where $\hat{a}(\omega)$ was found by direct calculation of Fourier transform.
(18) immediately shows that the spectrum of the original ODE system is
parametrized by $\hat{a}(\omega)$. In fact, this result is well-known, since
the system (1) on an infinite line belongs to the class of Laurent systems,
whose spectrum is known to be (18), see [22].

Now let us calculate the spectrum of the continualized system (14). By another
property of the Fourier transform, if the function $\rho(x)$ is sufficiently
smooth and its derivatives are integrable, we can recover their Fourier images
by

$\widehat{\left(\frac{\partial^{k}\rho}{\partial
x^{k}}\right)}\left(\omega\right)=\left(i\omega\right)^{k}\hat{\rho}(\omega).$

Therefore (14) is read in frequency domain as

$\dot{\hat{\rho}}(\omega)=\hat{c}(\omega)\hat{\rho}(\omega),\qquad\hat{c}(\omega)=\sum\limits_{k=0}^{d}c_{k}\frac{\Delta
x^{k}}{k!}\left(i\omega\right)^{k}.$ (19)

Substituting (13), we can rewrite (19) as

$\hat{c}(\omega)=\sum\limits_{j=1}^{N}a_{j}\sum\limits_{k=0}^{d}\frac{(is_{j}\Delta
x\omega)^{k}}{k!}.$ (20)

Now, comparing (20) with (18), it is clear that (20) uses the first $d+1$
terms of the Taylor expansion of the exponential function in (18). In fact,
since the exponential function is analytic on the whole complex plane, we have
just proven the following result:

######  Theorem 2.

The spectrum of the PDE (14) converges to the spectrum of the original ODE (1)
pointwise as $d\to\infty$.

However, the convergence of spectrums is not uniform. Moreover, the spectrum
(18) is an image of the unit circle and thus is a compact set, while the
spectrum (20) for any $d$ is a polynomial and thus unbounded. It is still
possible to state several corollaries of Theorem 2:

######  Corollary 1.

For any bounded subset of frequencies the spectrum of the PDE (20) converges
to the spectrum of the ODE (18) uniformly.

###### Proof.

The result follows directly from Theorem 2 and the fact that the Taylor
expansion of any analytic function converges to the function itself uniformly
on any bounded domain. ∎

Figure 3: Spectrum for the Transport ODE (21) $e^{i\omega}-1$ (blue circle)
together with spectrums of the continuations up to the order 6, including
(22). As $d$ increases, spectrums converge to the blue circle, however for
some orders (such as 4 or 5) they can become unstable.

######  Corollary 2.

If the original ODE (1) is unstable, there exists $D\geqslant 0$ such that for
all $d\geqslant D$ the continualized system (14) will also be unstable.

###### Proof.

Since the original system is unstable, there exists $\omega_{0}$ such that
$\operatorname{Re}\hat{a}(\omega_{0})>0$. Now, by the definition of a limit
there exists $D\geqslant 0$ such that for all $d\geqslant D$
$\operatorname{Re}\hat{c}(\omega_{0})>0$. ∎

Unfortunately, the stability of the original system cannot be always recovered
by taking large enough $d$. The reason for this is in the unboundedness of the
polynomial spectrum (20). Even if the original system is stable, for some $k$
it is possible to introduce artificial instability on high frequencies. There
are however some guidelines which can be used to properly choose the order of
continuation $d$:

######  Corollary 3.

PDE (14) with an odd order of continuation $d$ has the same stability
properties as a PDE with the order of continuation $d-1$.

###### Proof.

All odd terms in the spectrum (20) are purely imaginary and thus have no
impact on the stability. ∎

######  Corollary 4.

Artificial instability is introduced when the last even term in the PDE (14)
has $c_{k}>0$ if $k=4m$ or $c_{k}<0$ if $k=4m+2$ for some
$m\in\mathbb{Z}^{+}$.

###### Proof.

Artificial instability comes if the term of the polynomial (20) with the
highest even power is positive, which leads to a positive real part of a
spectrum on high frequencies. Positivity of the highest even term is exactly
equivalent to the statement of the corollary since $i^{4m}=1$ and
$i^{4m+2}=-1$ for any $m\in\mathbb{Z}^{+}$. ∎

We will demonstrate the convergence of spectrums on the Transport ODE

$\dot{\rho}_{i}=\rho_{i+1}-\rho_{i}.$ (21)

Assuming $\Delta x=1$, the continuation of (21) is:

$\frac{\partial\rho}{\partial
t}=\sum\limits_{k=1}^{d}\frac{1}{k!}\frac{\partial^{k}\rho}{\partial x^{k}},$
(22)

Spectrum of (21) by formula (18) is given by $e^{i\omega}-1$, which is
depicted as a blue circle in Fig. 3 together with the spectrums of the
continuations up to the order $d=6$. It is clear that as the order increases,
the approximations become better.

The original Transport ODE is stable. Moreover, it has an intrinsic diffusion
in it, which can be captured by the continuation of the second order. However,
the continuation of order 4 is unstable. It happens because of an artificial
instability as described in Corollary 4, since $c_{4}=1>0$. In general all
stable continuations of the Transport ODE are given by the orders
$\\{1,2,3,...,4m+2,4m+3,...\\}$ for all $m\in\mathbb{Z}^{+}$.

##  III Method for nonlinear systems

Finite differences give us a complete tool for linear systems, but for
nonlinear systems they should be applied in composition with nonlinearities.
Using an additional concept of computational graph it is possible to elaborate
the case of general nonlinear ODE systems.

As in the previous case we assume without loss of generality that the nodes
are equally spaced along the 1D line, a node $i$ having a state $\rho_{i}$ and
a position $x_{i}$. Then the general nonlinear ODE with space dependence takes
form of

$\dot{\rho}_{i}=F(\rho_{i+s_{1}},\rho_{i+s_{2}},...,\rho_{i+s_{N}}).$ (23)

We further assume that the function $F$ is continuous.

###  III-A Computational graph

In 1957 Kholmogorov [23] showed that every multidimensional continuous
function can be written as a composition of functions of one variable and
additions. This work laid the basis for the neural networks function
approximation, which is now a major branch of modern machine learning.

Here we will use this idea and assume that the function $F$ is given in the
form of computational graph (see [20] for example). This is a directed acyclic
graph, every node of which represents a one-dimensional function, applied to a
weighted sum of inputs coming to this node. We assume that the leaves of this
graph are the states of the system $\rho_{i+s_{j}}$ and the root node computes
the resulting value of $F$.

As an example of the computation graph we will consider a system

$\dot{\rho}_{i}=\sin(\rho_{i+1}-\rho_{i})-\sin(\rho_{i}-\rho_{i-1})$ (24)

which is a system of Kuramoto oscillators coupled on a ring. The computational
graph for (24) is presented in Fig. 4.

resultsinsin$\rho_{i-1}$$\rho_{i}$$\rho_{i}$$\rho_{i+1}$-11-11-11$i-\frac{1}{2}$$i+\frac{1}{2}$
Figure 4: Computation graph for the system (24). Similar subgraphs are
outlined by dashed rectangles of the same color. Possible choices of sinus
subgraph’ positions are written in the corners of blue rectangles.

###  III-B Similar subgraphs and their positions

Now let us introduce an original notion of similar subgraphs. Subgraph is a
computational graph which computes subexpression of the original computational
graph. Every node in a computational graph serves as the root of a subgraph
computing expression defined in this node. The leaf nodes are also the
subgraphs ”computing” themselves.

######  Definition 4.

We call two subgraphs similar if

  1. 1.

they serve as an input to the same node,

  2. 2.

they differ only in the positions of the leaf nodes, and this difference can
be represented by a single shift.

This is an equivalence relation, therefore we can speak about equivalence
classes which we call sets of similar subgraphs.

For example, in Fig. 4 there are three sets of similar subgraphs:

  1. 1.

$\rho_{i-1}$ and $\rho_{i}$ for the left sinus node,

  2. 2.

$\rho_{i}$ and $\rho_{i+1}$ for the right sinus node,

  3. 3.

$\sin(\rho_{i}-\rho_{i-1})$ and $\sin(\rho_{i+1}-\rho_{i})$ for the root node,
because they differ by a single shift which equals 1.

The last thing which should be defined is a position of a subgraph:

######  Definition 5.

Position of a subgraph is defined as a coordinate in space where the
expression of this subgraph is calculated.

The leaf nodes by definition are the states of the system, thus they are
calculated at some positions on the line. For example the leaf node
$\rho_{i+1}$ in Fig. 4 is defined in the point $x_{i+1}$, thus we will say
that its position is $i+1$.

The root node by definition has a position $i$, since it is exactly the
position of the left-hand side term in (23). We will define the positions of
other subgraphs as the average of their leaves’ positions. Note that in
general there is some freedom in the definition of the subgraphs’ positions,
with the only constraint that similar subgraphs should differ by a single
shift, but we will omit this for simplicity.

Since the position of a subgraph represents a position on the line, it is
natural to have non-integer position values, although the leaf nodes and the
root have only integer positions. As an example, defining the position as an
average, in Fig. 4 the node $\sin(\rho_{i+1}-\rho_{i})$ has its position
$i+1/2$.

###  III-C Continuation to a nonlinear PDE

When system (23) is expressed in a form of computational graph with similar
subgraphs being found and their positions being defined, one can perform a
continuation procedure described in section II-C to obtain a PDE.

Continuation should be performed recursively, starting from the leaves. Each
set of similar subgraphs by definition is used in their common ancestor node
as a linear combination of equivalent elements shifted by some distance.
Continuation of this linear combination by (11) replaces the set of similar
subgraphs by a weighted sum of partial derivatives of subexpressions,
calculated at the position of the ancestor node.

Let $\Delta x=x_{i+1}-x_{i}$ be a distance between any two neighbouring nodes.
Elaborating example (24) and using order of accuracy 0, we perform the
continuation in three steps:

  1. 1.

$\sin(\rho_{i+1}-\rho_{i})\quad\rightarrow\quad\sin\left(\Delta
x\dfrac{\partial\rho}{\partial x}(x_{i+1/2})\right)$,

  2. 2.

$\sin(\rho_{i}-\rho_{i-1})\quad\rightarrow\quad\sin\left(\Delta
x\dfrac{\partial\rho}{\partial x}(x_{i-1/2})\right)$,

  3. 3.

$\sin_{i+1/2}-\sin_{i-1/2}\quad\rightarrow\quad\Delta
x\dfrac{\partial}{\partial x}\sin$.

which finally gives a nonlinear PDE representation of (24):

$\frac{\partial\rho}{\partial t}=\Delta x\dfrac{\partial}{\partial
x}\sin\left(\Delta x\dfrac{\partial\rho}{\partial x}\right).$ (25)

Speaking about the accuracy, the simplest possible PDE can be obtained by
choosing orders of accuracy equal to 0 for each set of similar subgraphs. For
more accurate equations it makes sense to specify the desired order of the
equation $d$ and then get rid of all the terms which consist of composition of
derivatives of combined order higher than $d$.

##  IV Extensions

Until now we discussed systems with nodes which were uniformly placed on the
infinite 1D line and which had common space-independent dynamics. The method
can be extended to include more classes of systems.

Spatially invariant systems [21] such as periodic ones can be tackled by
choosing different index spaces. In the periodic case we can assume that the
positions $x\in\mathbb{S}$ are placed on the unit circle and indices
$i\in\mathbb{Z}/n\mathbb{Z}$ form a ring of integers modulo $n$, where $n$ is
the number of states of the original ODE. Since any function on $\mathbb{S}$
can be mapped to a periodic function on $\mathbb{R}$, the analysis in Sections
II and III remain the same.

Further the time dependence can be introduced into system gains both in the
ODE and in the PDE, where the continuation is performed independently for
every fixed $t$. This allows to use the method for time-varying systems and
switching networks. Also systems whose state is vector-valued can be
continualized using the same finite differences based scheme, thus in the
following we will assume that the state of a system is scalar.

In the following subsections we will explore how the method can be extended to
include systems with several spatial dimensions, systems with space dependence
or nonuniform placing and systems with boundaries. Finally we introduce a
concept of PDE with index derivatives which can be applied to systems whose
states coincide with the positions in space, for example particle systems.

###  IV-A Multidimensional systems

Let a position of a node $\rho_{i}$ be described by $x_{i}\in\mathbb{R}^{n}$.
Moreover, a node $\rho_{i}$ is referenced by a multi-index
$i=(i_{1},...,i_{n})\in\mathbb{Z}^{n}$. We assume that the position difference
between two neighbour nodes $i=(i_{1},...,i_{k},...,i_{n})$ and
$i^{\prime}=(i_{1},...,i_{k}+1,...,i_{n})$ is

$x_{i^{\prime}}-x_{i}=(0,...,\Delta x_{k},...,0)\qquad\forall k\in\\{1..n\\},$

and that there exists a vector $\Delta x=(\Delta x_{1},...,\Delta
x_{k},...,\Delta x_{n}).$

We shall describe the extension for linear systems, as nonlinear systems can
be dealt with in an analogous way by the same computational graph concept as
in Section III. Thus we start from the linear system (1). For nonnegative
multi-index $h$ we define an absolute value $|h|=\sum_{k=1}^{n}h_{k}$.
Further, we define multi-index power $h$ of a vector $x$ as
$x^{h}=\prod_{k=1}^{n}x_{k}^{h_{k}}$, with an assumption $0^{0}=1$. By Taylor
series

$\rho_{i+s_{j}}=\rho_{i}+\sum\limits_{|h|=1}s_{j}^{h}\Delta
x^{h}\frac{\partial\rho_{i}}{\partial
x^{h}}+\sum\limits_{|h|=2}\frac{s_{j}^{h}}{2}\Delta
x^{h}\frac{\partial^{2}\rho_{i}}{\partial x^{h}}+...,$ (26)

and the continuation up to the order $d$ thus is

$\frac{\partial\rho}{\partial
t}=\sum\limits_{\begin{subarray}{c}h\in\mathbb{Z}^{n}_{+}\\\ |h|\leqslant
d\end{subarray}}c_{h}\frac{\Delta
x^{h}}{|h|!}\frac{\partial^{|h|}\rho}{\partial x^{h}},\qquad
c_{h}=\sum\limits_{j=1}^{N}a_{j}s_{j}^{h}.$ (27)

Using $\mathbb{S}$ instead of $\mathbb{R}$, more complex multidimensional
spaces can be covered such as torus or cylinder.

###  IV-B Space-dependent systems

Let us now look on the linear system (1) with one important difference: the
system gains $a_{j}$, the shifts $s_{j}$ and the number of neighbours $N$
become space-dependent:

$\dot{\rho}_{i}=\sum\limits_{j=1}^{N_{i}}a_{ij}\rho_{i+s_{ij}}.$ (28)

Notice that equation (28) describes in fact any linear system.

Now, choosing a unique $d$ such that $d+1\geqslant N_{i}$ for all $i$, one can
perform a continuation (13) at every point $x_{i}$ up to the order $d$ and
obtain a PDE (14) with space dependent gains $c_{ik}$. This means that we know
the gains $c_{ik}$ at the points with coordinates $x_{i}$, which can be seen
as a sampling of some function $c_{k}(x)$ at points $x_{i}$.

We can now perform either an interpolation or an approximation based on this
sampling. In the first case we seek for $c_{k}(x)$ such that
$c_{k}(x_{i})=c_{ik}$, while in the second case it is enough to satisfy this
relation approximately, for example by determining $c_{k}(x)$ via least
squares. In either case, the resulting continuation of (28) is given by

$\frac{\partial\rho}{\partial t}=\sum\limits_{k=1}^{d}c_{k}(x)\frac{\Delta
x^{k}}{k!}\frac{\partial^{k}\rho}{\partial x^{k}}.$ (29)

In the case of nonlinear systems the continuation can be performed if the
computational graphs for every node compute the same dynamics. We can
formalize it by stating the following property:

######  Definition 6.

We say that two computational graphs have the same structure if

  1. 1.

their root nodes compute the same expression,

  2. 2.

any child subgraph of the root node of the first graph has the same structure
with some child subgraph of the root node of the second graph and vice versa.

This definition, formulated through recursion, essentially means that the
order of nonlinearities which is hidden in two computational graphs should
coincide, see Fig. 5.

$F(\cdot)$$F(\cdot)$$G(\cdot)$$H(\cdot)$$H(\cdot)$$G(\cdot)$$H(\cdot)$$a_{11}$$a_{12}$$a_{13}$$a_{21}$$a_{22}$
Figure 5: Illustration of two computational graphs having the same structure.

Finally, a continuation of a nonlinear ODE system can be performed if all the
computational graphs computing the dynamics for all states $\rho_{i}$ have the
same structure. Indeed, in this case it is possible to perform a continuation
for any set of similar subgraphs for each node as in the linear case of
(28)-(29). Moreover, by Definition 6 these sets of similar subgraphs for
different positions serve as inputs to the same nonlinearities, therefore a
unique PDE with space-dependent coefficients can be obtained.

######  Remark 1.

In theory, it is possible to satisfy Definition 6 for any nonlinear system
formulated through computational graphs. Indeed, assume two computational
graphs have two different root node expressions, denoted as $F(\cdot)$ and
$G(\cdot)$ respectively. Then we can artificially create a new common root
node which will compute $1\cdot F(\cdot)+0\cdot G(\cdot)$ for the first graph
and $0\cdot F(\cdot)+1\cdot G(\cdot)$ for the second. Thus we can satisfy the
first condition of Definition 6, and recursively applying this idea one can
transform any pair of computational graphs into the pair which has the same
structure. However, if the computational graphs of the system are too
different in different points, it can make no sense to represent a system as a
PDE, since it means that the dynamics of different parts of the system has
nothing in common.

###  IV-C Unequally-spaced systems

In the original derivation of (11) we assumed that the nodes are separated by
constant $\Delta x$ in space. This was done for simplicity, and in general
systems with nonuniform spacing can be tackled completely in the same way as
space-dependent systems described in previous subsection. Indeed, assume the
linear system is given by
$\dot{\rho}_{i}=\sum_{j=1}^{N_{i}}a_{ij}\rho_{i+s_{ij}}$, where
$x_{i+1}-x_{i}\not=x_{j+1}-x_{j}$ for $i\not=j$ in general. Then for every
point the continuation can be performed by defining
$c_{ik}=\sum_{j=1}^{N_{i}}a_{ij}(x_{i+s_{ij}}-x_{i})^{k}$. Representing
obtained gains by functions $c_{k}(x)$ as it was described in previous
subsection one finishes with a PDE

$\frac{\partial\rho}{\partial
t}=\sum\limits_{k=1}^{d}c_{k}(x)\frac{1}{k!}\frac{\partial^{k}\rho}{\partial
x^{k}}.$ (30)

###  IV-D Boundary conditions

Now let us look at the Heat PDE:

$\frac{\partial\rho}{\partial t}=\frac{\partial^{2}\rho}{\partial x^{2}}.$
(31)

Imagine that this equation is defined on an interval $x\leavevmode\nobreak\
\in\leavevmode\nobreak\ [0,\leavevmode\nobreak\ +\infty)$, that is there is a
boundary in the point $x=0$.

There are two types of boundary conditions (or BC) which can be supplied to
provide a well-posed boundary value problem. For example for some
$a\in\mathbb{R}$,

1) $\displaystyle\textit{Dirichlet}\text{ BC: }\rho(0)=a,$ (32) 2)
$\displaystyle\textit{Neumann}\text{ BC: }\partial\rho/\partial x\>(0)=a.$

There can also exist a linear combination of these boundary conditions, called
Robin BC.

If the Heat Equation (31) is discretized in stencil points $\\{i-1,i,i+1\\}$,
the result is

$\dot{\rho}_{i}=\frac{1}{\Delta
x^{2}}\left(\rho_{i-1}-2\rho_{i}+\rho_{i+1}\right).$ (33)

Assume now that there exists $i_{0}=1$ such that $x_{i_{0}-1}=0$. Depending on
the type of boundary conditions, the equation for the state $\rho_{1}$ can be
obtained by the discretization of a boundary value problem (31)-(32) in two
ways:

1) $\displaystyle\text{Dirichlet BC:
}\dot{\rho}_{1}=\left(a-2\rho_{1}+\rho_{2}\right)/\Delta x^{2},$ (34) 2)
$\displaystyle\text{Neumann BC:
}\dot{\rho}_{1}=\left(\rho_{2}-\rho_{1}\right)/\Delta x^{2}-a/\Delta x.$

Now imagine the system (31) is obtained by the continuation process from the
system (33). We can notice that states of (33) are governed by the same
dynamics except for the boundary state $\rho_{1}$. The question is how to
recover the boundary conditions (32) for the PDE from the dynamics of
$\rho_{1}$ in (34).

This indeed can be done if one assumes that there exists a ”ghost cell”
$\rho_{0}$ such that it has no dynamics, but is algebraically connected with
adjacent states. With a proper definition of $\rho_{0}$ the equation for
$\dot{\rho}_{1}$ can be represented in the same way as for other states (33)
and thus has the same continuation (31). For example, algebraic equations for
$\rho_{0}$ representing (33)-(34) are

1) $\displaystyle\text{Dirichlet BC: }\rho_{0}=a,$ (35) 2)
$\displaystyle\text{Neumann BC: }\rho_{0}=\rho_{1}-a\Delta x.$

The ghost cell $\rho_{0}=a$ for the case of Dirichlet BC is depicted in Fig.
6. Notice that equations (35) can be directly continualized, obtaining (32).

$\rho_{1}$$\rho_{2}$$a$$x$$0$$\Delta x$$2\Delta x$ Figure 6: Boundary of the
system (33) with Dirichlet boundary condition (34), represented by a ghost
cell $\rho_{0}=a$.

This procedure can be generalized to any ODE system: once the states near
boundaries change their dynamics with respect to the general governing
equation, this change can be represented by ”ghost cells” with algebraic
dependences on the ”real” states. Continualizing these algebraic equations
leads to the boundary conditions for the obtained PDE.

###  IV-E PDE with index derivatives

Usually PDEs have derivatives written with respect to the time and space
variables, thus their physical meaning is in the function continuously varying
in time and space. However, in general no one prevents us from writing a PDE
with respect to some other variables.

Assume a physical system is given by a set of interacting agents, with agents
being indexed by an integer index $i\in\mathbb{Z}$ (a general multiindex space
$\mathbb{Z}^{n}$ can also be used). Let an agent $i$ have a state $\rho_{i}$.
The index variable $i$ is by definition discrete. However we can make an
assumption that in between of two agents with consecutive indexes $i$ and
$i+1$ there is a continuum of virtual agents having state varying from
$\rho_{i}$ to $\rho_{i+1}$. Denoting this continuously varying index by
$M\in\mathbb{R}$ we can say that the state of the system $\rho$ is a
continuous and smooth function $\rho(M)$ with the property $\rho(i)=\rho_{i}$.
This definition of $M$ coincides with the definition of Moskowitz function
used to describe the number of vehicles passed through a fixed point in
traffic modeling [24].

Once the index variable is continuous, we can think about it as a new space
variable. Thus it is possible to use a continuation described in previous
sections, where the distance between two consecutive agents is obviously
$\Delta M=1$. The derivatives of the state with respect to the index can be
obtained by continuation, for example
$\rho_{i+1}-\rho_{i}\;\to\;\partial\rho/\partial M$.

##  V Derivation of the Euler Equations

In the beginning of the XX century Hilbert posed his 6th problem, where he
suggested to develop a rigorous way leading from the atomistic view to the
laws of motion of continua. In particular, the problem can be formulated as a
derivation of the Euler equations for compressible fluids from the Newton’s
dynamics of individual particles.

For the most famous case of particles interacting through collision the
Boltzmann equation was developed, describing evolution of the joint position-
velocity probability distribution of particles. The method of how to transform
individual’s dynamics into Boltzmann equation is based on the Boltzmann-Grad
limit [25], assuming velocities of colliding particles being independent. The
following transformation from the Boltzmann equation to the Euler equations
uses either Hilbert or Chapman-Erskog expansions with space contration limits
[26, 27], Grad moments [28] or the method of invariant manifolds [29].

Another situation arises when the particles interact through long-range
forces. In this case the Vlasov equation can be used instead of the Boltzmann
equation to describe the joint position-velocity probability distribution. The
derivation of the Euler equations from the Vlasov equation was performed in
[30] using space-contracting limit. In particular it was shown that the
resulting system has zero temperature, i.e. the velocities of individual
particles coincide with the velocity field. However, due to the space
contraction the particular form of the potential function was lost and the
obtained pressure was just a square of the density.

Here we present a derivation of Euler equations directly from individual’s
dynamics using the continuation method described in previous sections.
Contrary to other works, we do not use any kind of limits and we use only one
assumption on the isotropy of the space. The assumption requires that for any
particle its nearest neighbours are distributed around uniformly in every
direction, which can be seen as a counterpart to the molecular chaos
hypothesis for the standard derivation of the Boltzmann equation.

###  V-A System of particles

It is assumed that the fluid consists of small particles interacting with each
other, with every particle following simple Newton laws. We will study the
system with $n$ space dimensions, and the particles are assumed to have unit
mass.

We further assume that there is an interaction between each pair of particles
which is given by a force

$F(x_{i}-x_{j})=\frac{x_{i}-x_{j}}{\|x_{i}-x_{j}\|}f(\|x_{i}-x_{j}\|)=(x_{i}-x_{j})\phi(\|x_{i}-x_{j}\|),$
(36)

thus the force acts along the line connecting two particles with the magnitude
$f$ depending only on the distance between particles. For simplicity we also
define a function $\phi(s)=f(s)/s$ representing the scaled magnitude. Since we
assume the infinite number of particles and an infinitely large space, the
magnitude of the force should satisfy

$\int\limits_{\varepsilon}^{+\infty}s^{n-1}f(s)ds<\infty\quad\forall\varepsilon>0$
(37)

such that the cumulative force on any particle is finite. Thus the interaction
should be fast-decaying.

We then need to enumerate all particles. For this we will use multiindex
$i\in\mathbb{Z}^{n}$ as written in section IV-A. Now let us write the dynamics
of a particle with multiindex $i$ using the second Newton’s Law:

$\left\\{\begin{aligned} \dot{x}_{i}&=v_{i},\\\
\dot{v}_{i}&=\sum\limits_{q\not=0}F(x_{i}-x_{i+q}),\end{aligned}\right.$ (38)

where the summation is performed among all multiindices $q$ in
$\mathbb{Z}^{n}\setminus\\{0\\}$, since all the particles interact with each
other. Both the position $x_{i}$ and the velocity $v_{i}$ are vectors in
$\mathbb{R}^{n}$.

###  V-B Derivation in the Euclidean space

Treating the coordinate $x_{i}$ as a state and using the idea written in
section IV-E we define a multiindex function $M(x,t)$ which is the inverse
function of the coordinate: $M(x_{i},t):=i$. Likewise, $x(i,t)=x_{i}(t)$ and
thus $x(M(x,t),t)\equiv x\;\;\forall x\in\mathbb{R}^{n}$.

Now let us write a property of inverse function of multiindex as
$M(x(M,t),t)\equiv M\;\;\forall M\in\mathbb{R}^{n}$, where the space for
multiindices is continuous by the assumption in section IV-E. Taking the time
and the index derivatives, we obtain the following very useful relations on
Jacobians:

$\frac{\partial M}{\partial t}+\frac{\partial M}{\partial x}\frac{\partial
x}{\partial t}=0,$ (39) $\frac{\partial M}{\partial x}\frac{\partial
x}{\partial M}=I.$ (40)

Equation (39) can be seen as a PDE where the function $M$ depends both on $x$
and $t$. Recalling that the multiindex in assumed to be continuous, we can
further utilize the first equation of (38) written in a form
$\dot{x}(M,t)=v(M,t)$, substitute it in (39) and obtain the following equation
on the multiindex evolution:

$\frac{\partial M}{\partial t}=-\frac{\partial M}{\partial
x}v(M(x,t))=-\frac{\partial M}{\partial x}u(x,t),$ (41)

where the velocity function $u(x,t)=v(M(x,t),t)$ is defined as a velocity of a
particle at some given point in space. Finally, taking the derivative with
respect to space, we obtain

$\frac{\partial}{\partial t}\left(\frac{\partial M}{\partial
x}\right)=-\frac{\partial}{\partial x}\left(\frac{\partial M}{\partial
x}u\right).$ (42)

The Jacobian matrix $\frac{\partial M}{\partial x}(x,t)$ represents a
compression tensor, which measures how close are neighbour particles with
respect to different directions in the euclidean space. Evolution of this
Jacobian in the euclidean space is described by the matrix PDE (42), which is
essentially a transport equation with flow velocity given by $u(x,t)$.

Now we approach the second equation in (38). It would be desirable to
transform it in such a way that we could obtain an evolution equation for the
flow velocity $u(x,t)$. First of all, let us rewrite the second equation of
(38) in a way more suitable for continuation, namely

$\dot{v}_{i}=-\sum\limits_{q>0}\left(F(x_{i+q}-x_{i})-F(x_{i}-x_{i-q})\right),$
(43)

where the summation is performed among all multiinidices which are greater
than zero in lexicographical order, i.e. the first nonzero element of $q$
should be positive.

Now we use the continuation of order of accuracy 0 on a multidimensional
system as described in IV-A such that

$x_{i+q}-x_{i}\;\rightarrow\;\frac{\partial x}{\partial
M}\left(x_{i+q/2}\right)q,\quad x_{i}-x_{i-q}\;\rightarrow\;\frac{\partial
x}{\partial M}\left(x_{i-q/2}\right)q,$

which means that (43) becomes

$\dot{v}_{i}=-\sum\limits_{q>0}\left(F\left(\frac{\partial x}{\partial
M}q\right)_{i+q/2}-F\left(\frac{\partial x}{\partial
M}q\right)_{i-q/2}\right).$

Applying the continuation further to the forces, we obtain

$F_{i+q/2}-F_{i-q/2}\quad\rightarrow\quad\frac{\partial F}{\partial
M}(x_{i})q.$

Thus (43) transforms into

$\frac{\partial v}{\partial t}=-\sum\limits_{q>0}\frac{\partial}{\partial
M}\left(\left[\frac{\partial x}{\partial
M}q\right]\phi\left(\left\|\frac{\partial x}{\partial
M}q\right\|\right)\right)q,$ (44)

where we used a definition of the force (36).

Now, we state the following result:

######  Proposition 1.

For any $q\in\mathbb{Z}^{n}$ and for any smooth scalar field $\phi(x)$ the
following identity holds:

$\displaystyle\left[\frac{\partial}{\partial M}\left(\phi\frac{\partial
x}{\partial M}q\right)q\right]^{T}=$ (45)
$\displaystyle\nabla\cdot\left(\phi\frac{\partial x}{\partial
M}qq^{T}\frac{\partial x}{\partial
M}^{T}\right)-\phi\left(\nabla\cdot\left(\frac{\partial x}{\partial
M}\right)qq^{T}\frac{\partial x}{\partial M}^{T}\right),$

where $\nabla$ denotes a row vector of derivatives with respect to $x$.

###### Proof.

First, for convenience denote the left-hand side as a vector $Q$:

$Q:=\frac{\partial}{\partial M}\left(\phi\frac{\partial x}{\partial
M}q\right)q=\frac{\partial}{\partial x}\left(\phi\frac{\partial x}{\partial
M}q\right)\frac{\partial x}{\partial M}q.$ (46)

Also define $h=(\partial x/\partial M)q$. Expanding $\partial(\phi h)/\partial
x$, we get

$Q=h\frac{\partial\phi}{\partial x}h+\frac{\partial h}{\partial
x}h\phi=hh^{T}\frac{\partial\phi}{\partial x}^{T}+\frac{\partial h}{\partial
x}h\phi.$ (47)

Now, for any $h\in\mathbb{R}^{n}$

$\nabla\cdot(hh^{T})=\begin{pmatrix}\sum\limits_{i}h_{1}\frac{\partial
h_{i}}{\partial x_{i}}+\sum\limits_{i}h_{i}\frac{\partial h_{1}}{\partial
x_{i}}&\cdots&\sum\limits_{i}h_{n}\frac{\partial h_{i}}{\partial
x_{i}}+\sum\limits_{i}h_{i}\frac{\partial h_{n}}{\partial
x_{i}}\end{pmatrix},$

which means that

$\left(\nabla\cdot(hh^{T})\right)^{T}=\frac{\partial h}{\partial
x}h+(\nabla\cdot h)h.$ (48)

Therefore the transpose of (47) is

$Q^{T}=\frac{\partial\phi}{\partial
x}hh^{T}+\phi\nabla\cdot(hh^{T})-\phi(\nabla\cdot h)h^{T}.$ (49)

Since for any matrix $J$ and for any scalar field $\alpha$

$\nabla\cdot(\alpha J)=\frac{\partial\alpha}{\partial x}J+\alpha\nabla\cdot
J,$ (50)

we can simplify (49) as

$Q^{T}=\nabla\cdot(\phi hh^{T})-\phi(\nabla\cdot h)h^{T}.$ (51)

The result of the proposition follows by substituting $h$ and noticing that
$\nabla\cdot((\partial x/\partial M)q)=(\nabla\cdot(\partial x/\partial M))q$.
∎

Proposition 1 allows us to rewrite (44) as being dependent only on the
euclidean space divergences and the inverse of the compression tensor
$\partial M/\partial x$. To finalize the derivation of a complete set of
equations, recall the definition of the velocity field $u(x,t)=v(M(x,t),t)$.
Taking the time derivative:

$\frac{\partial u}{\partial t}=\frac{\partial v}{\partial t}+\frac{\partial
v}{\partial M}\frac{\partial M}{\partial t},$

which by (41) is

$\frac{\partial u}{\partial t}=-\frac{\partial v}{\partial M}\frac{\partial
M}{\partial x}u+\frac{\partial v}{\partial t}.$

This equation can be simplified by $\partial u/\partial x=\partial v/\partial
M\cdot\partial M/\partial x$. Finally, substituting (44) and (45) and
combining the result with (42) we obtain a system

$\left\\{\begin{aligned} \frac{\partial}{\partial t}\left(\frac{\partial
M}{\partial x}\right)=&-\frac{\partial}{\partial x}\left(\frac{\partial
M}{\partial x}u\right),\\\ \frac{\partial u}{\partial t}=&-\frac{\partial
u}{\partial x}u-\sum\limits_{q>0}\Bigg{[}\nabla\cdot\left(\phi\frac{\partial
x}{\partial M}qq^{T}\frac{\partial x}{\partial M}^{T}\right)\\\
&-\phi\left(\nabla\cdot\left(\frac{\partial x}{\partial
M}\right)qq^{T}\frac{\partial x}{\partial
M}^{T}\right)\Bigg{]}^{T},\end{aligned}\right.$ (52)

where $\phi=\phi(\left\|(\partial x/\partial M)q\right\|)$.

The system (52) has 12 states in 3-dimensional space, 9 for $\partial
M/\partial x\>(x,t)$ and 3 for $u(x,t)$. It resembles the famous Grad
13-moment system [28], which extends the Euler equations by considering
directional-dependent pressure tensor. The last state of the Grad 13-moment
system is the inner energy, which does not appear in (52). The reason for this
is that we derive a continuous interaction term explicitly from the
interaction forces, which is possible only if the forces are defined by long-
range potentials. As it was shown in [30], expressing a system with long-range
potentials by the Euler equations leads to the solution with zero temperature,
therefore the inner energy becomes functionally dependent on the velocity
field and its evolution equation can be omitted.

###  V-C Dimensionality reduction

It appears that in some special cases it is possible to reduce the system (52)
by considering only one scalar characteristic of a compression in any space
point instead of the whole compression tensor.

Indeed, we define a density as a determinant of the compression tensor,
$\rho(x,t):=\det\left(\partial M/\partial x\right)(x,t)$. Not only the
compression tensor itself, but also its determinant satisfies (42). This
nontrivial fact holds because the compression tensor is the Jacobian, and the
proof is given in Lemma 1 in the Appendix. Therefore from (42)

$\frac{\partial\rho}{\partial t}=-\nabla\cdot\left(\rho u\right).$ (53)

This equation is the first of the complete set of Euler equations.

Unfortunately, the second equation of (52) depends on the whole compression
tensor and thus it cannot be described only by the means of density. This is
reasonable since in general the system can have different pressures in
different directions in response to different compressions. Therefore in order
to simplify the system we need to assume that the compression can be
represented by a single number, i.e. that it is compressed equally in all
directions.

######  Assumption 1 (Isotropy).

Compression tensor $\partial M/\partial x(x,t)$ is isotropic (equal in all
directions), thus it can be represented as a rotation matrix multiplied by a
scalar.

This assumption looks restricting at first glance, but for the infinitely
large system with infinitely many particles the system indeed ”looks the same”
in all directions at every point, thus we can say it is isotropic.

Assumption 1 has long-lasting implications. Define $l(x,t):=\lambda(\partial
x/\partial M(x,t))$, since all the eigenvalues are equal. This variable,
called specific distance, represents an average distance between two
neighbouring particles at point $x$. By definition of the density
$\rho=l^{-n}$. Further, $\left\|\frac{\partial x}{\partial
M}q\right\|=l\left\|q\right\|$. Breaking the summation in (52) in a sum of all
possible lengths $r$ of multiindex vectors, we can rewrite the summation term
as

$\displaystyle\sum\limits_{r^{2}\in\mathbb{N}}\Bigg{[}$
$\displaystyle\nabla\cdot\left(\phi(rl)\frac{\partial x}{\partial
M}\sum\limits_{\begin{subarray}{c}q>0\\\
\left\|q\right\|=r\end{subarray}}\left(qq^{T}\right)\frac{\partial x}{\partial
M}^{T}\right)-$ (54)
$\displaystyle-\phi(rl)\left(\nabla\cdot\left(\frac{\partial x}{\partial
M}\right)\sum\limits_{\begin{subarray}{c}q>0\\\
\left\|q\right\|=r\end{subarray}}\left(qq^{T}\right)\frac{\partial x}{\partial
M}^{T}\right)\Bigg{]}^{T}.$

######  Proposition 2.

Given $r$ such that $r^{2}\in\mathbb{N}$, the summation over all outer
products of multiindices of a length $r$ is proportional to the identity
matrix, i.e. there exists $\beta(r)$ such that

$\sum\limits_{\begin{subarray}{c}q>0\\\
\left\|q\right\|=r\end{subarray}}qq^{T}=\beta(r)I.$ (55)

###### Proof.

First of all, we will show that all nondiagonal elements in (55) are zero.
Indeed, for any positive $q$ its contribution to $kj$-th element of matrix
(55) is given by $q_{k}q_{j}$. But for any $k\not=j$ we can pick $\bar{q}$
such that it equals $q$ except $\bar{q}_{\max(k,j)}\leavevmode\nobreak\
=\leavevmode\nobreak\ -q_{\max(k,j)}$. In this case $\bar{q}$ is also positive
and thus is included into the summation, while the contribution to $kj$-th
element of (55) has opposite sign. Therefore all nondiagonal elements of (55)
are zero.

Further, all diagonal elements of (55) are equal. This can be proven by
analogous argument. Indeed, we can take a positive $q$ and look at the
elements $q_{k}^{2}$ and $q_{j}^{2}$. Then $\bar{q}$ which is equal to $q$
except for $\bar{q}_{k}=\operatorname{sgn}(q_{k})|q_{j}|$ and
$\bar{q}_{j}=\operatorname{sgn}(q_{j})|q_{k}|$ is also positive, but swaps the
contributions between $k$-th and $j$-th diagonal elements. Thus all the
contributions to the diagonal elements are equal. Finally,

$\operatorname{Tr}\sum\limits_{\begin{subarray}{c}q>0\\\
\left\|q\right\|=r\end{subarray}}qq^{T}=\sum\limits_{\begin{subarray}{c}q>0\\\
\left\|q\right\|=r\end{subarray}}q^{T}q=r^{2}\cdot\\#_{r}q=n\beta(r),$ (56)

where $\\#_{r}q$ denotes the number of positive multiindices $q$ with length
$r$ and we define $\beta(r)=r^{2}/n\cdot\\#_{r}q$. It is worth noticing that
by [31] the average approximate behaviour of the number of positive
multiindices $q$ with length $r$ is $\\#_{r}q\propto r^{n-1}$ as
$r\to+\infty$, thus $\beta(r)\propto r^{n}$. ∎

By Assumption 1

$\frac{\partial x}{\partial M}\frac{\partial x}{\partial M}^{T}=l^{2}I.$ (57)

Using Proposition 2 and (57), (54) becomes

$\sum\limits_{r^{2}\in\mathbb{N}}\beta(r)\Bigg{[}\nabla\cdot\left(\phi(rl)l^{2}I\right)-\phi(rl)\left(\nabla\cdot\left(\frac{\partial
x}{\partial M}\right)\frac{\partial x}{\partial M}^{T}\right)\Bigg{]}^{T}.$

The value inside of the square brackets can be simplified further. Indeed, by
(50) it is possible to inject density inside, which gives

$\displaystyle\frac{1}{\rho}\Bigg{[}$
$\displaystyle\nabla\cdot\left(\rho\phi(rl)l^{2}I\right)-\frac{\partial\rho}{\partial
x}\phi(rl)l^{2}I-\phi(rl)\rho\left(\nabla\cdot\left(\frac{\partial x}{\partial
M}\right)\frac{\partial x}{\partial M}^{T}\right)\Bigg{]}^{T}$
$\displaystyle=\frac{1}{\rho}\Bigg{[}\nabla\cdot\left(\rho\phi(rl)l^{2}I\right)-\phi(rl)\left(\nabla\cdot\left(\rho\frac{\partial
x}{\partial M}\right)\frac{\partial x}{\partial M}^{T}\right)\Bigg{]}^{T}.$

Finally, the second term in the square brackets appears to be zero, since
Lemma 2 in the Appendix proves that $\nabla\leavevmode\nobreak\
\cdot\leavevmode\nobreak\ \left(\rho\frac{\partial x}{\partial
M}\right)\leavevmode\nobreak\ \frac{\partial x}{\partial
M}^{T}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0$. Using this Lemma and the
fact that $\nabla\cdot(\rho\phi(rl)l^{2}I)=\nabla(\rho\phi(rl)l^{2})$, we can
define the pressure:

$P=\sum\limits_{r^{2}\in\mathbb{N}}\beta(r)\rho\phi(lr)l^{2}=\sum\limits_{r^{2}\in\mathbb{N}}\frac{\beta(r)}{r}l^{1-n}f(lr).$
(58)

Note that the pressure is well-defined since the sum is convergent by the
property (37). With this definition, the system (52) together with (53) turns
into the famous Euler equations:

$\left\\{\begin{aligned} \frac{\partial\rho}{\partial
t}&=-\nabla\cdot\left(\rho u\right),\\\ \frac{\partial u}{\partial
t}&=-\frac{\partial u}{\partial x}u-\frac{\nabla
P}{\rho}^{T}.\end{aligned}\right.$ (59)

Therefore the following theorem was proven:

######  Theorem 3.

There exists a valid continuation process which leads from the Newtonian
system (38) to the Euler equations (59) under the assumption that the system
is locally isotropic in every point in space.

######  Remark 2 (Non-complete interaction topologies).

In the original ODE system (38) we assumed that an interaction exists between
every pair of particles, i.e. that the topology of interactions is all-to-all.
In general in order to obtain (38) it would be sufficient to use any topology
for which the isotropy required in Assumption 1 is possible. The difference in
topologies would modify the definitions of density $P(x,t)$ in (58).

For example, for the grid topology with equations given by

$\left\\{\begin{aligned} \dot{x}_{i}&=v_{i},\\\
\dot{v}_{i}&=\sum\limits_{k=1}^{n}\left(F(x_{i}-x_{i-e_{k}})-F(x_{i}-x_{i+e_{k}})\right),\end{aligned}\right.$
(60)

where $e_{k}$ denotes the $k$-th basis vector of $\mathbb{R}^{n}$, the
continuation renders the same Euler equations (59) with the pressure given by
$P=f(l)/l^{n-1}$.

##  VI Control of Robotic Swarm

In this section we will demonstrate how the continuation method described
above can help in the analysis and design of control laws for large-scale
systems. We will do it by using an example of a robotic swarm, i.e. a
formation of robots whose goal is to follow some desired trajectory while
passing through obstacles and preserving relative agents’ positions.

Control of robotic formations is an extensively studied topic, see recent
reviews [32, 33]. However most of the methods rely on the graph-theoretic
properties of interaction topology and on simple linear controllers to provide
stability. A PDE approach was taken in [34] where the Euler PDE with diffusion
terms was used to model the flocks of birds. The authors proposed a PDE to
describe the behaviour of agents and analyzed it to study a symmetry breaking
which leads to a coherent movement of birds. Similar PDE model was used to
control 3D agent formation with 2D disc communication topology via
backstepping in [35]. Another interesting concept was used in [36, 37], namely
the Partial difference Equations (PdEs), which are an analogue of ordinary
Partial Differential Equations defined on graphs. With the help of Lyapunov
analysis within PdE formalism it was proven that the linear diffusion
controller stabilizes the formation.

Contrary to previous works we will base our analysis on the continuation
procedure, rigorously introducing a PDE to describe a formation of drones. We
will study this PDE and recover a nonlinear local control law which, being
applied to the agents, forces the whole formation to follow the desired
density profile.

###  VI-A Continuation and PDE Control

Let us start from a system of drones having double integrator dynamics:

$\ddot{x}_{i}=\tau_{i}.$ (61)

Here $x_{i}\in\mathbb{R}^{n}$ is a position of the $i$-th drone in
$n$-dimensional space and $\tau_{i}\in\mathbb{R}^{n}$ is a control we want to
design. The drones are enumerated with multiindices $i\in\mathbb{Z}^{n}$.
Define $v_{i}=\dot{x}_{i}$. Similarly to the previous section we introduce
multiindex function $M(x,t)$ such that $M(x_{i},t)\equiv i$ and then perform a
continuation. The resulting system is

$\left\\{\begin{aligned} \frac{\partial\rho}{\partial t}&=-\nabla\cdot(\rho
u),\\\ \frac{\partial u}{\partial t}&=-\frac{\partial u}{\partial
x}u+\tau(x,t),\end{aligned}\right.$ (62)

where $\tau(x,t)=\tau(M(x,t),t)$ is a continuation of the control $\tau_{i}$.

Now let us formulate the desired system which will be used as a reference the
real formation should converge to. Given a velocity profile $u_{d}(x)$, we
define the desired density $\rho_{d}(x,t)$ to follow this velocity profile.
Essentially this means ”desired agents” have single-integrator dynamics. Note
that in general $u_{d}$ can be dependent on time but we don’t consider it for
simplicity of writing.

Thus we assume the desired system is governed by

$\frac{\partial\rho_{d}}{\partial t}=-\nabla\cdot(\rho_{d}u_{d}).$ (63)

Our goal is to derive $\tau(x,t)$ such that $\rho\to\rho_{d}$. First, direct
calculations from (62) and (63) lead to the following systems in terms of
flows $(\rho u)$ and $(\rho_{d}u_{d})$:

$\displaystyle\frac{\partial(\rho u)}{\partial t}$
$\displaystyle=-\nabla\cdot(\rho u)u-\rho\frac{\partial u}{\partial
x}u+\rho\tau(x,t),$ (64) $\displaystyle\frac{\partial(\rho_{d}u_{d})}{\partial
t}$ $\displaystyle=-\nabla\cdot(\rho_{d}u_{d})u_{d}.$

Define the deviation from the desired density $\tilde{\rho}=\rho-\rho_{d}$.
Then the second-order equation for the deviation is

$\frac{\partial^{2}\tilde{\rho}}{\partial
t^{2}}=\nabla\cdot\left[\nabla\cdot(\rho
u)u-\nabla\cdot(\rho_{d}u_{d})u_{d}+\rho\frac{\partial u}{\partial
x}u-\rho\tau(x,t)\right].$

In order to cancel the nonlinear terms, define now the control $\tau$ as

$\displaystyle\tau=\frac{\partial u}{\partial x}u+\frac{1}{\rho}\Big{[}\nabla$
$\displaystyle\cdot(\rho u)u-\nabla\cdot(\rho_{d}u_{d})u_{d}+{}$ (65)
$\displaystyle{}+\alpha(\rho_{d}u_{d}-\rho
u)+\beta\nabla(\rho_{d}-\rho)^{T}\Big{]},$

where $\alpha$ and $\beta$ are some positive gains. Then the equation for the
density deviation transforms into

$\frac{\partial^{2}\tilde{\rho}}{\partial
t^{2}}=-\alpha\frac{\partial\tilde{\rho}}{\partial
t}+\beta\nabla^{2}\tilde{\rho}.$ (66)

This equation is a wave equation with damping and thus it is asymptotically
stable if $\tilde{\rho}=0$ on the boundary of the domain [38]. Choosing a
desired system such that $\rho_{d}=0$ on the boundary and using a continuation
of $\rho$ such that $\rho=0$ on the boundary ensures satisfaction of the
boundary condition.

###  VI-B Discretization of the control

Formula (65) for PDE (62) is local by its nature, but it should be discretized
to be implemented on every agent of the original ODE (61). One particular
discretization is described next.

First of all, for the agent $i$ define a matrix $G_{i}$ as a discretization of
the compression tensor:

$[G_{i}]_{j}=(x_{i+e_{j}}-x_{i-e_{j}})/2\approx\frac{\partial x}{\partial
M_{j}}(x_{i},t),$ (67)

where $e_{j}$ is the $j$-th unit basis vector and $[G_{i}]_{j}$ represent the
$j$-th column of $G_{i}$. The matrix $G_{i}$ depends on the positions of $2n$
neighbouring agents of the $i$-th agent. In the same way as $G_{i}$ we define
a matrix $W_{i}$ representing a velocity Jacobian:

$[W_{i}]_{j}=(v_{i+e_{j}}-v_{i-e_{j}})/2\approx\frac{\partial u}{\partial
M_{j}}(x_{i},t).$ (68)

Now we can write formulas for all terms inside of (65) depending on the real
system:

1). $\displaystyle\quad\frac{\partial u}{\partial x}u=\frac{\partial
u}{\partial M}\frac{\partial M}{\partial x}u\approx W_{i}G_{i}^{-1}v_{i},$
(69) 2). $\displaystyle\quad\nabla\cdot u=\sum\limits_{j=1}^{n}\frac{\partial
u_{j}}{\partial M}\frac{\partial M}{\partial
x_{j}}\approx\sum\limits_{j=1}^{n}[W_{i}^{T}]_{j}\cdot[G_{i}^{-1}]_{j},$ 3).
$\displaystyle\quad\rho\approx 1/\det G_{i},$ 4).
$\displaystyle\quad\nabla\rho\approx-\rho^{2}\nabla(\det
G_{i})\approx-\rho^{2}\frac{\partial(\det G_{i})}{\partial M}G_{i}^{-1},$

where the gradient of the determinant $\det G_{i}$ should be computed
according to the determinant formula, using second derivatives of the
positions discretized similarly to (67):

$\displaystyle\frac{\partial^{2}x}{\partial M_{j}\partial M_{k}}(x_{i},t)$
$\displaystyle\approx(x_{i+e_{j}+e_{k}}+x_{i-e_{j}-e_{k}}-x_{i+e_{j}-e_{k}}-x_{i-e_{j}+e_{k}})/4,$
$\displaystyle\frac{\partial^{2}x}{\partial M_{j}^{2}}(x_{i},t)$
$\displaystyle\approx x_{i+e_{j}}-2x_{i}+x_{i-e_{j}}.$

Since the gradient of the determinant depends on the second derivatives, in
total each agent requires information about the velocities of its $2n$
neighbouring agents and the positions of its $2n^{2}$ neighbouring agents,
including diagonal ones.

Finally, substituting (69) into (65), the formula for the control action
$\tau_{i}$ appears as

$\displaystyle\tau_{i}$
$\displaystyle=\left[W_{i}G_{i}^{-1}+\sum\limits_{j=1}^{n}[W_{i}^{T}]_{j}\cdot[G_{i}^{-1}]_{j}-\alpha\right]v_{i}+{}$
(70) $\displaystyle{}+\left[\beta I-v_{i}v_{i}^{T}\right]\frac{1}{\det
G_{i}}G_{i}^{-T}\frac{\partial(\det G_{i})}{\partial M}^{T}+{}$
$\displaystyle{}+\det G_{i}\Big{[}\alpha\rho_{d}u_{d}+(\beta
I-u_{d}u_{d}^{T})\nabla\rho_{d}^{T}-\rho_{d}(\nabla\cdot u_{d})u_{d}\Big{]}.$

###  VI-C Boundary conditions

For the system (66) to converge to zero proper boundary conditions should be
used. Namely, the continuation should be chosen such that $\rho=0$ outside of
the formation. As it was shown in IV-D, boundary conditions for PDE correspond
to ”ghost agents” in the ODE case. In particular, information about neighbour
agents is used in (67) and (68). Therefore specifying boundary conditions
means specifying positions $x_{i\pm e_{j}}$ and velocities $v_{i\pm e_{j}}$
for the nonexisting agents.

Assume agent $i-e_{j}$ is a ghost agent. One natural choice, which we will use
for velocities, is to take $v_{i-e_{j}}=2v_{i}-v_{i+e_{j}}$. Being substituted
in (68) this leads to an approximation of the velocity gradient based solely
on the $i$ and $i+e_{j}$ agents.

$\rho$$x$$x_{i-e_{j}}$$x_{i}$$x_{i+e_{j}}$$x_{i+2e_{j}}$$s=2l$$l$$l$$\rho_{i-e_{j}}=0$$\rho_{i}=\frac{2}{3}l^{-1}$$\rho_{i+e_{j}}=\rho_{i+2e_{j}}=l^{-1}$
Figure 7: Left boundary of the system (61) with control (70). Agent $i$ is on
the boundary, the position of the ”ghost agent” $i-e_{j}$ is chosen such that
$\rho$ linearly goes to zero at $x_{i-e_{j}}$.

This idea can’t be used for $x_{i-e_{j}}$ since in this case the compression
tensor (67) will not ”feel” that the drone $i$ is on the border. Instead we
will use such an approximation that the density near the border will linearly
diminish to zero, see Fig. 7. Namely, let us look at 1D case and fix $i$-th
agent to be on the left border. Assume further that the distance between each
pair of existing agents is constant and equal to $l$. Then
$\rho_{i+e_{j}}=l^{-1}$. Also we define $\rho_{i-e_{j}}=0$ such that the ghost
agent has its density zero. Define an unknown distance $s:=x_{i}-x_{i-e_{j}}$.
Then asking for a linear dependency of a density on position, we have
necessarily

$\rho_{i}=\frac{l\rho_{i-e_{j}}+s\rho_{i+e_{j}}}{l+s}=\frac{s}{l(l+s)}.$

But by (67) $\rho_{i}=2/(l+s)$, which immediately gives the answer $s=2l$, or

$x_{i-e_{j}}=x_{i}+2(x_{i}-x_{i+e_{j}})=3x_{i}-2x_{i+e_{j}}.$ (71)

This finalizes the formulation of the boundary conditions and thus the correct
implementation of (70).

###  VI-D Numerical Simulation

To demonstrate the control policy (70) we performed a numerical simulation of
a cubic formation of 512 drones in 3D space. The goal was to reach a cubic
formation, fly through a window and restore the cubic formation after the
maneuver.

Assume the center of the window is placed at the point $(x_{0},0,0)$, and the
formation should fly through it starting from the origin. The desired velocity
field $u_{d}(x,y,z)$ able to fulfill the task was constructed as

$u_{d_{x}}=1,\quad
u_{d_{y|z}}=0.05\operatorname{atan}(x-x_{0})e^{-\frac{(x-x_{0})^{2}}{100}}y|z,$

where $y|z$ denotes $y$ or $z$, see the left panel of Fig. 8 for the
streamlines projected on the $x$-$y$ plane. For simplicity the desired system
(63) was simulated by first-order integrators following the desired velocity
profile, and the density $\rho_{d}(x,t)$ was interpolated between agents.

Both the desired system (63) and the real system (61) were simulated for the
cubic formation of $8\times 8\times 8$ drones using Euler method. The initial
positions for the real system were multiplied by 2 in comparison to the
desired system and a uniform noise $U(-2,2)$ was added. The control gains were
chosen as $\alpha=3$ and $\beta=100$. The convergence of the real density to
the desired one is shown on the right panel of Fig. 8 and snapshots of the
simulation are presented in Fig. 9. It is clear that the real formation, being
heavily disturbed in the beginning, converges to the desired shape in less
than 5 seconds and then follows the desired pattern, successfully passing
through the window.

Figure 8: Left: streamlines of the desired velocity field $u_{d}(x,y)$. Right:
convergence of the $L_{2}$ norm of the density deviation. Figure 9: Simulation
of drones flying through window. Rows correspond to times
$t=\\{0s,1s,5s,12s,20s,30s,45s\\}$. Left column, reference: desired system
(63), governed by single integrators. Right column, actual: heavility
perturbed real system (61) with control (70) which converges to the desired
one.

##  VII Conclusion

We presented a general process of transformation of ODE systems into their PDE
counterparts, defining the continuation to be valid if the original ODE system
could be obtained from the PDE version by a correct discretization. We further
showed that the spectrum of PDE converges to the spectrum of ODE. The
continuation method was then elaborated for many classes of systems including
nonlinear, multidimensional and space- and time-varying. Based on this method,
new continuous models can be derived and further utilized for analysis and
control purposes.

As an example we used the continuation to show how the Euler equations for
compressible fluid can be derived from the newtonian particle interactions,
providing more intuition into Hilbert’s 6th problem. The same continuation was
then used to describe a robot formation flying through window. We developed a
control algorithm to stabilize a desired trajectory based on a continuous
representation of the formation. This algorithm is distributed as every robot
requires information only about neighbouring robots.

It would be desirable to study further the continuation method, namely to
provide quantitative measures of how close is a PDE solution compared to the
original ODE one.

######  Lemma 1.

Let $J(x,t)\in\mathbb{R}^{n\times n}$ be the Jacobian matrix of function
$M(x,t)$. Let $J(x,t)$ satisfies the dynamic equation

$\frac{\partial J}{\partial t}=-\frac{\partial(Ju)}{\partial x},$ (72)

where $u=u(x,t)$ is some vector field. Then the determinant $\det J$ satisfies
the same equation:

$\frac{\partial\det J}{\partial t}=-\frac{\partial}{\partial x}\cdot\left(\det
J\cdot u\right).$ (73)

###### Proof.

First of all let us rewrite (72) for one element $J_{ik}$ of the matrix $J$:

$\displaystyle\frac{\partial J_{ik}}{\partial t}$
$\displaystyle=-\frac{\partial(J_{i}\>u)}{\partial
x_{k}}=-\sum\limits_{j=1}^{n}\frac{\partial^{2}M_{i}}{\partial x_{k}\partial
x_{j}}u_{j}-\sum\limits_{j=1}^{n}J_{ij}\frac{\partial u_{j}}{\partial x_{k}}=$
(74) $\displaystyle=-\sum\limits_{j=1}^{n}\frac{\partial J_{ik}}{\partial
x_{j}}u_{j}-\sum\limits_{j=1}^{n}J_{ij}\frac{\partial u_{j}}{\partial x_{k}},$

where we used the fact that $J=\partial M/\partial x$.

Now let us recall the definition of the determinant: $\det
J=\sum_{\sigma}\operatorname{sgn}(\sigma)\prod_{i=1}^{n}J_{\sigma_{i},i}$,
where $\sigma$ is a permutation of the set $\\{1,2...n\\}$ and $\sum_{\sigma}$
is taken over all possible permutations, with $\operatorname{sgn}(\sigma)$
being the sign of the permutation. Let us take the time derivative and then
substitute (74):

$\displaystyle\frac{\partial\det J}{\partial
t}=\sum\limits_{\sigma}\operatorname{sgn}(\sigma)\sum\limits_{k=1}^{n}\frac{\partial
J_{\sigma_{k},k}}{\partial t}\prod\limits_{i=1,i\not=k}^{n}J_{\sigma_{i},i}=$
(75)
$\displaystyle=-\sum\limits_{j=1}^{n}\sum\limits_{\sigma}\operatorname{sgn}(\sigma)\sum\limits_{k=1}^{n}\left[\frac{\partial
J_{\sigma_{k},k}}{\partial x_{j}}u_{j}+J_{\sigma_{k},j}\frac{\partial
u_{j}}{\partial x_{k}}\right]\prod\limits_{i=1,i\not=k}^{n}J_{\sigma_{i},i}$

We will investigate two parts of (75), corresponding to the first and the
second terms inside the square brackets. For the first part we have

$\displaystyle-\sum\limits_{j=1}^{n}\sum\limits_{\sigma}\operatorname{sgn}(\sigma)\sum\limits_{k=1}^{n}\frac{\partial
J_{\sigma_{k},k}}{\partial
x_{j}}u_{j}\prod\limits_{i=1,i\not=k}^{n}J_{\sigma_{i},i}$ $\displaystyle=$
(76) $\displaystyle=-\sum\limits_{j=1}^{n}\frac{\partial\det J}{\partial
x_{j}}u_{j}$ $\displaystyle=-\frac{\partial\det J}{\partial x}u.$

The second part is a little bit more tricky:

$\displaystyle-\sum\limits_{j=1}^{n}\sum\limits_{\sigma}\operatorname{sgn}(\sigma)\sum\limits_{k=1}^{n}J_{\sigma_{k},j}\frac{\partial
u_{j}}{\partial x_{k}}\prod\limits_{i=1,i\not=k}^{n}J_{\sigma_{i},i}=-\det
J\sum\limits_{j=1}^{n}\frac{\partial u_{j}}{\partial x_{j}}-{}$
$\displaystyle{}-\sum\limits_{j=1}^{n}\sum\limits_{\sigma}\operatorname{sgn}(\sigma)\sum\limits_{k=1,k\not=j}^{n}J_{\sigma_{k},j}\frac{\partial
u_{j}}{\partial x_{k}}\prod\limits_{i=1,i\not=k}^{n}J_{\sigma_{i},i}.$

Here we split the summation over $k$ into the term with $k=j$ and all other
terms. The first one immediately gives the determinant multiplied by the
divergence of the vector field. It appears that the sum over all other terms
is zero. Indeed, imagine a permutation $\bar{\sigma}$ such that it is equal to
$\sigma$ except $\sigma_{j}$ and $\sigma_{k}$ are swapped. Then the sign of
$\bar{\sigma}$ is opposite to the sign of $\sigma$. Further, since the product
$J_{\sigma_{k},j}J_{\sigma_{j},j}$ is the only way in which $\sigma_{k}$ and
$\sigma_{j}$ enter the formula, the absolute value does not change with the
change of permutation. Therefore for each $j,k$ and for each permutation there
exists a permutation which cancels them out.

Finally, substitution of the nonzero term of the last equation and (76) into
(75) leads to (73). ∎

######  Lemma 2.

Let $\partial x/\partial M$ be isotropic, i.e. represented by a scalar
multiplied by a rotation matrix, and let $\rho\leavevmode\nobreak\
=\leavevmode\nobreak\ \det(\partial M/\partial x)$. Then

$\nabla\cdot\left(\rho\frac{\partial x}{\partial M}\right)\frac{\partial
x}{\partial M}^{T}=0.$ (77)

###### Proof.

Define $\lambda=\lambda(\partial M/\partial x)$, thus $\rho=\lambda^{n}$. By
isotropy

$\frac{\partial x}{\partial M}=\lambda^{-2}\frac{\partial M}{\partial x}^{T},$
(78)

therefore the left-hand side of (77) is

$\displaystyle\nabla\cdot\left(\lambda^{n-2}\frac{\partial M}{\partial
x}^{T}\right)\frac{\partial M}{\partial x}\lambda^{-2}=$ (79)
$\displaystyle=\lambda^{n-4}\nabla\cdot\left(\frac{\partial M}{\partial
x}^{T}\right)\frac{\partial M}{\partial
x}+(n-2)\lambda^{n-3}\frac{\partial\lambda}{\partial x},$

where we used (50) and (78). Now let us investigate the first term more
closely. Taking the divergence and looking at $j$-th element, we see that

$\left[\nabla\cdot\left(\frac{\partial M}{\partial x}^{T}\right)\frac{\partial
M}{\partial x}\right]_{j}=\sum\limits_{k=1}^{n}\frac{\partial^{2}M}{\partial
x_{k}^{2}}^{T}\frac{\partial M}{\partial x_{j}}.$ (80)

Now, By isotropy

$\frac{\partial M}{\partial x_{j}}^{T}\frac{\partial M}{\partial
x_{k}}=0\;\;\forall j\not=k,\quad\frac{\partial M}{\partial
x_{k}}^{T}\frac{\partial M}{\partial x_{k}}=\lambda^{2}.$ (81)

Taking the derivative of the multiplication of basis vectors:

$\frac{\partial}{\partial x_{j}}\left(\frac{\partial M}{\partial
x_{k}}^{T}\frac{\partial M}{\partial
x_{k}}\right)=2\frac{\partial^{2}M}{\partial x_{j}\partial
x_{k}}^{T}\frac{\partial M}{\partial x_{k}},$ (82)

but at the same time the value under the derivative is $\lambda^{2}$ by (81),
therefore

$\frac{\partial}{\partial x_{j}}\left(\frac{\partial M}{\partial
x_{k}}^{T}\frac{\partial M}{\partial
x_{k}}\right)=\frac{\partial\lambda^{2}}{\partial
x_{j}}=2\lambda\frac{\partial\lambda}{\partial x_{j}}.$ (83)

Then, taking the derivative of multiplication of different basis vectors with
$j\not=k$, by (81) we obtain zero:

$\frac{\partial}{\partial x_{k}}\left(\frac{\partial M}{\partial
x_{j}}^{T}\frac{\partial M}{\partial
x_{k}}\right)=\frac{\partial^{2}M}{\partial x_{j}\partial
x_{k}}^{T}\frac{\partial M}{\partial x_{k}}+\frac{\partial M}{\partial
x_{j}}^{T}\frac{\partial^{2}M}{\partial x_{k}^{2}}=0,$

which by equality of (82) and (83) means that for $j\not=k$

$\frac{\partial M}{\partial x_{j}}^{T}\frac{\partial^{2}M}{\partial
x_{k}^{2}}=-\frac{\partial^{2}M}{\partial x_{j}\partial
x_{k}}^{T}\frac{\partial M}{\partial
x_{k}}=-\lambda\frac{\partial\lambda}{\partial x_{j}}.$ (84)

In the case of $j=k$ by equality of (82) and (83) we have

$\frac{\partial^{2}M}{\partial x_{j}^{2}}^{T}\frac{\partial M}{\partial
x_{j}}=\lambda\frac{\partial\lambda}{\partial x_{j}}.$ (85)

Combination of (84) and (85) means that (80) is

$\left[\nabla\cdot\left(\frac{\partial M}{\partial x}^{T}\right)\frac{\partial
M}{\partial x}\right]_{j}=(2-n)\lambda\frac{\partial\lambda}{\partial x_{j}}.$
(86)

Finally, substituting (86) in (79) gives zero. ∎

## Acknowledgment

The Scale-FreeBack project has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement N 694209).

## References

  * [1] L. Euler. Principia motus fluidorum. Novi commentarii academiae scientiarum Petropolitanae, 271-311. 1761.
  * [2] J. C. Maxwell. A treatise on electricity and magnetism (Vol. 1). Clarendon press. 1873.
  * [3] J. G. Verwer and J. M. Sanz-Serna. Convergence of method of lines approximations to partial differential equations. Computing, 33(3-4), 297-313. 1984.
  * [4] W. F. Ames. Numerical methods for partial differential equations. Academic press. 2014.
  * [5] M. Aoki. Control of large-scale dynamic systems by aggregation. IEEE Transactions on Automatic Control, vol. 13, no. 3, pp. 246–253, 1968.
  * [6] M. U. B. Niazi, X. Cheng, C. Canudas-de-Wit and J. Scherpen. Structure-based Clustering Algorithm for Model Reduction of Large-scale Network Systems. CDC 2019 - 58th IEEE Conference on Decision and Control, Nice, France, Dec 2019.
  * [7] M. U. B. Niazi, D. Deplano, C. Canudas-de-Wit, and A. Y. Kibangou. Scale-free estimation of the average state in large-scale systems. IEEE Control Systems Letters, vol. 4, no. 1, pp. 211–216, Jan 2020.
  * [8] D. Nikitin, C. Canudas-de-Wit and P. Frasca. Control of Average and Deviation in Large-Scale Linear Networks. Submitted to IEEE Transactions on Automatic Control.
  * [9] J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort and R. Spigler. The Kuramoto model: A simple paradigm for synchronization phenomena. Reviews of modern physics, 77(1), 137. 2005.
  * [10] D. Q. Nykamp and D. Tranchina. A population density approach that facilitates large-scale modeling of neural networks: Analysis and an application to orientation tuning. Journal of computational neuroscience, 8(1), 19-50. 2000.
  * [11] H. Grabert. Projection operator techniques in nonequilibrium statistical mechanics. Springer Tracts in Modern Physics, vol. 95, Springer-Verlag, Berlin-New York, 1982.
  * [12] Y. Yang, D. V. Dimarogonas and X. Hu. Shaping up crowd of agents through controlling their statistical moments. 2015 European Control Conference (ECC), pp. 1017-1022, Linz, 2015.
  * [13] S. Zhang, A. Ringh, X. Hu and J. Karlsson. Modeling collective behaviors: A moment-based approach. IEEE Transactions on Automatic Control, 2020.
  * [14] C. Kuehn. Moment closure—a brief review. In Control of self-organizing nonlinear systems (pp. 253-271). Springer, Cham. 2016.
  * [15] D. Nikitin, C. Canudas-de-Wit, P. Frasca. Shape-based nonlinear model reduction for 1D conservation laws. 21st IFAC World Congress 2020, Berlin, Germany, July 11-17, 2020.
  * [16] D. Glasscock. a Graphon?. Notices of the AMS, 62(1). 2015.
  * [17] L. Lovász. Large networks and graph limits (Vol. 60). American Mathematical Soc.. 2012.
  * [18] S. Gao and P. E. Caines. Graphon Control of Large-scale Networks of Linear Systems. in IEEE Transactions on Automatic Control, 2019.
  * [19] R. Vizuete, P. Frasca and F. Garin. Graphon-based sensitivity analysis of SIS epidemics. IEEE Control Systems Letters, IEEE, 2020, 4 (3), pp.542 - 547.
  * [20] A. G. Baydin, B. A. Pearlmutter, A. A. Radul and J. M. Siskind. Automatic differentiation in machine learning: a survey. The Journal of Machine Learning Research, 18(1), 5595-5637. 2017.
  * [21] B. Bamieh, F. Paganini and M. A. Dahleh. Distributed control of spatially invariant systems. IEEE Transactions on automatic control, 47(7), 1091-1107. 2002.
  * [22] A. E. Frazho and W. Bhosri. Toeplitz and Laurent Operators. In: An Operator Perspective on Signals and Systems. Operator Theory: Advances and Applications (Linear Operators and Linear Systems), vol 204. Birkhäuser Basel. 2010.
  * [23] A. N. Kolmogorov. On the representation of continuous functions of several variables by superposition of continuous functions of one variable and addition. Dokl. Akad. Nauk SSSR 114, 953-956. 1957.
  * [24] G. F. Newell. A simplified theory of kinematic waves in highway traffic, part I: General theory. Transportation Research Part B: Methodological, 27(4), 281-287. 1993.
  * [25] I. Gallagher, L. Saint-Raymond and B. Texier. From Newton to Boltzmann: hard spheres and short-range potentials. European Mathematical Society. 2013.
  * [26] L. Saint-Raymond. Hydrodynamic limits of the Boltzmann equation. Lecture Notes in Mathematics, vol. 1971. Springer-Verlag, 2009.
  * [27] S. Chapman and T. G. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, Cambridge, 1970.
  * [28] H. Grad. On the kinetic theory of rarefied gases. Commun. Pure Appl. Math. 2, 325, 1949.
  * [29] A. N. Gorban and I. Karlin. Hilbert’s 6th Problem: exact and approximate hydrodynamic manifolds for kinetic equations. Bull. Amer. Math. Soc. 51 (2): 186–246, 2014.
  * [30] S. Caprino, R. Esposito, R. Marra and M. Pulvirenti. Hydrodynamic limits of the Vlasov equation. Communications in partial differential equations, 18(5-6), 805-820. 1993.
  * [31] R. Takloo-Bighash. How many lattice points are there on a circle or a sphere?. In: A Pythagorean Introduction to Number Theory. Undergraduate Texts in Mathematics. Springer, Cham. 2018.
  * [32] K. K. Oh, M. C. Park and H. S. Ahn. A survey of multi-agent formation control. Automatica, 53, 424-440. 2015.
  * [33] S. Chung, A. A. Paranjape, P. Dames, S. Shen and V. Kumar. A Survey on Aerial Swarm Robotics. IEEE Transactions on Robotics, vol. 34, no. 4, pp. 837-855, Aug. 2018.
  * [34] J. Toner and Y. Tu, Long-range order in a two-dimensional dynamical XY model: How birds fly together. Phys. Rev. Lett., vol. 75, no. 3, pp. 4326–4329, 1995.
  * [35] J. Qi, R. Vazquez and M. Krstic. Multi-agent deployment in 3-D via PDE control. IEEE Transactions on Automatic Control, 60(4), 891-906. 2014.
  * [36] G. Ferrari-Trecate, A. Buffa, and M. Gati. Analysis of coordination in multi-agent systems through partial difference equations. IEEE Trans. Autom. Control, vol. 51, no. 6, pp. 1058–1063, Jun. 2006.
  * [37] M. Ji, G. Ferrari-Trecate, M. Egerstedt, and A. Buffa. Containment control in mobile networks. IEEE Trans. Autom. Control, vol. 53, no. 8, pp. 1972–1975, Sep. 2008.
  * [38] G. B. Folland. Introduction to partial differential equations (Vol. 102). Princeton university press. 1995.

|  Denis Nikitin received both the B.Sc. and the M.Sc. degree in Mathematics
and Mechanics Faculty of Saint Petersburg State University, Saint Petersburg,
Russia, specializing on the control theory and cybernetics. He won several
international robotics competitions while being student and was a teacher of
robotics in the Math&Phys Lyceum 239 in Saint Petersburg. He is currently a
doctoral researcher at CNRS, GIPSA-Lab, Grenoble, France. His current research
mainly focuses on control of large-scale systems.
---|---
|  Carlos Canudas-de-Wit (F’16) was born in Villahermosa, Mexico, in 1958. He
received the B.S. degree in electronics and communications from the Monterrey
Institute of Technology and Higher Education, Monterrey, Mexico, in 1980, and
the M.S. and Ph.D. degrees in automatic control from the Department of
Automatic Control, Grenoble Institute of Technology, Grenoble, France, in 1984
and 1987, respectively. He is currently a Directeur de recherche (Senior
Researcher) with CNRS, Grenoble, where he is the Leader of the NeCS Team, a
joint team of GIPSA-Lab (CNRS) and INRIA, on networked controlled systems. Dr.
Canudas-de-Wit is an IFAC Fellow. He was an Associate Editor of the IEEE
Transactions on Automatic Control, the Automatica, the IEEE Transactions on
Control Systems Technology. He is an Associate Editor of the Asian Journal of
Control, and the IEEE Transactions on Control of Network Systems. He served as
the President of the European Control Association from 2013 to 2015, and a
member of the IEEE Board of Governors of the Control System Society from 2011
to 2014. He holds the ERC Advanced Grant Scale-FreeBack from 2016 to 2021.
---|---
|  Paolo Frasca (M’13–SM’18) received the Ph.D. degree from Politecnico di
Torino, Turin, Italy, in 2009. After Postdoctoral appointments with the CNR-
IAC, Rome, and in Torino, he has been an Assistant Professor with the
University of Twente, Enschede, The Netherlands, from 2013 to 2016. Since
October 2016, he has been a CNRS Researcher with GIPSA-lab, Grenoble, France.
His research interests include the theory of network systems and cyber-
physical systems, with applications to robotics, sensors, infrastructural, and
social networks. On these topics, he has (co)authored more than 60 journal and
conference papers and the book Introduction to Averaging Dynamics Over
Networks (Springer). He has been an Associate Editor for the Editorial Boards
of several conferences and journals, including the IEEE Control Systems
Letters.
---|---