On Cucker-Smale model with noise and delayThe research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ ERC grant agreement No. 239870.

# On Cucker-Smale model with noise and delay††thanks: The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (Fp7/2007-2013)/ Erc grant agreement No. 239870.

Radek Erban,    Jan Haškovec   and   Yongzheng Sun
###### Abstract

A generalization of the Cucker-Smale model for collective animal behaviour is investigated. The model is formulated as a system of delayed stochastic differential equations. It incorporates two additional processes which are present in animal decision making, but are often neglected in modelling: (i) stochasticity (imperfections) of individual behaviour; and (ii) delayed responses of individuals to signals in their environment. Sufficient conditions for flocking for the generalized Cucker-Smale model are derived by using a suitable Lyapunov functional. As a byproduct, a new result regarding the asymptotic behaviour of delayed geometric Brownian motion is obtained. In the second part of the paper results of systematic numerical simulations are presented. They not only illustrate the analytical results, but hint at a somehow surprising behaviour of the system - namely, that an introduction of intermediate time delay may facilitate flocking.

11footnotetext: Mathematical Institute, University of Oxford, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, United Kingdom; e-mail: erban@maths.ox.ac.uk. Radek Erban would like to thank the Royal Society for a University Research Fellowship and the Leverhulme Trust for a Philip Leverhulme Prize.22footnotetext: Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; e-mail: jan.haskovec@kaust.edu.sa.33footnotetext: School of Sciences, China University of Mining and Technology, Xuzhou 221116, PR China; e-mail: yzsung@gmail.com. Yongzheng Sun acknowledges financial support from the China Scholarship Council (Grant No. 201308320087) and the National Natural Science Foundation of China (Grant No. 61403393).

Keywords: Cucker-Smale system, flocking, asymptotic behaviour, noise, delay, geometric Brownian motion.

## 1 Introduction

Collective coordinated motion of autonomous self-propelled agents with self-organization into robust patterns appears in many applications ranging from animal herding to the emergence of common languages in primitive societies [26]. Apart from its biological and evolutionary relevance, collective phenomena play a prominent role in many other scientific disciplines, such as robotics, control theory, economics and social sciences [5, 31, 23]. In this paper we study the interplay of noise and delay on collective behaviour. We investigate a modification of the well known Cucker-Smale model [6, 7] with multiplicative noise and reaction delays.

We consider autonomous agents described by their phase-space coordinates , , where (resp. ) are time-dependent position (resp. velocity) vectors of the -th agent. The governing equations are given as the following system of delayed Itô stochastic differential equations

 dxi = vidt, (1) dvi = λNN∑j=1ψij(˜vj−˜vi)dt+σiNN∑j=1ψij(˜vj−˜vi)dBti, (2)

where the delayed velocity is given by and is the reaction delay. The parameters and , measure the alignment and noise strength, respectively, and , , are independent -dimensional white noise vectors. In general, the communication rates are functions of the mutual distances , however, in most of our paper we will consider them as given functions of time satisfying certain assumptions. The standard Cucker-Smale model [6, 7] is a special case of equations (1)–(2) for and . Our aim is to investigate equations (1)–(2) for general values of reaction delay and noise strength parameters ,

The Cucker-Smale model was introduced and studied in the seminal papers [6, 7], originally as a model for language evolution. Later the interpretation as a model for flocking in animals (birds) prevailed. In general, the term flocking refers to the phenomena where autonomous agents reach a consensus based on limited environmental information and simple rules. The Cucker-Smale model is a simple relaxation-type model that reveals a phase transition depending on the intensity of communication between agents. Using and in (1)–(2), we can write the standard Cucker-Smale model as the following system of ordinary differential equations

 ˙xi = vi, (3) ˙vi = λNN∑j=1ψij(vj−vi),for i=1,2,…,N, (4)

where the dots denote the time derivatives. We note that the scaling by in (4) is significant to obtain a Vlasov-type kinetic equation in the mean-field limit , see, for example [28]. The communication rates introduced in [6, 7] and most of the subsequent papers are of the form

 ψij=ψ(|xi−xj|)with ψ(s)=1(1+s2)β,whereβ≥0. (5)

If , then the model exhibits the so-called unconditional flocking, where for every initial configuration the velocities converge to the common consensus value as . On the other hand, with the flocking is conditional, i.e., the asymptotic behaviour of the system depends on the value of and on the initial configuration. This result was first proved in [6, 7] using tools from graph theory (spectral properties of graph Laplacian), and slightly later reproved in [28] by means of elementary calculus. Another proof has been provided in [14], based on bounding by a system of dissipative differential inequalities, and, finally, the proof of [4] is based on bounding the maximal velocity.

Various modifications of the generic model have been considered. For instance, the case of singular communication rates was studied in [14, 24]. Motsch and Tadmor [21] scaled the communication rate between the agents in terms of their relative distance, so that their model does not involve any explicit dependence on the number of agents. The dependence of the communication rate on the topological rather than metric distance between agents was introduced in [15]. The influence of additive noise in individual velocity measurements was studied in [13] and [29]. Stochastic flocking dynamics with multiplicative white noises was considered in [1]. Delays in information processing were considered in [18], however, their analysis only applies to the Motsch-Tadmor variant of the model.

In this paper, we are interested in studying the combined influence of noise and delays on the asymptotic behaviour of the Cucker-Smale system. In particular, we derive a sufficient condition in terms of noise intensities and delay length that guarantees flocking. Our analysis is based on a construction of a Lyapunov functional and an estimate of its decay rate. To prove our main results, we make an additional structural assumption about the matrix of communication rates which, loosely speaking, means that the communication between agents is strong enough.

The paper is organized as follows: In Section 2 we show how our model is derived from the Cucker-Smale model and define what is meant by flocking. Moreover, we consider a simplified version of the model to provide an intuitive understanding of what qualitative properties may be expected. In Section 3 we derive a sufficient condition for flocking in terms of the parameters , and , based on a micro-macro decomposition and construction of a Lyapunov functional. Moreover, as a byproduct of our analysis, we provide a new result about the asymptotic behaviour of delayed geometric Brownian motion. Section 4 is devoted to a systematic numerical study of the model. First, we focus on simulation of delayed geometric Brownian motion, in particular, we study the dependence of its asymptotic behaviour on the delay and noise levels. Then, we perform the same study for system (1)–(2). This leads to the interesting observation that, for weak coupling and small noise levels, an introduction of intermediate delays may facilitate flocking. A systematic study of this effect concludes the paper.

## 2 The stochastic Cucker-Smale model with delay

In order to make the generic model (3)–(4) more realistic, we amend it with two additional features. First, we note that measurements in the real world are subject to errors and imprecisions that are typically modeled in terms of white noise. In particular, we assume that the state (velocity) of agent measured by agent is given by the expression

 ωi;j=vj+κi(vj−vi)dBti,% fori,j=1,2,…,N, (6)

where represents the imprecision of ’s measurement device, and are independent identically distributed -dimensional Brownian motions with zero mean and the covariance relations

 E[Bti⋅Bsj]=dδijδtst,fori,j=1,2,…,N and t,s≥0,

with iff and otherwise, and similarly for .

Note that the multiplicative structure of the noise term ensures that . Substituting given by (6) for in (4) and defining , we obtain the following system of stochastic differential equations (SDEs) for velocities

 dvi=λNN∑j=1ψij(vj−vi)dt+σiNN∑j=1ψij(vj−vi)dBti,fori=1,2,…,N,

with a positive coupling strength.

The second amendment is the introduction of delays, motivated by the fact that agents react to information received from their surroundings with some time lag. However, we assume that information propagates instantaneously, so the delay does not depend on the physical distance between agents. For simplicity, we assume the reaction lag to be the same for all agents, so that at time they react to information perceived at time for a fixed .

###### Convention 1

Throughout the paper, we denote by the quantity evaluated at time , i.e., , and by the same quantity evaluated at time , i.e., . We will also write resp. for the vectors of locations resp. velocities of the agents.

In general, the communication rates may be functions of the mutual distances . However, our analysis is based on a certain structural assumption about the communication matrix and the particular form of the dependence on the mutual distances is irrelevant. Therefore, we consider the rates as given adapted stochastic processes, so that decouples from . Moreover, we assume that are uniformly bounded,

 0<ψij(t)=ψji(t)≤1,for t≥0,i,j=1,2,…,N,almost surely. (7)

Thus, we finally arrive at the stochastic system of delayed differential equations that we will study,

 dvi = λNN∑j=1˜ψij(˜vj−˜vi)dt+σiNN∑j=1˜ψij(˜vj−˜vi)dBti, (8)

which is supplemented with the deterministic constant initial datum ,

 vi(t)≡v0i, for t∈(−τ,0],i=1,2,…,N. (9)

Let us note that we interpret the noise term in in terms of the Itô calculus [22, 19].

###### Theorem 1

The stochastic delay differential system  with initial condition  admits a unique global solution on which is an adapted process with for all , i.e., a martingale.

Proof: The proof follows directly from Theorem 3.1 and the subsequent remark on p. 157 of [19]. Indeed, (8) is of the form

 dv(t)=F(t,v(t−τ))dt+G(t,v(t−τ))dBt

for suitable functions and . In particular, the right-hand side is independent of the present state , so that the solution can be constructed by the method of steps. The second order moment is bounded on because of the linear growth of the right-hand side of (8) in .

We now define the property of asymptotic flocking for the solutions of (8)–(9).

###### Definition 1

We say that the system exhibits asymptotic flocking if the solution for any initial condition satisfies

 limt→∞∣∣E[vi(t)]−E[vj(t)]∣∣=0{\rm for all} i,j=1,2,…,N,

where denotes the expected value of a stochastic process.

### 2.1 Simplified case with ψ≡1

To get an intuition of what qualitative properties we may expect from the solutions of , we consider the case when the communication rate is constant, i.e., ; in other words, we set in (5). We also assume that is equal to the same constant for all , i.e. , and, moreover, that for some and all . Then, defining , we obtain

 dVc=σNN∑i=1(˜Vc−˜vi)dBti.

Since, by assumption, for , we have for all . Consequently, decouples into copies of the delayed SDE

 dw=−λ˜wdt−σ˜wdBt, (10)

where we denoted for any . We are not aware of any results concerning the asymptotic behaviour of . The method developed in [3] suggests that

 limt→∞E[|w(t)|2]=0if and % only if∫∞0rλ(t)2dt<1σ2,

where is the fundamental solution of the delayed ODE

 ˙w=−λ˜w, (11)

i.e., formally, solves (11) subject to the initial condition for . The fundamental solution can be constructed by the method of steps [25], however, evaluation of its -norm is an open problem. From this point of view, the analysis carried out in Section 3 provides new and valuable information about the asymptotics of (10), see Section 3.4. Let us note that setting in the above criterion recovers the well-known result about geometric Brownian motion [22]: the mean square fluctuation tends to zero if and only if .

Finally, for the convenience of the reader, we give an overview of the qualitative behaviour of solutions to (11) with , subject to a constant nonzero initial datum (see, e.g., Chapter 2 of [25]):

• If , the solution monotonically converges to zero as , hence no oscillations occur.

• If , oscillations appear, however, with asymptotically vanishing amplitude.

• If , periodic solutions exist.

• If , the amplitude of the oscillations diverges as .

Hence, we conclude that the (over)simplified model , corresponding to the delayed Cucker-Smale system with and no noise, exhibits flocking if and only if . In the next Section we derive a sufficient condition for flocking for the general model .

## 3 Sufficient condition for flocking

In this section we derive a sufficient condition for flocking in according to Definition 1. Our analysis will be based on a construction of a Lyapunov functional that will imply decay of velocity fluctuations for suitable parameter values. However, we will have to adopt an additional structural assumption on the matrix of communication rates .

Before we proceed, let us shortly point out the mathematical difficulties that arise due to the introduction of delay and noise into the Cucker-Smale system. The “traditional” proofs of flocking of model (3)–(4), for instance [6, 7, 28, 14], rely on the monotone decay of the kinetic energy (velocity fluctations) of the form

 ddtN∑i=1|vi|2=−λNN∑i=1N∑j=1ψij|vi−vj|2≤0.

However, this approach fails if processing delays are introduced, since for without noise (i.e., all ), we have

 ddtN∑i=1|vi|2=−λNN∑i=1N∑j=1˜ψij(vi−vj)⋅(˜vi−˜vj).

One then expects the product to be nonnegative for small enough, however, it is not clear how to prove this hypothesis.

The introduction of noise leads to additional difficulties - in particular, the classical bootstrapping argument [6, 7, 14] for fluctuations in velocity fails in this case. Similarly as in [13], we circumvent this problem by adopting, in addition to the boundedness (7), a structural assumption about the matrix of communication rates. We define the Laplacian matrix by

 Aij:=−ψij, for i≠j,Aii:=∑j≠iψij, (12)

and note that is symmetric, diagonally dominant with non-negative diagonal entries, thus it is positive semidefinite and has real nonnegative eigenvalues. Due to its Laplacian structure, its smallest eigenvalue is zero [6]. Let us denote its second smallest eigenvalue (the Fiedler number) . Our structural assumption is that there exists an such that

 μ2(t)≥ℓ>0,for t>0,almost surely. (13)

This can be guaranteed for instance by assuming that the communication rates are uniformly bounded away from zero, , since there exists a constant such that , see Proposition 2 in [6].

Moreover, we assume that the matrix of communication rates is uniformly Lipschitz continuous in the Frobenius norm, in particular, there exists a constant such that

 ∥A(t)−A(t−τ)∥F≤Lτfor t≥0, (14)

where denotes the Frobenius matrix norm.

To ease the notation and without loss of generality, we will consider the one-dimensional setting , i.e., and , where is the number of agents. Then, with the definition (12), we put (8) into the form

 dvi=−λN(˜A˜v)idt+σiN(˜A˜v)idBti,i=1,2,…,N. (15)

Our main result is the following.

###### Theorem 2

Let be given by satisfying , and . Let the parameters and satisfy

 σ2max<λ, (16)

then there exists a critical delay , independent of , such that for every the system  exhibits flocking in the sense of Definition 1.

Moreover, if the matrix of communication rates is constant, i.e. holds with , then is of the form

 τc=1λ2(−σ2max+√σ4max+112(λ−σ2max)2). (17)
###### Remark 1

The system with constant communication matrix can be seen as a linearization of the system about the equilibrium for with some . Note that in this case the formula  for the critical delay does not depend on the particular value of in .

### 3.1 Micro-macro decomposition

We introduce a micro-macro decomposition [28, 13] which splits into two parts: macroscopic, that describes the coarse-scale dynamics, and microscopic, that describes the fine-scale dynamics. The macroscopic part for the solution is set to be the mean velocity ,

 Vc(t):=1NN∑i=1vi(t). (18)

The microscopic variables are then taken as the fluctuations around their mean values,

 wi(t):=vi(t)−Vc(t),fori=1,2,…,N. (19)

We denote . Then we have

 w(t)=v(t)−Vc(t)e,wheree:=(1,1,…,1)T∈RN. (20)

Since is the eigenvector of corresponding to the zero eigenvalue, we have . Then (15) can be rewritten as follows

 dwi=−λN(˜A˜w)idt+σiN(˜A˜w)idBti−dVc. (21)

The macroscopic variable satisfies the following lemma.

###### Lemma 1

Let be a solution of  subject to the deterministic constant initial datum . Then for and for all .

Proof: The boundedness of follows directly from the definition (18) and the martingale property of provided by Theorem 1. Using (12), we have

 N∑i=1(˜A˜v)i=eT˜A˜v=0.

Summing equations (15), , using (18) and , we obtain that the macroscopic dynamics is governed by the system

 dVc=1N2N∑i=1σi(˜A˜w)idBti. (22)

After integration in time this implies

Since is a martingale, we have (see Theorem 5.8 on p. 22 of [19]). Thus we obtain .

###### Remark 2

Note that and are expressed in terms of the -variables only and so they form a closed system, which is equivalent to .

Clearly, due to (19), we have . Consequently, it is natural to introduce the decomposition , where is given by (20). We then have for all .

###### Lemma 2

Let , , be the matrix defined in and assume that and hold. Then:

(a) The maximal eigenvalue of is bounded by .

(b) We have for any vector

(c) We have for any vector .

(d) For any vectors , and we have

Proof:

1. The claim follows from the Gershgorin circle theorem. Indeed, since , the diagonal entries satisfy , and for all .

2. The smallest eigenvalue of is zero with the corresponding eigenvector . The second smallest eigenvalue (Fiedler number) is assumed to be positive by (13). Thus, is a symmetric, positive operator on the space and there exists an orthonormal basis of composed of eigenvectors of corresponding to the positive eigenvalues . Then, every vector can be decomposed as

 u=(u⋅e)e|e|2+N∑i=2(u⋅ξi)ξi. (23)

Thus, due to the above bound on the eigenvalues , we have

 |Au|2=N∑i=2(u⋅ξi)2μ2i≤2(N−1)N∑i=2(u⋅ξi)2μi=2(N−1)uTAu.
3. If , then (23) implies

 w=N∑i=2(w⋅ξi)ξi,andwTAw=N∑i=2(w⋅ξi)2μi.

Since nonzero eigenvalues are bounded from below by (using (13)) and from above by part (a) of this lemma, we obtain

 ℓ|w|2=ℓN∑i=2(w⋅ξi)2≤N∑i=2(w⋅ξi)2μi≤2(N−1)N∑i=2(w⋅ξi)2=2(N−1)|w|2.
4. With the orthonormality of the basis and the positivity of the eigenvalues , we have by the Cauchy-Schwartz inequality

 uTAw = (N∑i=2(u⋅ξi)ξi)TA(N∑i=2(w⋅ξi)ξi)=N∑i=2(u⋅ξi)(w⋅ξi)μi ≤ (N∑i=2(u⋅ξi)2μi)12(N∑i=2(w⋅ξi)2μi)12=(uTAu)12(wTAw)12,

and with any ,

 (uTAu)12(wTAw)12≤12δuTAu+δ2wTAw.

### 3.2 Lyapunov functional

The proof of Theorem 2 relies on estimating the decay rate of the following Lyapunov functional for (21)–(22),

 L(t):=|w(t)|2 + qN2∫tt−τ|A(s)w(s)|2ds + pN2∫tt−τ∫tθ|A(s−τ)w(s−τ)|2ds,

where , are positive constants depending on , and .

###### Lemma 3

Let the assumptions of Theorem 2 be satisfied. Then there exist positive constants , and such that for every solution of the Lyapunov functional satisfies

 ddtE[L(t)]≤−εNE[wT˜Aw]. (25)

Proof: We apply the Itô formula to calculate . Note that the Itô formula holds in its usual form also for systems of delayed SDE, see page 32 in [12] and [17, 8, 20]. Therefore, we obtain

With the identity (formula (6.11) on p. 36 of [19]), we have

 1N4N∑j=1N∑k=1σjσk(˜A˜w)j(˜A˜w)kdBtjdBtk=1N4N∑j=1σ2j(˜A˜w)2jdt.

Consequently, summing over , using and the identity

 N∑i=1(σiN(˜A˜w)idBti)dVc=1N3N∑i=1σ2i(˜A˜w)2idt,

we obtain

Consequently, we have

 dL(t) = (−2λNwT˜A˜w+N−1N3N∑i=1σ2i(˜A˜w)2i)dt+2NN∑i=1σiwi(˜A˜w)idBti (26) +qN2(|Aw|2−|˜A˜w|2)dt+pN2(τ|˜A˜w|2−∫tt−τ|Aw(s−τ)|2ds)dt.

Our goal is to estimate from above. First of all, we note that by the elementary property of the Itô integral (Theorem 5.8 on p. 22 of [19]),

 E[1NN∑i=1σiwi(˜A˜w)idBti]=0.

For the first term of the right-hand side in , we write

 −2λNwT˜A˜w=−2λNwT˜Aw+2λNwT˜A(w−˜w)

and apply Lemma 2(d) with ,

 2λNwT˜A(w−˜w)≤λδ−1NwT˜Aw+λδN(w−˜w)T˜A(w−˜w).

Using Lemma 2(c), we have

 λδN(w−˜w)T˜A(w−˜w)≤2λδN−1N|w−˜w|2. (27)

Now we write for , componentwise, using ,

Thus, we have for the expectation of the square

An application of the Cauchy-Schwartz inequality and Fubini’s theorem for the first term of the right-hand side yields

 3E[λN∫tt−τ(A(s−τ)w(s−τ))ids]2≤3λ2N2τ∫tt−τE[|(A(s−τ)w(s−τ))i|2]ds.

For the second term we use the fundamental property of the Itô integral (Theorem 5.8 on p. 22 of [19]),

 3E[σiN∫tt−τ(A(s−τ)w(s−τ))idBsi]2=3σ2iN2∫tt−τE[|(A(s−τ)w(s−τ))i|2]ds.

Similarly, the third term is estimated as

Thus, we get from , estimating ,

 λδNE[(w−˜w)T˜A(w−˜w)]≤6λδN2(λ2τ+2σ2max)N∑i=1∫tt−τE[|(A(s−τ)w(s−τ))i|2]ds.

An application of Lemma 2(b) gives

 qN2E[|Aw|2]≤2q(N−1)N2E[wTAw]≤2qNE[wTAw].

To balance this term with , we use assumption and Lemma 2(c) in

Collecting all the terms in finally leads to

We set

 p=6λδ(λ2τ+2σ2max),q=σ2max+pτ, (28)

then the above expression simplifies to

 ddtE[L(t)]≤1N[−2λ+λδ−1+2q(Lℓ−1τ+1)]E[wT˜Aw]. (29)

We want . Substituting (28) into this inequality leads to a third order polynomial inequality in . This polynomial has all positive coefficients but the zero order one, which is . If is satisfied, then choosing such that

 δ−1<2λ−1(λ−σ2max)

makes negative. Consequently, there exists a such that for any ,

 −2λ+λδ−1+2q(Lℓ−1τ+1)=−ε<0.

This completes the proof of .

It remains to study the case when is a constant matrix, i.e. in . Then simplifies to

 ddtE[L(t)]≤1N(−2λ+λδ−1+2q)E[wT˜Aw]

and we have to find such that . Again, substituting (28) for and leads to

 τ<λ(2−δ−1)−2σ2max12λδ(λ2τ+2σ2max). (30)

The maximum value of the right hand side is obtained for which is positive because of the first inequality in (17). Substituting into (30), we obtain

 τ<(λ−σ2max)212λ2(λ2τ+σ2max).

 τ<1λ2(−σ2max+√σ4max+112(λ−σ2max)2).

If the above sharp inequality is satisfied, there exists an such that and, consequently, (25) holds.

### 3.3 Proof of Theorem 2

An integration of in time gives

An application of Lemma 2(c) gives then

 E[|w|2(t)]≤E[L(0)]−εℓN∫t0E[|w(s)|2]ds,

so that the last integral is convergent as and, consequently, for all