# Uncertainty Quantification in Control Problems for Flocking Models

Authors: Giacomo Albi; Lorenzo Pareschi; Mattia Zanella
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The optimal control of flocking models with random inputs is investigated from a numerical point of view. The effect of uncertainty in the interaction parameters is studied for a Cucker-Smale type model using a generalized polynomial chaos (gPC) approach. Numerical evidence of threshold effects in the alignment dynamic due to the random parameters is given. The use of a selective model predictive control permits steering of the system towards the desired state even in unstable regimes.


---

## Body

## 1. Introduction

The aggregate motion of a multiagent system is frequently seen in the real world. Common examples are schools of fishes, swarms of bees and herds of sheep, natural phenomena that inspired important applications in many fields such as biology, engineering and economy [1]. As a consequence, the significance of new mathematical models, for understanding and predicting these complex dynamics, is widely recognized. Several heuristic rules of flocking have been introduced as alignment, separation, and cohesion [2, 3]. Nowadays these mathematical problems, and their constrained versions, are deeply studied from both the microscopic viewpoint [4–7] and their kinetic and mean-field approximations [8–13]. We refer to [1] for a recent introduction on the subject.In an applicative framework a fundamental step for the study of such models is represented by the introduction of stochastic parameters reflecting the uncertainty due to wide range of phenomena, such as the weathers influence during an experiment, temperature variations, or even human errors. Therefore quantifying the influence of uncertainties on the main dynamics is of paramount importance to build more realistic models and to give better predictions of their behavior. In the modeling of self-organized system, different ways to include random sources have been studied and analyzed; see, for example, [3, 14–17]. In this paper we focus on the case where the uncertainty acts directly in the parameter characterizing the interaction dynamic between the agents.We present a numerical approach having roots in the numerical techniques for uncertainty quantification (UQ) and model predictive control (MPC). Among the most popular methods for UQ, the generalized polynomial chaos (gPC) has recently received deepest attentions [18]. Jointly with Stochastic Galerkin (SG) this class of numerical methods is usually applied in physical and engineering problems, for which fast convergence is needed. Applications of gPC-Galerkin schemes to flocking dynamics, and their controlled versions, are almost unexplored in the actual state of art.We give numerical evidence of threshold effects in the alignment dynamic due to the random parameters. In particular the presence of a negative tail in the distribution of the random inputs leads to the divergence of the expected values for the system velocities. The use of a selective model predictive control permits steering of the system towards the desired state even in such unstable regimes.The rest of the paper is organized as follows. In Section2 we introduce briefly a Cucker-Smale dynamic with interaction function depending on stochastic parameters and analyze the system behavior in the case of uniform interactions. The gPC approach is then summarized in Section 3. Subsequently, in Section 4 we consider the gPC scheme in a constrained setting and derive a selective model predictive approximation of the system. Next, in Section 5 we report several numerical experiments which illustrate the different features of the numerical method. Extensions of the present approach are finally discussed in Section 6.


## 2. Cucker-Smale Dynamic with Random Inputs

We introduce a Cucker-Smale type [10] differential system depending on a random variable θ
∈
Ω
⊆
R with a given distribution f
(
θ
). Let (
x
i
,
v
i
)
∈
R
2
d,  d
≥
1, evolving as follows: (1)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
,where K is a time-dependent random function characterizing the uncertainty in the interaction rates and H
(
·
,
·
) is a symmetric function describing the dependence of the alignment dynamic from the agents positions. A classical choice of space-dependent interaction function is related to the distance between two agents (2)
H
x
,
y
=
1
ζ
2
+
x
-
y
2
γ
,where γ
≥
0 and ζ
>
0 are given parameters. Mathematical results concerning the system behavior in the deterministic case (K
≡
1) can be found in [10]. In particular unconditional alignment emerges for γ
<
1
/
2. Let us observe that, even for the model with random inputs (1), the mean velocity of the system is conserved in time (3)
V
θ
,
t
=
1
N
∑
i
=
1
N
v
i
θ
,
t
,
d
d
t
V
θ
,
t
=
0
,since the symmetry of H implies (4)
∑
i
,
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
=
∑
i
,
j
=
1
N
H
x
i
,
x
j
v
i
θ
,
t
.Therefore for each t
≥
0 we have V
(
θ
,
t
)
=
V
(
θ
,
0
).


### 2.1. The Uniform Interaction Case

To better understand the leading dynamic let us consider the simpler uniform interaction case whenH
≡
1, leading to the following equation for the velocities: (5)
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
v
j
θ
,
t
-
v
i
θ
,
t
=
K
θ
,
t
V
θ
,
0
-
v
i
θ
,
t
.The differential equation (5) admits an exact solution depending on the random input θ. More precisely if the initial velocities are deterministically known we have that (6)
v
i
θ
,
t
=
V
+
v
i
0
-
V
exp
⁡
-
∫
0
t
K
θ
,
s
d
s
,where V
=
V
(
0
) is the mean velocity of the system. In what follows we analyze the evolution of (6) for different choices of K
(
θ
,
t
) and of the distribution of the random variable θ.Example 1.
Let us consider a random scattering rate written in terms of the following decomposition:(7)
K
θ
,
t
=
k
θ
h
t
,where h
(
t
) is a nonnegative function depending on t
∈
R
+. The expected velocity of the ith agent is defined by (8)
v
-
i
t
=
E
θ
v
i
θ
,
t
=
∫
Ω
v
i
θ
,
t
f
θ
d
θ
.Then each agent evolves its expected velocity according to (9)
v
-
i
t
=
∫
Ω
V
+
v
i
0
-
V
e
x
p
-
k
θ
∫
0
t
h
s
d
s
·
f
θ
d
θ
.
For example, let us choosek
(
θ
)
=
θ, where the random variable is normally distributed; that is, θ
~
N
(
μ
,
σ
2
). Then, for each i
=
1
,
…
,
N, we need to evaluate the following integral: (10)
V
+
v
i
0
-
V
2
π
σ
2
·
∫
R
e
x
p
-
θ
∫
0
t
h
s
d
s
e
x
p
-
θ
-
μ
2
2
σ
2
d
θ
.The explicit form is easily found through standard techniques and yields(11)
v
-
i
t
=
V
+
v
i
0
-
V
·
e
x
p
-
μ
∫
0
t
h
s
d
s
+
σ
2
2
∫
0
t
h
s
d
s
2
.From (11) we observe a threshold effect in the asymptotic convergence of the mean velocity of each agent toward V. It is immediately seen that if (12)
∫
0
t
h
s
d
s
>
2
μ
σ
2it follows that, for t
→
+
∞, the expected velocity v
-
i diverges. In particular, if h
(
s
)
≡
1 we have that the solution starts to diverge as soon as t
>
μ
/
σ
2. Note that this threshold effect is essential due to the negative tail of the normal distribution. In fact, if we now consider a random variable taking only nonnegative values, for example, exponentially distributed θ
~
E
x
p
(
λ
) for some positive parameter λ
>
0, from (9) we obtain (13)
v
-
i
t
=
V
+
v
i
0
-
V
∫
0
+
∞
e
-
θ
t
λ
e
-
λ
θ
d
θ
,which corresponds to (14)
v
-
i
t
=
V
+
v
i
0
-
V
λ
t
+
λ
,and therefore v
-
i
(
t
)
→
V as t
→
∞. Then independently from the choice of the rate λ
>
0 we obtain for each agent convergences toward the average initial velocity of the system. Finally, in case of a uniform random variable θ
~
U
(
[
a
,
b
]
) we obtain (15)
v
-
i
t
=
V
+
v
i
0
-
V
∫
a
b
1
b
-
a
e
-
θ
t
d
θ
;that is,(16)
v
-
i
t
=
V
+
v
i
0
-
V
b
-
a
e
-
a
t
t
-
e
-
b
t
t
,which implies the divergence of the system in time as soon as a
∈
R assumes negative values.Example 2.
Next we consider a random scattering rate with time-dependent distribution function; that is,(17)
K
θ
,
t
=
θ
twith θ
(
t
)
~
f
(
θ
,
t
). As an example we investigate the case of a normally distributed random parameter with given mean and time-dependent variance: θ
~
N
(
μ
,
σ
2
(
t
)
). It is straightforward to rewrite θ as a translation of a standard normal-distributed variable θ
~; that is, (18)
θ
=
μ
+
σ
t
θ
~
,where θ
~
~
N
(
0,1
). The expected velocities read (19)
v
-
i
t
=
V
+
v
i
0
-
V
2
π
·
∫
R
e
x
p
-
μ
t
-
θ
~
∫
0
t
σ
s
d
s
e
x
p
-
θ
~
2
2
d
θ
~
,which correspond to(20)
v
-
i
t
=
V
+
v
i
0
-
V
e
x
p
-
μ
t
+
1
2
∫
0
t
σ
s
d
s
2
.
Similarly to the case of a time-independent normal variable a threshold effect occurs for large times; that is, the following condition on the variance of the distribution(21)
∫
0
t
σ
s
d
s
2
>
2
μ
timplies the divergence of system (5). As a consequence instability can be avoided by assuming a variance decreasing sufficiently fast in time. The simplest choice is represented by σ
(
t
)
=
1
/
t
α for some α
∈
[
1
/
2,1
). Condition (21) becomes (22)
t
1
-
α
1
-
α
2
>
2
μ
t
.
For example, ifα
=
1
/
2 the previous condition implies that the system diverges for each μ
<
2.


## 2.1. The Uniform Interaction Case

To better understand the leading dynamic let us consider the simpler uniform interaction case whenH
≡
1, leading to the following equation for the velocities: (5)
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
v
j
θ
,
t
-
v
i
θ
,
t
=
K
θ
,
t
V
θ
,
0
-
v
i
θ
,
t
.The differential equation (5) admits an exact solution depending on the random input θ. More precisely if the initial velocities are deterministically known we have that (6)
v
i
θ
,
t
=
V
+
v
i
0
-
V
exp
⁡
-
∫
0
t
K
θ
,
s
d
s
,where V
=
V
(
0
) is the mean velocity of the system. In what follows we analyze the evolution of (6) for different choices of K
(
θ
,
t
) and of the distribution of the random variable θ.Example 1.
Let us consider a random scattering rate written in terms of the following decomposition:(7)
K
θ
,
t
=
k
θ
h
t
,where h
(
t
) is a nonnegative function depending on t
∈
R
+. The expected velocity of the ith agent is defined by (8)
v
-
i
t
=
E
θ
v
i
θ
,
t
=
∫
Ω
v
i
θ
,
t
f
θ
d
θ
.Then each agent evolves its expected velocity according to (9)
v
-
i
t
=
∫
Ω
V
+
v
i
0
-
V
e
x
p
-
k
θ
∫
0
t
h
s
d
s
·
f
θ
d
θ
.
For example, let us choosek
(
θ
)
=
θ, where the random variable is normally distributed; that is, θ
~
N
(
μ
,
σ
2
). Then, for each i
=
1
,
…
,
N, we need to evaluate the following integral: (10)
V
+
v
i
0
-
V
2
π
σ
2
·
∫
R
e
x
p
-
θ
∫
0
t
h
s
d
s
e
x
p
-
θ
-
μ
2
2
σ
2
d
θ
.The explicit form is easily found through standard techniques and yields(11)
v
-
i
t
=
V
+
v
i
0
-
V
·
e
x
p
-
μ
∫
0
t
h
s
d
s
+
σ
2
2
∫
0
t
h
s
d
s
2
.From (11) we observe a threshold effect in the asymptotic convergence of the mean velocity of each agent toward V. It is immediately seen that if (12)
∫
0
t
h
s
d
s
>
2
μ
σ
2it follows that, for t
→
+
∞, the expected velocity v
-
i diverges. In particular, if h
(
s
)
≡
1 we have that the solution starts to diverge as soon as t
>
μ
/
σ
2. Note that this threshold effect is essential due to the negative tail of the normal distribution. In fact, if we now consider a random variable taking only nonnegative values, for example, exponentially distributed θ
~
E
x
p
(
λ
) for some positive parameter λ
>
0, from (9) we obtain (13)
v
-
i
t
=
V
+
v
i
0
-
V
∫
0
+
∞
e
-
θ
t
λ
e
-
λ
θ
d
θ
,which corresponds to (14)
v
-
i
t
=
V
+
v
i
0
-
V
λ
t
+
λ
,and therefore v
-
i
(
t
)
→
V as t
→
∞. Then independently from the choice of the rate λ
>
0 we obtain for each agent convergences toward the average initial velocity of the system. Finally, in case of a uniform random variable θ
~
U
(
[
a
,
b
]
) we obtain (15)
v
-
i
t
=
V
+
v
i
0
-
V
∫
a
b
1
b
-
a
e
-
θ
t
d
θ
;that is,(16)
v
-
i
t
=
V
+
v
i
0
-
V
b
-
a
e
-
a
t
t
-
e
-
b
t
t
,which implies the divergence of the system in time as soon as a
∈
R assumes negative values.Example 2.
Next we consider a random scattering rate with time-dependent distribution function; that is,(17)
K
θ
,
t
=
θ
twith θ
(
t
)
~
f
(
θ
,
t
). As an example we investigate the case of a normally distributed random parameter with given mean and time-dependent variance: θ
~
N
(
μ
,
σ
2
(
t
)
). It is straightforward to rewrite θ as a translation of a standard normal-distributed variable θ
~; that is, (18)
θ
=
μ
+
σ
t
θ
~
,where θ
~
~
N
(
0,1
). The expected velocities read (19)
v
-
i
t
=
V
+
v
i
0
-
V
2
π
·
∫
R
e
x
p
-
μ
t
-
θ
~
∫
0
t
σ
s
d
s
e
x
p
-
θ
~
2
2
d
θ
~
,which correspond to(20)
v
-
i
t
=
V
+
v
i
0
-
V
e
x
p
-
μ
t
+
1
2
∫
0
t
σ
s
d
s
2
.
Similarly to the case of a time-independent normal variable a threshold effect occurs for large times; that is, the following condition on the variance of the distribution(21)
∫
0
t
σ
s
d
s
2
>
2
μ
timplies the divergence of system (5). As a consequence instability can be avoided by assuming a variance decreasing sufficiently fast in time. The simplest choice is represented by σ
(
t
)
=
1
/
t
α for some α
∈
[
1
/
2,1
). Condition (21) becomes (22)
t
1
-
α
1
-
α
2
>
2
μ
t
.
For example, ifα
=
1
/
2 the previous condition implies that the system diverges for each μ
<
2.


## 3. A gPC Based Numerical Approach

In this section we approximate the Cucker-Smale model with random inputs using a generalized polynomial chaos approach. For the sake of clarity we first recall some basic facts concerning gPC approximations.


### 3.1. Preliminaries on gPC Approximations

Let(
Ω
,
F
,
P
) be a probability space, that is, an ordered triple with Ω any set, F a σ-algebra, and P
:
F
→
[
0,1
] a probability measure on F, where we define a random variable (23)
θ
:
Ω
,
F
⟶
R
,
B
R
,with B
R the Borel set of R. Moreover let us consider S
⊂
R
d,  d
≥
1, and [
0
,
T
]
⊂
R certain spatial and temporal subsets. For the sake of simplicity we focus on real-valued functions depending on a single random input (24)
g
x
,
θ
,
t
:
S
×
Ω
×
T
⟶
R
d
,
g
∈
L
2
Ω
,
F
,
P
.In any case it is possible to extend the setup of the problem to ap-dimensional vector of random variables θ
=
(
θ
1
,
…
,
θ
p
); see [19]. Let us consider now the linear space of polynomials of θ of degree up to M, namely, P
M
(
θ
). From classical results in approximation theory it is possible to represent the distribution of random functions with orthogonal polynomials {
Φ
k
(
θ
)
}
k
=
0
M, that is, an orthogonal basis of L
2
(
Ω
,
F
,
P
): (25)
E
θ
Φ
h
θ
Φ
k
θ
=
E
θ
Φ
h
θ
2
δ
h
kwith δ
h
k as the Kronecker delta function. Assuming that the probability law P
(
θ
-
1
(
B
)
)
,
∀
B
∈
B
R, involved in the definition of the introduced function g
(
x
,
θ
,
t
), has finite second-order moment, then the complete polynomial chaos expansion of g is given by (26)
g
x
,
θ
,
t
=
∑
m
∈
N
g
^
m
x
,
t
Φ
m
θ
.According to the Cameron-Martin theorem and to the Askey scheme, results that pave a connection between random variables and orthogonal polynomials [18, 20, 21], we chose a set of polynomials which constitutes the optimal basis with respect to the distribution of the introduced random variable in agreement with Table 1.Table 1
The different gPC choices for the polynomial expansions.
Probability law ofθ
Expansion polynomials
Support
Gaussian
Hermite
(
-
∞
,
+
∞
)
Uniform
Legendre
[
a
,
b
]
Beta
Jacobi
[
a
,
b
]
Gamma
Laguerre
[
0
,
+
∞
)
Poisson
Charlier
NLet us consider now a general formulation for a randomly perturbed problem(27)
D
x
,
t
,
θ
;
g
=
f
x
,
t
,
θ
,where we indicated with D a differential operator. In general the randomness introduced in the problem by θ acts as a perturbation of D, of the function g, or occurs as uncertainty of the initial conditions. In this work we focus on the first two aspects assuming that initial positions and velocities are deterministically known.The generalized polynomial chaos method approximates the solutiong
(
x
,
θ
,
t
) of (27) with its Mth-order polynomial chaos expansion and considers the Galerkin projections of the introduced differential problem; that is, (28)
E
θ
D
x
,
t
,
θ
;
g
·
Φ
h
θ
=
E
θ
f
x
,
θ
,
t
·
Φ
h
θ
,
h
=
0,1
,
…
,
M
.Due to the Galerkin orthogonality between the linear spaceP
M and the error produced in the representation of g
(
x
,
θ
,
t
) with a truncated series, from (28) we obtain a set of M
+
1 purely deterministic equations for the expansion coefficients g
^
m
(
x
,
t
). These subproblems can be solved through classical discretization techniques. From the numerical point of view through a gPC-type method it is possible to achieve an exponential order of convergence to the exact solution of the problem, unlike Monte Carlo techniques for which the order is O
(
1
/
M
) where M is the number of samples.


### 3.2. gPC Approximation of the Alignment Model

We apply the described gPC decomposition to the solution of the nonhomogeneous differential equationv
i
(
θ
,
t
) in (6) and to the stochastic scattering rate K
(
θ
,
t
); that is, (29)
v
i
M
θ
,
t
=
∑
m
=
0
M
v
^
i
,
m
t
Φ
m
θ
,
K
M
θ
,
t
=
∑
l
=
0
M
K
^
l
t
Φ
l
θ
,where (30)
v
^
i
,
m
t
=
E
θ
v
i
θ
,
t
Φ
m
θ
K
^
l
t
=
E
θ
K
θ
,
t
Φ
l
θ
.We obtain the following polynomial chaos expansion:(31)
d
d
t
∑
m
=
0
M
v
^
i
,
m
Φ
m
θ
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
l
,
m
=
0
M
K
^
l
t
v
^
j
,
m
-
v
^
i
,
m
Φ
l
θ
Φ
m
θ
.Multiplying the above expression by an orthogonal element of the basisΦ
h
(
θ
) and integrating with respect to the distribution of θ
(32)
E
θ
∑
m
=
0
M
d
d
t
v
^
i
,
m
Φ
m
θ
Φ
h
θ
=
E
θ
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
l
,
m
=
0
M
K
^
l
t
v
^
j
,
m
-
v
^
i
,
m
Φ
l
θ
Φ
m
θ
Φ
h
θwe find an explicit system of ODEs(33)
d
d
t
v
^
i
,
h
t
=
1
Φ
h
2
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
·
∑
l
=
0
M
K
^
l
t
e
l
m
h
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
K
^
m
h
t
,where e
l
m
h
=
E
θ
[
Φ
l
(
θ
)
Φ
m
(
θ
)
Φ
h
(
θ
)
] and (34)
K
^
m
h
t
=
1
Φ
h
2
∑
l
=
0
M
K
^
l
t
e
l
m
h
.We recall that the gPC numerical approach preserves the mean velocity of the alignment model (5). In fact, from (33) follows(35)
∑
i
=
1
N
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
,
i
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
-
1
N
∑
j
,
i
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
i
,
m
=
0
,thanks to the symmetry of H. More generally it can be shown that if (36)
1
N
∑
i
=
1
N
v
i
θ
,
t
=
V
,where V is time-independent, then its gPC decomposition is also mean-preserving since(37)
1
N
∑
i
=
1
N
∑
m
=
0
M
E
θ
v
i
θ
,
t
Φ
m
θ
Φ
m
θ
=
∑
m
=
0
M
E
θ
1
N
∑
i
=
1
N
v
i
θ
,
t
Φ
m
θ
Φ
m
θ
=
V
∑
m
=
0
M
E
θ
1
·
Φ
m
θ
Φ
m
θ
=
V
.Remark 3.
The gPC approximation (33) can be derived equivalently without expanding the kernel function K
(
θ
,
t
). In this way one obtains (38)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
K
^
m
h
,where now (39)
K
^
m
h
t
=
1
Φ
h
2
E
θ
K
θ
,
t
Φ
m
Φ
h
.
Note that since in generalN
≫
M, the overall computational cost is O
(
M
N
2
).


## 3.1. Preliminaries on gPC Approximations

Let(
Ω
,
F
,
P
) be a probability space, that is, an ordered triple with Ω any set, F a σ-algebra, and P
:
F
→
[
0,1
] a probability measure on F, where we define a random variable (23)
θ
:
Ω
,
F
⟶
R
,
B
R
,with B
R the Borel set of R. Moreover let us consider S
⊂
R
d,  d
≥
1, and [
0
,
T
]
⊂
R certain spatial and temporal subsets. For the sake of simplicity we focus on real-valued functions depending on a single random input (24)
g
x
,
θ
,
t
:
S
×
Ω
×
T
⟶
R
d
,
g
∈
L
2
Ω
,
F
,
P
.In any case it is possible to extend the setup of the problem to ap-dimensional vector of random variables θ
=
(
θ
1
,
…
,
θ
p
); see [19]. Let us consider now the linear space of polynomials of θ of degree up to M, namely, P
M
(
θ
). From classical results in approximation theory it is possible to represent the distribution of random functions with orthogonal polynomials {
Φ
k
(
θ
)
}
k
=
0
M, that is, an orthogonal basis of L
2
(
Ω
,
F
,
P
): (25)
E
θ
Φ
h
θ
Φ
k
θ
=
E
θ
Φ
h
θ
2
δ
h
kwith δ
h
k as the Kronecker delta function. Assuming that the probability law P
(
θ
-
1
(
B
)
)
,
∀
B
∈
B
R, involved in the definition of the introduced function g
(
x
,
θ
,
t
), has finite second-order moment, then the complete polynomial chaos expansion of g is given by (26)
g
x
,
θ
,
t
=
∑
m
∈
N
g
^
m
x
,
t
Φ
m
θ
.According to the Cameron-Martin theorem and to the Askey scheme, results that pave a connection between random variables and orthogonal polynomials [18, 20, 21], we chose a set of polynomials which constitutes the optimal basis with respect to the distribution of the introduced random variable in agreement with Table 1.Table 1
The different gPC choices for the polynomial expansions.
Probability law ofθ
Expansion polynomials
Support
Gaussian
Hermite
(
-
∞
,
+
∞
)
Uniform
Legendre
[
a
,
b
]
Beta
Jacobi
[
a
,
b
]
Gamma
Laguerre
[
0
,
+
∞
)
Poisson
Charlier
NLet us consider now a general formulation for a randomly perturbed problem(27)
D
x
,
t
,
θ
;
g
=
f
x
,
t
,
θ
,where we indicated with D a differential operator. In general the randomness introduced in the problem by θ acts as a perturbation of D, of the function g, or occurs as uncertainty of the initial conditions. In this work we focus on the first two aspects assuming that initial positions and velocities are deterministically known.The generalized polynomial chaos method approximates the solutiong
(
x
,
θ
,
t
) of (27) with its Mth-order polynomial chaos expansion and considers the Galerkin projections of the introduced differential problem; that is, (28)
E
θ
D
x
,
t
,
θ
;
g
·
Φ
h
θ
=
E
θ
f
x
,
θ
,
t
·
Φ
h
θ
,
h
=
0,1
,
…
,
M
.Due to the Galerkin orthogonality between the linear spaceP
M and the error produced in the representation of g
(
x
,
θ
,
t
) with a truncated series, from (28) we obtain a set of M
+
1 purely deterministic equations for the expansion coefficients g
^
m
(
x
,
t
). These subproblems can be solved through classical discretization techniques. From the numerical point of view through a gPC-type method it is possible to achieve an exponential order of convergence to the exact solution of the problem, unlike Monte Carlo techniques for which the order is O
(
1
/
M
) where M is the number of samples.


## 3.2. gPC Approximation of the Alignment Model

We apply the described gPC decomposition to the solution of the nonhomogeneous differential equationv
i
(
θ
,
t
) in (6) and to the stochastic scattering rate K
(
θ
,
t
); that is, (29)
v
i
M
θ
,
t
=
∑
m
=
0
M
v
^
i
,
m
t
Φ
m
θ
,
K
M
θ
,
t
=
∑
l
=
0
M
K
^
l
t
Φ
l
θ
,where (30)
v
^
i
,
m
t
=
E
θ
v
i
θ
,
t
Φ
m
θ
K
^
l
t
=
E
θ
K
θ
,
t
Φ
l
θ
.We obtain the following polynomial chaos expansion:(31)
d
d
t
∑
m
=
0
M
v
^
i
,
m
Φ
m
θ
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
l
,
m
=
0
M
K
^
l
t
v
^
j
,
m
-
v
^
i
,
m
Φ
l
θ
Φ
m
θ
.Multiplying the above expression by an orthogonal element of the basisΦ
h
(
θ
) and integrating with respect to the distribution of θ
(32)
E
θ
∑
m
=
0
M
d
d
t
v
^
i
,
m
Φ
m
θ
Φ
h
θ
=
E
θ
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
l
,
m
=
0
M
K
^
l
t
v
^
j
,
m
-
v
^
i
,
m
Φ
l
θ
Φ
m
θ
Φ
h
θwe find an explicit system of ODEs(33)
d
d
t
v
^
i
,
h
t
=
1
Φ
h
2
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
·
∑
l
=
0
M
K
^
l
t
e
l
m
h
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
K
^
m
h
t
,where e
l
m
h
=
E
θ
[
Φ
l
(
θ
)
Φ
m
(
θ
)
Φ
h
(
θ
)
] and (34)
K
^
m
h
t
=
1
Φ
h
2
∑
l
=
0
M
K
^
l
t
e
l
m
h
.We recall that the gPC numerical approach preserves the mean velocity of the alignment model (5). In fact, from (33) follows(35)
∑
i
=
1
N
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
,
i
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
-
1
N
∑
j
,
i
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
i
,
m
=
0
,thanks to the symmetry of H. More generally it can be shown that if (36)
1
N
∑
i
=
1
N
v
i
θ
,
t
=
V
,where V is time-independent, then its gPC decomposition is also mean-preserving since(37)
1
N
∑
i
=
1
N
∑
m
=
0
M
E
θ
v
i
θ
,
t
Φ
m
θ
Φ
m
θ
=
∑
m
=
0
M
E
θ
1
N
∑
i
=
1
N
v
i
θ
,
t
Φ
m
θ
Φ
m
θ
=
V
∑
m
=
0
M
E
θ
1
·
Φ
m
θ
Φ
m
θ
=
V
.Remark 3.
The gPC approximation (33) can be derived equivalently without expanding the kernel function K
(
θ
,
t
). In this way one obtains (38)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
v
^
j
,
m
-
v
^
i
,
m
K
^
m
h
,where now (39)
K
^
m
h
t
=
1
Φ
h
2
E
θ
K
θ
,
t
Φ
m
Φ
h
.
Note that since in generalN
≫
M, the overall computational cost is O
(
M
N
2
).


## 4. Selective Control of the gPC Approximation

In order to stabilize the gPC approximation of the Cucker-Smale type model (1) with random inputs, we introduce an additional term which acts as control of the approximated dynamic. More specifically we modify the approximation of the alignment model (1) by introducing a control term u
^
h to the hth component of its gPC approximation(40)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
u
^
h
Q
v
^
i
,
h
.The controlu
^
h is a solution of (41)
u
^
h
=
arg
⁡
min
u
^
h
∈
U
b
⁡
1
2
∫
0
T
1
N
∑
i
=
1
N
v
^
i
,
h
t
-
v
^
d
,
h
2
d
t
+
ν
2
∫
0
T
u
^
h
t
2
d
t
,where we assume to have, for each h
=
1
,
…
,
M, a bounded control u
^
h having value in an admissible set U
b
⊂
R
d; for example, in the one-dimensional case we consider u
^
h
∈
[
u
h
,
L
,
u
h
,
R
]. Parameter ν
>
0 is a regularization term and (
v
^
d
,
0
,
v
^
d
,
1
,
…
,
v
^
d
,
M
) are the desired values for the gPC coefficients. For example, (42)
v
^
d
,
h
=
E
θ
v
d
Φ
h
θ
=
v
d
E
Φ
h
θ
=
v
d
h
=
0
0
h
=
1
,
…
,
M
,where v
d is a desired velocity.Moreover the controller action is weighted by a bounded function,(43)
Q
:
R
d
⟶
R
d
.Due to the dependence of the controller effect from the single agent velocity, we refer to this approach as selective control; see [22].In order to tackle numerically the above problem, whose direct solution is prohibitively expansive for large numbers of individuals, we make use of model predictive control (MPC) techniques, also referred to as receding horizon strategy or instantaneous control [23]. These techniques have been used in [8, 9, 22] in the case of deterministic alignment systems.


### 4.1. Selective Model Predictive Control

Let us split the time interval[
0
,
T
] in N
~ time intervals of length Δ
t with t
n
=
n
Δ
t with n
=
0
,
…
,
N
~
-
1. The basic idea of the model predictive control approach is to consider a piecewise constant control (44)
u
^
h
t
=
∑
n
=
0
N
~
-
1
u
-
^
h
n
χ
t
n
,
t
n
+
1
t
.In this way it is possible to determine the value of the controlu
^
h
n
∈
R
d, solving for a state v
^
-
i
,
h the (reduced) optimization problem(45)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
u
^
h
Q
v
^
i
,
h
t
v
^
i
,
h
t
n
=
v
^
-
i
,
h
,
u
^
h
n
=
a
r
g
m
i
n
u
^
h
∈
U
b
∫
t
n
t
n
+
1
1
N
·
∑
i
=
1
N
1
2
v
^
i
,
h
t
-
v
^
d
,
h
2
+
ν
2
u
^
h
t
2
d
s
.Given the controlu
^
h
n on the time interval [
t
n
,
t
n
+
1
], we let v
^
i
,
h evolve according to the dynamics (46)
d
d
t
v
^
i
,
h
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
u
^
h
n
Q
v
^
i
,
h
n
tin order to obtain the new state v
^
-
i
,
h
=
v
^
i
,
h
(
t
n
+
1
). We again solve (45) to obtain u
^
h
n
+
1 with the modified initial data and we repeat this procedure until we reach N
~
Δ
t
=
T.The reduced optimization problem implies a reduction of the complexity of the initial problem since it depends only on the single real-valued variableu
^
h
n. On the other hand the price to pay is that in general the solution of the problem is suboptimal with respect to the full one ((40)-(41)).The quadratic cost and a suitable discretization of (46) allow an explicit representation of u
^
h
n in terms of v
^
-
i
,
h and v
^
i
,
h
n
+
1, as a feedback controlled system, as follows: (47)
v
^
i
,
h
n
+
1
=
v
^
i
,
h
n
+
Δ
t
N
∑
j
=
1
N
H
i
j
n
∑
m
=
0
M
K
^
m
h
n
v
^
j
,
m
n
-
v
^
i
,
m
n
+
Δ
t
u
^
h
n
Q
i
,
h
n
,
v
^
i
,
h
n
=
v
^
-
i
,
h
,
u
^
h
n
=
-
Δ
t
ν
N
∑
i
=
1
N
v
^
i
,
h
n
+
1
-
v
^
d
,
h
Q
i
,
h
n
,where H
i
j
n
≡
H
(
x
i
n
,
x
j
n
) and Q
i
,
h
n
≡
Q
(
v
^
i
,
h
n
). Note that since the feedback control u
^
h
n in (47) depends on the velocities at time n
+
1, the constrained interaction at time n is implicitly defined. The feedback controlled system in the discretized form results in(48)
v
^
i
,
h
n
+
1
=
v
^
i
,
h
n
+
Δ
t
N
∑
j
=
1
N
H
i
j
n
∑
m
=
0
M
K
^
m
h
n
v
^
j
,
m
n
-
v
^
i
,
m
n
-
Δ
t
2
ν
N
∑
j
=
1
N
v
^
j
,
h
n
+
1
-
v
^
d
,
h
Q
j
,
h
n
Q
i
,
h
n
,
v
^
i
,
h
n
=
v
^
-
i
,
h
.Again the action of the control is substituted by an implicit term representing the relaxation toward the desired component of the velocityv
^
d
,
h, and it can be inverted in a fully explicit system.Considering the scaling for the regularization parameterν
=
κ
Δ
t, the previous scheme is a consistent discretization of the following continuous system: (49)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
1
κ
N
∑
j
=
1
N
v
^
d
,
h
-
v
^
j
,
h
t
Q
v
^
j
,
h
t
Q
v
^
i
,
h
t
.Now the control is explicitly embedded in the dynamic of thehth component of the gPC approximation as a feedback term, and the parameter κ
>
0 determines its strength.


### 4.2. Choice of the Selective Control

For the specific choice of weight functionQ
(
·
)
≡
1 we refer in general to nonselective control. Note that in this case the action of the control is not strong enough to control the velocity of each agent; indeed in this case we are able only to control the mean velocity of the system. In fact the control term is reduced to (50)
1
κ
v
^
d
,
h
-
V
^
h
,where V
^
h is the hth coefficient of the expansion of V; that is, (51)
V
^
h
=
1
N
∑
j
=
1
N
v
^
j
,
h
t
.Then, only the projections of the mean velocity are steered toward the respective components of the target velocity; that is, as soon asκ
→
0 it follows that V
^
h
=
v
^
d
,
h. Therefore, the choice of the selective function Q
(
·
) is of paramount importance to ensure the action of the control on the single agent.In principle one can address directly the control problem on the original system (1) as (52)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
u
Q
v
i
θ
,
t
,where the control term u is solution of (53)
u
=
a
r
g
m
i
n
u
∈
U
b
1
2
∫
0
T
1
N
∑
i
=
1
N
v
i
θ
,
t
-
v
d
2
d
t
+
ν
2
∫
0
T
u
t
2
d
t
.Herev
d
∈
R
d is a target velocity, ν
>
0 a regularization parameter, and U
b the set of admissible control. Similarly to previous subsection, through the approach presented in [8, 9, 22], we can derive the time-continuous MPC formulation which explicitly embed the control term in the dynamic as follows:(54)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
1
κ
N
∑
j
=
1
N
v
d
-
v
j
θ
,
t
Q
v
j
t
,
θ
Q
v
i
t
,
θ
.Now the gPC approximation of (54) can be obtained as in Section 3 and leads to the set of ODEs (55)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
1
κ
N
∑
j
=
1
N
R
h
v
i
M
,
v
j
M
,where (56)
R
h
v
i
M
,
v
j
M
=
1
Φ
h
2
E
θ
v
d
-
v
j
M
Q
v
i
M
θ
,
t
·
Q
v
j
M
θ
,
t
Φ
h
θ
.In general systems (55) and (49), without further assumptions on the selective function Q
(
·
), are not equivalent. In addition to the nonselective case, there exists at least one choice of selective control that makes the two approaches totally interchangeable. In fact, taking (57)
Q
v
i
=
v
d
-
v
i
1
/
N
∑
j
=
1
N
v
d
-
v
j
2and Q
(
v
i
)
≡
0 if v
j
=
v
d
,
∀
j
=
1
,
…
,
N, we have that Q
(
·
) is bounded and the control term in (55) takes the following form: (58)
1
κ
N
∑
j
=
1
N
R
h
=
1
κ
Φ
h
2
E
θ
v
d
-
v
i
M
θ
,
t
Φ
h
θ
=
1
κ
v
^
d
,
h
-
v
^
i
,
h
.Similarly the control term in (49) reduces to (59)
1
κ
N
∑
j
=
1
N
v
^
d
,
h
-
v
^
j
,
h
t
Q
v
^
j
,
h
t
Q
v
^
i
,
h
t
=
1
κ
v
^
d
,
h
-
v
^
i
,
h
,and therefore system (55) coincides with (49). Note that as κ
→
0 both systems are driven towards the controlled state v
^
i
,
h
=
v
^
d
,
h which implies a strong control over each single agent.In Figure1 we summarize the two equivalent approaches. In the case of nonselective control and of selective function given by (57) the constrained gPC system can be obtained from our initial unconstrained model (1) through two different but equivalent methods. The first approximates the solution of the Cucker-Smale type model via the gPC projection and then introduces a control on the coefficients of the decomposition through a MPC approach in order to steer each component to (
v
^
d
,
0
,
v
^
d
,
1
,
…
,
v
^
d
,
M
), whereas the second method considers a constrained Cucker-Smale problem (52), introduces its continuous MPC approximation, and then computes the gPC expansion of the resulting system of constrained differential equations.Figure 1
Sketch of the two numerical approaches to solve the control problem with uncertainty, combining MPC and gPC. In both cases, of nonselective control, that is,Q
(
·
)
≡
1, and of selective control with Q
(
·
) defined in (57), the two approaches are equivalent.Remark 4.
We remark that the choice ofQ
(
·
) stated in (57), for which the two approaches sketched in Figure 1 are identical, is equivalent to consider the constrained dynamic (52), modified as follows: (60)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
u
i
,where the control term, u
i, for each agent i
=
1
,
…
,
N, is given by the minimization of the following functional: (61)
J
v
1
,
…
,
v
N
;
u
1
,
…
,
u
N
=
1
2
∫
0
T
1
N
∑
i
=
1
N
v
i
θ
,
t
-
v
d
2
+
ν
2
u
i
t
2
d
t
.Since the functional is strictly convex, applying the (MPC) procedure on a single time interval for the discretized dynamic of (60)-(61), we obtain u
i in terms of feedback control (62)
u
i
=
1
κ
v
d
-
v
i
,
i
=
1
,
…
,
N
.
Thus the same considerations on the equivalence of the approaches hold.


## 4.1. Selective Model Predictive Control

Let us split the time interval[
0
,
T
] in N
~ time intervals of length Δ
t with t
n
=
n
Δ
t with n
=
0
,
…
,
N
~
-
1. The basic idea of the model predictive control approach is to consider a piecewise constant control (44)
u
^
h
t
=
∑
n
=
0
N
~
-
1
u
-
^
h
n
χ
t
n
,
t
n
+
1
t
.In this way it is possible to determine the value of the controlu
^
h
n
∈
R
d, solving for a state v
^
-
i
,
h the (reduced) optimization problem(45)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
·
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
u
^
h
Q
v
^
i
,
h
t
v
^
i
,
h
t
n
=
v
^
-
i
,
h
,
u
^
h
n
=
a
r
g
m
i
n
u
^
h
∈
U
b
∫
t
n
t
n
+
1
1
N
·
∑
i
=
1
N
1
2
v
^
i
,
h
t
-
v
^
d
,
h
2
+
ν
2
u
^
h
t
2
d
s
.Given the controlu
^
h
n on the time interval [
t
n
,
t
n
+
1
], we let v
^
i
,
h evolve according to the dynamics (46)
d
d
t
v
^
i
,
h
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
u
^
h
n
Q
v
^
i
,
h
n
tin order to obtain the new state v
^
-
i
,
h
=
v
^
i
,
h
(
t
n
+
1
). We again solve (45) to obtain u
^
h
n
+
1 with the modified initial data and we repeat this procedure until we reach N
~
Δ
t
=
T.The reduced optimization problem implies a reduction of the complexity of the initial problem since it depends only on the single real-valued variableu
^
h
n. On the other hand the price to pay is that in general the solution of the problem is suboptimal with respect to the full one ((40)-(41)).The quadratic cost and a suitable discretization of (46) allow an explicit representation of u
^
h
n in terms of v
^
-
i
,
h and v
^
i
,
h
n
+
1, as a feedback controlled system, as follows: (47)
v
^
i
,
h
n
+
1
=
v
^
i
,
h
n
+
Δ
t
N
∑
j
=
1
N
H
i
j
n
∑
m
=
0
M
K
^
m
h
n
v
^
j
,
m
n
-
v
^
i
,
m
n
+
Δ
t
u
^
h
n
Q
i
,
h
n
,
v
^
i
,
h
n
=
v
^
-
i
,
h
,
u
^
h
n
=
-
Δ
t
ν
N
∑
i
=
1
N
v
^
i
,
h
n
+
1
-
v
^
d
,
h
Q
i
,
h
n
,where H
i
j
n
≡
H
(
x
i
n
,
x
j
n
) and Q
i
,
h
n
≡
Q
(
v
^
i
,
h
n
). Note that since the feedback control u
^
h
n in (47) depends on the velocities at time n
+
1, the constrained interaction at time n is implicitly defined. The feedback controlled system in the discretized form results in(48)
v
^
i
,
h
n
+
1
=
v
^
i
,
h
n
+
Δ
t
N
∑
j
=
1
N
H
i
j
n
∑
m
=
0
M
K
^
m
h
n
v
^
j
,
m
n
-
v
^
i
,
m
n
-
Δ
t
2
ν
N
∑
j
=
1
N
v
^
j
,
h
n
+
1
-
v
^
d
,
h
Q
j
,
h
n
Q
i
,
h
n
,
v
^
i
,
h
n
=
v
^
-
i
,
h
.Again the action of the control is substituted by an implicit term representing the relaxation toward the desired component of the velocityv
^
d
,
h, and it can be inverted in a fully explicit system.Considering the scaling for the regularization parameterν
=
κ
Δ
t, the previous scheme is a consistent discretization of the following continuous system: (49)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
1
κ
N
∑
j
=
1
N
v
^
d
,
h
-
v
^
j
,
h
t
Q
v
^
j
,
h
t
Q
v
^
i
,
h
t
.Now the control is explicitly embedded in the dynamic of thehth component of the gPC approximation as a feedback term, and the parameter κ
>
0 determines its strength.


## 4.2. Choice of the Selective Control

For the specific choice of weight functionQ
(
·
)
≡
1 we refer in general to nonselective control. Note that in this case the action of the control is not strong enough to control the velocity of each agent; indeed in this case we are able only to control the mean velocity of the system. In fact the control term is reduced to (50)
1
κ
v
^
d
,
h
-
V
^
h
,where V
^
h is the hth coefficient of the expansion of V; that is, (51)
V
^
h
=
1
N
∑
j
=
1
N
v
^
j
,
h
t
.Then, only the projections of the mean velocity are steered toward the respective components of the target velocity; that is, as soon asκ
→
0 it follows that V
^
h
=
v
^
d
,
h. Therefore, the choice of the selective function Q
(
·
) is of paramount importance to ensure the action of the control on the single agent.In principle one can address directly the control problem on the original system (1) as (52)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
u
Q
v
i
θ
,
t
,where the control term u is solution of (53)
u
=
a
r
g
m
i
n
u
∈
U
b
1
2
∫
0
T
1
N
∑
i
=
1
N
v
i
θ
,
t
-
v
d
2
d
t
+
ν
2
∫
0
T
u
t
2
d
t
.Herev
d
∈
R
d is a target velocity, ν
>
0 a regularization parameter, and U
b the set of admissible control. Similarly to previous subsection, through the approach presented in [8, 9, 22], we can derive the time-continuous MPC formulation which explicitly embed the control term in the dynamic as follows:(54)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
1
κ
N
∑
j
=
1
N
v
d
-
v
j
θ
,
t
Q
v
j
t
,
θ
Q
v
i
t
,
θ
.Now the gPC approximation of (54) can be obtained as in Section 3 and leads to the set of ODEs (55)
d
d
t
v
^
i
,
h
t
=
1
N
∑
j
=
1
N
H
x
i
,
x
j
∑
m
=
0
M
K
^
m
h
t
v
^
j
,
m
t
-
v
^
i
,
m
t
+
1
κ
N
∑
j
=
1
N
R
h
v
i
M
,
v
j
M
,where (56)
R
h
v
i
M
,
v
j
M
=
1
Φ
h
2
E
θ
v
d
-
v
j
M
Q
v
i
M
θ
,
t
·
Q
v
j
M
θ
,
t
Φ
h
θ
.In general systems (55) and (49), without further assumptions on the selective function Q
(
·
), are not equivalent. In addition to the nonselective case, there exists at least one choice of selective control that makes the two approaches totally interchangeable. In fact, taking (57)
Q
v
i
=
v
d
-
v
i
1
/
N
∑
j
=
1
N
v
d
-
v
j
2and Q
(
v
i
)
≡
0 if v
j
=
v
d
,
∀
j
=
1
,
…
,
N, we have that Q
(
·
) is bounded and the control term in (55) takes the following form: (58)
1
κ
N
∑
j
=
1
N
R
h
=
1
κ
Φ
h
2
E
θ
v
d
-
v
i
M
θ
,
t
Φ
h
θ
=
1
κ
v
^
d
,
h
-
v
^
i
,
h
.Similarly the control term in (49) reduces to (59)
1
κ
N
∑
j
=
1
N
v
^
d
,
h
-
v
^
j
,
h
t
Q
v
^
j
,
h
t
Q
v
^
i
,
h
t
=
1
κ
v
^
d
,
h
-
v
^
i
,
h
,and therefore system (55) coincides with (49). Note that as κ
→
0 both systems are driven towards the controlled state v
^
i
,
h
=
v
^
d
,
h which implies a strong control over each single agent.In Figure1 we summarize the two equivalent approaches. In the case of nonselective control and of selective function given by (57) the constrained gPC system can be obtained from our initial unconstrained model (1) through two different but equivalent methods. The first approximates the solution of the Cucker-Smale type model via the gPC projection and then introduces a control on the coefficients of the decomposition through a MPC approach in order to steer each component to (
v
^
d
,
0
,
v
^
d
,
1
,
…
,
v
^
d
,
M
), whereas the second method considers a constrained Cucker-Smale problem (52), introduces its continuous MPC approximation, and then computes the gPC expansion of the resulting system of constrained differential equations.Figure 1
Sketch of the two numerical approaches to solve the control problem with uncertainty, combining MPC and gPC. In both cases, of nonselective control, that is,Q
(
·
)
≡
1, and of selective control with Q
(
·
) defined in (57), the two approaches are equivalent.Remark 4.
We remark that the choice ofQ
(
·
) stated in (57), for which the two approaches sketched in Figure 1 are identical, is equivalent to consider the constrained dynamic (52), modified as follows: (60)
x
˙
i
θ
,
t
=
v
i
θ
,
t
v
˙
i
θ
,
t
=
K
θ
,
t
N
∑
j
=
1
N
H
x
i
,
x
j
v
j
θ
,
t
-
v
i
θ
,
t
+
u
i
,where the control term, u
i, for each agent i
=
1
,
…
,
N, is given by the minimization of the following functional: (61)
J
v
1
,
…
,
v
N
;
u
1
,
…
,
u
N
=
1
2
∫
0
T
1
N
∑
i
=
1
N
v
i
θ
,
t
-
v
d
2
+
ν
2
u
i
t
2
d
t
.Since the functional is strictly convex, applying the (MPC) procedure on a single time interval for the discretized dynamic of (60)-(61), we obtain u
i in terms of feedback control (62)
u
i
=
1
κ
v
d
-
v
i
,
i
=
1
,
…
,
N
.
Thus the same considerations on the equivalence of the approaches hold.


## 5. Numerical Tests

We present some numerical experiments of the behavior of the flocking model in the case of a Hermite polynomial chaos expansion. This choice corresponds to the assumption of a normal distribution for the stochastic parameter in the Cucker-Smale type equation (1) and in its constrained behavior (49). Numerical results show that the introduced selective control with the weight function (57) is capable of driving the velocity to a desired state even in case of a dynamic dependent on a normally distributed random input, with fixed or time-dependent variance. In the uniform interaction case, since the effect of agents’ positions does not influence the alignment we report only the results of the agents’ velocities.


### 5.1. Unconstrained Case

In Figures2 and 3 we present numerical results for the convergence of the error using the gPC scheme described in (33) for H
≡
1 and solved through a 4
t
h-order Runge-Kutta method. In particular Figure 2 shows the behavior of the error with respect to increasing terms of the gPC decomposition. Here we considered the average in time of the error for the mean and the variance at time t
>
0 in the L
1 norm (63)
E
v
-
t
=
1
N
∑
i
=
1
N
v
-
i
t
-
v
-
i
M
t
v
-
i
t
E
σ
-
2
t
=
1
N
∑
i
=
1
N
σ
-
i
2
t
-
σ
-
i
2
,
M
t
σ
-
i
2
t
,where (64)
σ
-
i
2
t
=
E
θ
v
i
θ
,
t
-
v
-
i
t
2with v
i
(
θ
,
t
) and v
-
i
(
t
) defined in (6) and (8). Observe that if the scattering rate K
(
θ
,
t
) is of the from described in (7) with h
(
·
)
≡
1 and k
(
θ
)
~
N
(
μ
,
σ
2
) then, in addition to the explicit evolution for the expected velocity as in (9), we can obtain the exact version for the evolution of the variance of the ith agent (65)
σ
-
i
2
t
=
v
i
0
-
V
2
·
e
x
p
-
2
μ
t
+
2
σ
2
t
2
-
e
x
p
-
2
μ
t
+
σ
2
t
2
.Figure 2
Error convergence for increasing number of polynomials in the gPC decomposition approximation. (a) Convergence of the mean error at two fixed timesT
=
1 and T
=
5. (b) Convergence of the variance error. In both cases we considered a random time-independent scattering K
(
θ
,
t
)
=
θ, where the random variable θ is normally distributed N
(
2,1
/
2
). The system of ODEs is solved through a 4
t
h-order Runge-Kutta with Δ
t
=
1
0
-
5.
(a)
(b)Figure 3
Evolution of the variance errorE
σ
-
2
(
t
) defined in (63) for the gPC decomposition for the unconstrained model (5) with K
(
θ
,
t
)
=
θ
~
N
(
2,1
/
2
) over the time interval [
0
,
T
] with T
=
5 and time step Δ
t
=
1
0
-
5.
(a)
(b)In (63) we indicated with σ
-
i
2
,
M
(
t
) the approximated variance (66)
σ
-
i
2
,
M
t
=
∑
h
=
0
M
v
^
i
,
h
2
t
E
θ
Φ
h
θ
2
-
v
^
i
,
0
2
t
.It is easily seen how the error decays spectrally for increasing value ofM; however the method is not capable of going above a certain accuracy and therefore for large M a threshold effect is observed. This can be explained by the large integration interval we have considered in the numerical computation and by the well-known loss of accuracy of gPC for large times [19]. In the case of the error of the variance (Figure 3) the gPC approximation exhibits a slower convergence with respect to the convergence of the mean. Next in Figure 4 we see how for large times the solution of the differential equation (5) diverges and the numerical approximation is capable of describing accurately its behavior only through an increasing number of Hermite polynomials.Figure 4
(a) Sixth-order Hermite gPC decomposition solved through a4
t
h-order Runge-Kutta. (b) 10
t
h-order Hermite gPC decomposition solved through a 4
t
h-order Runge-Kutta. In both cases the final time considered is T
=
6, with time step Δ
t
=
1
0
-
2.
(a)
(b)


### 5.2. Constrained Uniform Interaction Case

In Figure5 we show different scenarios for the uniform interaction dynamic with constraints. In the first row we represent the solution for N
=
10 agents, whose dynamic is described by (49) with v
d
=
1; different values of κ originate different controls on the average of the system, which however do not prevent the system from diverging. In the second row we show the action of selective control (57). It is evident that, with this choice, we are able to control the system also in the case with higher variance.Figure 5
Evolution of the uniform interaction alignment model (49) with N
=
10 agents, at t
=
0 distributed around V
=
2 with unitary variance, depending on a normal random parameter. Left column: θ
~
N
(
2,1
). Right column: θ
~
N
(
2,0.5
). The control term shows its ability to steer the system towards desired velocity v
d
=
1, with different intensities κ
=
1 and κ
=
0.1; when κ
=
∞ the control has no influence. First row shows the action of the control acting just on the average velocity, Q
≡
1. Second row shows the action of selective control with Q
(
·
) as in (57).
(a)
(b)
(c)
(d)Observe that the numerical results are coherent with the explicit solution of the controlled equation. Let us consider the time-independent scattering rateK
(
θ
,
t
)
=
θ
~
N
(
μ
,
σ
2
); then from the equation (67)
d
d
t
v
i
θ
,
t
=
θ
V
-
v
i
θ
,
t
+
1
κ
v
d
-
v
i
θ
,
twe can compute the exact solution given v
i
(
θ
,
0
)
=
v
i
(
0
): (68)
v
i
θ
,
t
=
κ
V
θ
+
v
d
κ
θ
+
1
+
v
i
0
-
κ
V
θ
+
v
d
κ
θ
+
1
e
x
p
-
θ
+
1
κ
t
.The asymptotic behavior of the expected value of (68) can be studied similarly to what we did in Section 2.1. In other words in order to prevent the divergence of the leading term of the controlled expected exact solution we might study (69)
e
x
p
-
μ
+
1
κ
t
+
σ
2
t
2
2
,which diverge if (70)
t
>
2
σ
2
μ
+
1
κ
.Then for each fixed time we could select a regularization parameterκ
>
0 so as to avoid the divergence of (68). Moreover we can observe that in the limit κ
→
0 the introduced selective control is capable of correctly driving the system (67) for each t
>
0.In Figure6 we consider the system with random time-dependent scattering rate θ
~
N
(
μ
,
σ
2
(
t
)
). The dynamic shows how, for the choice of time-dependent variance described within Example 2 in Section 2.1, that is, σ
(
t
)
=
1
/
t
α with α
=
1
/
2, the convergence depends on the mean value of the random input. In particular numerical experiments highlight the threshold effect for μ
=
2 which we derived in Section 2. In the second figure we show that the action of the selective control (57), with desired velocity v
d
=
V, is capable of stabilizing the system and driving the velocities towards the desired state.Figure 6
Solution of the uniform interaction case with time-dependent random parameterθ distributed accordingly to a normal distribution N
(
μ
,
σ
2
(
t
)
), with a time-dependent standard deviation σ
(
t
)
=
1
/
t
α, and α
=
1
/
2. (a) We see the threshold for different values of μ; that is, for μ
<
2 the system diverges. (b) Solution of the constrained model with κ
=
0.1; observe that we are able to steer the system to the desired velocity v
d
=
V, that is, the initial mean velocity of the system, using the selective control described in (57).
(a)
(b)


### 5.3. Constrained Space-Dependent Case

Next let us consider the full space nonhomogeneous constrained problem (1) with the interaction function defined in (2). In this case we assume that K
(
θ
)
=
θ with θ
~
N
(
μ
,
σ
2
). In Figures 7 and 8 we consider a system of N
=
100 agents with Gaussian initial position with zero mean and with variance 2 and Gaussian initial velocities clustered around ±
5 with variance 1
/
10. The numerical results for (33) have been performed through a 10
t
h-order gPC expansion. The dynamic has been observed in the time interval [
0,5
] with Δ
t
=
1
0
-
2. In Figure 8 we see how the selective control is capable of driving the velocity of each agent to the desired state v
d. In fact in case of no control (see Figure 7) we have that the velocities of the system naturally diverge.Figure 7
Numerical solution of (49), with γ
=
0.05
<
1
/
2, ζ
=
0.01, through a 10
t
h-order gPC Hermite decomposition (33) with κ
=
∞ with time step Δ
t
=
1
0
-
2. The random input is normally distributed θ
~
N
(
2,1
).
(a)
t
=
0
(b)
t
=
1
(c)
t
=
2
(d)
t
=
3
(e)
t
=
4
(f)
t
=
5Figure 8
Numerical solution of (49), with γ
=
0.05
<
1
/
2, ζ
=
0.01, through a 6th-order gPC Hermite decomposition for the selective control (33) with time step Δ
t
=
1
0
-
2. Here we considered a normally distributed random input θ
~
N
(
2,1
); the desired velocity is v
d
=
0 and the control parameter is κ
=
1.
(a)
t
=
0
(b)
t
=
1
(c)
t
=
2
(d)
t
=
3
(e)
t
=
4
(f)
t
=
5


## 5.1. Unconstrained Case

In Figures2 and 3 we present numerical results for the convergence of the error using the gPC scheme described in (33) for H
≡
1 and solved through a 4
t
h-order Runge-Kutta method. In particular Figure 2 shows the behavior of the error with respect to increasing terms of the gPC decomposition. Here we considered the average in time of the error for the mean and the variance at time t
>
0 in the L
1 norm (63)
E
v
-
t
=
1
N
∑
i
=
1
N
v
-
i
t
-
v
-
i
M
t
v
-
i
t
E
σ
-
2
t
=
1
N
∑
i
=
1
N
σ
-
i
2
t
-
σ
-
i
2
,
M
t
σ
-
i
2
t
,where (64)
σ
-
i
2
t
=
E
θ
v
i
θ
,
t
-
v
-
i
t
2with v
i
(
θ
,
t
) and v
-
i
(
t
) defined in (6) and (8). Observe that if the scattering rate K
(
θ
,
t
) is of the from described in (7) with h
(
·
)
≡
1 and k
(
θ
)
~
N
(
μ
,
σ
2
) then, in addition to the explicit evolution for the expected velocity as in (9), we can obtain the exact version for the evolution of the variance of the ith agent (65)
σ
-
i
2
t
=
v
i
0
-
V
2
·
e
x
p
-
2
μ
t
+
2
σ
2
t
2
-
e
x
p
-
2
μ
t
+
σ
2
t
2
.Figure 2
Error convergence for increasing number of polynomials in the gPC decomposition approximation. (a) Convergence of the mean error at two fixed timesT
=
1 and T
=
5. (b) Convergence of the variance error. In both cases we considered a random time-independent scattering K
(
θ
,
t
)
=
θ, where the random variable θ is normally distributed N
(
2,1
/
2
). The system of ODEs is solved through a 4
t
h-order Runge-Kutta with Δ
t
=
1
0
-
5.
(a)
(b)Figure 3
Evolution of the variance errorE
σ
-
2
(
t
) defined in (63) for the gPC decomposition for the unconstrained model (5) with K
(
θ
,
t
)
=
θ
~
N
(
2,1
/
2
) over the time interval [
0
,
T
] with T
=
5 and time step Δ
t
=
1
0
-
5.
(a)
(b)In (63) we indicated with σ
-
i
2
,
M
(
t
) the approximated variance (66)
σ
-
i
2
,
M
t
=
∑
h
=
0
M
v
^
i
,
h
2
t
E
θ
Φ
h
θ
2
-
v
^
i
,
0
2
t
.It is easily seen how the error decays spectrally for increasing value ofM; however the method is not capable of going above a certain accuracy and therefore for large M a threshold effect is observed. This can be explained by the large integration interval we have considered in the numerical computation and by the well-known loss of accuracy of gPC for large times [19]. In the case of the error of the variance (Figure 3) the gPC approximation exhibits a slower convergence with respect to the convergence of the mean. Next in Figure 4 we see how for large times the solution of the differential equation (5) diverges and the numerical approximation is capable of describing accurately its behavior only through an increasing number of Hermite polynomials.Figure 4
(a) Sixth-order Hermite gPC decomposition solved through a4
t
h-order Runge-Kutta. (b) 10
t
h-order Hermite gPC decomposition solved through a 4
t
h-order Runge-Kutta. In both cases the final time considered is T
=
6, with time step Δ
t
=
1
0
-
2.
(a)
(b)


## 5.2. Constrained Uniform Interaction Case

In Figure5 we show different scenarios for the uniform interaction dynamic with constraints. In the first row we represent the solution for N
=
10 agents, whose dynamic is described by (49) with v
d
=
1; different values of κ originate different controls on the average of the system, which however do not prevent the system from diverging. In the second row we show the action of selective control (57). It is evident that, with this choice, we are able to control the system also in the case with higher variance.Figure 5
Evolution of the uniform interaction alignment model (49) with N
=
10 agents, at t
=
0 distributed around V
=
2 with unitary variance, depending on a normal random parameter. Left column: θ
~
N
(
2,1
). Right column: θ
~
N
(
2,0.5
). The control term shows its ability to steer the system towards desired velocity v
d
=
1, with different intensities κ
=
1 and κ
=
0.1; when κ
=
∞ the control has no influence. First row shows the action of the control acting just on the average velocity, Q
≡
1. Second row shows the action of selective control with Q
(
·
) as in (57).
(a)
(b)
(c)
(d)Observe that the numerical results are coherent with the explicit solution of the controlled equation. Let us consider the time-independent scattering rateK
(
θ
,
t
)
=
θ
~
N
(
μ
,
σ
2
); then from the equation (67)
d
d
t
v
i
θ
,
t
=
θ
V
-
v
i
θ
,
t
+
1
κ
v
d
-
v
i
θ
,
twe can compute the exact solution given v
i
(
θ
,
0
)
=
v
i
(
0
): (68)
v
i
θ
,
t
=
κ
V
θ
+
v
d
κ
θ
+
1
+
v
i
0
-
κ
V
θ
+
v
d
κ
θ
+
1
e
x
p
-
θ
+
1
κ
t
.The asymptotic behavior of the expected value of (68) can be studied similarly to what we did in Section 2.1. In other words in order to prevent the divergence of the leading term of the controlled expected exact solution we might study (69)
e
x
p
-
μ
+
1
κ
t
+
σ
2
t
2
2
,which diverge if (70)
t
>
2
σ
2
μ
+
1
κ
.Then for each fixed time we could select a regularization parameterκ
>
0 so as to avoid the divergence of (68). Moreover we can observe that in the limit κ
→
0 the introduced selective control is capable of correctly driving the system (67) for each t
>
0.In Figure6 we consider the system with random time-dependent scattering rate θ
~
N
(
μ
,
σ
2
(
t
)
). The dynamic shows how, for the choice of time-dependent variance described within Example 2 in Section 2.1, that is, σ
(
t
)
=
1
/
t
α with α
=
1
/
2, the convergence depends on the mean value of the random input. In particular numerical experiments highlight the threshold effect for μ
=
2 which we derived in Section 2. In the second figure we show that the action of the selective control (57), with desired velocity v
d
=
V, is capable of stabilizing the system and driving the velocities towards the desired state.Figure 6
Solution of the uniform interaction case with time-dependent random parameterθ distributed accordingly to a normal distribution N
(
μ
,
σ
2
(
t
)
), with a time-dependent standard deviation σ
(
t
)
=
1
/
t
α, and α
=
1
/
2. (a) We see the threshold for different values of μ; that is, for μ
<
2 the system diverges. (b) Solution of the constrained model with κ
=
0.1; observe that we are able to steer the system to the desired velocity v
d
=
V, that is, the initial mean velocity of the system, using the selective control described in (57).
(a)
(b)


## 5.3. Constrained Space-Dependent Case

Next let us consider the full space nonhomogeneous constrained problem (1) with the interaction function defined in (2). In this case we assume that K
(
θ
)
=
θ with θ
~
N
(
μ
,
σ
2
). In Figures 7 and 8 we consider a system of N
=
100 agents with Gaussian initial position with zero mean and with variance 2 and Gaussian initial velocities clustered around ±
5 with variance 1
/
10. The numerical results for (33) have been performed through a 10
t
h-order gPC expansion. The dynamic has been observed in the time interval [
0,5
] with Δ
t
=
1
0
-
2. In Figure 8 we see how the selective control is capable of driving the velocity of each agent to the desired state v
d. In fact in case of no control (see Figure 7) we have that the velocities of the system naturally diverge.Figure 7
Numerical solution of (49), with γ
=
0.05
<
1
/
2, ζ
=
0.01, through a 10
t
h-order gPC Hermite decomposition (33) with κ
=
∞ with time step Δ
t
=
1
0
-
2. The random input is normally distributed θ
~
N
(
2,1
).
(a)
t
=
0
(b)
t
=
1
(c)
t
=
2
(d)
t
=
3
(e)
t
=
4
(f)
t
=
5Figure 8
Numerical solution of (49), with γ
=
0.05
<
1
/
2, ζ
=
0.01, through a 6th-order gPC Hermite decomposition for the selective control (33) with time step Δ
t
=
1
0
-
2. Here we considered a normally distributed random input θ
~
N
(
2,1
); the desired velocity is v
d
=
0 and the control parameter is κ
=
1.
(a)
t
=
0
(b)
t
=
1
(c)
t
=
2
(d)
t
=
3
(e)
t
=
4
(f)
t
=
5


## 6. Conclusions

We proposed a general approach for the numerical approximation of flocking models with random inputs through gPC. The method is constructed in two steps. First, the random Cucker-Smale system is solved by gPC. The presence of uncertainty in the interaction terms, which is a natural assumption in this kind of problems, leads to threshold effects in the asymptotic behavior of the system. Next a constrained gPC approximation is introduced and approximated though a selective model predictive control strategy. Relations under which the introduction of the gPC approximation and the model predictive control commute are also derived. The numerical examples illustrate that the assumption of positivity of the mean value of the random input is not sufficient for the alignment of the system but that a suitable choice of the selective control is capable of stabilizing the system towards the desired state. Extension of this technique to the case of a large number of interacting agents through mean-field and Boltzmann approximations is actually under study.



---

*Source: 850124-2015-06-01.xml*