Optimizing synchronization stability of the Kuramoto model in complex networks and power grids

Maintaining the stability of synchronization state is crucial for the functioning of many natural and artificial systems. In this study, we develop methods to optimize the synchronization stability of the Kuramoto model by minimizing the dominant Lyapunov exponent. Using the recently proposed cut-set space approximation of the steady states, we greatly simplify the objective function, and further derive its gradient and Hessian with respect to natural frequencies, which leads to an efficient algorithm with the quasi-Newton's method. The optimized systems are demonstrated to achieve better synchronization stability for the Kuramoto model with or without inertia in certain regimes. Hence our method is applicable in improving the stability of power grids. It is also viable to adjust the coupling strength of each link to improve the stability of the system. Various operational constraints can also be easily integrated into our scope by employing the interior point method in convex optimization. The properties of the optimized networks are also discussed.

Figure 1

${\varphi }_j-{\varphi }_i$ vs. $sin({\theta }_j^{{}^{*}}-{\theta }_i^{*})$. (a) Erdös-Rényi graph of 50 nodes (ER50), where $\omega$ is drawn from a Gaussian distribution and $K_{ij}=1$. Inset: root-mean-square error (RMSE) of estimator ${\varphi }_j-{\varphi }_i$ for $sin({\theta }_j^{*}-{\theta }_i^{*})$ among all the edges. Each data point is averaged over 100 samples. (b) IEEE reliability test system 96 (RTS96) [19], where $\omega$ is modified from the power injection data in the test system and $K_{ij}$ is defined to be the inverse of line reactance of edge $(i,j)$. Inset: RMSE of estimator ${\varphi }_j-{\varphi }_i$ for $sin({\theta }_j^{*}-{\theta }_i^{*})$ among all the edges. Each data point is averaged over 100 samples.

Figure 2

${\lambda }_2$ and ${\parallel \omega \parallel }_2$ through optimization for RTS96 power network. The initial natural frequency is modified from the power injection data in the test case. (a) Gradient ascent update. Both ${\lambda }_2(\tilde{L}\left(\varphi \right))$ and ${\lambda }_2(L\left({\theta }^{*}\right))$ increase gently in the later stage, and the natural frequency $\omega$ is approaching the optimal state $\omega =0$ very slowly due to the flat landscape. (b) Quasi-Newton update. Both ${\lambda }_2(\tilde{L}\left(\varphi \right))$ and ${\lambda }_2(L\left({\theta }^{*}\right))$ approach the optimum ${\lambda }_2(L({\theta }^{*}=0))=0.6889$ after six iterations.

Figure 3

(a) ${\lambda }_2$ through optimization for the RTS96 power network under Euclidean norm constraint. (b) The unoptimized and optimized system, where white square nodes have positive natural frequencies (generators) while gray circular nodes have nonpositive natural frequencies (loads or relay nodes). Edge color intensity encodes $cos({\theta }_i^{*}-{\theta }_j^{*})$. (c) Response of the RTS 96 power network governed by the first-order Kuramoto model. (d) Response of the RTS 96 power network governed by the second-order Kuramoto model with unit damping $D_i=1$ and small inertia $M_i=0.2$. In both (c) and (d), the disturbance, drawn from the Gaussian distribution with mean zero and standard deviation $0.05\text{rad}$, was applied to the steady state of phase oscillators at $t=0$.

Figure 4

Properties of the optimal system compared to systems with random frequencies. The networks are 100 realizations of ER random graphs with 50 nodes and edge connection probability $p=0.1$. (a) Correlation of increment of ${\lambda }_2(L\left({\theta }^{*}\right))$ and order parameter $r$. (b) Histogram of changes of phase angle differences among all edges $(i,j)$ in all realizations. (c) Average neighbor frequency ${\langle \omega \rangle }_i={\sum }_{j\in \partial i}{\omega }_j/d_i$ vs. natural frequency ${\omega }_i$. (d) Alignments of natural frequencies with graph Laplacian eigenvectors, i.e., $|\langle v_i|\omega \rangle {|}^2$ where $v_i$ is the normalized eigenvector corresponding to the $i\mathrm{th}$ smallest eigenvalue. The data was divided into 10 bins and $|\langle v_i|\omega \rangle {|}^2$ was first summed inside every bin for each sample, after which the sample mean and standard deviation of the bin summation quantity ${\sum }_{i\in \text{bin}}{\left|\langle v_i|\omega \rangle \right|}^2$ was calculated.

Figure 5

(a) Eigenvector $v_2$ corresponding to the second smallest eigenvalue ${\lambda }_2$ of the graph Laplacian matrix of a specific ER random graph, depicted on the network. The color indicates the sign of $v_{2j}$ on node $j$, i.e., white node corresponds to $v_{2j}>0$, while gray node corresponds to $v_{2j}\le 0$. The size of the node indicates the strength of $|v_{2j}|$ on that node. (b) Eigenvector $v_N$ corresponding to the largest eigenvalue ${\lambda }_N$ of the graph Laplacian matrix depicted on the network. (c) Frequencies ${\omega }^{\text{opt}}$ corresponding to the algebraic connectivity of the state-dependent Laplacian matrix evaluated at the optimal algebraic connectivity. Similarly, the white nodes correspond to ${\omega }^{\text{opt}}>0$, while the gray nodes correspond to ${\omega }^{\text{opt}}\le 0$.

Figure 6

Phase angles ${\theta }^{*}$ and state-dependent edge weights $W\left({\theta }^{*}\right)$ in a two-module network. In both cases, the $L^2$-norm of natural frequency is ${\parallel \omega \parallel }_2=4.26$. (a) Phases of the system depicted on the unit circle in Case 1. (b) The state-dependent edge weight $W{\left({\theta }^{*}\right)}_{ij}=K_{ij}cos({\theta }_i^{*}-{\theta }_j^{*})$ in Case 1. (c) Phases of the system depicted on the unit circle in Case 2. (d) The state-dependent edge weight $W{\left({\theta }^{*}\right)}_{ij}=K_{ij}cos({\theta }_i^{*}-{\theta }_j^{*})$ in Case 2.

Figure 7

(a) ${\lambda }_2$ and ${\parallel \omega \parallel }_1$ through optimization for the RTS96 power network under linear constraints at $\alpha =0.1$. The $L_1$-norm of natural frequency ${\parallel \omega \parallel }_1:={\sum }_i\left|{\omega }_i\right|$ is twice of total generation or total consumption. (b) ${\lambda }_2(L\left({\theta }^{*}\right))$ and ${\parallel \omega \parallel }_1$ of the optimal system as a function of $\alpha$ with variable demand (Problem 1) and fixed demand (Problem 2).

Figure 8

Optimizing the state algebraic connectivity by updating coupling strengths. (a) State algebraic connectivity ${\lambda }_2(L\left({\theta }^{*}\right))$, ${\lambda }_2(\tilde{L}\left(\varphi \right))$, and graph algebraic connectivity ${\lambda }_2\left(L\left[K\right]\right)$ through optimization. The initial state is the same as Case 1 in Sec. 4c. (b) The state-dependent edge weight $W{\left({\theta }^{*}\right)}_{ij}=K_{ij}cos({\theta }_i^{*}-{\theta }_j^{*})$ in the optimal state. Note the scale of color code is different from the cases of Fig. 6.