# Robust feedback stabilization of interacting multi-agent systems under
uncertainty

Giacomo Albi Department of Computer Science, University of Verona, Str. Le
Grazie 15, Verona, I-37134, Italy  Michael Herty IGPM, RWTH Aachen University,
Templergraben, 55, D-52062 Aachen, Germany  and Chiara Segala IGPM, RWTH
Aachen University, Templergraben, 55, D-52062 Aachen, Germany
MH and CS thank the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) for the financial support through 320021702/GRK2326,
333849990/IRTG-2379, HE5386/18-1,19-2,22-1,23-1 and under Germany’s Excellence
Strategy EXC-2023 Internet of Production 390621612. GA thanks the Italian
Ministry of Instruction, University and Research (MIUR) to support this
research with funds coming from PRIN Project 2017 (No. 2017KKJP4X entitled
“Innovative numerical methods for evolutionary partial differential equations
and applications”).

###### Abstract

We consider control strategies for large-scale interacting agent systems under
uncertainty. The particular focus is on the design of robust controls that
allow to bound the variance of the controlled system over time. To this end we
consider $\mathcal{H}_{\infty}$ control strategies on the agent and mean field
description of the system. We show a bound on the $\mathcal{H}_{\infty}$ norm
for a stabilizing controller independent on the number of agents. Furthermore,
we compare the new control with existing approaches to treat uncertainty by
generalized polynomial chaos expansion. Numerical results are presented for
one-dimensional and two-dimensional agent systems.

Keywords. Agent-based dynamics, mean-field equations, uncertainty
quantification, stochastic Galerkin, $\mathcal{H}_{\infty}$ control

AMS Classification. 49K15, 49M25, 93-10, 93D09, 35Q70, 35Q93

##  1 Introduction

We consider the mathematical modelling and control of phenomena of collective
dynamics under uncertainties. These phenomena have been studied in several
fields such as socio-economy, biology, and robotics where systems of
interacting particles are given by self-propelled particles, such as animals
and robots, see e.g. [1, 7, 15, 32, 24]. Those particles interact according to
a possibly nonlinear model, encoding various social rules as attraction,
repulsion, and alignment. A particular feature of such models is their rich
dynamical structure, which includes different types of emerging patterns,
including consensus, flocking, and milling [29, 50, 17, 23, 45]. Understanding
the impact of control inputs in such complex systems is of great relevance for
applications. Results in this direction allow to design optimized actions such
as collision-avoidance protocols for swarm robotics [14, 46, 48, 26],
pedestrian evacuation in crowd dynamics [16, 22], supply chain policies [34,
18], the quantification of interventions in traffic management [51, 31, 49] or
in opinion dynamics [27, 28]. Further, the introduction of uncertainty in the
mathematical modelling of real-world phenomena seems to be unavoidable for
applications, since often at most statistical information of the modelling
parameters is available. The latter has typically been estimated from
experiments or derived from heuristic observations [5, 8, 37]. To produce
effective predictions and to describe and understand physical phenomena, we
may incorporate parameters reflecting the uncertainty in the interaction
rules, and/or external disturbances [13].

Here, we are concerned with the robustness of controls influencing the
evolution of a collective motion of an interacting agent system. The controls
we are considering are aimed to stabilize the system’s dynamic under external
uncertainty. From a mathematical point of view, a description of self-
organized models is provided by complex system theory, where the overall
dynamics are depicted by a large-scale system of ordinary differential
equations (ODEs).

More precisely, we consider the control of high-dimensional dynamics
accounting $N$ agents with state
$v_{i}(t,\theta)\in\mathbb{R}^{d},\,i=1,\ldots,N$, evolving according to

$\frac{d}{dt}{v}_{i}(t,\theta)=\sum_{j=1}^{N}a_{ij}(v_{j}(t,\theta)-v_{i}(t,\theta))+u_{i}(t,\theta)+\sum_{k=1}^{Z}\theta_{k}\,,\qquad
v_{i}(0)=v_{i}^{0},$ (1.1)

where $A=[a_{ij}]\in\mathbb{R}^{N\times N}$ defines the nature of pairwise
interaction among agents, and
$\theta=(\theta_{1},\ldots,\theta_{Z})^{\top}\in\mathbb{R}^{Z\times d}$ is a
random input vector with a given probability density distribution on $Z$ as
$\rho\equiv\rho_{1}\otimes\ldots\otimes\rho_{Z}$. The control signal
$u_{i}(t,\theta)\in\mathbb{R}^{d}$ is designed to stabilize the state toward a
target state $\bar{v}\in\mathbb{R}^{N\times d}$, and its action is influenced
by the random parameter $\theta$. This is also due to the fact, that later we
will be interested in closed–loop or feedback controls on the state
$(v_{1}.\dots,v_{N})$ that in turn dependent on the unknown parameter
$\theta.$

Of particular interest will be controls designed via minimization of linear
quadratic (parametric) regulator functional such as

$\min_{u(\cdot,\theta)}{J}(u;v^{0}):=\int_{0}^{+\infty}\exp(-r\tau)\left[v^{\top}Qv+\nu
u^{\top}Ru\right]\,d\tau,$ (1.2)

with $Q$ positive semi-definite matrix of order $N$, $R$ positive definite
matrix of order $N$ and $r$ is a discount factor. In this case, the linear
quadratic dynamics allow for an optimal control $u^{*}$ stabilising the
desired state $v_{d}=0$, expressed in feedback form, and obtained by solving
the associated matrix Riccati -equations. Those aspects will be also addressed
in more detail below.

In order to assess the performances of controls, and quantify their robustness
we propose estimates using the concept of $\mathcal{H}_{\infty}$ control. In
this setting different approaches have been studied in the context of
$\mathcal{H}_{\infty}$ control and applied to first-order and higher-order
multiagent systems, see e.g. [43, 40, 44, 41, 42], in particular for an
interpretation of $\mathcal{H}_{\infty}$ as dynamic games we refer to [6].
Here we will study an approach based on the derivation of sufficient
conditions in terms of linear matrix inequalities (LMIs) for the
$\mathcal{H}_{\infty}$ control problem. In this way, consensus robustness will
be ensured for a general feedback formulation of the control action.
Additionally, we consider the large–agent limit and show that the robustness
is guaranteed independently of the number of agents.

Furthermore, we will discuss the numerical realization of system (1.1)
employing uncertainty quantification techniques. In general, at the numerical
level, techniques for uncertainty quantification can be classified into non-
intrusive and intrusive methods. In a non-intrusive approach, the underlying
model is solved for fixed samples with deterministic schemes, and statistics
of interest are determined by numerical quadrature, typical examples are
Monte-Carlo and stochastic collocation methods [19, 53]. While in the
intrusive case, the dependency of the solution on the stochastic input is
described as a truncated series expansion in terms of orthogonal functions.
Then, a new system is deduced that describes the unknown coefficients in the
expansion. One of the most popular techniques of this type is based on
stochastic Galerkin (SG) methods. In particular, generalized polynomial chaos
(gPC) gained increasing popularity in uncertainty quantification (UQ), for
which spectral convergence on the random field is observed under suitable
regularity assumptions [19, 35, 36, 53]. The methods, here developed, make use
of the stochastic Galerkin (SG) for the microscopic dynamics while in the
mean-field case we combine SG in the random space with a Monte Carlo method in
the physical variables.

The manuscript is organized as follows, in Section 2 we introduce the problem
setting and propose different feedback control laws; in Section 3 we
reformulate the problem in the setting of $\mathcal{H}_{\infty}$ control and
provide conditions for the robustness of the controls in the microscopic and
mean-field case. Section 4 is devoted to the description of numerical
strategies for the simulation of the agent systems, and to different numerical
experiments, which assess the performances and compare different methods.

##  2 Control of interacting agent system with uncertainties

The following notation is introduced with the control of high-dimensional
systems of interacting agents with random inputs. We consider the evolution of
$N$ agents with state $v(t,\theta)\in\mathbb{R}^{N\times d}$ as follows

$\frac{d}{dt}{v}_{i}(t,\theta)=\frac{1}{N}\sum_{j=1}^{N}\bar{p}(v_{j}(t,\theta)-v_{i}(t,\theta))+u_{i}(t,\theta)+\sum_{k=1}^{Z}\theta_{k}\,$
(2.1)

with deterministic initial data $v_{i}(0)=v_{i}^{0}$ for $i=1,\ldots,N$, and
where $\theta_{k}\in\Omega_{k}\subseteq\mathbb{R}^{d}$ for $k=1,\ldots,Z$ are
random inputs, distributed according to a compactly supported probability
density $\rho\equiv\rho_{1}\otimes\dots\otimes\rho_{Z}$, i.e.,
$\rho_{k}(\theta)\geq 0$ a.e., $\textrm{supp}(\rho_{k})\subseteq\Omega_{k}$
and $\int_{\Omega_{k}}\rho_{k}(\theta)\,d\theta=1$. For simplicity, we also
assume that the random inputs have zero average $\mathbb{E}[\theta_{k}]=0$.
The control signal $u(t,\theta)\in\mathbb{R}^{N\times d}$ is designed
minimizing the (parameterized) objective

$u^{*}(\cdot,\theta)=\arg\min_{u(\cdot,\theta)}{J}(u;v^{0}):=\int_{0}^{+\infty}\exp(-r\tau)\left(\frac{1}{N}\sum_{j=1}^{N}(|v_{j}(\tau,\theta)-\bar{v}|^{2}+\nu|u_{j}(\tau,\theta)|^{2})\right)\,d\tau,$
(2.2)

with $\nu>0$ being a penalization parameter for the control energy, the norm
$|\cdot|$ being the usual Euclidean norm in $\mathbb{R}^{d}$. The discount
factor $\exp(-r\tau)$ is introduced to have a well-posed integral.

We assume that $\bar{v}$ is a prescribed consensus point, namely, in the
context of this work we are interested in reaching a consensus velocity
$\bar{v}\in\mathbb{R}^{d}$ such that $v_{1}=\ldots=v_{N}=\bar{v}$, and
w.l.o.g. we can assume $\bar{v}=0$. Note that $\bar{v}=0$ is also the steady
state of the dynamics in absence of disturbances. Hence, we may view
$u(\cdot,\theta)$ as a stabilizing control of the zero steady state of the
system. Furthermore, we will be interested in feedback controls $u.$

Recall that the (deterministic) linear model (2.1), without uncertainties,
allows a feedback stabilization by solving the resulting optimal control
problem through a Riccati equations [33, 3, 2]. The functional $J$ in (2.2),
in absence of disturbances, reads as follows

$\displaystyle J(u;v^{0})=\int_{0}^{+\infty}\exp(-r\tau)\left(v^{\top}Qv+\nu
u^{\top}Ru\right)\,dt$

where $Q\equiv R=\frac{1}{N}\textrm{Id}_{N}$. In this case the controlled
dynamics (2.1) is reformulated in a matrix-vector notation

$\frac{d}{dt}v(t)=Av(t)+Bu(t),\qquad u(t)=-\frac{N}{\nu}Kv(t),$ (2.3)

with $B=\textrm{Id}_{N}$ the identity matrix of order $N$, and

$\displaystyle(A)_{ij}=\begin{cases}&a_{d}=\frac{\bar{p}(1-N)}{N},\qquad
i=j,\\\ &a_{o}=\frac{\bar{p}}{N},\qquad\qquad i\neq j.\\\ \end{cases}$ (2.4)

The matrix $K$ associated to feedback form of the optimal control has to
fulfilll the Riccati matrix-equation of the following form

$0=-rK+KA+A^{\top}K-\frac{N}{\nu}KK+Q.$ (2.5)

For a general linear system, we need to solve the $N\times N$ equations to
find $K$, which can be costly for large-scale agent-based dynamics. However,
we can use the same argument of [3] and exploit the symmetric structure of the
Laplacian matrix $A$ to reduce the algebraic Riccati equation. Unlike in [3],
where they investigate the case with finite terminal time, here we state the
following proposition for the infinite horizon case with discount factor $r$

######  Proposition 2.1 (Properties of the Algebraic Riccati Equation (ARE)).

For the linear dynamics (2.3), the solution of the Riccati equation (2.5)
reduces to the solution of

$\begin{split}0&=-rk_{d}-\frac{2\bar{p}(N-1)}{N}({k}_{d}-k_{o})-\frac{N}{\nu}\left({k}_{d}^{2}+(N-1){k}_{o}^{2}\right)+\frac{1}{N},\,\\\
0&=-rk_{o}+\frac{2\bar{p}}{N}(k_{d}-k_{o})-\frac{N}{\nu}\left(2{k}_{d}{k}_{o}+(N-2){k}_{o}^{2}\right).\end{split}$
(2.6)

The entries $(i,j)$ of the matrix $K$ of the algebraic Riccati equation (2.5)
is given by

$(K)_{ij}=\delta_{ij}k_{d}+(1-\delta_{ij})k_{o}.$

In order to allow the limit of infinitely many agents $N\to\infty$, we
introduce the following scalings

${k}_{d}\leftarrow Nk_{d},\quad{k}_{o}\leftarrow
N^{2}k_{o},\quad\alpha(N)=\frac{N-1}{N},$

and keeping the same notation also for the scaled variables $k_{d},k_{o}$, the
system (2.6) reads

$\begin{split}0&=-rk_{d}-2\bar{p}\alpha(N)\left(k_{d}-\frac{k_{o}}{N}\right)-\frac{1}{\nu}\left(k_{d}^{2}+\frac{\alpha(N)}{N}k_{o}^{2}\right)+1,\\\
0&=-rk_{o}+2\bar{p}\left(k_{d}-\frac{k_{o}}{N}\right)-\frac{1}{\nu}\left(2k_{d}k_{o}+\alpha(N)k_{o}^{2}-\frac{1}{N}k_{o}^{2}\right).\end{split}$
(2.7)

The previous considerations motivate to extend formula (2.3) to the parametric
case (2.1). Hence, in the presence of parametric uncertainty the feedback
control is written explicitly as follows

$u_{i}(t,\theta)=-\frac{1}{\nu}\left(Kv(t,\theta)\right)_{i}=-\frac{1}{\nu}\left(\left(k_{d}-\frac{k_{o}}{N}\right)v_{i}(t,\theta)+\frac{k_{o}}{N}\sum_{j=1}^{N}v_{j}(t,\theta)\right).$
(2.8)

The question arises if the given feedback is robust with respect to the
uncertainties $\theta$. In the following, we will provide a measure for the
robustness of control (2.8) in the framework of $\mathcal{H}_{\infty}$
control. Some additional remarks allow generalizing formula (2.8).

######  Remark 2.1 (Non-zero average).

In presence of general uncertainties with known expectations, we modify the
control (2.8) for model (2.1) including a correction factor given by the
expected values of the random inputs,

$u_{i}(t,\theta)=-\frac{1}{\nu}\left(\left(k_{d}-\frac{k_{o}}{N}\right)v_{i}(t,\theta)+\frac{k_{o}}{N}\sum_{j=1}^{N}v_{j}(t,\theta)\right)-\sum_{k=1}^{Z}\mu_{k},$
(2.9)

for $\mu_{k}=\mathbb{E}[\theta_{k}]$ for $k=1,\ldots,Z$.

######  Remark 2.2 (Averaged control).

In the case of a deterministic feedback control, we may consider the
expectation of the objective (2.2) subject to the noisy model (2.1)

$\bar{u}^{*}(\cdot)=\arg\min_{u(\cdot)}\mathbb{E}\left[\int_{0}^{+\infty}\frac{1}{2}(v^{\top}Qv+\nu
u^{\top}Ru)\ dt\right],$

where we introduce the matrices $Q=R=\frac{1}{N}\textrm{Id}_{N}$. In this
case, we have the following deterministic optimal feedback control is deduced

$\bar{u}_{i}(t)=-\frac{1}{\nu}\left({k}_{d}\mathbb{E}\left[v_{i}(t,\theta)\right]+\frac{{k}_{o}}{N}\sum_{j\neq
i}^{N}\mathbb{E}\left[v_{j}(t,\theta)\right]\right)-\sum_{k=1}^{Z}\mu_{k},$
(2.10)

where $\mu_{k}=\mathbb{E}[\theta_{k}]$ for $k=1,\ldots,Z$, ${k}_{d}$,
${k}_{o}$ satisfy equations (2.7). We refer to Appendix A for detail
computations for the synthesis of (2.10).

##  3 Robustness in the $\mathcal{H}_{\infty}$ setting

In the context of $\mathcal{H}_{\infty}$ theory, controllers are synthesized
to achieve stabilization with guaranteed performance. In this section, we
exploit the theory of Linear Matrix Inequality (LMI) to show the robustness of
the control. The introduction of LMI methods in control has dramatically
expanded the types and complexity of the systems we can control. In
particular, it is possible to use LMI solvers to synthesize optimal or
suboptimal controllers and estimators for multiple classes of state-space
systems and without giving a complete list of references, we refer to [47, 38,
9, 20].

Consider the linear system (2.1) with control (2.8) in the following
reformulation

$\displaystyle\frac{d}{dt}{v}(t)$
$\displaystyle=\widehat{A}v(t)+\widehat{B}\theta,$ (3.1)

where we consider the random input vector
$\theta=(\theta_{1},\ldots,\theta_{Z})^{\top}\in\mathbb{R}^{Z\times d}$, and
the matrices

$\widehat{A}=A-\frac{1}{\nu}{K},\qquad\widehat{B}=\mathds{1}_{N\times Z}.$

with $\mathds{1}$ a matrix of ones of dimension $N\times Z$. We introduce the
frequency transfer function
$\hat{G}(s):=(s\textrm{Id}_{N}-\widehat{A})^{-1}\widehat{B}$, such that
$\hat{G}\in R\mathcal{H}_{\infty}$. The latter is the set of proper rational
functions with no poles in the closed complex right half-plane, and the signal
norm $\|\cdot\|_{\mathcal{H}_{\infty}}$ measuring the size of the transfer
function in the following sense

$\|\hat{G}\|_{\mathcal{H}_{\infty}}=\text{ess}\sup_{\omega\in\mathbb{R}}\bar{\sigma}(\hat{G}(i\omega)),$
(3.2)

where for a given matrix $P$, $\bar{\sigma}(P)$ is the largest singular value
of $P$. The general $\mathcal{H}_{\infty}$-optimal control problem consists of
finding a stabilizing feedback controller $u=-\frac{1}{\nu}K$ which minimizes
the cost function (3.2) and we refer to Appendix B and to [21] for more
details.

However, the direct minimization of the cost
$\|\hat{G}\|_{\mathcal{H}_{\infty}}$ is in general a very hard task, and
possibly unfeasible by direct methods. To reduce the complexity, a possibility
consists in finding conditions for the stabilizing controller that achieves a
norm bound for a given threshold $\gamma>0$,

$\|\hat{G}\|_{\mathcal{H}_{\infty}}\leq\gamma.$ (3.3)

Hence, robustness of a given control $u=-\frac{1}{\nu}K$ is measured in terms
of the smallest $\gamma$ satisfying (3.3). In order to provide a quantitative
result, we can rely on the following result

######  Lemma 3.1.

Given the frequency transfer function $\hat{G}$, associated to (3.1), a
necessary and sufficient condition to guarantee the $\mathcal{H}_{\infty}$
bound (3.3) for $\gamma>0$, is to prove, that there exists a positive definite
square matrix $X$ of order $N$, $X>0$, such that the following algebraic
Riccati equation holds

$A^{\top}X+XA-\frac{1}{\nu}{K}^{\top}X-\frac{1}{\nu}X{K}+\frac{1}{\gamma}XX+\frac{1}{\gamma}\textrm{Id}_{N}=\mathcal{O}\left(\frac{1}{N}\right).$
(3.4)

For detailed proof of this result, we refer to Appendix B. Note that we are
interested in the large particle limit and hence may allow that the previous
equation is not exactly zero, but tends to zero at a rate $1/N.$

######  Theorem 3.2.

Consider system (3.1) with structure induced by (2.1), and consider a square
matrix $X$ with the following structure

$\displaystyle(X)_{ij}=\begin{cases}x_{d},\qquad i=j,\\\ x_{o},\qquad i\neq
j.\end{cases}\qquad$

Then for $N$ sufficiently large and finite, control $u=-\frac{1}{\nu}K$ given
by equation (2.8) is $\mathcal{H}_{\infty}$-robust for any $\gamma$ and
$c_{N}$ such that $\gamma\geq\frac{1}{c_{N}}$, $c_{N}>0$, where

$\displaystyle c_{N}=\bar{p}+\frac{1}{\nu}(k_{d}-\frac{k_{o}}{N}).$ (3.5)

###### Proof.

Under the hypothesis of the theorem, (3.4) reduces to the following system of
equations

$\displaystyle\mathcal{O}\left(\frac{1}{N}\right)$
$\displaystyle=\frac{2\bar{p}(1-N)}{N}x_{d}+\frac{2\bar{p}(N-1)}{N}x_{o}-\frac{2}{\nu}{k}_{d}x_{d}-\frac{2(N-1)}{\nu
N}{k}_{o}x_{o}$ (3.6)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{\gamma}x_{d}^{2}+\frac{N-1}{\gamma}x_{o}^{2}+\frac{1}{\gamma},$
(3.7) $\displaystyle\mathcal{O}\left(\frac{1}{N}\right)$
$\displaystyle=\frac{2\bar{p}(1-N)}{N}x_{o}+\frac{2\bar{p}}{N}x_{d}+\frac{2\bar{p}(N-2)}{N}x_{o}-\frac{2}{\nu}{k}_{d}x_{o}-\frac{2}{\nu
N}{k}_{o}x_{d}$ (3.8)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{2(N-2)}{\nu
N}{k}_{o}x_{o}+\frac{2}{\gamma}x_{d}x_{o}+\frac{N-2}{\gamma}x_{o}^{2}.$ (3.9)

We scale the off-diagonal elements $x_{o}$ of $X$ according to

$\quad\tilde{x}_{o}=\sqrt{N}x_{o},$

which, as we will see later, is also the consistent scaling for a mean-field
description of $X.$

Then, the previous system reads

$\displaystyle\begin{split}\mathcal{O}\left(\frac{1}{N}\right)&=\frac{2\bar{p}(1-N)}{N}x_{d}+\frac{2\bar{p}(N-1)}{N\sqrt{N}}\tilde{x}_{o}-\frac{2}{\nu}k_{d}x_{d}-\frac{2(N-1)}{\nu
N\sqrt{N}}k_{o}\tilde{x}_{o}+\frac{1}{\gamma}x_{d}^{2}+\frac{N-1}{\gamma
N}\tilde{x}_{o}^{2}+\frac{1}{\gamma},\\\
\mathcal{O}\left(\frac{1}{N}\right)&=\frac{2\bar{p}(1-N)}{N\sqrt{N}}\tilde{x}_{o}+\frac{2\bar{p}}{N}x_{d}+\frac{2\bar{p}(N-2)}{N\sqrt{N}}\tilde{x}_{o}-\frac{2}{\nu\sqrt{N}}k_{d}\tilde{x}_{o}-\frac{2}{\nu
N}k_{o}x_{d}-\frac{2(N-2)}{\nu
N\sqrt{N}}k_{o}\tilde{x}_{o}\cr&\qquad\qquad+\frac{2}{\gamma\sqrt{N}}x_{d}\tilde{x}_{o}+\frac{N-2}{\gamma
N}\tilde{x}_{o}^{2}.\end{split}$ (3.10)

From these two equations of system (3.10), and setting

$c=\bar{p}+\frac{k_{d}}{\nu},\quad\alpha=\bar{p}-\frac{k_{o}}{\nu},\quad\beta=\frac{k_{d}+k_{o}}{\nu}.$
(3.11)

we obtain two second order equations for $x_{d}$ and $\tilde{x}_{o}$

$\begin{split}0&=x_{d}^{2}-2\gamma
cx_{d}+\tilde{x}_{o}^{2}+\frac{2\gamma\alpha}{\sqrt{N}}\tilde{x}_{o}+1+\mathcal{O}\left(\frac{1}{N}\right),\\\
0&=\tilde{x}_{o}^{2}-\frac{2\gamma}{\sqrt{N}}\Bigl{(}\beta-\frac{x_{d}}{\gamma}\Bigr{)}\tilde{x}_{o}+\mathcal{O}\left(\frac{1}{N}\right).\end{split}$
(3.12)

For $N$ sufficiently large, their solutions is given by

$\displaystyle x_{d}^{\pm}$ $\displaystyle=\gamma
c\pm\sqrt{\gamma^{2}c^{2}-1-\tilde{x}_{o}^{2}-\frac{2\gamma\alpha}{\sqrt{N}}\tilde{x}_{o}},\quad\tilde{x}_{o}^{-}=0,\quad\tilde{x}_{o}^{+}=\frac{2\gamma}{\sqrt{N}}\Bigl{(}\beta-\frac{x_{d}}{\gamma}\Bigr{)}.$

Hence, we write the matrix $X$ as

$X=\frac{\tilde{x}_{o}^{\pm}}{\sqrt{N}}\
\mathds{1}_{N}+(x_{d}^{\pm}-\frac{\tilde{x}_{o}^{\pm}}{\sqrt{N}})\ {Id}_{N},$

and the eigenvalues of $X$ are

$\lambda_{i}=\lambda=x_{d}^{\pm}-\frac{\tilde{x}_{o}^{\pm}}{\sqrt{N}}\quad\text{for}\quad
i=1,\dots,N-1\quad\text{and}\quad\lambda_{N}=x_{d}^{\pm}+(N-1)\frac{\tilde{x}_{o}^{\pm}}{\sqrt{N}}.$

If we subtract equations in (3.10) we find the following second order equation
in the variable $\lambda=(x_{d}-\frac{\tilde{x}_{o}}{\sqrt{N}})$

$\lambda^{2}-2\gamma\Bigl{(}\bar{p}+\frac{1}{\nu}({k}_{d}-\frac{k_{o}}{N})\Bigr{)}\lambda+1=0,$

with its solutions $\lambda^{\pm}=x_{d}-\frac{\tilde{x}_{o}}{\sqrt{N}}=\gamma
c_{N}\pm\sqrt{\gamma^{2}c_{N}^{2}-1},$ and for

$c_{N}=\bar{p}+\frac{1}{\nu}\left(k_{d}-\frac{k_{o}}{N}\right).$ (3.13)

Therefore, we have $\lambda_{i}=\lambda^{\pm}\quad\text{for}\quad
i=1,\dots,N-1$ and
$\lambda_{N}=\lambda^{\pm}+N\frac{\tilde{x}_{o}}{\sqrt{N}}.$ Hence, there
exists a matrix $X$ satisfying (3.4), for

$x_{o}=x_{o}^{-}=0.$

In this case, the eigenvalues are $\lambda_{i}=\lambda_{N}=\lambda^{\pm}$. In
order to ensure the positivity of the eigenvalue $\lambda^{\pm}=\gamma
c_{N}\pm\sqrt{\gamma^{2}c_{N}^{2}-1}$ needs to be non–negative. This can be
expressed in terms of the choice of the parameters $\gamma$ and $c_{N}:$ First
of all, we have $\gamma\geq\frac{1}{c_{N}}$ to ensure the existence of the
square root. In addition, provided that $c_{N}>0$ we obtain

$\displaystyle\lambda^{+}>0\mbox{ and }\lambda^{-}>0.$

This finishes the proof. ∎

######  Remark 3.1.

We observe that Theorem 3.2 quantifies the robustness of the feedback control
with a lower bound on the parameters of the model. In particular, we can
achieve minimal value of $\gamma$ for large values of $c_{N}$ in (3.5), for
example if we fix the values $\bar{p},k_{d},k_{o}$ and $N$ for decreasing
values of the penalization $\nu$ the control is more robust.

###  3.1 Mean-field estimates for $\mathcal{H}_{\infty}$ control

Large system of interacting agents can be efficiently represented at different
scales to reduce the complexity of the microscopic description. Of particular
interest are models able to describe the evolution of a density of agents and
its moments [10, 12, 30].

In this section, we analyse robustness of controls in the case of a large
number of agents, i.e. $N\gg 1$, by means of the mean-field limit of the
interacting system. Hence, we consider the density distribution of agents
$f=f(t,v,\theta)$ to describe the collective behaviour of the ensemble of
agents. The empirical joint probability distribution of agents for the system
(2.1), is given by

$f^{N}(t,v,\theta)=\frac{1}{N}\sum_{i=1}^{N}\delta(v-v_{i}(t,\theta)),$

where $\delta(\cdot)$ is a Dirac measure over the trajectories
$v_{i}(t,\theta)$ dependent on the stochastic variable
$\theta=(\theta_{1},\ldots,\theta_{Z})$.

Hence, assuming enough regularity assuming that agents remain in a fixed
compact domain for all $N$ and in the whole time interval $[0,T]$, the mean-
field limit of dynamics (2.1) is obtained formally as

$\partial_{t}f(t,v,\theta)=-\nabla_{v}\cdot\left(f(t,v,\theta)\left(\left(\bar{p}-\dfrac{k_{o}}{\nu}\right)m_{1}[f](t,\theta)-\left(\bar{p}+\frac{k_{d}}{\nu}\right)v+\sum_{k=1}^{Z}\theta_{k}\right)\right),$

with initial data $f(0,v,\theta)=f^{0}(v,\theta)$.

The latter is obtained as the limit of $f^{N}(0,v,\theta)$ in the Wasserstein
distance given a sequence of initial agents, [11, 10]. The quantity $m_{1}[f]$
denotes the first moment of $f$ with respect to $v$

$\displaystyle m_{1}[f](t,\theta)=\int_{\mathbb{R}^{d}}vf(t,v,\theta)dv.$

For the many-particle limit, we recover a mean-field estimate of the
$\mathcal{H}_{\infty}$ condition similarly as in 3.2. Indeed, for
$N\rightarrow\infty$ the nonlinear system (3.10) yields

$\displaystyle 0$
$\displaystyle=-2\bar{p}x_{d}-\frac{2}{\nu}k_{d}x_{d}+\frac{1}{\gamma}x_{d}^{2}+\frac{1}{\gamma}x_{o}^{2}+\frac{1}{\gamma},\;0=\frac{1}{\gamma}x_{o}^{2}.$

Hence, for any fixed finite $N$ the matrix $X$ is diagonal with the entry

$x_{d}^{\pm}=\gamma\Bigl{(}\bar{p}+\frac{k_{d}}{\nu}\Bigr{)}\pm\sqrt{\gamma^{2}\Bigl{(}\bar{p}+\frac{k_{d}}{\nu}\Bigr{)}^{2}-1}+O(\frac{1}{N}).$

To ensure that $X$ is positive definite, we only have to assume that
$\gamma\geq\dfrac{1}{c},$ where $\gamma$ is the bound of the
$\mathcal{H}_{\infty}$ norm, and $c=\bar{p}+k_{d}/\nu+O(\frac{1}{N})$
corresponds to the value defined by equation (3.13). This shows that for any
$N$ there exists a positive definite matrix that guarantees robust
stabilization. Note that the condition for any $N$ is the limit of the finite-
dimensional conditions of the previous Lemma 3.1 and Theorem 3.2.

######  Remark 3.2.

For explicit values of the Riccati coefficients we can characterize the
previous estimates more precisely. In particular, for $N\rightarrow\infty$
system (2.7) reduces to

$\displaystyle 0=\frac{k_{d}^{2}}{\nu}+\left(2\bar{p}+r\right)k_{d}-1,\qquad
0=\frac{k_{o}^{2}}{\nu}+\left(\frac{2}{\nu}k_{d}+r\right)k_{o}-2\bar{p}k_{d},$

with solutions

$k_{d}^{\pm}=-\nu\left(\bar{p}+\frac{r}{2}\right)\pm\nu\sqrt{\left(\bar{p}+\frac{r}{2}\right)^{2}+\frac{1}{\nu}},\quad\quad
k_{o}^{\pm}=-\left(k_{d}+\frac{\nu r}{2}\right)\pm\sqrt{\left(k_{d}+\frac{\nu
r}{2}\right)^{2}+2\nu\bar{p}k_{d}}.$

In this particular case, the condition of Theorem 3.2 becomes

$\gamma\geq\frac{\sqrt{\nu}}{{\sqrt{\left(\bar{p}+\frac{r}{2}\right)^{2}\nu+1}}-\frac{r\sqrt{\nu}}{2}}.$
(3.14)

In Figure 3.1 we depict the lower bound of $\gamma$ for $r=0$ for different
values of $\nu$ and $\bar{p}$. As expected, smaller values of $\gamma$, hence
more robustness, is obtained if the penalization factor $\nu$ is small or when
$\bar{p}$ is large. The latter corresponds to a stronger attraction between
agents.

Figure 3.1: Numerical computation of the lower bound value of $\gamma$ in
(3.14) as function of $\nu$ and $\bar{p}$ when $r=0$, i.e.
$\gamma(\nu,\bar{p})=\sqrt{{\nu}/(\bar{p}^{2}\nu+1)}.$

##  4 Numerical approximation of the uncertain dynamics

In this section, we present numerical tests based on linear microscopic and
mean-field equations in presence of uncertainties. In particular, we give
numerical evidence of the robustness of the feedback control (2.8) and
illustrate a comparison with the averaged control (A.3). For the numerical
approximation of the random space, we employ the Stochastic Galerkin (SG)
method belonging to the class of generalized polynomial chaos (gPC) ([39,
53]). In the mean-field setting, the evolution of the density distribution is
approximated with a Monte Carlo (MC) method, in a similar spirt of particle
based gPC techniques developed in [13].

###  4.1 SG approximation for robust constrained interacting agent systems

We approximate the dynamics using a stochastic Galerkin approach applied to
the interacting particle system with uncertainties [4, 13]. Polynomial chaos
expansion provides a way to represent a random variable with finite variance
as a function of an $M$-dimensional random vector using a polynomial basis
that is orthogonal to the distribution of this random vector. Depending on the
distribution, different expansion types are distinguished, as shown in Table
4.1.

Probability law of $\theta$ | Expansion polynomials | Support
---|---|---
Gaussian | Hermite | $(-\infty,+\infty)$
Uniform | Legendre | $[a,b]$
Beta | Jacobi | $[a,b]$
Gamma | Laguerre | $[0,+\infty)$
Poisson | Charlier | $\mathbb{N}$
Table 4.1: The different choices for the polynomial expansions.

We recall first some basic notions on gPC approximation techniques and for the
sake of simplicity we consider a one-dimensional setting for the dynamical
state $v_{i},$ i.e., $d=1$. Let $(\Omega,\mathcal{F},P)$ be a probability
space where $\Omega$ is an abstract sample space, $\mathcal{F}$ a
$\sigma-$algebra of subsets of $\Omega$ and $P$ a probability measure on
$\mathcal{F}$. Let us define a random variable

$\theta:(\Omega,\mathcal{F})\rightarrow(I_{\Theta},\mathcal{B}(\mathbb{R}^{Z}))$
(4.1)

where $I_{\Theta}\in\mathbb{R}^{Z}$ is the range of $\theta$ and
$\mathcal{B}(\mathbb{R}^{Z})$ is the Borel $\sigma$-algebra of subsets of
$\mathbb{R}^{Z}$, we recall that $Z$ is the dimension of the random input
$\theta=(\theta_{1},\ldots,\theta_{Z})$ and where we assume that each
component is independent. We consider the linear spaces generated by
orthogonal polynomials of $\theta_{j}$ with degree up to $M$:
$\\{\Phi_{k_{j}}^{(j)}(\theta)\\}_{k_{j}=0}^{M}$, with $j=1,\dots,Z$. Assuming
that the probability law for the function $v_{i}(t,\theta)$ has a finite
second order moment, the complete polynomial chaos expansion of $v_{i}$ is
given by

$v_{i}(t,\theta)=\sum_{k_{1},\ldots,k_{Z}\in\mathbb{N}}\hat{v}_{i,k_{1}\dots
k_{Z}}(t)\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j}),$

where the coefficients $\hat{v}_{i,k_{1}\dots k_{Z}}(t)$ are defined as

$\hat{v}_{i,k_{1}\dots
k_{Z}}(t)=\mathbb{E}_{\theta}\left[v_{i}(t,\theta)\prod_{j=1}^{Z}\Phi^{(j)}_{k_{j}}(\theta_{j})\right]$

where the expectation operator $\mathbb{E}_{\theta}$ is computed with respect
to the joint distribution $\rho=\rho_{1}\otimes\ldots\otimes\rho_{Z}$, and
where $\\{\Phi^{(j)}_{k}(\theta_{j})\\}_{k}$ is a set of polynomials which
constitute the optimal basis with respect to the known distribution
$\rho(\theta_{j})$ of the random variable $\theta_{j}$, such that

$\displaystyle\mathbb{E}_{\theta_{j}}\left[\Phi^{(j)}_{k}(\theta_{1})\Phi^{(j)}_{h}(\theta_{j})\right]$
$\displaystyle=\mathbb{E}_{\theta_{j}}\left[\Phi^{(j)}_{h}(\theta_{j})^{2}\right]\delta_{hk},$

with $\delta_{hk}$ the Kronecker delta. From the numerical point of view, we
may have an exponential order of convergence for the SG series expansion,
unlike Monte Carlo techniques for which the order is $\mathcal{O}(1/\sqrt{M})$
where $M$ is the number of samples. Considering the noisy model 2.1 with
control $u_{i}(t)$ in 2.8, we have

$\dot{v}_{i}(t,\theta)=\frac{\bar{p}}{N}\sum_{j=1}^{N}\Bigl{(}v_{j}(t,\theta)-v_{i}(t,\theta)\Bigr{)}-\frac{1}{N\nu}\sum_{j=1}^{N}\Bigl{(}k_{o}(t)v_{j}(t,\theta)+k_{d}(t)v_{i}(t,\theta)\Bigr{)}+\sum_{j=1}^{Z}\theta_{j}.$
(4.2)

We apply the SG decomposition to the solution of the differential equation
$v_{i}(t,\theta)$ in (4.2) and to the stochastic variable $\theta_{k}$, and
for $i=1,\dots,N,\ l=1,\dots,Z$, we have

$\begin{split}v_{i}^{M}(t,\theta)&=\sum_{k_{1},\ldots,k_{Z}=0}^{M}\hat{v}_{i,k_{1}\dots
k_{Z}}(t)\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j}),\\\
\theta_{l}^{M}(\theta)&=\sum_{k_{1},\ldots,k_{Z}=0}^{M}\hat{\theta}_{l,k_{1}\dots
k_{Z}}(t)\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j}),\end{split}$ (4.3)

with

$\hat{\theta}_{l,k_{1}\dots
k_{Z}}=\mathbb{E}_{\theta}\left[\theta_{l}\prod_{j=1}^{Z}\Phi^{(j)}_{k_{j}}(\theta_{j})\right]=\mathbb{E}_{\theta_{l}}\left[\theta_{l}\Phi^{(l)}_{k_{l}}(\theta_{l})\right]\prod_{j=1,j\neq
l}^{Z}\mathbb{E}_{\theta_{j}}\left[\Phi^{(j)}_{k_{j}}(\theta_{j})\right].$

Then we obtain the following polynomial chaos expansion

$\begin{split}&\frac{d}{dt}\sum_{k_{1},\ldots,k_{Z}=0}^{M}\hat{v}_{i,k_{1}\dots
k_{Z}}\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j})=\\\
&=\frac{1}{N}\sum_{h=1}^{N}\sum_{k_{1},\ldots,k_{Z}=0}^{M}\left[\left(\bar{p}-\frac{k_{o}}{\nu}\right)\hat{v}_{h,k_{1}\dots
k_{Z}}-\left(\bar{p}+\frac{k_{d}}{\nu}\right)\hat{v}_{i,k_{1}\dots
k_{Z}}\right]\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j})\\\
&\quad+\sum_{l=1}^{Z}\sum_{k_{1},\ldots,k_{Z}=0}^{M}\hat{\theta}_{l,k_{1}\dots
k_{Z}}(t)\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j}).\end{split}$ (4.4)

Multiplying by $\prod_{j=1}^{Z}\Phi_{k_{j}}^{(j)}(\theta_{j})$ and integrating
with respect to the distribution $\rho(\theta)$, we end up with

$\displaystyle\frac{d}{dt}\hat{v}_{i,k_{1}\dots k_{Z}}=\frac{1}{N}$
$\displaystyle\sum_{h=1}^{N}\left[\left(\bar{p}-\frac{k_{o}}{\nu}\right)\hat{v}_{h,k_{1}\dots
k_{Z}}-\left(\bar{p}+\frac{k_{d}}{\nu}\right)\hat{v}_{i,k_{1}\dots
k_{Z}}\right]$ (4.5)
$\displaystyle+\sum_{l=1}^{Z}\mathbb{E}_{\theta_{l}}\left[\theta_{l}\Phi^{(l)}_{k_{l}}(\theta_{l})\right]\prod_{j=1,j\neq
l}^{Z}\mathbb{E}_{\theta_{l}}\left[\Phi^{(j)}_{k_{j}}(\theta_{j})\right].$
(4.6)

For the numerical tests, we approximate the integrals using quadrature rules.

######  Remark 4.1.

For model 2.10, where the control is averaged with respect to the random
sources,

$u_{i}=-\frac{1}{\nu}\left({k}_{d}\mathbb{E}_{\theta}\left[v_{i}\right]+\frac{{k}_{o}}{N}\sum_{h\neq
i}^{N}\mathbb{E}_{\theta}\left[v_{h}\right]+{s}\sum_{j=1}^{Z}\mathbb{E}_{\theta}\left[\theta_{j}\right]\right),$
(4.7)

the SG approximation is given by

$\begin{split}\frac{d}{dt}\hat{v}_{i,k_{1}\dots
k_{Z}}&=\frac{\bar{p}}{N}\sum_{h=1}^{N}\left(\hat{v}_{h,k_{1}\dots
k_{Z}}-\hat{v}_{i,k_{1}\dots k_{Z}}\right)\\\
-&\frac{\prod_{j=1}^{Z}\mathbb{E}_{\theta_{j}}\left[\Phi^{(j)}_{k_{j}}(\theta_{j})\right]}{\prod_{j=1}^{Z}\mathbb{E}_{\theta_{j}}\left[\left(\Phi^{(j)}_{k_{j}}(\theta_{j})\right)^{2}\right]}\left({k}_{d}\hat{v}_{i,00\dots
0}+\frac{{k}_{o}}{N}\sum_{h\neq i}^{N}{k}_{o}\hat{v}_{h,00\dots
0}+s\sum_{j=1}^{Z}\mu_{j}\right)\\\
+&\sum_{l=1}^{Z}\mathbb{E}_{\theta_{l}}\left[\theta_{l}\Phi^{(l)}_{k_{l}}(\theta_{l})\right]\prod_{j=1,j\neq
l}^{Z}\mathbb{E}_{\theta_{j}}\left[\Phi^{(j)}_{k_{j}}(\theta_{j})\right].\end{split}$
(4.8)

We recover the mean and the variance of the random variable $v(\theta)$ as

$\displaystyle\mathbb{E}_{\theta}[v_{i}(\theta)]$
$\displaystyle=\int_{\mathbb{R}^{2}}v_{i}(\theta)d\rho=\hat{v}_{i,00\dots 0},$
$\displaystyle\mathbb{V}_{\theta}[v_{i}(\theta)]$
$\displaystyle=\int_{\mathbb{R}^{2}}(v_{i}(\theta)-\hat{v}_{i,00\dots
0})^{2}d\rho=\sum_{k_{1},\ldots,k_{Z}=0}^{M}\hat{v}_{i,k_{1}\dots
k_{Z}}^{2}-\hat{v}_{i,00\dots 0}^{2}.$

###  4.2 Numerical tests

In this section we present different numerical tests on microscopic and mean-
field dynamics, to compare the robustness of controls described in sections 2.
We analyze one- and two-dimensional dynamics, for every test we consider the
attractive case with $\bar{p}=1$. The initial distribution of agents $v_{0}$
is chosen such that consensus to the target $\bar{v}=0$ would not be reached
without control action. We implement the SG approximations in (4.4) and (4.8)
and we perform the time integration until the final time $T=1$ of the
resulting system through a 4th order Runge-Kutta method.

We are taking into account a dynamics with $Z=2$ additive uncertanties,
$\theta_{1}$ a random variable with gaussian density distribution
$\rho_{1}\sim\mathcal{N}(\mu,\sigma^{2})$, and $\theta_{2}$ with uniform
density distribution $\rho_{2}\sim\mathcal{U}(a,b)$. This assumption of normal
and uniform distributions for the stochastic parameter corresponds to the case
of Hermite and Legendre polynomial chaos expansions, respectively, as shown in
Table 4.1. For every test, we have $M=10$ terms of the SG decomposition.

####  4.2.1 Test 1: one-dimensional consensus dynamics

In the one-dimensional microscopic case we take $N=100$ agents, and a uniform
initial distribution of agents, $v_{0}\sim\mathcal{U}(10,20)$. Figure 4.2
shows means, as continuous and dashed lines, and confidence regions of the two
noisy dynamics for different distributions $\rho_{1},\ \rho_{2}$. The shaded
region is computed as the region between the values

$\frac{1}{N}\sum_{i=1}^{N}\left(\mathbb{E}_{\theta}[v_{i}(\theta)]\right)\
\mp\
\max_{i=1,\dots,N}\left(\sqrt{\mathbb{V}_{\theta}[v_{i}(\theta)]}\right).$

Figure 4.2: Test 1. Comparison between the two controls in (2.8) and (2.10)
applied to the multiagent system with uncertainties, in terms of dynamics mean
and variance for different $\rho_{1},\rho_{2}$.

Numerical results show that both introduced controls are capable to drive the
agents to the desired state even in the case of a dynamic dependent on random
inputs. Moreover, we can observe that, with the $\mathcal{H}_{\infty}$
control, the variance of the uncertain dynamics is stabilized over time, while
in the case of a averaged control the variance keeps growing. This is because
the averaged control has information only on the mean value of the state and
uncertainty, while the $\mathcal{H}_{\infty}$ feedback control directly
depends on the state, and as a consequence on the randomness of the dynamics.
This is also expected given the robustness estimate on the feedback control.

####  4.2.2 Test 2: two-dimensional consensus dynamics

In the two-dimensional microscopic case, we take $N=100$ agents, and an
initial configuration is uniformly distributed on a 2D disc, as shown in
Figure 4.3. 2D means and confidence regions of the two noisy dynamics can be
seen in Figure 4.4, for different values of the penalization factor $\nu$ and
different distributions $\rho_{1},\ \rho_{2}$.

Figure 4.3: Test 2: Initial distribution of the agents in the two-dimensional
setting, $v_{0}\in\mathbb{R}^{2N}$.

We recall that, for the control $u$ in Eq. (2.8), the size of the transfer
function $\hat{G}$ related to the state-space system (3.1) in terms of the
$\mathcal{H}_{\infty}$ signal norm is

$\|\hat{G}\|_{\mathcal{H}_{\infty}}\leq\gamma.$

From Theorem 3.2 we know that the $\mathcal{H}_{\infty}$ control $u$ is robust
with a constant $\gamma>\frac{1}{c_{N}}$, $c_{N}>0$. We compute the value
$c_{N}$ for the two cases in Figure 4.4, and we have $c_{N}=14.29$ for a
penalization factor $\nu=0.01$, while $c_{N}=4.55$ for $\nu=0.1$. As expected,
we observe smaller regions for a smaller value of $\nu$, interpreted as the
control cost.

Figure 4.4: Test 2. Two-dimensional case. Comparison between the two controls
in (2.8) and (2.10) applied to the uncertain model, in terms of dynamics mean
and variance for different values of $\nu,\rho_{1},\rho_{2}$.

####  4.2.3 Test 3: mean-field consensus dynamics

In the mean-field limit, the Monte Carlo (MC) method is employed for the
approximation of the distribution function $f(t,v,\theta)$ in the phase space
whereas the random space at the particle level is approximated through the SG
technique.

Considering this MC-SG scheme, we work on an agent system using Monte Carlo
sampling with $N_{s}=10^{4}$ agents, then we consider the SG scheme at the
microscopic level. The probability density $f(t,v,\theta)$ is then
reconstructed as the histogram of $v(t,\theta)$. The reconstruction step of
the mean density has been done with $50$ bins. The mean and the variance of
the statistical quantity are computed as follows

$\begin{split}\mathbb{E}_{\theta}[f(t,v,\theta)]&=\int\int
f(t,v,\theta)d\rho_{1}(\theta_{1})d\rho_{2}(\theta_{2})=\sum_{l,h=1}^{L}f(t,v,\theta^{lh})\rho_{1}(\theta_{1}^{l})\rho_{2}(\theta_{2}^{h})\omega_{1}^{l}\omega_{2}^{h},\\\
\mathbb{V}_{\theta}[f(t,v,\theta)]&=\int\int
f(t,v,\theta)^{2}d\rho_{1}(\theta_{1})d\rho_{2}(\theta_{2})-\left(\mathbb{E}_{\theta}[f(t,v,\theta)]\right)^{2}\\\
&=\sum_{l,h=1}^{L}f(t,v,\theta^{lh})^{2}\rho_{1}(\theta_{1}^{l})\rho_{2}(\theta_{2}^{h})\omega_{1}^{l}\omega_{2}^{h}-\left(\mathbb{E}_{\theta}[f(t,v,\theta)]\right)^{2},\end{split}$
(4.9)

where for $l,h=1,\dots,L$ and $\hat{v}_{k,j}\in\mathbb{R}^{N_{s}}$,
$f(t,v,\theta^{lh})$ is reconstructed as the histogram of the data
$v(\theta^{lh})=\sum_{k,j=0}^{M}\hat{v}_{k,j}\Phi_{k}(\theta_{1}^{l})\Psi_{j}(\theta_{2}^{h}).$

(a) $\mathbb{E}_{\theta}[f(t,v,\theta)]$

(b) $\sqrt{\mathbb{V}_{\theta}[f(t,v,\theta)]}$

Figure 4.5: Test 3. Mean-field one-dimensional case. (a) Mean and (b) standard
deviation over time for the $\mathcal{H}_{\infty}$ control in (2.8), with
parameters
$\bar{p}=1,\nu=0.01,\rho_{1}\sim\mathcal{N}(0,5),\rho_{2}\sim\mathcal{U}(-5,5)$.

(a) $\mathbb{E}_{\theta}[f(t,v,\theta)]$

(b) $\sqrt{\mathbb{V}_{\theta}[f(t,v,\theta)]}$

Figure 4.6: Test 3. Mean-field one-dimensional case. (a) mean and (b) standard
deviation over time for the averaged control in (2.10), with parameters
$\bar{p}=1,\nu=0.01,\rho_{1}\sim\mathcal{N}(0,5),\rho_{2}\sim\mathcal{U}(-5,5)$.

We consider the same parameters as in Test 1 for the one dimensional
microscopic case, where $\theta_{1},\theta_{2}$ are uncertainties respectively
with Gaussian distribution $\mathcal{N}(0,5)$ and uniform distribution
$\mathcal{U}(-5,5)$. Hence we approximate the integrals in (4.9) using,
respectively for $\theta_{1}$ and $\theta_{2}$, a Gauss-Hermite and Legendre-
Gauss quadrature rules with $L=40$ quadrature points. Figures 4.5 and 4.6 show
a similar behavior with respect to the microscopic case in left plot of Figure
4.2, in particular we observe that less dispersion of the density for control
of type (2.8).

##  5 Conclusions

The introduction of uncertainties in multiagent systems is of paramount
importance for the description of realistic phenomena. Here we focused on the
mathematical modelling and control of collective dynamics with random inputs
and we investigated the robustness of controls proposing estimates based on
$\mathcal{H}_{\infty}$ theory in the linear setting. Reformulating the control
problem as a robust $\mathcal{H}_{\infty}$ control problem, we derived
sufficient conditions in terms of linear matrix inequalities (LMIs) to ensure
the control performance, independently of the type of random inputs. Moreover,
a robustness analysis is provided also in a mean-field framework, showing
consistent results with the microscopic scale. Different numerical tests were
proposed to compare the $\mathcal{H}_{\infty}$ control with control
synthesized minimizing the expectation of a function with respect to the
random inputs. The numerical methods here developed make use of the stochastic
Galerkin (SG) expansion for the microscopic dynamics while in the mean-field
case we combine an SG expansion in the random space with a Monte Carlo method
in the physical variables. The numerical experiments show that both controls
are capable to drive the average particle trajectories towards a consensus
state considering multiple sources of randomness and in different dimensions.
We further observe that, in the $\mathcal{H}_{\infty}$ setting, the variance
is stabilized over time, this is not surprising since the
$\mathcal{H}_{\infty}$ control accounts for the random state in a feedback
form, whereas in the noiseless control setting the uncertainty is averaged
out. Nonetheless, these results confirm the goodness of the estimates for the
control robustness for the uncertain dynamics. Further analysis is needed to
extend these results in the $\mathcal{H}_{\infty}$ setting to non-linear
dynamics with uncertainities. This can be studied for example introducing the
so-called Hamilton-Jacobi-Isaacs equation, whose solution can be extremely
challenging due to the high-dimensionality of multi-agent systems.

##  Appendix A Averaged control

In this section we consider a control by minimizing the expectation of the
cost functional (2.2) subject to the noisy model (2.1). Hence, we consider the
expected value of the quadratic cost

$\bar{u}^{*}(\cdot)=\arg\min_{u(\cdot)}\mathbb{E}\left[\int_{0}^{T}\frac{1}{2}(v^{\top}Qv+\nu
u^{\top}Ru)\ dt\right],$

where we introduce the matrices $Q=R=\frac{1}{N}\textrm{Id}_{N}$. We claim
that in this case an optimal feedback control is obtained as follows

$\displaystyle\bar{u}(t)=-\frac{N}{\nu}\Bigl{(}{K}\
\mathbb{E}_{\theta}\left[v\right]+{S}\
\mathbb{E}_{\theta}\left[\theta\right]\Bigr{)},$ (A.1)

where ${K}\in\mathbb{R}^{N\times N}$ and ${S}\in\mathbb{R}^{N\times Z}$
fulfill the Riccati matrix-equations

$\begin{cases}-\dot{{K}}={K}{A}+{A}^{\top}{K}-\frac{N}{\nu}K^{2}+\frac{1}{N}I_{N},\quad{K}(T)=0_{N},\\\
\ \\\ -\dot{{S}}=KB+A^{\top}S-\frac{N}{\nu}{K}{S},\quad{S}(T)=0_{N\times
M}.\end{cases}$ (A.2)

######  Theorem A.1.

Assume matrices ${K}$ and $S$ have the following structures

$\displaystyle({K})_{ij}=\begin{cases}{k}_{d},\qquad i=j,\\\ {{k}_{o}},\qquad
i\neq j,\end{cases}\quad({S})_{ij}=s\cdot\mathds{1}_{N\times Z}.$

Matrices $K$ and $S$ are defined by $2$ and $1$ elements respectively.

Then the $i-th$ component of the control $u$ is given by

$u_{i}=-\frac{1}{\nu}\left({k}_{d}\mathbb{E}\left[v_{i}\right]+\frac{{k}_{o}}{N}\sum_{j\neq
i}^{N}\mathbb{E}\left[v_{j}\right]+{s}\sum_{j=1}^{Z}\mathbb{E}\left[\theta_{j}\right]\right),$
(A.3)

where ${k}_{d}$ and ${k}_{o}$, after a scaling, turn out to be the same as in
equations (2.7), while ${s}$ satisfies

$-\dot{{s}}=k_{d}+\alpha(N)k_{o}-\frac{s}{\nu}\left(k_{d}+\alpha(N)k_{o}\right),\qquad
s(T)=0,$ (A.4)

with $\alpha(N)=\frac{N-1}{N}$.

###### Proof.

From Proposition 2.1 of [3], we can prove that ${k}_{d}$, ${k}_{o}$ satisfy
equations (2.7). Given the structure of the matrices $S$, $K$ and $B$, and
solving the second equation in (A.2) componentwise leads to the following
identities:

$\displaystyle(KB)_{ij}$ $\displaystyle=k_{d}+(N-1)k_{o},$
$\displaystyle(A^{\top}S)_{ij}$
$\displaystyle=s\left(a_{d}+(N-1)a_{o}\right),$ $\displaystyle(KS)_{ij}$
$\displaystyle=s\left(k_{d}+(N-1)k_{o}\right).$

We can further simplify the Riccati-matrix system (A.2) using the dependency
of coefficients $a_{d},a_{o}$. And since we are interested in the dynamics of
large number of agents, we introduce the following scalings

$s\leftarrow Ns,\quad{k}_{d}\leftarrow Nk_{d},\quad{k}_{o}\leftarrow
N^{2}k_{o},\quad\alpha(N)=\frac{N-1}{N}.$ (A.5)

For the sake of simplicity, we keep the same notation also for the scaled
variables $s,k_{d},k_{o}$. Under this scaling, the second equation in (A.2)
reads

$-\dot{s}=k_{d}+\alpha(N)k_{o}-\frac{s}{\nu}\left(k_{d}+\alpha(N)k_{o}\right),\qquad
s(T)=0\,,$

and the Riccati feedback law (A.1) is given by Eq. (A.3) The coefficients
$s(t),k_{d}(t),k_{o}(t)$, have to be determined integrating backwards in time.
∎

######  Remark A.1.

In the infinite horizon case and without discount factor, (A.4) reduces to

$\displaystyle
0=k_{d}+\alpha(N)k_{o}-\frac{s}{\nu}\left(k_{d}+\alpha(N)k_{o}\right),$

hence s = $\nu$.

##  Appendix B $\mathcal{H}_{\infty}$ control setting

Define a state-space system $G:L^{2}\rightarrow L^{2}$ by $y=G\theta$, with
$\theta$ a random input as defined in (4.1), if

$\displaystyle\frac{d}{dt}{v}(t)$
$\displaystyle=\widetilde{A}v(t)+\widetilde{B}\theta,$ $\displaystyle y(t)$
$\displaystyle=Cv(t)+D\theta,$

where $v(t)$ is the system state, and $y(t)$ is the observed output. It is
proved (see e.g. [21, 9]) that for any stable state-space system, $G$, there
exists a frequency transfer function $\hat{G}\in R\mathcal{H}_{\infty}$ such
that

$\hat{G}(s)=\left[\begin{array}[]{@{}c|c@{}}\widetilde{A}&\widetilde{B}\\\
\hline\cr C&D\end{array}\right]=D+C(sI_{N}-\widetilde{A})^{-1}\widetilde{B}.$
(B.1)

where $s$ is a complex number and $R\mathcal{H}_{\infty}$ is the set of proper
rational functions with no poles in the closed right half-plane, in particular
$R\mathcal{H}_{\infty}=R\cap\mathcal{H}_{\infty}$, where $R$ is the space of
rational functions and $\mathcal{H}_{\infty}$ is a signal space of “transfer
functions” for linear time-invariant systems, we refer to [9, 20, 25] for
further theoretical details.

State space $\widetilde{A},\widetilde{B},C,D$ or the transfer function is a
representation of a system and these formats uses matrices or complex-valued
functions (a signal) to parameterize the representation. The signal norm
$\|\cdot\|_{\mathcal{H}_{\infty}}$ measures the size of the transfer function
in a certain sense and the $\mathcal{H}_{\infty}$-optimal control problem
consists of finding a stabilizing controller $u=Ky$ which minimizes the cost
function

$\|\hat{G}\|_{\mathcal{H}_{\infty}}=\text{ess}\sup_{\omega\in\mathbb{R}}\bar{\sigma}(\hat{G}(i\omega)).$

The direct minimization of the cost $\|\hat{G}\|_{\mathcal{H}_{\infty}}$ turns
out to be a very hard problem, and it is therefore not feasible to tackle it
directly. Instead, it is much easier to construct conditions that state
whether there exists a stabilizing controller which achieves the norm bound

$\|\hat{G}\|_{\mathcal{H}_{\infty}}\leq\gamma,$

for a given $\gamma>0$.

The history of LMIs in the analysis of dynamical systems begins in about 1890,
when Lyapunov published his seminal work introducing what we now call Lyapunov
theory. One of the major next major contributions that we use in this work
came in the early 1960’s, when Yakubovich, Popov, Kalman, and other
researchers succeeded in reducing the solution of the LMIs to what we now call
the positive-real (PR) lemma, that shows how LMIs can be used to constrain the
eigenvalues of a system [54, 55].

######  Lemma B.1.

Given the frequency transfer function $\hat{G}$, the following are equivalent:

  * •

$\|\hat{G}\|_{\mathcal{H}_{\infty}}\leq\gamma$.

  * •

$\exists$ a positive definite square matrix of order $N$, $X>0$ s.t. (B.2)
holds.

$\begin{bmatrix}\widetilde{A}^{\top}X+X\widetilde{A}&X\widetilde{B}\\\
\widetilde{B}^{\top}X&-\gamma
I_{Z}\end{bmatrix}+\frac{1}{\gamma}\begin{bmatrix}C^{\top}\\\
D^{\top}\end{bmatrix}\begin{bmatrix}C&D\end{bmatrix}<0.$ (B.2)

In [52] the following Lemma with equivalent characterization through a Riccati
equation as been established:

######  Lemma B.2.

The following are equivalent:

  * •

$\exists X>0$ s.t. (B.2) holds.

  * •

$\exists X>0$ s.t. (B.3) holds.

$\widetilde{A}^{\top}X+X\widetilde{A}-(X\widetilde{B}+C^{\top}D)(-\gamma
I_{Z}+\frac{1}{\gamma}D^{\top}D)^{-1}(\widetilde{B}^{\top}X+D^{\top}C)+\frac{1}{\gamma}C^{\top}C=0.$
(B.3)

###### Proof.

The structure of of Eq. (B.2) is

$\begin{bmatrix}\hat{A}&\hat{B}\\\
\hat{B}^{\top}&\hat{D}\end{bmatrix}\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix},$

this system can be solved using the Schur-complement theory. Provided that
$\widehat{D}^{-1}$ exists, we have

$\begin{cases}\hat{A}x_{1}+\hat{B}x_{2}=0,\\\
\hat{B}^{\top}x_{1}+\hat{D}x_{2}=0,\end{cases}$

$\rightarrow\quad
x_{2}=-\hat{D}^{-1}\hat{B}^{\top}x_{1}\quad\rightarrow\quad\left(\hat{A}-\hat{B}\hat{D}^{-1}\hat{B}^{\top}\right)x_{1}=0.$

Hence, provided that exists $X$ such that (B.3) has a solution, then for all
$\xi$

$\xi^{\top}\left(\widetilde{A}^{\top}X+X\widetilde{A}-(X\widetilde{B}+C^{\top}D)(-\gamma
I_{Z}+\frac{1}{\gamma}D^{\top}D)^{-1}(\widetilde{B}^{\top}X+D^{\top}C)+\frac{1}{\gamma}C^{\top}C\right)\xi=0.$

Further, (B.3) is the Schur-complement of

$M_{sc}:=\begin{bmatrix}\widetilde{A}^{\top}X+X\widetilde{A}+\frac{1}{\gamma}C^{\top}C&\widetilde{B}^{\top}X+D^{\top}C\\\
X\widetilde{B}+C^{\top}D&\frac{1}{\gamma}D^{\top}D-\gamma I_{Z}\end{bmatrix}.$

Hence for $\eta=-\hat{D}\hat{B}^{\top}\xi$ and $\forall\xi$, we have that

$\begin{bmatrix}\xi&\eta\end{bmatrix}M_{sc}\begin{bmatrix}\xi\\\
\eta\end{bmatrix}=0.$

For a general vector $\left[\xi\quad\varrho\right]^{\top}$, we compute

$\begin{bmatrix}\xi\\\ \varrho\end{bmatrix}=\begin{bmatrix}\xi\\\
\eta\end{bmatrix}+\begin{bmatrix}0\\\ \varrho-\eta\end{bmatrix},$

and then

$\begin{bmatrix}\xi&\varrho\end{bmatrix}M_{sc}\begin{bmatrix}\xi\\\
\varrho\end{bmatrix}=0+\begin{bmatrix}0&\varrho-\eta\end{bmatrix}M_{sc}\begin{bmatrix}0\\\
\varrho-\eta\end{bmatrix}=\begin{bmatrix}\varrho-\eta\end{bmatrix}^{\top}\begin{bmatrix}\frac{1}{\gamma}D^{\top}D-\gamma
I_{Z}\end{bmatrix}\begin{bmatrix}\varrho-\eta\end{bmatrix}<0$

for $\gamma$ sufficiently large. ∎

## References

  * [1] G. Albi, N. Bellomo, L. Fermo, S.-Y. Ha, J. Kim, L. Pareschi, D. Poyato, and J. Soler.  Vehicular traffic, crowds, and swarms: from kinetic theory and multiscale methods to applications and research perspectives.  Math. Models Methods Appl. Sci., 29(10):1901–2005, 2019.
  * [2] G. Albi, S. Bicego, and D. Kalise.  Gradient-augmented supervised learning of optimal feedback laws using state-dependent Riccati equations.  IEEE Control Systems Letters, 6:836–841, 2021.
  * [3] G. Albi, M. Herty, D. Kalise, and C. Segala.  Moment-driven predictive control of mean-field collective dynamics.  SIAM Journal on Control and Optimization, 60(2):814–841, 2022.
  * [4] G. Albi, L. Pareschi, and M. Zanella.  Uncertainty quantification in control problems for flocking models.  Mathematical Problems in Engineering, 2015, 04 2015.
  * [5] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, et al.  Empirical investigation of starling flocks: a benchmark study in collective animal behaviour.  Animal behaviour, 76(1):201–215, 2008.
  * [6] T. Başar and P. Bernhard.  H-infinity optimal control and related minimax design problems: a dynamic game approach.  Springer Science & Business Media, 2008.
  * [7] N. Bellomo and J. Soler.  On the mathematical theory of the dynamics of swarms viewed as complex systems.  Math. Models Methods Appl. Sci., 22(suppl. 1):1140006, 29, 2012\.
  * [8] M. Bongini, M. Fornasier, M. Hansen, and M. Maggioni.  Inferring interaction rules from observations of evolutive systems i: The variational approach.  Mathematical Models and Methods in Applied Sciences, 27(05):909–951, 2017.
  * [9] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan.  Linear matrix inequalities in system and control theory.  SIAM, 1994.
  * [10] J. A. Carrillo, Y.-P. Choi, and M. Hauray.  The derivation of swarming models: mean-field limit and Wasserstein distances.  In Collective dynamics from bacteria to crowds, pages 1–46. Springer, 2014.
  * [11] J. A. Carrillo, M. R. D’Orsogna, and V. Panferov.  Double milling in self-propelled swarms from kinetic theory.  Kinet. Relat. Models, 2(2):363–378, 2009.
  * [12] J. A. Carrillo, M. Fornasier, G. Toscani, and F. Vecil.  Particle, kinetic, and hydrodynamic models of swarming.  In Mathematical modeling of collective behavior in socio-economic and life sciences, pages 297–336. Springer, 2010.
  * [13] J. A. Carrillo, L. Pareschi, and M. Zanella.  Particle based gPC methods for mean-field models of swarming with uncertainty.  Communications in Computational Physics, 25(2):508–531, 2018.
  * [14] Y.-P. Choi, D. Kalise, J. Peszek, and A. A. Peters.  A collisionless singular Cucker-Smale model with decentralized formation control.  SIAM J. Appl. Dyn. Syst., 18(4):1954–1981, 2019.
  * [15] S. Cordier, L. Pareschi, and G. Toscani.  On a kinetic model for a simple market economy.  J. Stat. Phys., 120(1-2):253–277, 2005.
  * [16] E. Cristiani, B. Piccoli, and A. Tosin.  Multiscale modeling of pedestrian dynamics, volume 12 of MS&A. Model. Simul. Appl. Springer, Cham, 2014.
  * [17] F. Cucker and S. Smale.  Emergent behavior in flocks.  IEEE Trans. Automat. Control, 52(5):852–862, 2007.
  * [18] P. Degond, S. Göttlich, M. Herty, and A. Klar.  A network model for supply chains with multiple policies.  Multiscale Model. Simul., 6(3):820–837, 2007.
  * [19] G. Dimarco, L. Pareschi, and M. Zanella.  Uncertainty quantification for kinetic models in socio–economic and life sciences.  In Uncertainty quantification for hyperbolic and kinetic equations, pages 151–191. Springer, 2017.
  * [20] G.-R. Duan and H.-H. Yu.  LMIs in control systems: analysis, design and applications.  CRC press, 2013.
  * [21] G. E. Dullerud and F. Paganini.  A course in robust control theory: a convex approach, volume 36.  Springer Science & Business Media, 2013.
  * [22] J. R. Dyer, A. Johansson, D. Helbing, I. D. Couzin, and J. Krause.  Leadership, consensus decision making and collective behaviour in humans.  Philos. Trans. Roy. Soc. B, 364(1518):781–789, 2009.
  * [23] M. R. D’Orsogna, Y.-L. Chuang, A. L. Bertozzi, and L. S. Chayes.  Self-propelled particles with soft-core interactions: patterns, stability, and collapse.  Phys. Rev. Lett., 96(10):104302, 2006.
  * [24] G. Estrada-Rodriguez and H. Gimperlein.  Interacting particles with Lévy strategies: limits of transport equations for swarm robotic systems.  SIAM J. Appl. Math., 80(1):476–498, 2020.
  * [25] G. F. Franklin, J. D. Powell, A. Emami-Naeini, and J. D. Powell.  Feedback control of dynamic systems, volume 4.  Prentice hall Upper Saddle River, 2002.
  * [26] G. Freudenthaler and T. Meurer.  PDE-based multi-agent formation control using flatness and backstepping: analysis, design and robot experiments.  Automatica, 115:108897, 13, 2020.
  * [27] J. Garnier, G. Papanicolaou, and T.-W. Yang.  Consensus convergence with stochastic effects.  Vietnam J. Math., 45(1-2):51–75, 2017.
  * [28] B. D. Goddard, B. Gooding, H. Short, and G. Pavliotis.  Noisy bounded confidence models for opinion dynamics: the effect of boundary conditions on phase transitions.  IMA Journal of Applied Mathematics, 87(1):80–110, 2022.
  * [29] J. Gómez-Serrano, C. Graham, and J.-Y. Le Boudec.  The bounded confidence model of opinion dynamics.  Math. Models Methods Appl. Sci., 22(2):1150007, 46, 2012.
  * [30] S.-Y. Ha and E. Tadmor.  From particle to kinetic and hydrodynamic descriptions of flocking.  Kinet. Relat. Models, 1(3):415–435, 2008.
  * [31] Y. Han, A. Hegyi, Y. Yuan, S. Hoogendoorn, M. Papageorgiou, and C. Roncoli.  Resolving freeway jam waves by discrete first-order model-based predictive control of variable speed limits.  Transportation Research Part C: Emerging Technologies, 77:405–420, 2017.
  * [32] M. Herty and L. Pareschi.  Fokker-Planck asymptotics for traffic flow models.  Kinet. Relat. Models, 3(1):165–179, 2010.
  * [33] M. Herty, L. Pareschi, and S. Steffensen.  Mean-field control and Riccati equations.  Netw. Heterog. Media, 10(3):699–715, 2015.
  * [34] M. Herty and C. Ringhofer.  Feedback controls for continuous priority models in supply chain management.  Comput. Methods Appl. Math., 11(2):206–213, 2011.
  * [35] J. Hu and S. Jin.  Uncertainty quantification for kinetic equations.  In Uncertainty quantification for hyperbolic and kinetic equations, pages 193–229. Springer, 2017.
  * [36] J. Hu, S. Jin, and D. Xiu.  A stochastic galerkin method for Hamilton–Jacobi equations with uncertainty.  SIAM Journal on Scientific Computing, 37(5):A2246–A2269, 2015.
  * [37] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin.  Inferring the structure and dynamics of interactions in schooling fish.  Proceedings of the National Academy of Sciences, 108(46):18720–18725, 2011.
  * [38] I. Khalil, J. Doyle, and K. Glover.  Robust and optimal control.  prentice hall, new jersey, 1996.
  * [39] O. Le Maître and O. M. Knio.  Spectral methods for uncertainty quantification: with applications to computational fluid dynamics.  Springer Science & Business Media, 2010.
  * [40] P. Lin and Y. Jia.  Robust H-infinity consensus analysis of a class of second-order multi-agent systems with uncertainty.  IET control theory & applications, 4(3):487–498, 2010.
  * [41] J. Liu, Y. Zhang, H. Liu, Y. Yu, and C. Sun.  Robust event-triggered control of second-order disturbed leader-follower mass: A nonsingular finite-time consensus approach.  International Journal of Robust and Nonlinear Control, 29(13):4298–4314, 2019.
  * [42] Y. Liu and Y. Jia.  Robust H-infinity consensus control of uncertain multi-agent systems with time delays.  International Journal of Control, Automation and Systems, 9, 12 2011\.
  * [43] Y. Luo and W. Zhu.  Event-triggered h-infinity finite-time consensus control for nonlinear second-order multi-agent systems with disturbances.  Advances in Difference Equations, 2021(1):1–19, 2021.
  * [44] L. P. Mo, H. Y. Zhang, and H. Y. Hu.  Finite-time H-infinity consensus of multi-agent systems with a leader.  In Applied Mechanics and Materials, volume 241, pages 1608–1613. Trans Tech Publ, 2013.
  * [45] S. Motsch and E. Tadmor.  Heterophilious dynamics enhances consensus.  SIAM review, 56(4):577–621, 2014.
  * [46] K.-K. Oh, M.-C. Park, and H.-S. Ahn.  A survey of multi-agent formation control.  Automatica, 53:424–440, 2015.
  * [47] M. M. Peet.  Lecture notes in LMI methods in optimal and robust control, 2020.
  * [48] A. A. Peters, R. H. Middleton, and O. Mason.  Leader tracking in homogeneous vehicle platoons with broadcast delays.  Automatica, 50(1):64–74, 2014.
  * [49] R. E. Stern, S. Cui, M. L. Delle Monache, R. Bhadani, M. Bunting, M. Churchill, N. Hamilton, H. Pohlmann, F. Wu, B. Piccoli, et al.  Dissipation of stop-and-go waves via control of autonomous vehicles: Field experiments.  Transp. Research Part C: Emerging Techn., 89:205–221, 2018.
  * [50] G. Toscani.  Kinetic models of opinion formation.  Commun. Math. Sci., 4(3):481–496, 2006.
  * [51] A. Tosin and M. Zanella.  Kinetic-controlled hydrodynamics for traffic models with driver-assist vehicles.  Multiscale Model. Simul., 17(2):716–749, 2019.
  * [52] J. Willems.  Least squares stationary optimal control and the algebraic Riccati equation.  IEEE Transactions on automatic control, 16(6):621–634, 1971.
  * [53] D. Xiu.  Numerical methods for stochastic computations.  Princeton university press, 2010.
  * [54] V. A. Yakubovich.  Solution of certain matrix inequalities encountered in non-linear control theory.  In Doklady Akademii Nauk, volume 156, pages 278–281. Russian Academy of Sciences, 1964.
  * [55] V. A. Yakubovich.  The method of matrix inequalities in the stability theory of nonlinear control systems, i, ii, iii.  Automation and remote control, 25(4):905–917, 1967.

E-mail address<EMAIL_ADDRESS>
E-mail address<EMAIL_ADDRESS>
E-mail address<EMAIL_ADDRESS>


# Fractional Hadamard transform with continuous variables in the context of
quantum optics††thanks: Project supported by the National Natural Science
Foundation of China (Grant Nos 10775097 and 10874174) and the Research
Foundation of the Education Department of Jiangxi Province of China.

Li-yun Hu, Xue-xiang Xu, and Shan-jun Ma
College of Physics and Communication Electronics, Jiangxi Normal University,
Nanchang 330022, China Corresponding author. E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We introduce the quantum fractional Hadamard transform with continuous
variables. It is found that the corresponding quantum fractional Hadamard
operator can be decomposed into a single-mode fractional operator and two
single-mode squeezing operators. This is extended to the entangled case by
using the bipartite entangled state representation. The new transformation
presents more flexibility to represent signals in the fractional Hadamard
domain with extra freedom provided by an angle and two-squeezing parameters.

Keywords: fractional Hadamard transform, fractional Hadamard operator, the
additivity of operator

PACC: 0367, 4250

##  1 Introduction

Fractional Fourier transform (FrFT) is a generalization of the ordinary
Fourier transform, which has been used in signal processing and image
manipulations [1, 2]. The concept of the FrFT was originally described by
Condon [3] and was later introduced for signal processing by Namias [4] as a
Fourier transform of fractional order. The 1-dimension FrFT of $\alpha$-order
is defined in Refs.[5, 6] as

$g\left(x^{\prime}\right)=\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\sin\alpha}}\int
e^{-\frac{\mathtt{i}\left(x^{\prime
2}+x^{2}\right)}{2\tan\alpha}+\frac{\mathtt{i}xx^{\prime}}{\sin\alpha}}f\left(x\right)\mathtt{d}x.$
(1)

The usual Fourier transform is a special case with order $\alpha=\pi/2$. On
the other hand, many orthogonal transform have been successfully used in
signal processing, such as discrete cosine transform [7], discrete Hartley
transform [8] and Hadamard transform.

Hadamard transform is not only an important tool in classical signal
processing, but also is of great importance for quantum computation
applications [9]. This transform, used to go from the position basis
$\left|x\right\rangle$ to the momentum basis, is defined as [10, 11]

$\mathcal{F}\left|x\right\rangle=\frac{1}{\sqrt{\pi}\sigma}\int_{-\infty}^{\infty}e^{2\mathtt{i}xy/\sigma^{2}}\left|y\right\rangle\mathtt{d}y,$
(2)

where $\sigma$ is the scale length (also makes the expression in the
exponential dimensionless), $\left|x\right\rangle$ and $\left|y\right\rangle$
are the eigenvector of coordinate operator $X.$ In Ref.[12], the explicit form
of $\mathcal{F}$ has been derived by using the technique of integration within
an ordered product (IWOP) of operators [13, 14, 15], and it is found that it
can be decomposed into a single-mode squeezing operator and a position-
momentum mutual transform operator, i.e., $\mathcal{F=}S_{1}^{-1}(-1)^{i\pi
a^{{\dagger}}a/2}$. In addition, the two-mode Hadamard transform with
continuous variables is also introduced by using the bipartite entangled state
representation, whose Hadamard operator involves a two-mode squeezing operator
and a mutual transform operator.

In this paper, we shall introduce the continuous fractional Hadamard transform
(CFrHT), which is a generalization of the usual Hadamard transform in Eq.(2).
The development of the CFrHT is based upon the same spirit of continuous
fractional Fourier transform (CFrFT). Then the CFrHT operator (CFrHTO) is
derived by using the IWOP technique, and its properties are analyzed. It is
found that the CFrHTO can be decomposed into a single-mode fractional operator
$e^{\mathtt{i}\alpha a^{\dagger}a}$ and two single-mode squeezing operators.
On the other hand, since the publication of the paper of Einstein, Podolsky
and Rosen (EPR) in 1935 [16], arguing the incompleteness of quantum mechanics,
the conception of entanglement has become more and more fascinating and
important as it plays a central role in quantum imformation and quantum
computation, we also shall introduce the two-mode CFrHO in bipartite entangled
state representation, which turns out to involve the two fractional operators
and two two-mode squeezing operators.

Our work is arranged as follows. In section 2, for the single-mode case, the
normally ordered fractional Hadamard operator is derived by using the IWOP
technique. The properties of fractional Hadamard operator is discussed in
section 3, such as the unitarity, the decomposition of the CFrHO and its
transform relation. Then the single-mode case is extended to two-mode case in
section 4 and some similar discussions to singled-mode case are presented.
Section 5 is devoted to exploring the measurements for the output states from
the CFrHT. Conclusions are involved in the last section.

##  2 Normally Ordered Fractional Hadamard Operator

In this section, we first introduce the continuous fractional Hadamard
transform (CFrHT), i.e.,

$\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|x\right\rangle=\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\mu\nu\sin\alpha}}\int_{-\infty}^{\infty}\exp\left\\{-\frac{\mathtt{i}\left(x^{2}/\mu^{2}+y^{2}/\nu^{2}\right)}{2\tan\alpha}+\frac{\mathtt{i}xy}{\mu\nu\sin\alpha}\right\\}\left|y\right\rangle\mathtt{d}y,$
(3)

where $\mu,\nu$ are the scale length (also make the expression in the
exponential dimensionless), $\alpha$ is an angle, and
$\mathcal{H}^{\alpha}\left(\mu,\nu\right)$ is called the CFrHT operator
(CFrHO). In particular, when $\alpha=\pi/2$ and $\mu=\nu=\sigma/\sqrt{2},$
Eq.(3) just reduces to Eq.(2).

In order to find the explicit expression of the CFrHO, multiplying Eq.(3) by
the bra $\int\mathtt{d}x\left\langle x\right|$ from the rights in two-side,
where $\left|y\right\rangle$ and $\left|x\right\rangle$ are coordinate
eigenvectors, $X\left|x\right\rangle=x\left|x\right\rangle$, and

$\left|x\right\rangle=\pi^{-1/4}\exp\left\\{-\frac{x^{2}}{2}+\sqrt{2}xa^{\dagger}-\frac{a^{\dagger
2}}{2}\right\\}\left|0\right\rangle,$ (4)

we can recast the CFrHO $\mathcal{H}_{\alpha}\left(\mu,\nu\right)$ into the
following integral form,

$\mathcal{H}_{\alpha}\left(\mu,\nu\right)=\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\mu\nu\sin\alpha}}\int_{-\infty}^{\infty}\exp\left\\{-\frac{\mathtt{i}\left(x^{2}/\mu^{2}+y^{2}/\nu^{2}\right)}{2\tan\alpha}+\frac{\mathtt{i}xy}{\mu\nu\sin\alpha}\right\\}\left|y\right\rangle\left\langle
x\right|\mathtt{d}x\mathtt{d}y.$ (5)

Then using the vacuum projector’s normal ordering form
$\left|0\right\rangle\left\langle 0\right|=\colon
e^{-a^{{\dagger}}a}\colon$(where the symbol $\colon$ $\colon$ denotes the
normally ordering) and the IWOP technique to directly perform the integration,
we finally obtain

$\displaystyle\mathcal{H}_{\alpha}\left(\mu,\nu\right)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\mu\nu\sin\alpha}}\int_{-\infty}^{\infty}\colon\exp\left\\{-\frac{A}{2\mu^{2}}x^{2}+\sqrt{2}xa+\frac{\mathtt{i}xy}{\mu\nu\sin\alpha}\right.$
(6)
$\displaystyle-\left.\frac{B}{2\nu^{2}}y^{2}+\sqrt{2}ya^{\dagger}-\frac{\left(a^{{\dagger}}+a\right)^{2}}{2}\right\\}\colon\mathtt{d}x\mathtt{d}y$
$\displaystyle=$ $\displaystyle\sqrt{\frac{2\mu\nu
e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{u\sin\alpha}}\exp\left\\{\left(\frac{\nu^{2}A}{u}-\frac{1}{2}\right)a^{{\dagger}2}\right\\}$
$\displaystyle\times\exp\left\\{a^{\dagger}a\ln\frac{\mathtt{i}2\mu\nu}{u\sin\alpha}\right\\}\exp\left\\{\left(\frac{\mu^{2}B}{u}-\frac{1}{2}\right)a^{2}\right\\},$

where we have set $\allowbreak
A=\mathtt{i}\cot\alpha+\mu^{2},B=\mathtt{i}\cot\alpha+\nu^{2},u=\csc^{2}\alpha+AB,$and
used the operator identity in the last step of Eq.(6),

$\exp\left\\{fa^{\dagger}a\right\\}=\colon\exp\left\\{\left(e^{f}-1\right)a^{\dagger}a\right\\}\colon,$
(7)

Eq.(6) is the normally ordered form of the CFrHO. In particular, when
$\alpha=\pi/2$ and $\mu=\nu=\sigma/\sqrt{2},$ leading to
$A=B=\sigma^{2}/2,u=1+\sigma^{4}/4,$ then Eq.(6) becomes

$\displaystyle\mathcal{H}_{\pi/2}\left(\sigma/\sqrt{2},\sigma/\sqrt{2}\right)$
$\displaystyle=$
$\displaystyle\frac{2\sigma}{\sqrt{\sigma^{4}+4}}\exp\left\\{\frac{\sigma^{4}-4}{\sigma^{4}+4}\frac{a^{\dagger
2}}{2}\right\\}$ (8)
$\displaystyle\times\exp\left\\{a^{\dagger}a\ln\frac{4\mathtt{i}\sigma^{2}}{\sigma^{4}+4}\right\\}\exp\left\\{\frac{\sigma^{4}-4}{\sigma^{4}+4}\frac{a^{2}}{2}\right\\},$

which is just the result Eq.(7) in Ref.[12].

##  3 Properties of Fractional Hadamard Operator

From Eq.(5) one can see that the CFrHO is a unitary one, i.e.,
$\mathcal{H}_{\alpha}\left(\mu,\nu\right)\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)=\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)\mathcal{H}_{\alpha}\left(\mu,\nu\right)=1.$
In fact, uisng Eq.(5) and the orthogonality of coordinate state, $\left\langle
x^{\prime}\right.\left|x\right\rangle=\delta\left(x-x^{\prime}\right)$, we
have

$\displaystyle\mathcal{H}_{\alpha}\left(\mu,\nu\right)\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi\mu\nu\sin\alpha}\int_{-\infty}^{\infty}\exp\left\\{\allowbreak\frac{\mathtt{i}\left(y^{\prime
2}-y^{2}\right)}{2\nu^{2}\tan\alpha}+\mathtt{i}x\frac{y-y^{\prime}}{\mu\nu\sin\alpha}\right\\}\left|y\right\rangle\left\langle
y^{\prime}\right|\mathtt{d}x\mathtt{d}y\mathtt{d}y^{\prime}$ (9)
$\displaystyle=$
$\displaystyle\frac{1}{\mu\nu\sin\alpha}\int_{-\infty}^{\infty}\delta\left(\frac{y-y^{\prime}}{\mu\nu\sin\alpha}\right)\exp\left\\{\allowbreak\frac{\mathtt{i}\left(y^{\prime
2}-y^{2}\right)}{2\nu^{2}\tan\alpha}\right\\}\left|y\right\rangle\left\langle
y^{\prime}\right|\mathtt{d}y\mathtt{d}y^{\prime}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left|y\right\rangle\left\langle
y\right|\mathtt{d}y=\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)\mathcal{H}_{\alpha}\left(\mu,\nu\right)=1.$

In order to see clearly its transform relation under the CFrHO, next we
examine its decomposition. Performing the change of variables,
$x/\mu\rightarrow x,$ $y/\nu\rightarrow y$, we can be recast Eq.(5) into the
following form,

$\mathcal{H}_{\alpha}\left(\mu,\nu\right)=\sqrt{\frac{\mu\nu
e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\sin\alpha}}\int_{-\infty}^{\infty}\exp\left\\{-\frac{\mathtt{i}\left(x^{2}+y^{2}\right)}{2\tan\alpha}+\frac{\mathtt{i}xy}{\sin\alpha}\right\\}\left|\nu
y\right\rangle\left\langle\mu x\right|\mathtt{d}x\mathtt{d}y.$ (10)

By noticing that the single-mode squeezing operator $S_{1}$ [17] has its
natural expression in coordinate representation [13], i.e.,

$S_{1}\left(\mu\right)=\frac{1}{\sqrt{\mu}}\int_{-\infty}^{\infty}\mathtt{d}x\left|\frac{x}{\mu}\right\rangle\left\langle
x\right|,$ (11)

which leads to $\left|\nu
y\right\rangle=\frac{1}{\sqrt{\nu}}S_{1}^{-1}\left(\nu\right)\left|y\right\rangle,$
$\left\langle\mu x\right|=\frac{1}{\sqrt{\mu}}\left\langle
x\right|S_{1}\left(\mu\right),$ so Eq.(10) can be decomposed into

$\mathcal{H}_{\alpha}\left(\mu,\nu\right)=S_{1}^{-1}\left(\nu\right)\mathcal{F}_{\alpha}S_{1}\left(\mu\right)=S_{1}^{-1}\left(\nu\right)e^{\mathtt{i}\alpha
a^{\dagger}a}S_{1}\left(\mu\right),$ (12)

where $\mathcal{F}_{\alpha}$ is given by

$\mathcal{F}_{\alpha}\equiv\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\sin\alpha}}\int_{-\infty}^{\infty}e^{-\frac{\mathtt{i}(x^{2}+y^{2})}{2\tan\alpha}+\frac{\mathtt{i}xy}{\sin\alpha}}\left|y\right\rangle\left\langle
x\right|\mathtt{d}x\mathtt{d}y=e^{\mathtt{i}\alpha a^{\dagger}a},$ (13)

this integral result can be obtained by using a similar way to deriving
Eq.(6). Thus we see that the CFrHO can be decomposed as a fractional operator
and two-single-mode squeezing operators.

Using the decomposition of the CFrHO in Eq.(12), and noticing that
$S_{1}\left(\mu\right)XS_{1}^{-1}\left(\mu\right)=\mu X,$
$S_{1}\left(\mu\right)PS_{1}^{-1}\left(\mu\right)=P/\mu,$ and
$e^{\mathtt{i}\alpha a^{\dagger}a}ae^{-\mathtt{i}\alpha
a^{\dagger}a}=ae^{-\mathtt{i}\alpha},$ which leads to

$e^{\mathtt{i}\alpha a^{\dagger}a}Xe^{-\mathtt{i}\alpha
a^{\dagger}a}=\allowbreak X\cos\alpha+P\sin\alpha,\text{ }e^{\mathtt{i}\alpha
a^{\dagger}a}Pe^{-\mathtt{i}\alpha a^{\dagger}a}=\allowbreak
P\cos\alpha-X\sin\alpha,$ (14)

thus we have

$\displaystyle\mathcal{H}_{\alpha}\left(\mu,\nu\right)X\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)$
$\displaystyle=$ $\displaystyle\mu
S_{1}^{-1}\left(\nu\right)\left(X\cos\alpha+P\sin\alpha\right)S_{1}\left(\nu\right)$
(15) $\displaystyle=$ $\displaystyle\frac{\mu}{\nu}X\cos\alpha+\mu\nu
P\sin\alpha,$
$\displaystyle\mathcal{H}_{\alpha}\left(\mu,\nu\right)P\mathcal{H}_{\alpha}^{\dagger}\left(\mu,\nu\right)$
$\displaystyle=$
$\displaystyle\frac{\allowbreak\nu}{\mu}P\cos\alpha-\frac{X}{\mu\nu}\sin\alpha,$
(16)

from which we see that the CFrHO plays the role of combining the coordinate
operator $X$ and momentum operator $P$ in a certain way (15)-(16), i.e.,
including the squeezing and the rotation. In paticular, when
$\alpha=\frac{\pi}{2},$ Eqs.(15)-(16) become

$\mathcal{H}_{\frac{\pi}{2}}\left(\mu,\nu\right)X\mathcal{H}_{\frac{\pi}{2}}^{\dagger}\left(\mu,\nu\right)=\mu\nu
P,\text{
}\mathcal{H}_{\frac{\pi}{2}}\left(\mu,\nu\right)P\mathcal{H}_{\frac{\pi}{2}}^{\dagger}\left(\mu,\nu\right)=-\frac{X}{\mu\nu},$
(17)

i.e., the mutual exchanging of coordinate-momentum operators.

On the other hand, there is a most important feature of the FrFT is that the
FrFT obeys the additivity rule, i.e., two successive FrFT of order $\alpha$
and $\beta$ makes up the FrFT of order $\alpha+\beta$. Then a question
naturally arises: Is the two successive CFrHOs still a CFrHO? To answer this
question, we examine the direct product
$\mathcal{H}_{\alpha}\left(\mu,\nu\right)\otimes\mathcal{H}_{\beta}\left(\mu^{\prime},\nu^{\prime}\right)$.
Using Eq.(12) it is easily seen that when $\mu=$ $\nu^{\prime}$ there is an
additivity of operator as follows

$\displaystyle\mathcal{H}_{\alpha}\left(\mu,\nu\right)\otimes\mathcal{H}_{\beta}\left(\mu^{\prime},\mu\right)$
$\displaystyle=$ $\displaystyle S_{1}^{-1}\left(\nu\right)e^{\mathtt{i}\alpha
a^{\dagger}a}S_{1}\left(\mu\right)S_{1}^{-1}\left(\mu\right)e^{\mathtt{i}\beta
a^{\dagger}a}S_{1}\left(\mu^{\prime}\right)$ (18) $\displaystyle=$
$\displaystyle\mathcal{H}_{\alpha+\beta}\left(\mu^{\prime},\nu\right),$

which can be seen as the additivity property of the CFrHOs. Here it should be
pointed out that the condition of additivitive operator for the CFrHOs is that
the parameter $\mu$ of the prior cascade operator should be equal to the
parameter $\nu^{\prime}$ of the next one, i.e., $\mu=\nu^{\prime}.$ This can
be clearly seen from the viewpoint of classical optics transform.

##  4 Two-mode CFrHT

Next, we shall extend the single-mode CFrHT to two-mode case by using the
entangled state representation [18],

$\left|\eta\right\rangle=\exp\left\\{-\frac{1}{2}\left|\eta\right|^{2}+\eta
a_{1}^{\dagger}-\eta^{\ast}a_{2}^{\dagger}+a_{1}^{\dagger}a_{2}^{\dagger}\right\\}\left|00\right\rangle,$
(19)

where $\left|\eta=\eta_{1}+\mathtt{i}\eta_{2}\right\rangle$ is the common
eigenvector of two-particle’s relative coordinate $X_{1}-X_{2}$ and total
momentum $P_{1}+P_{2},$

$\left(X_{1}-X_{2}\right)\left|\eta\right\rangle=\sqrt{2}\eta_{1}\left|\eta\right\rangle,\left(P_{1}+P_{2}\right)\left|\eta\right\rangle=\sqrt{2}\eta_{2}\left|\eta\right\rangle,$
(20)

and $\left|\eta\right\rangle$ possesses the completeness and the
orthogonality,

$\int_{-\infty}^{\infty}\frac{\mathtt{d}^{2}\eta}{\pi}\left|\eta\right\rangle\left\langle\eta\right|=1,\text{
}\left\langle\eta\right|\left.\eta^{\prime}\right\rangle=\pi\delta\left(\eta-\eta^{\prime}\right)\delta\left(\eta^{\ast}-\eta^{\prime\ast}\right).$
(21)

In a similar way to introducing Eq.(3), we examine the following transform,

$\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)\left|\eta\right\rangle=\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\mu\nu\sin\alpha}\int\frac{\mathtt{d}^{2}\eta^{\prime}}{\pi}e^{-\frac{\mathtt{i}\left(\allowbreak\left|\eta^{\prime}\right|^{2}/\nu^{2}+\left|\eta\right|^{2}/\mu^{2}\allowbreak\right)}{2\tan\alpha}+\frac{\mathtt{i}\left(\eta^{\prime}{}^{\ast}\allowbreak\eta+\eta^{\ast}\allowbreak\eta^{\prime}\right)}{2\mu\nu\sin\alpha}}\left|\eta^{\prime}\right\rangle.$
(22)

Using Eq.(21), one can see that $\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)$
is a unitary operator, i.e.,
$\mathcal{H}_{\alpha}^{C}\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}=\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}\mathcal{H}_{\alpha}^{C}=1$.
Here we should emphasize that, the exponential item in the right hand side of
Eq.(22) can be decomposed into a direct product of two exponential items in
the right hand side of Eq.(3), but $\left|\eta\right\rangle$ is an entangled
state (not the direct product of two single-mode coordinate states, which can
be seen clearly from its Schmidt decomposition [19]), thus this is a
nontrivial extension from single-mode case to two-mode case.

Performing a similar procedure to single-mode case, i.e., noticing that the
two-mode squeezing operator has its natural expression in the entangled state
representation,

$S_{2}\left(\mu\right)=\exp\left[\left(a_{1}^{\dagger}a_{2}^{\dagger}-a_{1}a_{2}\right)\ln\mu\right]=\frac{1}{\mu}\int\frac{d^{2}\eta}{\pi}\left|\frac{\eta}{\mu}\right\rangle\left\langle\eta\right|,$
(23)

which leads to
$\frac{1}{\nu}\left|\frac{\eta^{\prime}}{\nu}\right\rangle=S_{2}\left(\nu\right)\left|\eta^{\prime}\right\rangle,\frac{1}{\mu}\left|\frac{\eta}{\mu}\right\rangle=S_{2}\left(\mu\right)\left|\eta\right\rangle,$
then using the completeness of $\left|\eta\right\rangle$ and
$\left|00\right\rangle\left\langle 00\right|=\colon
e^{-a^{{\dagger}}a-b^{{\dagger}}b}\colon$ and the orthogonality in Eq.(21) we
can further decompose the operator
$\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)$ into the following form,

$\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)=S_{2}^{\dagger}\left(\nu\right)\mathcal{F}_{\alpha}^{C}S_{2}\left(\mu\right)=S_{2}^{\dagger}\left(\nu\right)\exp\left\\{\mathtt{i}\alpha\left(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2}\right)\right\\}S_{2}\left(\mu\right),$
(24)

where the operatror $\mathcal{F}_{\alpha}^{C}$ is given by

$\displaystyle\mathcal{F}_{\alpha}^{C}$ $\displaystyle=$
$\displaystyle\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\sin\alpha}\int\frac{\mathtt{d}^{2}\eta^{\prime}\mathtt{d}^{2}\eta}{\mu^{2}\nu^{2}\pi}e^{-\frac{\mathtt{i}\left(\allowbreak\left|\eta^{\prime}\right|^{2}/\nu^{2}+\left|\eta\right|^{2}/\mu^{2}\allowbreak\right)}{2\tan\alpha}+\frac{\mathtt{i}\left(\eta^{\prime}{}^{\ast}\allowbreak\eta+\eta^{\ast}\allowbreak\eta^{\prime}\right)}{2\mu\nu\sin\alpha}}\left|\frac{\eta^{\prime}}{\nu}\right\rangle\left\langle\frac{\eta}{\mu}\right|$
(25) $\displaystyle=$
$\displaystyle\exp\left\\{\mathtt{i}\alpha\left(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2}\right)\right\\}.$

Thus we see that the two-mode CFrHO can be decomposed into the form in
Eq.(24), i.e., two fractional operators and two two-mode squeezing operators.

This is a convient expression for further deriving the transforms and the
condition of additivitive operator. In fact, using Eqs.(24), (14) and
Eqs.(20), (21) leading to

$\displaystyle
S_{2}\left(\mu\right)\left(X_{1}-X_{2}\right)S_{2}^{\dagger}\left(\mu\right)$
$\displaystyle=$ $\displaystyle\mu\left(X_{1}-X_{2}\right),\text{
}S_{2}\left(\mu\right)\left(P_{1}+P_{2}\right)S_{2}^{\dagger}\left(\mu\right)=\mu\left(P_{1}+P_{2}\right),$
(26) $\displaystyle
S_{2}\left(\mu\right)\left(X_{1}+X_{2}\right)S_{2}^{\dagger}\left(\mu\right)$
$\displaystyle=$ $\displaystyle\frac{1}{\mu}\left(X_{1}+X_{2}\right),\text{
}S_{2}\left(\mu\right)\left(P_{1}-P_{2}\right)S_{2}^{\dagger}\left(\mu\right)=\frac{1}{\mu}\left(P_{1}-P_{2}\right),$
(27)

we have

$\displaystyle\mathcal{H}_{\alpha}^{C}\left(X_{1}-X_{2}\right)\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}$
$\displaystyle=$ $\displaystyle\mu
S_{2}^{\dagger}\left(\nu\right)\mathcal{F}_{\alpha}^{C}\left(X_{1}-X_{2}\right)\left[\mathcal{F}_{\alpha}^{C}\right]^{{\dagger}}S_{2}\left(\nu\right)$
(28) $\displaystyle=$ $\displaystyle\mu
S_{2}^{\dagger}\left(\nu\right)\left(\left(X_{1}-X_{2}\right)\cos\alpha+\left(P_{1}-P_{2}\right)\sin\alpha\right)S_{2}\left(\nu\right)$
$\displaystyle=$
$\displaystyle\frac{\mu}{\nu}\left(X_{1}-X_{2}\right)\cos\alpha+\mu\nu\left(P_{1}-P_{2}\right)\sin\alpha,$
$\displaystyle\mathcal{H}_{\alpha}^{C}\left(X_{1}+X_{2}\right)\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}$
$\displaystyle=$
$\displaystyle\frac{\nu}{\mu}\left(X_{1}+X_{2}\right)\cos\alpha+\frac{1}{\mu\nu}\left(P_{1}+P_{2}\right)\sin\alpha,$
(29)

and

$\displaystyle\mathcal{H}_{\alpha}^{C}\left(P_{1}-P_{2}\right)\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}$
$\displaystyle=$
$\displaystyle\frac{\nu}{\mu}\left(P_{1}-P_{2}\right)\cos\alpha-\frac{1}{\mu\nu}\left(X_{1}-X_{2}\right)\sin\alpha,$
(30)
$\displaystyle\mathcal{H}_{\alpha}^{C}\left(P_{1}+P_{2}\right)\left[\mathcal{H}_{\alpha}^{C}\right]^{{\dagger}}$
$\displaystyle=$
$\displaystyle\frac{\mu}{\nu}\left(P_{1}+P_{2}\right)\cos\alpha-\mu\nu\left(X_{1}+X_{2}\right)\sin\alpha.$
(31)

From Eqs.(28)-(31) it is easy to see that when $\alpha=\pi/2,$ the role of
$\mathcal{H}_{\pi/2}^{C}$ is just exchanging $\left(X_{1}-X_{2}\right)$ and
$\left(P_{1}-P_{2}\right)$, $\left(X_{1}+X_{2}\right)$ and
$\left(P_{1}+P_{2}\right);$ while for $\alpha=\pi,$ $\mathcal{H}_{\pi}^{C}$
can be seen as an identity operator.

In addition, from the decomposition (24) one can see that the direct product
$\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)\otimes\mathcal{H}_{\beta}^{C}\left(\mu^{\prime},\nu^{\prime}\right)$
satisfies the additivity rule when $\mu=$ $\nu^{\prime}$, i.e.,

$\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)\otimes\mathcal{H}_{\beta}^{C}\left(\mu^{\prime},\nu^{\prime}\right)=\mathcal{H}_{\alpha+\beta}^{C}\left(\mu^{\prime},\nu\right).$
(32)

##  5 Measurements for the output states from the CFrHT

The measurement for quantum state plays an important role in quantum
computation and quantum imfromation. When a quantum state
$\left|f\right\rangle$ is transformed by the CFrHO, then what is the
measurement result with continuous orthogonal basis? For single-mode case, the
output state from the CFrHT is
$\left|g\right\rangle_{out}=\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|f\right\rangle$.
The measurement basis is choosen as a coordiante eigenvector, then the
measurement result is given by

$\displaystyle\left\langle x\right.\left|g\right\rangle_{out}$
$\displaystyle=$ $\displaystyle\left\langle
x\right|\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|f\right\rangle$ (33)
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\mathtt{d}x^{\prime}\left\langle
x\right|\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|x^{\prime}\right\rangle\left\langle
x^{\prime}\right.\left|f\right\rangle$ $\displaystyle=$
$\displaystyle\sqrt{\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\pi\mu\nu\sin\alpha}}\int_{-\infty}^{\infty}f\left(x^{\prime}\right)e^{-\frac{\mathtt{i}\left(x^{\prime
2}/\mu^{2}+x^{2}/\nu^{2}\right)}{2\tan\alpha}+\frac{\mathtt{i}x^{\prime}x}{\mu\nu\sin\alpha}}\mathtt{d}x^{\prime},$

which just corresponds to a generalized fractional Fourier transform of wave
function $f\left(x^{\prime}\right)=\left\langle
x^{\prime}\right.\left|f\right\rangle.$

For two-mode case, the measurement result by two-mode entangled state Bell
basis is

$\displaystyle\left\langle\eta^{\prime}\right.\left|g\right\rangle_{out}$
$\displaystyle=$
$\displaystyle\left\langle\eta^{\prime}\right|\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)\left|f\right\rangle$
(34) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{\mathtt{d}^{2}\eta}{\pi}\left\langle\eta^{\prime}\right|\mathcal{H}_{\alpha}^{C}\left(\mu,\nu\right)\left|\eta\right\rangle\left\langle\eta\right.\left|f\right\rangle$
$\displaystyle=$
$\displaystyle\frac{e^{\mathtt{i}\left(\frac{\pi}{2}-\alpha\right)}}{2\mu\nu\sin\alpha}\int_{-\infty}^{\infty}\frac{\mathtt{d}^{2}\eta}{\pi}e^{-\frac{\mathtt{i}\left(\allowbreak\left|\eta^{\prime}\right|^{2}/\nu^{2}+\left|\eta\right|^{2}/\mu^{2}\allowbreak\right)}{2\tan\alpha}+\frac{\mathtt{i}\left(\eta^{\prime}{}^{\ast}\allowbreak\eta+\eta^{\ast}\allowbreak\eta^{\prime}\right)}{2\mu\nu\sin\alpha}}f\left(\eta\right),$

which is just a generalized complex fractional Fourier transform, and the wave
function $f\left(\eta\right)$ is the projection of quantum state
$\left|f\right\rangle$ on $\left\langle\eta\right|$. From Eqs.(33) and (34) we
can clearly see that the generalized FrFT of the wavefunction for any quantum
state $\left|f\right\rangle$ in coordinate/entangled state corresponds to the
wavefunction of Hadamard-transformed
($\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|f\right\rangle$) in
coordinate/entangled state. In other words, the generalized FrFT of
wavefunction is just the matrix element of CFrHO in $\left\langle x\right|$
($\left\langle\eta^{\prime}\right|$) and $\left|f\right\rangle$.

##  6 Conclusion

Based on quantum Hadamard transform, we have introduced the quantum fractional
Hadamard transform with continuous variables. It is found that the
corresponding quantum fractional Hadamard operator can be decomposed into a
single-mode fractional operator $e^{\mathtt{i}\alpha a^{\dagger}a}$ and two
single-mode squeezing operators. The two-mode fractional Hadamard transform is
also introduced by using the bipartite entangled state representation. It is
shown that the corresponding transform operatror involves two single-mode
fractional operators and two two-mode squeezing operators. For any quantum
state vector $\left|f\right\rangle$, the measurement results for the
transformed quantum state (for instance
$\left|g\right\rangle_{out}=\mathcal{H}_{\alpha}\left(\mu,\nu\right)\left|f\right\rangle)$
by continuous coordinate state $\left|x\right\rangle$ (or bipartite entangled
state $\left|\eta\right\rangle)$ just corresponds to a generalized (complex)
fractional Fourier transform. In addition, the new transformation gives us
more flexibility to represent signals in the fractional Hadamard domain with
extra freedom provided by an angle, and two-squeezing parameters. For more
discussions about the optical transforms in the context of quantum optics and
the discrete fractional Hadmard transform, we refer to Refs.[20, 21, 22].

## References

  * [1] L. B. Almeida 1994 IEEE _Trans. Signal Process_. 42 3084
  * [2] Tao R, Deng B, Zhang W Q, et. al. 2008 _IEEE Trans. Signal Process._ 56 158
  * [3] Condon E U 1937 Proc. Nat. Acad. Sci. 23 158
  * [4] Namias V 1980 _J. Inst. Math. Appl._ 25 241
  * [5] McBride A C and Kerr F H 1987 _IMA J. Appl. Math._ 39 159
  * [6] Kerr F H 1988 _J. Math. Anal. Appl._ 136 404
  * [7] Pei S C, Yeh M H 2001 _IEEE Trans. Signal Process._ 49 1198
  * [8] Pei S C, Tseng C C, Yeh M H and Shyu J J 1998 _IEEE Trans. Circuit Ssytem II_ 45 665
  * [9] Nielsen M A and Chuang I L 2000 The Quantum Computation and Quantum Information, Cambridge University Press, Cambridge
  * [10] Braunstein S L 1998 _Phys. Rev. Lett._ 80 4084
  * [11] Parker S, Bose S and Plenio M B 2000 _Phys. Rev. A._ 61 032305
  * [12] Fan H Y and Guo Q 2008 _Commun. Theor. Phys._ 49 859
  * [13] Fan H Y, H. R. Zaidi and J. R. Klauder 1987 Phys. Rev. D 35 1831
  * [14] Fan H Y 2003 J. Opt. B: Quantum. Semiclass. Opt. 5 R147
  * [15] Hu L Y and Fan H Y 2009 _Commun. Theor. Phys._ (Beijing, China) 52 1071; Hu L Y and Fan H Y 2009 Chin. Phys. B 18 4657
  * [16] Einstein A, Podolsky B and Rosen N 1935 _Phys. Rev._ 47 777
  * [17] Mandel L and Wolf E 1995 Optical Coherence and Quantum Optics, Cambridge Press, (London)
  * [18] Fan H Y and Klauder J R 1994 Phys. Rev. A 49 704; Fan H Y and Fan Y 1996 Phys. Rev. A 54 958
  * [19] Preskill J 1998 Quantum Information and Computation, California Institute of Technology; Hu L Y and Lu H L. 2007 Chin. Phys. 16 2200
  * [20] Hu L Y and Fan H Y 2008 _J. Mod. Opt._ 55 1835; Fan H Y and Hu L Y 2009 _Chin. Phys. B_ 18 0611
  * [21] Tao R, Lang J, Wang Y 2009 Opt. Commun. 282 1531
  * [22] Hu L Y and Fan H Y 2009 Int. J. Theor. Phys. DOI 10.1007/s10773-009-0008-z