## Abstract

The control of complex systems and network-coupled dynamical systems is a topic of vital theoretical importance in mathematics and physics with a wide range of applications in engineering and various other sciences. Motivated by recent research into smart grid technologies, we study the control of synchronization and consider the important case of networks of coupled phase oscillators with nonlinear interactions—a paradigmatic example that has guided our understanding of self-organization for decades. We develop a method for control based on identifying and stabilizing problematic oscillators, resulting in a stable spectrum of eigenvalues, and in turn a linearly stable synchronized state. The amount of control, that is, number of oscillators, required to stabilize the network is primarily dictated by the coupling strength, dynamical heterogeneity, and mean degree of the network, and depends little on the structural heterogeneity of the network itself.

- control
- synchronization
- Complex Networks
- nonlinear dynamics
- power grid

## INTRODUCTION

Complex networks and complex systems describe the physical, biological, and social structures that connect our world and host the dynamical processes vital to our lives (*1*–*3*). The failure of such large-scale systems to operate in the desired way can thus lead to catastrophic events such as power outages (*4*, *5*), extinctions (*6*, *7*), and economic collapses (*8*, *9*). Thus, the development and design of efficient and effective control mechanisms for such systems is not only a question of theoretical interest to mathematicians but also has a wide range of important applications in physics, chemistry, biology, engineering, and the social sciences (*10*, *11*).

The roots of modern linear and nonlinear control reach back several decades, but recently, research in this direction has seen a revival in physics and engineering communities. For instance, the concept of “structural controllability,” which is based on the paradigm of linear homogeneous dynamical systems, was first introduced by Lin (*12*) and more recently investigated by Liu *et al*. (*13*) and Yuan *et al*. (*14*). These advances have enabled further progress related to structural controllability such as centrality (*15*), energy (*16*), effect of correlations (*17*), emergence of bimodality (*18*), transition and nonlocality (*19*), the specific role of individual nodes (*20*), target control (*21*), and control of edges in switchboard dynamics (*22*). Significant advances have also been made in the control of nonlinear systems, for instance, the control of chaotic systems using unstable periodic orbits (*23*), control via pinning (*24*–*26*), control and rescue of networks using compensatory perturbations (*27*, *28*), and control via structural adaptation (*29*). Implicit in all such network control problems are the following questions: (i) What form(s) of control should one choose? (ii) How much effort is needed to attain a desired state (*30*)?

Motivated by ongoing studies on the stability and function of power grids (*31*, *32*), we study the control of heterogeneous coupled oscillator networks (*33*, *34*). Recent research into smart grid technologies has shown that certain power grid networks called “microgrids” evolve and can be treated as networks of Kuramoto phase oscillators (*35*). A microgrid consists of a relatively small number of localized sources and loads that, while typically operating in connection to a larger central power grid, can disconnect itself and operate autonomously as may be necessitated by physical or economical constraints. In particular, by means of a method known as “frequency drooping,” the dynamics of microgrids become equivalent to Kuramoto oscillator networks—a class of system for which a large body of literature detailing various dynamical phenomena exists (*36*). Here, we develop a control mechanism for such coupled oscillator networks, thus providing a solution with potentially direct application to the control of certain power grids.

Our goal is to induce a synchronized state (that is, consensus) in a given coupled oscillator network and guarantee asymptotic stability by applying as few control gains to the network as possible. Our method is based on calculating the Jacobian of the desired synchronized state and studying its spectrum, by which we identify the oscillators in the network that contribute to unstable eigenvalues and thus destabilize the synchronized state. Our method not only identifies which oscillators require control but also the required strength of each control gain. We find that the control required to stabilize a network is dictated by the coupling strength, dynamical heterogeneity, and mean degree of the network, and depends little on the structural heterogeneity of the network. That is, the number of nodes requiring control depends surprisingly little on the network topology and degree distribution and is more sensitive to the average connectivity of the network and the dynamical parameters. Although Kuramoto oscillator networks serve as our primary system of interest due to both its specific correlation with microgrids and its rich body of literature, we note that our method can be applied to a much wider set of oscillator networks, provided that their linearized dynamics take a certain form. Moreover, because Kuramoto and other oscillator network models have served as a paradigmatic example for modeling and studying synchronization in various contexts, we hypothesize that our results may shed light more generally on the control of synchronization processes and could potentially give insight into other important applications such as the termination of cardiac arrhythmias (*37*) and treatments for pathological brain dynamics (*38*).

## RESULTS

### The Kuramoto model

We consider the famous Kuramoto model for the entrainment of many coupled dissipative oscillators (*39*). The Kuramoto model consists of *N* phase oscillators θ_{i} for *i* = 1, …, *N* that, when placed on a network dictating their pairwise interactions, evolve according to:(1)Each oscillator *i* has a unique nature frequency ω_{i} that describes its preferred angular velocity in the absence of interactions, which is typically drawn randomly from a distribution *g*(ω). Furthermore, the global coupling strength *K* describes the influence that oscillators have on one another via the network connectivity, which is encoded in the adjacency matrix [*A*_{ij}]. Here, we focus on the simple case of an undirected, unweighted network (*A*_{ij} = 1 if oscillators *i* and *j* are connected by a link and *A*_{ij} = 0 otherwise), but we note that all results presented here easily generalize to directed and weighted networks. We also assume that the network is connected, that is, irreducible. Over the last few decades, the Kuramoto model has proven to be very useful for modeling real-world systems (*36*, *40*), uncovering the mechanisms behind emergent collective behavior (*41*, *42*), exploring additional effects such as time delays (*43*) and community structure (*44*), and finding optimal network structure (*45*).

Depending on the coupling strength *K*, as well as the frequency vector ω and the network topology, the steady-state dynamics of Eq. 1 can attain many different states that included complete incoherence, partial synchronization, and full synchronization. The latter is characterized by and is also referred to as full phase locking, frequency synchronization, or consensus. The fully synchronized state (henceforth called the synchronized state) typically displays a large degree of phase synchronization *r* ≈ 1, where is the standard Kuramoto order parameter. In Fig. 1A, we illustrate a synchronized state in a group of five oscillators, each moving with an angular velocity of ω. The order parameter is illustrated as vector of length *r* with an offset angle ψ from the positive real axis.

### Stability and control

Here, we address the problem of control by first assuming that because of the system parameters (that is, the coupling strength, natural frequency sequence, or network topology), the steady-state dynamics of Eq. 1 are at least partially incoherent, that is, one or more oscillators remain desynchronized. In the example of the power grid, a single desynchronized oscillator represents a single power failure but can have further damaging effects, in particular, triggering a cascade of additional failures and ultimately a power outage (*5*). Thus, our goal is to find a synchronized state and stabilize it. If all oscillators are initially synchronized, then our goal is trivially realized; however, our method can be used to make this state more robustly stable. In a synchronized state such as the one we seek here, we expect the oscillators to be clustered in a single reasonably tight cluster such that |θ_{j} − θ_{i}| ≪ 1, and thus, Eq. 1 can be linearized to:(2)where *L* is the network Laplacian matrix with entries *L*_{ij} = δ_{ij} ∑_{l} *A*_{il} − *A*_{ij}, and δ_{ij} is the Kronecker delta. A straightforward analysis yields a “target” synchronized state (within the rotating reference frame θ ↦ θ + 〈ω〉*t*) given by the vector θ* = *K*^{−1}*L*^{†}ω, where *L*^{†} is the pseudo-inverse of the Laplacian (*46*). (We summarize a derivation of this result in Materials and Methods.) We note that, because the system is assumed to be partially incoherent, the fixed point θ = θ* either does not exist or is unstable. However, we take θ = θ* to represent the closest synchronized fixed point for the given parameter values, and therefore, we use it as a target. We also note that, although Eq. 2 was directly obtained from linearizing Eq. 1, other systems of more general forms yield equivalent linearizations and therefore can also be controlled using the method we provide here. We present an example of such a general system with arbitrary coupling function [for example, see (*47*)] in Materials and Methods.

The stability of θ = θ* is dictated by the Jacobian matrix whose entries are defined , and is stable if all the eigenvalues of *DF*|_{θ*} are nonpositive. In our case, we have that:(3)

We note that each row (and column) of *DF* sums to zero, that is, satisfies *DF*_{ii} = −∑_{j≠i} *DF*_{ij}. This is a particularly convenient property for using the Gershgorin circle theorem (*48*), which implies that the eigenvalues of *DF* lie within the union of closed discs *D*_{i} for *i* = 1, …, *N*, which are each centered at *DF*_{ii} and have radius *R*_{i}, where *R*_{i} = ∑_{j≠i}|*DF*_{ij}|. (The full theorem is given in Materials and Methods.) In particular, if all the off-diagonal entries of *DF* are nonnegative, then it follows that each Gershgorin disc is contained in the left-half plane, implying that all eigenvalues are nonpositive and the solution is stable. An illustration of this case is presented in Fig. 1B. If, however, one or more nondiagonal entries of *DF* are negative, then each Gershgorin disc corresponding to a row with a negative off-diagonal entry enters the right-half plane, admitting the possibility for one or more positive eigenvalues and thus destabilization. Thus, the oscillators that require control can be easily identified as those whose corresponding rows have one or more negative off-diagonal entries.

We aim to stabilize the synchronized solution by adding one or more control gains to the system, as illustrated in Fig. 1C. In the following recent literature, we will refer to oscillators to which we apply control as “driver nodes,” and to oscillators to which we do not apply control as “free nodes.” We choose the control gains to take the form *f*_{i}(*t*) = *F*_{i} sin(ϕ_{i} − θ_{i}) where *F*_{i} is the strength of the *i*th control gain and ϕ_{i} is a target phase that can, in principle, depend on either local or global information, and vary in time. Here, we focus on the choice of target phase and discuss other possibilities below. Because the control gain depends on the current state of the system, this can be thought of as a form of feedback control. The new dynamics are then given by:(4)where we take *F*_{i} = 0 for free nodes. Whereas the off-diagonal entries of *DF* remain unaltered, the new diagonal entries are given by . Thus, we set coupling gain strength of each driver node *i* such that it satisfies . This ensures that all Gershgorin discs are contained in the left-half plane, implying that (up to the linear approximation of θ*) all eigenvalues are nonpositive and the synchronized state is stable. (In the case of a directed network, this implies that all eigenvalues have nonpositive real part and the synchronized state is stable.)

We now briefly comment on the choice of target phases ϕ_{i} in the control gains. In the method outlined above, we have set the target phase equal to the steady-state phase, . This is a convenient choice for the derivation described above. Additionally, we find that, in practice, other choices also yield positive results. In particular, one choice that tends to yield slightly better results is to force each driver node toward the center of the synchronized cluster, that is, , where we assume that the cluster is centered at the angle . Target phases can also be chosen according to global or distributed control strategies. In particular, given the standard Kuramoto order parameter or the set of local order parameters , the choices ϕ_{i} = ψ and ϕ_{i} = ψ_{i} correspond to global and distributed control strategies that typically yield favorable results. We also note before presenting examples that because the method outlined above depends on the approximation of the steady-state solution θ* ≈ *K*^{−1} *L*^{†}ω, in practice, we add a buffer margin when identifying unstable oscillators, looking for nondiagonal entries of *DF* that are not necessarily negative, but rather less than some ϵ > 0. We find that the choice ϵ = 0.2 is sufficient, and is what we use in the examples below.

### Control of random networks

We now demonstrate our approach by considering two types of random networks: Erdős-Rényi (ER) (*49*) networks and scale-free (SF) networks. Each ER network is constructed using a fixed link probability *p*, and each SF network is built using the configuration model (*50*) with a degree sequence drawn from the distribution *P*(*k*) ∝ *k*^{−γ} with γ = 3 and enforced minimum degree *k*_{0}. To tune the mean degree 〈*k*〉 of each network, we set either *p* = 〈*k*〉/(*N* − 1) or *k*_{0} = 〈*k*〉/(γ − 1). Figure 2 illustrates our results with an example of each type of network, where we have used networks of size *N* = 1000 with mean degree 〈*k*〉 = 6, set the coupling strength *K* = 0.4, and used natural frequencies drawn from a uniform distribution with zero mean and unit variance. The top and bottom rows display the results for the ER and SF networks, respectively, displaying the time series θ(*t*) of a randomly selected 10% of the oscillators without control (Fig. 2, A and D) and with control (Fig. 2, B and E), after discarding a long transient. The difference between no control to control is quite drastic, with a large fraction of desynchronized oscillators without control, and full synchronization with control. Driver nodes constituted 37.1 and 44.5% of the ER and SF networks, respectively, to attain the synchronized state. In Fig. 2 (C and F), we present the degree of phase synchronization, plotting the order parameter *r*(*t*) for both solutions with and without control. Unlike the solution without control, which fluctuates significantly at a relative low value, the solution with control reaches a steady value near *r*(*t*) ≈ 1.

Next, we investigate the properties of driver and free nodes by revisiting our method. This is an essential question because of the heterogeneity of the oscillators in terms of both network structure (that is, degree distribution) and local dynamics (that is, natural frequencies). An oscillator *i* is a driver node if one of its off-diagonal entries is negative, and therefore, oscillators with large (small) steady-state values tend to be driver (free) nodes. Furthermore, we find that these values scale approximately linearly with the ratio of the natural frequencies to degrees, that is, . We illustrate this in Fig. 3, where we plot the relationship [*L*^{†}ω]_{i} versus ω_{i}/*k*_{i} for the example ER and SF networks presented above (Fig. 3, A and B, respectively), denoting driver nodes with red crosses and free nodes with blue circles. These results show the important role that dynamics, in addition to network structure, plays in dictating controlling the system. In particular, driver nodes of the system tend to balance a large ratio (in absolute value) of natural frequencies to degrees.

Finally, we quantify the overall effort required for consensus by studying how the fraction of driver nodes, denoted *n*_{D} = *N*_{D}/*N*, where *N*_{D} is the total number of driver nodes, depends on both the system’s dynamical and structural parameters. Presenting our results in Fig. 4, we first explore how the fraction of driver nodes depends on the coupling strength by plotting in Fig. 4A *n*_{D} versus *K* for both ER and SF networks with mean degrees 〈*k*〉 = 4, 8, and 12 (blue circles, red triangles, and green squares, respectively). Results for ER and SF networks are plotted with unfilled and filled symbols, respectively, and each curve represents an average of over 100 network realizations, each averaged over 100 random natural frequency realizations. Although it is expected that *n*_{D} decreases monotonically with *K*, the curves’ dependence on network topology and mean degree is nontrivial. In particular, the shape of *n*_{D} versus *K* depends more sensitively on the mean degree than the topology, suggesting that network heterogeneity has little effect on overall control in comparison to average connectivity. In light of the significant dependence of overall control on the coupling strength, we investigate the coupling strength required to synchronize a network if a limited amount of control is available. To this end, we calculate for each family of networks the required coupling strengths *K*_{5%}, *K*_{10%}, and *K*_{20%} for which, on average, a fraction *n*_{D} = 0.05, 0.1, and 0.2 will achieve synchronization as a function of the average degree 〈*k*〉. We plot the results in Fig. 4B. We point out again that ER and SF networks behave very similarly on average, and that with a larger mean degree, a smaller coupling strength is required to achieve synchronization.

## DISCUSSION

The theoretical and practical aspects of the control of dynamical processes remain important and ongoing areas of interdisciplinary research at the intersection between mathematics, physics, biology, chemistry, engineering, and the social sciences. The control of complex networks and complex systems is particularly important because, together, they comprise most of the world we live in (*51*); however, the nonlinear nature of realistic dynamical processes and the complex network topologies of real networks represent challenges for the scientific community. Building on concepts from classical linear control theory (*12*), recent work has made significant advances in understanding structural controllability (*13*, *14*), and significant progress has been made in the development of control mechanisms for networks of nonlinear systems (*23*, *24*, *28*). Nonetheless, because of the problem-sensitive nature of most real-world problems and applications requiring control techniques, further progress in designing and implementing efficient and effective control mechanisms for a wide range of problems with practical constraints remains an important avenue of research.

Here, we have focused on the control of synchronization (that is, consensus) in coupled oscillator networks. Our primary inspiration has been advances in the research of power grid networks (*52*, *53*). In particular, recent studies have shown that certain power grids known as microgrids can be treated as Kuramoto oscillator networks (*35*, *39*). Here, we have presented a control method that can easily be applied to Kuramoto networks and other phase oscillator networks, thus providing a control framework with potentially direct application to these new technologies. Our method is based on identifying and stabilizing a synchronized state for a given network via spectral properties of the Jacobian matrix, and we have demonstrated its effectiveness on both ER and SF networks. We have observed that driver nodes, that is, oscillators that require control, tend to balance (in absolute value) large natural frequencies with small degrees. Furthermore, the overall amount of control required to achieve synchronization decreases with both coupling strength and mean degree, whereas the total effort required to attain a synchronized state depends sensitively on the average connectivity of the network and the dynamical parameters, but surprisingly little on the network topology and degree distribution. These results enhance our understanding of and ability to understand, optimize, and ultimately control synchronization in power grid networks [see, in particular, (*5*, *32*)], and more generally complement important work on the control of network-coupled nonlinear dynamical systems (*26*, *28*, *29*).

Although our central inspiration and target application are in the area of power grid technology, synchronization phenomenon plays a vital role in a variety of complex processes that occur in both natural and manmade systems, including healthy cardiac behavior (*54*), functionality of cell circuits (*55*), stability of pedestrian bridges (*56*), and communications security (*57*). Given this broad range of applications, we hypothesize that our findings here may potentially shed some light on the control of synchronization in other contexts, such as cardiac physiology and neuroscience. For instance, a large amount of research has recently been devoted to the development of cardiac arrhythmia treatments that require minimal shock to knock out fatal asynchronous behavior such as cardiac fibrillation (*58*) and the promotion of normal brain oscillations (*59*) while repressing disorders such as Parkinson’s disease, which are associated with abnormal oscillations (*60*).

## MATERIALS AND METHODS

### Steady-state solution

To derive the steady-state solution θ* = *K*^{−1}*L*^{†}ω, we begin with Eq. 2, which represents the linearized dynamics of Eq. 1. Recall that this linearization requires that we are searching for a synchronized state where all oscillators are tightly packed in a single cluster, so we expect that |θ_{j} − θ_{i}| ≪ 1. We also note that the mean frequency of all oscillators is given by the mean natural frequency 〈ω〉. For simplicity, we enter the rotating frame θ ↦ θ + 〈ω〉*t*, effectively setting the mean frequency to zero. It is then convenient to write Eq. 2 in vector form, that is:(5)where *L* is the network Laplacian whose entries are defined *L*_{ij} = δ_{ij} ∑_{l} *A*_{il} − *A*_{ij}. Although *L* has a zero eigenvalue, denoted λ_{1} = 0, rendering it noninvertible, it does have a pseudoinverse defined using its other eigenvalues (which are nonzero provided that the network is connected) and corresponding eigenvectors, (*46*). Each eigenvector is normalized such that forms an orthonormal basis for the space of vectors in ℝ^{N} with zero mean. Thus, both *L* and *L*^{†} share a null space, which is spanned by the eigenvector *v*^{1} ∝ 1, and therefore map vectors onto the space of zero-mean vectors in ℝ^{N}. With the pseudoinverse in hand, we can finally obtain the desired steady-state solution by setting θ = 0 and solving for θ, which yields the solution θ* = *K*^{−1}*L*^{†}ω, as desired.

### General oscillator networks

Here, we present an example of a more general oscillator network than that in Eq. 1 that can be controlled using the same method detailed above. In particular, we generalize to account for an arbitrary coupling function *H*(θ), yielding:(6)We assume that *H*(θ) is 2π periodic and at least once continuously differentiable. *H* need not satisfy *H*(0) = 0, and thus, coupling between neighboring oscillators can be “frustrated” (*47*), denoting that even when two oscillators are exactly equal, their interaction term does not vanish. Provided that the coupling frustration is not too large, for example, , a tightly clustered synchronized state is attainable, and linearizing Eq. 6 yields:(7)By defining the quantities and , it is easy to see that the linearized dynamics of Eq. 7 are of the same form as Eq. 2, and therefore, the control method we present above can be readily applied.

### Gershgorin circle theorem

Definition. (Gershgorin discs) Let *M* be an *N* × *N* complex matrix. For *i* = 1, …, *N*, let *R*_{i} = ∑_{j≠i} |*M*_{ij}| be the sum of absolute values of nondiagonal elements of row *i*, and define *D*(*M*_{ii}, *R*_{i}) closed disc of radius *R*_{i} centered at *M*_{ii}. *D*_{i} = *D*(*M*_{ii}, *R*_{i}) is the *i*th Gershgorin disc.

Theorem. (Gershgorin) All eigenvalues of the matrix *M* lie within the union of Gershgorin discs.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Funding:**This work was supported by the James S. McDonnell Foundation (P.S.S. and A.A.), Spanish DGICYT grant FIS2012-38266 (A.A.), and FET Project No. MULTIPLEX (317532) (A.A.).

**Author contributions:**P.S.S. and A.A. conceived the study. P.S.S. performed the numerical experiments and analyzed the results. P.S.S. and A.A. wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

- Copyright © 2015, The Authors