Analysis of the Dynamics and Topology Dependencies of Small Perturbations in Electric Transmission Grids
{Luiscarlos~A.~Torres-S\'anchez
Giuseppe T. Freitas de Abreu
~\IEEEmembership{Senior Member~IEEE
We study the dynamics of disturbances in unregulated electric transmission grids through a Synchronous Motor load model. We start our analysis with linearized system equations, which we represent as complex Fourier series to find their eigenmodes and eigenfrequencies. This reduces the problem to the diagonalization of a finite dimensional matrix, which depends on the stationary phase solutions of the grid, and is thereby inherently conditioned by the topology and the power distribution. This matrix is found to belong to the generalized Laplacian matrices, which relates the analysis of perturbation dynamics to a graph theory problem. We consider three networks: Small-world, Random and German transmission grid. We find that the density of eigenfrequencies highly depends on the topology. For the Random network, it resembles the Mar\u{chenko-Pastur distribution. For the Small-world network and the German transmission grid, we find strongly peaked densities with long tails. Moreover, we find that the algebraic connectivity of the generalized Laplacian matrix determines the long-time transient behavior of perturbations, whereas the intensity of its eigenvector, the generalized Fiedler vector, discloses information of the modes localization. Finally, we find strong indications that the degree of localization tends to increase with a decrease of the network's average clustering coefficient.
bibliography.tex
\begin{thebibliography}{1}
%1-5

\bibitem{bna}
Versorgungsqualit\"at-SAIDI-Werte 2006-2014, Bundesnetzagentur, 2015.


\bibitem{meier}
   A. von Meier, Electric Power Systems: A Conceptual Introduction. John Wiley Sons, Inc, 2006.
   

\bibitem{kundur_ieee} P. Kundur, J. Paserba and V. Ajjarapu, ``Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions'', {\it IEEE Transactions on  Power Systems}, vol. 19, no. 3, pp. 1387 - 1401, Aug. 2004.



\bibitem {ulbig}
	A. Ulbig, T. Borsche and G. Andersson, arXiv:1312.6435 {\bf[math.OC]} (2014).

  \bibitem{cuadra}
 L. Cuadra, S. Salcedo-Sanz, J. del Ser and Z. Woo Geem. A critical review of robustness in power grids using complex networks concepts. Energies {\bf 8}, 9211-9265 (2015).

 \bibitem{timme}
 %  Self-Organized Synchronization in Decentralized Power Grids,
M. Rohden, A. Sorge, M. Timme and D. Witthaut, Phys. Rev. Lett.
{\bf 109}, 064101 (2012).


\bibitem{jung}
M. Rohden,  D. Jung, S. Tamrakar and S. Kettemann, Cascading failures in ac electricity grids, Phys. Rev. {\bf E 94}, 032209 (2016).  


%6-10

\bibitem{machowski} J. Machowski, J. Bialek, J. Bumby, Power System Dynamics: Stability and Control,  Wiley, 2008. 

 
 \bibitem{coletta}
 T. Coletta, P. Jacquod,
  Linear Stability and the Braess Paradox in Coupled-Oscillator Networks and Electric Power Grids, Phys Rev E {\bf 93}, 032222 (2016).
 



\bibitem{hill} A. Bergen and D. Hill, ``A structure preserving model for power system stability analysis'', {\it IEEE Trans. on Power Apparatus and Systems}, vol. PAS-100, no.1, pp. 25-35, Jan. 1981.


\bibitem{filatrella}  G. Filatrella, A. Nielsen, and N. Pedersen, Eur Phys. J B {\bf 61}, 485 (2008).


 %11-15

\bibitem{schmietendorf} K. Schmietendorf, J. Peinke, R. Friedrich and O. Kamps, 
 Eur. Phys. J. Spec. Top. {\bf 223}, 2577 (2014). 
 
 
 \bibitem{heuck} K. Heuck, K.-D. Dettmann and D. Schulz,  
 {\it Elekrtische Energieversorgung}, Springer, Wiesbaden (2013). 
 

  \bibitem{menck} P. J. Menck, J. Heitzig, J. Kurths and H. Schellnhuber,
%How dead ends undermine power grid stability
Nature Communications, 5, 3969 (2014).



 \bibitem{motter}
T. Nishikawa, A. Motter, New J. Phys. {\bf 17}, 015012 (2015).
 
 \bibitem{kundur_book} P. Kundur, Power System Stability and Control, Mc Graw Hill (1994). 

 %16-20

\bibitem{kettemann}

S. Kettemann, Delocalization of disturbances and the stability of AC electricity grids, Phys. Rev. E {\bf 94}, 062311 (2016).


  \bibitem{myers}
 A. Motter, S. Myers, M. Anghel and Takashi Nishikawa, Nature Physics {\bf 9}, 191-197 (2013).
 
 \bibitem{overbye}
T. Overbye, X. Cheng, Y. Sun, ``A comparison of the AC and DC power flow models for LMP Calculations'', {\it Proceedings of the 37th Annual Hawaii International Conference on System Sciences}, 2004.

 
   \bibitem{watts_strogatz}
 D. Watts and S. Strogatz, Collective dynamics of small-world networks, Nature {\bf 393}, 440-442 (1998).

 

 \bibitem{barabasi}
 R. Albert and A. Barab\'asi, Statistical mechanics of complex networks, Reviews of Modern Physics, vol. 74, no.1, 47-97 (2002).


 %21-25
 
  \bibitem{brenton} 
  
  B. Prettejohn, M. Berryman and M. McDonnell, Methods for generating complex networks with selected structural properties for simulations: a review and tutorial for neuroscientists, Frontiers in Computational Neuroscience, vol. 5, 5-11 (2011).
 
 
   \bibitem{hu}
 L. Fu, W. Huang, S. Xiao, Y. Li and S. Guo. Vulnerability assessment for power grid based
on small-world topological model. {\it Pacific Power and Energy Engineering Conference (APPEEC), 2010 Asia-Pacific}, 2010.

  
   \bibitem{pagani}

G. Pagani and M. Aiello. Power grid complex network evolutions for the smart grid. Phys.
A Stat. Mech. Its Appl. , 396, 248-266, 2014.

  \bibitem{scigrid}
 W. Medjroubi, C. Matke and D. Kleinhans, SciGRID an open source reference model for the European transmission network, version 0.1, http://www.scigrid.de/ (2015).
 
  \bibitem{dorfler}
 F. Dorfler, M. Chertkov and F. Bullo. Synchronization in complex oscillator networks and
smart grids, {\it Proceedings of the National Academy of Sciences of the United States of America}, vol. 110, no. 6, pp. 2005-2010 (2013).


  %26-30
  \bibitem{hogben}
  
  L. Hogben, Spectral graph theory and the inverse eigenvalue problem of a graph, Electronic Journal of Linear Algebra, ISSN 1081-3810, vol. 14,  pp. 12-31, Jan. 2005. 
   
   
 \bibitem{fiedler}
M. Fiedler, Algebraic connectivity of graphs. Czechoslovak Math. J. 23 (1973), 298-305.


 %31


\bibitem{kirkland}

S. Kirkland, I. Rocha and V. Trevisan, Algebraic connectivity of k-connected graphs, Czechoslovak Mathematical Journal, vol.65, no.1, pp. 219-236, Apr. 2015.


\bibitem{bermann}
A. Bermann, X-D. Zhang, Lower bounds for the eigenvalues of Laplacian matrices, Linear Algebra and its Applications, vol. 316, 1-3, pp. 13-20, Sept. 2000.

 \bibitem{ministry}
 
 German Federal Ministry for Economic Affairs and Energy. Available in https://www.bmwi.de/Redaktion/EN/Artikel/Energy/strommarkt-der-zukunft-zahlen-und-fakten.html (2017).


  \bibitem{birron}
 
 D. Kothari and I. Nagrath, Power System Engineering 2nd edition, The McGraw Hill Companies (2008).
 
 
 \bibitem{buzna}
L. Buzna, S.Lozano and A. D\'iaz-Guilera. Synchronization in symmetric bipolar population
network, Physical Review E {\bf 80}, 066120, 2009.

  \bibitem{dyson} F. J. Dyson, Phys. Rev.  {\bf 92}, 1331 (1953).
  
  \bibitem{alexander}  S. Alexander, J. Bernasconi, W. R. Schneider, R. Orbach, Rev. Mod. Phys. {\bf 53}, 175 (1981).
   
    \bibitem{ziman}
   T. A.  L. Ziman, Phys. Rev. Lett. {\bf 49}, 337 (1982).
   
    \bibitem{wegner}
 F. J. Wegner in   E. Abrahams, 50 years of Anderson Localization, World Scientific (2010). 
 
  \bibitem{john}
 S. John, H. Sompolinsky, M.J. Stephen, Phys. Rev. B {\bf 27}, 5592 (1983).
 
   \bibitem{anderson}
 P. W. Anderson, Phys. Rev. {\bf 109}, 1492 (1958).

 \bibitem{derong}
 D. Liu, P. Antsaklis, Stability and Control of Dynamical Systems with Application, p. {\bf 328}, Springer, 2003.

  \bibitem{sari}
 A. T. Sari\' c, A. M. Stankovi\' c, "Rapid small-signal stability
assessment and enhancement following changes in topology", {\it IEEE Transactions on  Power Systems}, vol. 30, no. 3, pp. 1155–1163, May, 2015.

 \bibitem{wang}
 W. Wang, G. M. Huang, Impacts of smart grid topology control on power system stability with renewable integration, {\it North American Power Symposium (NAPS)}, 2016.
 
 



%\bibitem{verdu}

%Verdu



\end{thebibliography}

% biography section
% 
% If you have an EPS/PDF photo (graphicx package needed) extra braces are
% needed around the contents of the optional argument to biography to prevent
% the LaTeX parser from getting confused when it sees the complicated
% \includegraphics command within an optional argument. (You could create
% your own custom macro containing the \includegraphics command to make things
% simpler here.)
%\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{mshell}}]{Michael Shell}
% or if you just want to reserve a space for a photo:

%\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.5in,clip,keepaspectratio]{lats_photo}}]{Luiscarlos A. Torres-S\'anchez} was born in Managua, Nicaragua, in 1993. He received the B.Eng in electrical and computer engineering with specialization in the analysis and control of dynamical systems, and wireless communications, from Jacobs University, Bremen, Germany, in 2016. He worked as a research intern in the area of non-bipedal robot stability for the New Mexico Institute of Mining and Technology, Socorro, United States, and in the area of optimal scheduling of power system generation for the Fraunhofer Institute for Factory Operation and Automation, Magdeburg, Germany. During his undergraduate studies, he was simultaneously part of the ConDyNet project under Prof. Stefan Kettemann's supervision and the Wireless Communications, Signal Processing Research and Development group, under Prof. Giuseppe T. F. Abreu's guidance. His research interests mainly include hybrid approaches in power system analysis, comprising graph theory and electromagnetic transients.  He currently works for decon international GmbH, offering engineering consultancy services for international projects in power transmission, distribution, and renewable energy systems.
%\end{IEEEbiography}


%\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{ab_photo}}]{Giuseppe T. F. de Abreu} (S '99, M '04, SM '09) received the B.Eng. degree in electrical engineering and a specialization (Latu Sensu) degree in telecommunications engineering from the Universidade Federal da Bahia (UFBa), Salvador, Brazil, in 1996 and 1997, respectively, and the M.Eng. degree in physics and the Ph.D. degree in electrical and computer engineering from the Yokohama National University, Yokohama, Japan, in 2001 and 2004, respectively. He was a Postdoctoral Fellow and Adjunct Professor (Docent) of statistical signal processing and communications theory with the Department of Electrical and Information Engineering, University of Oulu, Oulu, Finland, from 2004 to 2006 and 2006 to 2011, respectively. Since 2011, he has been a Professor of electrical engineering with Jacobs University Bremen, Bremen, Germany, and since 2015, he has been a Professor of electrical engineering with Ritsumeikan University, Kyoto, Japan. His research interests include communications and signal processing, such as communications theory, estimation theory, statistical modelling, wireless localization, cognitive radio, wireless security, MIMO systems, energy harvesting networks, random networks, connected vehicles network, and many other topics. He served as an Associate Editor of the IEEE Transactions on Wireless Communications from 2009 to 2014 and currently serves as an Associate Editor of the IEEE Transactions on Communications. He coauthored articles short-listed for the Best Paper Awards at the Asilomar Conference on Signals Systems and Computers in 2009, 2012, and 2014, and was the corecipient of the Best Paper Award at the Workshop on Positioning Navigation and Communications in 2012 and 2014. He was also the recipient of the prestigious JSPS, Heiwa Nakajima and NICT Fellowships in 2010, 2013, and 2015, respectively and the Uenohara Award by Tokyo University in 2000 for his Graduate work.
%\end{IEEEbiography}


%\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{sk_photo}}]{Stefan Kettemann} received his PhD degree in physics at Max Planck Institute for Solid State Physics from the University of Stuttgart, Germany, in 1994. After Postdoctoral positions at the Weizmann Institute in Rehovot, Israel, and the Max Planck Institute for Complex Systems in Dresden, Germany, he became a lecturer at the University of Hamburg, Germany, in 2001. In 2008, he became a lecturer at Jacobs University, Bremen, Germany. In 2009, he was appointed Associate Professor at Postech, Pohang, South Korea. He is currently a Professor of physics at Jacobs University. His  research has concentrated on the application of spectral analysis methods to  disordered and complex systems, including doped semiconductors, quantum networks, spintronic materials and electric power grids. He is a member of the German Physical Society. 
%\end{IEEEbiography}
conclusion.tex
\section{Conclusion}

We  analyze disturbances in three different 
topologies  by mapping them into the problem of  
  diagonalizing a generalized Laplacian matrix,
   solving then for their  eigenmodes and eigenfrequencies. 
      We find that the density of eigenfrequencies
        depends strongly on the topology. For the Random network, it resembles the Mar\u{c}henko-Pastur distribution,
          expected for  uncorrelated random matrices. For the Small-world network and the 
           German transmission grid, we find, however, a 
            strongly peaked density of eigenfrequencies  with long tails. We identify that the 
      long-time transient behavior of perturbations is 
      determined by the generalized algebraic connectivity of the coupling matrix. We find that for the Random network, the intensity 
       of the generalized Fiedler vector is strongly localized.
          We  also find strong indications  that  the degree of localization tends to increase with 
          a decrease of the network's average clustering coefficient $C$. 
 The observation 
      that the eigenmode is strongly localized in the Random network 
      may have important consequences for the design of stable electricity grids: 
                       if  the phase perturbation  is initially in  a  state localized 
                       around a node $r_0$, 
                        with localization length $\xi_k$, then that disturbance  remains 
                         localized there and it 
                        decays exponentially  in time \cite{kettemann}. Thus, lesser meshed grids  may help to localize 
           disturbances more strongly.  Finally, in order to analyze the degree of localization 
           and its connection to the relaxation rate more quantitatively 
           as function of the system topology and its parameters 
            such as the inertia, 
                         we plan to  study this  dependence 
                         systematically in future studies using the inverse participation ratio.  


introduction.tex

\section{Introduction}
 Electric
power grids   provide a highly reliable electrical service
  to billions of customers.
   In fact,  the average outage time experienced by a consumer 
 has kept  decreasing in recent  years, reaching a record low of $12.5$ minutes in Germany, in $2014$ \cite{bna}. However, the energy transition from a centralized power production with unilateral power flow towards an increased supply of decentralized and more volatile  renewable energy resources with bidirectional flow, might become harmful for the stability of electricity grids in the future.
   In  the currently existing grids, the synchronous 
    generators  and synchronous motors
    provide, with their rotating masses, high inertia to the system, which automatically 
     reacts to disturbances \cite{meier,kundur_ieee}.     
     For instance,  an abrupt  increase in  load demand 
    can be momentarily balanced by a change of the kinetic energy of  
     rotating synchronous generators, 
       causing some generators to slow down and deviate from 
        the grid frequency,  but ensuring the overall stability of the network.
       With an increasing share of renewable energy,  this  buffer 
        for the electrical energy  is expected to decrease 
         since solar cells and conventional wind turbines  do not  provide such  inertia to the system \cite{ulbig}. 
          Therefore, it will be increasingly important to  obtain a deeper 
           understanding of
            how fast disturbances decay and  spread  in the grid and how this 
              depends on the topological connectivity and the system parameters, 
              in order 
               to maintain a reliable control of the network.
              

Many authors have studied the role of system topology for the robustness of power grids against large disturbances, such as intentional and random removals of nodes and edges \cite{cuadra, timme,jung}. For small disturbances, on the other hand, the small-disturbance rotor angle stability has been properly defined \cite{kundur_ieee} and thoroughly studied by assessing the solution of the system swing equations and its conditions of stability. In fact, extensive attention has been given to the eigenvalues and eigenvectors of the 
 stability matrix of multiple-machine systems \cite{timme,machowski, coletta} to, for instance, optimize the parameters and grid location of Power System Stabilizers (PSS) \cite{machowski}. Nonetheless, little attention has been given to the 
  propagation  of small disturbances and 
   how the latter depends on   grid topology and the distribution of  system parameters. 
   In order to study the decay and propagation 
    of disturbances, we  implement a hybrid approach  to combine graph theory tools with electric parameters of inductive grids \cite{cuadra} and consider only undirected graphs to depict the smart grid concept, according to which consumers could rapidly become producers and exchange the existing hierarchical power transmission into a bidirectional system.
    
  We specifically study small perturbations in unregulated AC electricity grids by adopting a Synchronous Motor load model. Moreover, since the perturbations are assumed to be small so that the system remains close to its steady state, we linearize the equations in the perturbations. After solving the resulting discrete wave equation on 
   the grid, we find its eigenmodes and eigenfrequencies from a complex Fourier series representation. Then, we study different topologies and randomize the power distribution to analyze how the network arrangement and system paramenters influence the behavior of these modes. Finally, we introduce some potential applications of our research to power system planning and grid control.


 
  
%  %... in which the rotor angle dynamics of generators and loads is described by the same swing equation \cite{meier, motter}. 
%The article is split/formed into sections.. The first section...
mainfinal2610.tex

\documentclass[journal]{IEEEtran}
%\usepackage{blindtext}
\usepackage{graphicx}
%\usepackage{epsfig} 
\usepackage{amsmath}
%\usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
\usepackage[caption=false,font=normalsize]{subfig}
\usepackage{float}

\newcommand{\valuegiv}{0.2675} 



\begin{document}

\title{Analysis of the Dynamics and Topology Dependencies of Small Perturbations in Electric Transmission Grids}


\author{{Luiscarlos~A.~Torres-S\'anchez, Giuseppe T. Freitas de Abreu,~\IEEEmembership{Senior Member~IEEE}~and~Stefan Kettemann} %<-this % stops a space
\thanks{Luiscarlos A. Torres-S\'anchez is with the Department of Energy Division, decon international GmbH, 61352 Bad Homburg, Germany (e-mail: l.torres@decon.de).}% <-this % stops a space
\thanks{G. T. Freitas de Abreu is with the Department of Computer Science and Electrical Engineering, Jacobs University
Bremen, 28759 Bremen, Germany (e-mail: g.abreu@jacobs-university.de), and also with the Department of Electrical
and Electronic Engineering, Ritsumeikan University, Kusatsu 525-8577, Japan (e-mail: g-abreu@fc.ritsumei.ac.jp).}% <-this % stops a space
\thanks{Stefan Kettemann is with the Department of Physics and Earth Science, Jacobs University
Bremen, 28759 Bremen, Germany, and also with the Division of Advanced
Materials Science, Pohang University of Science and Technology, Pohang 790-784, South Korea (e-mail: s.kettemann@jacobs-university.de). }}


%\markboth{Journal of \LaTeX\ Class Files,~Vol.~6, No.~1, January~2007}%
%{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for Journals}





% If you want to put a publisher's ID mark on the page you can do it like
% this:
%\IEEEpubid{0000--0000/00\$00.00~\copyright~2007 IEEE}
% Remember, if you use this you must call \IEEEpubidadjcol in the second
% column for its text to clear the IEEEpubid mark.



% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}




% make the title area
\maketitle


\begin{abstract}
 
We study the dynamics of disturbances in unregulated electric transmission grids through a Synchronous Motor load model. We start our analysis with linearized system equations, which we represent as complex Fourier series to find their eigenmodes and eigenfrequencies. This reduces the problem to the diagonalization of a finite dimensional matrix, which depends on the stationary phase solutions of the grid, and is thereby inherently conditioned by the topology and the power distribution. This matrix is found to belong to the generalized Laplacian matrices, which relates the analysis of perturbation dynamics to a graph theory problem. We consider three networks: Small-world, Random and German transmission grid. We find that the density of eigenfrequencies highly depends on the topology. For the Random network, it resembles the Mar\u{c}henko-Pastur distribution. For the Small-world network and the German transmission grid, we find strongly peaked densities with long tails. Moreover, we find that the algebraic connectivity of the generalized Laplacian matrix determines the long-time transient behavior of perturbations, whereas the intensity of its eigenvector, the generalized Fiedler vector, discloses information of the modes localization. Finally, we find strong indications that the degree of localization tends to increase with
 a decrease of  the network's average clustering coefficient.
\end{abstract}

% Note that keywords are not normally used for peerreview papers.
%\begin{IEEEkeywords}
%Synchronous machines, Synchronous generators, Syncrhonous motors, Graph theory, Spectral analysis, Power system simulation, Topology, Eigenvalues and eigenfunctions, EHV Transmission, Electromagnetic transient propagation, Frequency stability, Power system dynamic stability, Power system control, Power system transient, Power system transient stability, Power transmission control, Stability


%Small Perturbations, Synchronous Motor Model, Generalized Laplacian Matrices, Eigenmode, Eigenfrequency, Localization, Fiedler vector, algebraic connectivity, Small-World Network, Random Network, German Grid, Marchenko-Pastur Distribution.
%\end{IEEEkeywords}


% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
%\IEEEpeerreviewmaketitle

\input{introduction.tex}

\input{phasedyn.tex}

\input{transgrids.tex} %Path Length, Clustering Coefficient 
					%Small-world


         
\input{spectralanalysis.tex}   %Random
							%German
  
    
\input{simulation.tex}

\input{plancontrol.tex}

\input{conclusion.tex}

\section{Acknowledgment}


  We gratefully acknowledge the support of BMBF
   CoNDyNet FK. 03SF0472A. 

%\pagebreak
\input{bibliography.tex}

 
                
   
    
    
       
% needed in second column of first page if using \IEEEpubid
%\IEEEpubidadjcol

% An example of a floating figure using the graphicx package.
% Note that label must occur AFTER (or within) \caption.
% For figures, \caption should occur after the \includegraphics.
% Note that IEEEtran v1.7 and later has special internal code that
% is designed to preserve the operation of \label within \caption
% even when the captionsoff option is in effect. However, because
% of issues like this, it may be the safest practice to put all your
% \label just after \caption rather than within \caption{}.
%
% Reminder: the "draftcls" or "draftclsnofoot", not "draft", class
% option should be used if it is desired that the figures are to be
% displayed while in draft mode.
%
%\begin{figure}[!t]
%\centering
%\includegraphics[width=2.5in]{myfigure}
% where an .eps filename suffix will be assumed under latex, 
% and a .pdf suffix will be assumed for pdflatex; or what has been declared
% via \DeclareGraphicsExtensions.
%\caption{Simulation Results}
%\label{fig_sim}
%\end{figure}

% Note that IEEE typically puts floats only at the top, even when this
% results in a large percentage of a column being occupied by floats.


% An example of a double column floating figure using two subfigures.
% (The subfig.sty package must be loaded for this to work.)
% The subfigure \label commands are set within each subfloat command, the
% \label for the overall figure must come after \caption.
% \hfil must be used as a separator to get equal spacing.
% The subfigure.sty package works much the same way, except \subfigure is
% used instead of \subfloat.
%
%\begin{figure*}[!t]
%\centerline{\subfloat[Case I]\includegraphics[width=2.5in]{subfigcase1}%
%\label{fig_first_case}}
%\hfil
%\subfloat[Case II]{\includegraphics[width=2.5in]{subfigcase2}%
%\label{fig_second_case}}}
%\caption{Simulation results}
%\label{fig_sim}
%\end{figure*}
%
% Note that often IEEE papers with subfigures do not employ subfigure
% captions (using the optional argument to \subfloat), but instead will
% reference/describe all of them (a), (b), etc., within the main caption.


% An example of a floating table. Note that, for IEEE style tables, the 
% \caption command should come BEFORE the table. Table text will default to
% \footnotesize as IEEE normally uses this smaller font for tables.
% The \label must come after \caption as always.
%
%\begin{table}[!t]
%% increase table row spacing, adjust to taste
%\renewcommand{\arraystretch}{1.3}
% if using array.sty, it might be a good idea to tweak the value of
% \extrarowheight as needed to properly center the text within the cells
%\caption{An Example of a Table}
%\label{table_example}
%\centering
%% Some packages, such as MDW tools, offer better commands for making tables
%% than the plain LaTeX2e tabular which is used here.
%\begin{tabular}{|c||c|}
%\hline
%One & Two\\
%\hline
%Three & Four\\
%\hline
%\end{tabular}
%\end{table}


% Note that IEEE does not put floats in the very first column - or typically
% anywhere on the first page for that matter. Also, in-text middle ("here")
% positioning is not used. Most IEEE journals use top floats exclusively.
% Note that, LaTeX2e, unlike IEEE journals, places footnotes above bottom
% floats. This can be corrected via the \fnbelowfloat command of the
% stfloats package.










% if have a single appendix:
%\appendix[Proof of the Zonklar Equations]
% or
%\appendix  % for no appendix heading
% do not use \section anymore after \appendix, only \section*
% is possibly needed

% use appendices with more than one appendix
% then use \section to start each appendix
% you must declare a \section before using any
% \subsection or using \label (\appendices by itself
% starts a section numbered zero.)
%


%\appendices
%\section{Proof of the First Zonklar Equation}
%Some text for the appendix.

% use section* for acknowledgement



% Can use something like this to put references on a page
% by themselves when using endfloat and the captionsoff option.
%\ifCLASSOPTIONcaptionsoff
 % \newpage
%\fi



% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section

% can use a bibliography generated by BibTeX as a .bbl file
% BibTeX documentation can be easily obtained at:
% http://www.ctan.org/tex-archive/biblio/bibtex/contrib/doc/
% The IEEEtran BibTeX style support page is at:
% http://www.michaelshell.org/tex/ieeetran/bibtex/
%\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
%\bibliography{IEEEabrv,../bib/paper}
%
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)




\end{document}



phasedyn.tex

\section{Phase Dynamics Analysis}

\subsection{Mathematical Model}

Phase  dynamics in AC electricity grids have been modeled by active power balance equations with  additional  terms  describing the dynamics of rotating 
    machines \cite{timme, hill,filatrella,schmietendorf,heuck,menck}. We specifically assume loads to be synchronous motors  whose rotor angle $\varphi_i$ dynamics can be modeled by the swing equation of synchronous generators \cite{motter}. This second-order differential equation describes the inertia to changes in kinetic energy through the moment of inertia $J$, and the machine damping through the damping coefficient $\gamma$. Adding these terms to the active power balance equations yields, for purely inductive transmission lines \cite{timme, hill,filatrella,heuck,menck},
    %
    \begin{equation} \label{dynamic0}
        P_i =\left( \frac{J}{2}  \frac{d}{dt}  + \gamma \right)  \left(\frac{d \varphi_i}{d t} \right)^2  + \sum_j K_{ij} \sin 	(\varphi_i-\varphi_j),
    \end{equation}
%
where $K_{ij}$ is the power line capacity. Eq.(\ref{dynamic0}) is analog to an unregulated generator, where primary and secondary 
 frequency control and  voltage regulator actions are disregarded or considered to have large time constants. Therefore, we can assume the mechanical power $P_i$ and excitation voltages to be constant in time \cite{machowski}. Eq.(\ref{dynamic0}), with $K_{ij}=K A_{ij}$, where $A$ is the Adjacency matrix of the network,  corresponds to a homogeneous system, in which all generators, motors, and transmission lines, have identical inertia, damping and power line capacity parameters. 
 
 A homogeneous system allows to focus on the influence of network topology \cite{motter} and can also be used to model a grid with low level of inertia (produced by the integration of renewable energy sources), whose frequency deviations are controlled by a simplistic consideration of fast primary control (e.g. Battery Energy Storage Systems). This control response can  be simply modeled as an additional damping term \cite{ulbig}. Here, we consider fixed voltages (i.e at $V_i= 1$ p.u), which eliminates dynamic terms in the reactive power balance equation
     as they only appear in higher order
      when voltage dynamics- in addition to phase dynamics- 
        are considered \cite{machowski,schmietendorf, kundur_book}.
        
The rotor angle is expressed as $\varphi_i(t) =  \omega   t + \theta_i(t)$, where $\omega$ is the grid frequency and $\theta_i(t)$ is a phase shift. By assuming that $\dot{\theta_i}<<\omega$ and that the rate at which energy is stored in the kinetic term is much less than the rate at which energy is dissipated by friction (i.e  $|J\ddot{\theta}_i|<<2\gamma\omega$), Eq.(\ref{dynamic0}) can be simplified as \cite{filatrella}: 
       
        \begin{equation} \label{dynamic}
        P_i =J \omega \ddot{\theta_i}+2\gamma\omega\dot{\theta_i}+ K\sum_j A_{ij} \sin{(\theta_i-\theta_j)}.
       \end{equation}
       
       
      \subsection{Dynamics of Disturbances in the Grid:}
      
      
 In order to study the propagation of disturbances, we set  $ \varphi_i(t) =  \omega   t + \theta^0_i +  \alpha_i(t)$ 
 with steady state phases  $\theta^0_i$, the
solutions of Eq.(\ref{dynamic}). 
   The dynamics of   the phase disturbances   $\alpha_i(t)$ are governed by: 
   %
        \begin{equation} \label{alphaharmonic3}
        \begin{split}
             \partial_t^2   \alpha_i  + 
                  2 \Gamma  \partial_t   \alpha_i   =& \frac{P_{i}}{J \omega}\\ &- \sum_j \frac{K}{J \omega} A_{ij} \sin   (\theta^0_i
       - \theta^0_j + \alpha_i -\alpha_j),
       \end{split}
          \end{equation}       
          %
          where $\Gamma = \gamma/J$ is the local relaxation rate. Since the steady-state natural or inherent stability of a system can be analyzed via a linearized unregulated condition \cite{machowski}, we consider  small  perturbations from  the  stationary state, as typically done in small-signal stability analyses, and expand  Eq.(\ref{alphaharmonic3}) in ($\alpha_i -\alpha_j$),
        which yields  linear equations on the grid \cite{kettemann}:
       % 
    \begin{equation} \label{alphaharmonic4}
    \partial_t^2   \alpha_i  + 2 \Gamma  \partial_t   \alpha_i 
                   = -
       \sum_j   t_{ij}    (\alpha_i -\alpha_j),
          \end{equation}
         %
         with  coupling amplitude $t_{ij} = \frac{K}{J\omega} A_{ij}\cos(\theta^0_i - \theta^0_j)$ \cite {kettemann}. We let $E$ to be the \textit{coupling matrix} whose entries are formed as $E_{ii}=\sum_{j} t_{ij}$, and $E_{ij}=-t_{ij}$. This matrix has been previously identified in synchronization studies of coupled-oscillator networks, with possibly including ohmic losses \cite{myers}, as well as in linear stability studies for purely inductive grids \cite{coletta}, under the name of  \textit{stability matrix}. 
         
          With $N$ as the number of nodes and the perturbation at node $i$ expressed as a complex Fourier series,  $\alpha_i (t)   =  \sum_{k=1}^N 
               c_{ik} \exp {(-j \epsilon_k t)}$, we plug into
               Eq. (\ref{alphaharmonic4}) to obtain the relation:
                %
                 \begin{eqnarray} \label{randommatrix1}
                    (\epsilon_k^2 +  j2\Gamma \epsilon_k)c_{ik}  
                   =      
       \sum_j   t_{ij}    (c_{ik} -c_{jk}),
          \end{eqnarray}               
               %  
or $(\epsilon^2+j2\Gamma\epsilon)I\vec{c}=E\vec{c}$ in matrix form. From the eigenvalue equation $E\vec{c}= \Lambda\vec{c}$, we calculate the eigenmodes $\vec{c}$ and subsequently find the roots of $(\epsilon^2 +  j2\Gamma\epsilon)I=\Lambda$ to retrieve the eigenfrequencies, $\epsilon$. $c_1=0$ and $\epsilon_1=0$  correspond to the stationary solution. From the real symmetry of $E$, we expect real eigenvalues and real eigenmodes. For $\Gamma =0$, we get $\epsilon^2 I= \Lambda$
                   with eigenfrequencies $\epsilon_k^2= \Lambda_{kk}$ and eigenmodes $c_k$. For $\Gamma \neq 0$ we obtain
               the same eigenmodes $c_k$ and
                complex eigenfrequencies $\Omega_k =  -j \Gamma
                \pm j \sqrt{\Gamma^2 - \epsilon_k^2 }$. For $\epsilon_k \geq \Gamma$,  $\Im (\Omega_{k})=-\Gamma$. For $\epsilon_k < \Gamma$,  we obtain $\Im (\Omega_{k^-})<-\Gamma$, which produces the fastest amplitude decay; but we also obtain $\Im (\Omega_{k^+})>-\Gamma$, which produces the slowest amplitude decay, creating long-lasting perturbations. 
               Since this may increase the impact of 
                    disturbances on the  power system stability,
                it  is highly important to 
                find out 
          the conditions   for  such  slow amplitude decays. 
                   
                          
           Eq.(\ref{alphaharmonic3}), together with the linearized dynamics Eq.(\ref{alphaharmonic4}), contains information about other perturbations in the system. For instance, the change in normalized power flow through the transmission line connecting nodes $i$ and $j$ is given by $\delta F_{ij} (t)= \frac{K}{J\omega}[\sin(\theta_i^0- \theta_j^0 +\alpha_i-\alpha_j)-\sin(\theta_i^0- \theta_j^0)]$. From Eq.(\ref{alphaharmonic3}), $-\sum_j \delta F_{ij} (t)=  \partial_t^2   \alpha_i  + 
                   2 \Gamma  \partial_t   \alpha_i$; which results in $\sum_j \delta F_{ij} (t)= \sum_{k=1}^N c_{ik}(\epsilon_k^2+j2\Gamma \epsilon_k) \exp (-j \epsilon_k t)$; or simply: $\sum_j \delta F_{ij} (t)= \sum_{k=1}^N c_{ik} \Lambda_{kk} \exp (-j \epsilon_k t)$. This is the modal representation of the sum of normalized power flow changes among all the lines connected to node $i$ as well as the anti-phase change of normalized mechanical power at the node, i.e $\frac{\delta P_i(t)}{J\omega}=-\sum_j \delta F_{ij}(t)$.
       
       
       % \subsection{Stationary State:} 
       
       Moreover, the stationary state of Eq.(\ref{dynamic}) can be written in matrix form, using the graph oriented incidence matrix, $B$:
         
           \begin{eqnarray}\label{steady_nonlin}
           P= KB\sin(B^T\theta).
           \end{eqnarray} 
         
        
          To solve for the steady state phases, we choose a DC approximation $\left(\text{i.e}~ B^T \theta<< \frac{\pi}{2}\right)$, which reduces the computational time considerably \cite{overbye} when compared to other more accurate methods such as solving the coupled nonlinear swing equations Eq.(\ref{dynamic}), or solving Eq.(\ref{steady_nonlin}) via a root-finding algorithm. The oriented incidence matrix of size $(N,\kappa)$, where $\kappa$ is the number of edges in the grid, is related to the Laplacian matrix $L$, $L=BB^T$, and to the Adjacency matrix $A$ by $L=D-A$; where $D$ is a diagonal matrix with the node-degree values. In a DC approximation, this leads to: $P=KBB^T\theta$ or $P=KL\theta$. Since we only consider fully connected graphs, $L$ has a zero eigenvalue, which makes the matrix non-invertible. Therefore, we solve for $\theta$ by applying the Moore-Penrose pseudoinverse:
         %
         
          \begin{equation}\label{steady_lin}
           \theta= \frac{1}{K}L^+P,
           \end{equation} 
         %
where $L^+=(L^TL)^{-1}L^T$. This approximation is accurate enough as long as $P_i<< d_i K$, where $d_i$ is the node degree. The reason is intuitive: 
           $P_{avg_i}=\frac{P_i}{d_i}$  is the average mechanical power generated or consumed at node $i$, which leaves or enters the node in the form of electric power through the transmission lines that are connected to it. If we consider that the electric power in each of these lines does not deviate much from the average, then $P_{avg_i}\approx K\sin(\theta_i-\theta_j$), and since the condition $|\theta_i-\theta_j|<< \frac{\pi}{2}$ is needed for linearization, then $\frac{P_{avg_i}}{K}<<1$ must be fulfilled. The statement $P_i<< d_i K$ follows.
           
          
plancontrol.tex

\section{Potential Contributions to Grid Planning and Control}
 
Power System Planning: The addition of one node or one edge should lower the network's clustering coefficient in order to improve the grid's robustness against cascading failure, to generate a greater generalized algebraic connectivity of the coupling matrix and a much more confined eigenfrequency spectrum (e.g Fig. \ref{fig:rand2}). This would also allow the system to strongly localize the perturbation eigenmodes in a given area, which would improve the controllability of the disturbance. \\

Power System Stabilizer: The main method used today to guard power systems against small-signal instability is the off-line tuning of PSS. These are local controllers to the generators for which significant disadvantages have been found \cite{derong}. This locality can be overcome by new controllers that use remote signal inputs from the network or adaptively change their setpoints according to the power system operating conditions while bypassing the large analytical calculations needed to obtain the results \cite{derong}. Analyzing disturbance modes from a topology-based approach could provide a default state of the system under disturbances, requiring infrequent tunning of the control devices in real time. 
%%
Let us take, for instance, the PSS whose input parameters are the end-of-schaft speed and the change in generated real power. This PSS calculates an equivalent speed deviation  $\Delta \omega_{eq}$ to then output a leading voltage signal to control the perturbation modes. The $\Delta \omega_{eq}$ requires the installation of a mesuring system to obtain the end-of-schaft speed. It also requires the calculation of two integrals to first obtain the change in mechanical power and then the equivalent speed deviation \cite{machowski}. If, however, we represent the latter in terms of known modes, we could simplify the system into    $\Delta \omega_{eq} (t)= \dot{\alpha_i}(t)$.\\


%%-\Sigma_{k=2}^N j\epsilon_{k}c_{ik}\exp{(-j\epsilon_{k}t)}$. 



%\textbf {Secondary Frequency Control}:
%Authors in [] establish an optimization problem in which secondary controllers should react to power imbalance in the grid in the most economically viable way. In our case, we could consider the use of secondary controllers to optimally react to the global maximum of the change in normalized mechanical power, that is to react to the value $\max\Sigma_i \delta P_i(t)/(J\omega)$ at $t_{max}$ in order to drive the system faster to its steady-state.\\


Dynamic Transmission Topology Control (TC): Although scalable practical solutions have not been yet achieved, there is an ongoing interesest in the research community for this emerging control technology \cite{sari} for its potential to manage the uncertainty of power sources and flow patterns on a grid with high penetration of renewable energy \cite{wang}. The control strategy consists of switching lines on or off to relieve voltage and line flow violations \cite{wang}, with a considerable impact on small-signal and transient stability \cite{sari, wang}. Based on our analysis, any TC action that modifies the structure of the grid should firstly guarantee phase-cohesiveness and secondly that $\sqrt{a(L')}$, where $L'$ is the Laplacian matrix of the modified topology, is never less than the damping rate of the system; otherwise, there will surely exist at least one mode that decays slowly in time. 

simulation.tex

 
\section{Simulation}

\subsection{Selection of Grid Parameters}
 
    We assign values of power to each node from a bipolar distribution, i.e $P_i=\pm P$ in Watts. $P>0$ for generators and $P<0$ for motors (consumers). Eq.(\ref{dynamic}) synchronizes at a frequency $\dot{\theta}_{\text{synch}}= \sum_i^N \frac{P_i}{\gamma}$ \cite{dorfler}, which implies that the condition $\sum_i^N P_i=0$ must be fulfilled at all times for the system to reach steady state. This is a realistic consideration since power generation must constantly match load demand. We take as reference the German installed capacity of $199.2$ GW as per November $10^{th}$, $2015$ \cite{ministry}, and consider half of the nodes-for an even number- to be generators and the remaining half to be consumers. For the $500$-node complex networks (Small-world and Random) this results in $P_i= \pm 796.80$ MW. For the $489$-node German grid, we have on average $P_i\approx \pm 814.72$ MW. We choose the grid angular frequency $\omega= 2\pi$($50$ Hz) and moment of inertia $J=10^5$ kgm$^2$. This $J$ is, for instance, for a generator working  at $\omega$, with inertia constant $H=3 \frac{MJ}{MVA}$, rated at $100$ MVA. This $H$ resembles those of high-speed and slow-speed water-wheel generators and non-condensing turbine generators \cite{birron}. We apply Eq.(\ref{syn_cond}) to the complex networks and retrieve $K_{sworld}\approx 5.24$ GW and $K_{rand} \approx  335.14$ MW. For the German transmission grid, we assign $K_{G}=10$ GW, which is strong enough to keep small angular differences. Finally, we select the damping rate $\Gamma=1$ Hz for all topologies considered.
    %
    \subsection{Density of Eigenfrequencies} 
    
    We study the eigenfrequency density $\rho_m(\epsilon_k)$ for different arrangements of generators and consumers by randomizing, within the bipolar distribution, the $P_i$ of each node. We perform $R=1500$ iterations to obtain the average density, $\rho(\epsilon_k)= \frac{1}{R}\sum_{m=1}^{R}\rho_m(\epsilon_k)$. The results are shown in Figs. (\ref{fig:sw2},\ref{fig:rand2},\ref{fig:ger2}). The stationary solution corresponds to $\epsilon_1= 0$ Hz and it is not shown. We observe that: 
    \begin{itemize}
    \item The nonzero eigenfrequencies for all networks exceed, for the chosen parameters, the damping rate $\Gamma$, $\epsilon_k> \Gamma$, so that disturbances decay exponentially fast with relaxation rate $\Gamma$. However, under the current high integration of renewable energy, the reduction of $J$ \cite{ulbig} changes the power system parameters and may thereby reduce the eigenfrequencies $\epsilon_k$ in comparison to $\Gamma$, leading to slower relaxation rates. 
    \item Although highly distributed, a significant peak of the German grid eigenfrequency density is located around $\sqrt{\frac{K_{G}}{J\omega}} \approx 17.84$ Hz. This could be attributed to the influence of the 149 one-degree nodes in the grid whose diagonal entries in $E$ are, with very small angular differences in the system, $\sqrt{\frac{K_{G}}{J\omega} \cos(\theta_i-\theta_j)} \approx \sqrt{\frac{K_{G}}{J\omega}}$.
    \item The average eigenfrequency density of the Random network resembles the Mar\u{c}henko-Pastur distribution, expected for uncorrelated random matrices. 
    
    \end{itemize}
    
    
    \subsection{The Effect of Clustering}
    
  The values $K_{sworld}$ and $K_{rand}$ were obtained for the vectors $P_{sworld}$ and $P_{rand}$ assigned to Figs. (\ref{fig:sw},\ref{fig:rand}), in which there are visible clusters of generators and consumers. If $P$ is randomized, reducing the size of the clusters, smaller values of $K$ can be found. This effect was studied in \cite{buzna} for a bipolar distribution of frequencies (power in our case), and it was shown that synchronization is enhanced when adjacent nodes have opposite frequencies, resulting in a diminished frequency similarity throughout the grid. This simply means that synchronization is enhanced when generators are surrounded by consumers and vice versa. Moreover, critical effects, such as cascading failures, are less likely to be triggered if a greater frequency dissimilarity prevails. This was demonstrated statistically in \cite {jung}, where authors claim that the existence of large clusters of generators and consumers turns the grid vulnerable against cascading failures, since the likelihood for a whole cluster to disconnect at once appears to increase with increasing cluster size. 
 
 Figs. (\ref{fig:sw1},\ref{fig:rand1}) provide an insight into the effect of randomization. The squares are those machines that have switched their power in comparison to Figs. (\ref{fig:sw},\ref{fig:rand}). It is clear that clusters are reduced, resulting in smaller power line capacities (i.e $K_{sworld} \approx 1.12$ GW and $K_{rand} \approx 309.34$ MW), but also in  smaller nonzero eigenfrequencies. For the Small-world network, we obtained for Fig. (\ref{fig:sw}), $\epsilon_2=5.5470$ Hz, whereas for Fig. (\ref{fig:sw1}), $\epsilon_2=2.7198$ Hz. For the Random network, we obtained for Fig. (\ref{fig:rand}), $\epsilon_2= 5.8986$ Hz, whereas for Fig. (\ref{fig:rand1}), $\epsilon_2= 5.6250$ Hz. In fact, out of the $1500$ iterations for each complex network, no single value of $\epsilon_2$ was greater than the ones from Figs. (\ref{fig:sw},\ref{fig:rand}). 
 This implies that whereas clustering is detrimental to grid stability and to cascading outage prevention, the larger power capacity needed to ensure stability in the presence of clusters results in an increment of the smallest nonzero eigenfrequency, leading in fact to faster mode damping rates and thereby to a greater resilience of the power system to small perturbations.

 
  \subsection{Spatial Distribution of the Eigenmode Intensity}
   In Figs. (\ref{fig:sw3},\ref{fig:rand3},\ref{fig:ger3}), we show the intensity, $|c_{i2}|^2$, of the generalized Fiedler vector; the eigenmode corresponding to the smallest nonzero eigenfrequency. The intensity at each node is 
    divided by $\max_i |c_{i2}|^2$.
    We observe that in the Random network, Fig. (\ref{fig:rand3}),
     the eigenmode is strongly localized with most of its intensity on a single node. 
     In the Small-world network,  Fig.  (\ref{fig:sw3}), the eigenmode intensity is 
   spread over many nodes, which are far away from each other. In the German transmission grid, 
       the intensity is spread over most of the grid, with greater intensity in the Southwestern and Northwestern part of the system, see Fig. (\ref{fig:ger3}). 
       To understand this behavior, we can refer to the fact that the discrete wave equation 
         Eq.(\ref{alphaharmonic4}) was derived first  
         for the problem of randomly coupled atoms
                 in  harmonic approximation and has been studied intensively  
                  for various random distributions of the coupling $t_{ij}$ \cite{dyson,alexander,ziman,wegner}.            
        For nonzero eigenfrequency  $\epsilon_k$, 
               the eigenmodes were found,  for a random chain of nodes,
                to be localized with  localization length  
                $ \xi(\epsilon_k) \sim 1/\epsilon_k$ \cite{dyson,alexander,ziman,wegner,john},
                 due to the random scattering of waves along the chain. This is an example of the 
                  so-called
                  Anderson localization, which  
                  is enhanced when the amplitude of  randomness is increased \cite{anderson}.
                 In grids with higher degree $d$, the localization length is typically found to be larger.                 
                     Moreover, 
                      the localization length is typically smallest in tree-like grids, whereas it becomes 
                       larger the more meshed the grid is; in which case the eigenmodes 
                         can become even delocalized  \cite{wegner, john}. On the other hand, the average clustering coefficient $C$ is a measure of 
                         how strongly meshed a grid is. 
                     We observe that  the Random network, Fig. (\ref{fig:rand3}),
                       has a  very small average clustering coefficient $C(1) =0.016$, 
                       explaining the fact that its eigenmode is strongly localized, 
                        whereas the  German transmission grid, shown in Fig. (\ref{fig:ger3}), is  meshed
                      with  $C= 0.2021$  and  the Small-world network in 
                       Fig.   (\ref{fig:sw3}) is more strongly meshed  with $C = 0.61841$, explaining that the eigenmode intensity 
                        in these grids is more delocalized and spread over  many nodes.     
                                
     
                  
 
 % Fig. SW, gave the epsilon=5.5470, and for the random network epsilon=5.8986. Out of the 1500 iterations, none of them was bigger thant this one..
 
 %On the other hand, Eq.(\ref{syn_cond}) does not necessarily hold for the German grid; thus, we choose a common value for power line capacity, $K=10$ GW [reference].
 
 

 


 
 %\include{mp}

spectralanalysis.tex


  \section{Spectral Analysis}
   
          \subsection{Generalized Laplacian Matrix} 
         
   % \subsection{Sufficient Power Line Capacity} 
   
      Authors in  \cite{dorfler} showed that a sufficient condition for Small-world and Random networks to reach cohesive phases (that is, that all angular distances $|\theta_i- \theta_j|$ are bounded $|\theta_i- \theta_j| \leq \zeta < \frac{\pi}{2}$,
       where the upper bound is known from the power-angle curve of a synchronous generator connected to an infinite busbar), is given by $||B^{T}L^{+}P||_{\infty}\leq K\sin(\zeta)$, where $||\cdot||_{\infty}$ is the infinity-norm. In the limit $\zeta -> \frac{\pi}{2}$, we find:

 \begin{equation} \label{syn_cond}
  ||B^{T}L^{+}P||_{\infty}< K.
 \end{equation}
 
 If we compare Eq.(\ref{steady_lin}) and Eq.(\ref{syn_cond}), we observe that $||B^{T}\theta||_{\infty}< 1$. This implies that $\max_{i,j}|\theta_i-\theta_j|<1$ (approx. $57.30^o$). If Eq.(\ref{syn_cond}) is imposed, then $E$ is a diagonally dominant matrix. A diagonally dominant matrix $W$ must satisfy $|W_{ii}|\geq\Sigma_{j\neq i} |W_{ij}|$. Since $\max_{i,j}|\theta_i-\theta_j|<1$, the equality $|- \Sigma_{j \neq i} E_{ij}|=\Sigma_{j \neq i} |E_{ij}|$ holds, as the left-hand side, which represents how the diagonal of $E$ is formed with respect to the nonzero off-diagonal entries, only contains negative terms for each $E_{ij}$. This makes $E$ a diagonally dominant matrix with positive diagonal entries. Therefore, $E$ is a positive semi-definite matrix, just as the undirected-unweighted graph Laplacian $L$. In conclusion, the coupling matrix $E$ is a positive semi-definite and real-symmetric matrix; thus, its eigenvalues are real and non-negative, indicating that the system always reaches steady-state under phase-cohesiveness. Moreover, $E$ has positive entries along the diagonal, negative entries for adjacent nodes and zeros for nonadjacent nodes; therefore, it can be considered as a generalized Laplacian matrix. These matrices, which include the graph Laplacian and are found in the inverse eigenvalue problem of a graph \cite{hogben}, fulfill precisely the same conditions for the off-diagonal entries, but have no restrictions on the diagonal entries. Beyond that, the coupling matrix can be considered as nothing else than the Laplacian matrix of a positively weighted graph. Then, all properties of $L$ also apply to $E$, including the decomposition into the product of oriented incidence matrices, i.e $E= QQ^T$, where $Q_{ij}=\pm\sqrt{\frac{K}{J\omega}\cos(\theta_i-\theta_j)}.$



%%\begin{equation}
%Q_{ij}=
%\begin{cases} 
 %     -\sqrt{\frac{K}{J\omega}\cos(\theta_i-\theta_j)} & \text{if}~ B_{ij}=-1, \\
  %    \hfil \sqrt{\frac{K}{J\omega}\cos(\theta_i-\theta_j)} & \text{if}~  B_{ij}=1, \\
   % \hfil 0 & \text{otherwise}, \\ 
%\end{cases}
%\end{equation}
%\[
%Q_{ij}=
 %\begin{cases} 
  %    -\sqrt{\frac{K}{J\omega}\cos(\theta_i-\theta_j)} & \text{if}~ B_{ij}=-1, \\
   %   \hfil \sqrt{\frac{K}{J\omega}\cos(\theta_i-\theta_j)} & \text{if}~  B_{ij}=1, \\
    %\hfil 0 & \text{otherwise}, \\ 
   %\end{cases}
%\]

 %Recall that the analysis in this subsection is only valid as long as $B^T\theta<\frac{\pi}{2}$, which is anyhow fulfilled by Eq.(\ref{syn_cond}).

 \subsection{Algebraic Connectivity}
 
 The smallest nonzero eigenvalue of the Laplacian matrix is called the algebraic connectivity, and its corresponding eigenvector, the Fiedler vector \cite{fiedler}. We refer to this algebraic connectivity of a graph $G$ as $a(G)$. 
 Since the coupling matrix depends on the angular differences, which are inherently related to the power distribution, it would be highly convenient to know lower and upper bounds 
 for its respective connectivity to avoid performing an eigenvalue decomposition every time $P$ changes. The algebraic connectivity $a_L(G)$ of the matrix $\frac{K}{J\omega} L$, can be set as the upper bound for the generalized  algebraic connectivity $a_E(G)$ of 
  the  coupling matrix $E$. The smaller the angular differences, the more the coupling matrix approaches the scaled Laplacian $\left(\text{i.e}~ E \rightarrow \frac{K}{J\omega}L \right)$, and the closer $a_E(G)$ gets to $a_L(G)$. The matrix $\frac{K}{J\omega} L$ corresponds to the coupling matrix of a network with $P=0$, where the only power in the grid is that of the perturbation itself. It would also be very useful to provide
  lower bounds for the generalized algebraic connectivity of $E$. Some lower bounds have been derived for weighted graph Laplacian matrices \cite{kirkland, bermann}, but the inclusion of, for instance, the weighted isoperimetric number \cite{bermann}, makes their calculation computational expensive  in comparison to 
   the explicit solution of the  eigenvalue problem, which  we perform in the following section. 
 

         
 
%Recall that the analysis in this subsection is only valid as long as $B^T\theta<\frac{\pi}{2}$, which is anyhow fulfilled by Eq.(\ref{syn_cond}).

 	% 1---L
    % 2---D
    % 3---R
    % 4--- U   
 
 



transgrids.tex
 \section{Electric Power Transmission Grid Models}
        
          
         %\subsection{Small-world and Random Network:} 
         Authors in \cite{watts_strogatz} proposed a model that interpolates between a lattice and a random graph based on a rewiring probability $p$ for each edge. For a certain range of $p$, there is a coexistence of small Average Path Length, $V(p)$, and high Average Clustering Coefficient, $C(p)$, forming the Small-world network, which mimics many real-world networks that contain small average path lengths, but also have unusually large clustering coefficients \cite{barabasi}. The Average Path Length is defined as: $V= \frac{1}{N(N-1)}\sum_{i,j}l_{ij}$, where $l_{ij}$ is the shortest distance between node $i$ and $j$ for an undirected graph, and the sum runs over all $i$-$j$ combinations. The Clustering Coefficient ($CC$) measures how nodes tend to cluster together by analyzing how neighbors (one link separation) of a given node $i$ are connected among themselves. Mathematically, $CC_i= \frac{h_i}{\frac{1}{2}(d_i)(d_i-1)}$, where $h_i$ is the number of links shared among the $d_i$ neighbors of $i$. The Average Clustering Coefficient is the $CC$ averaged over all nodes, $C= \frac{1}{N}\sum_i CC_i$. 
         
         A rewiring probability $p=1$ generates a Random network, which may not be necessarily similar to the Erd\"os-R\' enyi random network commonly referred to in the literature. To be more precise, despite similar clustering coefficients and average path lengths, a Watts-Strogatz network with $p=1$ is not identical to an  Erd\"os-R\' enyi random network with same size and same mean degree $d_{avg}$, since for example, the Watts-Strogatz algorithm does not allow nodes to exist with degree smaller than $\frac{d_{avg}}{2}$, whereas Erd\"os-R\' enyi does \cite{brenton}.
         
         Here, we study Small-world and Random networks. Firstly, there is a considerable amount of transmission grids that present similar characteristics to the former: Sweden, Finland, Norway, part of Denmark, U.S Western States, Shanghai, Italy, France, Spain \cite{cuadra} and Northern China \cite{hu}. Secondly, Small-world networks have  economical and structural feasible features for electricity distribution in smart grids, as proven by using real data from the Dutch power grid \cite{pagani}. On the other hand, Random networks are proven to be more robust than multiple networks against intentional attacks \cite{cuadra}, which makes their inclusion also important for our study.

     % 1---Left
    % 2---Down
    % 3---Right
    % 4--- Up
    
    
     \begin{figure}[h]
   \includegraphics[totalheight=\valuegiv\textheight, width=\columnwidth]{apl_acc.eps}
   \caption{Normalized $V(p)$ and $C (p)$, with $V(0)=25.40$ and $C(0)\approx 0.67$,  for the Watts-Strogatz model, with $N=500$ and $d_{avg}=10$. Twenty-five realizations per $p$.}  
   \label{fig:apl_acc}
    \end{figure}
         
         To create our grids, we select $N=500$ and $d_{avg}=10$. For the Small-world network we fulfill the condition $N>> d_{avg} >> \ln (N) >> 1$, to have a sparse but connected graph \cite{barabasi}. Then, we find $V(p)$ and $C(p)$, for different values of $p$. For each $p$, we average $V(p)$ and $C(p)$ over 25 realizations. We compare both normalized parameters in Fig. (\ref{fig:apl_acc}), and then retrieve the $p$ that yields the greatest difference between them. This rewiring probability is $p=3.42*10^{-2}$, with normalized parameters $V(p)/V(0)=0.198$ and $C(p)/C(0)=0.923$, which results in the Small-world network shown in Fig. (\ref{fig:sw}).  For the Random network, we set $p=1$, which results in the grid shown in Fig. (\ref{fig:rand}), with $V(p)/V(0)=0.1154$ and $C(p)/C(0)=0.0239$. The parameters for $p=0$ are $V(0)=25.40$ and $C(0)\approx 0.67$.
         

% 1---Left
    % 2---Down
    % 3---Right
    % 4--- Up
    
    
   \begin{figure*}[htp]
 \centering
 \subfloat[Generators are represented as circles and consumers as crosses. Power line capacity of the network, $K_{sworld}\approx 5.24$ GW. Smallest nonzero eigenfrequency, $\epsilon_2= 5.5470$ Hz.]{\includegraphics[trim = 28.2mm 13.5mm 18mm 22mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{small_world.eps}\label{fig:sw}}\hfill
   \subfloat[Randomization of generator and consumer clusters. The squares are those machines that have switched to the opposite power in comparison to Fig. (\ref{fig:sw}). Power line capacity of the network, $K_{sworld} \approx 1.12$ GW. Smallest nonzero eigenfrequency, $\epsilon_2 \approx 2.7198$ Hz.  ]{\includegraphics[trim = 28.2mm 13.5mm 18.2mm 22mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{small_world1.eps}\label{fig:sw1}}\hfill
  \subfloat[Bars size, $0.125$ Hz. Smallest value of $\epsilon_2$ found after $1500$ iterations, $\epsilon_2=2.3764$ Hz.]{\includegraphics[trim = 0mm 0mm 0mm 0mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{small_world2.eps}\label{fig:sw2}} \hfill
   \subfloat[Intensity of the generalized Fiedler vector components (i.e $|c_{i2}|^2$) for the grid in Fig. (\ref{fig:sw1}).]{\includegraphics[trim = 10.5mm 13mm 18.2mm 19mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth, clip, totalheight=0.30\textheight, width=\columnwidth]{small_world3.eps}\label{fig:sw3}} \hfill
  \caption{Small-world Network from the Watts-Strogatz Model with Parameters: $N=500$, $d_{avg}=10$, $p=3.42*10^{-2}$, $V(p)/V(0)=0.198$, $C(p)/C(0)=0.923$, $V(0)=25.40$ and $C(0)\approx 0.67$.}
\end{figure*}


 % 1---Left
    % 2---Down
    % 3---Right
    % 4--- Up
         
 \begin{figure*}[!t]
 \centering
 \subfloat[Generators are represented as circles and consumers as crosses. Power line capacity of the network, $K_{rand}\approx 335.14$ MW. Smallest nonzero eigenfrequency, $\epsilon_2= 5.8986$ Hz.]{\includegraphics[trim = 35mm 15mm 15mm 16mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{random.eps}\label{fig:rand}}\hfill
   \subfloat[Randomization of generator and consumer clusters. The squares are those machines that have switched to the opposite power in comparison to Fig. (\ref{fig:rand}). Power line capacity of the network, $K_{rand} \approx 309.34$ MW. Smallest nonzero eigenfrequency, $\epsilon_2 \approx 5.6250$ Hz.]{\includegraphics[trim = 35mm 15mm 15mm 16mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{random1.eps}\label{fig:rand1}}\hfill
     \subfloat[Bars size, $0.125$ Hz. Smallest value of $\epsilon_2$ found after $1500$ iterations, $\epsilon_2=5.4880$ Hz.]{\includegraphics[trim = 0mm 0mm 0mm 0mm, clip,totalheight=\valuegiv\textheight, width=\columnwidth]{random2.eps}\label{fig:rand2}} \hfill
   \subfloat[Intensity of the generalized Fiedler vector components (i.e $|c_{i2}|^2$) for the grid in Fig. (\ref{fig:rand1}).]{\includegraphics[trim = 13mm 13mm 15mm 12mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{random3.eps}\label{fig:rand3}} \hfill
  \caption{Random Network from the Watts-Strogatz Model with Parameters: $N=500$, $d_{avg}=10$, $p=1$, $V(p)/V(0)=0.1154$, $C(p)/C(0)=0.0239$, $V(0)=25.40$ and $C(0)\approx 0.67$.}
\end{figure*}

    

      % \subsection {German Transmission Grid:}
        %The data for the Extra-high-AC Voltage (380 kV and 220 kV) fully connected graph can be found in \cite{scigrid}. The grid, consisting of 489 nodes, is shown in Fig. (\ref{fig:ger}), with parameters $V=9.9384$ and $C=0.2021$.
        
        We also consider the fully connected graph of the Extra-high-AC Voltage (380 kV and 220 kV) German transmission grid, which can be found in \cite{scigrid}. The grid, consisting of 489 nodes, is shown in Fig. (\ref{fig:ger}), with parameters $V=9.9384$ and $C=0.2021$.
       
 \begin{figure*}[!t]
 \centering
 \subfloat[Generators are represented as circles and consumers as crosses. Power line capacity of the network, $K_{G}= 10$ GW.]{\includegraphics[trim = 25mm 14mm 25mm 12mm, clip, totalheight=0.36\textheight, width=\columnwidth]{german.eps}\label{fig:ger}}\hfill
 \subfloat[Intensity of the  generalized Fiedler vector components (i.e $|c_{i2}|^2$) for the grid in Fig.(\ref{fig:ger}).]{\includegraphics[trim =25mm 14mm 12mm 12mm, clip, totalheight=0.36\textheight, width=\columnwidth]{german3.eps}\label{fig:ger3}}\hfill
     \subfloat[Bars size, $0.5$ Hz. Smallest value of $\epsilon_2$ found after $1500$ iterations, $\epsilon_2=1.7343$ Hz. Peak density around $\rho_{peak}= \sqrt{\frac{K_G}{J\omega}}= 17.84$ Hz.]{\includegraphics[trim = 0mm 0mm 0mm 0mm, clip, totalheight=\valuegiv\textheight, width=\columnwidth]{german2.eps}\label{fig:ger2}} \hfill
  \caption{Extra-high-AC Voltage (380 kV and 220 kV) German Transmission Grid with Parameters: $N=489$, $d_{avg}=2.71$, $V=9.9384$, $C=0.2021$.}
\end{figure*}
       
      % {\includegraphics[trim = 25mm 14mm 25mm 12mm, clip, totalheight=0.375\textheight, width=\columnwidth]{german.eps}\label{fig:ger}}\hfill

 % 1---Left
    % 2---Down
    % 3---Right
    % 4--- Up
         
1706.10130.extracted_bib
{bna Versorgungsqualit\"at-SAIDI-Werte 2006-2014, Bundesnetzagentur, 2015.
1706.10130.extracted_bib
{meier A. von Meier, Electric Power Systems: A Conceptual Introduction. John Wiley Sons, Inc, 2006.
1706.10130.extracted_bib
{kundur_ieee P. Kundur, J. Paserba and V. Ajjarapu, ``Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions'', {IEEE Transactions on Power Systems, vol. 19, no. 3, pp. 1387 - 1401, Aug. 2004.
1706.10130.extracted_bib
{ulbig A. Ulbig, T. Borsche and G. Andersson, arXiv:1312.6435 {[math.OC] (2014).
1706.10130.extracted_bib
{cuadra L. Cuadra, S. Salcedo-Sanz, J. del Ser and Z. Woo Geem. A critical review of robustness in power grids using complex networks concepts. Energies {8, 9211-9265 (2015).
1706.10130.extracted_bib
{timme % Self-Organized Synchronization in Decentralized Power Grids, M. Rohden, A. Sorge, M. Timme and D. Witthaut, Phys. Rev. Lett. {109, 064101 (2012).
1706.10130.extracted_bib
{jung M. Rohden, D. Jung, S. Tamrakar and S. Kettemann, Cascading failures in ac electricity grids, Phys. Rev. {E 94, 032209 (2016). %6-10
1706.10130.extracted_bib
{machowski J. Machowski, J. Bialek, J. Bumby, Power System Dynamics: Stability and Control, Wiley, 2008.
1706.10130.extracted_bib
{coletta T. Coletta, P. Jacquod, Linear Stability and the Braess Paradox in Coupled-Oscillator Networks and Electric Power Grids, Phys Rev E {93, 032222 (2016).
1706.10130.extracted_bib
{hill A. Bergen and D. Hill, ``A structure preserving model for power system stability analysis'', {IEEE Trans. on Power Apparatus and Systems, vol. PAS-100, no.1, pp. 25-35, Jan. 1981.
1706.10130.extracted_bib
{filatrella G. Filatrella, A. Nielsen, and N. Pedersen, Eur Phys. J B {61, 485 (2008). %11-15
1706.10130.extracted_bib
{schmietendorf K. Schmietendorf, J. Peinke, R. Friedrich and O. Kamps, Eur. Phys. J. Spec. Top. {223, 2577 (2014).
1706.10130.extracted_bib
{heuck K. Heuck, K.-D. Dettmann and D. Schulz, {Elekrtische Energieversorgung, Springer, Wiesbaden (2013).
1706.10130.extracted_bib
{menck P. J. Menck, J. Heitzig, J. Kurths and H. Schellnhuber, %How dead ends undermine power grid stability Nature Communications, 5, 3969 (2014).
1706.10130.extracted_bib
{motter T. Nishikawa, A. Motter, New J. Phys. {17, 015012 (2015).
1706.10130.extracted_bib
{kundur_book P. Kundur, Power System Stability and Control, Mc Graw Hill (1994). %16-20
1706.10130.extracted_bib
{kettemann S. Kettemann, Delocalization of disturbances and the stability of AC electricity grids, Phys. Rev. E {94, 062311 (2016).
1706.10130.extracted_bib
{myers A. Motter, S. Myers, M. Anghel and Takashi Nishikawa, Nature Physics {9, 191-197 (2013).
1706.10130.extracted_bib
{overbye T. Overbye, X. Cheng, Y. Sun, ``A comparison of the AC and DC power flow models for LMP Calculations'', {Proceedings of the 37th Annual Hawaii International Conference on System Sciences, 2004.
1706.10130.extracted_bib
{watts_strogatz D. Watts and S. Strogatz, Collective dynamics of small-world networks, Nature {393, 440-442 (1998).
1706.10130.extracted_bib
{barabasi R. Albert and A. Barab\'asi, Statistical mechanics of complex networks, Reviews of Modern Physics, vol. 74, no.1, 47-97 (2002). %21-25
1706.10130.extracted_bib
{brenton B. Prettejohn, M. Berryman and M. McDonnell, Methods for generating complex networks with selected structural properties for simulations: a review and tutorial for neuroscientists, Frontiers in Computational Neuroscience, vol. 5, 5-11 (2011).
1706.10130.extracted_bib
{hu L. Fu, W. Huang, S. Xiao, Y. Li and S. Guo. Vulnerability assessment for power grid based on small-world topological model. {Pacific Power and Energy Engineering Conference (APPEEC), 2010 Asia-Pacific, 2010.
1706.10130.extracted_bib
{pagani G. Pagani and M. Aiello. Power grid complex network evolutions for the smart grid. Phys. A Stat. Mech. Its Appl. , 396, 248-266, 2014.
1706.10130.extracted_bib
{scigrid W. Medjroubi, C. Matke and D. Kleinhans, SciGRID an open source reference model for the European transmission network, version 0.1, http://www.scigrid.de/ (2015).
1706.10130.extracted_bib
{dorfler F. Dorfler, M. Chertkov and F. Bullo. Synchronization in complex oscillator networks and smart grids, {Proceedings of the National Academy of Sciences of the United States of America, vol. 110, no. 6, pp. 2005-2010 (2013). %26-30
1706.10130.extracted_bib
{hogben L. Hogben, Spectral graph theory and the inverse eigenvalue problem of a graph, Electronic Journal of Linear Algebra, ISSN 1081-3810, vol. 14, pp. 12-31, Jan. 2005.
1706.10130.extracted_bib
{fiedler M. Fiedler, Algebraic connectivity of graphs. Czechoslovak Math. J. 23 (1973), 298-305. %31
1706.10130.extracted_bib
{kirkland S. Kirkland, I. Rocha and V. Trevisan, Algebraic connectivity of k-connected graphs, Czechoslovak Mathematical Journal, vol.65, no.1, pp. 219-236, Apr. 2015.
1706.10130.extracted_bib
{bermann A. Bermann, X-D. Zhang, Lower bounds for the eigenvalues of Laplacian matrices, Linear Algebra and its Applications, vol. 316, 1-3, pp. 13-20, Sept. 2000.
1706.10130.extracted_bib
{ministry German Federal Ministry for Economic Affairs and Energy. Available in https://www.bmwi.de/Redaktion/EN/Artikel/Energy/strommarkt-der-zukunft-zahlen-und-fakten.html (2017).
1706.10130.extracted_bib
{birron D. Kothari and I. Nagrath, Power System Engineering 2nd edition, The McGraw Hill Companies (2008).
1706.10130.extracted_bib
{buzna L. Buzna, S.Lozano and A. D\'iaz-Guilera. Synchronization in symmetric bipolar population network, Physical Review E {80, 066120, 2009.
1706.10130.extracted_bib
{dyson F. J. Dyson, Phys. Rev. {92, 1331 (1953).
1706.10130.extracted_bib
{alexander S. Alexander, J. Bernasconi, W. R. Schneider, R. Orbach, Rev. Mod. Phys. {53, 175 (1981).
1706.10130.extracted_bib
{ziman T. A. L. Ziman, Phys. Rev. Lett. {49, 337 (1982).
1706.10130.extracted_bib
{wegner F. J. Wegner in E. Abrahams, 50 years of Anderson Localization, World Scientific (2010).
1706.10130.extracted_bib
{john S. John, H. Sompolinsky, M.J. Stephen, Phys. Rev. B {27, 5592 (1983).
1706.10130.extracted_bib
{anderson P. W. Anderson, Phys. Rev. {109, 1492 (1958).
1706.10130.extracted_bib
{derong D. Liu, P. Antsaklis, Stability and Control of Dynamical Systems with Application, p. {328, Springer, 2003.
1706.10130.extracted_bib
{sari A. T. Sari\' c, A. M. Stankovi\' c, "Rapid small-signal stability assessment and enhancement following changes in topology", {IEEE Transactions on Power Systems, vol. 30, no. 3, pp. 1155–1163, May, 2015.
1706.10130.extracted_bib
{wang W. Wang, G. M. Huang, Impacts of smart grid topology control on power system stability with renewable integration, {North American Power Symposium (NAPS), 2016. %
1706.10130.extracted_bib
{verdu %Verdu