Uncertainty Quantification in Control Problems for Flocking Models

In this paper 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 to steer the system towards the desired state even in unstable regimes.


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][5][6][7] and their kinetic and mean-field approximations [8][9][10][11][12][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][15][16][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 Section 2 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 2 Mathematical Problems in Engineering 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.

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 ( ). Let ( , V ) ∈ R 2 , ≥ 1, evolving as follows: where is a time-dependent random function characterizing the uncertainty in the interaction rates and (⋅, ⋅) 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 ( , ) = 1 where ≥ 0 and > 0 are given parameters. Mathematical results concerning the system behavior in the deterministic case ( ≡ 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 since the symmetry of implies Therefore for each ≥ 0 we have V( , ) = V( , 0).

The Uniform Interaction Case.
To better understand the leading dynamic let us consider the simpler uniform interaction case when ≡ 1, leading to the following equation for the velocities: 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 where V = V(0) is the mean velocity of the system. In what follows we analyze the evolution of (6) for different choices of ( , ) and of the distribution of the random variable .
Example 1. Let us consider a random scattering rate written in terms of the following decomposition: where ℎ( ) is a nonnegative function depending on ∈ R + . The expected velocity of the th agent is defined by Then each agent evolves its expected velocity according to For example, let us choose ( ) = , where the random variable is normally distributed; that is, ∼ N( , 2 ). Then, for each = 1, . . . , , we need to evaluate the following integral: The explicit form is easily found through standard techniques and yields 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 it follows that, for → +∞, the expected velocity V diverges. In particular, if ℎ( ) ≡ 1 we have that the solution starts to diverge as soon as > / 2 . Note that this threshold effect is essential due to the negative tail of the normal distribution. Engineering   3 In fact, if we now consider a random variable taking only nonnegative values, for example, exponentially distributed ∼ Exp( ) for some positive parameter > 0, from (9) we obtain

Mathematical Problems in
which corresponds to and therefore V ( ) → V as → ∞. 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 ∼ ([ , ]) we obtain that is, which implies the divergence of the system in time as soon as ∈ R assumes negative values.
Example 2. Next we consider a random scattering rate with time-dependent distribution function; that is, with ( ) ∼ ( , ). As an example we investigate the case of a normally distributed random parameter with given mean and time-dependent variance: ∼ N( , 2 ( )). It is straightforward to rewrite as a translation of a standard normal-distributed variablẽ; that is, wherẽ∼ N(0, 1). The expected velocities read which correspond to 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 implies 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 ( ) = 1/ for some ∈ [1/2, 1). Condition (21) becomes For example, if = 1/2 the previous condition implies that the system diverges for each < 2.

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.

Preliminaries on gPC
Approximations. Let (Ω, F, ) be a probability space, that is, an ordered triple with Ω any set, F a -algebra, and : F → [0, 1] a probability measure on F, where we define a random variable with B R the Borel set of R. Moreover let us consider ⊂ R , ≥ 1, and [0, ] ⊂ R certain spatial and temporal subsets. For the sake of simplicity we focus on real-valued functions depending on a single random input In any case it is possible to extend the setup of the problem to a -dimensional vector of random variables = ( 1 , . . . , ); see [19]. Let us consider now the linear space of polynomials of of degree up to , namely, P ( ). From classical results in approximation theory it is possible to represent the distribution of random functions with orthogonal polynomials {Φ ( )} =0 , that is, an orthogonal basis of 2 (Ω, F, ): with ℎ as the Kronecker delta function. Assuming that the probability law ( −1 ( )), ∀ ∈ B R , involved in the definition of the introduced function ( , , ), has finite second-order moment, then the complete polynomial chaos expansion of is given by 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.
Let us consider now a general formulation for a randomly perturbed problem 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 , 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 solution ( , , ) of (27) with its th-order polynomial chaos expansion and considers the Galerkin projections of the introduced differential problem; that is, Due to the Galerkin orthogonality between the linear space P and the error produced in the representation of ( , , ) with a truncated series, from (28) we obtain a set of + 1 purely deterministic equations for the expansion coef-ficientŝ( , ). 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 (1/ √ ) where is the number of samples.

gPC Approximation of the Alignment Model.
We apply the described gPC decomposition to the solution of the nonhomogeneous differential equation V ( , ) in (6) and to the stochastic scattering rate ( , ); that is, We obtain the following polynomial chaos expansion: Hermite Charlier N Multiplying the above expression by an orthogonal element of the basis Φ ℎ ( ) and integrating with respect to the distribution of we find an explicit system of ODEŝ We recall that the gPC numerical approach preserves the mean velocity of the alignment model (5). In fact, from (33) follows thanks to the symmetry of . More generally it can be shown that if Mathematical Problems in Engineering 5 where V is time-independent, then its gPC decomposition is also mean-preserving since Remark 3. The gPC approximation (33) can be derived equivalently without expanding the kernel function ( , ).
In this way one obtainŝ where noŵℎ Note that since in general ≫ , the overall computational cost is ( 2 ).

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̂ℎ to the ℎth component of its gPC approximation The control̂ℎ is a solution of where we assume to have, for each ℎ = 1, . . . , , a bounded control̂ℎ having value in an admissible set U ⊂ R ; for example, in the one-dimensional case we consider ℎ ∈ [ ℎ, , ℎ, ]. Parameter ] > 0 is a regularization term and (V ,0 ,V ,1 , . . . ,V , ) are the desired values for the gPC coefficients. For example, where V is a desired velocity. Moreover the controller action is weighted by a bounded function, 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.

Selective Model Predictive Control.
Let us split the time interval [0, ] iñtime intervals of length Δ with = Δ with = 0, . . . ,̃− 1. The basic idea of the model predictive control approach is to consider a piecewise constant control In this way it is possible to determine the value of the control̂ℎ ∈ R , solving for a stateV ,ℎ the (reduced) optimization problem Given the control̂ℎ on the time interval [ , +1 ], we let V ,ℎ evolve according to the dynamicŝ 6 Mathematical Problems in Engineering in order to obtain the new stateV ,ℎ =V ,ℎ ( +1 ). We again solve (45) to obtain̂+ 1 ℎ with the modified initial data and we repeat this procedure until we reach̃Δ = .
The reduced optimization problem implies a reduction of the complexity of the initial problem since it depends only on the single real-valued variablêℎ. 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̂ℎ in terms ofV ,ℎ andV +1 ,ℎ , as a feedback controlled system, as follows: where ≡ ( , ) and ,ℎ ≡ (V ,ℎ ). Note that since the feedback control̂ℎ in (47) depends on the velocities at time + 1, the constrained interaction at time is implicitly defined. The feedback controlled system in the discretized form results in Again the action of the control is substituted by an implicit term representing the relaxation toward the desired component of the velocityV ,ℎ , and it can be inverted in a fully explicit system.
Considering the scaling for the regularization parameter ] = Δ , the previous scheme is a consistent discretization of the following continuous system: Now the control is explicitly embedded in the dynamic of the ℎth component of the gPC approximation as a feedback term, and the parameter > 0 determines its strength.

Choice of the Selective Control.
For the specific choice of weight function (⋅) ≡ 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 whereV ℎ is the ℎth coefficient of the expansion of V; that is, 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 thatV ℎ =V ,ℎ . Therefore, the choice of the selective function (⋅) 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) aṡ where the control term is solution of Here V ∈ R is a target velocity, ] > 0 a regularization parameter, and U 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: Mathematical Now the gPC approximation of (54) can be obtained as in Section 3 and leads to the set of ODEŝ where In general systems (55) and (49), without further assumptions on the selective function (⋅), 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 and (V ) ≡ 0 if V = V , ∀ = 1, . . . , , we have that (⋅) is bounded and the control term in (55) takes the following form: Similarly the control term in (49) reduces to and therefore system (55) coincides with (49). Note that as → 0 both systems are driven towards the controlled statê V ,ℎ =V ,ℎ which implies a strong control over each single agent.
In Figure 1 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 ,0 ,V ,1 , . . . ,V , ), 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.

Remark 4.
We remark that the choice of (⋅) 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: where the control term, , for each agent = 1, . . . , , is given by the minimization of the following functional: Since the functional is strictly convex, applying the (MPC) procedure on a single time interval for the discretized dynamic of (60)-(61), we obtain in terms of feedback control Thus the same considerations on the equivalence of the approaches hold.

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 Mathematical Problems in Engineering influence the alignment we report only the results of the agents' velocities. Figures 2 and 3 we present numerical results for the convergence of the error using the gPC scheme described in (33) for ≡ 1 and solved through a 4th-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 > 0 in the 1 norm

Unconstrained Case. In
with V ( , ) and V ( ) defined in (6) and (8). Observe that if the scattering rate ( , ) is of the from described in (7) with ℎ(⋅) ≡ 1 and ( ) ∼ 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 th agent In (63) we indicated with 2, ( ) the approximated variance It is easily seen how the error decays spectrally for increasing value of ; however the method is not capable of going above a certain accuracy and therefore for large 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 5 we show different scenarios for the uniform interaction dynamic with constraints. In the first row we represent the solution for = 10 agents, whose dynamic is described by (49) with V = 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. Observe that the numerical results are coherent with the explicit solution of the controlled equation. Let us consider the time-independent scattering rate ( , ) = ∼ N( , 2 ); then from the equation

Constrained Uniform Interaction Case. In
we can compute the exact solution given V ( , 0) = V (0):  Figure 6: Solution of the uniform interaction case with time-dependent random parameter distributed accordingly to a normal distribution N( , 2 ( )), with a time-dependent standard deviation ( ) = 1/ , 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 = V, that is, the initial mean velocity of the system, using the selective control described in (57).
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 which diverge if 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 > 0.
In Figure 6 we consider the system with random timedependent scattering rate ∼ N( , 2 ( )). The dynamic shows how, for the choice of time-dependent variance described within Example 2 in Section 2.1, that is, ( ) = 1/ 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 = V, is capable of stabilizing the system and driving the velocities towards the desired state.

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 ( ) = with ∼ N( , 2 ). In Figures 7 and 8 we consider a system of = 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 10th-order gPC expansion. The dynamic has been observed in the time interval [0, 5] with Δ = 10 −2 . In Figure 8 we see how the selective control is capable of driving the velocity of each agent to the desired state V . In fact in case of no control (see Figure 7) we have that the velocities of the system naturally diverge.

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 Here we considered a normally distributed random input ∼ N(2, 1); the desired velocity is V = 0 and the control parameter is = 1.