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

Gemma Rosell-Tarragó<EMAIL_ADDRESS>Departament de Física de la
Matèria Condensada, Universitat de Barcelona, Martí Franquès 1, 08028
Barcelona, Spain  Universitat de Barcelona Institute of Complex Systems
(UBICS), Martí Franquès, 1, 08028 Barcelona, Spain  Albert Díaz-Guilera
<EMAIL_ADDRESS>Departament de Física de la Matèria Condensada, Universitat
de Barcelona, Martí Franquès 1, 08028 Barcelona, Spain  Universitat de
Barcelona Institute of Complex Systems (UBICS), Martí Franquès, 1, 08028
Barcelona, Spain

###### Abstract

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 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 $\pi/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 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, $\theta_{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, $\mathcal{G}(\mathcal{V},\mathcal{E})$. The
coupling strength describes the intensity of such pair-wise interactions,
$K_{ij}$. The set of nodes of the network, $\mathcal{V}(\mathcal{G})$ consists
of all the oscillators, while the set of edges, $\mathcal{E}(\mathcal{G})$, is
made of the links between them. In his original work, Kuramoto assumed
homogeneous interactions, i.e, $K_{ij}=K\ \forall(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]:

$\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j}A_{ij}\sin(\theta_{j}-\theta_{i})\
i=1,...,N\ j\in\Gamma_{i}$ (1)

where $\Gamma_{i}$ is the set of neighbors of node $i$ and $\omega_{i}$ is the
natural frequency of each unit.

We consider that two nodes are phase synchronized when their phases have the
same value,

$\theta_{i}(t)-\theta_{j}(t)=0\ \forall t>t_{0}$

When the phase difference has a constant value, that is,
$\theta_{i}(t)-\theta_{j}(t)=c\ \forall 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:

$\frac{d\theta_{i}}{dt}-\frac{d\theta_{j}}{dt}=0\ \forall t>t_{0}$

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:

$\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j}A_{ij}\sin(\theta_{j}-\theta_{i}-\alpha)\
i=1,...,N\ j\in\Gamma_{i}$ (2)

where $\alpha$ 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 $|\alpha|<\pi/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, $\omega_{i}=\omega_{0}\
\forall i$, the frustration parameter, $\alpha$, 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:

$\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j}A_{ij}\sin(\theta_{j}-\theta_{i}-\alpha_{i})\
i=1,...,N\ j\in\Gamma_{i}$ (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:

$\displaystyle\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j}W_{ij}(\theta_{j}-\theta_{i}-\alpha_{i})=$
$\displaystyle=\omega_{i}-K\sum_{j}L_{ij}\theta_{j}-K\alpha_{i}s_{i}$ (4)

where $W_{ij}$ is the value of the weight of the edge between node $i$ and
node $j$, $s_{i}\equiv\sum_{j}W_{ij}$ is the weighted degree of the $i$th node
and $L$ is the weighted Laplacian matrix defined as
$L_{ij}\equiv\delta_{ij}s_{i}-W_{ij}$. In the stable regime, a synchronized
frequency is achieved and, for all oscillators $\dot{\theta}_{i}=\Omega$. We
can derive the value of the common frequency oscillation, $\Omega$, summing
Eq.(III) over $i$:

$\sum_{i}\Omega=\sum_{i}\omega_{i}-K\sum_{i}\sum_{j}L_{ij}\theta_{j}^{*}-K\sum_{i}\alpha_{i}s_{i}$
(5)

Taking into account the steady state $\dot{\theta}_{i}=\Omega\ \forall i$ and
arranging summations:

$N\Omega=\sum_{i}\omega_{i}-K\sum_{j}\theta_{j}^{*}\sum_{i}L_{ij}-K\sum_{i}\alpha_{i}s_{i}.$
(6)

and finally,

$\Omega=\left\langle\omega\right\rangle-K\left\langle\alpha s\right\rangle.$
(7)

where we have used the Laplacian matrix property: $\sum_{i}L_{ij}=0$ and
defined the averages $\sum_{i}\alpha_{i}s_{i}/N=\left\langle\alpha
s\right\rangle$ and $\sum_{i}\omega_{i}/N=\left\langle\omega\right\rangle$.
Now we can plug expression Eq.(7) to Eq.(III) to get the stable phases of
oscillators, $\theta_{i}^{*}$:

$\sum_{j}L_{ij}\theta_{j}^{*}=\frac{\omega_{i}}{K}-\frac{\left\langle\omega\right\rangle}{K}+\left\langle\alpha
s\right\rangle-\alpha_{i}s_{i}\hskip 28.45274pt\forall i$ (8)

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
$\\{\theta_{i}^{*}\\}=f_{i}(t)$, but with the phase differences with respect
to a reference node, once the stationary state is achieved,

$\phi_{i}\equiv\theta_{i}-\theta_{R}$ (9)

In this way, we work with time independent values. In this situation,
$\phi_{R}=0$, by definition, as $\phi_{R}\equiv\theta_{R}-\theta_{R}=0$.

On the other hand, the contribution $\left\langle\alpha
s\right\rangle-\alpha_{i}s_{i}$ of the right-hand side of Eq.(8) can be
written in matrix form as:

$\displaystyle-\begin{pmatrix}\frac{N-1}{N}&-\frac{1}{N}&-\frac{1}{N}&...\\\
-\frac{1}{N}&\frac{N-1}{N}&-\frac{1}{N}&...\\\ ...&...&...&...\\\
-\frac{1}{N}&-\frac{1}{N}1&...&\frac{N-1}{N}\end{pmatrix}\cdot\begin{pmatrix}s_{0}&0&...&0\\\
0&s_{1}&...&0\\\ 0&...&...&0\\\ 0&...&0&s_{N-1}\end{pmatrix}\cdot$
$\displaystyle\cdot\begin{pmatrix}\alpha_{0}\\\ \alpha_{1}\\\ ...\\\
\alpha_{N-1}\end{pmatrix}=(-M\cdot D_{s})\vec{\alpha}$

where we have defined:
$M\equiv\begin{pmatrix}\frac{N-1}{N}&-\frac{1}{N}&-\frac{1}{N}&...\\\
-\frac{1}{N}&\frac{N-1}{N}&-\frac{1}{N}&...\\\ ...&...&...&...\\\
-\frac{1}{N}&-\frac{1}{N}1&...&\frac{N-1}{N}\end{pmatrix}$ and
$D_{s}\equiv\begin{pmatrix}s_{0}&0&...&0\\\ 0&s_{1}&...&0\\\ 0&...&...&0\\\
0&...&0&s_{N-1}\end{pmatrix}$. We write Eq.(8) in matrix as

$L\vec{\theta}^{*}=\frac{1}{K}\vec{\Delta\omega}-M\cdot D_{s}\vec{\alpha}$
(10)

where $\Delta\omega_{i}\equiv\omega_{i}-\left\langle\omega\right\rangle$.
Finally, we obtain the set of unknowns $\\{\alpha_{i}\\}$:

$M\cdot D_{s}\vec{\alpha}=\frac{1}{K}\vec{\Delta\omega}-L\vec{\theta}^{*}$
(11)

Equation (11), however, does not have a solution, because of the singular
nature of $M\cdot D_{s}$ matrix. $M$ matrix is singular too, and hence, its
inverse matrix does not exist. Mathematically, $det(M\cdot D_{s})=det(M)\cdot
det(D_{s})=0$

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:

$\kappa_{i}\equiv\alpha_{i}-\alpha_{C}$ (12)

where $\alpha_{C}$ is the value of the control node. In this situation,
$\kappa_{C}=0$, by definition, as $\kappa_{C}\equiv\alpha_{C}-\alpha_{C}=0$.

To easily write the matrix expressions, we define the selection matrix
$J_{(n,m)}$, which is, in general, an $(N-1)\times(N-1)$ identity matrix after
the removal of the $mth$ row and the $nth$ column.

$L\vec{\theta}^{*}$ turns to $\tilde{L}(k,R)\vec{\phi}^{*}$, 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,
$\tilde{L}(k,R)=J_{(,k)}\cdot L\cdot J_{(R,)}\equiv\tilde{L}$.

Similarly,
$\vec{\tilde{\phi}}(k)=J_{(,k)}\cdot\vec{\phi}\equiv\vec{\tilde{\phi}}$, where
we have removed the $kth$ row.

In an equivalent way as the definition of the reduced Laplacian:

$\tilde{MD_{s}}(k,C)=J_{(,k)}\cdot MD_{s}\cdot J_{(C,)}\equiv\tilde{MD_{s}}$

where $\tilde{MD_{s}}$ is $MD_{s}$ without the $kth$ row and the $Cth$ column.

Similarly,
$\vec{\tilde{\kappa}}(k)=J_{(,k)}\cdot\vec{\kappa}\equiv\vec{\tilde{\kappa}}$
and
$\vec{\tilde{\Delta\omega}}(k)=J_{(,k)}\cdot\vec{\Delta\omega}\equiv\vec{\tilde{\Delta\omega}}$,
where we have removed the $kth$ row.

Considering all the previous definitions and remarks, Eq.(10) can be rewritten
as:

$\displaystyle\tilde{MD_{s}}\vec{\tilde{\kappa}}=\frac{1}{K}\vec{\tilde{\Delta\omega}}-\tilde{L}\vec{\tilde{\phi}}^{*}-\alpha_{C}\cdot
J_{(,k)}\sum_{j}^{\rightarrow}[MD_{s}]_{ij}$ (13)

and finally,

$\displaystyle\vec{\tilde{\kappa}}=\left(\tilde{MD_{s}}\right)^{-1}\left(\frac{1}{K}\vec{\tilde{\Delta\omega}}-\tilde{L}\vec{\tilde{\phi}}^{*}-\alpha_{C}\cdot\tilde{M}\vec{s}\right)$
(14)

where we have used $J_{(,k)}\vec{\sum}_{j}[MD_{s}]_{ij}=\tilde{M}\vec{s}$.
Notice that $MD_{s}$ matrix is singular, but the row sum is not zero, although
it is so for the column sum. Hence, we need to set $\alpha_{C}=0$ if we want
to avoid extra constant arrays in the final expression. In this particular
case:

$\displaystyle\vec{\tilde{\kappa}}=\left(\tilde{MD_{s}}\right)^{-1}\left(\frac{1}{K}\vec{\tilde{\Delta\omega}}-\tilde{L}\vec{\phi}^{*}\right)$
(15)

and keep in mind that $\kappa_{C}=0$.

The obtained values of $\vec{\alpha}$ depend on both the chosen control node,
$C$, and the value we set for its frustration parameter, $\alpha_{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-row-sum equal to zero of
$MD_{s}$ matrix, the differences between the obtained values are dependent of
the control node choice. Mathematically,
$\alpha_{i}-\alpha_{j}(C=l)\neq\alpha_{i}-\alpha_{j}(C=k)$ if $l\neq k$. 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., $\omega_{i}=\omega\ \forall i$. In this case,
Eq.(14) turns to:

$\vec{\tilde{\kappa}}=\left(\tilde{MD_{s}}\right)^{-1}\left(-\tilde{L}\vec{\tilde{\phi}}^{*}-\alpha_{C}\cdot\tilde{M}\vec{s}\right)$

which in the case of the network depicted in Figure 1,

Figure 1: Network of seven nodes.

leads to the solution

$\begin{pmatrix}\kappa_{0}\\\ \kappa_{2}\\\ \kappa_{3}\\\ \kappa_{4}\\\
\kappa_{5}\\\
\kappa_{6}\end{pmatrix}=\begin{pmatrix}\frac{-2\alpha_{1}+3\phi_{1}+\phi_{3}+3\phi_{6}}{4}\\\
\frac{3(\phi_{1}-\phi_{2})}{2}\\\
\frac{2\phi_{1}-\phi_{2}-2\phi_{3}+\phi_{4}}{2}\\\
\frac{2\phi_{1}-\phi_{2}+\phi_{3}-2\phi_{4}+\phi_{5}}{2}\\\
\frac{2\phi_{1}-\phi_{2}+\phi_{4}-2\phi_{5}-\phi_{6}}{2}\\\
\frac{2\phi_{1}\phi_{2}+\phi_{5}-2\phi_{6}}{2}\end{pmatrix}$ (16)

where we have chosen $\kappa_{1}=0$ and $\phi_{0}=0$. Hence, the results are
written as a function of the value $\alpha_{1}$ and $\phi_{i}\ i\neq 0$.
Therefore, we can achieve any phase configuration, given by the set
$\\{\phi_{i}\\}$ by tuning the frustration parameters set $\\{\alpha\\}$,
where $\alpha_{i}=\kappa_{i}+\alpha_{C}$.

To illustrate how we obtain the final values, let us consider the following
phase configuration:

$\vec{\tilde{\phi}}_{(R=0)}=(0.1,0.2,0.25,-0.2,-0.1,0.0)$ (17)

In the general case where $\alpha_{C}=\alpha_{1}\neq 0$:

$\vec{\tilde{\kappa}}_{(C=1)}=(0.1375-\frac{\alpha_{1}}{2},-0.15,-0.35,0.275,0.0,-0.05).$

If we choose $\alpha_{C}=0$, then $\alpha_{i}=\kappa_{i}$, we can include the
value of the control node $C=1$:

$\vec{\tilde{\alpha}}=(0.1375,0.0,-0.15,-0.35,0.275,0.0,-0.05).$

Alternatively, we can choose whatever value we need regarding the control
node. For instance, if $\alpha_{C}=\alpha_{1}=0.1$:

$\vec{\tilde{\alpha}}=(0.1875,0.1,-0.05,-0.25,0.375,0.1,0.05)$

and the phases configuration is the same. Importantly, we recover the same
phase differences using the nonlinear model with the tuned $\alpha$’s, up to
an error. For this last example and using the frustration parameters obtained
by setting $\alpha_{1}=0.1$, the nonlinear model leads to final phases vector

$\begin{split}\vec{\tilde{\phi}}_{(R=0)}=(0.09969,0.19944,0.25097,\\\
-0.19798,-0.09897,0.00012)\end{split}$ (18)

which represents $\sim 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:

$e_{T}(C)\equiv\sum_{i}|\alpha_{i}(C)|$ (19)

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,
$\alpha_{C}$.

The optimal configuration is given by the solution of the minimization problem

$\min_{C,x}e_{T}(C,x)=\min_{C,x}\sum_{i}^{N}\alpha_{i}(C,x)$ (20)

where the $x$ variable is not yet specified. Depending on the problem we are
interested in we would set it either to $\omega_{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, $\alpha_{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 $\alpha_{C}$. In
Figure 2 we present an example, following with the network presented in
Section III and choosing different values of $\alpha_{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.

Figure 2: (Color online) Implied cost to achieve the phase configuration in
Eq.(17) as a function of the chosen control node, $C$, for the network in
Figure 1 and considering five different values of $\alpha_{C}$. Notice that
the minimum cost is given, in this case, by $\alpha_{C}=0$ and $C=1$ or $C=5$.

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, $\alpha_{i}=\alpha_{h}\ \forall i$. For our
purposes, we consider $\alpha_{h}>0$. In this situation, the trivial solution
of the frustration parameters, $\alpha_{i}=\alpha_{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.

###  V.1 Optimal cost tuning when $\alpha_{C}=0$

Let us firstly consider the case where $\alpha_{C}=0$ and homogeneous natural
frequencies $\omega_{i}=\omega_{h}\ \forall i$. In the particular case of the
symmetric configuration, that is, the phase configuration given by
$\alpha_{i}=\alpha_{h}\ \forall i$ the solution of the frustration parameters
is given by:

$\vec{\tilde{\kappa}}=\left(\tilde{MD_{s}}\right)^{-1}\left(\frac{1}{K}\vec{\tilde{\Delta\omega}}-\tilde{L}\vec{\tilde{\phi}}^{*}\right)=-\left(\tilde{MD_{s}}\right)^{-1}\tilde{L}\vec{\tilde{\phi}}^{*}$
(21)

But $\vec{\tilde{\phi}}^{*}$ corresponds to the symmetric case. Hence (see
Section II),

$\vec{\tilde{\phi}}^{*}=\alpha_{h}\tilde{L}^{-1}\vec{\tilde{\Delta s}}$ (22)

where $\vec{\tilde{\Delta s}}_{i}\equiv\left\langle s\right\rangle-s_{i}$ and
the tilde touches on $kth$ row removal.

Plugging Eq.(22) into Eq.(21):

$\vec{\tilde{\kappa}}=-\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\tilde{L}\tilde{L}^{-1}\vec{\tilde{\Delta
s}}=-\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\vec{\tilde{\Delta s}}$

But $\vec{\tilde{\Delta s}}$ can be written as:

$\vec{\tilde{\Delta s}}=-\tilde{M}\vec{s}$ (23)

Putting it all together:

$\vec{\tilde{\kappa}}=-\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\tilde{L}\tilde{L}^{-1}\vec{\tilde{\Delta
s}}=\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\tilde{M}\vec{s}\\\ \
\vec{\tilde{\kappa}}$ (24)

which in vector form is written as:

$\vec{\tilde{\kappa}}=\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\tilde{M}\vec{s}\\\
\ \vec{\tilde{\kappa}}=\alpha_{h}\begin{pmatrix}1-\frac{s_{C}}{s_{0}}\\\
1-\frac{s_{C}}{s_{1}}\\\ \cdots\\\ 1-\frac{s_{C}}{s_{N-1}}\\\ \end{pmatrix}$
(25)

And considering the relation between $\alpha$ and $\kappa$, in Eq.(12):

$\begin{split}\vec{\alpha}=\alpha_{h}\begin{pmatrix}1-\frac{s_{C}}{s_{0}}\\\
1-\frac{s_{C}}{s_{1}}\\\ \cdots\\\ 0\text{ ($C$ node)}\\\ \cdots\\\
1-\frac{s_{C}}{s_{N-1}}\\\ \end{pmatrix}\end{split}$ (26)

Equation (25) gives us the tuned values of the frustration parameters as a
function of the chosen control node, $C$, when $\alpha_{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):

$e_{T}(C)=\left\lvert\alpha_{h}\right\rvert\sum_{i}^{N-1}\left\lvert
1-\frac{s_{C}}{s_{i}}\right\rvert=|\alpha_{h}|\sum_{i}^{N-1}\Big{\lvert}\frac{s_{i}-s_{C}}{s_{i}}\Big{\rvert}$
(27)

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}\ll s_{i}$ or
$s_{C}\gg 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.

For example, the network in Figure 3(a), with $\vec{s}=(1,6,2,1,2,2,2,2)$ has
the set of unique degrees $\vec{s}_{unique}=(1,2,6)$ and hence three possible
values of the cost, shared by some nodes. If $C=\\{0,3\\}$, $s_{C}=1$:

$\displaystyle
e_{T}(C)=|\alpha_{h}|\left(|1-\frac{1}{1}|+5|1-\frac{1}{2}|+|1-\frac{1}{6}|\right)=$
$\displaystyle=|\alpha_{h}|\left(\frac{5}{2}+\frac{5}{6}\right)=\frac{10}{3}|\alpha_{h}|$

If $C=\\{2,4,5,6,7\\}$, $s_{C}=2$:

$\displaystyle
e_{T}(C)=|\alpha_{h}|\left(2|1-\frac{2}{1}|+4|1-\frac{2}{2}|+|1-\frac{2}{6}|\right)=$
$\displaystyle=|\alpha_{h}|\left(2+\frac{2}{3}\right)=\frac{8}{3}|\alpha_{h}|$

And, finally, if $C=1$, $s_{C}=6$:

$\displaystyle
e_{T}(C)=|\alpha_{h}|\left(2|1-\frac{6}{1}|+5|1-\frac{6}{2}|\right)=$
$\displaystyle=|\alpha_{h}|\left(10+10\right)=20|\alpha_{h}|$

The minimum value of the energy is $\frac{8}{3}|\alpha_{h}|$, corresponding to
the choice $C\in\\{2,4,5,6,7\\}$ with $s_{C}=2$.

Notice that the optimal choice of the control node (or nodes) does not depend
on the value of $\alpha_{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):

$\min_{C,x}e_{T}(C,x)=\min_{C,x}\sum_{i}^{N}\alpha_{i}(C,x)$

which, when considering the symmetric configuration case, turns to

$\min_{s_{C}}|\alpha_{h}|\sum_{i}^{N-1}\Big{\lvert}1-\frac{s_{C}}{s_{i}}\Big{\rvert}=|\alpha_{h}|\min_{s_{C}}\sum_{i}^{N}\Big{\lvert}\frac{s_{C}-s_{i}}{s_{i}}\Big{\rvert}$
(28)

Equation (28) is equivalent to the minimization of the absolute value of the
relative error of the degree:

$|\alpha_{h}|\min_{s_{C}}\sum_{i}^{N}|\mathcal{E}_{i}|$ (29)

where
$\mathcal{E}_{i}=\Big{\lvert}\displaystyle\frac{s_{C}-s_{i}}{s_{i}}\Big{\rvert}$.

The most general minimization problem of the relative error of a variable [23]
can be written as

$\min_{d}\sum_{i=1}^{N}w_{i}|x_{i}-d|\ ;d>0$ (30)

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

$\displaystyle d=x_{m}\text{ , where }m\equiv\min\\{i\ |\
\sum_{k=1}^{i}w_{k}\geq\sum_{k=i}^{n}w_{k}\\}$ $\displaystyle
i\in\\{1,...,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 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|\sum_{k=1}^{i}w_{k}\geq\sum_{k=i}^{n}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
$\vec{s}=(1,6,2,1,2,2,2,2)$, let us compute the optimal value of $s_{C}$ by
using Eq.(V.1).

$\displaystyle\text{sorted}(\vec{s})=(1,1,2,2,2,2,2,6)$
$\displaystyle\vec{w}=\left(1,1,\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{6}\right)$

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.

$1+1+\frac{1}{2}+\frac{1}{2}=3\geq
2.17=\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+\frac{1}{6}$

which corresponds to $s_{C}=2$, in agreement with the location of the minimum
for $\alpha_{C}=0$ in Figure 3(b) corresponding to the network in Figure 3(a).

Figure 3: (Color online) Implied cost to achieve the symmetric configuration
as a function of the degree corresponding to different choices of the control
node, $s_{C}$, for a network of 8 nodes (upper panels) and 9 nodes (lower
panels). The distinct colors and markers correspond to different values of
$\alpha_{C}$. The symmetric configuration is generated by a value of
$\alpha_{h}=0.1$.

###  V.2 Optimal cost tuning when $\alpha_{C}\neq 0$

We next ask which is the optimal choice of the control node in the case we let
$\alpha_{C}\neq 0$ and $\omega_{i}=\omega_{h}\ \forall i$. In this case, we
should look at Eq.(13) and set $\vec{\tilde{\Delta\omega}}=0$. Making use of
the analytical solution of the symmetric configuration in Eq.(22):

$\tilde{MD_{s}}\vec{\tilde{\kappa}}=-\alpha_{h}\tilde{L}\tilde{L}^{-1}\vec{\tilde{\Delta
s}}-\alpha_{C}\cdot\sum_{j}^{\rightarrow}[MD_{s}]_{ij}$

Using the properties $\tilde{L}\tilde{L}^{-1}=\mathbb{I}$ and
$\vec{\tilde{\Delta s}}=\tilde{M}\vec{s}$,

$\vec{\tilde{\kappa}}=\alpha_{h}\left(\tilde{MD_{s}}\right)^{-1}\left(\tilde{M}\vec{s}-\alpha_{C}\left(\tilde{MD_{s}}\right)^{-1}\sum_{j}^{\rightarrow}[MD_{s}]_{ij}\right)$

Finally, in vector form,

$\displaystyle\vec{\tilde{\kappa}}=\alpha_{h}\begin{pmatrix}1-\frac{s_{C}}{s_{0}}\\\
1-\frac{s_{C}}{s_{1}}\\\ \cdots\\\ 1-\frac{s_{C}}{s_{N-1}}\\\
\end{pmatrix}-\alpha_{C}\begin{pmatrix}1-\frac{s_{C}}{s_{0}}\\\
1-\frac{s_{C}}{s_{1}}\\\ \cdots\\\ 1-\frac{s_{C}}{s_{N-1}}\\\ \end{pmatrix}$
(32)

Using Eq.(12),

$\displaystyle\vec{\alpha}=(\alpha_{h}-\alpha_{C})\begin{pmatrix}1-\frac{s_{C}}{s_{0}}\\\
1-\frac{s_{C}}{s_{1}}\\\ \cdots\\\ 0\text{ ($C$ node)}\\\ \cdots\\\
1-\frac{s_{C}}{s_{N-1}}\\\ \end{pmatrix}+\alpha_{C}$ (33)

where we have used the result in Eq.(25) and the relation

$\displaystyle\left(\tilde{MD_{s}}\right)^{-1}\sum_{j}^{\rightarrow}[MD_{s}]_{ij}=\left(\tilde{MD_{s}}\right)^{-1}\tilde{M}\vec{s}=$
$\displaystyle=\begin{pmatrix}1-\frac{\alpha_{C}}{\alpha_{0}}\\\
1-\frac{\alpha_{C}}{\alpha_{1}}\\\ \cdots\\\
1-\frac{\alpha_{C}}{\alpha_{N-1}}\\\ \end{pmatrix}$

In the particular case that $\alpha_{C}=\alpha_{h}$ we recover the trivial
initial configuration $\alpha_{i}=\alpha_{h}\ \forall 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):

$e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}(\alpha_{h}-\alpha_{C})\left(1-\frac{s_{C}}{s_{i}}\right)+\alpha_{C}\Big{\rvert}$
(34)

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
$\alpha_{C}=0$. Equation (34) can be rearranged as

$\displaystyle
e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{s_{i}-\left(1-\frac{\alpha_{C}}{\alpha_{h}}\right)s_{C}}{s_{i}/\alpha_{h}}\Big{\rvert}$
(35)

and thereby can be easily mapped to the solution of the minimization problem
defined and solved in Eq.(30) and Eq.(V.1), respectively. Looking at Eq.(34),
we should choose $x_{i}=s_{i}$, $w_{i}=\alpha_{h}/s_{i}$ and
$d=\left(1-\alpha_{C}/\alpha_{h}\right)s_{C}$. With this choice, the value of
$d$ that minimizes Eq.(34) corresponds to the weighted median of the set
$\\{s_{i}\\}$ with weights $\alpha_{h}/s_{i}$. Therefore, the value of $d$ is
the same as the solution of the case $\alpha_{C}=0$, but $d\neq 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 $\alpha_{h}>0$:

  * •

$\alpha_{C}>0$. In this case we inspect Eq.(35) and distinguish two more
cases:

    * o

$\alpha_{C}>\alpha_{h}$: in this case, the prefactor of $s_{C}$ is negative,
and we can write:

$\displaystyle
e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{s_{i}+\Big{\lvert}1-\frac{\alpha_{C}}{\alpha_{h}}\Big{\rvert}s_{C}}{s_{i}/\alpha_{h}}\Big{\rvert}=$
$\displaystyle=\sum_{i=0}^{N-1}\Big{\lvert}\alpha_{h}+M\alpha_{h}\frac{s_{C}}{s_{i}}\Big{\rvert}$

where $M\equiv\Big{\lvert}1-\frac{\alpha_{C}}{\alpha_{h}}\Big{\rvert}>0$ is a
positive number. Hence, as the cost function increases with increasing
$s_{C}$, the minimum is achieved when $s_{C}=\min(s_{i})$ (See Figure 3 at
$\alpha_{C}=0.2$).

    * o

$\alpha_{C}<\alpha_{h}$: in this case, the prefactor of $s_{C}$ is positive,
and we can write:

$e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{s_{i}-\Big{\lvert}1-\frac{\alpha_{C}}{\alpha_{h}}\Big{\rvert}s_{C}}{s_{i}/\alpha_{h}}\Big{\rvert}$

taking into account that
$d=\Big{\lvert}1-\frac{\alpha_{C}}{\alpha_{h}}\Big{\rvert}s_{C}$ and
considering that, in this case, $0<\alpha_{C}<\alpha_{h}$ and hence
$0\leq\Big{\lvert}1-\frac{\alpha_{C}}{\alpha_{h}}\Big{\rvert}\leq 1$ and the
weighted mean is bounded by $\min(s_{i})\leq d\leq\max(s_{i})$, the optimal
value of $s_{C}$ falls in the range $d\leq s_{C}\leq\max(s_{i})$. Hence, the
optimal value of $s_{C}$ is always larger than the weighted median, $d$ (see
Figure 3 at $\alpha_{C}=0.05$).

  * •

$\alpha_{C}<0$. In this case we can rewrite Eq.(35) as

$\displaystyle
e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{s_{i}-\left(1+\frac{|\alpha_{C}|}{\alpha_{h}}\right)s_{C}}{s_{i}/\alpha_{h}}\Big{\rvert}$

and distinguish two more cases:

    * o

$|\alpha_{C}|>|\alpha_{h}|$: In this case, the prefactor of $s_{C}$ is
positive and bounded by
$2\leq\left(1+\frac{|\alpha_{C}|}{\alpha_{h}}\right)\leq\infty$. In this case,
$d=\left(1+\frac{|\alpha_{C}|}{\alpha_{h}}\right)s_{C}$ and hence $0\leq
s_{C}\leq 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
$\alpha_{C}=-0.2$).

    * o

$|\alpha_{C}|<|\alpha_{h}|$:n this case, the prefactor of $s_{C}$ is positive
and bounded by $1\leq\left(1+\frac{|\alpha_{C}|}{\alpha_{h}}\right)\leq 2$. In
this case, $d=\left(1+\frac{|\alpha_{C}|}{\alpha_{h}}\right)s_{C}$ and hence
$d/2\leq s_{C}\leq d$. Hence, the optimal value of $s_{C}$ is always smaller
than the weighted median, $d$ (see Figure 3 at $\alpha_{C}=-0.05$).

  * •

$\alpha_{C}=0$: This case is explored in Section V.1. Equation (35) turns to

$\displaystyle
e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{s_{i}-s_{C}}{s_{i}/\alpha_{h}}\Big{\rvert}$

. The optimal value of $s_{C}$ is the same as the weighted median, $d$,
without any further transformation (see Figure 3 at $\alpha_{C}=0.0$).

  * •

$\alpha_{C}=\alpha_{h}$: This case is discussed in the introduction of the
present section. Eq.(35) turns to

$e_{T}(C)=\sum_{i=0}^{N-1}\alpha_{h}=N\alpha_{h}$

and hence the value of the cost is the same constant value for all nodes (see
Figure 3 at $\alpha_{C}=0.1$).

Amid all the cases considered concerning the value of $\alpha_{C}$, the global
minimum cost is given by $\alpha_{C}=0$, as shown in Figure 3. This result can
be proved by considering a simplified version of Eq.(35), defined as

$f(x)=\Big{\lvert}\frac{a-(1-x/b)c}{a/b}\Big{\rvert}$ (36)

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$, $\alpha_{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 $\alpha_{C}=0$.

Summing up, in order to obtain the optimal $\\{\alpha_{i}\\}$ parameters’ set
in order to achieve the symmetric phase configuration with the minimum implied
cost in the Kuramoto-Sakaguchi model, we should set $\alpha_{C}=0$,
independently of the value of $\alpha_{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_{i}^{-1}$).

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, $\vec{\phi}^{*}=\vec{0}$. If we set, as in Section V,
$\omega_{i}=\omega_{h}\ \forall i$, we end up with the trivial solution
$\alpha_{i}=0\ \forall 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:

$\displaystyle\vec{\tilde{\kappa}}=\left(\tilde{MD_{s}}\right)^{-1}\left(\frac{1}{K}\tilde{M}\vec{\omega}-\alpha_{C}\cdot\tilde{M}\vec{s}\right)$
(37)

and in vector form,

$\displaystyle\vec{\tilde{\kappa}}=\begin{pmatrix}\frac{\alpha_{C}(s_{C}-s_{0})-(\omega_{C}-\omega_{0})/K}{s_{0}}\\\
\frac{\alpha_{C}(s_{C}-s_{1})-(\omega_{C}-\omega_{1})/K}{s_{1}}\\\ \cdots\\\
\frac{\alpha_{C}(s_{C}-s_{N-1})-(\omega_{C}-\omega_{N-1})/K}{s_{N-1}}\end{pmatrix}\hskip
14.22636pt$ (38)

where we have used: $\vec{\tilde{\Delta\omega}}=\tilde{M}\vec{\omega}$.

Finally, from the $\vec{\kappa}$ in Eq.(38) we can obtain $\vec{\alpha}$:

$\vec{\alpha}=\begin{pmatrix}\frac{\alpha_{C}s_{C}-(\omega_{C}-\omega_{0})/K}{s_{0}}\\\
\frac{\alpha_{C}s_{C}-(\omega_{C}-\omega_{1})/K}{s_{1}}\\\ \cdots\\\
\alpha_{C}\\\ \cdots\\\
\frac{\alpha_{C}s_{C}-(\omega_{C}-\omega_{N-1})/K}{s_{N-1}}\end{pmatrix}$ (39)

Similarly as the result of the symmetric configuration, given in Eq.(33), the
solution of the fully synchronized configuration concerning $\vec{\alpha}$ is
a continuous spectrum of values, depending on the choice of the control node,
$C$, the value of its frustration parameter $\alpha_{C}$, which is a free
parameter, and the natural frequencies of the oscillators. In Sections VI.1
and VI.2 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.3).

###  VI.1 Optimal cost tuning when $\alpha_{C}=0$

Using the definition of cost in Eq.(19) and the general solution of the
frustration parameters in Eq.(39) we get:

$e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{\alpha_{C}s_{C}-(\omega_{C}-\omega_{i})/K}{s_{i}}\Big{\rvert}$
(40)

In the particular choice $\alpha_{C}=0$:

$e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{\omega_{C}-\omega_{i}}{Ks_{i}}\Big{\rvert}$
(41)

Equation (41) shows that the relevant piece of information regarding the
control node is given by its natural frequency, $\omega_{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):

$\min_{\omega_{C}}\Big{\lvert}\frac{\omega_{C}-\omega_{i}}{Ks_{i}}\Big{\rvert}=\frac{1}{K}\min_{\omega_{C}}\Big{\lvert}\frac{\omega_{C}-\omega_{i}}{s_{i}}\Big{\rvert}$
(42)

The optimization problem is equivalent to the most general problem, described
in Eq.(30), with solution given by Eq.(V.1). In this case, $d=\omega_{C}$,
$x_{i}=\omega_{i}$ and the weight $w_{i}=s_{i}^{-1}$. Accordingly, and in a
similar way as in Section V, the solution of $\omega_{C}$ corresponds to the
weighted median of the set $\\{\omega_{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.1). 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.1, for the network in Figure
3(a), with degree sequence $\vec{s}=(1,6,2,1,2,2,2,2)$, let us compute the
optimal value of $\omega_{C}$ by using Eq.(V.1). Consider the following
natural frequencies

$\vec{\omega}=(0.1,0.2,0.05,0.45,0.3,0.4,0.25,0.15)$ (43)

which lead to

$\text{sorted}(\vec{\omega})=(0.05,0.1,0.15,0.2,0.25,0.3,0.4,0.45)$ (44)

and the corresponding weights

$\vec{w}=\left(\frac{1}{2},1,\frac{1}{2},\frac{1}{6},\frac{1}{2},\frac{1}{2},\frac{1}{2},1\right)$
(45)

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.

$\frac{1}{2}+1+\frac{1}{2}+\frac{1}{6}+\frac{1}{2}=2.67\geq
2.5=\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+1$

Therefore, the optimal value of natural frequency corresponds to the choice
$C=6$ [see $\alpha_{C}=0$ line in Figure 4(a)], with $\omega_{C}=0.25$ [see
$\alpha_{C}=0$ line in Figure 4(b)] and a degree of $s_{C}=2$.

Figure 4: (Color online) Implied cost to achieve the fully synchronized
configuration as a function of the chosen control node, $C$ (upper panel) and
natural frequencies of nodes (bottom panel) for the network in Figure 3(a).
Five different values of $\alpha_{C}$ are considered (marked colored lines).
Natural frequencies are set as the example in Eq.(43).

###  VI.2 Optimal cost tuning when $\alpha_{C}\neq 0$

The cost corresponding to the fully synchronized configuration case is given
by Eq.(40). In the general case where $\alpha_{C}\neq 0$, we can minimize the
cost with respect to $\omega_{C}$ or to $s_{C}$. If we minimize with respect
to $\omega_{i}$, we first have to rewrite Eq.(40) as

$\displaystyle
e_{T}(C)=\sum_{i=0}^{N-1}\Big{\lvert}\frac{\alpha_{C}s_{C}-(\omega_{C}-\omega_{i})/K}{s_{i}}\Big{\rvert}$
$\displaystyle\frac{1}{K}\sum_{i=0}^{N-1}\Big{\lvert}\frac{\omega_{i}-(\omega_{C}-\alpha_{C}s_{C}K)}{s_{i}}\Big{\rvert}$
(46)

Again, the problem and the solution of Eq.(VI.2) can be taken from Eq.(30) and
Eq.(V.1), choosing $d\equiv\omega_{C}-\alpha_{C}s_{C}K$, $w_{i}\equiv
1/Ks_{i}$ and $x_{i}\equiv\omega_{i}$.

Hence, the value $d$ that minimizes the cost is the weighed median considering
the same weight as in Section VI.1, $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 $\alpha_{C}$, maintaining
$\omega_{C}$ and $s_{C}$ constant:

  * •

$\omega_{C}>\alpha_{C}s_{C}K$ or $\alpha_{C}<\frac{\omega_{C}}{Ks_{C}}$: We
can write

$\displaystyle\sum_{i}^{N}\Big{\lvert}\frac{\omega_{i}-|\omega_{C}-\alpha_{C}s_{C}K|}{Ks_{i}}\Big{\rvert}$

The value which minimizes cost is given by $d=\omega_{k}$, corresponding to
the weighted median. However this is not directly the value of $\omega_{C}$,
as $d=|\omega_{C}-\alpha_{C}s_{C}K|$ in this case. The real values of the pair
$\\{\omega_{C},s_{C}\\}$ are given by
$\min_{C}(\omega_{k}-(\omega_{C}-\alpha_{C}s_{C}K))$. Following with the
example in Section VI.1, 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 $\alpha_{C}=0.1$,
for instance, we find that, $|\omega_{C}-0.1s_{c}|=0.25$. In Figure 4 we see
that the optimal choice is given by $\omega_{C}=0.4$, which corresponds to
$C=5$ and $s_{C}=2$.

  * •

$\omega_{C}<\alpha_{C}s_{C}K$ or $\alpha_{C}>\frac{\omega_{C}}{Ks_{C}}$: We
can write

$\displaystyle\sum_{i}^{N}\Big{\lvert}\frac{\omega_{i}+|\omega_{C}-\alpha_{C}s_{C}K|}{Ks_{i}}\Big{\rvert}$

Hence, as the function increases with increasing
$(\omega_{i}-\alpha_{C}s_{i}K)$, the minimum is achieved by
$\min_{C}(\omega_{C}-\alpha_{C}s_{C}K)$.

###  VI.3 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:

$\dot{\theta}_{i}=\omega_{i}+K\sum_{j}W_{ij}\sin(\theta_{j}-\theta_{i}-\alpha_{i})$
(47)

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 $\vec{\alpha}$ for
a given configuration of random $\vec{\omega}$ 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 $\omega_{i}\in[-1,1]$.

Figure 5: (Color online) Average order parameter, $\left\langle
r\right\rangle$, as a function 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 homogeneous node degree $s_{i}=4$ and $N=100$
nodes. Natural frequencies are obtained from a uniform random distribution in
the range $\omega_{i}\in[-1,1]$. Each data point represents an average over
ten optimized configurations. We compare three types of tuning for the set of
frustration parameters, $\\{\alpha_{i}\\}$: the original Kuramoto dynamics or
$\alpha_{i}=0\ \forall i$ (spotted continuous red line); type II [21] KS
dynamics with the frustration parameters set to
$\sin(\alpha_{i})=-\omega_{i}/(Ks_{i})$ if $|\omega_{i}|<Ks_{i}$ and
$\sin(\alpha_{i})=\pm 1$ otherwise (squared dashed green line); and KS
dynamics with the frustration parameters determined by the derived linear
approximation in Eq.(39))(squared discontinuous purple line).

From Figure 5(a), the linear case of the Kuramoto-Sakaguchi model [see
Eq.(III)], 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, $\Omega$, 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 $\Omega$, 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 degree-frequency framework.

##  VII Conclusions

The Kuramoto-Sakaguchi model adds to the original Kuramoto model a homogeneous
phase lag, $\alpha$, between nodes which promotes a phase shift between
oscillators. We consider a more general framework, in which the phase lag or
the frustration parameter, $\alpha_{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, $\\{\alpha_{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, $\\{\omega_{i}\\}$,
frustration parameters $\\{\alpha_{i}\\}$, and phases in the steady state
$\phi^{*}_{i}$, are coupled by an equation that allows to tune them for a
desired configuration. While the set $\phi^{*}_{i}$ is uniquely determined,
the set $\alpha_{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 $\\{\alpha_{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 frequencies-frustration 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 $\\{\alpha_{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 phase-lag 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.

###### Acknowledgements.

The authors acknowledge financial support from MINECO via Project No.
PGC2018-094754-B-C22 (MINECO/FEDER,UE) and Generalitat de Catalunya via Grant
No. 2017SGR341. G.R.-T. also acknowledges MECD for Grant No. FPU15/03053

##  Appendix A Step-by-step derivation of the example

Consider the network in Figure 1, with its Laplacian matrix:

$L=\begin{pmatrix}4&-1&-1&-1&0&0&-1\\\ -1&2&-1&0&0&0&0\\\ -1&-1&2&0&0&0&0\\\
-1&0&0&2&-1&0&0\\\ 0&0&0&-1&2&-1&0\\\ 0&0&0&0&-1&2&1\\\ -1&0&0&0&0&-1&2\\\
\end{pmatrix}$

We develop the equation step by step:

$\displaystyle\sum_{j}L_{ij}\theta_{j}^{*}=\frac{\omega_{i}}{K}-\frac{\left\langle\omega\right\rangle}{K}+\left\langle\alpha
s\right\rangle-\alpha_{i}s_{i}\hskip 28.45274pt\forall i$
$\displaystyle\begin{pmatrix}4&-1&-1&-1&0&0&-1\\\ -1&2&-1&0&0&0&0\\\
-1&-1&2&0&0&0&0\\\ -1&0&0&2&-1&0&0\\\ 0&0&0&-1&2&-1&0\\\ 0&0&0&0&-1&2&1\\\
-1&0&0&0&0&-1&2\\\ \end{pmatrix}\cdot\begin{pmatrix}\theta_{0}^{*}\\\
\theta_{1}^{*}\\\ \theta_{2}^{*}\\\ \theta_{3}^{*}\\\ \theta_{4}^{*}\\\
\theta_{5}^{*}\\\ \theta_{6}^{*}\end{pmatrix}=$
$\displaystyle\begin{pmatrix}\frac{\omega_{0}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{1}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{2}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{3}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{4}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{5}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{6}-\left\langle\omega\right\rangle}{K}\\\
\end{pmatrix}-\begin{pmatrix}\frac{6}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&\frac{6}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&-\frac{1}{7}&\frac{6}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&\frac{6}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&\frac{6}{7}&-\frac{1}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&\frac{6}{7}&-\frac{1}{7}\\\
-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&-\frac{1}{7}&\frac{6}{7}\end{pmatrix}\cdot$
$\displaystyle\cdot\begin{pmatrix}4&0&0&0&0&0&0\\\ 0&2&0&0&0&0&0\\\
0&0&2&0&0&0&0\\\ 0&0&0&2&0&0&0\\\ 0&0&0&0&2&0&0\\\ 0&0&0&0&0&2&0\\\
0&0&0&0&0&0&2\\\ \end{pmatrix}\cdot\begin{pmatrix}\alpha_{0}\\\ \alpha_{1}\\\
\alpha_{2}\\\ \alpha_{3}\\\ \alpha_{4}\\\ \alpha_{5}\\\ \alpha_{6}\\\
\end{pmatrix}=\begin{pmatrix}\frac{\omega_{0}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{1}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{2}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{3}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{4}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{5}-\left\langle\omega\right\rangle}{K}\\\
\frac{\omega_{6}-\left\langle\omega\right\rangle}{K}\\\ \end{pmatrix}+$
$\displaystyle\begin{pmatrix}-\frac{24}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}\end{pmatrix}\cdot\begin{pmatrix}\alpha_{0}\\\
\alpha_{1}\\\ \alpha_{2}\\\ \alpha_{3}\\\ \alpha_{4}\\\ \alpha_{5}\\\
\alpha_{6}\\\ \end{pmatrix}$

If we set all natural frequencies to the same value: $\omega_{i}=\omega\
\forall i$:

$\displaystyle\begin{pmatrix}4&-1&-1&-1&0&0&-1\\\ -1&2&-1&0&0&0&0\\\
-1&-1&2&0&0&0&0\\\ -1&0&0&2&-1&0&0\\\ 0&0&0&-1&2&-1&0\\\ 0&0&0&0&-1&2&1\\\
-1&0&0&0&0&-1&2\\\ \end{pmatrix}\cdot\begin{pmatrix}\theta_{0}^{*}\\\
\theta_{1}^{*}\\\ \theta_{2}^{*}\\\ \theta_{3}^{*}\\\ \theta_{4}^{*}\\\
\theta_{5}^{*}\\\ \theta_{6}^{*}\end{pmatrix}=$
$\displaystyle=\begin{pmatrix}-\frac{24}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}\end{pmatrix}\cdot\begin{pmatrix}\alpha_{0}\\\
\alpha_{1}\\\ \alpha_{2}\\\ \alpha_{3}\\\ \alpha_{4}\\\ \alpha_{5}\\\
\alpha_{6}\\\ \end{pmatrix}$

Now we choose $R=0$ and $C=1$, i.e., all $\phi_{i}=\theta_{i}-\theta_{0}$ and
$\kappa_{i}=\alpha_{i}-\alpha_{1}$:

Let us write explicitly the change of variables (the red and the blue columns
are ones we can remove due to the change of variables, as they do not affect
to the system of equations anymore):

$\displaystyle\begin{pmatrix}\color[rgb]{1,0,0}4&-1&-1&-1&0&0&-1\\\
\color[rgb]{1,0,0}-1&2&-1&0&0&0&0\\\ \color[rgb]{1,0,0}-1&-1&2&0&0&0&0\\\
\color[rgb]{1,0,0}-1&0&0&2&-1&0&0\\\ \color[rgb]{1,0,0}0&0&0&-1&2&-1&0\\\
\color[rgb]{1,0,0}0&0&0&0&-1&2&1\\\ \color[rgb]{1,0,0}-1&0&0&0&0&-1&2\\\
\end{pmatrix}\cdot\begin{pmatrix}\color[rgb]{1,0,0}\phi_{0}^{*}=0\\\
\phi_{1}^{*}\\\ \phi_{2}^{*}\\\ \phi_{3}^{*}\\\ \phi_{4}^{*}\\\
\phi_{5}^{*}\\\ \phi_{6}^{*}\end{pmatrix}=$
$\displaystyle=\begin{pmatrix}-\frac{24}{7}&\color[rgb]{0,0,1}\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}\\\
\frac{4}{7}&\color[rgb]{0,0,1}\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}\end{pmatrix}\cdot$
$\displaystyle\cdot\begin{pmatrix}\kappa_{0}\\\
\color[rgb]{0,0,1}\kappa_{1}=0\\\ \kappa_{2}\\\ \kappa_{3}\\\ \kappa_{4}\\\
\kappa_{5}\\\ \kappa_{6}\\\
\end{pmatrix}+\begin{pmatrix}\frac{-12}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\
\end{pmatrix}$

If we look carefully at Eq.(A), we see that although the left-hand side and
the right hand-side matrices are both singular, the first one has both column
and row sums equal to zero, while the second one has only column sum equal to
zero. This is reflected in the additional constant term that appears when
doing the change of variables regarding $\alpha_{i}$, which can be written as:

$b_{i}=\sum_{j}[M\cdot D_{s}]_{ij}\ \neq 0\text{ in general}$ (48)

We can choose whatever row to remove from either sides. We choose row 0:

$\displaystyle\begin{pmatrix}2&-1&0&0&0&0\\\ -1&2&0&0&0&0\\\ 0&0&2&-1&0&0\\\
0&0&-1&2&-1&0\\\ 0&0&0&-1&2&1\\\ 0&0&0&0&-1&2\\\
\end{pmatrix}\cdot\begin{pmatrix}\phi_{1}^{*}\\\ \phi_{2}^{*}\\\
\phi_{3}^{*}\\\ \phi_{4}^{*}\\\ \phi_{5}^{*}\\\ \phi_{6}^{*}\end{pmatrix}=$
$\displaystyle=\begin{pmatrix}\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}\end{pmatrix}\cdot\begin{pmatrix}\kappa_{0}\\\
\kappa_{2}\\\ \kappa_{3}\\\ \kappa_{4}\\\ \kappa_{5}\\\ \kappa_{6}\\\
\end{pmatrix}+\begin{pmatrix}\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \end{pmatrix}$

In this situation, we can solve for the set $\vec{\tilde{\kappa}}$:

$\displaystyle\begin{pmatrix}\kappa_{0}\\\ \kappa_{2}\\\ \kappa_{3}\\\
\kappa_{4}\\\ \kappa_{5}\\\ \kappa_{6}\\\
\end{pmatrix}=\begin{pmatrix}\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}&\frac{2}{7}\\\
\frac{4}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&\frac{2}{7}&-\frac{12}{7}\end{pmatrix}^{-1}\cdot$
$\displaystyle\cdot\left[\begin{pmatrix}2&-1&0&0&0&0\\\ -1&2&0&0&0&0\\\
0&0&2&-1&0&0\\\ 0&0&-1&2&-1&0\\\ 0&0&0&-1&2&1\\\ 0&0&0&0&-1&2\\\
\end{pmatrix}\cdot\begin{pmatrix}\phi_{1}^{*}\\\ \phi_{2}^{*}\\\
\phi_{3}^{*}\\\ \phi_{4}^{*}\\\ \phi_{5}^{*}\\\
\phi_{6}^{*}\end{pmatrix}-\begin{pmatrix}\frac{2}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\
\frac{2}{7}\alpha_{1}\\\ \frac{2}{7}\alpha_{1}\\\ \end{pmatrix}\right]$

Which leads to the result:

$\begin{pmatrix}\kappa_{0}\\\ \kappa_{2}\\\ \kappa_{3}\\\ \kappa_{4}\\\
\kappa_{5}\\\
\kappa_{6}\end{pmatrix}=\begin{pmatrix}\frac{-2\alpha_{1}+3\phi_{1}+\phi_{3}+3\phi_{6}}{4}\\\
\frac{3(\phi_{1}-\phi_{2})}{2}\\\
\frac{2\phi_{1}-\phi_{2}-2\phi_{3}+\phi_{4}}{2}\\\
\frac{2\phi_{1}-\phi_{2}+\phi_{3}-2\phi_{4}+\phi_{5}}{2}\\\
\frac{2\phi_{1}-\phi_{2}+\phi_{4}-2\phi_{5}-\phi_{6}}{2}\\\
\frac{2\phi_{1}\phi_{2}+\phi_{5}-2\phi_{6}}{2}\end{pmatrix}$ (49)

Note that $\kappa_{1}=0$ and $\phi_{0}=0$.

Consider the following configuration:

$\vec{\tilde{\phi}}_{(R=0)}=(0.1,0.2,0.25,-0.2,-0.1,0.0)$

In this case:

$\vec{\tilde{\kappa}}_{(C=1)}=(0.1375-\alpha_{1}/2,-0.15,-0.35,0.275,0.0,-0.05)$

Note the definition
$\kappa_{i}\equiv\alpha_{i}-\alpha_{C}\Rightarrow\alpha_{i}=\kappa_{i}+\alpha_{C}$.

If we choose $\alpha_{C}=0\Rightarrow\alpha_{i}=\kappa_{i}$, then we can
include the value of the control node $C=1$:

$\vec{\tilde{\alpha}}=(0.1375,0.0,-0.15,-0.35,0.275,0.0,-0.05)$

Alternatively, we can choose whatever value we need regarding the control
node. For instance, if $\alpha_{C}=\alpha_{1}=0.1$:

$\vec{\tilde{\alpha}}=(0.1875,0.1,-0.05,-0.25,0.375,0.1,0.05)$

and the phases configuration is the same, as shown in Figure 6

Figure 6: For the network in Fig 1, phases obtained after tuning the set of
frustration parameters to the symmetric configuration (four distinct
symmetries).

##  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.

###  B.1 Symmetric configuration case

Considering the symmetric configuration case and choosing $\alpha_{C}=0$, the
continuous optimization problem can be written as

$\displaystyle\frac{\partial e_{T}(C)}{\partial
s_{C}}=\frac{\partial}{\partial
s_{C}}|\alpha_{h}|\sum_{i}^{N-1}\Big{\lvert}1-\frac{s_{C}}{s_{i}}\Big{\rvert}=$
$\displaystyle|\alpha_{h}|\sum_{i}^{N-1}\frac{\text{sgn}(s_{C}-s_{i})}{s_{i}}$
(50)

Equation (50) is based on the function

$f(x)=\Big{\lvert}\frac{a-x}{a}\Big{\rvert}\ a,x>0$ (51)

which is depicted in Figure 7 for different values of $a$ and the sum of all
of them.

Figure 7: Three examples of the general function
$f(x)=\lvert\frac{a-x}{a}\rvert$, with $a=2$, $a=3$ and $a=4$, and the
resulting sum of them.

Regardless of the set of $a_{i}$ values, the sum function $\sum_{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.(51):

$\frac{df(x)}{dx}=\frac{\text{sgn}(x-a)}{a}$ (52)

and hence, $d\sum_{i}f(x,a_{i})/dx=\sum_{i}\text{sgn}(x-a_{i})/a_{i}$, which
is depicted in Figure 8.

Figure 8: Derivative of the function
$f(x)=\lvert\frac{2-x}{2}\rvert+\lvert\frac{3-x}{3}\rvert+\lvert\frac{4-x}{4}\rvert$
defined in Figure 7. Red dashed line at $y=0$.

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(50) behaves equivalently as the function defined in Eq.(51)
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.(50) is equivalent to the minimization of the absolute value of the
relative error:

$\sum_{i}^{N-1}\Big{\lvert}1-\frac{s_{C}}{s_{i}}\Big{\rvert}=\sum_{i}^{N}\Big{\lvert}\frac{s_{C}-s_{i}}{s_{i}}\Big{\rvert}=\sum_{i}^{N}|\mathcal{E}_{i}|$
(53)

The general problem can be written as [23]:

$\displaystyle\min_{d}\sum_{i=1}^{N}w_{i}|x_{i}-d|\ ;d>0$

with the solution:

$\displaystyle d=x_{m}\text{ where }m\equiv\min\\{i\ |\
\sum_{k=1}^{i}w_{k}\geq\sum_{k=i}^{n}w_{k}\\}$ $\displaystyle
i\in\\{1,...,n\\}$ (54)

In other words, the $d$ value that minimizes Eq.(B.1) 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|\sum_{k=1}^{i}w_{k}\geq\sum_{k=i}^{n}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.

Following the example provided in Figure 7, $\\{x\\}=\\{2,3,4\\}$ and
$\\{w\\}=\\{1/2,1/3,1/4\\}$.

The weighted median is achieved for $k=2$, corresponding to $x_{2}=3$ and
weight $w_{2}=1/3$ as $1/2+1/3=5/6>1/3+1/4=7/12$. Conversely, if we let $k=1$,
and hence $x_{1}=2$ and $w_{1}=1/2$, the condition on the weights will not be
true: $1/2\ngtr 1/2+1/3+1/4$.

## References

  * Nicolis and Nicolis [2007] G. Nicolis and C. Nicolis,  _Foundations of Complex Systems: Nonlinear Dynamics, Statistical Physics, Information and Prediction_ (World Scientific, 2007).
  * Barrat  _et al._ [2008] A. Barrat, M. Barthelemy, and A. Vespignani,  _Dynamical Processes on Complex Networks_ (Cambridge University Press, 2008).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A universal concept in nonlinear sciences_ (Cambridge University Press, Cambridge, UK, 2001).
  * Osipov  _et al._ [2007] G. Osipov, J. Kurths, and C. Zhou,  _Synchronization in Oscillatory Networks_, Springer Series in Synergetics (Springer Berlin Heidelberg, 2007).
  * Arenas  _et al._ [2008] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Synchronization in complex networks, Physics Reports 469, 93 (2008).
  * Dörfler and Bullo [2014] F. Dörfler and F. Bullo, Synchronization in complex networks of phase oscillators: A survey, Automatica 50, 1539 (2014).
  * Boccaletti  _et al._ [2018] S. Boccaletti, A. N. Pisarchik, C. I. del Genio, and A. Amann,  _Synchronization_ (Cambridge University Press, 2018).
  * [8] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in  _International Symposium on Mathematical Problems in Theoretical Physics_ (Springer-Verlag) pp. 420–422.
  * Sakaguchi and Kuramoto [1986] H. Sakaguchi and Y. Kuramoto, A Soluble Active Rotator Model Showing Phase Transitions via Mutual Entrainment, Progress of Theoretical Physics 76, 576 (1986).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Reviews of Modern Physics 77 (2005).
  * Liu and Barabási [2016] Y.-Y. Liu and A.-L. Barabási, Control principles of complex systems, Rev. Mod. Phys. 88, 035006 (2016).
  * Frolov  _et al._ [2020] N. Frolov, V. Maksimenko, S. Majhi, S. Rakshit, D. Ghosh, and A. Hramov, Chimera-like behavior in a heterogeneous kuramoto model: The interplay between attractive and repulsive coupling, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 081102 (2020).
  * Kuramoto [2003] Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Dover Publications, 2003) p. 156.
  * Nicosia  _et al._ [2013] V. Nicosia, M. Valencia, M. Chavez, A. Díaz-Guilera, and V. Latora, Remote Synchronization Reveals Network Symmetries and Functional Modules, Physical Review Letters 110, 174102 (2013).
  * Nishikawa and Motter [2016] T. Nishikawa and A. E. Motter, Symmetric states requiring system asymmetry, Physical Review Letters 117, 114101 (2016).
  * Molnar  _et al._ [2020] F. Molnar, T. Nishikawa, and A. E. Motter, Network experiment demonstrates converse symmetry breaking, Nature Physics 16, 351 (2020).
  * Zhang and Motter [2020] Y. Zhang and A. E. Motter, Symmetry-independent stability analysis of synchronization patterns, SIAM Rev. 62, 817 (2020).
  * Rosell-Tarragó and Díaz-Guilera [2020] G. Rosell-Tarragó and A. Díaz-Guilera, Functionability in complex networks: Leading nodes for the transition from structural to functional networks through remote asynchronization, Chaos 30, 062315 (2020).
  * Lei  _et al._ [2001] X. Lei, X. Li, and D. Povh, A nonlinear control for coordinating TCSC and generator excitation to enhance the transient stability of long transmission systems, Electric Power Systems Research 59, 103 (2001).
  * Acharya  _et al._ [2007] U. R. Acharya, K. P. Joseph, N. Kannathal, L. C. Min, and J. S. Suri, Heart Rate Variability, in  _Advances in Cardiac Signal Processing_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 121–165.
  * Brede and Kalloniatis [2016] M. Brede and A. C. Kalloniatis, Frustration tuning and perfect phase synchronization in the Kuramoto-Sakaguchi model, Physical Review E 93, 062315 (2016).
  * Skardal  _et al._ [2014] P. S. Skardal, D. Taylor, and J. Sun, Optimal synchronization of complex networks, Physical Review Letters 113, 144101 (2014).
  * Semsar-Kazerooni and Khorasani [2009] E. Semsar-Kazerooni and K. Khorasani, Multi-agent team cooperation: A game theory approach, Automatica 45, 2205 (2009).
  * D’Souza  _et al._ [2019] R. M. D’Souza, J. Gómez-Gardeñes, J. Nagler, and A. Arenas, Explosive phenomena in complex networks, Advances in Physics 68, 123 (2019).
  * Gómez-Gardeñes  _et al._ [2011] J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Phys. Rev. Lett. 106, 128701 (2011).