Data-driven Discovery of Emergent Behaviors in Collective Dynamics

2019·arXiv

Abstract

1 Introduction

Emergent behavior in collective dynamics, such as clustering of opinions [2, 3, 4, 5], flocking of birds [6, 7, 8], milling of fish [9, 10, 11, 12], and concentric trajectories of planetary motion [13], is among one of the most interesting phenomena in macroscopic and microscopic scale systems. It occurs in systems used across many disciplines, including biology, social science, particle physics, astronomy, economics, and many more. Extensive studies have been conducted in order to understand the mechanism behind such intricate and yet geometrically simple behaviors. As shown in [14, 15, 7, 16, 17, 18, 5], these emergent behaviors are steady-states of various types of collective dynamics, and they can be qualitatively studied when the governing equations are known beforehand. However, if only the short-time trajectories of the dynamics are observed, it may be challenging to make accurate prediction about the emergent behaviors of the observed dynamics without prior knowledge of the governing equations. We offer a learning approach to overcome this difficulty by first discovering the governing equations from the observational data, and then use the estimated equations for large-time prediction.

Research on discovering governing equations of dynamical systems has enjoyed a long history in the science and engineering community; it can be traced back to the earlier work of Lagrange, Laplace and Gauss [19]. Among the many inspiring studies, the lengthy discovery of gravity had immense impact. In 1605, Kepler announced his first law of planetary motion, from his work on showing Mars’ elliptical orbit based on Tycho Brahe’s observational data. Based on Kepler’s first law and the assumption that gravity has a parametric form, namely , Newton formulated his law of universal gravitation, i.e., that gravity has the form 1, in 1687. Our learning approach can re-discover the 1form of the law of universal gravitation in a highly efficient and precise manner without the assumption of gravitation having a parametric form and planetary motion being elliptical, for details see Sec. 7.

Our work is focused on discovering governing equations of collective dynamics (also known as self-organized dynamics), a special kind of interacting particle- and agent-based dynamical systems. These are a distinguished subset of general autonomous systems of ODEs

In this general case, given (0 ), system identification consists in inferring f from X(t) and ˙X(t) observed at various t’s. Classical regression techniques (e.g. [20, 21, 22, 23]) have recently been brought to bear on this problem (e.g. [24, 25, 26, 27]). However, lack of independence among the observation data and the curse of dimensionality due to D typically being large, are all obstructions to finding the desired f effectively and efficiently (see [1] for an extended discussion) of these points. The works in [28, 1, 29] proposed a learning approach that exploits the special form of collective dynamics to overcome the difficulties mentioned above. It began with first order systems of the form

The 1-dimensional function is referred to as the interaction kernel. Here and in what follows we assume, with possibly abuse of notation, that the term = i in the sum in the r.h.s. is 0, even in cases where may not be defined at 0 (e.g. ) = 1). The aforementioned works consider the problem of estimating given trajectory observations, in terms of positions and velocities of the agents at various times, along one or multiple trajectories (with different initial conditions (ICs), e.g. sampled at random from some probability distribution on the state space). In [28, 1] a nonparametric learning approach to construct an estimator for is considered, that exploits the governing structure of the dynamics in (2), which is a special (yet ubiquitous) case of the general equations (1). The work [28] considered a first-order model of homogeneous agents (derived from gradient flow), and studied the convergence to its mean field limit and the inference of the mean-field limit interaction kernel from observations of trajectories of the system with a finite and yet increasing number of agents. The work [1] extended the approach in [28] to the situation where the number of agents is fixed, but the number of observations increases, showing that the nonparametric estimators for the interaction kernel converge at the near optimal rate for regression in one dimension, in particular independent of the dimension of the state space. It generalized the estimators to first and second-order of heterogeneous agents with 1-dimensional interaction kernels based on pairwise distances, providing substantial numerical evidence of the performance of these generalizations. The work [29] analyzes in detail the estimators for first order heterogeneous agent-based systems, generalizing the theoretical results of [1] to that case, while sharpening some of the constructions.

Here we extend the application of these approaches to rather general classes of agent-based systems, driven by first- and second-order dynamics, interaction kernels depending not only on pairwise distances but also on other pairwise data, which also depends on the states between the agents. We show that the estimators can be generalized to these settings, and measure their performance at both approximating the interaction kernel (in a suitable, dynamics-driven weighted -distance) and in predicting trajectories from the same initial conditions (ICs). These estimators may be constructed in a memory-efficient way, i.e. scalable to large data sets with a large number of agents. The estimated interaction kernels, inferred from short-time trajectory data, can provide very good approximations to the original unknown interaction kernels, and yield predictions from new initial conditions. Furthermore, the estimated interaction kernels can also provide insight into discovering the correct emergent behaviors at large time, as we will demonstrate in several examples in section 5 and section 6. We also extend these estimators to consider not only a single system, but a family of systems, governed by a family of interaction kernels . We consider the case of gravity in section 7, and show that we can discover both the “common structure” of the interaction kernel, namely the 1dependency on pairwise distance, and the dependency on mass.

The structure of our paper is as follows. In section 2 we discuss in detail the two models, of first and second-order systems respectively, which we are considering. In section 3 we outline the learning algorithm for each model. These algorithms are efficient and scalable. They enjoy highly favorable performance in terms of computational complexity as described in section 3.1. Section 4 defines the various performance measures, confusion matrices, and pattern indicator scores. It also discuses how we set up the numerical experiments. Section 5 to Section 7 provide a detailed study of five fundamental dynamical systems that vary across order, interaction kernel form, and agent characteristics, as well as learning of interaction kernels that involve parameters. Finally, we conclude the paper and discuss various future research directions stimulated by these results in section 8.

2 Model Description

The main focus of this work is to numerically investigate the capability of the estimators at predicting emergent behaviors of various collective dynamics using our extended learning approach. Our extended learning covers dynamical systems much more elaborate than those considered in [1]. These systems have a complex balance between non-collective and collective forces, include interaction laws depending on more than one variable, and allow for parametric families of interaction laws. We motivate these extensions with families of dynamical systems exhibiting these more intricate governing structures, which were motivated by various applications and whose dynamical properties have been studied in the scientific literature.

Here we consider particle- and agent-based systems that model rather general complex systems, beyond those considered in [28, 1, 29]. The first-order models are governed by the following system of coupled ODEs

These systems contain heterogeneous agents: the agents are partitioned into K different types, with containing the indices of the agents of type k, for k = 1, . . . , K. Table 1 shows the definitions of variables in (3).

Table 1: Notation for first-order models

Remark 2.1. Compared to Eqn. (8) in [1], our new equation (3) has the following additions: the new variable, non-collective forces , and 2-dimensional interaction laws (as opposed to only single-variable, pairwise-distance-based interactions).

for i = 1. Table 2 shows the definitions of variables in (4). Natural regularity and growth assumptions on the functions in the right-hand side are made so that the system has a unique solution for all times. For example assuming that the functions involved are at least Lipschitz and decay sufficiently rapidly at infinity would suffice.

Table 2: Notation for second-order models

Remark 2.2. Compared to Eqn. (11) in [1], our new equation (3) has the following additions: slightly different non-collective forces, , and 2-dimensional interaction laws.

We are given observation data, namely (((=

for second order systems) at time instances . In the case of missing derivative data, namely ˙, we will approximate ˙using appropriate finite difference schemes. The observation data is generated from M initial conditions (ICs), (0)), which are i.i.d samples from a (typically unknown) probability distribution (= for first order and = for second order). The unknowns in these systems are the interaction laws and the distribution of the initial conditions, while everything else is assumed known. We construct estimators for (resp. for second-order systems) that are close to the true interaction laws with high probability. Moreover, such estimators yield approximate systems, whose dynamics are approximations to the dynamics of the original system within the training time interval [], but can also provide approximations for emergent behaviors of collective dynamics, ranging from first-order opinion dynamics to second-order gravitational dynamics governing the planetary movement in our solar system. A key component of evaluating the emergent dynamics are appropriate measures of the presence of a specific emergent behavior, which will be discussed in 4.

3 Learning Algorithm

Similar to the the algorithm presented in [1], the learning algorithm which we use for the more complex dynamics considered here starts from the introduction of suitable cost functions whose minimizers, over a suitable approximation space, determine the estimators. Equation (3) can be rewritten in a more compact form: ) + )

Here ; for the interaction kernels, we use the vectorized notations, = and , and are the collection of the

corresponding right hand side terms in (3) respectively. Lastly, ) is defined as the vectorization of

the non-collective forces and ) is defined as the vectorization of the non-collective

where the norm is defined as

for with each or 1). Here is the same norm used in (3) and (4); = and are finite-dimensional hypothesis spaces. We choose each of the hypothesis space be to a finite dimensional function space of piece-wise polynomials of degree p, with p = 0 or 1 (polynomials of higher degree can be used and other type of basis functions are also possible, e.g., clamped B-splines, see [1]), with polynomial pieces supported on intervals that form a uniform partition of the observed range of variables []. Hence, each can be expressed in terms of the linear combination of the basis functions as follows

Similar definitions are used for each . Substituting this expression into the functionals above, the mini- mization becomes a set of linear equations,

Remark 3.1. In the case of missing ˙) (for first order system) or ¨) (for second order system), we will approximate it using an appropriate finite difference scheme. See Sec. 4.4 for details on how we setup the examples with or without derivative information.

In the case of the second-order dynamics described in (4), we introduce a new variable ) = ˙and let , a compact form of (4) is given as follows,

Here = . By choosing appropriate finite dimensional hypothesis spaces for and , e.g., piece-wise polynomials of degree p, we can simplify the least square problems down to the following linear systems which we solve to generate the necessary coefficients:

Here, with being the collection of ’s and being the collection of ’s.

Remark 3.2. For further details regarding the construction of the learning matrices, , and the right hand side vectors, , we refer the reader to Section 2 : Algorithm of the Supplementary Information of [1].

3.1 Computational Complexity

The learning approach, which is described in Sec. 3, can be easily parallelized in the m (number of initial conditions) variable. Although it takes MLD double-precision floating-point numbers (D = Nd + N for a first-order system, and D = 2Nd + N for a second-order system) to store the discrete trajectory data, each computing core j only needs to store floating-point numbers, with . Furthermore, each computing core does not need to hold all of the trajectory data in memory, since the assembly of the learning matrix and the right-hand-side vector needs only LD floating-point numbers (one system trajectory at a time). The sizes for the learning matrix and right hand side vector are: and 1 (or for a first-order system and + or for a second-order system), respectively. Since we have , which makes solving for our estimators extremely memory efficient. At each time instance, we have to compute the various pairwise variables, requiring ) distance calculations, hence the algorithm performs a total of ) computations of pairwise variables. In solving the linear system, it performs ) operations (or log(n)), we take the worst cases scenario since we use the built-in pseudoinverse routine in MATLAB to avoid any possible issues with numerical stability). The total computational complexity is + ). Online learning can be built into our learning approach: as trajectory data from different initial conditions comes in, one can simply average the estimators from previous trajectory data with the estimators from the new trajectory data to obtain a better approximation.

4 Performance Measures

We consider three different kinds of performance measures: how close the estimated interaction kernel(s) are to the true one(s), how well the trajectories of the system driven by the estimated interaction kernel(s) approximate the trajectories of the original system, and finally how well emergent patterns are reproduced/predicted in the system driven by the estimated interaction kernels. We use appropriate dynamics adapted measures, specified in section 4.1 and the Appendix.

4.1 Estimation error of interaction kernels

Following the definitions in [1], we introduce a set of probability measures to calculate the learning error between

(for continuous trajectory) and (for discrete trajectory) are only used in the theoretical setting; in practice, we use (with large M and L) for actual implementations and applications. These measures depend on the dynamical system and the distribution of initial conditions, weighting the areas of pairwise distances (the variable r) and of variables based on how often trajectories of the system explore them.

Table 3 explains the definitions of the variables in (5).

’s, Definition of the Variables

In the case of = 0, we define the corresponding ) to a zero function.

We measure the error of the interaction kernel estimators, , using the dynamics-induced weighted

However since is not calculable, we use instead. The weight, , comes from the governing structure of (3). Theoretical guarantees such as those in [1, 29] bound these errors, with high probability, as M grows. Extending those bounds to the general types of systems considered here will be investigated in future work. The results of our numerical experiments suggest that the learning rate, i.e. the rate of decrease of the error in (6), as a function of M is independent of the dimension of the state space of the system, and only depends crucially on the number of variables in the interaction kernel. The curse of dimensionality (of the state space) is therefore avoided.

4.2 Trajectory errors

We consider another performance measure, which might be estimated from data, especially when the true interaction kernel is not known, that quantifies the prediction capability of our estimators, by comparing the observed trajectories to the estimated trajectories evolved from the same initial conditions but using the estimated interaction laws. We will consider both X(t) =(and ) =for ]. Let , then the following norm is used

Here ˆis the estimated trajectory using our estimators with the same initial condition as in . The scaling by maxenables us to compare trajectory errors for different kinds of dynamics, especially those with large . Similarly,

For performance measures defined for and the second order systems, please see sec. A in the appendix.

4.3 Confusion Matrix and Pattern Indicator Scores

When a system is highly sensitive on small perturbations, or even chaotic, it is hopeless to expect that the estimated system will produce trajectories that are accurate approximations of the trajectories of the original system, except perhaps for very small times. However we have observed that certain large-time aspects of the dynamics of the system, such as certain emergent behavior including flocking or milling or clustering, are preserved in our estimated system, even when the trajectory-wise errors are relatively large. On the one hand, this may seem surprising, as at no point do we inject any knowledge about such emergent behaviors, into our estimator; on the other hand if such emergent behaviors are thought of as being “structurally robust” to perturbations of the system, and even perturbations of the laws of the system (the interaction kernels), then it becomes reasonable to expect that our estimated systems should preserve, at least to some degree, such emergent behaviors. We therefore introduce a way of measuring quantitatively the presence of such emergent behaviors, and quantify the performance in reproducing them in our estimated systems.

In order to accurately describe the capability of our estimators to predict the correct emergent behaviors at large time , we consider confusion matrices and “pattern indicator scores”. These are defined differently for each dynamical system to measure its unique emergent behavior.

Several aspects of the emergent behaviors that we are interested in are observables (i.e. functions defined on the state variables of the system). We define various emergent behavior scores, such as the flocking score, the milling score, etc., and choose a target range for the score to be in as an indicator of occurrence of the emergent behavior. For example, if the flocking score is within (0.99, 1], then flocking occurs. We calculate these scores on the true and estimated systems (systems with the same initial conditions as the true systems but evolved using the learned interaction law(s)). From this indicator of whether the emergent behavior occurred in the true/estimated system, we construct a confusion matrix, given as follows (in the case of learning flocking systems),

Table 4: General form of a confusion matrix. Each shows a probability (represented in percentages) of the combination, e.g., is the probability of the predicted system (evolved using the estimated interaction laws) showing non-flocking behavior given that the true system shows non-flocking behavior with the same initial conditions.

It is used to present the probability of the occurrence of the desired emergent behaviors in the true and estimated systems. Namely, if the true systems exhibit flocking with high probability, then the estimated systems should ideally show flocking with similar probability.

In order to provide deeper insight about the prediction of emergent behavior via confusion matrices, we also consider the following statistics from the confusion matrix.

Table 5: Definition of the statistical terms related to the confusion matrix.

Next we use a more refined measurement, a so-called ‘pattern indicator score’, to further demonstrate the capabilities of the estimated system at predicting emergent behaviors. Besides the emergent behavior scores, there are other quantitative descriptions of the emergent behaviors, such as the center-of-mass velocity in flocking, the common rotational axis in milling, the conservation of total energy in concentric trajectories, etc. The pattern indicator scores use these, sometimes together with the previously defined emergent behavior scores, to measure how well the estimated systems are predicting these observables compared to the true systems. Details of the definition of the confusion matrices and pattern indicator scores for each dynamics are in Sec. 5 to Sec. 7.

4.4 Setup of the Numerical Experiments

Here we describe the general setup for the subsequent sections of experiments. The various dynamical systems we consider exhibit a wide variety of emergent behaviors: clustering, flocking, milling, synchronization, and concentric trajectories. Different forms of interaction kernels are also considered, i.e., ) and ), where P is an unknown vector of parameters. These dynamics range from first-order dynamics of homogeneous agents to second-order dynamics of heterogeneous agents. We arrange the examples in three major sections based on the different types of the interaction laws.

The experiments are setup as follows: we first run different initial conditions generated i.i.d from the probability measure for initial condition, and evolvethe dynamics from 0 to T: the dynamics observed in [] is used to compute the probability measures ’s, which are empirical approximations to the prob- ability measures ’s. We do this only to compute and report the ) approximation errors; in practice this step is not required nor needed. Next, we generate another set of M random initial conditions and corresponding trajectories of the dynamics for [0, T], with each dynamics observed at L equidistant times , producing the observation data, i.e. , without the corresponding derivative information (i.e., ˙is not given, except for Synchronized Oscillator Dynamics and Gravitational Dynamics), as input to our estimation procedure. We construct the hypothesis spaces, where the estimators are found, on the learning intervals, e.g. [], derived from the observation data; the numbers of basis functions, as well as their degrees, are reported in each section. We report the ) errors between the estimated and true interaction kernels, as well as the trajectory errors based on the statistics over the training set and over a testing set (with new initial conditions), in the form of (mean value)(standard deviation). Then we consider the emergent behavior of the true dynamics and the predicted dynamics at , and evaluate “pattern indicator scores” and confusion matrices corresponding to the various kinds of emergent behaviors. The parameters used by all experiments are reported in table 6.

Table 6: Common Parameters

Each section/subsection is presented in a similar manner: we introduce the model and discuss why such model interesting for our learning approach; then, we relate the model equation to the learning paradigm presented in (3) and (4). Next, we present our learning results in figures and tables, in terms of approximation error, trajectory error, confusion matrix and pattern indicators. We end with a brief discussion of the learning results.

5 Emergent Behaviors Induced by φ(r)

We consider here three prototypical types of emergent behavior: clustering, flocking and milling. We examined four different examples of collective dynamics in order to comprehensively explore all three types of emergent behaviors, with very different dynamical behaviors.

5.1 Opinion Dynamics

The opinion dynamics (OD) model, first introduced in [2], is a prototypical first-order model of homogeneous agents which describes the interaction of people’s opinions through time, see details and extensions in [3, 4, 5, 30, 31, 32]. These models have gained popularity in modeling human’s social behavior, and they can be used to predict interesting social phenomena, namely, clustering/consensus of opinions.

The governing equations (being a vector of opinions) are:

Here 0 for all 0. With the interaction kernels giving attractive influences only, these models are bound to have clusters of opinions at large time. Table 7a shows how this dynamical system is mapped to the general form (3). The parameters used for setting up the experiment used are shown in table 7b.

Table 7: Opinion Dynamics (OD)

We consider the following interaction law,

Piece-wise constant polynomials with = 99 basis functions are used to approximate . The comparison of the true and the estimated ˆis shown in Fig.1.

Figure 1: (OD) Comparison of and ˆ, with the relative error being 11010(calculated using (6)). The true interaction kernel is shown in black solid line, whereas the mean estimated interaction kernel is shown in blue solid line with its confidence interval shown in blue dotted lines. Shown in the background is the comparison of approximated versus the empirical .

As it is shown in Fig. 1, not only can our estimator detect the discontinuity in the , but also can it detect the compact support of . Meanwhile, there is higher uncertainty in learning the interaction kernel at r = 0 (the information of (0) is lost since it is weighted by corresponding ) and at those discontinuity points. Since is non-negative, the agents in the system would eventually converge to clusters, this decreases the effective number of pairwise distance data for inferring . However, we are still able to provide an accurate estimator of by the continuity of the estimator. The comparison of a trajectory driven by the true versus the other one driven by the estimated ˆis shown in Fig. 2: there is no major visual difference between the true and predicted trajectories (generated from the training initial condition); the differences are quantified in table 8.

Figure 2: (OD) Comparison of X and ˆX, with the errors reported in table 8. The first row of trajectories are generated from an initial condition taken from the observation data. The second row of trajectories are generated from another randomly chosen initial condition. The first column of trajectories are generated from the true interaction kernel, whereas the second column of trajectories are generated from our estimated kernel with the same initial conditions. The color of the trajectory indicates the flow of time, from deep blue (at ) to light green (at ).

Table 8: (OD) Trajectory Errors: Initial Conditions (ICs) used in the training set (first two rows), new ICs randomly drawn from (second set of two rows). The trajectory errors is calculated using (7).

The confusion matrix and pattern indicator scores used to examine the capability of our estimators predicting the proper emergent behaviors associated with the Opinion Dynamics model are defined as follows. First, a confusion matrix is used to show the accuracy of our estimator to display the same clustering behavior as the true systems, see the results in table 9.

Table 9: (OD) Confusion Matrix: ICs used in the training set (first two rows), new ICs randomly drawn from (second set of two rows). It is generated using table 4.

We provide more statistics about the confusion matrix in order to understand our prediction of clustering better in table 10.

Table 10: (OD) Confusion Matrix Statistics: ICs used in the training set, new ICs randomly drawn from . It is generated using table 5.

Next, when the true system has clustering, we want to know if the predicted system can have the same number of clusters as the true system has. Hence, we assign a score of 1 when the predicted system shows the same number of clusters as the true systems; and a score of 0 when it predicts the wrong number of clusters. PIis the average of those scores over M trials. Lastly, we want to compare the clusters between the true and predicted systems. Let C contain the centers of the clusters at time from the true system, ˆC contain the centers of clusters from the estimated system; we shall use Hausdorff distance to calculate the distance between C and ˆC. PIis the average of M trials of such distances. See table 11 for details.

Table 11: (OD) Pattern Indicator Scores: ICs used in the training set (first two rows), new ICs randomly drawn from (second set of two rows).

The smaller PIis, the better we are at predicting the number of clusters right. Meanwhile, not only could we predict the number of clusters with high confidence, but also could we predict the actual location of the clusters.

5.2 Cucker-Smale Dynamics

Modeling how animals (or other living agents) move in a cohesive group formation has been a challenging and well-studied problem [6, 33, 34], [9, 35, 12], [36, 37, 15]. There are different degress of cohesion in a collective system: flocking (where each agent shares a common velocity), milling (where each agent rotates around the same axis or about the same oint), and swarming (transition state between flocking and milling). We first consider the simplest cohesion in a collective system, namely flocking (see detailed work in [38, 39, 8, 40, 41, 34, 42], its mean field limit in [43, 44, 45], and extension to a stochastic system in [46] and references therein), and investigate the learnability of these flocking systems.

The Cucker-Smale (CS) dynamics is one of the prototypical examples of flocking agents. Its governing equations are

Here ) = (1 + where are chosen parameters. Table 12 shows how this dynamical system is mapped to the general form (4).

Table 13: (CS) Parameters for Experiment Setup

With certain choices of H and , the Cucker-Smale system is guaranteed to produce flocking (where all agents have the same final velocity) see [38]. For example, when , the system is guaranteed to have flocking regardless of initial conditions; when , the system has conditional flocking depending on the initial configuration of velocities; when , the system has conditional flocking depending on the initial configuration of both positions and velocities.

We consider the following interaction law,

With this interaction kernel, the agents are guaranteed to flock (see theorem 2, 3 in [38]). We use the following parameters in table 13 to set up the experiment. Piece-wise linear polynomials with = 100 basis functions are used to approximate . The comparison of the true and the estimated ˆis shown in Fig.3.

Figure 3: (CS) Comparison of and ˆtogether with a plot of versus , with the relative error being 41010(calculated using (19)). The true interaction kernel is shown by in a black solid line, whereas the mean estimated interaction kernel is shown in a blue solid line with its confidence interval shown in blue dotted lines. Shown in the background is the comparison of approximated versus the empirical .

As it is shown in Fig. 3, our learning approach produce faithful approximation to , especially capturing the tail behavior of the original interaction law, notwithstanding the scarcity of samples in that region of pairwise distances and speeds; towards r = 0, the estimated kernel is also close to the true kernel. The comparison of true trajectory X(t) and learned ˆX(t) is shown in Fig. 4.

Figure 4: (CS) Comparison of X and ˆX, with the errors reported in table 14. The first row of trajectories are generated from an initial condition taken from the observation data. The second row of trajectories are generated from another randomly chosen initial condition. The first column of trajectories are generated from the true interaction kernel, whereas the second column of trajectories are generated from our estimated kernel with the same initial conditions. The color of the trajectory indicates the flow of time, from deep blue (at ) to light green (at ).

Fig. 4 shows no visual difference between the true trajectories and the learned trajectories (for the training initial condition and a randomly chosen initial condition), we provide a quantitative insight into the difference between trajectories in table 14.

Table 14: (CS) Trajectory Errors: Initial Conditions (ICs) used in the training set (first two rows), new ICs randomly drawn from (second set of two rows). The trajectory errors in x/v is calculated using (7)/(8).

We consider the Flocking score (at ) taken from [6],

When = 0, perfect flocking occurs; however we use 1 to indicate flocking. The prediction capability of the estimated systems in the form of confusion matrix is reported in table 15.

Table 15: (CS) Confusion Matrix: ICs used in the training set (first two rows), new ICs randomly drawn from (second set of two rows). It is generated using table 4.

With = , the true system is guaranteed to show flocking. Since we have no control over when the flocking would occur, we use to help us capture the essence of flocking behavior, i.e. = 0 (or close to 0). And our estimated system can capture the same behavior (flocking or not) with high probability. We provide more statistics about the confusion matrix in order to understand our prediction of clustering better in table 16.

Table 16: (CS) Confusion Matrix Statistics: ICs used in the training set, new ICs randomly drawn from . It is generated using table 5.

In order for us to provide more quantitative insight into the predication capability of our estimator for the case of flocking, we consider two different pattern indicator scores. First, PIis the relative error of between true and predicted systems, averaged over M trials. Second, we consider another important quantity, the center of mass velocity, . It is given by =

for i = 1= ). Then we define PIto be the relative error of the predicted center of mass velocity and true center of mass velocity averaged over M trials. The scores are reported in table 17.

Table 17: (CS) Pattern Indicator Scores: ICs used in the training set (first two rows), new ICs randomly drawn from (second set of two rows).

As it is shown in table 17, our estimated system can predict with relatively high accuracy. Surprisingly, our estimated system can reproduce down to numerical accuracy.

5.3 Fish Milling in 2 dimensions

Next we consider a more complicated cohesive collective system: a dynamical system which produces milling patterns, where each agent rotates around the same axis or about the same point. The models we consider have been proposed in [9, 12] (see references therein for a variety of sources for the biological roots of these models). Useful background references for the two-dimensional models are [47, 10] as well as the primer [48]. Further theoretical study of models of this type has been done in [49, 50, 11].

The governing equations of the Fish Milling Dynamics in (FM2D) of [9] are,

Here, is the Morse potential describing the interaction of the agent with the other agents in the system, defined as follows

Here are the attraction/repulsion strengths and are the effective attraction/repulsion lengths. Table 18 shows how the FM2D dynamics fits into the framework of (4).

Table 19: (FM2D) Parameters for Experiment Setup

The delicate balance between the self-propelling force produced by ) and the collective force induced by the energy kernel can create a wide range of patterns for different initial conditions. Milling patterns (single or double milling) are one of the most interesting ones. Unlike the well-understood Cucker-Smale model, necessary and sufficient conditions on the interaction kernels and ICs that guarantee the existence milling patterns seem to be unknown. These milling patterns result from the balance of the non-collective force and the collective force induced by the energy kernel (especially when is not H-stable, double-milling would occur, see Fig. 1 in [9]), and are therefore rather sensitive to the selection of parameters: relatively small differences in the interaction laws can correspond to dynamical systems with very different dynamical patterns. The estimator error between the true and estimated interaction kernel may therefore offer little insight information into how well our estimated dynamics can re-produce milling patterns at large time. The confusion matrix and pattern indicator scores are finer indicators of performance in this case.

We consider the following interaction law,

With this setup, a double-milling pattern appears 100% of the time (see [12]). The other parameters are reported in table 19. We use piecewise constant polynomials with = 122 basis functions to approximate . The comparison of the true and the estimated ˆis shown in Fig.5.

Figure 5: (FM2D) Comparison of and ˆ, with the relative error being 6(calculated using (6)). The true interaction kernel is shown in black solid line, whereas the mean estimated interaction kernel is shown in blue solid line with its confidence interval shown in blue dotted lines. Shown in the background is the comparison of approximated versus the empirical .

As it is shown in Fig. 5, our estimator closely resembles , however when r is close to 0, there is a sharp drop of to , the availability of r data close to 0 becomes scarcer, and since we are using a uniform basis to approximate , the difference between and ˆis apparent in this range. The comparison of the true trajectory X(t) and learned ˆX(t) is shown in Fig. 6.

Figure 6: (FM2D) Comparison of X and ˆX, with the errors reported in table 20. The first row of trajectories are generated from an initial condition taken from the observation data. The second row of trajectories are generated from another randomly chosen initial condition. The first column of trajectories are generated from the true interaction kernel, whereas the second column of trajectories are generated from our estimated kernel with the same initial conditions. The color of the trajectory indicates the flow of time, from deep blue (at ) to light green (at ).

Our predicted system can still estimate the position/velocity of the agents in large time, i.e., for , with relatively small error, around 10. Moreover, when the dynamics enters its milling state, our predicted system is also in the same milling state. We provide a quantitative insight into the difference between trajectories in table 20.

Table 20: (FM2D) Trajectory Errors: Initial Conditions (ICs) used in the training set (first two rows), new ICs randomly drawn from (second set of two rows). The trajectory errors in x/v is calculated using (7)/(8).

We are getting 10relative accuracy for estimating positions/velocities within the learning time interval (i.e. [0, T]); however, as time goes on, we can see roughly linear growth of the errors (= 5T).

We consider the center of mass position ) =

) = ). We consider the indicator score (at , where are taken from [9]. Here, the flocking score is defined as,

Again = 1 when perfect flocking occurs. Then we consider the the milling score as follows,

When = 1 when perfect milling(around the same axis) occurs; meanwhile = 0 if = 1. Therefore 1]. As suggested by the thresholds in [12], we use the case, 5, to indicate milling. The true systems always show milling pattern (in fact, it shows double milling), and 100% of our estimated systems also show milling.

Next for the pattern indicator scores, let PIbe the relative error of over M trials. And PIis the relative error between the pair () (in norm) over M trials. The scores are reported in table 21.

Table 21: (FM2D) Pattern Indicator Scores: ICs used in the training set (first two rows), new ICs randomly drawn from (second set of two rows).

Milling patterns in dynamics are very delicate. The intricate balance and the H-stability of decides the appearance of milling in a dynamics, especially when is not H-stable for double milling patterns. In the case of the true dynamics showing milling (to be exact, double milling), our predicted systems can capture the same behavior (with high accuracy) both in terms of and the pair, ().

5.4 Fish Milling in 3 dimensions

Next, we consider a cohesive collective dynamics in 3D of self-propelled particles within a fluid environment, introduced in [12]. It is a more complicated 3D extension of the FM2D model, where agents could experience self-propelling force in a fluid.

The governing equations of the Fish Milling Dynamics in (FM3D) are,

Here, u is the lab-frame fluid velocity generated at position ( ˙)) gives the drag force (0), ( ˙)) represents the self-propelling motility force, and is the agent-to-agent interaction force on agent i, and the energy potential is the same Morse potential defined in sec. 5.3. is defined as follows