Sparse Network Optimization for Synchronization

Spontaneous synchronization happens frequently in nature: common examples are fireflies flashing, crickets chirping, planets orbiting and neurons firing. All these phenomena consist of a set of agents that exhibit a cyclic behaviour: these agents are called oscillators [1]. Two or more oscillators are said to be coupled if they influence each other by some physical or chemical process. In our setting, each oscillator is associated with the node of a graph, and the whole graph represents a network. The result of such interactions is often synchrony, which means that all entities end up with the same frequency after a short period of time. The Kuramoto model [2] was originally motivated whereby a system of coupled oscillating nodes for sufficient coupling lock on to a common frequency, in this case the mean of the ensemble of natural frequencies of the oscillators, despite the variance in those frequencies. For a recent literature review of the model see [3].

In the following we summarise salient aspects of the literature relevant to network optimization. The original model, as in [2], was posed for all-to-all coupling in a complete graph. Since then the pattern of synchronization has been explored for a variety of regular and irregular graphs. At one end, sparse regular graphs are the poorest in synchronization, by which we mean they require high coupling strength to achieve synchronization. For example, synchronization on regular structures such as trees [4] and rings [5, 6] admit a degree of solubility for critical coupling, or bounds thereon. These are in some respects balanced structures: sparse, with uniform degree distributions. However, as they scale in size they require increasingly higher values of coupling to achieve synchrony. Irregular or complex graphs, such as Erdös-Rényi [7], small world [8], and the scale-free Barabasi-Alberts [9] networks synchronize for significantly smaller coupling, with small world generally the worst of these [10]. However, their complexity leads to very uneven distribution of degrees which may incur some cost, for example imbalances in work or communication loads depending on the physical setting of the synchronization process. In our case, our interest in such graphs is for the purpose of synchronized decision-making in organizations [11] where span of control and cognitive overload require sparseness and balanced degree distribution.

In terms of approaches that explicitly use a form of optimization, [12] employed a stochastic hill-climbing approach based on random network link rewirings and acceptance or rejection on the criterion of improving the Kuramoto order parameter $r$, to be defined below. Contemporaneously, [13] applied gradient descent to the order parameter to be able to more efficiently compute optimal weighted networks. The authors of [14] follow an evolution process similar to that of [12], but, to avoid hubs in the network, the hill-climbing trajectory only accepts a rewiring that preserves the initial degree distribution. In the same work, the authors propose the use of an additional energy-like function in the selection process.

A Markov Monte Carlo replica exchange approach was developed by [15] to optimise network synchronization against noise, but the method becomes computationally demanding beyond $N=15$ size graphs. In a second order generalization of the Kuramoto model, relevant to electrical power systems, [16] used convex optimization to determine optimal placements of frequencies or link weights within an existing network structure. The examples considered in that work were of similarly small size.

A number of optimization approaches have been applied to a related, but simpler, coupled oscillator model in which the graph Laplacian figures explicitly in the interaction function, and to which the Kuramoto model approximates under linearization. Using the master stability function of [17], the authors of [18] observed that graphs with lower Laplacian spectral ratio exhibited improved synchronization. In this spirit, a stochastic rewiring approach based on improving the spectral ratio was explored by [19, 20] yielding heterogeneous and typically dense graphs. Complementary to the spectral ratio approach is that of expander graphs which have the property that every partition of the nodes into two subsets has a number of boundary edges between them that scales with the size of the smallest subset. Here [21] conducted computational search over instances of graphs with such properties, with results that have excellent expansion and synchronization properties but that increase in degree as the graphs are scaled in size.

To conclude this summary we mention that [22] have provided a construction for the Kuramoto model that yields modified trees – which are therefore sparse and well-balanced – that are provably expanders, with enhanced synchronization. They also showed that a computational search over random variations of the links on the leaf nodes of this construction (maintaining the degree distribution) to find instances with lowest possible spectral ratio did not necessarily give graphs of better synchronization.

For more details on similar techniques, which progressively search for a network so as to promote its fitness for synchronization, see [3, Section 8].

In most of these approaches, however, optimization takes the form of a hill-climbing evolution, where, according to certain rules and merit functions, a sequence of “mutations” of an initial network, is constructed. After a prescribed number of steps (usually stipulated a priori), the last mutation accepted will exhibit some “optimal” properties, in the sense that the merit functions used in the process have the best values among all the previous mutations.

In the present paper, we propose, implement and solve a new optimization model for designing graphs that (i) are sparse, (ii) can achieve synchronization in the context of the Kuramoto model [3] and (iii) have a relatively narrow range of small degrees. Conditions (i) and (iii) are helpful for ease of implementation of a network, as well as for promoting balanced sharing of communication- or work-loads.

We mathematically incorporate these three features into our optimization model. Indeed, our optimization variables are the set of all weighted (connected) undirected networks (expressed as $N\times N$ symmetric matrices with nonnegative coefficients), and a set of $N$ functions, namely the phases $\{\theta_{i}\}$ for $i=1,\ldots,N$ in the Kuramoto system. Our objective function involves a penalty term ensuring sparsity (given by the $\ell_{1}$ norm of the matrix) and two other terms promoting the small variance of the degrees (see the second and third terms in the objective function in (6)).

To the authors knowledge, there is no available approach that proposes a mathematical optimization model that performs a “universal search” among all undirected networks which are able to synchronize by means of the Kuramoto model.

The initial version of the optimization model we construct is an infinite-dimensional mixed-integer programming problem, which is in the form of an optimal control problem. By relaxing the integer variables, we obtain a continuous-variable model, which is nonsmooth and still infinite-dimensional, and therefore not tractable yet. We then discretize the problem, so that optimization software can be employed, as it is done for optimal control problems. The discretization produces a large-scale problem, which, in general, prohibits any use of nonsmooth solvers. However, in our particular model, non-smoothness appears in a peculiar form that can be transformed into a smooth, i.e., differentiable, form. We use tricks and techniques from optimization to convert the problem into a smooth one so that powerful differentiable solvers can be employed to get a solution, even though the problem is still a large-scale and non-convex one.

One should note that a solution to our discrete and smooth optimization model does not represent a graph immediately. A key procedure we propose, Algorithm 1, uses solutions of our discrete and smooth optimization model so as to construct increasingly sparse graphs while keeping the spectral ratio of the Laplacian relatively small. In other words, Algorithm 1 helps us

• (i)

  find a sparse adjacency matrix of the graph such that the system defined by the equations in (1) below is synchronized, and

• (ii)

  ensure that the resulting graph exhibits a small spectral ratio (defined in (5) below) and

• (iii)

  has a possibly small and narrow range of node degrees.

We present extensive numerical experiments by using this algorithm for graphs with $N=20$ and $N=127$. We illustrate the fact that, when adding a moderate number of edges, a great improvement in synchronization properties of the graph can be achieved. We also make comparisons with certain graphs from the literature, in particular from [22], where graphs are obtained by appending at the bottom level of a hierarchical tree, additional matchings linking the leaf nodes. The graphs found in this way in [22] enhance synchronization compared to the original tree. We emphasise the key result: that we can generate graphs with optimal sparseness, balanced load, and good synchronization properties starting from a random frequency and initial phase input. Because the graphs we generate show low degree heterogeneity, the resultant synchronization is not so sensitive to the input frequency choices.

The paper is organized as follows. In Section 2, we recall classical concepts and properties of the Kuramoto model and the spectral ratio of a graph. In Section 3 we introduce our mathematical model that promotes synchronization as well as pose its relaxation. In Section 4, we discuss discretization of the time-dependent functions in the model and smoothing. In Section 5 we present our numerical experiments and provide a detailed interpretation of the numerical results. In Section 6, we provide concluding remarks and comment on future work.

In this section we briefly recall the mathematical concepts we will use in our models, including the Tree-expander Construction approach presented in [22].

As mentioned in Section 1, our aims are (i) to find a sparse $A$, and to achieve (ii) synchronization and (iii) a narrow distribution of degrees. These aims, when realised, are expected to yield a relatively small $Q_{L}$. To obtain such a matrix $A$, we propose optimization models with decision variables $(A,\theta)$, where $A\in{\rm{I\ \kern-5.39993ptR}}^{N\times N}$ and $\theta$ is the vector function defined in (1). When the graph is assumed to be undirected, the matrix $A$ is symmetric. Similar models can be posed for the directed case. We fix a random vector $\omega$, and a time $t_{c}$ after which we impose synchronization among $\theta_{i}$’s.

As it stands in Section 3.2, the $\ell_{1}$-model is nondifferentiable. On the other hand, it is well-known that a smooth re-formulation of the $\ell_{1}$-model can be obtained by using standard nonlinear programming techniques—see e.g. [23, 24]. This allows us to use smooth optimization solvers for solving the discretized $\ell_{1}$-model.

With the elements $\widetilde{A}_{ij}$s of the matrix $\widetilde{A}$ interpreted as constant control functions and $\theta_{i}$s interpreted as the state variables, the model in (10) is an optimal control problem. We have discretized the $\ell_{1}$-model, using the Trapezoidal rule, in a fashion similar to the way optimal control problems are discretized; see, for example, [25, 26]. Challenging optimal control problems have successfully been solved using direct discretization previously [28, 30, 27, 29].

Suppose that the ODE in Model (10) is defined over the time horizon $[0,T]$, with $T=t_{c}$. Let

$$f(\theta_{1}(t),\ldots,\theta_{N}(t),\widetilde{A}):=\omega_{i}+\sum_{{j=1}}^{{N}}\widetilde{A}_{ij}\sin(\theta_{j}(t)-\theta_{i}(t)),\quad\forall t\in[0,T]\,.$$

(11)

Consider a regular partition of $[0,T]$, such that $0=t_{0}<t_{1}<\ldots<t_{M}=T$, with $t_{k+1}:=t_{k}+h$ and $h:=T/M$. Let $\theta_{i}^{k}$ be an approximation of $\theta_{i}(t_{k})$, $k=0,1,\ldots,M-1$. Then the ODE can be discretized by applying the trapezoidal rule:

$\displaystyle\theta_{i}^{k+1}$

$\displaystyle=$

$\displaystyle\theta_{i}^{k}+\frac{h}{2}\left(f(\theta_{1}^{k},\ldots,\theta_{N}^{k},\widetilde{A})+f(\theta_{1}^{k+1},\ldots,\theta_{N}^{k+1},\widetilde{A})\right),$

(12)

$k=0,1,\ldots,M-1$. After this discretization, and the replacement of the ODE with this discretization, the optimization model in (10) becomes a finite-dimensional optimization problem. With the number of optimization variables being $[(N-1)NM/2]$, the discretized problem is large-scale even with a moderate graph size $N$ and a moderate partition size $M$. Moreover, the objective function of the problem is nondifferentiable, prohibiting efficient use of smooth optimization methods and software.

Techniques from the optimization literature can be employed to transform the nonsmooth objective function into a smooth one as follows. Nonsmoothness of the first term (the $\ell_{1}$-norm) is caused by the moduli, which can be avoided by defining two new variables, $B_{ij}^{1}$ and $B_{ij}^{2}$ such that (see e.g. [24])

$$\widetilde{A}_{ij}:=B_{ij}^{1}-B_{ij}^{2}\,,\quad B_{ij}^{1},B_{ij}^{2}\geq 0\,.$$

(13)

Then one can show that

$$|\widetilde{A}_{ij}|=B_{ij}^{1}+B_{ij}^{2}\,.$$

(14)

The RHS of the ODEs in (11) then becomes

$$f(\theta_{1}(t),\ldots,\theta_{N}(t),B_{ij}^{1},B_{ij}^{2}):=\omega_{i}+\sum_{{j=1}}^{{N}}(B_{ij}^{1}-B_{ij}^{2})\,\sin(\theta_{j}(t)-\theta_{i}(t))\,,\quad\forall t\in[0,T].$$

(15)

To tackle the second and third terms in the objective function, one can define the new optimization variables (see e.g. [23])

$$\alpha:=\widetilde{d}_{\max}\quad\mbox{and}\quad\beta:=\widetilde{d}_{\min}\,.$$

Then

$$\min\frac{\widetilde{d}_{\max}}{\widetilde{d}_{\min}}\equiv\min\frac{\alpha}{\beta}$$

with

	$\displaystyle\sum_{j=1}^{N}\widetilde{A}_{ij}\leq\alpha\,,\quad\mbox{for }i=1,\ldots,N\,,$		(16)
	$\displaystyle\sum_{j=1}^{N}\widetilde{A}_{ij}\geq\beta\,,\quad\mbox{for }i=1,\ldots,N\,.$		(17)

The synchrony constraint in Model (10) is also nonsmooth, which can be replaced by a smooth version given by

$\displaystyle\sum_{i\neq j}\left(\theta_{i}^{M}-\theta_{j}^{M}\right)^{2}\leq\varepsilon\,.$

(18)

Incorporation of (11)–(18) transforms Model (10) into the smooth optimization problem below.

$$\left\{\begin{array}[]{ll}&\hbox{\em objective}:\\ &\displaystyle\min_{B^{1},B^{2},\theta^{k},\alpha,\beta}\,\,c_{1}\sum_{i,j}(B_{ij}^{1}+B_{ij}^{2})+c_{2}\,\dfrac{\alpha}{\beta}+c_{3}\,\left(\dfrac{\alpha}{\beta}-1\right)\\[17.07164pt] &\hbox{\em subject to the following constraints.}\\ &\\ &\hbox{\em symmetricity}:\ \ B_{ij}^{1}-B_{ij}^{2}=B_{ji}^{1}-B_{ji}^{2}\,,\ \ \forall i,j,\\ &\\ &\hbox{\em nonnegativity}:\ \ B_{ij}^{1}-B_{ij}^{2}\geq 0\,,\ \ \forall i,j,\\ &\\ &\hbox{\em discretization}:\\[5.69054pt] &\theta_{i}^{k+1}=\theta_{i}^{k}+\displaystyle\frac{h}{2}\left(f(\theta_{1}^{k},\ldots,\theta_{N}^{k},B_{ij}^{1},B_{ij}^{2})+f(\theta_{1}^{k+1},\ldots,\theta_{N}^{k+1},B_{ij}^{1},B_{ij}^{2})\right),\\[5.69054pt] &\hskip 199.16928ptk=0,1,\ldots,M-1\,,\ \ \theta_{i}^{0}=\theta_{i,0}\,,\ \ \forall i,\\[5.69054pt] &\hbox{\em synchrony}:\\[2.84526pt] &\displaystyle\sum_{i\neq j}\left(\theta_{i}^{M}-\theta_{j}^{M}\right)^{2}\leq\varepsilon\,,\\ &\\ &\hbox{\em auxiliary}:\ \ B_{ij}^{1},B_{ij}^{2}\geq 0\,,\forall i,j,\\ &\\ &\hskip 56.9055pt\sum_{j=1}^{N}(B_{ij}^{1}-B_{ij}^{2})\leq\alpha\,,\ \ \forall i\,,\\ &\\ &\hskip 56.9055pt\sum_{j=1}^{N}(B_{ij}^{1}-B_{ij}^{2})\geq\beta\,,\ \ \forall i\,.\end{array}\right.$$

(19)

For the above problem a standard differentiable optimization software can now be used to find a solution.

We have proposed an optimization algorithm (Algorithm 1) for designing sparse undirected networks, i.e., graphs, which synchronize well. The algorithms involve solution of a discretized relaxed mathematical model which maximize sparsity and minimize the spread of degrees (in some sense) at the same time. Synchrony of the network is posed as a constraint by means of the well-known Kuramoto system of coupled oscillators. By means of carefully chosen examples, from small to large scale, we illustrate the working of our algorithm and optimization model. The outcome is a method to generate sparse well balanced graphs with good synchronization from a random set of frequencies and initial phases. These graphs retain these properties when other frequency choices are used. In other words, in the resultant graphs there is not a strong correlation between the frequencies and degrees of the nodes (with a cubic generated in one case) underlying this insensitivity.

We have illustrated that our algorithm is successful in designing sparse synchronizing graphs with a relatively narrow distribution of degrees and small Laplacian eigenvalue ratio $Q_{L}$. By allowing more edges in the graph, but without sacrificing the sparsity much, we managed to obtain impressive synchronization behaviour in terms of the expected order parameter curve $\operatorname*{\mbox{E}}[\bar{r}]$ and its variance $\operatorname*{\mbox{Var}}[\bar{r}]$. These graphs synchronize in a robust manner in that for relatively small synchronizing values of $\bar{r}$ the variance of $\bar{r}$ under randomized data is small.

As pointed out in Remark 1, Algorithm 1 might be used as a valuable tool in identifying the so-called vital nodes and the critical set of links of the optimized sparse networks.

We have not explicitly incorporated the spectral ratio $Q_{L}$ into our optimization model. However, the numerical experiments have demonstrated that the additional terms involving the minimum and maximum degrees in the objective function promote smaller $Q_{L}$. Having said this, a measure of $Q_{L}$ might still be directly/explicitly included in a future model by using the Rayleigh quotient and the representation properties of the Laplacian using quadratic functions, such as the one in [35, Equation (14)]. Similar expressions that can be investigated can be found in [36, Remark 4.2].

Another interesting feature of a future optimization model would be constraints imposed on the degrees which the nodes of a graph should have. This would result in more targeted design of networks with certain desired load distributions.

The authors would like to acknowledge valuable discussions with Subhra Dey at the initiation of this project and with Richard Taylor during its later stages. This research was a collaboration under the auspices of the Modelling Complex Warfighting initiative between the Commonwealth of Australia (represented by the Defence Science and Technology Group) and the University of South Australia through a Defence Science Partnerships agreement.