Locating and counting equilibria of the Kuramoto model with rank one coupling

TL;DR

This paper introduces an efficient algorithm for locating real equilibria of the Kuramoto model with rank one coupling, significantly reducing computation time and providing rigorous guarantees, applicable to systems with modest oscillator counts.

Contribution

The paper presents a novel algebraic geometric algorithm that locates only real equilibria of the Kuramoto model, improving efficiency and enabling precise approximation, with theoretical and practical validation.

Findings

Algorithm finds all and only real equilibria.

Computational time is reduced by several orders of magnitude.

Maximum number of equilibria grows at the same rate as complex solutions.

Abstract

The Kuramoto model describes synchronization behavior among coupled oscillators and enjoys successful application in a wide variety of fields. Many of these applications seek phase-coherent solutions, i.e., equilibria of the model. Historically, research has focused on situations where the number of oscillators, $n$, is extremely large and can be treated as being infinite. More recently, however, applications have arisen in areas such as electrical engineering with more modest values of $n$. For these, the equilibria can be located by finding the real solutions of a system of polynomial equations utilizing techniques from algebraic geometry. However, typical methods for solving such systems locate all complex solutions even though only the real solutions give equilibria. In this paper, we present an algorithm to locate only the real solutions of the model, thereby shortening computation time by several orders of magnitude in certain situations. This is accomplished by choosing specific equilibria representatives and the consequent algebraic decoupling of the system. The correctness of the algorithm (that it finds only and all the equilibria) is proved rigorously. Additionally, the algorithm can be implemented using interval methods so that the equilibria can be approximated up to any given precision without significantly more computational effort. We also compare this solving approach to other computational algebraic geometric methods. Furthermore, analyzing this approach allows us to prove, asymptotically, that the maximum number of equilibria grows at the same rate as the number of complex solutions of a corresponding polynomial system. Finally, we conjecture an upper bound on the maximum number of equilibria for any number of oscillators which generalizes the known cases and is obtained on a range of explicitly provided natural frequencies.

Tables

Table 1. Table 1: Comparison of various solving methods

$n$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
# real	$2$	$2$	$4$	$4$	$4$	$4$	$4$	$4$	$8$	$8$
# complex	$6$	$12$	$28$	$56$	$118$	$238$	$486$	$976$	$1972$	$3958$
Macaulay2 degree	s	s	$0.1$s	$1.1$s	$7.0$s	$72.6$s	$716.5$s	$10783.7$s	$149578.0$s	$‡$
Bertini regeneration	$0.3$s	$1.2$s	$3.4$s	$13.4$s	$45.1$s	$116.6$s	$210.1$s	$486.2$s	$1493.1$s	$3443.5$s
Bertini parameter	s	s	$0.2$s	$0.4$s	$1.1$s	$2.2$s	$6.9$s	$15.0$s	$36.9$s	$116.8$s
Algorithm 3.15	s	s	s	s	s	s	s	s	s	s

Table 2. Table 2: Generic root counts for symmetric coupling matrices of various ranks

$n$	rank $1$	rank $2$	rank $3$	rank $4$	rank $5$
2	2	2	2	2	2
3	6	6	6	6	6
4	14	20	20	20	20
5	30	70	70	70	70
6	62	232	252	252	252
7	126	714	924	924	924
8	254	2056	3362	3432	3432
9	510	5646	11,860	12,870	12,870
10	1022	14,864	40,136	48,368	48,620

Equations

$$\frac{d{\theta }_{\nu }}{dt}={\omega }_{\nu }-\frac{K}{n}\sum _{\mu =1}^n\sin ({\theta }_{\nu }-{\theta }_{\mu }),for\nu =1,\dots ,n$$

$$\frac{d{\theta }_{\nu }}{dt}={\omega }_{\nu }-\frac{1}{n}\sum _{\mu =1}^n{k}_{\nu }{k}_{\mu }\sin ({\theta }_{\nu }-{\theta }_{\mu }),for\nu =1,\dots ,n.$$

$${\omega }_{\nu }=\frac{1}{n}\sum _{\mu =1}^n{k}_{\nu }{k}_{\mu }\sin ({\theta }_{\nu }-{\theta }_{\mu }),for\nu =1,\dots ,n.$$

$${\omega }_{\nu }=\frac{1}{n}\sum _{\mu =1}^n{k}_{\nu }{k}_{\mu }\left(s_{\nu }{c}_{\mu }-s_{\mu }{c}_{\nu }\right),$$

$${\omega }_{\nu }=\frac{1}{n}\sum _{\mu =1}^n{\kappa }_{\nu \mu }\sin ({\theta }_{\nu }-{\theta }_{\mu }),for\nu =1,\dots ,n,$$

$$\sum _{\nu =1}^n\frac{d{\theta }_{\nu }}{dt}=\sum _{\nu =1}^n{\omega }_{\nu }$$

$$\begin{array}[]{ccccc}\sin(\theta_{1}-\theta_{2})+\sin(\theta_{1}-\theta_{3})+\sin(\theta_{1}-\theta_{4})&=&\sin(\theta_{2}-\theta_{1})+\sin(\theta_{2}-\theta_{3})+\sin(\theta_{2}-\theta_{4})&=&\\ \sin(\theta_{3}-\theta_{1})+\sin(\theta_{3}-\theta_{2})+\sin(\theta_{3}-\theta_{4})&=&\sin(\theta_{4}-\theta_{1})+\sin(\theta_{4}-\theta_{2})+\sin(\theta_{4}-\theta_{3})&=&0.\end{array}$$

$${\theta }_3={\theta }_1+\pi and{\theta }_4={\theta }_2+\pi .$$

$${\Theta }_{\omega ,k}$$

$${\Theta }_{\omega ,k,\sigma }$$

$${\mathcal{R}}_{\omega ,k,\sigma }$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=\frac{1}{n}\sum _{\mu =1}^n{k}_{\nu }{k}_{\mu }\sin \left({\theta }_{\nu }-{\theta }_{\mu }\right)\text{and }OC1\right\}.$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=\frac{1}{n}{k}_{\nu }\mathrm{Im}{e}^{i{\theta }_{\nu }}\sum _{\mu =1}^n{k}_{\mu }{e}^{-i{\theta }_{\mu }}\text{and }OC1\right\}.$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=\frac{1}{n}{k}_{\nu }\mathrm{Im}{e}^{i{\theta }_{\nu }}\overline{\sum _{\mu =1}^n{k}_{\mu }{e}^{i{\theta }_{\mu }}}\text{and }OC1\right\}.$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{r\in {\mathbb{R}}_{\ge 0}}{\exists }r=\frac{1}{n}\sum _{\mu =1}^n{k}_{\mu }{e}^{i{\theta }_{\mu }}\text{and }\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=r{k}_{\nu }\mathrm{Im}{e}^{i{\theta }_{\nu }}\right\}.$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{r\in {\mathbb{R}}_{>0}}{\exists }r=\frac{1}{n}\sum _{\mu =1}^n{k}_{\mu }{e}^{i{\theta }_{\mu }}\text{and }\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=r{k}_{\nu }\sin {\theta }_{\nu }\right\}.$$

$${\Theta }_{\omega ,k}=\left\{\theta \in U:\underset{r\in {\mathbb{R}}_{>0}}{\exists }r=\frac{1}{n}\sum _{\mu =1}^n{k}_{\mu }\cos {\theta }_{\mu }\text{and }\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }{\omega }_{\nu }=r{k}_{\nu }\sin {\theta }_{\nu }\right\}.$$

$${\Theta }_{\omega ,k}$$

$${\Theta }_{\omega ,k,\sigma }$$

$${\Theta }_{\omega ,k,\sigma }=\left\{\theta \in U:\underset{r\in {\mathbb{R}}_{>0}}{\exists }{r}^2=\frac{1}{n}\sum _{\mu =1}^n{\sigma }_{\mu }\sqrt{k_{\mu }^2{r}^2-{\omega }_{\mu }^2}\text{and }\underset{\nu \in \left\{1,\dots ,n\right\}}{\forall }\sin {\theta }_{\nu }=\frac{{\omega }_{\nu }}{k_{\nu }r}\text{and }\mathrm{sign}\cos {\theta }_{\nu }={\sigma }_{\nu }\right\}.$$

$$\sum _{\mu =1}^n{k}_{\mu }{e}^{i{\theta }_{\mu }}\in {\mathbb{R}}_{>0}$$

$$\sum _{\mu =1}^n{k}_{\mu }{e}^{i({\psi }_{\mu }+\varphi )}=e^{i\varphi }\sum _{\mu =1}^n{k}_{\mu }{e}^{i{\theta }_{\mu }}\in {\mathbb{R}}_{>0}.$$

$$f_{\sigma }(R)=-R+\frac{1}{n}\sum _{\mu =1}^n{\sigma }_{\mu }\sqrt{k_{\mu }^{2}R-{\omega }_{\mu }^2}.$$

$$\sin {\theta }_{\nu }=\frac{{\omega }_{\nu }}{k_{\nu }\sqrt{R}}and\mathrm{sign}\cos {\theta }_{\nu }={\sigma }_{\nu }for\nu =1,\dots ,n.$$

$$\left[{\left(\frac{{\omega }_n}{k_n}\right)}^2,{\left(\frac{1}{n}\sum _{\mu \in {\sigma }_{+}}{k}_{\mu }\right)}^2\right].$$

$$R\ge {\left(\frac{{\omega }_{\nu }}{k_{\nu }}\right)}^{2}for\nu =1,\dots ,n.$$

$$0\le \sqrt{k_{\mu }^{2}R-{\omega }_{\mu }^2}\le k_{\mu }\sqrt{R}.$$

$$R$$

$$=\left(\frac{1}{n}\sum _{\mu \in {\sigma }_{+}}\sqrt{k_{\mu }^{2}R-{\omega }_{\mu }^2}\right)-\left(\frac{1}{n}\sum _{\mu \in {\sigma }_{-}}\sqrt{k_{\mu }^{2}R-{\omega }_{\mu }^2}\right)$$

$$\le \frac{1}{n}\sum _{\mu \in {\sigma }_{+}}{k}_{\mu }\sqrt{R}.$$

$$s_{\ell }=\sum _{\mu =1}^{\ell }{\sigma }_{\mu }{k}_{\mu }\le 0,$$

$$R\ge {\left(\frac{{\omega }_n}{k_n}\right)}^2\ge {\left(\frac{{\omega }_{n-1}}{k_{n-1}}\right)}^2\ge \cdots \ge {\left(\frac{{\omega }_1}{k_1}\right)}^2$$

$$\sqrt{R-{\left(\frac{{\omega }_1}{k_1}\right)}^2}\ge \cdots \ge \sqrt{R-{\left(\frac{{\omega }_n}{k_n}\right)}^2}\ge 0.$$

Full text

In this paper, we present an algorithm to locate only the real solutions of the model, thereby shortening computation time by several orders of magnitude in certain situations. This is accomplished by choosing specific equilibria representatives and the consequent algebraic decoupling of the system. The correctness of the algorithm (that it finds only and all the equilibria) is proved rigorously. Additionally, the algorithm can be implemented using interval methods so that the equilibria can be approximated up to any given precision without significantly more computational effort. We also compare this solving approach to other computational algebraic geometric methods.

Furthermore, analyzing this approach allows us to prove, asymptotically, that the maximum number of equilibria grows at the same rate as the number of complex solutions of a corresponding polynomial system. Finally, we conjecture an upper bound on the maximum number of equilibria for any number of oscillators which generalizes the known cases and is obtained on a range of explicitly provided natural frequencies.

Keywords. Kuramoto model, equilibria, univariate solving, homotopy continuation, numerical algebraic geometry

AMS Subject Classification. 65H10, 68W30, 14Q99

1 Introduction

Oscillatory dynamics characterize many important systems. For such systems, it is important to understand the synchronization behavior of coupled oscillators, especially when conducting stability assessments. Synchronization behavior is characterized by the equilibria of the associated dynamic model. This paper is concerned with locating and counting equilibria of a certain generalization of the Kuramoto model [22], which we call a rank-one coupled Kuramoto model.

Kuramoto model:

The standard Kuramoto model for $n\geq 2$ oscillators has all-to-all and uniform coupling among the oscillators. It is formulated as the following system of coupled first-order ordinary differential equations:

[TABLE]

where $K>0$ is the uniform coupling strength, and each parameter $\omega_{\nu}$ and variable $\theta_{\nu}$ denote the natural frequency and phase angle of the $\nu^{\rm th}$ oscillator, respectively. There is a large body of literature for the Kuramoto model (1) and its many variants, e.g., non-uniform coupling among oscillators and allowance for second-order dynamics. The wide variety of applications of the Kuramoto model in modeling oscillatory behavior include electrical engineering [14, 15, 39], biology [36], and chemistry [3, 23, 33]. See [1, 15, 37] and the references therein for a more detailed survey of the relevant literature and extensive applications.

Rank-one coupled Kuramoto model:

In this paper, we consider a slight generalization with a non-uniform coupling between the oscillators described by a symmetric rank-one matrix. In particular, for $k=(k_{1},\dots,k_{n})\in\mathbb{R}_{>0}^{n}$, the $\nu^{\rm th}$ and $\mu^{\rm th}$ oscillators are coupled with strength $k_{\nu}k_{\mu}$ yielding the model

[TABLE]

The standard Kuramoto case (1) corresponds with $k=(\sqrt{K},\dots,\sqrt{K})$. We are concerned with the equilibria of the rank-one coupled Kuramoto model (2), which are the real solutions to the system of nonlinear equations resulting from setting $\dfrac{d\theta_{\nu}}{dt}$ equal to [math] in (2), namely

[TABLE]

That is, we aim to compute the values of the variables $\theta_{1},\dots,\theta_{n}$ such that (3) holds for given values of the parameters $n$, $k_{1},\dots,k_{n}$, and $\omega_{1},\dots,\omega_{n}$. This generalization was originally motivated by applications where the coupling is non-uniform, such as in a power flow model [14, 15, 39] in which the coupling matrix could be of arbitrary rank. However, as demonstrated in Ex. 3.11, with a lossless power system and uniform line susceptances, the equilibria of the power flow equations correspond to the equilibria of the rank-one coupled Kuramoto model (3). Hence, (3) can be viewed as an initial generalization (rank one) toward the full generalization (arbitrary rank).

We address two natural problems: locating all equilibria and counting them.

Locating all equilibria:

In [28, 29, 34], homotopy continuation and numerical algebraic geometry [4, 35] were applied to the standard Kuramoto model and various non-uniform coupling generalizations by converting the corresponding system describing the equilibria into a polynomial system. For example, with $s_{\nu}=\sin(\theta_{\nu})$ and $c_{\nu}=\cos(\theta_{\nu})$, (3) corresponds to the polynomial system

[TABLE]

Even though all complex solutions were computed, only the real solutions are physically meaningful, i.e., correspond to equilibria, so that a post-processing step is necessary to filter out the non-real solutions. In other words, homotopy continuation expends computational effort to compute all complex solutions when only the real solutions are relevant. Using parallel computing techniques, such a method has been applied to problems with $n\leq 18$ [29].

In [25], a specialized continuation method was proposed which computes only equilibria so that the computational cost scales with the number of real solutions of the corresponding polynomial system rather than the number of complex solutions. Moreover, this continuation method is applicable to a more general class of problems (the power flow equations) which include (3) as a special case. However, the robustness proof showing that it locates all equilibria for this more general class of problems was shown to be flawed [9] with a counterexample presented in [30]. In [24], a modification of the method based on an elliptical reformulation of equations was shown to have improved robustness. There currently does not exist a robustness proof for this modification or a known counterexample, so the capabilities of this method remain to be fully characterized.

In summary, despite significant progress, the aforementioned approaches either quickly become intractable as $n$ increases or are not proven to find all equilibria. One of the main contributions of this paper is to provide a new algorithm that can handle much larger values of $n$ which is also rigorously proved to find all equilibria. For instance, in Section 3.3, we demonstrate our approach on an example with $n=60$ which computes all equilibria in under a second.

Counting equilibria:

The second problem is to determine the maximum number of equilibria (up to trivial shifts – see Section 2). Existing upper bounds on the number of equilibria are based on bounds for the number of complex solutions to (4). In [2], an upper bound on the number of equilibria of the Kuramoto model with an arbitrary coupling matrix $\kappa\in\mathbb{R}^{n\times n}$, i.e., the equilibria satisfy

[TABLE]

is $\binom{2n-2}{n-1}$. This bound is sharp for $n=2$ and $n=3$. It is an open question (first posed in [2, Question 5.1]) whether the upper bound of $\binom{2n-2}{n-1}$ can be achieved for $n\geq 4$.

Other research [10, 11, 31] has produced tighter upper bounds on the number of complex solutions to (5) when the oscillators are not completely connected, i.e., “topologically dependent” bounds.

The number of equilibria for the standard Kuramoto model has been studied for small values of $n$. In the standard Kuramoto setting, i.e., $k=(\sqrt{K},\dots,\sqrt{K})$, there are at most $2$ equilibria satisfying (3) when $n=2$. For $n=3$ and $n=4$, elimination theory was used in [40] to produce a degree six and degree fourteen univariate polynomial, respectively, yielding bounds of at most $6$ and $14$ equilibria. Morse Theory was used to derive similar results for the $n=3$ and $n=4$ cases in [2]. These aforementioned bounds are tight for $n=2$ and $n=3$, but it is currently unknown whether the upper bound of $14$ can be achieved for $n=4$. The authors of [40] find a maximum of $10$ equilibria in the $n=4$ case, which is smaller than the upper bound of $14$. Since this maximum was obtained via a computational experiment which gridded the parameter space, they conjecture that $10$ is indeed the maximum number of equilibria when $n=4$.

In summary, despite significant progress, there remains several open questions regarding the number of equilibria to the rank-one coupled Kuramoto model. First, for the polynomial system (4), the generic root count, which is the number of solutions for generic values of the parameters, is unknown. Clearly, this is bounded above by $\binom{2n-2}{n-1}$ which is the generic root count for the corresponding polynomial system in the arbitrary coupling case (5). Moreover, the quality of the relationship between the maximum number of equilibria and the generic root count has not been explored. Three contributions of this paper are to provide such a generic root count for (4), count the number of equilibria in particular cases, and use these cases to analyze the asymptotic behavior of the ratio between the maximum number of equilibria and the generic root count for (4).

Approach:

This paper locates and counts equilibria for arbitrary $n$ by reformulating (3) into a family of decoupled univariate radical equations. (This reformulation is similar in spirit but different than the approach in [40]. Further, the proposed reformulation is not limited to $n=2,3,4$.) This reformulation enables the development of both new theoretical results and computational tools. Our solving algorithm exploits this reformulation together with new results regarding cases where equilibria cannot exist. Computational experiments demonstrate that this algorithm can be several orders of magnitude faster than the more general computational algebraic geometric methods [4, 17, 28, 29, 34, 35] and elliptical continuation [24, 25] algorithms when applied to (4).

This reformulation allows us to count the number of equilibria for (3) where the parameters are carefully chosen to have many equilibria. These results extend a conjecture from [40] that the maximum number of equilibria for the standard Kuramoto model when $n=4$ is $10$. Moreover, the particular cases allow us to show that the maximum number of equilibria and the generic root count for (4) have the same asymptotic scaling. This suggests that algorithms which only compute equilibria to (3) will, in the worst-case, computationally scale similar to algorithms that compute all the complex solutions to (4). However, it may be the case that algorithms which compute only the real solutions to (4) have significant computational advantages for many practical problems. For instance, typical operating conditions of power flow problems are expected to have few equilibria relative to the number of complex solutions [34].

The rest of the paper is organized as follows. Section 2 presents the decoupling approach. Section 3 describes an algorithm that uses this reformulation to compute all equilibria satisfying (3) and compares the computational performance of this algorithm with other methods. Section 4 counts the number of equilibria in particular cases and compares the maximum number of equilibria with the generic root count of (4). A short conclusion is provided in Section 5.

2 Decoupled reformulation

One approach to solving a multivariate system of equations involves first decoupling the system. In this section, we take this approach and decouple the system (3). A standard method for decoupling is to apply computational tools from elimination theory (i.e., multivariate resultants and Gröbner basis techniques [6, 7, 8, 12, 13, 16, 26, 38]) to the polynomial system (4), thereby obtaining several univariate polynomials, say $f$ and $g_{1},\ldots,g_{2n}$, such that each solution of (4) is the value of the polynomial system $g=(g_{1},\dots,g_{2n})$ evaluated at a root of $f$. However, the major drawback of this approach is that the obtained polynomials $f$ and $g_{1},\dots,g_{2n}$ are of very high degree (exponential in $n$) with no naturally discernible structure. Thus solving with this method is very time-consuming, even for relatively small $n$. In the following, we will instead use an alternate method to decouple the system (3) adapted from Kuramoto’s approach [22, 23, §5.4]. This method allows us to exploit the inherent structure of the equations to obtain explicit radical expressions that are quickly solvable by standard univariate solvers.

Using (2), it is easy to see that

[TABLE]

so that equilibria can only exist when $\sum_{\nu=1}^{n}\omega_{\nu}=0$. Therefore, we only need to consider solving such cases, which we list as our first input condition (IC):

1. IC1:

$\omega_{1}+\cdots+\omega_{n}=0$.

If $\theta=(\theta_{1},\dots,\theta_{n})$ is a solution of (3), then shifting all angles by $\phi$, i.e., $(\theta_{1}+\phi,\dots,\theta_{n}+\phi)$, is also a solution. Thus, we want to both compute and count equilibria modulo shift. One approach, e.g., as used in [29], is to set one of the angles, say $\theta_{n}$, to be zero. A second approach is to fix the “weighted average angle” which is the following output condition (OC):

1. OC1:

$\sum_{\mu=1}^{n}k_{\mu}e^{i\theta_{\mu}}\in\mathbb{R}_{\geq 0}$

In paritcular, OC1 is equivalent to selecting $\theta$ so that $\sum_{\mu=1}^{n}k_{\mu}\sin(\theta_{\mu})=0$ and $\sum_{\mu=1}^{n}k_{\mu}\cos(\theta_{\mu})\geq 0$. It is a natural extension of the condition used by Kuramoto [22, 23, §5.4].

When each $\omega_{i}=0$, the following shows that there can be infinitely many equilibria modulo shift.

Example 2.1

For $n=4$ with $\omega_{i}=0$ and $k_{i}=1$, (3) is equivalent to

[TABLE]

This system has infinitely many equilibria modulo shift which can be seen, for example, by taking

[TABLE]

Since we aim to enumerate all equilibria modulo shift, we will not consider the case when every $\omega_{i}=0$ and will leave this positive-dimensional case as a possible future research direction. This forms our second input condition:

1. IC2:

$\omega\neq\left(0,\ldots,0\right)$

Finally, we want to consider rank-one coupled Kuramoto models that are fully coupled, i.e., $k_{\nu}>0$ for $\nu=1,\dots,n$, so that each oscillator is positively impacted by every other oscillator. This is a natural extension of the classical Kuramoto model (1) for which the uniform coupling strength $K$ is positive. With this assumption, we can, without loss of generality, adjust the indexing to order the input parameters based on $|\omega_{\nu}/k_{\nu}|$. This forms our third input condition:

1. IC3:

$k_{1},\ldots,k_{n}>0$ which are ordered so that $\left|\dfrac{\omega_{1}}{k_{1}}\right|\leq\left|\dfrac{\omega_{2}}{k_{2}}\right|\leq\cdots\leq\left|\dfrac{\omega_{n}}{k_{n}}\right|$

With this setup, we are now ready to decouple the multivariate system of equations (3).

Theorem 2.2** (Decoupled Reformulation)**

Suppose that $\omega\in\mathbb{R}^{n}$ and $k\in\mathbb{R}_{>0}^{n}$ satisfy IC1, IC2, and IC3. If $\Theta_{\omega,k}$ is the set of all equilibria described via OC1 satisfying (3), then

[TABLE]

Proof: For $U=(-\pi,\pi]^{n}$, we have

[TABLE]

Since $\sin\eta=\operatorname{Im}e^{i\eta}$, IC3 and factoring yields

[TABLE]

From $e^{-i\alpha}=\overline{e^{i\alpha}}$ and IC3, we have

[TABLE]

Since OC1 is equivalent to $\underset{r\in\mathbb{R}_{\geq 0}}{\exists}\ r=\frac{1}{n}\sum_{\mu=1}^{n}k_{\mu}e^{i\theta_{\mu}}$, we have

[TABLE]

If $r=0,$ then $\omega=\left(0,\ldots,0\right),$ contradicting IC2. Hence, $r\neq 0$ so that

[TABLE]

From IC1, we have $0=\sum_{\mu=1}^{n}\omega_{\mu}=\sum_{\mu=1}^{n}rk_{\mu}\sin\theta_{\mu}=r\sum_{\mu=1}^{n}k_{\mu}\sin\theta_{\mu}.$ Since $r\neq 0$, we know $\sum_{\mu=1}^{n}k_{\mu}\sin\theta_{\mu}=0$. Thus,

[TABLE]

Since $\cos\alpha=\pm\sqrt{1-\sin^{2}\alpha}$ and $\sin\theta_{\mu}=\dfrac{\omega_{\mu}}{k_{\mu}r}$, we have

[TABLE]

Since $r,k_{\mu}>0$, we can simplify to

[TABLE]

Since $r>0$, for $R=r^{2}>0$, we have $r=\sqrt{R}$ yielding the result. $\Box$

Remark 2.3

The proof of Theorem 2.2 shows that we could update OC1 to be

[TABLE]

which yields that a unique representative is computed for each equilibria modulo shift. In particular, if $\psi=(\psi_{1},\dots,\psi_{n})$ is such that $\sum_{\mu=1}^{n}k_{\mu}e^{i\psi_{\mu}}\neq 0$, then there is a unique $\phi\in(-\pi,\pi]$ such that

[TABLE]

3 Locating the equilibria

Theorem 2.2 immediately yields an algorithm for locating all equilibria satisfying (3). After some improvements, we compare the resulting method with other approaches.

3.1 Basic algorithm

For each $\sigma\in\{-1,+1\}^{n}$, the first step to utilize Theorem 2.2 for locating all equilibria is to find the positive roots of

[TABLE]

The following algorithm depends upon a root finding method that returns the set of all positive roots of $f_{\sigma}$, denoted $\mathbf{Solve}(f_{\sigma},+)$. Our implementation uses an interval Newton method [18, Chap. 6], which allows each positive root of $f_{\sigma}$ to be approximated up to any given precision. For each positive root $R$ of $f_{\sigma}$, the second step from Theorem 2.2 is to compute the equilibria via

[TABLE]

This is summarized in the following algorithm.

Algorithm 3.1** (Basic)**

In:

$\omega\in\mathbb{R}^{n}$* and $k\in\mathbb{R}_{>0}^{n}$ satisfying IC1, IC2, and IC3.*

Out:

$\Theta$, the set of equilibria satisfying OC1.

$\Theta\leftarrow\{\}$** 2. 2.

For $\sigma\in{\left\{-1,+1\right\}^{n}}$ do

1. (a)

$f_{\sigma}\leftarrow-R+\frac{1}{n}\sum_{\mu=1}^{n}\sigma_{\mu}\sqrt{k_{\mu}^{2}R-\omega_{\mu}^{2}}$** 2. (b)

$\mathcal{R}\leftarrow\mathbf{Solve}(f_{\sigma},+)$** 3. (c)

For $R\in\mathcal{R}$ do

1. i.

Compute $\theta\in(-\pi,\pi]^{n}$ such that $\sin\theta_{\nu}=\frac{\omega_{\nu}}{k_{\nu}\sqrt{R}}\hbox{~{}and~{}}\mathrm{sign}\cos\theta_{\nu}=\sigma_{\nu}$ for $\nu=1,\dots,n$. 2. ii.

Add $\theta$ to $\Theta$

Example 3.2

To illustrate for $n=2$, consider $\omega=(4,-4)$ and $k=(5,2)$. There are $4$ sign patterns $\sigma$ to consider:

• •

$\sigma=(-1,-1):$**

• $\circ$

$f_{\sigma}(R)=-R+\frac{1}{2}\left(-\sqrt{25R-16}-\sqrt{4R-16}\right)$* has no positive roots.*

• •

$\sigma=(-1,+1):$**

• $\circ$

$f_{\sigma}(R)=-R+\frac{1}{2}\left(-\sqrt{25R-16}+\sqrt{4R-16}\right)$* has no positive roots.*

• •

$\sigma=(+1,-1):$**

• $\circ$

$f_{\sigma}(R)=-R+\frac{1}{2}\left(\sqrt{25R-16}-\sqrt{4R-16}\right)$* has one positive root, namely $R=4.25$.*

• $\circ$

This yields the equilibrium $\theta=(0.3985,-1.8158)$.

• •

$\sigma=(+1,+1):$**

• $\circ$

$f_{\sigma}(R)=-R+\frac{1}{2}\left(\sqrt{25R-16}+\sqrt{4R-16}\right)$* has one positive root, namely $R=10.25$.*

• $\circ$

This yields the equilibrium $\theta=(0.2526,-0.6747)$.

In summary, there are two equilibria satisfying (3).

3.2 Optimizations

In Algorithm 3.1, $\mathbf{Solve}(f_{\sigma},+)$, which computed all positive roots of $f_{\sigma}$, was called for all $2^{n}$ sign patterns. This exponential scaling in the number of oscillators is not much better than the previous approaches discussed earlier. As such, the goal of this section is to prune out sign patterns $\sigma$ for which $f_{\sigma}$ as in (6) has no positive roots. This improvement allows for the optimized algorithm to essentially scale based on the number of equilibria, provided an extra condition is satisfied, yielding much shorter computation times.

Throughout this section, we assume $\omega\in\mathbb{R}^{n}$ and $k\in\mathbb{R}_{>0}^{n}$ satisfy IC1 – IC3, $\sigma\in\{-1,+1\}^{n}$, and $f_{\sigma}$ as in (6). With this setup, the following provides an interval containing all positive roots of $f_{\sigma}$.

Proposition 3.3

If $\sigma_{+}=\{\mu:\sigma_{\mu}=+1\}$, then every positive root of $f_{\sigma}$ is contained in the interval

[TABLE]

Proof: Suppose that $R$ is a positive root of $f_{\sigma}$. Since each $\frac{\omega_{\nu}^{2}}{k_{\nu}^{2}R}=\sin^{2}\theta_{\nu}\leq 1$ by Theorem 2.2,

[TABLE]

Hence, IC3 shows that $R\geq\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2}$.

Moreover,

[TABLE]

For $\sigma_{-}=\{\mu:\sigma_{\mu}=-1\}$, we have

[TABLE]

This is equivalent to $R\leq\left(\frac{1}{n}\sum_{\mu\in\sigma_{+}}k_{\mu}\right)^{2}$. $\Box$

Example 3.4

With the setup from Ex. 3.2, we consider the four cases:

• •

$\sigma=(-1,-1):$**

• $\circ$

no positive roots since Prop. 3.3 provides the “interval” $[4,0]$.

• •

$\sigma=(-1,+1):$**

• $\circ$

no positive roots since Prop. 3.3 provides the “interval” $[4,1]$.

• •

$\sigma=(+1,-1):$**

• $\circ$

Prop. 3.3 provides the interval $[4,6.25]$, which contains the positive root $R=4.25$.

• •

$\sigma=(+1,+1):$**

• $\circ$

Prop. 3.3 provides the interval $[4,12.25]$, which contains the positive root $R=10.25$.

As shown in Ex. 3.4, Prop. 3.3 can exclude sign patterns $\sigma$ for which $f_{\sigma}$ has no positive roots. The following provides another such test.

Proposition 3.5

If, for all $\ell=1,2,\ldots,n$,

[TABLE]

then $f_{\sigma}$ has no positive roots.

Proof: Suppose that $R$ is a positive root of $f_{\sigma}$. Then, by Prop. 3.3 and IC3,

[TABLE]

so that

[TABLE]

Thus, for every $\ell=2,\dots,n$, we have $s_{\ell}=s_{\ell-1}+\sigma_{\ell}k_{\ell}\leq 0$ by definition and

[TABLE]

Hence, combining with $f_{\sigma}(R)=0$, we have

[TABLE]

This is a contradiction since $R>0$. $\Box$

Example 3.6

With the setup from Ex. 3.2, Prop. 3.5 shows that $f_{\sigma}$ can have no positive roots for $\sigma=(-1,-1)$ and $\sigma=(-1,+1)$.

The following will be used to show additional conditions for which $f_{\sigma}$ has no positive roots.

Lemma 3.7

$f_{\sigma}$* has no positive roots if and only if $f_{\sigma}<0$ on $\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\infty\right)$.*

Proof: Let $I=\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\infty\right)$. If $R\in I$, then

[TABLE]

so that

[TABLE]

Since $f_{\sigma}$ is continuous on $I$, we must have $f_{\sigma}<0$ on $I$ when $f_{\sigma}$ has no positive roots. Furthermore, the interval $\left[\left(\frac{\omega_{n}}{k_{n}}\right)^{2},\;\;\left(\frac{1}{n}\sum_{\mu\in\sigma_{+}}k_{\mu}\right)^{2}~{}\right]$ contains all the positive roots of $f_{\sigma}$ by Prop. 3.3 and is contained in $I$. Hence, if $f_{\sigma}<0$ on $I$, then $f_{\sigma}$ has no positive roots. $\Box$

If $f_{\sigma}$ has no positive roots, the following yields additional cases which also have no positive roots.

Lemma 3.8

Let $\mu$ be such that $\sigma_{\mu}=+1$. Let $\sigma^{\prime}$ such that $\sigma_{\mu}^{\prime}=-1$ and $\sigma^{\prime}_{\nu}=\sigma_{\nu}$ for $\nu\neq\mu$. If $f_{\sigma}$ has no positive roots, then $f_{\sigma^{\prime}}$ also has no positive roots.

Proof: Since $f_{\sigma}$ has no positive roots, Lemma 3.7 shows that $f_{\sigma}<0$ on $\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\infty\right)$. Since $f_{\sigma^{\prime}}\leq f_{\sigma}$ on $\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\infty\right)$, $f_{\sigma^{\prime}}$ also does not have any positive roots by Lemma 3.7. $\Box$

Example 3.9

With the setup from Ex. 3.2, since $f_{\sigma}$ for $\sigma=(-1,+1)$ has no positive roots, $f_{\sigma^{\prime}}$ also has no positive roots for $\sigma^{\prime}=(-1,-1)$.

The following excludes additional cases by swapping entries of $\sigma$.

Lemma 3.10

Suppose $\mu$ and $\nu$ are such that $\sigma_{\mu}=+1$ and $\sigma_{\nu}=-1$. Let $\sigma^{\prime}$ be the same as $\sigma$ except that $\sigma_{\mu}^{\prime}=\sigma_{\nu}=-1$ and $\sigma_{\nu}^{\prime}=\sigma_{\mu}=+1$. If $f_{\sigma}$ has no positive roots where

[TABLE]

then $f_{\sigma^{\prime}}$ also has no positive roots.

Proof: Given (8) and $R\in\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\left(\displaystyle\frac{1}{n}\sum_{\mu=1}^{n}k_{\mu}\right)^{2}~{}\right]$, we have

[TABLE]

Rearranging gives

[TABLE]

Hence, $f_{\sigma^{\prime}}(R)\leq f_{\sigma}(R)$. Therefore, the result follows from Lemma 3.7. $\Box$

A natural way to order the sign patterns $\sigma$ is to present them using binary representations of the numbers in base $10$ from [math] to $2^{n}-1$ where “0” in binary represents $-1$ and “1” in binary represents $+1$. For example, $\sigma=(+1,-1)$ corresponds to the binary number $10_{2}$, so we can say $\sigma$ corresponds to the number $2$ in base $10$. We demonstrate this on a concrete application (power flow analysis) from electrical engineering [15] and apply all the previous results.

Example 3.11** (Power flow model: 4-bus system)**

Figure 1 depicts a lossless four-bus power system with active power injections $P_{1},\ldots,P_{4}$, voltage magnitudes $\left|V_{1}\right|,\ldots,\left|V_{4}\right|$, and line susceptances $b_{12}=b_{13}=b_{14}=b_{23}=b_{24}=b_{34}=-1$. The equilibria of the power flow equations correspond to the equilibria of the rank-one coupled Kuramoto model, namely the solutions of (3) where $\omega=\left(P_{1},\ldots,P_{4}\right)$, $k=(2|V_{1}|,2|V_{2}|,2|V_{3}|,2|V_{4}|)$, and $\theta=\left(\theta_{1},\theta_{2},\theta_{3},\theta_{4}\right)$ are the voltage angles.

Let us consider the case with $P=(1.00,-1.25,2.00,-1.75)$ and $|V|=(1.10,0.93,1.05,0.90)$. That is, we aim to solve (3) where $\omega=(1.00,-1.25,2.00,-1.75)$ and $k=(2.20,1.86,2.10,1.80)$. By taking the $16$ possible sign patterns as the numbers

[TABLE]

some possibilities can immediately be ruled out:

• •

[math], $1$, $2$, $4$, $5$, $8$ by Prop. 3.3;

• •

[math], $1$, $2$, $3$, $4$, $5$ by Prop. 3.5.

We now consider the remaining possibilities starting from the largest:

• •

One equilibria resulting from each of the following: $15$, $14$, $13$, $12$, $11$, $10$;

• •

No equilibria resulting from $9\equiv(+1,-1,-1,+1)$;

• •

Two equilibria resulting from $7\equiv(-1,+1,+1,+1)$;

• •

No equilibria resulting from $6\equiv(-1,+1,+1,-1)$.

For example, since $9\equiv(+1,-1,-1,+1)$ yields no equilibria, Lemma 3.8 provides that $8$, $1$, and [math] also yield no equilibria while Lemma 3.10 provides that $5$ yields no equilibria. In summary, this particular case has a total of eight equilibria satisfying (3).

We now turn to consider a special case for which we can provide further optimizations:

1. IC4:

$k_{1}\geq k_{2}\geq\cdots\geq k_{n}$

We note that IC4 is independent of the implicitly assumed conditions IC1–IC3 so we will explicitly state when this condition is also required. With IC4, we provide a simplification of Lemma 3.10.

Lemma 3.12

Suppose that IC4 is satisfied and $\mu$ and $\nu$ are such that $\sigma_{\mu}=+1$ and $\sigma_{\nu}=-1$. Let $\sigma^{\prime}$ be the same as $\sigma$ except that \sigma_{\mu}^{\prime}=\sigma_{\nu}=-1\ and $\sigma_{\nu}^{\prime}=\sigma_{\mu}=+1$. If $\mu<\nu$ and $f_{\sigma}$ has no positive roots, then $f_{\sigma^{\prime}}$ also has no positive roots.

Proof: From IC4, and $\mu<\nu$, we have

[TABLE]

For $R\geq\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2}$,

[TABLE]

Hence, $f_{\sigma^{\prime}}(R)\leq f_{\sigma}(R)$. Therefore, the result follows from Lemma 3.7. $\Box$

Writing $\sigma$ as a binary number, Lemma 3.8 allows changing a “1” to a “0.” With IC4, Lemma 3.12 allows swapping a “0” and a “1” provided the “0” is on the right of “1.” Thus, with this understanding and ordering, we state the main optimization result.

Theorem 3.13

Assume that IC4 is satisfied.

Suppose $\sigma=(+1,\dots,+1)$ and $f_{\sigma}$ has no positive roots. Then, for every $\sigma^{\prime}\in\{-1,+1\}^{n}$, $f_{\sigma^{\prime}}$ has no positive roots. 2. 2.

Suppose $\sigma$ has exactly one entry which is $-1$ and that $f_{\sigma}$ has no positive roots. Then, $f_{\sigma^{\prime}}$ also has no positive roots for every $\sigma^{\prime}\in\{-1,+1\}^{n}$ which is smaller than $\sigma$ using the aforementioned binary representation. 3. 3.

Suppose that $\sigma$ has at least two entries equal to $-1$ and $f_{\sigma}$ has no positive roots. Let $\ell$ be the penultimate entry of a $-1$ in $\sigma$. Let $\sigma=(\rho_{1},\rho_{2})$ where $\rho_{1}=(\sigma_{1},\dots,\sigma_{\ell-1},-1)$ and $\rho_{2}=(\sigma_{\ell+1},\dots,\sigma_{n})$. Then, for every $\sigma^{\prime}=(\rho_{1},\rho_{2}^{\prime})\in\{-1,+1\}^{n}$ such that $\rho_{2}^{\prime}$ is smaller than $\rho_{2}$ using the aforementioned binary representation, $f_{\sigma^{\prime}}$ also has no positive roots.

Proof: We prove the three cases as follows.

This case follows immediately by repeated application of Lemma 3.8. 2. 2.

This case follows by alternately applying Case 1 to parts of $\sigma$ and Lemma 3.12. 3. 3.

This case follows by applying Case 2 to $\rho_{2}$.

$\Box$

The main benefit of this theorem is that it allows one to skip sequential sign cases by directly computing the next case that needs to be checked from the current case.

Example 3.14

To illustrate, suppose the input parameters satisfy IC4 and $f_{\sigma}$ has no positive roots for $\sigma=\left(+1,+1,-1,+1,-1,+1\right)\equiv 110101_{2}=53$. Theorem 3.13 shows that $f_{\sigma^{\prime}}$ also has no positive roots for the following sequential sign patterns $\sigma^{\prime}$:

[TABLE]

Furthermore, 48 can be immediately calculated from 53 by zeroing out everything from the next to last 0 onward, so that the five listed cases do not need to be considered at all.

The following utilizes these previous results assuming IC4 to more efficiently compute the set of all equilibria to (3). This depends on two algorithms: a root finding method that returns the set of all roots of $f_{\sigma}$ in an interval $I$, denoted $\mathbf{Solve}(f_{\sigma},I)$, and a method that returns a sign pattern in $\{-1,+1\}^{n}$ given a number $0\leq\iota\leq 2^{n}-1$, denoted $\mathbf{Convert}(\iota)$.

Algorithm 3.15** (Optimized)**

In:

$\omega\in\mathbb{R}^{n}$* and $k\in\mathbb{R}_{>0}^{n}$ satisfying IC1–IC4.*

Out:

$\Theta$, the set of equilibria satisfying OC1.

$\Theta\leftarrow\{\}$** 2. 2.

$\iota\leftarrow 2^{n}-1$** 3. 3.

While $\iota\geq 0$ do

1. (a)

$\sigma\leftarrow\mathbf{Convert}(\iota)$** 2. (b)

$I\leftarrow\left[\left(\dfrac{\omega_{n}}{k_{n}}\right)^{2},\;\;\left(\displaystyle\frac{1}{n}\sum_{\mu\in\sigma_{+}}k_{\mu}\right)^{2}~{}\right]$** 3. (c)

If $I$ is empty, then

1. i.

Decrement $\iota$ according to Theorem 3.13 2. ii.

Continue (go back to the start of Step 3) 4. (d)

If $\sum_{\mu=1}^{\ell}\sigma_{\mu}k_{\mu}\leq 0$ for all $\ell=1,2,\ldots,n$, then

1. i.

Decrement $\iota$ according to Theorem 3.13. 2. ii.

Continue (go back to the start of Step 3) 5. (e)

$f_{\sigma}\leftarrow-R+\frac{1}{n}\sum_{\mu=1}^{n}\sigma_{\mu}\sqrt{k_{\mu}^{2}R-\omega_{\mu}^{2}}$** 6. (f)

$\mathcal{R}\leftarrow\mathbf{Solve}(f_{\sigma},I)$** 7. (g)

If $\mathcal{R}=\emptyset$, then

1. i.

Decrement $\iota$ according to Theorem 3.13 2. ii.

Continue (go back to the start of Step 3) 8. (h)

For $R\in\mathcal{R}$ do

1. i.

Compute $\theta\in(-\pi,\pi]^{n}$ such that $\sin\theta_{\nu}=\frac{\omega_{\nu}}{k_{\nu}\sqrt{R}}\hbox{~{}and~{}}\mathrm{sign}\cos\theta_{\nu}=\sigma_{\nu}$ for $\nu=1,\dots,n$. 2. ii.

Add $\theta$ to $\Theta$ 9. (i)

$\iota\leftarrow\iota-1$**

Remark 3.16

In Algorithm 3.15, Steps 3b and 3d follow from Prop. 3.3 and 3.5, respectively.

Example 3.17

To illustrate, we apply Algorithm 3.15 to the setup from Ex. 3.2.

• •

$\iota=3$* yielding $\sigma=(+1,+1)$:*

• $\circ$

$I=[4,12.25]$**

• $\circ$

One positive root of $f_{\sigma}(R)=-R+\frac{1}{2}\left(\sqrt{25R-16}+\sqrt{4R-16}\right)$ on $I$, namely $R=10.25$.

• $\circ$

This yields the equilibrium $\theta=(0.2526,-0.6747)$.

• •

$\iota=2$* yielding $\sigma=(+1,-1)$:*

• $\circ$

$I=[4,6.25]$**

• $\circ$

One positive root of $f_{\sigma}(R)=-R+\frac{1}{2}\left(\sqrt{25R-16}-\sqrt{4R-16}\right)$ in $I$, namely $R=4.25$.

• $\circ$

This yields the equilibrium $\theta=(0.3985,-1.8158)$.

• •

$\iota=1$* yielding $\sigma=(-1,+1)$:*

• $\circ$

$I=[4,1]$* is empty*

• $\circ$

Theorem 3.13 removes the $\iota=0$ case.

In summary, there are two equilibria satisfying (3).

3.3 Performance

We implemented both Algorithms 3.1 and 3.15 in C++ using the C-XSC library [21] with the univariate solver being an interval Newton method [18, Chap. 6]. The implementation is available at http://dx.doi.org/10.7274/R09W0CDP. In this section, we benchmark the performance of this with the following methods for computing all equilibria to (3):

• •

solve (4) using Gröbner basis techniques in Macaulay2 [17];

• •

solve (4) using homotopy continuation in Bertini [4] as in [29];

• •

compute equilibria for (3) using elliptical continuation from [24].

We end with an example having $n=60$ that is easily solvable using Algorithm 3.15.

Comparison with computational algebraic geometry:

We use the following setup from [29] to compare with solving (4) using Macaulay2 and Bertini with serial computations. For each $n=3,\dots,12$, the natural frequencies are equidistant, namely $\omega_{\mu}=-1+(2\mu-1)/n$ for $\mu=1,\dots,n$, with uniform coupling $k=(\sqrt{1.5},\dots,\sqrt{1.5})$. To simplify the algebraic geometry computations using Macaulay2 and Bertini, we compute the equilibria as in [29] by setting $\theta_{n}=0$ ($s_{n}=0$ and $c_{n}=1$) with the results summarized in Table 1.

With Macaulay2, we simply computed the total number of complex solutions, i.e., the degree of the ideal generated by the polynomials in (4) when $s_{n}=0$ and $c_{n}=1$. Thus, one would need to perform additional computations to compute the number of real solutions. The symbol $\ddagger$ means that the computation did not complete within 48 hours.

With Bertini, we performed two different computations. The first was to directly solve (4) using regeneration [19] and the second utilized a parameter homotopy [32]. Both of these computations provide all real and non-real solutions to (4).

Although Bertini is parallelized and Algorithm 3.15 is parallelizable, we again note that the data in Table 1 is based on using serial processing. Nonetheless, this shows the advantage of using Algorithm 3.15 to compute all equilibria without needing to compute the non-real solutions of (4).

Comparison with elliptical continuation:

We next compare Algorithm 3.15 with the elliptical continuation method proposed in [24]. While having the advantage of being applicable to a more general setting of power flow equations, the elliptical continuation method in [24] comes with both theoretical and computational drawbacks relative to Algorithm 3.15 when considered in the context of the Kuramoto model. In contrast to Algorithm 3.15, there currently is no theoretical guarantee that the elliptical continuation method in [24] will compute all equilibria. Moreover, the computational speed of Algorithm 3.15 can be several orders of magnitude faster than the elliptical continuation method in [24]. Consider, for instance, a test case with $n=18$, $k=(1,\dots,1)$, and

[TABLE]

When interpreted as a power flow problem, this test case represents a power system composed of $18$ buses with fixed, unity voltage magnitudes and specified active power injections given by $\omega$ in normalized “per unit” values. The buses are completely connected by lines with unity reactance and zero resistance. While this is a very special example of a power system network, the corresponding test case enables comparison between Algorithm 3.15 and the elliptical continuation method in [24] in the context of the Kuramoto model.

A serial implementation of the elliptical continuation method in [24] in Matlab yielded $8538$ equilibria satisfying (3) in $1.935\times 10^{5}$ seconds ($53.77$ hours). For a fair comparison, we used a serial implementation of Algorithm 3.15 in Matlab which computed $8538$ equilibria in $13.9$ seconds. Hence, the implementation of Algorithm 3.15 in Matlab is roughly four orders of magnitude faster than the Matlab implementation of [24] for this example. We note that the C++ implementation of Algorithm 3.15 took $6.6$ seconds.

An example with $n=60$:

We conclude with an example solved by Algorithm 3.15 for $n=60$ having $k=(60,\dots,60)$ and

[TABLE]

This example has 2 equilibria satisfying (3) with the total computation time using the C++ implementation of Algorithm 3.15 taking under a second. For comparison, the elliptical continuation method as described in the previous example took $5609$ seconds ($93.5$ minutes). This example is simply too large for current methods that compute all complex roots. Generally, problems with $\left(\dfrac{w_{n}}{k_{n}}\right)^{2}$ near $\left(\dfrac{1}{n}\displaystyle\sum_{\mu=1}^{n}k_{\mu}\right)^{2}$ will be solved quickly by Algorithm 3.15 as a consequence of Prop. 3.3 and Theorem 3.13.

4 Counting equilibria

After reviewing known information, we compute the generic root count for (4) which bounds the number of equilibria to (3). By analyzing the number of equilibria in particular cases, we can asymptotically compare the maximum number of equilibria to the generic root count of (4).

4.1 Summary of known results

As mentioned in the Introduction, the arbitrary coupling case (5) has at most $\binom{2n-2}{n-1}$ equilibria [2] and, for $n\geq 4$, it is currently unknown if this bound can be achieved. The minimum number of equilibria is easily observed to be [math].

There are results regarding the number of equilibria for the standard Kuramoto model that apply to the rank-one coupled Kuramoto model as well. When $n=2$, it is easy to see that the maximum number of equilibria satisfying (3) is $2$. By IC1, we have $\omega_{2}=-\omega_{1}\neq 0$, so, without loss of generality, we assume $\omega_{1}>0$. With $k=(1,1)$, one can verify:

• •

$2$ equilibria if $0<\omega_{1}<\dfrac{1}{2}$;

• •

$1$ equilibrium (of “multiplicity 2”) if $\omega_{1}=\dfrac{1}{2}$;

• •

[math] equilibria if $\omega_{1}>\dfrac{1}{2}$.

For $n=3$, the maximum number of equilibria is $6$ [2, 40]. When $k=(1,1,1)$, Prop. 3.3 shows that equilibria can only occur if each $|\omega_{\nu}|\leq 1$. By taking $\omega_{3}=-\omega_{1}-\omega_{2}$ due to IC1, Figure 2 plots the regions having [math], $2$, $4$, and $6$ distinct equilibria for $\omega_{1},\omega_{2}\in[-1,1]$. Such a plot has appeared previously, e.g., [5, 20].

For $n=4$, the maximum number of equilibria is $14$ [2, 40] and it is an open problem to determine if this bound is sharp. A recent experiment [40] applied to the standard Kuramoto model computed all equilibria for selected values of $\omega\in\mathbb{R}^{4}$ in a relevant compact parameter space based on a grid with step-size $1/2000$. Since this experiment attained a maximum of $10$ equilibria, they conjecture that the maximum number of equilibria satisfying (3) when $n=4$ and $k=(\sqrt{K},\sqrt{K},\sqrt{K},\sqrt{K})$ is $10$, which is strictly smaller than the upper bound of $14$. We revisit this case in Ex. 4.8 and 4.10.

4.2 Bounding the number of equilibria

As summarized in Section 4.1, the maximum number of equilibria to (3) is $2,6,14$ for $n=2,3,4$, respectively. Theorem 4.3 shows that $2^{n}-2$ bounds the number of equilibria with Corollary 4.4 showing that $2^{n}-2$ is actually the generic root count for the polynomial system (4) modulo shift.

Let $\omega\in\mathbb{R}^{n}$ and $k\in\mathbb{R}_{>0}^{n}$ satisfy IC1-IC3. The following shows that the function

[TABLE]

is actually a reducible polynomial.

Proposition 4.1

The univariate function $g$ in (9) is a polynomial of degree $2^{n}$. Moreover, there exists a polynomial $h(R)$ of degree $2^{n}-2$ with

[TABLE]

Proof: Since $g$ is a product over all $2^{n}$ conjugates, it immediately follows that $g$ is a polynomial with leading term $(-R)^{2^{n}}$ showing that $g$ is a polynomial of degree $2^{n}$.

In order to show that $R^{2}$ is a factor of $g$, we simply need to show that $g(0)=g^{\prime}(0)=0$. To that end, consider $\sigma_{\omega}={\rm sign~{}}\omega\in\{-1,+1\}^{n}$ where ${\rm sign~{}}\omega_{i}=1$ if $\omega_{i}\geq 0$, otherwise ${\rm sign~{}}\omega_{i}=-1$. Then,

[TABLE]

by IC1. By (9), this immediately shows that $g(0)=0$ since one of the terms in the product is [math].

By a similar argument as above, $f_{-\sigma_{\omega}}(0)=0$ by IC1. This shows that at least two terms in the product defining $g$ in (9) are zero. Hence, the product rule for differentiation shows that $g^{\prime}(0)=0$. $\Box$

Example 4.2

For $n=2$, we have

[TABLE]

which is indeed a polynomial of degree $2^{2}=4$. Moreover, IC1 implies $\omega_{2}=-\omega_{1}$ so that

[TABLE]

Proposition 4.1 immediately provides the following upper bound.

Theorem 4.3

If $\omega\in\mathbb{R}^{n}$ and $k\in\mathbb{R}_{>0}^{n}$ satisfy IC1–IC3, then there are at most $2^{n}-2$ equilibria satisfying (3).

Proof: This follows from Theorem 2.2 since $g(R)$ in (9) has at most $2^{n}-2$ positive roots. $\Box$

Corollary 4.4

The generic root count modulo shift to (4) is $2^{n}-2$.

Proof: Reviewing the proof of Theorem 2.2 shows that $2^{n}-2$ also bounds the number of complex solutions to (4). For $g$ in (9), $g^{\prime\prime}(0)\neq 0$ for generic values of the parameters yielding that there are generically $2^{n}-2$ nonzero roots of $g$. Hence, $2^{n}-2$ is the generic root count of (4). $\Box$

Example 4.5

Table 1 shows that the polynomial system (4) for $n=4$, $\omega=(-3/4,-1/4,1/4,3/4)$, and $k=(\sqrt{1.5},\sqrt{1.5},\sqrt{1.5},\sqrt{1.5})$ has $12$ complex roots modulo shift, which is less than the generic root count of $2^{4}-2=14$. In fact, as in the proof of Prop. 4.1, this is due to the following four quantities being equal to zero:

[TABLE]

Hence, $g$ in (9) has $g(0)=g^{\prime}(0)=g^{\prime\prime}(0)=g^{\prime\prime\prime}(0)=0$, namely

[TABLE]

Theorem 4.3 provides an upper bound of $2^{n}-2$ when the symmetric coupling matrix has rank one while [2] provides an upper bound of $\binom{2n-2}{n-1}$ in the general case. By Stirling’s formula,

[TABLE]

showing the bound in Theorem 4.3 for the rank-one case is roughly the square root of the general purpose bound from [2]. Due to this difference, we computed the generic root counts for the corresponding polynomial system associated with (5) when the coupling matrix $\kappa$ is a symmetric matrix of various ranks for $n=2,\dots,10$ using Bertini [4]. The results are presented in Table 2. This data, for selected values of $r$ and $n$, shows that the generic root counts for a symmetric coupling matrix of rank $r$ and rank $r+1$ are equal whenever $n\leq 2r+1$ and differ when $n\geq 2r+2$. In fact, the difference between the generic root counts for rank $r$ and rank $r+1$ symmetric coupling matrices when $n=2r+2$ is equal to $\binom{2r+2}{r+1}=\binom{n}{n/2}$. We leave it as an open problem to fully understand the behavior for all choices of $r$ and $n$.

4.3 Counting equilibria for particular cases

Motivated by [40], we use Theorem 2.2 to analyze the number of equilibria satisfying (3) for particular cases when $n$ is even (Theorem 4.6 and Corollary 4.9) and when $n$ is odd (Theorem 4.12).

Theorem 4.6

Suppose that $n\geq 2$ is even and $q>0$. For $\omega=\left(nq,\ldots,nq,-nq,\ldots,-nq\right)$ and $k=(n,\ldots,n)$, there are exactly

[TABLE]

equilibria satisfying (3) counting multiplicity. Hence, the number of equilibria changes precisely at the integers $q=1,2,\dots,n/2$.

Proof: Since $k_{\mu}^{2}=n^{2}$ and $\omega_{\mu}^{2}=n^{2}q^{2}$, Theorem 2.2 shows that we need to compute all $R>0$ where

[TABLE]

with $S=\sum_{\mu=1}^{n}\sigma_{\mu}$ and $\sigma\in\{-1,+1\}^{n}$.

If $S\leq 0$, then (11) has no positive solutions. Since $n$ is even, the remaining cases have $S\geq 2$. Thus, the positive solutions of (11) must satisfy

[TABLE]

This yields three cases:

$2\leq S<2q$: (11) has no positive solutions; 2. 2.

$S=2q\geq 2$: (11) has one positive solution of multiplicity 2, namely $R=S^{2}/2$; 3. 3.

$S>2q$ with $S\geq 2$: (11) has two distinct positive solutions.

Suppose that $q$ is not an integer. Since $S$ is even, we have $S\neq 2q$. Hence, the number of equilibria is exactly

[TABLE]

Since $\binom{n}{n/2+\ell}=\binom{n}{n/2-\ell}$ and $2^{n}=\sum_{\ell=0}^{n}\binom{n}{\ell}$, the number of equilibria when $q$ is not an integer is

[TABLE]

When $q$ is an integer, we need to add in the case when $S=2q$ yielding

[TABLE]

$\Box$

Example 4.7

For $n=2$ and $q>0$, the case of $\omega=(2q,-2q)$ and $k=(2,2)$ corresponds to $\omega=(q/2,-q/2)$ and $k=(1,1)$ as considered in Section 4.1. Hence, counting multiplicity, there are two equilibria for $q\leq 1$ and no equilibria for $q>1$ in agreement with Theorem 4.6.

Example 4.8

For $n=4$ and $q>0$, the case of $\omega=(4q,4q,-4q,-4q)$ and $k=(4,4,4,4)$ corresponds to $\omega=(q/4,q/4,-q/4,-q/4)$ and $k=(1,1,1,1)$ as considered in Section 4.1. Figure 3(a) plots the regions based on the number of equilibria when $k=(1,1,1,1)$ such that $\omega_{3}=\omega_{4}=-(\omega_{1}+\omega_{2})/2$. With this setup, $\omega_{1}=\omega_{2}=q/4$ implies $\omega_{3}=\omega_{4}=-q/4$. Since the sign is arbitrary, the plot in Figure 3(b) incorporates the line $\omega_{1}=\omega_{2}=q/4$. By Theorem 4.6, there are $10$ equilibria for $0<|q|<1$, $2$ equilibria for $1<|q|<2$, and no equilibria for $|q|>2$.

Theorem 4.6 immediately yields the following.

Corollary 4.9

Suppose that $n\geq 2$ is even and $q>0$. The maximum number of distinct equilibria satisfying (3) when $\omega=(nq,\dots,nq,-nq,\dots,-nq)$ and $k=(n,\dots,n)$ is

[TABLE]

which occurs for all $0<q<1$.

Example 4.10

For $n=4$, Corollary 4.9 provides a maximum of $2^{4}-\binom{4}{2}=10$ distinct equilibria which matches the computational results in [40] as discussed in Section 4.1.

Before considering the odd case, we first define the constants

[TABLE]

and prove an inequality regarding them.

Lemma 4.11

For $0<q<q_{o}$, $R_{o}+\sqrt{R_{o}}-2\sqrt{R_{o}-q^{2}}<0$ where $q_{o}$ and $R_{o}$ as defined in (13).

Proof: Since $q<q_{o}$ and $R_{o}-q^{2}>R_{o}-q_{o}^{2}>0$, we have

[TABLE]

$\Box$

With Lemma 4.11, we now consider the case when $n$ is odd.

Theorem 4.12

Suppose that $n\geq 3$ is odd and let $0<q<q_{o}$ where $q_{o}$ is defined by (13). For $\omega=(nq,\dots,nq,-nq,\dots,-nq,0)$ and $k=(n,\ldots,n)$, the number of equilibria satisfying (3) is

[TABLE]

Proof: Since $k_{\mu}^{2}=n^{2}$ for $\mu=1,\dots,n$, $\omega_{\nu}^{2}=n^{2}q^{2}$ for $\nu=1,\dots,n-1$, and $\omega_{n}=0$, Theorem 2.2 shows that we need to compute all $R>0$ with

[TABLE]

where $S=\sum_{\mu=1}^{n-1}\sigma_{\mu}$ and $\sigma\in\{-1,+1\}^{n}$. Define $p_{\sigma}(R)=R-\sigma_{n}\sqrt{R}-S\sqrt{R-q^{2}}$.

Since $n-1$ is even, we know that $S$ is also even. This yields three cases to consider.

$S<0$:

Rewriting (15) as

[TABLE]

shows that the right-hand size is non-positive. Hence, to have a solution, we need $\sigma_{n}=+1$ and $R\in(q^{2},1)$. Since $p_{\sigma}(q^{2})=q^{2}-q<0$ and $p_{\sigma}(1)=-S\sqrt{1-q^{2}}>0$, we know that there is at least one root in $(q^{2},1)$. In fact, since $S\leq-2$, it is easy to see that $p_{\sigma}$ is a strictly increasing function on $(q^{2},1)$ since

[TABLE]

for all $R\in(q^{2},1)$. Thus, this case yields one equilibrium for each $\sigma\in\{-1,+1\}^{n}$ such that $\sigma_{n}=+1$ and $S<0$ for a total of

[TABLE]

$S=0$:

Since (15) becomes $R=\sigma_{n}\sqrt{R}$, this case requires $\sigma_{n}=+1$ and $R=1$. The total number of equilibria for this case is thus

[TABLE]

$S>0$:

We split this into two cases based on the value of $\sigma_{n}$.

$\sigma_{n}=+1$:

Rewriting (15) as

[TABLE]

shows that the right-hand size is nonnegative. Hence, to have a solution, we need $R>1$. Since $p_{\sigma}(1)=-S\sqrt{1-q^{2}}<0$ and $\lim_{R\rightarrow\infty}p_{\sigma}(R)=\infty$, we know that there is at least one root in $(1,\infty)$. In fact, the root is unique since the graph of $p_{\sigma}$ is concave up due to

[TABLE]

for $R>1$. Hence, the total number of equilibria for this case is

[TABLE]

$\sigma_{n}=-1$:

We need to compute the number of roots of $p_{\sigma}$ for $R>q^{2}$. Since $S\geq 2$ and $R^{3/2}>(R-q^{2})^{3/2}$ for all $R>q^{2}$, it follows that

[TABLE]

when $R>q^{2}$. Hence, $p_{\sigma}$ is concave up on $R>q^{2}$ with $p_{\sigma}(q^{2})=q^{2}+q>0$ and $\lim_{R\rightarrow\infty}p_{\sigma}(R)=\infty$. Thus, the number of roots depends on the sign of the minimum value of $p_{\sigma}$ on $R>q^{2}$. Since increasing $S$ makes $p_{\sigma}$ more negative and Lemma 4.11 shows that $p_{\sigma}(R_{o})<0$ when $S=2$, there are always two roots with $R>q^{2}$. Hence, the total number of equilibria for this case is

[TABLE]

The result is obtained by simply summing the number of equilibria from all of these cases. $\Box$

Example 4.13

For $n=3$, Theorem 4.12 shows that the number of equilibria for $\omega=(3q,-3q,0)$ and $k=(3,3,3)$ is $2^{3}-\binom{2}{1}=6$ whenever $0<q<q_{o}$ with $q_{o}$ defined in (13). This is equivalent to the case when $\omega=(q/3,-q/3,0)$ and $k=(1,1,1)$ for $0<q<q_{o}$. Since the ordering of the elements in $\omega$ is arbitrary, Figure 4 is an enhanced version of Figure 2 that plots, in red, the corresponding three segments within the region having $6$ equilibria:

• •

$\{(\alpha,0,-\alpha)~{}:~{}|\alpha|<q_{0}\}$* is the horizontal segment,*

• •

$\{(0,\alpha,-\alpha)~{}:~{}|\alpha|<q_{0}\}$* is the vertical segment, and*

• •

$\{(\alpha,-\alpha,0)~{}:~{}|\alpha|<q_{0}\}$* is the diagonal segment.*

The following suggests an upper bound on the maximum number of equilibria.

Conjecture 4.14

For $n\geq 2$, the maximum number of equilibria satisfying (3) with $n$ oscillators is

[TABLE]

which are achieved in Corollary 4.9 and Theorem 4.12, respectively.

As summarized in Section 4.1, this conjecture matches the known cases of $n=2$ and $n=3$, and agrees with the conjecture for $n=4$ provided in [40] for the standard Kuramoto model.

4.4 Asymptotic behavior

Even though we can only conjecture an upper bound on the number of equilibria, the results from Corollary 4.9 and Theorem 4.12 provide the following result: there can asymptotically be as many equilibria satisfying (3) as the number of complex solutions to (4) modulo shift.

Theorem 4.15

As $n\rightarrow\infty$, the ratio of the maximum number of equilibria satisfying (3) and the generic root count to (4) limits to $1$.

Proof: For each $n\geq 2$, let $\Omega(n)$ denote this ratio. Theorems 4.3 and 4.12 together with Corollaries 4.4 and 4.9 show that, for every $\ell\geq 1$,

[TABLE]

Stirling’s formula yields

[TABLE]

so that

[TABLE]

Similarly, Stirling’s formula yields

[TABLE]

so that

[TABLE]

Therefore, $\Omega(n)\rightarrow 1$ as $n\rightarrow\infty$. $\Box$

5 Conclusion

The Kuramoto model is a standard model used to describe the behavior of coupled oscillators which has proven to be useful in many applications, e.g., electrical engineering [15, 39], biology [36], and chemistry [3, 23, 33]. When the coupling matrix is a symmetric matrix of rank one, which is a slight generalization of the standard Kuramoto model (1), the reformulation (Theorem 2.2) permits all equilibria to be computed efficiently and effectively (Section 3.3) without the need to compute all complex solutions to a corresponding polynomial system. Moreover, this reformulation is also useful for computing an upper bound on the number of equilibria (Theorem 4.3), computing the exact number of equilibria for particular cases (Theorem 4.6 and Theorem 4.12), and understanding the asymptotic behavior of the maximum number of equilibria (Theorem 4.15).

Even with the broad use of the Kuramoto model and the new results presented in this paper regarding the equilibria, many questions still remain. One prominent question is how to compute the maximum number of equilibria when the coupling matrix has rank one, which we have conjectured (Conjecture 4.14) is strictly smaller than the upper bound of $2^{n}-2$ for all $n\geq 4$, an extension of the computational results for the standard Kuramoto when $n=4$ from [40]. It is also unknown how to extend our decoupling method to higher rank models. One final question regards the relationship between the rank of the coupling matrix, the number of oscillators, and the number of equilibria (Table 2), which may yield new approaches for computing all equilibria when the coupling matrix has rank $r>1$.

Acknowledgment

We would like to thank Dhagash Mehta for helpful discussions regarding the Kuramoto model, and Bernard Lesieutre and Dan Wu for sharing a Matlab implementation of their elliptical continuation method proposed in [24].