Evolution of clusters in large-scale dynamical networks

Recent tremendous progress in electronics, complexity theory and network science provides new opportunities for intellectual control of complex large-scale systems operating in turbulent environment via networks of interconnected miniature devices, serving as actuators, sensors and data processors. Actual dynamics of the resulting control systems are too sophisticated to be examined controlled by traditional methods, which primarily deal with ordinary differential equations. However, their complexity can be dramatically reduced by fast processes, organizing the elementary units of the system (called agents) into relatively small number of clusters. The clusters emerge and deteriorate in response to changes in the environment, and the processes of their formation and destruction are very short in time. During the periods of the clusters’ existence, the system’s dynamics is essentially low-dimensional due to synchronization between the agents in each cluster. An enormously complicated system is thus reduced to a finite-dimensional model with time-varying structure of the state vector. The low-dimensionality of the reduced model allows to control it by using classical methods, e.g. model-predictive or adaptive control.
This philosophy of complex systems control is illustrated on an experimental setup, called the “airplane with feathers”. The wings of this airplane are equipped with arrays of microsensors, microcomputers, and microactuators (“feathers”). The feathers self-organize into clusters by using a multi-agent consensus protocol; the aim of this coordination is to reduce the perturbing forces, affecting the airplane in a turbulent flow.


Introduction
Whereas interconnected systems have been studied for many decades, recently it has been realized that many of them are governed by similar principles and obey similar mathematical models and thus can be examined by the same methods. This has given rise to new paradigms in study of complex structures and their dynamics, manifested by Network Science [Barabasi, 2002, Newman et al., 2006, Easley and Kleinberg, 2010, Multi-Agent Systems (MAS) theory [Shoham andLeyton-Brown, 2008, Mesbahi andEgerstedt, 2010] and the theory of complex systems (or "complexity") [Bar-Yam, 1997]. Attracting enormous attention of researchers from various fields, these new emerging areas "break down barriers between physics, chemistry and biology ... and the soft sciences of psychology, sociology, economics,and anthropology" [Bar-Yam, 1997, Wooldridge, 2002]. Interest to system-theoretic properties of complex networks has been motivated by necessity to develop robust and resilient algorithms to control large-scale and spatially distributed systems such as e.g. continental power grids or smart infrastructures, functioning in the information rich world [Murray, 2003, Samad and Annaswamy, 2011, Annaswamy et al., 2017. To describe complexity effects, two types of models can be used: agent-based (networked, microscopic) and statistical (known under numerous alternative names, see below).

Statistical (Macroscopic) Modeling
Historically, statistical models of complex dynamics had appeared and been examined long before the emer-gence of network science. Since elementary parts of a complex system (e.g. individuals in a society, cars in dense urban traffic or molecules of some gas in a vessel) are numerous and their trajectories can be chaotic, the direct measurement of each part's characteristics is practically impossible. A natural idea is thus to describe the distribution of some variable of interest (e.g. the voters' preferences before an election or referendum, cars' or molecules' velocities etc.) or even some special metrics of such a distribution (e.g. temperature or average kinetic energy of the molecules). For example, models of fluid dynamics and continuum mechanics [Batchelor, 1967, Mase, 1967 describe the evolution of density and velocity of flows, where individual particles are replaced by indistinguishable infinitesimal elements of the flow. Similar statistical ("macroscopic") models naturally describe flows in transportation networks [Helbing, 1995].
In biology and social sciences, statistical (compartmental) models date back to early works on mathematical epidemics [Kermack and McKendrick, 1927] and mathematical bioscience [Rashevsky, 1935, Rashevsky, 1938. Behaviors between social or biological agents in such models are portrayed by interactions of multiple "species" or compartments [Jacquez, 1985].
Compartments are indecomposable entities that can grow or decline. For instance, in SIR models of epidemic spread the three compartments stand for the proportions of susceptible (S), infectious (I) and recovered (R) individuals; in the classical Lotka-Volterra model compartments are populations of predators and preys. Similar statistical models, arising in evolutionary game theory [Maynard Smith, 1982] and sociodynamics [Weidlich, 1971], are used to describe behaviors of social and economic agents. Such models can e.g. describe how the distribution of opinions evolves (paying no attention to the opinion of a particular social actor) and are sometimes also referred as Eulerian or "continuum-agent" [Como and Fagnani, 2011, Canuto et al., 2012, Mirtabatabaei and Bullo, 2012, Hendrickx and Olshevsky, 2016.

Agent-Based (Microscopic) Models
The networked or agent-based approach considers a complex system as a union of simple autonomous units, referred to as agents. 1 This bottom-up approach is broadly used in sociology since the pioneering works on consensus in social groups and evolution of social power [French Jr., 1956, Harary, 1959, Abelson, 1967, DeGroot, 1974, modeling large-scale biological populations (e.g. herds, flocks and schools [Reynolds, 1987]), ensembles of particles [Vicsek et al., 1995] and coupled oscillators [Kuramoto, 1991, Mirollo andStrogatz, 1990]. Statistical models can be thought as limits of agentbased models as the number of agents tends to infinity (in practice, the system becomes very large) and distinctions between individual agents are neglected. Usually the agents are mapped to nodes of some network, or interaction graph. The graph's arcs represent connections, or couplings between the agents. The coupling can stand for a physical interaction or presence of information flow between the two agents. In most interesting situations, the connections are sparse: each agent interacts with only a few number of neighbors, adjacent to its node in the graph, and has no information about the remaining part of the network. The topology of the graph can be static or time-varying, in the latter case, it can evolve independently of the agents or co-evolve with them [Belykh et al., 2014a]. It should be noted that structural and statistical properties of complex graphs, representing e.g. large-scale power grids or human relations in social groups, lead to a number of difficult problems on the borderline between statistical physics and graph theory [Watts and Strogatz, 1998, Barabási and Albert, 1999, Barabási et al., 2000, Albert and Barabási, 2002, Newman, 2003.
Many problems in this area has, in fact, been anticipated by extensive studies in social network analysis [Wasserman andFaust, 1994, Freeman, 2004]. However, the most intriguing problems are concerned with dynamics over complex networks.

Synchronization and Consensus
Perhaps, the most studied phenomena in dynamic complex networks are various effects of "self-organization" and "spontaneous order" [Néda et al., 2000, Strogatz, 2003] emerging from local interactions among the agents. The term "spontaneous order" has appeared in social philosophy and economics [Barry, 1982] and means self-organization of social or economic self-interested actors, established as a result of their independent actions rather than governmental regulations. In mathematical theory of dynamic networks, this term is typically understood as a kind of synchronization [Blekhman, 1988a, Blekhman et al., 1997, Strogatz, 2003, Fradkov, 2007, Pikovsky et al., 2001, Arenas et al., 2008, Wu, 2007, Fradkov, 2017, Boccaletti et al., 2018 between the agents' trajectories. Studies on self-synchronization among dynamical systems date back to the seminal experiments by Huygens in 17th century; the efforts to reproduce his experiments has led to a number of new results and open problems [Bennett et al., 2002, Czolczynski et al., 2011, Pena Ramirez et al., 2016. Several decades after Huygens' experiments, de Mairan [de Mairan, 1729] discovered the effect of forced synchronization (or entrainment) of plants' circadian clocks with an environmental periodic signals. Studies on entrainment of biological rhythms gave birth to the field of chronobiology and opened up the new perspectives in control of limit-cycle oscillators [Johnson, 1999, Roenneberg et al., 2003, Efimov, 2011, Sacre and Sepulchre, 2014, Proskurnikov and Cao, 2017c, Medvedev et al., 2018. The seminal work [Pecora and Carroll, 1998b] on master-slave synchronization of chaotic system has attracted enormous attention of physicists and engineers. The idea of chaotic synchronization has found numerous applications in telecommunications; this and other aspects of chaos "taming" can be found in surveys [Andrievsky andFradkov, 2004, Fradkov andEvans, 2005].
Special cases of synchronization are asymptotic coordinate synchronization and consensus [Olfati-Saber et al., 2007a, Ren and Beard, 2008, Ren and Cao, 2011] that take their origin in sociological studies Tempo, 2018].

From Synchronous to Clustered Behaviors
Synchronization is, however, not the only possible form of a complex system's organization. Much more complicated (and less studied) are various forms of clustering 2 , also called partial or cluster synchronization. Special case of clustering includes the situation, where the states or outputs of all agents should be separated by a positive distance, as in the problems of collision avoidance in swarms and mobile sensor deployment [Ny et al., 2012, Schwager et al., 2011. Much more interesting are situations where the number of clusters is much smaller than the number of agents; as will be discussed, such a clustering greatly enables efficient control of the system's dynamics using a few control inputs. In terms of statistical models with infinite number of agents, synchronization typically corresponds to a unimodal distribution with a "weak" tail, where clustering corresponds to a bimodal (polarization) or multimodal distribution. Paradoxically, models providing agents' splitting into a few clusters are much less studied than synchronization and consensus algorithms; the structural and dynamical properties of the networks leading to formation and deterioration of stable clusters still remain, to a great extent, unexplored. The models explaining cluster synchronization can be divided into several classes.
One can expect that clustering may be caused by the agents' heterogeneity. One can expect, for instance, that in a network constituted by agents of K ≥ 2 different types, agents of the same type should reach synchrony, constituting thus K different clusters. Such models of clustering have been examined e.g. in Aeyels, 2009, Lu et al., 2010], see also recent works [Lu et al., 2017, Du et al., 2018 and references therein. Heterogeneity does not necessarily implies different dynamics of agents. Heterogeneous dynamics of agents can be caused by presence of non-identical exogenous signals, that can be considered as control inputs, disturbances or exclusive information [Wu et al., 2009, Aeyels and Smet, 2011, Xia and Cao, 2011. In social networks, agents differ in their susceptibility to external influence, being fully or partially closed to assimilation of the peers' opinions , Friedkin, 2015, Parsegov et al., 2017. Besides agents, heterogeneity may occur in couplings, e.g. agents can split into communities or clusters such that couplings within a single community are stronger, in some sense, than couplings between different communities; different clusters may even repulse [Xia and Cao, 2011, Fläche and Macy, 2011, Altafini, 2012, Altafini, 2013, Ruiz-Silva and Gonzalo Barajas-Ramírez, 2018.
Homogeneous networks, however, also exhibit cluster synchronization. One of the reasons for that can be the symmetry of couplings, enabling stable invariant manifolds of partial synchrony [Belykh et al., 2000, Pogromsky et al., 2002, Pogromsky, 2008; similar effects have been earlier revealed in oscillator networks [Okuda, 1993]. Desynchronization and cluster formation in oscillator networks may be caused by presence of common noise [Gil et al., 2009].
Recently, it has been demonstrated that clustering is possible even if the agents and couplings are identical during special structure of the coupling functions that can e.g. vanish as the distance between agents grows (the principle of bounded confidence models [Hegselmann and Krause, 2002, Deffuant et al., 2001, Proskurnikov and Tempo, 2018 or the agents approach some stationary points [Proskurnikov and Cao, 2016].
The distribution of agents between clusters can be uniform or highly unequal, e.g. one huge cluster including all opinions and few clusters of 1-2 elements each.
In oscillator networks, such dispropor-tional structures of clusters are known as "chimera states" [Abrams and Strogatz, 2004].

Clusters and Control of Networks
Besides portraying important phenomena, observed in natural, social and robotic networks, models of synchronization and clustering play an important role in control of complex networks. It has been shown in [Zubarev, 1974] that no finite number of variables (degrees of freedom) can be sufficient to describe a non-equilibrium process with a variable number of degrees of freedom. At the same time, it is known that complex phenomena such as turbulence and spatiotemporal chaos [Xiao et al., 1998, Grigoriev et al., 1997 can be very controlled in a very "parsimonious" way. For instance, to suppress turbulence efficiently 3 , one needs a finite number of control inputs that is much less than the number of positive Lyapunov exponents [Xiao et al., 1998], provided that these controls are injected in proper points of the system. This seminal idea is known as pinning control [Grigoriev et al., 1997, Wu et al., 2009, Yu et al., 2014, Fu et al., 2017a.
In multiagent control, this idea is illustrated by algorithms of leader-following synchronization Beard, 2008, Lewis et al., 2014]: a group of agents applied a usual consensus algorithm in order to follows a single leader agent (influencing the remaining systems directly or indirectly), which can be externally controlled. Similar in spirit algorithms of containment control [Ren and Cao, 2011] gather a group of agents inside a convex polytope, spanned by the states of several leaders; controlling the leaders' trajectories, one then can control the motion of the whole swarm without precise synchronization of their trajectories with the leaders. A similar role in social networks is played by stubborn agents [Yildiz et al., 2013, Parsegov et al., 2017. The clustered structure of a network substantially simplifies its structure as a control system and enables relatively simple (low-dimensional) controllers, which can operate in uncertain conditions, as exemplified by modern methods of robust and stochastic model predictive control (MPC) [Mayne et al., 2000, Mesbah, 2018, Seron et al., 2018, adaptive control [Fradkov et al., 1999, Fradkov, 2007, data-driven and learning-based control [Benosman, 2018].

Forced Clustering and Slow-Fast Dynamics
Summarizing the discussion from the previous subsections, it can be said that "spontaneous order" effects such as synchronization and clustering in complex networks are not only crucial for understanding natural and social phenomena, but also important from control-theoretic prospective and enable efficient control of very complex dynamics. In spite of enormous progress in understanding these effects, as well as controllability properties of complex systems, most of the existing results are limited in several aspects.
First, the existing results are primarily focused on understanding the roles of couplings between the agents and connectivity properties of the interconnection topology in reaching full or partial synchronization. Much less understood is the role of environment, which can exert heterogeneous forces on one or several agents (one can imagine viscous friction forces, arising when agents move through some medium, or beams of particles that are bombarding agents). In turn, the environment can depend on the agents' states, e.g. neural cells can secrete neurotransmitter peptides, changing conductivity of the intercellular medium [Gonze et al., 2005, Vasalou andHenson, 2011], moving charged particles create a time-varying electromagnetic field etc. The arising clusters can be forced by the environmental "cues" and thus adapt to the changing environment.
Vortex-wave turbulence structures in fluid, gas, multiphase and plastic streams Khantuleva, 2015, Khantuleva andMeshcheryakov, 2016] can also be described in terms of birth and dissolution of clusters. Although cluster synchronization dramatically reduces the dimension of a control system, this dimension (and structure of the state vector) is time-varying. One paradigm to cope with such systems is hybrid systems theory [Goebel et al., 2009], imposing however a number of restrictions on the system's dynamics. In statistical physics, studies on cluster dynamics in fact date to Boltzmann's works [Gabrielov et al., 2008].
Second, most of the existing results consider systems with a single time scale, whereas in reality some processes within an agent or unfolding in the whole network can be much faster than others, as exemplified by spiking (very fast) and rest (relatively slow) dynamics in neuron ensembles [Pereira et al., 2007]. One way to describe fast processes is to consider them as pulses; such networks can be studied by methods of impulsive systems theory [Gelig and Churilov, 1998] or hybrid systems [Goebel et al., 2009]. Most studied systems of this type are pulse-coupled oscillator (PCO) networks that portray behavior of biological neurons and have found many applications in electronics [Mirollo and Strogatz, 1990, Kuramoto, 1991, Proskurnikov and Cao, 2017d, Núñez et al., 2015. Another approach is the method of singular perturbations, broadly used in physics and mechanics [Cole, 1968, Dyke, 1964. Networks with singularly perturbed agents have not been thoroughly studied, most of them stipulate linear dynamics of the agents [Yang et al., 2017, Rejeb et al., 2018 or cope with very special type of networks, e.g. Kuramoto-like oscillator ensembles [Dörfler and Bullo, 2012]. In practice, external control can be rather slow compared to the internal processes of cluster formation and deterioration. Understanding of interrelation between fast processes of clustering and slow-input controllability of the network remains an unsolved problem.

The Paper's Contribution and Organization
In this paper, we propose a new general "philosophy" of complex processes modeling, exploiting multi-agent dynamical models with variable structures 4 of clusters. The relevant mathematical concepts and brief literature survey are given in Section 2. This philosophy has been, to a great extent, inspired by a practical problem of reducing an aircraft's vibration in a turbulent flow. In Section 3, we present a novel approach to this problem, employing a special equipment ("feathers", consisting of miniature sensors, computers and actuators) installed on the wing surface and used to equalize the air pressure. Section 4 concludes the paper, discussing some open problems and directions for future research.

Clusters Formation and Related Problems
As has been already mentioned, many complex phenomena studied in natural, engineering and social sciences can be described as emergence and deterioration of clusters in complex dynamical networks. The clusters can emerge due to some endogenous structural and/or dynamical properties of the system (that is, the topology and structure of coupling), but can be also influenced by external factors (including e.g. control inputs). As a result of fast external processes, the clusters deteriorate, and the agents reassemble into a different configuration of clusters. Supposing that the processes of clusters' formation and destruction are sufficiently fast, so that the configuration of clusters remains stable for sufficiently long time, one can replace each cluster by a single "averaged" system and thus dramatically reduce the state space of the overall complex systems. This general program decomposes into a number of difficult problems, in particular: 1. why and how do the agents split into clusters? 2. how to predict the structure and number of clusters, given their previous configuration? 3. how fast is the cluster formation process? 4. how to control a cluster in a non-destructive way, preserving consensus of constituent agents? 5. which properties of the original complex system are inherited by its reduced-order model, obtained by "gluing together" agents from the same cluster and treating them as a single "macro-agent"?
The existing literature on cluster synchronization mainly focuses on the first two problems, which, how-4 The theory of cluster dynamics should not be confused with variable-structure systems (VSS) theory, focusing on sliding-mode control of uncertain systems [Emelyanov, 1967, Utkin, 1992. ever, are also far from being thoroughly studied. Problems 3 and 4 are mainly considered for special situations (e.g. linear dynamics of the agents). Problem 5 relates to the advanced theory of singularly perturbed systems, studying coupled dynamics on different time scales (dynamics of the cluster restructuring are very fast, compared to periods of the clusters' co-existence).

Multi-Agent Networks
Henceforth we confine ourselves to the special situation, where the system is constituted by a finite 5 (yet large) number of agents N ≫ 1, indexed 1 through N ; we use N = {1, . . . , N } to denote the set of agents. Agent i is described by the differential equationṡ (1) Here x i (t) ∈ R ni stands for the state vector of agent i ∈ N . The joint state vector of the system X(t) = [x ⊤ 1 , . . . , x ⊤ N ] ⊤ ∈ X = Rn thus belongs to the space of very high dimensionn = ∑ i∈N n i . Each agent can depend on the (finite-or infinite-dimensional) uncertain vector θ i ∈ Θ i , generally, time-varying (in some situations, it is convenient to consider θ i as a stochastic variable). The components of this vector can represent parameters of the agent or external signals, which can stand for e.g. reference trajectories or disturbances. The functions f i and θ i should, of course, be "good" enough to guarantee a solution's existence; both of them, however, may become discontinuous due to rapid changes in the structure of the agent and its environment. The relevant existence theorems can be found in classical monographs on functional equations and differential inclusions [Coddington and Levinson, 1955, Hale, 1977, Gelig et al., 2004. Along with the continuous-time model (1), discretetime agents may be considered whose states evolve as Continuous-time dynamics naturally describe physical processes, and discrete-time equations arise in models in digital devices and software components. Many "cyberphysical" systems arising in engineering include both physical and digital components, and hence are described by impulsive [Gelig and Churilov, 1998] or hybrid models [Goebel et al., 2009], where some state components evolve in continuous-time as in (1), and the others obey recurrent equations as in (2). Unless otherwise stated, we consider continuous-time agents (1). Each agent has two control inputs, which we refer to as "fast" and "slow" controls respectively. The first of them, u i (t) is used for coupling between the agents that provides global consensus or cluster synchronization (such algorithms will be considered below in detail). The second input U i (t) is used to control the agent. It can be e.g. identical for each agent in a cluster, treating thus clusters as "macro-agents", and depend on some averaged characteristics such as the average of agents in a cluster (as the number of agents N → ∞, such an average can be approximated by an integral over some domain). Such an example will be considered below in Subsection 2.5. Another paradigm (known as "pinning control") is to influence a single agent or few agents in each cluster, in this situation, U i (t) = 0 for most agents except for a few "pinned" leaders. Notice that in general the agents may be not only heterogeneous, but even have state vectors of different dimensions (more generally, belonging to different vector spaces). To describe clustering of agents, some "distance" between the agents has to be introduced. For this reason, the output z i (t) is introduced, containing some characteristics of interest. We always assume that which allows to measure distances between these outputs. In examples considered below, usually z i = x i . Besides z i , the agent has another output y i that is used for the interactions, or coupling, among the agents. At each time t, the set N i (t) ⊆ N is defined (possibly, empty), whose members we refer to as neighbors of agent i. Decision of agent i at time t is based on the outputs y j (t), "displayed" by the neighbors. 6 The set of outputs {y j (t)} j∈N i can be measured by the agent (e.g. a mobile robot can detect positions and velocities of its neighbors) or communicated by the neighbors. As a result of interaction with neighbors, each agent takes a decision or computes the control command Here F i , i ∈ N are functions, defined on appropriate spaces. The coupling algorithm (3) (also called protocol) is distributed in the sense that each agent uses only "local" information Y i (t) from its neighbors. Notice that F i can be random functions; also, in the case of constant sets N i one can consider more general functionals, depending on the whole trajectory {Y (s) : 0 ≤ s ≤ t}. Notice that the coupling (3) can stand for some physical connection between the agents. For instance, agents can represent mechanical oscillators (pendulums or vibrating using), installed on the same platform [Blekhman, 1988b].
In many situations, it is convenient to depict the relations j ∈ N i (t) (agent j influences agent i at time t) by directed arcs i − → j. A collection of such arcs E(t) along with the set of nodes N constitutes the directed interaction graph G(t) = (N , E(t)), referred also to as the topology of the network (1),(3).
For the subsequent definition, it is convenient to introduce the deviations between the outputs Here ∥ · ∥ stands for some predefined norm in the space of outputs. On finite-dimensional spaces, all norms are equivalent, however, from technical viewpoint computation of some norm may appear simpler. Most typically, ∥z∥ stands for one of the norms Definition 2. We say that agents i and j are (output) asymptotically) if any pair i, j ∈ M is synchronized in this sense. If this holds for all agents (M = N ), the system is said to be globally synchronized. The same definitions apply to ε-synchronization.
It should be noticed that synchronization from Definition 2 is a special case of coordinate synchronization. In this paper, we do not consider general concept of synchronization in dynamic systems; the interested readers are referred to monographs and surveys [Blekhman, 1988b, Blekhman et al., 1997, Fradkov, 2017. Notice that synchronization of two agents at time t * , in general, does not imply synchronization at any t ≥ t * : as the system evolve, synchronization can be destroyed. However, usually algorithms used to reach consensus provide the invariance of global synchronization: if z i (t * ) = z j (t * ) holds for any i, j, this relation retains its validity for t ≥ t * (more formally, it holds until the system is "singularly" disturbed, see below).
Some clusters M j may consist of only one agent. A (0, 0)-clustering is henceforth referred to simply as clustering (or cluster synchronization) of the agents.
Notice that in the second inequality of (5) we use the lower limit rather than the upper one, thus, for t * being large, agents from different clusters are δ-separated for any t ≥ t * . Obviously, ε-synchronization is a special case of (ε, δ)-clustering with the only cluster M 1 = N (in this case, δ ≥ 0 can be arbitrary). The most interesting situation of (ε, δ)-clustering is the case where, first, δ ≥ ε (the inter-cluster distance is greater than the clusters' diameters) and, second, the number of clusters k is relatively small compared to N = |N |. Effects of full or cluster synchronization (as well as other "spontaneous order" phenomena that are beyond the scope of this paper) among the agents are explained by presence of couplings among them. A natural question arises how a network synchronizes or splits into clusters. This problem has attracted enormous attention of different research communities in the recent decades. In spite of a substantial progress, the behaviors of quite simple models of dynamic networks (such as e.g. Kuramoto oscillators [Dörfler and Bullo, 2014a] or models of opinion dynamics [Proskurnikov and Tempo, 2018]) remain, in many extent, mysterious. In the next subsections, we consider several most studied situations.

The French-DeGroot and Abelson Models
One of the first models, explaining how consensus among the agents can arise, was suggested by French [French Jr., 1956]. Consider a group of N social actors, whose opinions at the stage t = 0, 1, 2, . . . are expressed by some scalar values x 1 (t), . . . , x N (t). The pattern of interpersonal relations is described by a graph G = (V, E), whose nodes are in one-to-one correspondence with the actors and each arc corresponds to direct communication between the actor at the head of the arc and that at its tail. At each step agents simultaneously calculate the means of their own opinions with opinions, communicated by their neighbors. Introducing the adjacency matrix A = (a ij ) of the graph, the dynamics of the ith opinion is given by (6) This model was investigated by F. Harary [Harary, 1959], who established the first graphtheoretic criterion of asymptotic consensus ("unanimity") which means that all of the opinions x i (t) converge to the same limit 7 : The iteration (6) is a special case of more general dynamics, referred to as linear opinion pooling [DeGroot, 1974]. DeGroot's model involves a rowstochastic 8 matrix W = (w ij ) and is given by stands for the stack vector of opinions. Establishing of consensus among the opinions x j (t) under any initial conditions X(0) means that the limit W * = lim t→∞ W t exists, and W * = 7 Consensus in this sense may seem a more restrictive condition than asymptotic synchrony from Definition 2, but in fact the convergence of the opinions to the same limit (7) and vanishing of their deviations appear to be equivalent in the case of linear algorithms. The corresponding fact remains valid for time-varying matrices W (t) and is known as the equivalence of weak and strong ergodicity of Markov chains (or infinite matrix products) [Seneta, 1981].
8 A square matrix W = (w ij ) is row-stochastic, if all of its entries are non-negative w ij ≥ 0 and the sum in each row ∑ j=1 w ij = 1. Any such matrix is the transition matrix of a stationary Markov chain [Norris, 1998], where w ij is a probability of transition from the state i to state j.
(1, 1, . . . , 1) ⊤ π ⊤ * , where π ⊤ * = (π 1 , . . . , π N ) is a nonnegative row vector. Obviously, W * W = W * , and hence π ⊤ * W = π * . Thus π * is the Perron-Frobenius left eigenvector of W , corresponding to the maximal eigenvalue 1. Since W * is row-stochastic, one also has ∑ i π i = 1. It can be easily shown that π * is the only eigenvector with such properties. Notice also that the consensus opinion, on which the agents converge, is Conditions ensuring the existence of vector π * have been thoroughly investigated in the theory of Markov chains; the aforementioned consensus property is equivalent to the ergodicity of the Markov chain with the transition matrix W , whose unique stationary distribution is π * . Starting at any initial distribution π 0 = (π 0 1 , . . . , π 0 N ), the chain converges to the stationary distribution π * = π 0 W * = lim t→∞ π 0 W t . A necessary and sufficient condition on W , guaranteeing the ergodicity, is being a SIA matrix (stochastic, indecomposable, aperiodic) [Ren and Beard, 2008]. From the algebraic viewpoint, this means that all eigenvalues of W except for λ 1 = 1 are strictly stable (|λ j | < 1, j = 2, . . . , N ), and λ 1 has multiplicity 1 (such matrices are also called fully regular [Gantmacher, 2000].) In terms of the Markov chain, this means that all recurrent states [Norris, 1998] of the Markov chain communicate to each other and aperiodic. A more convenient graphical condition can be obtained, introducing the interaction graph of the multi-agent system. The set of the ith agent's neighbors is N i = {j : w ji > 0} (obviously, the ith agent is affected by the jth one if and only if w ij ̸ = 0). Communicating recurrent classes correspond to the only source strongly connected component [Harary et al., 1965], from which all other components can be reached by paths. The existence of such a component is equivalent to the existence of a directed (out-branching) spanning tree in a graph [Harary et al., 1965, Ren and Beard, 2008, Blondel et al., 2005, Cao et al., 2008. Consensus is reached if and only if this source component is aperiodic, that is, greatest common divisor of all cycles' lengths is 1. In the simplest situation where w ii > 0 ∀i, aperiodicity automatically holds, and the system reaches consensus if and only if the graph has a directed spanning tree [Ren and Beard, 2008, Blondel et al., 2005, Cao et al., 2008. The mulit-agent system (8) can be easily rewritten in the general form (2), (3), for instance, Hence, the evolution of an agent's opinion is driven by its deviations y j − y i from the neighbors' outputs 9 . This coupling mechanism is often called diffusive [Pogromsky and Nijmeijer, 2001]. A continuous-time counterpart of (10), replacing the opinion's increment by the instantaneous "velocity", iṡ Here L = (l ij ) is a so-called Laplacian matrix [Ren and Beard, 2008, Mesbahi and Egerstedt, 2010, Olfati-Saber and Murray, 2004a, its off-diagonal entries l ij = −ε ij are non-positive (i.e. ε ij ≥ 0) and the rows sum to zero: It is remarkable that the algorithm (11) first appeared and was studied in [Abelson, 1964], long before the recent "boom" in multi-agent and networked control. Consensus in (11) holds if and only of the interaction graph, corresponding to the adjacency matrix (ε ij ), has a directed spanning tree. This criterion has been formulated in [Abelson, 1964], however, the proof was incomplete. The first correct proof appeared in [Agaev and Chebotarev, 2000]. Similarly to the discrete-time case, the consensus value can be represented as (9), where π * ≥ 0 is the (unique) left eigenvector of the Laplacian matrix L, such that π ⊤ * L = 0, ∑ j π j = 1. The detailed history and analysis of properties of the algorithms (8),(11) is available in . The consensus algorithm (11) is in fact a special case of speed gradient (or gradient descent) dynamics [Fradkov et al., 1999, Fradkov, 2005, Khantuleva and Shalymov, 2017. Consider the "cost function" that penalizes disagreement between the agents A straightforward computation shows [Olfati-Saber andMurray, 2004b, Olfati-Saber et al., 2007b] that (41) is nothing else than the conventional speed-gradient algorithṁ In other words, each agent is constantly moving in the 9 The agent's equation could also be written as x i (t+1) = u i (t), from control viewpoint, (10) is more natural (the state usually cannot be controlled directly) and stands for the discretized integrator. direction of fastest descent of the cost function. Since and the graph is connected, the set whereQ = 0 (equivalently,ẋ i = 0) coincides with the synchronous manifold, i.e. the set of all X such that x i = x j ∀i, j. An important property of the algorithms (8) and (11) is the robustness of synchronization against external disturbances. For instance, the algorithṁ where θ(t) ∈ R N is the vector of bounded disturbances, provides ε-synchronization with the accuracy ε proportional to sup t≥0 max i,j |θ i (t) − θ j (t)|.
In the case of vanishing at infinity θ(·), asymptotic synchronization is preserved. Mathematical results on consensus robustness for continuous-time and discrete-time cases are available in [Shi and Johansson, 2013a]. Besides this, consensus algorithms are robust against delays in measurements , Münz et al., 2011, Tian and Liu, 2008, Blondel et al., 2005, Proskurnikov, 2012, Shi and Johansson, 2013a, Amelina et al., 2015, Seuret et al., 2008, Proskurnikov and Matveev, 2018. Another important property of the consensus protocol has been established in [Fu et al., 2017b]. Assuming that the initial vector x(0) is a continuous random variable with a finite entropy h(x(0)), the same holds for any x(t) (also being random). Furthermore, the entropy is non-increasing along the trajectories, in particular, h(x(t)) ≤ h(x(0)). Another relation between the entropy and consensus will be discussed in Section 2.3.
In other words, the asymptotic convergence rate depends on the second largest eigenvalue of W .
Similarly, in the case of continuous-time dynamics (11), one has [Olfati-Saber et al., 2007a, Ren andBeard, 2008] where λ 2 stands for the eigenvalue of L with the second minimal real part 10 . For practical reasons, it is often more important to know the convergence time of the algorithm, namely, how long will it take the agents to reach ε-synchrony. Instead of the maximal deviation max i,j |x i − x j |, it is often more convenient to evaluate the deviation ∥x(t) − x * ∥; the choice of norm here can be arbitrary, however, for technical reasons the norm ∥ · ∥ = | · | ∞ is usually considered Tsitsiklis, 2011, Nedić and.
Using (15), one has Obviously, the agents reach ε-consensus in no longer than T (ε/∥x(0) − x * ∥) iterations. A question arises how to estimate ρ and T (r) in a general situation 11 . The most studied is the case of bidirectional graphs (i ∈ N j ⇔ j ∈ N i ) with self-loops (i ∈ N i ), where each agent assigns equal influence weights to its neighbors. w ij = 1 |N i | ∀j ∈ N i . The result by [Olshevsky and Tsitsiklis, 2011] states that in such a situation, where γ > 0 is a constant, independent of N or any specific graph. Moreover, this asymptotical estimate cannot be improved: there exist graphs, on which ρ ≥ 1 −γN −3 . Using this result, it can be shown [Olshevsky and Tsitsiklis, 2011] that however, the tightness of this estimate remains an open problem (it is known, however, that T (r) = Ω(N 3 log r) for an appropriate graph with N nodes.
The problem of convergence rate's estimation becomes even more challenging for randomized algorithms, where either matrix W (t) is random (as e.g. in gossiping algorithms, surveyed in [Proskurnikov and Tempo, 2018]) or the system is disturbed by random noises. In such situations, synchronization is often understood in meansquare sense; however, one may also try to estimate (random) time of reaching ε-consensus, provided that disturbances are sufficiently small compared to ε.

Extensions and Applications
The properties of SIA matrices [Cao et al., 2008, Jadbabaie et al., 2003, Ren and Beard, 2008 and Lyapunov-based arguments [Moreau, 2005, Blondel et al., 2005, Münz et al., 2011 allow to obtain consensus criteria under non-stationary interactions among the agents, where the rowstochastic matrix W = W (t) in (8) or the Laplacian matrix L = L(t) in (11) evolves; the relevant results can be found in Beard, 2008, Proskurnikov andTempo, 2018]. These criteria require the existence of spanning tree in the union of interaction graphs over sufficiently long period and in general are only sufficient for reaching agreement. Necessary and sufficient conditions are known only for special types of interaction graphs [Matveev et al., 2013, Hendrickx and Tsitsiklis, 2013, Shi and Johansson, 2013b.
It should be noticed that dynamics similar to (8) and (11) naturally arise in many applications. Consensus algorithms (called also protocols) similar to (8) were independently introduced by J. Tsitsiklis [D. Bertsekas andTsitsiklis, 1989, Tsitsiklis, 1984] in some distributed algorithms for multi-processor clusters. Similar dynamics were suggested for modeling of alignment in ensembles of moving particles [Vicsek et al., 1995, Helbing, 2001, that can be considered as a deterministic analogue of phase transitions in thermodynamics.

From Consensus to Clustering
While the phenomenon of consensus has been studied up to certain exhaustiveness, opinions of social actors often do not reach any agreement but rather form highly irregular fractions or clusters of different sizes. Along with the consensus problem, the problem of describing clustering mechanisms has arisen in mathematical sociology [Abelson, 1967], where it is now referred to as the community cleavage problem [Friedkin, 2015] or Abelson's diversity puzzle [Kurahashi-Nakamura et al., 2016]. Since convexcombination mechanisms (8),(11) have been widely adopted in opinion dynamics modeling, a question has arisen which factors lead to clustering of opinions under such a protocol. The proposed explanations can be (roughly) divided into several major classes.

Clustering Due to Lack of Connectivity
The most obvious explanation for absence of consensus is the absence of connectivity in the network, e.g. if the graph is composed of two or more disjoint components (equivalently, W is a block-diagonal matrix), then the opinions in these components may become disparate. Less obvious is the situation where the graph is weakly connected yet does not have a directed spanning tree, that is, its condensation [Harary et al., 1965] has several source components, and any other component is reachable from them (equivalently, the graph is covered by a spanning forest of several trees). As a simple example, one may consider the system (8) with several stubborn agents. Agent i is stubborn if w ii = 1, whereas w ij = 0 ∀j ̸ = i. Such an agent is anchored at its initial opinion x i (t) ≡ x i (0), being insusceptible to social influence; at the same time, stubborn agents can influence the others.
In presence of several stubborn agents with heterogeneous opinions, the group does not reach consensus, but rather demonstrates clustering behavior. Typically, the opinions of normal ("open-minded") agents get into the convex hull, spanned by the stubborn agents' opinions. This containment property of consensus dynamics is explored in problems of containment control with multiple static or dynamic leaders Cao, 2011, Proskurnikov and. Notice, however, that there is no simple relation between the number of stubborn agents and number of clusters. The structure of clusters, formed by consensus-like algorithms over weakly connected graphs, depends on the matrix of spanning forests, whose computation for large-scale graph is challenging Agaev, 2002, Chebotarev andAgaev, 2014].

External Inputs and Information
The idea of stubborn agents can be further extended, replacing the (static) radical opinions by additional external inputs, static or dynamic, that influence some of the agents. Consider, for instance, the system (14), where θ i (t) = 0 for all agents except for i ∈ I 0 , where I 0 is a set of several "informed" leaders, whose inputs {θ j (t) : j ∈ I 0 }, are mutually different. One can expect that such a group also exhibits cluster synchronization; the relevant conditions have been found in [Smet and Aeyels, 2009, Xia and Cao, 2011, Aeyels and Smet, 2011. A similar idea lies in the heart of the so-called Friedkin-Johnsen (FJ) model [Friedkin and Johnsen, 1999], squarely based on the DeGroot dynamics (8), Similar to (8), here W is a row-stochastic matrix, A is a non-negative diagonal matrix, 0 ≤ a ii ≤ 1, and u is a constant vector. In the original model from [Friedkin and Johnsen, 1999], u i = x i (0), and the diagonal entry a ii is treated as an agent's susceptibility to opinion assimilation. If a ii = 0, the agent is totally stubborn x i (t) ≡ x i (0). A maximally susceptible agent with a ii = 1 assimilates opinions of the other agents in the same way as in (8) x An agent with a ii ∈ (0, 1) is partially "anchored" at its initial opinion u i and factors it into every iteration of the opinion. As discussed in , in general u i may be considered as a "prejudice" of agent i, influencing dynamics of its opinion. In spite of its simplicity, many properties of the model (16), including graph-theoretic conditions for its stability, convergence rate and the structure of final opinions have been studied only recently [Parsegov et al., 2017.
Recently, extensions of the Friedkin-Johnsen models have been proposed, allowing to describe dynamics of multidimensional opinions, representing individual positions on several logically related topics [Parsegov et al., 2017]. (8) and (11), allowing the agents' splitting into two clusters with opposite opinions ("polarization") is known as the Altafini model [Altafini, 2012, Altafini, 2013, Liu et al., 2017, Xia et al., 2016. The discrete-time Altafini model is coincident to (8) with the only difference is that W may have negative off-diagonal entries w ij , i ̸ = j, however, the diagonal weights are non-negative and the matrix of absolute values (|w ij |) i,j is row-stochastic. The continuous-time Altafini model is coincident to (11), however, L is a signed Laplacian matrix, which means that l ij , i ̸ = j, may be positive or negative and l ii = ∑ j |l ij |. Note that the Altafini model steps away from the convex combination mechanism, in particular, the convex hull of the opinions can expand.

Balance and Negative Weights An elegant extension of the consensus algorithms
In the case of non-negative coefficients (respectively, w ij or ε ij = (−l ij )), the Altafini model reduces to the usual consensus algorithm. Otherwise, the positive sign of a coupling gains can be considered as "trust" or "friendship" between the corresponding agents, whereas the negative gain means "distrust" or "enmity". If the graph has a directed spanning tree and is structurally balanced [Harary, 1953], i.e. agents split into two "hostile camps" such that members of the same camp are "friends" (w ij ≥ 0), and the members of different camps are "enemies" (w ij ≤ 0), then the Altafini model provides polarization of the opinions: members of the two camps agree on the opposite values x 0 and −x 0 , where x 0 depends on the initial conditions and the graph's structure. If these conditions do not hold, the model can provide a more complicated type of clustering, called "interval bipartite consensus" [Meng et al., 2016]. A degenerate situation is also possible, where the system is exponentially stable and all opinions converge to 0.
For the detailed analysis of Altafini models over static and time-varying graphs the reader is referred to , Meng et al., 2016, Xia et al., 2016, Proskurnikov and Cao, 2017a, Proskurnikov and Cao, 2017b, Liu et al., 2017.

Nonlinear Couplings
It can be noticed that the aforementioned models of clustering stipulate some heterogeneity of the agents (e.g. presence of heterogeneous disturbances, different levels of attachment to the initial opinions) or their couplings (e.g. positive/negative). Also, the number of clusters in such models is usually very large (in general, each agent can form its own cluster), except for the Altafini model, which can lead to polarized opinions (with only two clusters). It appears, however, that the cluster synchronization (with arbitrary number of clusters, which can be greater than two but much smaller than the number of agents) can be provided by special structure of couplings among the agents. Unlike (8) and (11), the relevant models of clustering are essentially nonlinear. The nonlinearity of couplings is, however, is not sufficient for clustering: consensus is guaranteed by many nonlinear algorithms, preserving the convex hull [Moreau, 2005, Lin et al., 2007b, Matveev et al., 2013, Proskurnikov, 2013. A number of clustering models can be written aṡ where φ(·) is a (nonlinear) coupling maps and ε ij ≥ 0 are coupling gains, encoding the interaction graph. In the case where φ(ξ, ζ)(ξ − ζ) > 0 for any ξ ̸ = ζ and φ is continuous, this protocol guarantees consensus [Lin et al., 2007b, Matveev et al., 2013. It is remarkable that algorithm (17) naturally arises in some models of statistical mechanics, e.g. considering x as a discrete probability distribution ( ∑ i x i = 1) and as-suming that φ(x i , x j ) = log(x j /x i ) = log x j − log x i and ε ij = ε ji , one obtains the speed-gradient dynamics maximizing the conventional Shannon entropy . Consensus at the uniform distribution corresponds to the maximum of the entropy. Remarkably, similar entropy-maximizing (MAXENT) algorithms for other types of entropies arising in statistical physics lead to non-uniform distributions, corresponding to clustering behavior Fradkov, 2018].
One way to provide clustering is to require clustering is to assume that φ(x) = 0 for |x| ≥ R, which is in harmony with effects of homophily and biased assimilation [Dandekar et al., 2013] in social dynamics: a social actor readily accepts opinions of like-minded individuals (agents j such that |x j (t) − x i (t)| < R), examining the deviant opinions with discretion or rejecting them. The discrete-time counterpart of this model was first proposed in [Krause, 2000] and is now referred to as the Hegselmann-Krause model [Hegselmann and Krause, 2002] x The model (18) is a modification of DeGroot's dynamics (8), where the graph is state-dependent and evolves together with the agents' opinions (such networks are called co-evolutionary or coevolving).
Bounded confidence models (in continuous and discrete time) have been extensively studied in the literature [Hegselmann and Krause, 2002, Deffuant et al., 2001, Lorenz, 2007, Blondel et al., 2009, Blondel et al., 2010, Ceragioli and Frasca, 2012, Motsch and Tadmor, 2013, Jabin and Motsch, 2014, Frasca et al., 2016, however, their nonlinear structure makes them very difficult for analysis; in particular, the results. Some results on stability analysis of fixed points in the state space (standing for possible configurations of opinion clusters) are given in [Blondel et al., 2009, Blondel et al., 2010, Ceragioli and Frasca, 2012; it was noticed, in particular, that some trajectories can converge to unstable equilibria. A novel model of opinion dynamics with bounded confidence has been proposed in [Pilyugin and Campi, 2018]. In the latter model, agents' attitudes belong to a symmetric interval, e.g. [−1, 1], and evolve as follows Here h > 0 is a constant, and N i (x) is the same as in (18). Two main differences with the model (18) are the opinion saturation at the extremal values and the effect of opinion reinforcement: if N i (x(t)) = {i}, agent i gets no counter-argument against its attitude and tends to "strengthen" it: It has been shown in [Pilyugin and Campi, 2018] that this model leads to polarization of opinions at the extremal values ±1; the structure of arising clusters can be determined. Unlike the model from [Ceragioli et al., 2016] with similar properties, the model from [Pilyugin and Campi, 2018] does not employ negative couplings between the agents. There are other models with nonlinear couplings, admitting clusters. For instance, a model proposed in [Fläche and Macy, 2011] allows repulsion between distant opinions, whereas close opinions attract 12 . The structure of clusters in this model still remains a challenging question; in the experiments reported in [Fläche and Macy, 2011], no more than 5 clusters are reported. Another model of cluster dynamics has been proposed in [Proskurnikov and Cao, 2016] Here α > 0 is a positive constant and h(x) is a coupling function, controlling the structure of clusters. In the simplest situation of linear coupling maps h(x) = kx, k ̸ = 0, the model (19) becomes a special case of the consensus algorithm (11), where l ij = −|k|a ij . Consensus is the most typical outcome (under the assumption of the network's connectivity) whenever h is strictly monotone. Otherwise, the dynamics (19) allows clustered equilibria that can be classified [Proskurnikov and Cao, 2016] and are determined, to a great extent, by the stationary points of h (i.e. the points where h ′ (x i ) = 0). Notice that an agent starting at such a point is usually 13 stubborn: (19). Under some conditions (e.g. the all-to-all communication case where a ij = 1 ∀i, j), it has been proved [Proskurnikov and Cao, 2016] that any of the solutions converges to one of the equilibria; the general case is a topic of ongoing research.
In this section, we consider only a few general classes of agents, for which synchronization problems have been most studied.

Linear Stationary Networks
Recall that system (11) may be considered as a group of agentṡ coupled via the following protocol A natural extension of the first-order integrator agent (20) is a general linear time-invariant agenṫ where x j (t) stands for the state vector of the agent, u j (t) ∈ R m is its control or input, and y j (t) ∈ R l is some output. The state vectors are generally unavailable for measurement, and the only information agents can use comes to their own and neighbors' outputs, or even deviation among the outputs. Introducing joint stack vectors and U (t) = [u 1 (t) ⊤ , . . . , u N (t) ⊤ ], a natural extension of the protocol (21) is as follows Beard, 2008, Olfati-Saber et al., 2007a]: Here K is a m × l-matrix and ⊗ stands for the Kronecker product of two matrices [Laub, 2005]: if A = (a ij ) anf B are arbitrary matrices, their product is a 12 B . . . a 1n B  a 21 B a 22 B . . . a 2n As we have discussed, the matrix L always has a zero eigenvalue λ 1 = 0, corresponding to the eigenvector (1, 1, . . . , 1) ⊤ ∈ R N . This is the only eigenvector at 0 if and only if the interaction graph of the network has a directed spanning tree. In this case, all the remaining eigenvalues of λ 2 , . . . , λ N lie in the open right halfplane of C: Re λ j > 0 [Agaev and Chebotarev, 2005]. We are interested in asymptotic synchronization of the agents' states (the synchronization outputs are z i = x i ); henceforth we refer to is as synchronization among the agents (22) for brevity. The following fundamental result [Fax and Murray, 2004, Olfati-Saber et al., 2007a, Li et al., 2010 establishes criterion for such a synchronization.
Theorem 1. The protocol (23) establishes asymptotic synchronization among the agents (22) for any if and only if two conditions hold: 1. the interaction graph has a directed spanning tree (and hence λ 1 = 0 is a simple eigenvalue of L); 2. for other eigenvalues λ 2 , . . . , λ N of L the matrices A − λ j BKC are Hurwitz. 14 Theorem 1 can be extended to more sophisticated models of networks, including those with time-delays [Olfati-Saber andMurray, 2004a, Michiels et al., 2009b].
Another approach to analysis of linear time-invariant networks is the frequency-domain technique [Tian and Liu, 2008, Tian and Zhang, 2012, Tian, 2012,Lestas and Vinnicombe, 2010,Münz, 2010, operating with Fourier or Laplace transforms of the solutions x j (t). Some extensions to heterogeneous agents and time-varying topologies can be found in [Wieland et al., 2011, Ren and Cao, 2011, Trentelman et al., 2013.

Nonstationary Networks with Static
Graphs and MSF Approach Whereas Theorem 1 appeared (in a modified form) only in [Fax and Murray, 2004] (see also [Wu, 2007]), its non-stationary extension, based on the Master Stability Function (MSF) approach has been widely known in physical literature since the seminal paper [Pecora and Carroll, 1998a], offering a synchronization criterion for the agentṡ applying a non-stationary distributed algorithm It should emphasized that the topology of the network here remains fixed, whereas the agents and coupling gains may evolve. The system (24),(25) arises in [Pecora and Carroll, 1998a] from linearization 15 of a time-invariant nonlinear network around a non-equilibrium solution (e.g. a cycle) For any z ∈ C consider a system of ordinary differential equationṡ

ξ(t) = [A(t) − zB(t)K(t)C(t)]ξ(t),
and let λ max (z) be its maximal Lyapunov exponent. The function λ max is referred to as the master stability function (MSF) of the system. In the stationary case λ max (z) is simply the maximal real part of the eigenvalues of the matrix A − zBKC and hence λ max (z) < 0 if and only if this matrix is Hurwitz. Condition 2 in Theorem 1 can be now reformulated as follows: λ max (λ j ) < 0 for j = 2, . . . , N . In this form it appears to be valid in the non-stationary case [Pecora and Carroll, 1998a]: Theorem 2. The protocol (25) establishes asymptotic synchronization among the agents (24) if and only if two conditions hold: 1. the interaction graph has a directed spanning tree (thus λ 1 = 0 is a simple eigenvalue of L); 2. for other eigenvalues λ 2 , . . . , λ N of L one has λ max (λ j ) < 0.
In the case where the coefficients are periodic, the Lyapunov exponents are negative if and only if the Floquet multipliers are less than 1 in modulus; this condition can be validated by e.g. solving the periodic Riccati equation [Yakubovich et al., 2007]. Notice that Theorem 2, applied to (26), guarantees only local conditions for synchronization; the estimation of the relevant basin of attraction remains a non-trivial problem. It is remarkable that a relaxation of the conditions from Theorems 1 and 2, in some situations, leads to cluster synchronization. This happens when the graph has some symmetries, which means that SL = LS for some permutation matrix S. Existence of symmetries enables the existence of clusters, whose structure is determined by the graph's symmetry group [Pecora et al., 2014]. These clusters are stable (locally in the case of nonlinear agents) if the conditions of Theorem 2 hold for some spectral points of the Laplacian matrix. It is remarkable that results on global stability of such clusters are available only in special situations [Pogromsky et al., 2002, Pogromsky, 2008.

Synchronization among
General Nonlinear Agents In spite of significant progress in recent years, see e.g. [Arenas et al., 2008, Belykh et al., 2014b, the problem of reaching synchronization or other type of regular behavior in a general complex network where both nodes and couplings may be nonlinear, and the topology may evolve, still remains a challenge. A wide class of such networks can be described by the following mathematical model. Let the dynamics of nodes bė coupled via the following non-stationary protocol Here φ jk (·) are some mappings (in general, nonlinear), referred to as coupling maps or couplings, and ε jk (t) are coupling gains that characterize the couplings' "strengths" and define, in particular, the interaction topology (generally time-varying and disconnected). The goal is to disclose conditions, providing asymptotic state synchronization between the agents. As was discussed in the foregoing, for stationary networks (ε jk = const) synchronization can be proved locally, under some technical assumptions, via the MSF approach [Pecora and Carroll, 1998a]. More subtle techniques, applicable to time-varying topologies, were proposed in [Lü andChen, 2005, Lu et al., 2007].
These methods, however, allow to estimate neither the basin of attraction for a synchronous solution nor the convergence speed. The results, establishing global synchronization in a general network (27),(28), mostly employ either contraction principles [DeLellis et al., 2011] or Lyapunov functions. Many result explicitly or implicitly employ "incremental" Lyapunov function of the formV (X) = where V stands for the positive definite quadratic form or another radially unbounded function. Often one can prove the protocol (28) to satisfy an "incremental" quadratic constraint: Suppose, for instance, that the interaction graph is undirected and the couplings are anti-symmetric assume also that φ jk (y) ⊤ y ≥ 0. A straightforward computation yields that which is a special case of (29). More subtle costraints (29) arise when the coupling obey general sector inequalities [Proskurnikov, 2014]. Assuming that (29) holds, to construct an incremental Lyapunov function V (X), it suffices to find such a function V (x) that for any solutions x 1 , u 1 , y 1 and x 2 , u 2 , y 2 of (27) one haṡ This principle is applicable to many networks with linear couplings φ jk (y) = Ky [Hamadeh et al., 2012, Stan and Sepulchre, 2007, Belykh et al., 2006b, Belykh et al., 2006a, Belykh et al., 2007, Belykh et al., 2005.
The condition (31) is a special case of incremental dissipativity [Stan and Sepulchre, 2007] which extends the conventional Willems dissipativity [Willems, 1972] has proved to be an important tool in network analysis [Liu et al., 2011.
In the case of linear agents (22) and nonlinear couplings (28), considered in [Proskurnikov, 2014, Proskurnikov andMatveev, 2015], incremental Lyapunov functions enable to derive counterparts of the seminal circle and Popov's crite-ria [Gelig et al., 2004]; in the discrete-time, counterparts of the Tsypkin and Jury-Lee criteria are obtained [Proskurnikov and Matveev, 2018]. Notice that up to now we have dealt with homogeneous (identical) agents.
In the case where agents are passive in Willems sense, some strong synchronization criteria, applicable to heterogeneous agents as well, were obtained in [Arcak, 2007, Chopra and Spong, 2006, Proskurnikov and Mazo, 2017. In general, heterogeneity of the agents is recognized as a potential reason for clustering, e.g. K different types of systems can be clustered into K synchronous groups [Lu et al., 2010, Lu et al., 2017, Du et al., 2018.

Consensus-preserving Control and Slow-Fast
Dynamics in Networks Up to now, we have considered a group of agents, reaching a consensus at some point; as follows from (9), this point depends on the initial condition and the structure of coupling matrix; the latter dependence is highly nonlinear. A natural question arises whether it is possible to control the multi-agent system (via additional inputs U i (t), as in (1)) in such a way that, first, synchrony among the agents is not destroyed and, second, the consensus value is driven towards the desired trajectory. The answer to this question is affirmative, and the construction used to solve this problem is known as the mean-field control [Kadanoff, 2009]. To keep things simpler, we confine ourselves to the case of balanced Laplacian matrix in (11), which means that ⊤ N L = 0, and therefore (9) holds with π * = N −1 N (average consensus) [Olfati-Saber et al., 2007a]. Notice first that an identical control input, added to all agents in (11), does not destroy consensus: Here we used the fact that L N = 0 and therefore e −Lt N = N . Since, by assumption, e −Lt x(0) → x * and x * = x 0 N , the asymptotic synchrony is preserved (furthermore, agents need the same time to reach εsynchrony as in the uncontrolled system (11)). Now it can be noticed that the consensus value x 0 = N −1 N x(t) obeys the integrator dynamicṡ If one is capable to measure the average value x 0 (t) by some sensor (e.g. by organizing an opinion poll in the case of social network), it becomes possible to control it by applying the proper control input U (t) to all of the agents. For instance, one is able to drive it to a predefined value x 0 des (constant or slowly changing) by using e.g. the conventional PI controller Notice the properties of such a control algorithm: -each agent directly interacts only with neighboring agents from N i = {j : l ij ̸ = 0}; -implicitly, each agent is also influenced by a (common) control input, preserving consensus; -the external control input depends on an averaged characteristics of the cluster and may be considered as a "mean field", depending on all the agents; -the resulting motion is a superposition of the two motions, caused respectively by the local couplings and the mean-field control; -the "mean field" control depends only on some averaged characteristics; A similar idea can be applied in a more general situation of multi-dimensional agents. Consider a group of linear agents, having a common input U i (t) = U (t) (33) Assume that the agents are coupled via the protocol (28), obeying symmetry conditions (30). Introducing the new variablesx i (t) = x i (t) − B 0 ∫ t 0 U (s) ds, the system (33),(28) reduces to (22),(28), which means that if the algorithm (28) establishes consensus, this consensus is not destroyed by the meanfield control U (t). Introducing the function x 0 (t) = N −1 N x(t), which is no longer constant, and noticing In other words, similar to the first-order case, the dynamics of the (single) cluster's center is decoupled from the "centripetal" consensus dynamics of individual agents. This separation principle opens up the perspective to control the large-scale formation of agents efficiently by controlling the formation's center x 0 (t), being an averaged characteristics of the multi-agent system. A similar principle can be applied to continuum agents (N → ∞) under the assumption of coupling symmetry [Hendrickx and Olshevsky, 2016]; in the latter case, the average value has to be replaced by the integral over the cluster of agents. As the dynamics of an agent (and the cluster's center) become highdimensional, partially uncertain and various constraints occur that have to be taken into account, for controlling U (t) efficient methods of advanced control theory have to be employed, such as e.g. MPC and adaptive control.
The idea of mean-field control thus appears to be a simple yet efficient technique to control the dynamics of clusters without destroying them; its practical implementation in the case of multiple clusters (discussed in the next sections), heterogeneous dynamics of agents or asymmetric coupling remains an open problem. A disadvantage of this approach is the necessity to create a "homogeneous" field, equally influencing all agents. An alternative approach is to control a single agent in the community; it can be expected that (under some conditions) the couplings between the agents make the whole group follow this dedicated "leader". The relevant algorithms are known as leader-following or reference-tracking consensus protocols. In the case of non-trivial desired trajectory of the leader one, however, either needs to estimate its derivative or use sliding-mode control. We do not consider such algorithms in this paper and refer the reader to Beard, 2008, Lewis et al., 2014].
Another conjecture, inspired by the aforementioned decoupling property, is that general complex systems can be decoupled into low-dimensional dynamics and "slow" of a few cluster's centers (average values, or integrals of the states over the clusters) and high-dimensional and "fast" dynamics of individual agents, forming a cluster due to some coupling law. Whereas in general situation simple separation of these dynamics, similar to (34), does not seem to be possible, methods of averaging and singular perturbations could be employed, allowing to separate slow and fast processes in dynamical systems [Mitropolsky, 1967, Blekhman, 2000, Kokotovic et al., 1986, O'Malley, 1991, Cole, 1968. It should be noticed, however, that the existing theory of singularly perturbed systems is mainly focused on two problems. The first problem, pioneered in [Tikhonov, 1948] is the convergence of the solutions of the "perturbed" system to "unperturbed" ones (such a convergence appears to be uniform only on compact intervals, which does not allow to investigate asymptotic behaviors of the perturbed systems) [Lizama andPrado, 2006, Parand andRad, 2011].
The second problem developed in [Klimushev and Krasovskii, 1961, Khalil, 1981, Kokotovic et al., 1986 is the global asymptotic stability of singularly perturbed systems with a unique equilibrium.
Behavior of other types of attractors under singular perturbations has been studied only in special situations [Fenichel, 1979, Krupa and Wechselberger, 2010, Jardon-Kojakhmetov et al., 2016.
Recently, a progress has been achieved in examining pendulumlike systems with multiple equilibria under singular perturbations; it has been shown, in particular, that "gradient-like" behavior (convergence of all solutions) is usually preserved for small values of the parameter, as well as some characteristics of the solutions such as the number of slipped cycles [Smirnova et al., 2015b, Smirnova et al., 2015a, Smirnova et al., 2016, Smirnova and Proskurnikov, 2017.
The analysis of the aforementioned slow-fast dynamics effects in multi-agent networks remains, however, a difficult problem that is a subject of ongoing research.

Airplane with "Feathers" in a Turbulent Airflow
In this section, we give a detailed exposition of the problem formulated in [Granichin and Khantuleva, 2017].

Flight Stabilization in a Turbulence Region
The problem of reducing the aerodynamical forces influence force on the hull arises for an airplane that has fallen into a turbulence zone, where it is subject to rapidly changing forces and their moments that lead to potentially dangerous ride. Birds can partially control the influence of forces and moments with their feathers. Suppose we are able to construct an airplane with a large number of "feathers" that would smooth the influence of turbulent wind flows. From a traditional point of view we need a central fastest computer which performs very quickly complex computations collect based on collected data from all "feathers". In parallel with computation it feeds to all "feathers" computed control actions. This can be avoided if we take into account the possibility for self-organization between "feathers". For this purpose, at the first, we need to organize the local interaction between neighbor "feathers" by transmitting electric signals that by exchanging data regarding the forces acting on each "feather". At the second, we suppose that the "feathers" are able to turn in such a way that equalizes the forces between them induced by the turbulent flow. In the process of equating the induced forces, dynamic structures will form on the surface in the form of clusters with virtually identical values of induces forces. With time, these clusters will grow and merge with each other, reducing their number on the airplane's surface until they form a single cluster. By smoothing force fields acting on the airplane in a turbulence region, we can make its flow mode close to laminar. Let an airplane with the total mass M is flying at a velocity V(t) due to the sum of forces acting on it where the forces F e , F d , F g , F l are engine thrust, air drag force, gravity and lift respectively. If the regime of the flight is stationary all the forces are counterbalanced F e − F d − F g + F l = 0, the plane is moving along the straight line at a constant velocity V 0 = col(V 0 , 0, 0).
The gas flow over the plane considers being laminar and known for the given body's form. It means that the distribution of the forces over the body's surface is also known. When the airplane gets into the turbulent wind flow the forces distribution changes, the total balance of forces is broken, the flight becomes non-stationary and jolting Besides deviations from the given trajectory due to the force F 1 applied to the resultant force point r 1 (t), moment of the force M 1 (t) = [r 1 (t) × F 1 (t)] causes irregular rotation of the plane around its center of mass changing angular momentum of the plane L 1 (t) = [r 1 (t)×mV 1 (t)] according to the law dL1(t) dt = M 1 (t). Hereinafter [·×·] denotes the vector product. This rotation leads not only to deviations from a given trajectory but also to a ride due to inertia and rapid turns around the center of mass in different directions. To a large extent the irregular rotation of the plane is connected to the effect of jolting. It is impossible to completely avoid the influence of a turbulent flow, but it is desirable to make it less chaotic.
In order to reduce jolting it would be desirable to compensate the forces acting on the plane from the turbulent flow by additional forces. To achieve this goal the special constructive elements like the bird's feathers on the body's surface may be used. Assume that the surface of an airplane wing consists of N similar elements ("feathers") a 1 , . . . , a N (see. Fig 1) Here i ∈ N = {1, . . . , N } stands for the number of "feather". Denote r i is a vector from the center of mass of airplane to the center of "feather" a i . Each element is equipped with a pressure sensor (force-sensitive resistor), allowing to measure pressure force y i , and an actuator u i with two servomechanisms, allowing to change the air pressure. Each feather can rotate in two directions, parameterized by tilt angle α i ∈ [α − , α + ] in the vertical plane and by rotation angle β i ∈ [β − , β + ] in horizontal plane. The tilt and rotation angle determine, respectively, the vertical and cross components of the air drag force.
The integral air force acting on i-th element with area S has three components and three degrees of freedom (normal reaction of the plate and two angles α i and rotation β i ) are needed to control the feathers where n i (t) is the normal vector to i-th plate, C d (α i (t), β i (t)) is the aerodynamic coefficient of i-th plate defined by two angles α i (t), β i (t), v i (t) is the velocity magnitude of the flow over the i-th plate on the plane's surface at instant t. The force d i (t) has three projections: one component along x-axis defines air resistance, vertical component determines the lift and the third is lateral yaw component. Summarizing the activities of all the feathers one gets a total force acting on the plane and changing its trajectory. For rotations of agent i we can use micro-control action u i (t) The macro control action U i (t) for agent i is defined by airplane engine power. Suppose that airplane moves in a laminar stationary air flow from the initial instant t 0 till some instant t 1 , during this time period [t 0 , t 1 ] the integral air force acting on i-th element ("feathers") is constant d 0 i and all "feathers" are lying on the surface: α i = 0 and β i = 0 (see Fig. 2). Let x i (t) be the state vector of agent ("feathers") i. It consists of 11 components x (j) i (t), j = 1, 2, . . . , 5: 3 coordinates, 3 velocities, 3 rotations, and two angles. We assume that y i (t) = ∥d i (t)∥ without loss of generality. Let z i (t) = d i (t) − d 0 i be the deviation of the integral air force acting on i-th element from initial vector d 0 i in laminar wind flow. For the forces F e , F d , F g , F l we have in laminar wind flow where we introduce additional agent N + 1 located in the center of mass with r N +1 = 0 and m N +1 = M − ∑ i∈N θ In the laminar wind flow, when the plane is moving along the straight line at a constant velocity V 0 , dynamics of agent i is described by the differential equationṡ i (t) = 0, j = 2, 3, 6, . . . , 11, x (4) where constant e is determined by a constant engine power.
In turbulent flow "feathers" begin to rise and turn. Each feather is a plate on which the force is acting depending on its position relative to the flow direction. Fig. 3 shows the different forces (different colors) for different units of airplane in case of a turbulent wind flow when all feathers remain the initial (equal) orientations. The forces deviations z i (t) are induced by the laminar flow small perturbations during the plates' turns and by the turbulent flow around them. We consider "intellectual feather" which are able to arrange a self-regulation in the system. We denote set N i of feather i neighbors and assume that each feather i gets the information about d j (t) and angles α j (t), β j (t) for j ∈ {i} ∪ N i . This process should be much faster and can be repeated at any change of the turbulent flow. The system of such miniature intellectual agents has the ability to selforganization. As a result of the forces synchronization, groups of agents get approximately equal components of the forces deviations acting on agents and form clusters on the surface of the plane (see Fig. 4). Such structuring is provided by the goal-directed collective behavior of the system. Suppose that with respect to the airplane's axis of motion the wings are large in the direction perpendicular to the motion and small lengthwise, while the size of the airplane hull in the perpendicular direction is also small. Under these assumptions, among disturbing forces z i (t) with amplitudes of approximately the same order of magnitude acting on the entire surface of the airplane the only vertical and lateral projections of the forces deviations are important for the problem of reducing the jolting, and the total effect of all other moments is small. For the mentioned above disturbing force F 1 (t) we have We can use micro-control actions u i (t) to equalize deviation forces z i (t) if the changing of turbulent wind flow is sufficiently slow so that we have no any significant changing before instant t 2 . Assume that we can do it during a small time interval τ : t 1 + τ < t 2 , and z i (t) = z 1 . By virtue (39) we obtain that F 1 (t) = N z 1 = const during time interval [t 1 + τ, t 2 ]. Hence, we have for the rotation moment that For homogeneous distribution of feathers on the symmetrical body ∑ n i=1 r i = 0 and we have M 1 = 0. So, the most disturbing factor connected to jolting disappeared during the process of equalization of the forces acting on different feathers. When the forces deviations induced by the turbulent wind flow are equalized the total force correction on account of the feather's work tends to F 1 = const and the flow over the body becomes almost laminar again until the turbulent wind will not change. Because of the small and constant force F 1 the trajectory of the plane is slowly changing in one direction and can be corrected by the macro control U i . More precisely, during time interval [t 1 + τ, t 2 ] we can use the following model of agent dynamicṡ where model parameters θ i (t) are determined by force F 1 , mass, position, orientation, and initial parameters of motion in the laminar wind flow.

Consensus Algorithm
We are interested in synchronization of the agents with respect to the outputs z i (t), i = 1, . . . , N , introduced above. If it were possible to control directly the velocitiesż i , the vector-valued modification of the protocol (11) could be useḋ . . . , N. (41) Here B = (b ij ) is a non-negative weighted adjacency matrix of the communication graph, N i = {j : b ij ̸ = 0} is the set of neighboring nodes for node i and γ > 0 is the control gain, introduced for convenience. Algorithm (41) coincides with (11), where ε ij = γb ij . Henceforth we always assume that communication is bidirectional (b ij = b ji ) and connected, which, as has been discussed, guarantees the exponential convergence of (41) to consensus. We have also seen that (41) in the case of undirected graph is nothing else than an algorithm of gradient descent, applied to the "energy" function The algorithm (41), however, cannot be directly exploited in our situation since the actual control inputs are the ratesα i andβ i . Consider now the control algorithṁ Here γ > 0 is a control gain, and we assume thatα i = β i = 0 if y i = 0. A straightforward computation shows that the algorithm (43),(44) is equivalent to (41) and hence exponentially converges to consensus. More important, it is possible to estimate time needed to reach an approximate consensus in the sense that |Q(Z(t))| ≤ ε.
The estimate involves ε, γ, initial conditions (that is, the values of Z(t 1 +) and characteristics of the graph. The exact formulation is available in [Granichin and Khantuleva, 2017].

Conclusion and Future Work
The problem of determination of the constants θ i (t)) in equations (40) is not considered in this paper. This parameters remain constants during a finite time horizon. Traditional asymptotic methods cannot be used and we plan to apply for this case extended LSCR (Leave-out Sign-dominant Correlation Regions) approach [Amelin and Granichin, 2016], which was proposed earlier by M. Campi with co-authors to increase the effectiveness of adaptive control based on finite (not large) set of experimental data only [Campi et al., 2009, Campi andWeyer, 2010]. Another useful SPS (Sign-Perturbed Sums) method is proposed in [Care et al., 2017].