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.

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 [47] and their kinetic and mean-field approximations [813]. 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, 1417]. 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 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 with a given distribution . Let ,  , 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 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 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: . It is straightforward to rewrite as a translation of a standard normal-distributed variable ; that is, where . The expected velocities read which correspond toSimilarly 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 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 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 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 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 ; 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 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 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 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 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 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 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 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 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 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 , a bounded control having value in an admissible set ; for example, in the one-dimensional 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 real-valued 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 time-continuous 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 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 , 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 , 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 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 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 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 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

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 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.

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.

Observe that the numerical results are coherent with the explicit solution of the controlled equation. Let us consider the time-independent 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 time-dependent scattering rate . The dynamic shows how, for the choice of time-dependent 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.

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 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.

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.

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 ERC-Starting Grant High-Dimensional Sparse Optimal Control (HDSPCONTR-306274).