Optimal cost tuning of frustration: Achieving desired states in the Kuramoto-Sakaguchi model
OOptimal cost tuning of frustration: Achieving desired states in theKuramoto-Sakaguchi model
Gemma Rosell-Tarrag´o ∗ and Albert D´ıaz-Guilera † Departament de F´ısica de la Mat`eria Condensada,Universitat de Barcelona, Mart´ı Franqu`es 1, 08028 Barcelona, Spain andUniversitat de Barcelona Institute of Complex Systems (UBICS), Mart´ı Franqu`es, 1, 08028 Barcelona, Spain (Dated: January 29, 2021)There are numerous examples of studied real-world systems that can be described as dynamicalsystems characterized by individual phases and coupled in a network like structure. Within theframework of oscillatory models, much attention has been devoted to the Kuramoto model, whichconsiders a collection of oscillators interacting through a sinus function of the phase differences.In this paper, we draw on an extension of the Kuramoto model, called the Kuramoto-Sakaguchimodel, which adds a phase lag parameter to each node. We construct a general formalism thatallows to compute the set of lag parameters that may lead to any phase configuration within alinear approximation. In particular, we devote special attention to the cases of full synchronizationand symmetric configurations. We show that the set of natural frequencies, phase lag parametersand phases at the steady state is coupled by an equation and a continuous spectra of solutions isfeasible. In order to quantify the system’s strain to achieve that particular configuration, we definea cost function and compute the optimal set of parameters that minimizes it. Despite considering alinear approximation of the model, we show that the obtained tuned parameters for the case of fullsynchronization enhance frequency synchronization in the nonlinear model as well.
I. INTRODUCTION
Emergence is one of the key concepts in the analysisof complex systems [1]. Collective properties emerge as aconsequence of irregular interactions among its elementalconstituents [2]. One of the most paradigmatic examplesof emergence is synchronization [3, 4], because the inter-play between populations of oscillatory units gives riseto a variety of global states, ranging from perfect syn-chronization or phase locked stationary configurations tochimera states [5–7]. Among the different models thathave been used to understand such collective behavior,a lot of effort has been devoted to the Kuramoto model(KM), in which phase oscillators interact continuouslywith other units through a sine function of the phasedifference [8–10].In the past few years there has been a growing inter-est in the concept of controllability, which quantifies thefeasibility to achieve a desired final state of a given dy-namical system [11]. As stated above, the KM can giverise to a wide variety of stationary (phase or frequencysynchronized) or not stationary states, being chimeras anunexpected mixture of both types of behaviors [12]. Inthis context, controllability can be understood as a tun-ing of the internal parameters of the oscillators to reachspecific phase configurations. The most simple settingsstand for a collection of identical oscillators interactingthrough a sinus function of the phase differences. In thiscase it is quite intuitive to see that the final state is a per-fectly synchronized one in which all oscillators have ex-actly the same phase and frequency (the same frequency ∗ [email protected] † [email protected] than the intrinsic one). It is the existence of a distribu-tion of frequencies that gives rise to a transition, in termsof the strength of the coupling, from an incoherent stateto a coherent one [13]; such a transition is robust in thesense that the introduction of a lag term, a phase addedto the argument of the sinus function, does not changethe behavior, as far as it is kept below π/ a r X i v : . [ phy s i c s . s o c - ph ] J a n actors for the network, because they enable more vari-ability in the states of the system, but also potentiallydangerous ones, as tiny perturbations can produce cas-cade like effects that completely changes the network dy-namics.As stated, the addition of a phase lag parameter en-ables a richer configuration state. However, it is clearthat a tuning of a single parameter will not be enough togenerate the wide variety of stationary states that a pop-ulation of Kuramoto oscillators can achieve. Notwith-standing, the question that arises is whether a fine tuningof a set of individual parameters can make it possible. Inthis paper, this is our proposal, we construct a generalformalism that allows, within a linear approximation, tocompute the set of lag parameters that may lead to anyphase configuration for a fixed set of intrinsic frequencies.The problem can also be posed the other way around.Namely, given a set of frequencies, we may derive theconfiguration of phases that is produced by a set of lagparameters.There are numerous examples of real-world systemsthat can be described as dynamical systems characterizedby individual phases and which functioning are object ofinvestigation. Some examples are the brain functionalnetworks arising from temporal correlation patterns, acpower in power grids[19], heartbeats[20], multiprocessorsand multicore processors, or traffic signaling. Not onlythe synchronization between their constituents may beintended or prevented, but also other particular configu-rations may be of relevant interest. For this reason, wepropose a mechanism for tuning the intrinsic parametersof the system to achieve any desired phase configuration.A previous work proposes a methodology to enhancefrequency synchronization for the nonlinear Kuramoto-Sakaguchi model (extension of the Kuramoto model witha node phase lag parameter) [21]. Another work sug-gests that an unstable synchronized state becomes sta-ble when, and only when, the oscillator parameters aretuned to nonidentical values [15]. We highlight the workdone in Reference [22], where the particular configurationof perfect synchronization is studied and the synchronyalignment function is defined in order to minimize the or-der parameter of the system considering different topolo-gies and frequency scenarios. We address the most gen-eral question, following a similar path to that pursued bythem, forcing the system to achieve any particular con-figuration for the linear case of the Kuramoto-Sakaguchimodel by means of a fine tuning of the phase lag or frus-tration parameter set. Despite considering a linear ap-proximation of the model, we show that the obtainedtuned parameters for the case of full synchronization en-hance frequency synchronization in the nonlinear modelas well.The structure of the paper is the following. In SectionII we present the Kuramoto-Sakaguchi model (a variationof the Kuramoto model) as the proper framework for ourpurposes. Next, in Section III, we derive the analytic ex-pression of the fine tuning of the frustration parameters so as to achieve any phase configuration. Then, in Sec-tion IV, we define a cost function to assess the expenseof achieving a particular phase configuration by its corre-sponding tuning and derive the analytic solution for thecases of symmetric and fully synchronized configurations,in Sections V and VI, respectively, as well as commenton the nonlinear validity of our results. We concludein Section VII. An easy-to-follow example and furthermathematical derivations can be found in the Appendix. II. THE KURAMOTO-SAKAGUCHI MODEL
In 1975, Kuramoto suggested one of the best-knowndynamical equation to model interacting oscillatorysystems[8]: a set of N phase oscillators characterized bytheir phase, θ i , and coupled between each other by thesine of their phase differences. Each unit is influenceddirectly by the set of its nearest neighbours via the adja-cency matrix of the network corresponding to the system, G ( V , E ). The coupling strength describes the intensity ofsuch pair-wise interactions, K ij . The set of nodes of thenetwork, V ( G ) consists of all the oscillators, while the setof edges, E ( G ), is made of the links between them. In hisoriginal work, Kuramoto assumed homogeneous interac-tions, i.e, K ij = K ∀ ( i, j ) [8, 10]. Taking into accountthe connectivity or topology of the network and the os-cillatory dynamics, the dynamics of the system can bewritten as a system of differential equations[5]: dθ i dt = ω i + K (cid:88) j A ij sin( θ j − θ i ) i = 1 , ..., N j ∈ Γ i (1)where Γ i is the set of neighbors of node i and ω i is thenatural frequency of each unit.We consider that two nodes are phase synchronizedwhen their phases have the same value, θ i ( t ) − θ j ( t ) = 0 ∀ t > t When the phase difference has a constant value, that is, θ i ( t ) − θ j ( t ) = c ∀ t > t , we say there is a phase lockingbetween nodes i and j . Similarly, we consider that twonodes are frequency synchronized when their frequencieshave the same value: dθ i dt − dθ j dt = 0 ∀ t > t We say that two nodes are fully synchronized when theyare phase synchronized, because this implies frequencysynchronization.As long as the distribution of natural frequencies ishomogeneous, namely, all units have the same naturalfrequency, there is only one attractor of the dynamics:the fully synchronized state. It can be shown that, ifthe distribution of natural frequencies is unimodal, thesystem becomes frequency synchronized as long as thecoupling strength is larger than a threshold value[10].In 1986, Kuramoto together with Sakaguchi presenteda similar model which incorporated a constant phase lagbetween oscillators[9] which can be written as follows: dθ i dt = ω i + K (cid:88) j A ij sin( θ j − θ i − α ) i = 1 , ..., N j ∈ Γ i (2)where α is a homogeneous phase lag parameter. Thismodel has also become well-known and several variationsof the model have been studied. It has been shown that,as long as | α | < π/
2, the system is not chaotic and athreshold value for the coupling strength exists abovewhich the system becomes synchronized to a resultingfrequency [9]. In the particular case that considers ho-mogeneous natural frequencies, i.e, ω i = ω ∀ i , the frus-tration parameter, α , forces the system to break the oth-erwise original fully synchronized state. However, partialsynchronization is conserved for symmetric nodes in thenetwork[14, 15]. As the frustration increases the asyn-chronous groups’ phase move away from each other.We are interested here in a more general case, wherethe frustration parameter is not homogeneous but an in-trinsic property of each unit: dθ i dt = ω i + K (cid:88) j A ij sin( θ j − θ i − α i ) i = 1 , ..., N j ∈ Γ i (3)In this context, a recent work studies a particular effect ofthis frustration parameter by defining functionability, anew centrality measure of the nodes in a network, in orderto address the issue of which nodes, when perturbed,move the system from a synchronized state to one thatis more asynchronous in the sense that it enhances thephase differences between all pairs of oscillators [18]. III. ANALYTIC EXPRESSION OF THEFRUSTRATION PARAMETERS TUNING
We address the most general problem, which considersthe same dynamics as in Eq.(3), while allowing the edgesof the network to be weighted, a more realistic scenariofor real-world networks.For small values of the frustration parameters andphases close to each other, which is the case in frequencysynchronization, we can linearize Eq.(3) as follows: dθ i dt = ω i + K (cid:88) j W ij ( θ j − θ i − α i ) == ω i − K (cid:88) j L ij θ j − Kα i s i (4)where W ij is the value of the weight of the edge betweennode i and node j , s i ≡ (cid:80) j W ij is the weighted degreeof the i th node and L is the weighted Laplacian matrixdefined as L ij ≡ δ ij s i − W ij . In the stable regime, asynchronized frequency is achieved and, for all oscillators ˙ θ i = Ω. We can derive the value of the common frequencyoscillation, Ω, summing Eq.(4) over i : (cid:88) i Ω = (cid:88) i ω i − K (cid:88) i (cid:88) j L ij θ ∗ j − K (cid:88) i α i s i (5)Taking into account the steady state ˙ θ i = Ω ∀ i and ar-ranging summations: N Ω = (cid:88) i ω i − K (cid:88) j θ ∗ j (cid:88) i L ij − K (cid:88) i α i s i . (6)and finally, Ω = (cid:104) ω (cid:105) − K (cid:104) αs (cid:105) . (7)where we have used the Laplacian matrix property: (cid:80) i L ij = 0 and defined the averages (cid:80) i α i s i /N = (cid:104) αs (cid:105) and (cid:80) i ω i /N = (cid:104) ω (cid:105) . Now we can plug expression Eq.(7)to Eq.(4) to get the stable phases of oscillators, θ ∗ i : (cid:88) j L ij θ ∗ j = ω i K − (cid:104) ω (cid:105) K + (cid:104) αs (cid:105) − α i s i ∀ i (8)The solution of Eq.(8) regarding phases is undetermineddue to the singular nature of the Laplacian matrix.Hence, Eq.(8) is, in general, an undetermined system oflinear equations, that is, there is one free phase, which weshould use as a reference value for the solution. Nonethe-less, we do not work directly with the functional form ofphases because they are time dependent { θ ∗ i } = f i ( t ),but with the phase differences with respect to a refer-ence node, once the stationary state is achieved, φ i ≡ θ i − θ R (9)In this way, we work with time independent values. Inthis situation, φ R = 0, by definition, as φ R ≡ θ R − θ R = 0.On the other hand, the contribution (cid:104) αs (cid:105) − α i s i of theright-hand side of Eq.(8) can be written in matrix formas: − N − N − N − N ... − N N − N − N ...... ... ... ... − N − N ... N − N · s ... s ... ... ... ... s N − ·· α α ...α N − = ( − M · D s ) (cid:126)α where we have defined: M ≡ N − N − N − N ... − N N − N − N ...... ... ... ... − N − N ... N − N and D s ≡ s ... s ... ... ... ... s N − . We write Eq.(8) in matrixas L(cid:126)θ ∗ = 1 K (cid:126) ∆ ω − M · D s (cid:126)α (10)where ∆ ω i ≡ ω i − (cid:104) ω (cid:105) . Finally, we obtain the set ofunknowns { α i } : M · D s (cid:126)α = 1 K (cid:126) ∆ ω − L(cid:126)θ ∗ (11)Equation (11), however, does not have a solution, be-cause of the singular nature of M · D s matrix. M matrixis singular too, and hence, its inverse matrix does not ex-ist. Mathematically, det ( M · D s ) = det ( M ) · det ( D s ) = 0Similarly as we did for phases [18], we solve the singu-larity problem by setting a reference node, which we call control node, regarding frustration parameters, i.e., wewould not obtain the value for each of the parameters,but a relation between them: κ i ≡ α i − α C (12)where α C is the value of the control node. In this situa-tion, κ C = 0, by definition, as κ C ≡ α C − α C = 0.To easily write the matrix expressions, we define theselection matrix J ( n,m ) , which is, in general, an ( N − × ( N −
1) identity matrix after the removal of the mth rowand the nth column.
L(cid:126)θ ∗ turns to ˜ L ( k, R ) (cid:126)φ ∗ , where we have removed the kth row and the Rth column. The result does not depend onwhich row we remove, hence we can choose any k . Usingthe selection matrix, ˜ L ( k, R ) = J ( ,k ) · L · J ( R, ) ≡ ˜ L .Similarly, (cid:126) ˜ φ ( k ) = J ( ,k ) · (cid:126)φ ≡ (cid:126) ˜ φ , where we have removedthe kth row.In an equivalent way as the definition of the reducedLaplacian:˜ M D s ( k, C ) = J ( ,k ) · M D s · J ( C, ) ≡ ˜ M D s where ˜ M D s is M D s without the kth row and the Cth column.Similarly, (cid:126) ˜ κ ( k ) = J ( ,k ) · (cid:126)κ ≡ (cid:126) ˜ κ and (cid:126) ˜∆ ω ( k ) = J ( ,k ) · (cid:126) ∆ ω ≡ (cid:126) ˜∆ ω , where we have removed the kth row.Considering all the previous definitions and remarks,Eq.(10) can be rewritten as:˜ M D s (cid:126) ˜ κ = 1 K (cid:126) ˜∆ ω − ˜ L(cid:126) ˜ φ ∗ − α C · J ( ,k ) → (cid:88) j [ M D s ] ij (13)and finally, (cid:126) ˜ κ = (cid:16) ˜ M D s (cid:17) − (cid:18) K (cid:126) ˜∆ ω − ˜ L(cid:126) ˜ φ ∗ − α C · ˜ M(cid:126)s (cid:19) (14)where we have used J ( ,k ) (cid:126) (cid:80) j [ M D s ] ij = ˜ M(cid:126)s . Notice that
M D s matrix is singular, but the row sum is not zero,although it is so for the column sum. Hence, we need toset α C = 0 if we want to avoid extra constant arrays inthe final expression. In this particular case: (cid:126) ˜ κ = (cid:16) ˜ M D s (cid:17) − (cid:18) K (cid:126) ˜∆ ω − ˜ L(cid:126)φ ∗ (cid:19) (15)and keep in mind that κ C = 0. The obtained values of (cid:126)α depend on both the chosencontrol node, C , and the value we set for its frustrationparameter, α C . Notice, therefore, that there is a contin-uous spectrum of values for the frustration parameter inorder to achieve a particular phase configuration.Moreover and more importantly, due to the non-row-sum equal to zero of M D s matrix, the differences betweenthe obtained values are dependent of the control nodechoice. Mathematically, α i − α j ( C = l ) (cid:54) = α i − α j ( C = k )if l (cid:54) = k . This property will lead us to the definition ofa cost for the system to move to the final configuration,which will depend on both the control node and the valueof its frustration parameter.We provide an example of a toy network for the caseof a homogeneous natural frequencies distribution, i.e., ω i = ω ∀ i . In this case, Eq.(14) turns to: (cid:126) ˜ κ = (cid:16) ˜ M D s (cid:17) − (cid:16) − ˜ L(cid:126) ˜ φ ∗ − α C · ˜ M(cid:126)s (cid:17) which in the case of the network depicted in Figure 1,
012 3 456
FIG. 1. Network of seven nodes. leads to the solution κ κ κ κ κ κ = − α +3 φ + φ +3 φ φ − φ )22 φ − φ − φ + φ φ − φ + φ − φ + φ φ − φ + φ − φ − φ φ φ + φ − φ (16)where we have chosen κ = 0 and φ = 0. Hence, theresults are written as a function of the value α and φ i i (cid:54) = 0. Therefore, we can achieve any phase config-uration, given by the set { φ i } by tuning the frustrationparameters set { α } , where α i = κ i + α C .To illustrate how we obtain the final values, let us con-sider the following phase configuration: (cid:126) ˜ φ ( R =0) = (0 . , . , . , − . , − . , .
0) (17)In the general case where α C = α (cid:54) = 0: (cid:126) ˜ κ ( C =1) = (0 . − α , − . , − . , . , . , − . . If we choose α C = 0, then α i = κ i , we can include thevalue of the control node C = 1: (cid:126) ˜ α = (0 . , . , − . , − . , . , . , − . . Alternatively, we can choose whatever value we need re-garding the control node. For instance, if α C = α = 0 . (cid:126) ˜ α = (0 . , . , − . , − . , . , . , . α ’s, up to an error. For this lastexample and using the frustration parameters obtainedby setting α = 0 .
1, the nonlinear model leads to finalphases vector (cid:126) ˜ φ ( R =0) = (0 . , . , . , − . , − . , . ∼ .
3% of relative error with respect tothe initial Eq.(17). See the full derivation of the analyt-ical solution in Appendix A.
IV. OPTIMAL COST TUNING OFFRUSTRATION
As pointed out in Section III, there is a continuousspectrum of values for the choice of the frustration pa-rameters that enables the system access a particularphase configuration. The following question arises natu-rally: Among all the possible solutions, which is the onethat makes the system achieve a particular phase config-uration with the minimum required cost?This question is of particular relevance when we con-sider the plausible real nature of the system. If a realsystem needs to access a particular phase configuration,which may be associated with a precise function, it willtend to minimize the effort or cost to do so.In order to quantify the required cost, we define it asfollows: e T ( C ) ≡ (cid:88) i | α i ( C ) | (19)Henceforth, the cost associated to each node is given bythe absolute value of the required frustration parameter.The absolute value operator allows for a sign-free contri-bution of each node, a very convenient choice in the case that the system is not beforehand specified, and a gen-eral definition is proposed instead. Furthermore, unlikeother nonlinear cost functions such as the square sum ofthe parameters, no extra weight is given to larger values,besides the corresponding to a linear function.As previously remarked, e T ( C ) will depend both on thechosen control node, C , as well as the particular choiceof its frustration parameter, α C .The optimal configuration is given by the solution ofthe minimization problemmin C,x e T ( C, x ) = min
C,x N (cid:88) i α i ( C, x ) (20)where the x variable is not yet specified. Depending onthe problem we are interested in we would set it eitherto ω i , s i or any other combination of the parameters ofthe model. The minimal value of the cost will depend onthe proper choice of the control node, C , in addition ofthe particular value of its frustration parameter, α C , asthe free parameter left to be set. In Sections V and VIwe provide a thorough analysis of it.The cost required to achieve a particular phase config-uration depends on that configuration, the control nodeand the chosen value of α C . In Figure 2 we present anexample, following with the network presented in SectionIII and choosing different values of α C , we compute nu-merically the values of the required cost using Eq.(19) toachieve the phase configuration given in Eq.(17). Noticethat the global minimum depends on the control nodeand its frustration parameter. FIG. 2. (Color online) Implied cost to achieve the phase con-figuration in Eq.(17) as a function of the chosen control node, C , for the network in Figure 1 and considering five differentvalues of α C . Notice that the minimum cost is given, in thiscase, by α C = 0 and C = 1 or C = 5. In Section III we have derived the general analyticalsolution of the frustration parameters as a function ofa particular choice for the phase configuration. In thissection we have defined a cost function in order to assessthe optimal choice of such configuration.Depending on the phase configuration one is interestedin achieving, results will vary and the analytical expres-sions will have different features.In the following sections we will focus on two particularconfigurations, due to its intrinsic importance, in order toobtain and discuss the analytical solution of Eq.(20): Theconfiguration given by the symmetries of the network[14]and the fully synchronized state.
V. SYMMETRIC PHASE CONFIGURATION
As explained in Section II, a particular example of theKuramoto-Sakaguchi model is the symmetric case, ob-tained by a homogeneous distribution of frustration pa-rameters, i.e, α i = α h ∀ i . For our purposes, we consider α h >
0. In this situation, the trivial solution of the frus-tration parameters, α i = α h , is another one of the valuesout of the continuous spectrum. That is, we can recoverthe landscape given by the symmetric configuration inmany different ways. We are however, interested in com-puting the analytical expression of the cost function inorder to select the one corresponding to the minimumcost. A. Optimal cost tuning when α C = 0 Let us firstly consider the case where α C = 0 and ho-mogeneous natural frequencies ω i = ω h ∀ i . In the par-ticular case of the symmetric configuration, that is, thephase configuration given by α i = α h ∀ i the solution ofthe frustration parameters is given by: (cid:126) ˜ κ = (cid:16) ˜ M D s (cid:17) − (cid:18) K (cid:126) ˜∆ ω − ˜ L(cid:126) ˜ φ ∗ (cid:19) = − (cid:16) ˜ M D s (cid:17) − ˜ L(cid:126) ˜ φ ∗ (21)But (cid:126) ˜ φ ∗ corresponds to the symmetric case. Hence (seeSection II), (cid:126) ˜ φ ∗ = α h ˜ L − (cid:126) ˜∆ s (22)where (cid:126) ˜∆ s i ≡ (cid:104) s (cid:105) − s i and the tilde touches on kth rowremoval.Plugging Eq.(22) into Eq.(21): (cid:126) ˜ κ = − α h (cid:16) ˜ M D s (cid:17) − ˜ L ˜ L − (cid:126) ˜∆ s = − α h (cid:16) ˜ M D s (cid:17) − (cid:126) ˜∆ s But (cid:126) ˜∆ s can be written as: (cid:126) ˜∆ s = − ˜ M(cid:126)s (23)Putting it all together: (cid:126) ˜ κ = − α h (cid:16) ˜ M D s (cid:17) − ˜ L ˜ L − (cid:126) ˜∆ s = α h (cid:16) ˜ M D s (cid:17) − ˜ M(cid:126)s (cid:126) ˜ κ (24) which in vector form is written as: (cid:126) ˜ κ = α h (cid:16) ˜ M D s (cid:17) − ˜ M(cid:126)s (cid:126) ˜ κ = α h − s C s − s C s · · · − s C s N − (25)And considering the relation between α and κ , in Eq.(12): (cid:126)α = α h − s C s − s C s · · · C node) · · · − s C s N − (26)Equation (25) gives us the tuned values of the frustrationparameters as a function of the chosen control node, C ,when α C = 0. Notice that the result depends nonlinearlyonly on the ratio between the degree of each node and thecontrol node. This informs us that nodes with the samedegree would be tuned to the same value or, in otherwords, the tuning depends only on the degree sequenceof the network.Once we have computed the analytical solution of thefrustration parameters, we derive the expression of therequired cost to achieve such state with the particularchoice of C . Using the definition in Eq.(19): e T ( C ) = | α h | N − (cid:88) i (cid:12)(cid:12)(cid:12)(cid:12) − s C s i (cid:12)(cid:12)(cid:12)(cid:12) = | α h | N − (cid:88) i (cid:12)(cid:12)(cid:12) s i − s C s i (cid:12)(cid:12)(cid:12) (27)Before we provide the mathematical solution to theminimization problem defined in Eq.(20) for this partic-ular case, let us gain an intuitive understanding of it.Looking at Eq.(27) we see that the contribution of the ith node to the cost increment depends on | s C − s i | and,hence, if the chosen control node, C , has an extremevalue, i.e, s C (cid:28) s i or s C (cid:29) s i , the contribution willbe larger. On the contrary, if the degree of the controlnode is similar to that of the remaining nodes, then theincrease in cost will be smaller.For example, the network in Figure 3(a), with (cid:126)s =(1 , , , , , , ,
2) has the set of unique degrees (cid:126)s unique =(1 , ,
6) and hence three possible values of the cost, sharedby some nodes. If C = { , } , s C = 1: e T ( C ) = | α h | (cid:18) | − | + 5 | − | + | − | (cid:19) == | α h | (cid:18)
52 + 56 (cid:19) = 103 | α h | If C = { , , , , } , s C = 2: e T ( C ) = | α h | (cid:18) | − | + 4 | − | + | − | (cid:19) == | α h | (cid:18) (cid:19) = 83 | α h | And, finally, if C = 1, s C = 6: e T ( C ) = | α h | (cid:18) | − | + 5 | − | (cid:19) == | α h | (10 + 10) = 20 | α h | The minimum value of the energy is | α h | , correspondingto the choice C ∈ { , , , , } with s C = 2.Notice that the optimal choice of the control node (ornodes) does not depend on the value of α h in the sym-metric configuration, but only on the degree sequenceof the network. Moreover, this example illustrates thatthe degree of the control node corresponds to an inter-mediate value within the degree sequence of the networkand not an extreme value. A more detailed inspectionof Eq.(27) discloses that the proper choice of the controlnode (or control nodes) corresponds to the minimizationof the relative error of degrees. In order to find the par-ticular value of the degree that the control node musthave we should solve the minimization problem definedin Eq.(20): min C,x e T ( C, x ) = min
C,x N (cid:88) i α i ( C, x )which, when considering the symmetric configurationcase, turns tomin s C | α h | N − (cid:88) i (cid:12)(cid:12)(cid:12) − s C s i (cid:12)(cid:12)(cid:12) = | α h | min s C N (cid:88) i (cid:12)(cid:12)(cid:12) s C − s i s i (cid:12)(cid:12)(cid:12) (28)Equation (28) is equivalent to the minimization of theabsolute value of the relative error of the degree: | α h | min s C N (cid:88) i |E i | (29)where E i = (cid:12)(cid:12)(cid:12) s C − s i s i (cid:12)(cid:12)(cid:12) .The most general minimization problem of the relativeerror of a variable [23] can be written asmin d N (cid:88) i =1 w i | x i − d | ; d > d is the variable one is interested in and w i is theweight corresponding to each x i variable. The solutionof Eq.(30) is given by d = x m , where m ≡ min { i | i (cid:88) k =1 w k ≥ n (cid:88) k = i w k } i ∈ { , ..., n } (31)In other words, the value of d that minimizes Eq.(30)corresponds to the weighted median of the variable x or,equivalently, the 50% weighted percentile. The weightedmedian of a set n distinct ordered elements x , x , ..., x n with positive weights w , w , ..., w n , is the element x k satisfying min { i | (cid:80) ik =1 w k ≥ (cid:80) nk = i w k } . In other words,the solution is given by x k , the value such that the sumof the weights at each side of the pivot, k , are as even aspossible.The particular case defined in Eq.(28) can be mappedto the most general problem defined in Eq.(30), choosing w i = 1 /s i , x i = s i and d = s C . Accordingly, the solutionof s C corresponds to the weighted median of the set { s i } ,with weights given by the inverse of the node degree.Following the example of the network in Figure 3(a),with degree sequence (cid:126)s = (1 , , , , , , , s C by using Eq.(31).sorted( (cid:126)s ) = (1 , , , , , , , (cid:126)w = (cid:18) , , , , , , , (cid:19) To find the weighted median, we have to find the mini-mum value such that the sum of the weights at each sideof the pivot are as even as possible.1 + 1 + 12 + 12 = 3 ≥ .
17 = 12 + 12 + 12 + 12 + 16which corresponds to s C = 2, in agreement with the lo-cation of the minimum for α C = 0 in Figure 3(b) corre-sponding to the network in Figure 3(a). FIG. 3. (Color online) Implied cost to achieve the symmetricconfiguration as a function of the degree corresponding todifferent choices of the control node, s C , for a network of 8nodes (upper panels) and 9 nodes (lower panels). The distinctcolors and markers correspond to different values of α C . Thesymmetric configuration is generated by a value of α h = 0 . B. Optimal cost tuning when α C (cid:54) = 0 We next ask which is the optimal choice of the controlnode in the case we let α C (cid:54) = 0 and ω i = ω h ∀ i . Inthis case, we should look at Eq.(13) and set (cid:126) ˜∆ ω = 0.Making use of the analytical solution of the symmetricconfiguration in Eq.(22):˜ M D s (cid:126) ˜ κ = − α h ˜ L ˜ L − (cid:126) ˜∆ s − α C · → (cid:88) j [ M D s ] ij Using the properties ˜ L ˜ L − = I and (cid:126) ˜∆ s = ˜ M(cid:126)s , (cid:126) ˜ κ = α h (cid:16) ˜ M D s (cid:17) − ˜ M(cid:126)s − α C (cid:16) ˜ M D s (cid:17) − → (cid:88) j [ M D s ] ij Finally, in vector form, (cid:126) ˜ κ = α h − s C s − s C s · · · − s C s N − − α C − s C s − s C s · · · − s C s N − (32)Using Eq.(12), (cid:126)α = ( α h − α C ) − s C s − s C s · · · C node) · · · − s C s N − + α C (33)where we have used the result in Eq.(25) and the relation (cid:16) ˜ M D s (cid:17) − → (cid:88) j [ M D s ] ij = (cid:16) ˜ M D s (cid:17) − ˜ M(cid:126)s == − α C α − α C α · · · − α C α N − In the particular case that α C = α h we recover the trivialinitial configuration α i = α h ∀ i , as expected from themodel.Once we have computed the analytical solution of thefrustration parameters, we derive the expression of theimplied cost to achieve such state with the particularchoice of C . Using the definition in Eq.(19): e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) ( α h − α C ) (cid:18) − s C s i (cid:19) + α C (cid:12)(cid:12)(cid:12) (34)We next derive the analytical solution of the optimalchoice of the control node and finally proof that the globalminimum corresponds to a value of α C = 0. Equation(34) can be rearranged as e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) s i − (cid:16) − α C α h (cid:17) s C s i /α h (cid:12)(cid:12)(cid:12) (35) and thereby can be easily mapped to the solution of theminimization problem defined and solved in Eq.(30) andEq.(31), respectively. Looking at Eq.(34), we shouldchoose x i = s i , w i = α h /s i and d = (1 − α C /α h ) s C .With this choice, the value of d that minimizes Eq.(34)corresponds to the weighted median of the set { s i } withweights α h /s i . Therefore, the value of d is the same asthe solution of the case α C = 0, but d (cid:54) = s C and thus wemust apply a transformation in order to obtain the op-timal choice of s C . We have to distinguish several cases,considering α h > • α C >
0. In this case we inspect Eq.(35) and distin-guish two more cases:o α C > α h : in this case, the prefactor of s C isnegative, and we can write: e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) s i + (cid:12)(cid:12)(cid:12) − α C α h (cid:12)(cid:12)(cid:12) s C s i /α h (cid:12)(cid:12)(cid:12) == N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) α h + M α h s C s i (cid:12)(cid:12)(cid:12) where M ≡ (cid:12)(cid:12)(cid:12) − α C α h (cid:12)(cid:12)(cid:12) > s C , the minimum is achieved when s C = min( s i ) (See Figure 3 at α C = 0 . α C < α h : in this case, the prefactor of s C ispositive, and we can write: e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) s i − (cid:12)(cid:12)(cid:12) − α C α h (cid:12)(cid:12)(cid:12) s C s i /α h (cid:12)(cid:12)(cid:12) taking into account that d = (cid:12)(cid:12)(cid:12) − α C α h (cid:12)(cid:12)(cid:12) s C andconsidering that, in this case, 0 < α C < α h and hence 0 ≤ (cid:12)(cid:12)(cid:12) − α C α h (cid:12)(cid:12)(cid:12) ≤ s i ) ≤ d ≤ max( s i ),the optimal value of s C falls in the range d ≤ s C ≤ max( s i ). Hence, the optimal value of s C is always larger than the weighted median, d (see Figure 3 at α C = 0 . • α C <
0. In this case we can rewrite Eq.(35) as e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) s i − (cid:16) | α C | α h (cid:17) s C s i /α h (cid:12)(cid:12)(cid:12) and distinguish two more cases:o | α C | > | α h | : In this case, the prefactor of s C ispositive and bounded by 2 ≤ (cid:16) | α C | α h (cid:17) ≤ ∞ .In this case, d = (cid:16) | α C | α h (cid:17) s C and hence0 ≤ s C ≤ d/
2. Hence, the optimal valueof s C is always smaller than half the valueof the weighted median, d (see Figure 3 at α C = − . | α C | < | α h | :n this case, the prefactor of s C is positive and bounded by 1 ≤ (cid:16) | α C | α h (cid:17) ≤
2. In this case, d = (cid:16) | α C | α h (cid:17) s C and hence d/ ≤ s C ≤ d . Hence, the optimal value of s C is always smaller than the weighted median, d (see Figure 3 at α C = − . • α C = 0: This case is explored in Section V A. Equa-tion (35) turns to e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) s i − s C s i /α h (cid:12)(cid:12)(cid:12) . The optimal value of s C is the same as theweighted median, d , without any further transfor-mation (see Figure 3 at α C = 0 . • α C = α h : This case is discussed in the introductionof the present section. Eq.(35) turns to e T ( C ) = N − (cid:88) i =0 α h = N α h and hence the value of the cost is the same constantvalue for all nodes (see Figure 3 at α C = 0 . α C ,the global minimum cost is given by α C = 0, as shownin Figure 3. This result can be proved by considering asimplified version of Eq.(35), defined as f ( x ) = (cid:12)(cid:12)(cid:12) a − (1 − x/b ) ca/b (cid:12)(cid:12)(cid:12) (36)The minimum value of Eq.(36) is achieved when x = 0,as long as a > b > c >
0. This conditions areequivalent to s i > α h > s C >
0, and are truefor all the summation terms in Eq.(35). Therefore, theminimum value is given by setting α C = 0.Summing up, in order to obtain the optimal { α i } pa-rameters’ set in order to achieve the symmetric phaseconfiguration with the minimum implied cost in theKuramoto-Sakaguchi model, we should set α C = 0, inde-pendently of the value of α h . The remaining parametershave to be tuned using Eq.(33). Moreover, the optimalchoice of the control node (or nodes) corresponds to thatwith s C located at the weighted median of { s i } (withweight equal to s − i ).Notice also that nodes are grouped by degree regardingthe tuned values of its frustration parameters. In otherwords, there may be different potential control nodes, aslong as they share the same degree. VI. FULLY SYNCHRONIZED PHASECONFIGURATION
Another particular phase configuration is given by thephase synchronization of nodes, that is, (cid:126)φ ∗ = (cid:126)
0. If we set,as in Section V, ω i = ω h ∀ i , we end up with the trivialsolution α i = 0 ∀ i . In the case of full synchronizationwe want to recover the completely in-phase state froma phase dispersion produced by a distribution of naturalfrequencies, which we consider to be positive. Hence,applying Eq.(13) to this case: (cid:126) ˜ κ = (cid:16) ˜ M D s (cid:17) − (cid:18) K ˜ M (cid:126)ω − α C · ˜ M(cid:126)s (cid:19) (37)and in vector form, (cid:126) ˜ κ = α C ( s C − s ) − ( ω C − ω ) /Ks α C ( s C − s ) − ( ω C − ω ) /Ks · · · α C ( s C − s N − ) − ( ω C − ω N − ) /Ks N − (38)where we have used: (cid:126) ˜∆ ω = ˜ M (cid:126)ω .Finally, from the (cid:126)κ in Eq.(38) we can obtain (cid:126)α : (cid:126)α = α C s C − ( ω C − ω ) /Ks α C s C − ( ω C − ω ) /Ks · · · α C · · · α C s C − ( ω C − ω N − ) /Ks N − (39)Similarly as the result of the symmetric configuration,given in Eq.(33), the solution of the fully synchronizedconfiguration concerning (cid:126)α is a continuous spectrum ofvalues, depending on the choice of the control node, C ,the value of its frustration parameter α C , which is a freeparameter, and the natural frequencies of the oscillators.In Sections VI A and VI B we will make a in-depth anal-ysis of the problem, as well as comment on the nonlinearexpansion of the Kuramoto-Sakaguchi model and the va-lidity of our approach in this case (Section VI C). A. Optimal cost tuning when α C = 0 Using the definition of cost in Eq.(19) and the generalsolution of the frustration parameters in Eq.(39) we get: e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) α C s C − ( ω C − ω i ) /Ks i (cid:12)(cid:12)(cid:12) (40)In the particular choice α C = 0: e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) ω C − ω i Ks i (cid:12)(cid:12)(cid:12) (41)0Equation (41) shows that the relevant piece of informa-tion regarding the control node is given by its natural fre-quency, ω C . Similarly to the minimization problem posedin Section V, and in order to find the optimal choice ofthe control node we need to solve Eq.(20) considering thesolution of Eq.(41):min ω C (cid:12)(cid:12)(cid:12) ω C − ω i Ks i (cid:12)(cid:12)(cid:12) = 1 K min ω C (cid:12)(cid:12)(cid:12) ω C − ω i s i (cid:12)(cid:12)(cid:12) (42)The optimization problem is equivalent to the most gen-eral problem, described in Eq.(30), with solution givenby Eq.(31). In this case, d = ω C , x i = ω i and the weight w i = s − i . Accordingly, and in a similar way as in Sec-tion V, the solution of ω C corresponds to the weightedmedian of the set { ω i } , with weights given by the inverseof the node degree. Notice that the optimal choice ofthe control node is in general different to that given inSection V A). This is due to the fact that the weightsof the weighted median have to be sorted according todescending order of natural frequencies instead of nodedegree.Following with the example provided in Section V A,for the network in Figure 3(a), with degree sequence (cid:126)s = (1 , , , , , , , ω C by using Eq.(31). Consider the following naturalfrequencies (cid:126)ω = (0 . , . , . , . , . , . , . , .
15) (43)which lead tosorted( (cid:126)ω ) = (0 . , . , . , . , . , . , . , .
45) (44)and the corresponding weights (cid:126)w = (cid:18) , , , , , , , (cid:19) (45)To find the weighted median, we have to find the mini-mum value such that the sum of the weights at each sideof the pivot are as even as possible.12 + 1 + 12 + 16 + 12 = 2 . ≥ . C = 6 [see α C = 0 line in Figure4(a)], with ω C = 0 .
25 [see α C = 0 line in Figure 4(b)]and a degree of s C = 2. B. Optimal cost tuning when α C (cid:54) = 0 The cost corresponding to the fully synchronized con-figuration case is given by Eq.(40). In the general casewhere α C (cid:54) = 0, we can minimize the cost with respect to ω C or to s C . If we minimize with respect to ω i , we first FIG. 4. (Color online) Implied cost to achieve the fully syn-chronized configuration as a function of the chosen controlnode, C (upper panel) and natural frequencies of nodes (bot-tom panel) for the network in Figure 3(a). Five differentvalues of α C are considered (marked colored lines). Naturalfrequencies are set as the example in Eq.(43). have to rewrite Eq.(40) as e T ( C ) = N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) α C s C − ( ω C − ω i ) /Ks i (cid:12)(cid:12)(cid:12) K N − (cid:88) i =0 (cid:12)(cid:12)(cid:12) ω i − ( ω C − α C s C K ) s i (cid:12)(cid:12)(cid:12) (46)Again, the problem and the solution of Eq.(46) can betaken from Eq.(30) and Eq.(31), choosing d ≡ ω C − α C s C K , w i ≡ /Ks i and x i ≡ ω i .Hence, the value d that minimizes the cost is theweighed median considering the same weight as in Sec-tion VI A, w i = 1 /s i (notice, however, that the orderingis determined by natural frequencies and not degrees).Let us analyze the different possibilities regarding thevalues of α C , maintaining ω C and s C constant:1 • ω C > α C s C K or α C < ω C Ks C : We can write N (cid:88) i (cid:12)(cid:12)(cid:12) ω i − | ω C − α C s C K | Ks i (cid:12)(cid:12)(cid:12) The value which minimizes cost is given by d = ω k ,corresponding to the weighted median. Howeverthis is not directly the value of ω C , as d = | ω C − α C s C K | in this case. The real values of the pair { ω C , s C } are given by min C ( ω k − ( ω C − α C s C K )).Following with the example in Section VI A, thevalue of the weighted median is d = 0 .
25. In thecase we are considering, however, this is not theoptimal choice of the parameters for the controlnode. We must shift the values considering therelation between d and the other parameters. Ifwe choose α C = 0 .
1, for instance, we find that, | ω C − . s c | = 0 .
25. In Figure 4 we see that theoptimal choice is given by ω C = 0 .
4, which corre-sponds to C = 5 and s C = 2. • ω C < α C s C K or α C > ω C Ks C : We can write N (cid:88) i (cid:12)(cid:12)(cid:12) ω i + | ω C − α C s C K | Ks i (cid:12)(cid:12)(cid:12) Hence, as the function increases with increas-ing ( ω i − α C s i K ), the minimum is achieved bymin C ( ω C − α C s C K ). C. Non-linear expansion of theKuramoto-Sakaguchi model
The results obtained in Section VI are based on a linearapproximation of the Kuramoto-Sakaguchi model. Wehave derived the results based on the phase synchroniza-tion requirement, and assuming that frequency synchro-nization is already achieved in the steady state. Nev-ertheless, when measuring the order parameter with alarge dispersion of natural frequencies or low couplingconstant, we do not expect such steady state. However,we ask to which extend the proposed values of the ob-tained frustration parameters are also able to enhancefrequency synchronization considering the original non-linear Kuramoto-Sakaguchi model:˙ θ i = ω i + K (cid:88) j W ij sin( θ j − θ i − α i ) (47)We compare the results from Ref. [21] considering its Type II frustration parameters tuning for both the lin-ear and the nonlinear Kuramoto model and we find that,despite our approach does not consider the enhancementof frequency synchronization on the nonlinear regime, itis able to improve the value of the order parameter, in asimilar fashion as in Ref. [21]. This work considers the nonlinear Kuramoto-Sakaguchi model and seeks to im-prove the number of nodes that fall into the recruitmentcondition so as to achieve the same common oscillatoryfrequency. The considered network class is the same asthe mentioned paper, as well as the statistics study.We make use of the expression in Eq.(39) to tune theset of (cid:126)α for a given configuration of random (cid:126)ω and studythe effect on the synchronization of the system for differ-ent values of the coupling strength.We consider two cases: the linear and the nonlinearmodel with natural frequencies obtained from a uniformdistribution ω i ∈ [ − , VII. CONCLUSIONS
The Kuramoto-Sakaguchi model adds to the originalKuramoto model a homogeneous phase lag, α , betweennodes which promotes a phase shift between oscillators.We consider a more general framework, in which thephase lag or the frustration parameter, α i , is an intrin-sic property of each node. A very relevant question inoscillatory models is finding the conditions of networksynchronization. In the present work, we bring forwarda methodology not only to obtain the desired synchro-nized state, but any convenient phase configuration inthe steady state, by means of a fine tuning of the phaselag or frustration parameters, { α i } . We feature the ana-lytical solution of frustration parameters so as to achieveany phase configuration, by linearizing the most generalmodel. The three intrinsic parameters of the nodes in themodel, natural frequencies, { ω i } , frustration parameters2 FIG. 5. (Color online) Average order parameter, (cid:104) r (cid:105) , as afunction of the coupling strength, K , for the linear [in panel(a)] and the nonlinear [in panel (b)] Kuramoto-Sakaguchi(KS) model on regular random graphs with homogeneousnode degree s i = 4 and N = 100 nodes. Natural frequenciesare obtained from a uniform random distribution in the range ω i ∈ [ − , { α i } : the original Ku-ramoto dynamics or α i = 0 ∀ i (spotted continuous red line);type II [21] KS dynamics with the frustration parameters setto sin( α i ) = − ω i / ( Ks i ) if | ω i | < Ks i and sin( α i ) = ± { α i } , and phases in the steady state φ ∗ i , are coupled byan equation that allows to tune them for a desired con-figuration. While the set φ ∗ i is uniquely determined, theset α i has a continuous spectrum of solutions.A main result is that a given phase configuration can beaccess via a continuous spectrum of frustration param-eters, i.e, one phase and one frustration parameter areleft as free parameters. The nodes we choose their valuesconcerning phase and frustration parameter, are namedreference and control nodes, respectively. Once the frus-tration parameters are tuned so as to obtain the desiredstate, we define a cost function to assess the overheadthat the system requires to achieve such parameters’ con-figuration. Among all possible tuning solutions of { α i } ,we request those which minimize the cost to obtain them.We develop the analytical solution of the cost function forthe cases of symmetric configuration and fully synchro-nized state and discuss them.A key result is the solution to the minimization costproblem: For the case of symmetric configuration, thenodes which are to be set as control nodes are thosewhose degree is the weighed median of the sample, witha weight equal to the inverse of its degree. On the otherhand, for the case of fully synchronized state, controlnodes are those whose natural frequency is the weightedmedian of the sample, with a weight equal to the inverseof its degree. An extensive analysis of several cases isdone in the text and a detailed example of a toy networkis provided. We highlight the connection made with thenonlinear Kuramoto-Sakaguchi model. Despite our anal-ysis being based on the linear version of the model, weshow that the proposed parameters’ tuning is also able toenhance frequency synchronization, as done in Ref. [21].We stress the fact that the question ‘among all the pos-sible solutions, which is the one that makes the systemachieve a particular phase configuration with the mini-mum required cost?’ is of particular relevance when weconsider the plausible real nature of the system. If a realsystem needs to access a particular phase configuration,which may be associated with a singular function, then itwill tend to minimize the effort or cost to do so. Furtherwork can be done within this framework by doing realexperiments on measuring the energy needed to access aparticular configuration. Moreover, other nonlinear os-cillatory models can be analyzed and compared with theKuramoto-Sakaguchi model.Other questions regarding the model are left open. Wehave considered the coupled trio of natural frequencies-frustration parameters-steady state phases. A naturalextension to this would be to inspect the possibility toalso tune the weights of the network edges in order toaccess a particular configuration. The higher dimensionof the latter with respect to the vectors of parameterswould require further assumptions about the model orthe network structure, such as positive weights or partic-ular distributions or topologies. Another research venuewould be to consider the effect of removing a node of thenetwork and the { α i } set needed to minimize the effect3on the removal on the whole network.Despite we provide the analytical solution to the opti-mal choice of parameters in order to minimize the costof achieving both the symmetric and the fully synchro-nized configurations, the access to all nodes’ parametersrequirement may not be feasible in real-world networks.Our methodology is quite general and the optimizationprocedure refers to a set of parameters to be tuned. Inparticular, a finite subset of nodes with accessible phase-lag parameter could be chosen (the choice could be re-stricted to any subset of nodes), holding all other nodesunaltered. This would provide a nonoptimal global con-dition but a restricted and approximated one that coulddeal with a subset of available nodes. A meaningful anal-ysis would be to identify which subset of nodes is the onethat enables to get a closest approximate solution, andrelate those nodes with their topological properties, al-though this question is beyond the goal of this work. ACKNOWLEDGMENTS
The authors acknowledge financial support fromMINECO via Project No. PGC2018-094754-B-C22(MINECO/FEDER,UE) and Generalitat de Catalunyavia Grant No. 2017SGR341. G.R.-T. also acknowledgesMECD for Grant No. FPU15/03053
Appendix A: Step-by-step derivation of the example
Consider the network in Figure 1, with its Laplacianmatrix: L = − − − − − − − − − − − − − − − We develop the equation step by step: (cid:88) j L ij θ ∗ j = ω i K − (cid:104) ω (cid:105) K + (cid:104) αs (cid:105) − α i s i ∀ i − − − − − − − − − − − − − − − · θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ = ω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) K − − − − − − − −
17 67 − − − − − − −
17 67 − − − − − − −
17 67 − − − − − − −
17 67 − − − − − − −
17 67 − − − − − − −
17 67 ·· · α α α α α α α = ω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) Kω −(cid:104) ω (cid:105) K + −
247 27 27 27 27 27 2747 −
127 27 27 27 27 2747 27 −
127 27 27 27 2747 27 27 −
127 27 27 2747 27 27 27 −
127 27 2747 27 27 27 27 −
127 2747 27 27 27 27 27 − · α α α α α α α ω i = ω ∀ i : − − − − − − − − − − − − − − − · θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ θ ∗ == −
247 27 27 27 27 27 2747 −
127 27 27 27 27 2747 27 −
127 27 27 27 2747 27 27 −
127 27 27 2747 27 27 27 −
127 27 2747 27 27 27 27 −
127 2747 27 27 27 27 27 − · α α α α α α α Now we choose R = 0 and C = 1, i.e., all φ i = θ i − θ and κ i = α i − α :Let us write explicitly the change of variables (the redand the blue columns are ones we can remove due to thechange of variables, as they do not affect to the systemof equations anymore): − − − − − − − − − − − − − − − · φ ∗ = 0 φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ == −
247 27 27 27 27 27 2747 −
127 27 27 27 27 2747 27 −
127 27 27 27 2747 27 27 −
127 27 27 2747 27 27 27 −
127 27 2747 27 27 27 27 −
127 2747 27 27 27 27 27 − ·· κ κ = 0 κ κ κ κ κ + − α α α α α α α If we look carefully at Eq.(A1), we see that although theleft-hand side and the right hand-side matrices are bothsingular, the first one has both column and row sumsequal to zero, while the second one has only column sumequal to zero. This is reflected in the additional constantterm that appears when doing the change of variablesregarding α i , which can be written as: b i = (cid:88) j [ M · D s ] ij (cid:54) = 0 in general (A1) We can choose whatever row to remove from either sides.We choose row 0: − − − − − − − · φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ ==
47 27 27 27 27 2747 −
127 27 27 27 2747 27 −
127 27 27 2747 27 27 −
127 27 2747 27 27 27 −
127 2747 27 27 27 27 − · κ κ κ κ κ κ + α α α α α α In this situation, we can solve for the set (cid:126) ˜ κ : κ κ κ κ κ κ =
47 27 27 27 27 2747 −
127 27 27 27 2747 27 −
127 27 27 2747 27 27 −
127 27 2747 27 27 27 −
127 2747 27 27 27 27 − − ·· − − − − − − − · φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ φ ∗ − α α α α α α Which leads to the result: κ κ κ κ κ κ = − α +3 φ + φ +3 φ φ − φ )22 φ − φ − φ + φ φ − φ + φ − φ + φ φ − φ + φ − φ − φ φ φ + φ − φ (A2)Note that κ = 0 and φ = 0.Consider the following configuration: (cid:126) ˜ φ ( R =0) = (0 . , . , . , − . , − . , . (cid:126) ˜ κ ( C =1) = (0 . − α / , − . , − . , . , . , − . κ i ≡ α i − α C ⇒ α i = κ i + α C .If we choose α C = 0 ⇒ α i = κ i , then we can includethe value of the control node C = 1: (cid:126) ˜ α = (0 . , . , − . , − . , . , . , − . α C = α =0 . (cid:126) ˜ α = (0 . , . , − . , − . , . , . , .
045 36 1 2Node 0Node 1Node 2Node 3Node 4Node 5Node 6
FIG. 6. For the network in Fig 1, phases obtained after tuningthe set of frustration parameters to the symmetric configura-tion (four distinct symmetries).
Appendix B: Mathematical solution of the costoptimization problem and intuitive insights
In order to gain a more intuitive understating of theanalytical expression and solution of the considered costfunction, we consider the analysis of the continuous case.
1. Symmetric configuration case
Considering the symmetric configuration case andchoosing α C = 0, the continuous optimization problemcan be written as ∂e T ( C ) ∂s C = ∂∂s C | α h | N − (cid:88) i (cid:12)(cid:12)(cid:12) − s C s i (cid:12)(cid:12)(cid:12) = | α h | N − (cid:88) i sgn( s C − s i ) s i (B1)Equation (B1) is based on the function f ( x ) = (cid:12)(cid:12)(cid:12) a − xa (cid:12)(cid:12)(cid:12) a, x > a and the sum of all of them. Regardless of the set of a i values, the sum function (cid:80) i f ( x, a i ) (see the example inthe black line in Figure 7), is a concave function and hasa unique minimum, which corresponds to one of the a i values.In order to assess the value of a i where the minimumis located, we compute the derivative of Eq.(B2): df ( x ) dx = sgn( x − a ) a (B3) x f ( x ) = | x | g ( x ) = | x | h ( x ) = | x | f ( x ) + g ( x ) + h ( x ) FIG. 7. Three examples of the general function f ( x ) = | a − xa | ,with a = 2, a = 3 and a = 4, and the resulting sum of them. x f ( x ) + g ( x ) + h ( x ) f ( x ) + g ( x ) + h ( x ) FIG. 8. Derivative of the function f ( x ) = | − x | + | − x | + | − x | defined in Figure 7. Red dashed line at y = 0. and hence, d (cid:80) i f ( x, a i ) /dx = (cid:80) i sgn( x − a i ) /a i , which isdepicted in Figure 8. Notice that, despite the derivativeof the function is not defined at the values x = a i , thederivative changes its sign when moving from x < x > a i .To conclude, Eq(B1) behaves equivalently as the func-tion defined in Eq.(B2) and hence, displays only one min-imum, which is achieved at the s i where there is a change6of sign in the derivative.Alternatively and as explained in the main text, we canunderstand the minimization problem as part of a generalframework. The minimization of Eq.(B1) is equivalentto the minimization of the absolute value of the relativeerror: N − (cid:88) i (cid:12)(cid:12)(cid:12) − s C s i (cid:12)(cid:12)(cid:12) = N (cid:88) i (cid:12)(cid:12)(cid:12) s C − s i s i (cid:12)(cid:12)(cid:12) = N (cid:88) i |E i | (B4)The general problem can be written as [23]:min d N (cid:88) i =1 w i | x i − d | ; d > d = x m where m ≡ min { i | i (cid:88) k =1 w k ≥ n (cid:88) k = i w k } i ∈ { , ..., n } (B5)In other words, the d value that minimizes Eq.(B5) cor-responds to the weighted median of the variable x , or the50% weighted percentile. Weighted median:
For n distinct ordered ele-ments x , x , ..., x n with positive weights w , w , ..., w n ,the weighted median is the element x k satisfyingmin { i | (cid:80) ik =1 w k ≥ (cid:80) nk = i w k } Therefore, the solution is given by x k , the value suchthat the sum of the weights at each side of the pivot, k ,are as even as possible.Our problem is a special case of the discrete weightedmedians with weights 1 /s i , which are a special case ofthe medians of a measure.Following the example provided in Figure 7, { x } = { , , } and { w } = { / , / , / } .The weighted median is achieved for k = 2, corre-sponding to x = 3 and weight w = 1 / / / / > / / /
12. Conversely, if we let k = 1, andhence x = 2 and w = 1 /
2, the condition on the weightswill not be true: 1 / ≯ / / / [1] G. Nicolis and C. Nicolis, Foundations of Complex Sys-tems: Nonlinear Dynamics, Statistical Physics, Informa-tion and Prediction (World Scientific, 2007).[2] A. Barrat, M. Barthelemy, and A. Vespignani,
DynamicalProcesses on Complex Networks (Cambridge UniversityPress, 2008).[3] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchroniza-tion: A universal concept in nonlinear sciences (Cam-bridge University Press, Cambridge, UK, 2001).[4] G. Osipov, J. Kurths, and C. Zhou,
Synchronizationin Oscillatory Networks , Springer Series in Synergetics(Springer Berlin Heidelberg, 2007).[5] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, andC. Zhou, Synchronization in complex networks, PhysicsReports , 93 (2008).[6] F. D¨orfler and F. Bullo, Synchronization in complex net-works of phase oscillators: A survey, Automatica ,1539 (2014).[7] S. Boccaletti, A. N. Pisarchik, C. I. del Genio,and A. Amann, Synchronization (Cambridge UniversityPress, 2018).[8] Y. Kuramoto, Self-entrainment of a population of cou-pled non-linear oscillators, in
International Sympo-sium on Mathematical Problems in Theoretical Physics (Springer-Verlag) pp. 420–422.[9] H. Sakaguchi and Y. Kuramoto, A Soluble Active Ro-tator Model Showing Phase Transitions via MutualEntrainment, Progress of Theoretical Physics , 576(1986).[10] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ri-tort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Reviews ofModern Physics (2005).[11] Y.-Y. Liu and A.-L. Barab´asi, Control principles of com-plex systems, Rev. Mod. Phys. , 035006 (2016).[12] N. Frolov, V. Maksimenko, S. Majhi, S. Rakshit,D. Ghosh, and A. Hramov, Chimera-like behavior in aheterogeneous kuramoto model: The interplay betweenattractive and repulsive coupling, Chaos: An Interdisci-plinary Journal of Nonlinear Science , 081102 (2020).[13] Y. Kuramoto, Chemical oscillations, waves, and turbu-lence (Dover Publications, 2003) p. 156.[14] V. Nicosia, M. Valencia, M. Chavez, A. D´ıaz-Guilera,and V. Latora, Remote Synchronization Reveals NetworkSymmetries and Functional Modules, Physical ReviewLetters , 174102 (2013).[15] T. Nishikawa and A. E. Motter, Symmetric states re-quiring system asymmetry, Physical Review Letters ,114101 (2016).[16] F. Molnar, T. Nishikawa, and A. E. Motter, Networkexperiment demonstrates converse symmetry breaking,Nature Physics , 351 (2020).[17] Y. Zhang and A. E. Motter, Symmetry-independent sta-bility analysis of synchronization patterns, SIAM Rev. , 817 (2020).[18] G. Rosell-Tarrag´o and A. D´ıaz-Guilera, Functionabilityin complex networks: Leading nodes for the transitionfrom structural to functional networks through remoteasynchronization, Chaos , 062315 (2020).[19] X. Lei, X. Li, and D. Povh, A nonlinear control for coor-dinating TCSC and generator excitation to enhance thetransient stability of long transmission systems, Electric Power Systems Research , 103 (2001).[20] U. R. Acharya, K. P. Joseph, N. Kannathal, L. C. Min,and J. S. Suri, Heart Rate Variability, in Advances inCardiac Signal Processing (Springer Berlin Heidelberg,Berlin, Heidelberg, 2007) pp. 121–165.[21] M. Brede and A. C. Kalloniatis, Frustration tuningand perfect phase synchronization in the Kuramoto-Sakaguchi model, Physical Review E , 062315 (2016).[22] P. S. Skardal, D. Taylor, and J. Sun, Optimal synchro-nization of complex networks, Physical Review Letters , 144101 (2014).[23] E. Semsar-Kazerooni and K. Khorasani, Multi-agentteam cooperation: A game theory approach, Automat-ica , 2205 (2009).[24] R. M. D’Souza, J. G´omez-Garde˜nes, J. Nagler, andA. Arenas, Explosive phenomena in complex networks,Advances in Physics , 123 (2019).[25] J. G´omez-Garde˜nes, S. G´omez, A. Arenas, andY. Moreno, Explosive synchronization transitions inscale-free networks, Phys. Rev. Lett.106