Research Article  Open Access
Uncertainty Quantification in Control Problems for Flocking Models
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 CuckerSmale 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.
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 meanfield 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 selforganized 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 gPCGalerkin 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 CuckerSmale 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. CuckerSmale Dynamic with Random Inputs
We introduce a CuckerSmale type [10] differential system depending on a random variable with a given distribution . Let , , evolving as follows: where is a timedependent 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 spacedependent interaction function is related to the distance between two agents where and are given parameters. Mathematical results concerning the system behavior in the deterministic case () can be found in [10]. In particular unconditional alignment emerges for . 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 we have .
2.1. The Uniform Interaction Case
To better understand the leading dynamic let us consider the simpler uniform interaction case when , 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 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 . 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, . Then, for each , 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 . It is immediately seen that if it follows that, for , the expected velocity diverges. In particular, if we have that the solution starts to diverge as soon as . 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 for some positive parameter , from (9) we obtain which corresponds to and therefore as . Then independently from the choice of the rate 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 assumes negative values.
Example 2. Next we consider a random scattering rate with timedependent distribution function; that is, with . As an example we investigate the case of a normally distributed random parameter with given mean and timedependent variance: . It is straightforward to rewrite as a translation of a standard normaldistributed variable ; that is, where . The expected velocities read which correspond toSimilarly to the case of a timeindependent 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 for some . Condition (21) becomes For example, if the previous condition implies that the system diverges for each .
3. A gPC Based Numerical Approach
In this section we approximate the CuckerSmale 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 be a probability space, that is, an ordered triple with any set, a algebra, and a probability measure on , where we define a random variable with the Borel set of . Moreover let us consider , , and certain spatial and temporal subsets. For the sake of simplicity we focus on realvalued 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 ; see [19]. Let us consider now the linear space of polynomials of of degree up to , namely, . From classical results in approximation theory it is possible to represent the distribution of random functions with orthogonal polynomials , that is, an orthogonal basis of : with as the Kronecker delta function. Assuming that the probability law , involved in the definition of the introduced function , has finite secondorder moment, then the complete polynomial chaos expansion of is given by
According to the CameronMartin 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 a differential operator. In general the randomness introduced in the problem by acts as a perturbation of , 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 thorder polynomial chaos expansion and considers the Galerkin projections of the introduced differential problem; that is,
Due to the Galerkin orthogonality between the linear space and the error produced in the representation of with a truncated series, from (28) we obtain a set of purely deterministic equations for the expansion coefficients . These subproblems can be solved through classical discretization techniques. From the numerical point of view through a gPCtype 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 where 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 equation in (6) and to the stochastic scattering rate ; that is, where
We obtain the following polynomial chaos expansion:
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 ODEswhere and
We recall that the gPC numerical approach preserves the mean velocity of the alignment model (5). In fact, from (33) followsthanks to the symmetry of . More generally it can be shown that if where is timeindependent, then its gPC decomposition is also meanpreserving since
Remark 3. The gPC approximation (33) can be derived equivalently without expanding the kernel function . In this way one obtains where now Note that since in general , the overall computational cost is .
4. Selective Control of the gPC Approximation
In order to stabilize the gPC approximation of the CuckerSmale 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 , a bounded control having value in an admissible set ; for example, in the onedimensional case we consider . Parameter is a regularization term and are the desired values for the gPC coefficients. For example, where 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.
4.1. Selective Model Predictive Control
Let us split the time interval in time intervals of length with with . 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 , solving for a state the (reduced) optimization problem
Given the control on the time interval , we let evolve according to the dynamics in order to obtain the new state . We again solve (45) to obtain 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 realvalued variable . 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 of and , as a feedback controlled system, as follows: where and . Note that since the feedback control in (47) depends on the velocities at time , 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 velocity , 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 determines its strength.
4.2. Choice of the Selective Control
For the specific choice of weight function 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 where is the th coefficient of the expansion of ; 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 it follows that . 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) as where the control term is solution of
Here is a target velocity, a regularization parameter, and the set of admissible control. Similarly to previous subsection, through the approach presented in [8, 9, 22], we can derive the timecontinuous MPC formulation which explicitly embed the control term in the dynamic as follows:
Now the gPC approximation of (54) can be obtained as in Section 3 and leads to the set of ODEs 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 if , 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 both systems are driven towards the controlled state 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 CuckerSmale 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 , whereas the second method considers a constrained CuckerSmale 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 , 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.
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 CuckerSmale 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 timedependent 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 Figures 2 and 3 we present numerical results for the convergence of the error using the gPC scheme described in (33) for and solved through a order RungeKutta 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 in the norm where with and defined in (6) and (8). Observe that if the scattering rate is of the from described in (7) with and 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
(a)
(b)
(a)
(b)
In (63) we indicated with 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 wellknown 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.
(a)
(b)
5.2. Constrained Uniform Interaction Case
In Figure 5 we show different scenarios for the uniform interaction dynamic with constraints. In the first row we represent the solution for agents, whose dynamic is described by (49) with ; 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.
(a)
(b)
(c)
(d)
Observe that the numerical results are coherent with the explicit solution of the controlled equation. Let us consider the timeindependent scattering rate ; then from the equation we can compute the exact solution given :
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 so as to avoid the divergence of (68). Moreover we can observe that in the limit the introduced selective control is capable of correctly driving the system (67) for each .
In Figure 6 we consider the system with random timedependent scattering rate . The dynamic shows how, for the choice of timedependent variance described within Example 2 in Section 2.1, that is, with , the convergence depends on the mean value of the random input. In particular numerical experiments highlight the threshold effect for which we derived in Section 2. In the second figure we show that the action of the selective control (57), with desired velocity , is capable of stabilizing the system and driving the velocities towards the desired state.
(a)
(b)
5.3. Constrained SpaceDependent 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 . In Figures 7 and 8 we consider a system of agents with Gaussian initial position with zero mean and with variance and Gaussian initial velocities clustered around with variance . The numerical results for (33) have been performed through a order gPC expansion. The dynamic has been observed in the time interval with . In Figure 8 we see how the selective control is capable of driving the velocity of each agent to the desired state . In fact in case of no control (see Figure 7) we have that the velocities of the system naturally diverge.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
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 CuckerSmale 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 meanfield and Boltzmann approximations is actually under study.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Giacomo Albi acknowledges the support of the ERCStarting Grant HighDimensional Sparse Optimal Control (HDSPCONTR306274).
References
 L. Pareschi and G. Toscani, Interacting Multiagent Systems: Kinetic Equations & Monte Carlo Methods, Oxford University Press, New York, NY, USA, 2013.
 C. W. Reynolds, “Flocks, herds and schools: a distributed behavioral mode,” Computer Graphics, vol. 21, no. 4, pp. 25–34, 1987. View at: Publisher Site  Google Scholar
 T. Vicsek, A. Czirók, E. BenJacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of selfdriven particles,” Physical Review Letters, vol. 75, no. 6, pp. 1226–1229, 1995. View at: Publisher Site  Google Scholar
 G. Albi, D. Balagué, J. A. Carrillo, and J. von Brecht, “Stability analysis of flock and mill rings for second order models in swarming,” SIAM Journal on Applied Mathematics, vol. 74, no. 3, pp. 794–818, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 J. A. Carrillo, M. Fornasier, J. Rosado, and G. Toscani, “Asymptotic flocking dynamics for the kinetic CuckerSmale model,” SIAM Journal on Mathematical Analysis, vol. 42, no. 1, pp. 218–236, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 S. Motsch and E. Tadmor, “Heterophilious dynamics enhances consensus,” SIAM Review, vol. 56, no. 4, pp. 577–621, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 H. Zhou, W. Zhou, and W. Zeng, “Flocking control of multiple mobile agents with the rules of avoiding collision,” Mathematical Problems in Engineering. In press. View at: Google Scholar
 G. Albi, M. Herty, and L. Pareschi, “Kinetic description of optimal control problems in consensus modeling,” Communications in Mathematical Sciences, vol. 13, no. 6, pp. 1407–1429, 2015. View at: Google Scholar
 G. Albi, L. Pareschi, and M. Zanella, “Boltzmanntype control of opinion consensus through leaders,” Philosophical Transactions of the Royal Society A, vol. 372, no. 2028, 2014. View at: Publisher Site  Google Scholar
 F. Cucker and S. Smale, “Emergent behavior in flocks,” IEEE Transactions on Automatic Control, vol. 52, no. 5, pp. 852–862, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 P. Degond, J.G. Liu, and C. Ringhofer, “Largescale dynamics of meanfield games driven by local Nash equilibria,” Journal of Nonlinear Science, vol. 24, no. 1, pp. 93–115, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 P. Degond, M. Herty, and J. G. Liu, “Meanfield games and model predictive control,” http://arxiv.org/abs/1412.7517. View at: Google Scholar
 S. Motsch and E. Tadmor, “A new model for selforganized dynamics and its flocking behavior,” Journal of Statistical Physics, vol. 144, no. 5, pp. 923–947, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 S. M. Ahn and S.Y. Ha, “Stochastic flocking dynamics of the CuckerSmale model with multiplicative white noises,” Journal of Mathematical Physics, vol. 51, no. 10, Article ID 103301, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 S.Y. Ha, K. Lee, and D. Levy, “Emergence of timeasymptotic flocking in a stochastic CuckerSmale system,” Communications in Mathematical Sciences, vol. 7, no. 2, pp. 453–469, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 Y. Sun, “Mean square consensus for uncertain multiagent systems with noises and delays,” Abstract and Applied Analysis, vol. 2012, Article ID 621060, 18 pages, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 C. A. Yates, R. Erban, C. Escudero et al., “Inherent noise can facilitate coherence in collective swarm motion,” Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 14, pp. 5464–5469, 2009. View at: Publisher Site  Google Scholar
 D. Xiu, Numerical Methods for Stochastic Computations, Princeton University Press, Princeton, NJ, USA, 2010. View at: MathSciNet
 M. Gerritsma, J.B. van der Steen, P. Vos, and G. Karniadakis, “Timedependent generalized polynomial chaos,” Journal of Computational Physics, vol. 229, no. 22, pp. 8333–8363, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 R. H. Cameron and W. T. Martin, “Transformations of Wiener integrals under translations,” Annals of Mathematics, vol. 45, no. 2, pp. 386–396, 1944. View at: Publisher Site  Google Scholar  MathSciNet
 D. Xiu and G. E. Karniadakis, “The WienerAskey polynomial chaos for stochastic differential equations,” SIAM Journal on Scientific Computing, vol. 24, no. 2, pp. 619–644, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 G. Albi and L. Pareschi, “Meanfield selective control of flocking dynamics,” In preparation. View at: Google Scholar
 H. Michalska and D. Q. Mayne, “Robust receding horizon control of constrained nonlinear systems,” IEEE Transactions on Automatic Control, vol. 38, no. 11, pp. 1623–1633, 1993. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2015 Giacomo Albi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.