Optimal cost tuning of frustration: Achieving desired states in the Kuramoto-Sakaguchi model

There are numerous examples of studied real-world systems that can be described as dynamical systems characterized by individual phases and coupled in a network like structure. Within the framework of oscillatory models, much attention has been devoted to the Kuramoto model, which considers 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-Sakaguchi model, which adds a phase lag parameter to each node. We construct a general formalism that allows to compute the set of lag parameters that may lead to any phase configuration within a linear approximation. In particular, we devote special attention to the cases of full synchronization and symmetric configurations. We show that the set of natural frequencies, phase lag parameters and phases at the steady state is coupled by an equation and a continuous spectra of solutions is feasible. In order to quantify the system's strain to achieve that particular configuration, we define a cost function and compute the optimal set of parameters that minimizes it. Despite considering a linear approximation of the model, we show that the obtained tuned parameters for the case of full synchronization enhance frequency synchronization in the nonlinear model as well.


I. INTRODUCTION
Emergence is one of the key concepts in the analysis of complex systems [1]. Collective properties emerge as a consequence of irregular interactions among its elemental constituents [2]. One of the most paradigmatic examples of emergence is synchronization [3,4], because the interplay between populations of oscillatory units gives rise to a variety of global states, ranging from perfect synchronization or phase locked stationary configurations to chimera states [5][6][7]. Among the different models that have been used to understand such collective behavior, a lot of effort has been devoted to the Kuramoto model (KM), in which phase oscillators interact continuously with other units through a sine function of the phase difference [8][9][10].
In the past few years there has been a growing interest in the concept of controllability, which quantifies the feasibility to achieve a desired final state of a given dynamical system [11]. As stated above, the KM can give rise to a wide variety of stationary (phase or frequency synchronized) or not stationary states, being chimeras an unexpected mixture of both types of behaviors [12]. In this context, controllability can be understood as a tuning of the internal parameters of the oscillators to reach specific phase configurations. The most simple settings stand for a collection of identical oscillators interacting through a sinus function of the phase differences. In this case it is quite intuitive to see that the final state is a perfectly synchronized one in which all oscillators have exactly the same phase and frequency (the same frequency * gemmaroselltarrago@gmail.com † albert.diaz@ub.edu than the intrinsic one). It is the existence of a distribution of frequencies that gives rise to a transition, in terms of the strength of the coupling, from an incoherent state to a coherent one [13]; such a transition is robust in the sense that the introduction of a lag term, a phase added to the argument of the sinus function, does not change the behavior, as far as it is kept below π/2 [9]. However, the introduction of this lag term for identical oscillators changes completely the structure of the, in principle, synchronized state. In Reference [14] it was shown that, for small and common values of the lag parameters, the synchronized state breaks into partially synchronized groups of oscillators, being symmetry the reason for the phase synchronization of the oscillators. When increasing this common lag parameter the system enters into a incoherent chaotic state. Actually, there has been an increasing interest in the last months on the role that symmetries plays in the synchronization of oscillatory units and how the lack of homogeneity in some of the parameters can be compensated by other choices [15][16][17].
In a previous work we introduced the concept of "functionability" as a measure of the ability of a given node to change the state of the system by just tuning one internal variable, the node lag in the argument of the sinus function of the interaction [18]. Being an intrinsic property of the node, its change produces a global change in the phases of the system of oscillators that can be measured. Functionability stands for the reaction of the whole phase distribution to a small change in a node. The analytical expression of functionability reports its quadratic dependence on the node degree and the node lag value, but also a structural term, such that the most peripheral nodes in the network have also larger contributions to functionability centrality measure. The nodes with higher functionability values may represent positive arXiv:2101.12171v1 [physics.soc-ph] 28 Jan 2021 actors for the network, because they enable more variability in the states of the system, but also potentially dangerous ones, as tiny perturbations can produce cascade like effects that completely changes the network dynamics.
As stated, the addition of a phase lag parameter enables a richer configuration state. However, it is clear that a tuning of a single parameter will not be enough to generate the wide variety of stationary states that a population of Kuramoto oscillators can achieve. Notwithstanding, the question that arises is whether a fine tuning of a set of individual parameters can make it possible. In this paper, this is our proposal, we construct a general formalism that allows, within a linear approximation, to compute the set of lag parameters that may lead to any phase 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 the configuration of phases that is produced by a set of lag parameters.
There are numerous examples of real-world systems that can be described as dynamical systems characterized by individual phases and which functioning are object of investigation. Some examples are the brain functional networks arising from temporal correlation patterns, ac power in power grids [19], heartbeats [20], multiprocessors and multicore processors, or traffic signaling. Not only the synchronization between their constituents may be intended or prevented, but also other particular configurations may be of relevant interest. For this reason, we propose a mechanism for tuning the intrinsic parameters of the system to achieve any desired phase configuration.
A previous work proposes a methodology to enhance frequency synchronization for the nonlinear Kuramoto-Sakaguchi model (extension of the Kuramoto model with a node phase lag parameter) [21]. Another work suggests that an unstable synchronized state becomes stable when, and only when, the oscillator parameters are tuned to nonidentical values [15]. We highlight the work done in Reference [22], where the particular configuration of perfect synchronization is studied and the synchrony alignment function is defined in order to minimize the order parameter of the system considering different topologies and frequency scenarios. We address the most general question, following a similar path to that pursued by them, forcing the system to achieve any particular configuration for the linear case of the Kuramoto-Sakaguchi model by means of a fine tuning of the phase lag or frustration parameter set. Despite considering a linear approximation of the model, we show that the obtained tuned parameters for the case of full synchronization enhance frequency synchronization in the nonlinear model as well.
The structure of the paper is the following. In Section II we present the Kuramoto-Sakaguchi model (a variation of the Kuramoto model) as the proper framework for our purposes. Next, in Section III, we derive the analytic expression of the fine tuning of the frustration parameters so as to achieve any phase configuration. Then, in Section IV, we define a cost function to assess the expense of achieving a particular phase configuration by its corresponding tuning and derive the analytic solution for the cases of symmetric and fully synchronized configurations, in Sections V and VI, respectively, as well as comment on the nonlinear validity of our results. We conclude in Section VII. An easy-to-follow example and further mathematical derivations can be found in the Appendix.

II. THE KURAMOTO-SAKAGUCHI MODEL
In 1975, Kuramoto suggested one of the best-known dynamical equation to model interacting oscillatory systems [8]: a set of N phase oscillators characterized by their phase, θ i , and coupled between each other by the sine of their phase differences. Each unit is influenced directly by the set of its nearest neighbours via the adjacency matrix of the network corresponding to the system, G(V, E). The coupling strength describes the intensity of such pair-wise interactions, K ij . The set of nodes of the network, V(G) consists of all the oscillators, while the set of edges, E(G), is made of the links between them. In his original work, Kuramoto assumed homogeneous interactions, i.e, K ij = K ∀(i, j) [8,10]. Taking into account the connectivity or topology of the network and the oscillatory dynamics, the dynamics of the system can be written as a system of differential equations [5]: where Γ i is the set of neighbors of node i and ω i is the natural frequency of each unit. We consider that two nodes are phase synchronized when their phases have the same value, When the phase difference has a constant value, that is, θ i (t) − θ j (t) = c ∀t > t 0 , we say there is a phase locking between nodes i and j. Similarly, we consider that two nodes are frequency synchronized when their frequencies have the same value: We say that two nodes are fully synchronized when they are phase synchronized, because this implies frequency synchronization.
As long as the distribution of natural frequencies is homogeneous, namely, all units have the same natural frequency, there is only one attractor of the dynamics: the fully synchronized state. It can be shown that, if the distribution of natural frequencies is unimodal, the system becomes frequency synchronized as long as the coupling strength is larger than a threshold value [10].
In 1986, Kuramoto together with Sakaguchi presented a similar model which incorporated a constant phase lag between oscillators [9] which can be written as follows: (2) where α is a homogeneous phase lag parameter. This model has also become well-known and several variations of the model have been studied. It has been shown that, as long as |α| < π/2, the system is not chaotic and a threshold value for the coupling strength exists above which the system becomes synchronized to a resulting frequency [9]. In the particular case that considers homogeneous natural frequencies, i.e, ω i = ω 0 ∀i, the frustration parameter, α, forces the system to break the otherwise original fully synchronized state. However, partial synchronization is conserved for symmetric nodes in the network [14,15]. As the frustration increases the asynchronous groups' phase move away from each other.
We are interested here in a more general case, where the frustration parameter is not homogeneous but an intrinsic property of each unit: (3) In this context, a recent work studies a particular effect of this frustration parameter by defining functionability, a new centrality measure of the nodes in a network, in order to address the issue of which nodes, when perturbed, move the system from a synchronized state to one that is more asynchronous in the sense that it enhances the phase differences between all pairs of oscillators [18].

III. ANALYTIC EXPRESSION OF THE FRUSTRATION PARAMETERS TUNING
We address the most general problem, which considers the same dynamics as in Eq.(3), while allowing the edges of the network to be weighted, a more realistic scenario for real-world networks.
For small values of the frustration parameters and phases close to each other, which is the case in frequency synchronization, we can linearize Eq.(3) as follows: where W ij is the value of the weight of the edge between node i and node j, s i ≡ j W ij is the weighted degree of the ith node and L is the weighted Laplacian matrix defined as L ij ≡ δ ij s i − W ij . In the stable regime, a synchronized frequency is achieved and, for all oscillatorsθ i = Ω. We can derive the value of the common frequency oscillation, Ω, summing Eq.(4) over i: Taking into account the steady stateθ i = Ω ∀i and arranging summations: and finally, where we have used the Laplacian matrix property: i L ij = 0 and defined the averages i α i s i /N = αs and i ω i /N = ω . Now we can plug expression Eq. (7) to Eq.(4) to get the stable phases of oscillators, θ * i : The solution of Eq.(8) regarding phases is undetermined due to the singular nature of the Laplacian matrix. Hence, Eq.(8) is, in general, an undetermined system of linear equations, that is, there is one free phase, which we should use as a reference value for the solution. Nonetheless, we do not work directly with the functional form of phases because they are time dependent {θ * i } = f i (t), but with the phase differences with respect to a reference node, once the stationary state is achieved, In this way, we work with time independent values. In this situation, φ R = 0, by definition, as φ R ≡ θ R −θ R = 0. On the other hand, the contribution αs − α i s i of the right-hand side of Eq.(8) can be written in matrix form as: and D s ≡ where ∆ω i ≡ ω i − ω . Finally, we obtain the set of unknowns {α i }: Equation (11), however, does not have a solution, because of the singular nature of M · D s matrix. M matrix is singular too, and hence, its inverse matrix does not ex- Similarly as we did for phases [18], we solve the singularity problem by setting a reference node, which we call control node, regarding frustration parameters, i.e., we would not obtain the value for each of the parameters, but a relation between them: where α C is the value of the control node. In this situation, κ C = 0, by definition, as To easily write the matrix expressions, we define the selection matrix J (n,m) , which is, in general, an (N −1)× (N − 1) identity matrix after the removal of the mth row and the nth column.
L θ * turns toL(k, R) φ * , where we have removed the kth row and the Rth column. The result does not depend on which row we remove, hence we can choose any k. Using the selection matrix,L(k, In an equivalent way as the definition of the reduced Laplacian:M whereM D s is M D s without the kth row and the Cth column. Similarly, κ(k) = J (,k) · κ ≡ κ and ∆ ω(k) = J (,k) · ∆ω ≡ ∆ ω, where we have removed the kth row.
Considering all the previous definitions and remarks, Eq.(10) can be rewritten as: (13) and finally, where we have used 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 to set α C = 0 if we want to avoid extra constant arrays in the final expression. In this particular case: and keep in mind that κ C = 0.
The obtained values of α depend on both the chosen control node, C, and the value we set for its frustration parameter, α C . Notice, therefore, that there is a continuous spectrum of values for the frustration parameter in order to achieve a particular phase configuration.
Moreover and more importantly, due to the non-rowsum equal to zero of M D s matrix, the differences between the obtained values are dependent of the control node choice. Mathematically, This property will lead us to the definition of a cost for the system to move to the final configuration, which will depend on both the control node and the value of its frustration parameter.
We provide an example of a toy network for the case of a homogeneous natural frequencies distribution, i.e., ω i = ω ∀i. In this case, Eq.(14) turns to: which in the case of the network depicted in Figure  leads to the solution where we have chosen κ 1 = 0 and φ 0 = 0. Hence, the results are written as a function of the value α 1 and φ i i = 0. Therefore, we can achieve any phase configuration, given by the set {φ i } by tuning the frustration parameters set {α}, where α i = κ i + α C .
To illustrate how we obtain the final values, let us consider the following phase configuration: In the general case where α C = α 1 = 0: If we choose α C = 0, then α i = κ i , we can include the value of the control node C = 1: which represents ∼ 0.3% of relative error with respect to the initial Eq. (17). See the full derivation of the analytical solution in Appendix A.

IV. OPTIMAL COST TUNING OF FRUSTRATION
As pointed out in Section III, there is a continuous spectrum of values for the choice of the frustration parameters that enables the system access a particular phase configuration. The following question arises naturally: Among all the possible solutions, which is the one that makes the system achieve a particular phase configuration with the minimum required cost?
This question is of particular relevance when we consider the plausible real nature of the system. If a real system needs to access a particular phase configuration, which may be associated with a precise function, it will tend to minimize the effort or cost to do so.
In order to quantify the required cost, we define it as follows: Henceforth, the cost associated to each node is given by the absolute value of the required frustration parameter. The absolute value operator allows for a sign-free contribution of each node, a very convenient choice in the case that the system is not beforehand specified, and a general definition is proposed instead. Furthermore, unlike other nonlinear cost functions such as the square sum of the 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 the chosen control node, C, as well as the particular choice of its frustration parameter, α C .
The optimal configuration is given by the solution of the minimization problem where the x variable is not yet specified. Depending on the problem we are interested in we would set it either to ω i , s i or any other combination of the parameters of the model. The minimal value of the cost will depend on the proper choice of the control node, C, in addition of the particular value of its frustration parameter, α C , as the free parameter left to be set. In Sections V and VI we provide a thorough analysis of it. The cost required to achieve a particular phase configuration depends on that configuration, the control node and the chosen value of α C . In Figure 2 we present an example, following with the network presented in Section III and choosing different values of α C , we compute numerically the values of the required cost using Eq. (19) to achieve the phase configuration given in Eq. (17). Notice that the global minimum depends on the control node and its frustration parameter. In Section III we have derived the general analytical solution of the frustration parameters as a function of a particular choice for the phase configuration. In this section we have defined a cost function in order to assess the optimal choice of such configuration.
Depending on the phase configuration one is interested in achieving, results will vary and the analytical expressions will have different features.
In the following sections we will focus on two particular configurations, due to its intrinsic importance, in order to obtain and discuss the analytical solution of Eq.(20): The configuration 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 the Kuramoto-Sakaguchi model is the symmetric case, obtained by a homogeneous distribution of frustration parameters, i.e, α i = α h ∀i. For our purposes, we consider α h > 0. In this situation, the trivial solution of the frustration parameters, α i = α h , is another one of the values out of the continuous spectrum. That is, we can recover the landscape given by the symmetric configuration in many different ways. We are however, interested in computing the analytical expression of the cost function in order to select the one corresponding to the minimum cost.
A. Optimal cost tuning when αC = 0 Let us firstly consider the case where α C = 0 and homogeneous natural frequencies ω i = ω h ∀i. In the particular case of the symmetric configuration, that is, the phase configuration given by α i = α h ∀i the solution of the frustration parameters is given by: But φ * corresponds to the symmetric case. Hence (see Section II), where ∆ s i ≡ s − s i and the tilde touches on kth row removal. Plugging Eq. (22) into Eq. (21): But ∆ s can be written as: Putting it all together: which in vector form is written as: And considering the relation between α and κ, in Eq. (12): Equation (25) gives us the tuned values of the frustration parameters as a function of the chosen control node, C, when α C = 0. Notice that the result depends nonlinearly only on the ratio between the degree of each node and the control node. This informs us that nodes with the same degree would be tuned to the same value or, in other words, the tuning depends only on the degree sequence of the network.
Once we have computed the analytical solution of the frustration parameters, we derive the expression of the required cost to achieve such state with the particular choice of C. Using the definition in Eq. (19): Before we provide the mathematical solution to the minimization problem defined in Eq. (20) for this particular 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 extreme value, i.e, s C s i or s C s i , the contribution will be larger. On the contrary, if the degree of the control node is similar to that of the remaining nodes, then the increase in cost will be smaller.
Notice that the optimal choice of the control node (or nodes) does not depend on the value of α h in the symmetric configuration, but only on the degree sequence of the network. Moreover, this example illustrates that the degree of the control node corresponds to an intermediate value within the degree sequence of the network and not an extreme value. A more detailed inspection of Eq.(27) discloses that the proper choice of the control node (or control nodes) corresponds to the minimization of the relative error of degrees. In order to find the particular value of the degree that the control node must have we should solve the minimization problem defined in Eq. (20): which, when considering the symmetric configuration case, turns to Equation (28) is equivalent to the minimization of the absolute value of the relative error of the degree: where The most general minimization problem of the relative error of a variable [23] can be written as where d is the variable one is interested in and w i is the weight corresponding to each x i variable. The solution of Eq.(30) is given by 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 weighted median of a set n distinct ordered elements x 1 , x 2 , ..., x n with positive weights w 1 , w 2 , ..., w n , is the element x k satisfying min{i| i k=1 w k ≥ n k=i w k }. In other words, the solution is given by x k , the value such that the sum of the weights at each side of the pivot, k, are as even as possible.
The particular case defined in Eq.(28) can be mapped to the most general problem defined in Eq.(30), choosing w i = 1/s i , x i = s i and d = s C . Accordingly, the solution of 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 s = (1, 6, 2, 1, 2, 2, 2, 2), let us compute the optimal value of s C by using Eq.(31). sorted( s) = (1, 1, 2, 2, 2, 2, 2, 6) To find the weighted median, we have to find the minimum value such that the sum of the weights at each side of the pivot are as even as possible. which corresponds to s C = 2, in agreement with the location of the minimum for α C = 0 in Figure 3(b) corresponding to the network in Figure 3(a). We next ask which is the optimal choice of the control node in the case we let α C = 0 and ω i = ω h ∀i. In this case, we should look at Eq.(13) and set ∆ ω = 0. Making use of the analytical solution of the symmetric configuration in Eq. (22): Using the propertiesLL −1 = I and ∆ s =M s, Using Eq. (12), where we have used the result in Eq.(25) and the relation In the particular case that α C = α h we recover the trivial initial configuration α i = α h ∀i, as expected from the model.
Once we have computed the analytical solution of the frustration parameters, we derive the expression of the implied cost to achieve such state with the particular choice of C. Using the definition in Eq. (19): We next derive the analytical solution of the optimal choice of the control node and finally proof that the global minimum corresponds to a value of α C = 0. Equation (34) can be rearranged as and thereby can be easily mapped to the solution of the minimization problem defined and solved in Eq.(30) and Eq.(31), respectively. Looking at Eq.(34), we should choose With this choice, the value of d that minimizes Eq.(34) corresponds to the weighted median of the set {s i } with weights α h /s i . Therefore, the value of d is the same as the solution of the case α C = 0, but d = s C and thus we must apply a transformation in order to obtain the optimal choice of s C . We have to distinguish several cases, considering α h > 0: • α C > 0. In this case we inspect Eq.(35) and distinguish two more cases: o α C > α h : in this case, the prefactor of s C is negative, and we can write: o α C < α h : in this case, the prefactor of s C is positive, and we can write: α h s C and considering that, in this case, 0 < α C < α h and hence 0 ≤ 1 − α C α h ≤ 1 and the weighted mean is bounded by min(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.05).
• α C < 0. In this case we can rewrite Eq.(35) as In this case, d = 1 + |α C | α h s C and hence 0 ≤ s C ≤ d/2. Hence, the optimal value of s C is always smaller than half the value of the weighted median, d (see Figure 3 at o |α C | < |α h |:n this case, the prefactor of s C is positive and bounded by 1 ≤ 1 + |α C | α h ≤ 2. In this case, d = 1 + |α C | α h s C and hence d/2 ≤ s C ≤ d. Hence, the optimal value of s C is always smaller than the weighted median, d (see Figure 3 at α C = −0.05).
• α C = 0: This case is explored in Section V A. Equation (35) turns to The optimal value of s C is the same as the weighted median, d, without any further transformation (see Figure 3 at α C = 0.0).
• α C = α h : This case is discussed in the introduction of the present section. Eq.(35) turns to and hence the value of the cost is the same constant value for all nodes (see Figure 3 at α C = 0.1).
Amid all the cases considered concerning the value of α C , the global minimum cost is given by α C = 0, as shown in Figure 3. This result can be proved by considering a simplified version of Eq.(35), defined as The minimum value of Eq.(36) is achieved when x = 0, as long as a > 0, b > 0 and c > 0. This conditions are equivalent to s i > 0, α h > 0 and s C > 0, and are true for all the summation terms in Eq.(35). Therefore, the minimum value is given by setting α C = 0. Summing up, in order to obtain the optimal {α i } parameters' set in order to achieve the symmetric phase configuration with the minimum implied cost in the Kuramoto-Sakaguchi model, we should set α C = 0, independently of the value of α h . The remaining parameters have to be tuned using Eq.(33). Moreover, the optimal choice of the control node (or nodes) corresponds to that with s C located at the weighted median of {s i } (with weight equal to s −1 i ). Notice also that nodes are grouped by degree regarding the tuned values of its frustration parameters. In other words, there may be different potential control nodes, as long as they share the same degree.

VI. FULLY SYNCHRONIZED PHASE CONFIGURATION
Another particular phase configuration is given by the phase synchronization of nodes, that is, φ * = 0. If we set, as in Section V, ω i = ω h ∀i, we end up with the trivial solution α i = 0 ∀i. In the case of full synchronization we want to recover the completely in-phase state from a phase dispersion produced by a distribution of natural frequencies, which we consider to be positive. Hence, applying Eq.(13) to this case: and in vector form, where we have used: ∆ ω =M ω.
Finally, from the κ in Eq.(38) we can obtain α: Similarly as the result of the symmetric configuration, given in Eq.(33), the solution of the fully synchronized configuration concerning α is a continuous spectrum of values, depending on the choice of the control node, C, the value of its frustration parameter α C , which is a free parameter, and the natural frequencies of the oscillators. In Sections VI A and VI B we will make a in-depth analysis of the problem, as well as comment on the nonlinear expansion of the Kuramoto-Sakaguchi model and the validity 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 general solution of the frustration parameters in Eq.(39) we get: In the particular choice α C = 0: Equation (41) shows that the relevant piece of information regarding the control node is given by its natural frequency, ω C . Similarly to the minimization problem posed in Section V, and in order to find the optimal choice of the control node we need to solve Eq. (20) considering the solution of Eq.(41): The optimization problem is equivalent to the most general problem, described in Eq.(30), with solution given by Eq. (31). In this case, d = ω C , x i = ω i and the weight w i = s −1 i . Accordingly, and in a similar way as in Section V, the solution of ω C corresponds to the weighted median of the set {ω i }, with weights given by the inverse of the node degree. Notice that the optimal choice of the control node is in general different to that given in Section V A). This is due to the fact that the weights of the weighted median have to be sorted according to descending order of natural frequencies instead of node degree.
Following with the example provided in Section V A, for the network in Figure 3(a), with degree sequence s = (1, 6, 2, 1, 2, 2, 2, 2), let us compute the optimal value of ω C by using Eq.(31). Consider the following natural frequencies and the corresponding weights To find the weighted median, we have to find the minimum value such that the sum of the weights at each side of the pivot are as even as possible. Therefore, the optimal value of natural frequency corresponds to the choice C = 6 [see α C = 0 line in Figure  4 The cost corresponding to the fully synchronized configuration case is given by Eq.(40). In the general case where α C = 0, we can minimize the cost with respect to ω C or to s C . If we minimize with respect to ω i , we first have to rewrite Eq.(40) as Again, the problem and the solution of Eq.(46) can be taken from Eq.(30) and Eq.(31), choosing d ≡ ω C − α C s C K, w i ≡ 1/Ks i and x i ≡ ω i . Hence, the value d that minimizes the cost is the weighed median considering the same weight as in Section VI A, w i = 1/s i (notice, however, that the ordering is determined by natural frequencies and not degrees). Let us analyze the different possibilities regarding the values of α C , maintaining ω C and s C constant: The value which minimizes cost is given by d = ω k , corresponding to the weighted median. However this 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, the value of the weighted median is d = 0.25. In the case we are considering, however, this is not the optimal choice of the parameters for the control node. We must shift the values considering the relation between d and the other parameters. If we choose α C = 0.1, for instance, we find that, |ω C − 0.1s c | = 0.25. In Figure 4 we see that the optimal choice is given by ω C = 0.4, which corresponds to C = 5 and s C = 2. • Hence, as the function increases with increasing (ω i − α C s i K), the minimum is achieved by min C (ω C − α C s C K).

C. Non-linear expansion of the Kuramoto-Sakaguchi model
The results obtained in Section VI are based on a linear approximation of the Kuramoto-Sakaguchi model. We have derived the results based on the phase synchronization requirement, and assuming that frequency synchronization is already achieved in the steady state. Nevertheless, when measuring the order parameter with a large dispersion of natural frequencies or low coupling constant, we do not expect such steady state. However, we ask to which extend the proposed values of the obtained frustration parameters are also able to enhance frequency synchronization considering the original nonlinear Kuramoto-Sakaguchi model: We compare the results from Ref. [21] considering its Type II frustration parameters tuning for both the linear and the nonlinear Kuramoto model and we find that, despite our approach does not consider the enhancement of frequency synchronization on the nonlinear regime, it is able to improve the value of the order parameter, in a similar fashion as in Ref. [21]. This work considers the nonlinear Kuramoto-Sakaguchi model and seeks to improve the number of nodes that fall into the recruitment condition so as to achieve the same common oscillatory frequency. The considered network class is the same as the mentioned paper, as well as the statistics study. We make use of the expression in Eq.(39) to tune the set of α for a given configuration of random ω and study the effect on the synchronization of the system for different values of the coupling strength.
We consider two cases: the linear and the nonlinear model with natural frequencies obtained from a uniform distribution ω i ∈ [−1, 1].
From Figure 5(a), the linear case of the Kuramoto-Sakaguchi model [see Eq.(4)], our approach, derived from the analytic expression of the linear approximation, advances the analytic tuning of frustration parameters suggested by Ref. [21]. This is because they look for an enhancement in the number of nodes that are oscillating at the same frequency, Ω, but they do not worry about the exact values of the phases they achieve. On the contrary, we assume nodes are already synchronized (without setting the specific value of Ω, as they do) and we look for the full synchronization state.
In Figure 5(b), the linear tuning squared discontinuous purple line) approaches the type II (squared dashed green line) tuning in the case of the nonlinear Kuramoto-Sakaguchi model, even for small values of the coupling strength. Hence, despite the aim of our approach is not achieving frequency synchronization, the obtained tuning of the frustration parameters helps enhancing it as well. In principle this behavior is reminiscent of the so called explosive percolation (see Ref. [24] and references therein), since the transition to the synchronized state is abrupt, as it happens in a first order phase transition. We are adjusting the phase-lag parameter as a response to the frequencies, and then in some sense it is similar to the original proposal in Ref. [25], the correlated degreefrequency framework.

VII. CONCLUSIONS
The Kuramoto-Sakaguchi model adds to the original Kuramoto model a homogeneous phase lag, α, between nodes which promotes a phase shift between oscillators. We consider a more general framework, in which the phase lag or the frustration parameter, α i , is an intrinsic property of each node. A very relevant question in oscillatory models is finding the conditions of network synchronization. In the present work, we bring forward a methodology not only to obtain the desired synchronized state, but any convenient phase configuration in the steady state, by means of a fine tuning of the phase lag or frustration parameters, {α i }. We feature the analytical solution of frustration parameters so as to achieve any phase configuration, by linearizing the most general model. The three intrinsic parameters of the nodes in the model, natural frequencies, {ω i }, frustration parameters {α i }, and phases in the steady state φ * i , are coupled by an equation that allows to tune them for a desired configuration. While the set φ * i is uniquely determined, the set α i has a continuous spectrum of solutions.
A main result is that a given phase configuration can be access via a continuous spectrum of frustration parameters, i.e, one phase and one frustration parameter are left as free parameters. The nodes we choose their values concerning phase and frustration parameter, are named reference and control nodes, respectively. Once the frustration parameters are tuned so as to obtain the desired state, we define a cost function to assess the overhead that the system requires to achieve such parameters' configuration. 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 for the cases of symmetric configuration and fully synchronized state and discuss them.
A key result is the solution to the minimization cost problem: For the case of symmetric configuration, the nodes which are to be set as control nodes are those whose degree is the weighed median of the sample, with a weight equal to the inverse of its degree. On the other hand, for the case of fully synchronized state, control nodes are those whose natural frequency is the weighted median of the sample, with a weight equal to the inverse of its degree. An extensive analysis of several cases is done in the text and a detailed example of a toy network is provided. We highlight the connection made with the nonlinear Kuramoto-Sakaguchi model. Despite our analysis being based on the linear version of the model, we show that the proposed parameters' tuning is also able to enhance frequency synchronization, as done in Ref. [21]. We stress the fact that the question 'among all the possible solutions, which is the one that makes the system achieve a particular phase configuration with the minimum required cost?' is of particular relevance when we consider the plausible real nature of the system. If a real system needs to access a particular phase configuration, which may be associated with a singular function, then it will tend to minimize the effort or cost to do so. Further work can be done within this framework by doing real experiments on measuring the energy needed to access a particular configuration. Moreover, other nonlinear oscillatory models can be analyzed and compared with the Kuramoto-Sakaguchi model.
Other questions regarding the model are left open. We have considered the coupled trio of natural frequenciesfrustration parameters-steady state phases. A natural extension to this would be to inspect the possibility to also tune the weights of the network edges in order to access a particular configuration. The higher dimension of the latter with respect to the vectors of parameters would require further assumptions about the model or the network structure, such as positive weights or particular distributions or topologies. Another research venue would be to consider the effect of removing a node of the network and the {α i } set needed to minimize the effect on the removal on the whole network.
Despite we provide the analytical solution to the optimal choice of parameters in order to minimize the cost of achieving both the symmetric and the fully synchronized configurations, the access to all nodes' parameters requirement may not be feasible in real-world networks. Our methodology is quite general and the optimization procedure refers to a set of parameters to be tuned. In particular, a finite subset of nodes with accessible phaselag parameter could be chosen (the choice could be restricted to any subset of nodes), holding all other nodes unaltered. This would provide a nonoptimal global condition but a restricted and approximated one that could deal with a subset of available nodes. A meaningful analysis would be to identify which subset of nodes is the one that enables to get a closest approximate solution, and relate those nodes with their topological properties, although this question is beyond the goal of this work. Step-by-step derivation of the example Consider the network in Figure 1, with its Laplacian matrix: We develop the equation step by step: Which leads to the result:

Appendix B: Mathematical solution of the cost optimization problem and intuitive insights
In order to gain a more intuitive understating of the analytical expression and solution of the considered cost function, we consider the analysis of the continuous case.

Symmetric configuration case
Considering the symmetric configuration case and choosing α C = 0, the continuous optimization problem can be written as which is depicted in Figure 7 for different values of a and the sum of all of them. Regardless of the set of a i values, the sum function i f (x, a i ) (see the example in the black line in Figure 7), is a concave function and has a unique minimum, which corresponds to one of the a i values.
In order to assess the value of a i where the minimum is located, we compute the derivative of Eq.(B2): and hence, d i f (x, a i )/dx = i sgn(x−a i )/a i , which is depicted in Figure 8. Notice that, despite the derivative of the function is not defined at the values x = a i , the derivative changes its sign when moving from x < 3 to x > 3 and hence, the minimum is located at this value of a i .
To conclude, Eq(B1) behaves equivalently as the function defined in Eq.(B2) and hence, displays only one minimum, which is achieved at the s i where there is a change of sign in the derivative.
Alternatively and as explained in the main text, we can understand the minimization problem as part of a general framework. The minimization of Eq.(B1) is equivalent to the minimization of the absolute value of the relative error: The general problem can be written as [23]: In other words, the d value that minimizes Eq.(B5) corresponds to the weighted median of the variable x, or the 50% weighted percentile.
Weighted median: For n distinct ordered elements x 1 , x 2 , ..., x n with positive weights w 1 , w 2 , ..., w n , the weighted median is the element x k satisfying min{i| i k=1 w k ≥ n k=i w k } Therefore, the solution is given by x k , the value such that 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 weighted medians with weights 1/s i , which are a special case of the medians of a measure.