Document [original]

Chaos 28, 101102 (2018); https://doi.org/10.1063/1.5054850 28, 101102

Qualitative stability and synchronicity

analysis of power network models in port-

Hamiltonian form

Cite as: Chaos 28, 101102 (2018); https://doi.org/10.1063/1.5054850

Submitted: 04 September 2018 . Accepted: 01 October 2018 . Published Online: 18 October 2018

Volker Mehrmann, Riccardo Morandin , Simona Olmi , and Eckehard Schöll

ARTICLES YOU MAY BE INTERESTED IN

Collective mode reductions for populations of coupled noisy oscillators

Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 101101 (2018); https://

doi.org/10.1063/1.5053576

Dissipative solitons for bistable delayed-feedback systems

Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 101103 (2018); https://

doi.org/10.1063/1.5062268

Power grid stability under perturbation of single nodes: Effects of heterogeneity and internal

nodes

Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 103120 (2018); https://

doi.org/10.1063/1.5040689

CHAOS 28, 101102 (2018)

Qualitative stability and synchronicity analysis of power network models in

port-Hamiltonian form

Volker Mehrmann,1,a)Riccardo Morandin,1,b)Simona Olmi,2,3,4,c)and Eckehard Schöll2,d)

1Institut für Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, Germany

2Institut für Theoretische Physik, Sekr. EW 7-1, TU Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany

3INRIA Sophia Antipolis Méditerranée, 2004 Route des Lucioles, 06902 Valbonne, France

4CNR—Consiglio Nazionale delle Ricerche—Istituto dei Sistemi Complessi, 50019 Sesto Fiorentino, Italy

(Received 4 September 2018; accepted 1 October 2018; published online 18 October 2018)

In view of highly decentralized and diversified power generation concepts, in particular with renew-

able energies, the analysis and control of the stability and the synchronization of power networks is

an important topic that requires different levels of modeling detail for different tasks. A frequently

used qualitative approach relies on simplified nonlinear network models like the Kuramoto model

with inertia. The usual formulation in the form of a system of coupled ordinary differential equations

is not always adequate. We present a new energy-based formulation of the Kuramoto model with

inertia as a polynomial port-Hamiltonian system of differential-algebraic equations, with a quadratic

Hamiltonian function including a generalized order parameter. This leads to a robust representation

of the system with respect to disturbances: it encodes the underlying physics, such as the dissipa-

tion inequality or the deviation from synchronicity, directly in the structure of the equations, and it

explicitly displays all possible constraints and allows for robust simulation methods. The model is

immersed into a system of model hierarchies that will be helpful for applying adaptive simulations

in future works. We illustrate the advantages of the modified modeling approach with analytics and

numerical results. Published by AIP Publishing. https://doi.org/10.1063/1.5054850

To reach the goal of temperature reduction to limit the

climate change, as stipulated at the Paris Conference in

2015, it is necessary to integrate renewable energy sources

into the existing power networks. Wind and solar power

are the most promising ones, but the integration into the

electric power grid remains an enormous challenge due

to their variability that requires storage facilities, back-

up plants, and accurate control processing. The current

approach to describe the dynamics of power grids in terms

of simplified nonlinear models, like the Kuramoto model

with inertia, may not be appropriate when different con-

trol and optimization tasks are needed to be addressed.

Under this aspect, we present a new energy-based formu-

lation of the Kuramoto model with inertia that allows for

an easy extension if further effects have to be included

and higher fidelity is required for qualitative analysis. We

illustrate the new modeling approach with analytic results

and numerical simulations carried out for a semi-realistic

model of the Italian grid and indicate how this approach

can be generalized to models of finer granularity.

I. INTRODUCTION

The increasing usage of renewable energy sources,

such as wind and solar power, and the decentralization of

power generation make the stability and synchronization con-

trol of modern power systems very difficult. To address

a)Electronic mail: [email protected]

b)Electronic mail: [email protected]

c)Electronic mail: [email protected]

d)Electronic mail: [email protected]

different control and optimization tasks, there are different

modeling approaches for the various components of power

networks; we will briefly present a model hierarchy for

synchronous generators and transmission lines. The simplest

model in the hierarchy is the Kuramoto model with inertia,

which is often used for a qualitative analysis of the network

behavior.1,2

The usual formulation of the Kuramoto model with iner-

tia in the form of a coupled system of ordinary differential

equations as in Refs. 1–7is, however, not always appropriate,

because physical properties like the conservation of energy

and momentum are only implicitly represented in the equa-

tions, and thus, in numerical simulation or control, they may

be violated and lead to un-physical behavior. To prevent this,

we present a new energy-based formulation of the Kuramoto

model with inertia as port-Hamiltonian system of differential-

algebraic equations (pHDAEs). This formulation is closely

related to the one in Refs. 8and 9, but it uses a change of vari-

ables, similar to Ref. 10, to turn the system into a polynomial

system with quadratic Hamiltonian.

Port-Hamiltonian (pH) systems (Refs 11–13)have

become a paradigm for multi-physics systems modeling incor-

porated in automated modeling frameworks such as 20-SIM

http://www.20sim.com/;seeRef.14. To include interface con-

ditions or constraints, in Refs. 15 and 16, ordinary pH systems

have been extended to port-Hamiltonian DAEs (pHDAEs).

Although the discussed systems are nonlinear, we present

here, for simplicity, the definition in a linear time-varying

form, which describes the local behavior in the neighborhood

of solutions; see Ref. 16 for the nonlinear version. Denote for

a compact interval I⊆Rby Cj(I×Rn,Rm)the set of j-times

continuously differentiable functions from I×Rnto Rm.

1054-1500/2018/28(10)/101102/6/$30.00 28, 101102-1 Published by AIP Publishing.

101102-2 Mehrmann

et al

.Chaos 28, 101102 (2018)

Definition 1. A linear variable coefficient system

E˙x=[(J−R)Q−ET]x+(B−P)u,

y=(B+P)TQx +(S+N)u,(1)

with E,Q∈C1(I×Rn,Rn,n),J,R,T∈C0(I×Rn,Rn,n),

B,P∈C0(I×Rn,Rn,m),S,N∈C0(I×Rn,Rm,m),andS=

ST,N=−NTis called port-Hamiltonian differential-

algebraic system (pHDAE) if

(i) Along all solutions x:I→Rnof (1),wehave

(QTE)[t,x(t)]=(ETQ)[t,x(t)]∈C1(I,Rn,n)

and for all t∈I,wehave

dt(QTE)[t,x(t)]=QT(ET −JQ)

+(ET −JQ)TQ[t,x(t)];

(ii) the Hamiltonian function H=1

2xT[(QTE)(t,x)]x∈

C1(I×Rn,R)satisfies H[t,x(t)]≥h0uniformly for

some h0∈R,allt∈Iand all solutions x(t)of (1);

(iii) for all t∈I,W=WT≥0, where

W:=QTRQ QTP

PTQS

∈C0(I×Rn,Rn+m,n+m).(2)

pHDAEs satisfy the dissipation inequality

H[t1,x(t1)]−H[t0,x(t0)]≤t1

y(t)Tu(t)dt,(3)

which shows that (1) is a passive system (see Ref. 17), and

since H(x)defines a Lyapunov function, pH systems are

implicitly Lyapunov stable. An advantage of pHDAEs in the

context of power system modeling is that they are closed

under power-conserving interconnection, which allows one

to build-up models in a modularized way (see Ref. 18), and

Galerkin projection (see Ref. 19), which allows systematic

discretization and model reduction.

II. MODEL HIERARCHIES FOR POWER NETWORKS

A power network consists of many different types of

components, e.g., the synchronous generators, the loads, and

the transmission lines. Each of these components can be mod-

eled with a different degree of fidelity, depending on the state

of the system, the dynamics of the network, the task that has to

be completed, the accuracy requirements, or the allowed com-

putation time. The innate modularity of pH systems allows

one to interconnect different components of the network while

retaining energy-consistency, even between different models

in the hierarchy, and can be easily exploited by automated

modeling frameworks. It can also be combined with structure-

preserving model reduction methods that allow one to speed

up computations while retaining qualitative behavior.

A detailed synchronous generator model including the

electromagnetic phenomena internal to the generator (see,

e.g., Refs. 9and 20) can be drastically simplified by neglect-

ing the dynamics of the magnetic variables, leading to the

swing equation,9which can be modeled with different accu-

racy (see, e.g., Ref. 21).

FIG. 1. Model hierarchies for synchronous generators (on the left) and

transmission lines (on the right).

Also, the model of a transmission line, represented as a

circuit with resistive, inductive, and capacitive components

(see, e.g., Refs. 2and 20) can be further simplified: both the

shunt capacitance and the series resistance are often negligi-

ble when compared with the inductance, and when assuming

steady-state, the phasor model can be used, substituting the

equations for the current through the transmission line by a

linear relation. See Fig. 1for a summary of the two model

hierarchies.

For the purpose of this paper, we consider a power

network consisting of nsynchronous generators, connected

through transmission lines, using the conventional swing

equation for the synchronous generators, and the phasor

L-model for the transmission lines. The resulting system [see

Ref. 9(subsection 6.2)] can be written as

mj˙ωj=−djωj+



k=1

gjk sin(θk−θj)+Pj,(4)

for j=1,...,n,whereθjis the phase shift, ωj=˙

θjis the

deviation from a given (nominal) frequency ,mj>0isthe

inertia, dj>0 is the damping constant, Pj∈Ris the power

balance of the j-th generator, gjk =BjkEjEk≥0, where Bjk

is the susceptance of the transmission line between the j-th

and k-th node, and Ej,Ekare the respective amplitudes of

the voltage potentials (assumed constant). Note that we allow

the power balance to be negative, so that a generator can

also be interpreted as a load. This model can be written as

an m+n-dimensional ordinary port-Hamiltonian system (see

Refs. 8and 9) with variables q=[θk−θj]∈Rm(indexed as

the edges of the network graph) and p=[mjωj]∈Rn:

˙q=ATω,

M˙ω=−Dω−Asinq+P,(5)

with Hamiltonian

H(q,p):=1

2pTM−1p−1Tcosq,(6)

where M=diag(mj),D=diag(dj)∈Rn,n,=diag[gjk]∈

Rm,m,P=[Pj]∈Rn,A∈Rn,mis the incidence matrix of the

graph, and 1is the vector of all ones.

System (4) is a generalization of the canonical second-

order Kuramoto model with inertia, which consists of a

101102-3 Mehrmann

et al

.Chaos 28, 101102 (2018)

system of nfully-coupled oscillators satisfying

mj¨

θj+dj˙

θj=j+K



k=1

sin(θk−θj),(7)

for j=1,...,n,whereθj,mj,djare as before, jare the

natural frequencies, and K>0 is the coupling constant of

the system. We will refer to system (4) as the generalized

Kuramoto model (GKM). To study qualitatively the synchro-

nization of oscillatory networks, typically the complex order

parameter (COP)

reiφ:=1



j=1

eiθj∈{z∈C:|z|≤1}

is used where r=1 when the system is in a fully synchronized

state and 0 when it is completely desynchronized.

III. THE GENERALIZED ORDER PARAMETER

While the COP is linked to the fully-coupled canonical

Kuramoto model with inertia, (see Ref. 3)inthecaseofsparse

networks, it can show weird behavior (see, e.g., Fig. 3or

Ref. 22), because it does not take the topology of the system

into account, and so it fails to capture well the transitions in

the partially phase locked regime. To address this problem, we

introduce a generalized order parameter (GOP) that applies

to the GKM and is consistent with the COP. Treating (7) as a

special case of (4) with a complete graph and gjk ≡K, we can

write

r2=|reiφ|2=1

n2

j,k

ei(θk−θj)=1

n2

j,k

cos(θk−θj)

n2

j<k

cos(θk−θj)+1

n=2

Kn21Tcosq+1

thus, in this case, 1Tcos q=1

2Kn2r2−1

n. Since the pH

formulation does not change if we add a constant to its

Hamiltonian, we can replace (6) with

H1(q,p):=1

2pTM−1p+1

2Kn2[1 −r(q)2]≥0. (8)

To extend the definition of the order parameter to the gen-

eral case, without losing its connection to the Hamiltonian,

we define the GOP ξ(q):=a1Tcosq+bfor some network-

dependent constants a,b∈R, which should satisfy the

following conditions:

(1) In the case of system (7),wehaveξ(q)≡r2(q);

(2) the maximum value of ξ(q)is 1, and it corresponds to the

fully synchronized state;

(3) the minimum value of ξ(q)is greater or equal to a fixed

value ξ0∈R.

We can restate these conditions in terms of aand b:

(1) Since in (7) we have r2=2

Kn21Tcosq+1

n,forthat

network we must have a=2/Kn2and b=1/n.

(2) The fully synchronized state ¯qcorresponds to cos ¯q=

1and ξ(¯q)=a1T1+b,soξ(¯q)is the maximum if

and only if a≥0, and it is equal to 1 if and only if

b=1−a1T1.

(3) If we assume (2), then ξ(q)≥−a1T1+b=2b−1for

all q. Whether this value can actually be reached depends

on the network, however, there is no clear general formula

for minqξ(q). In any case, b≥1+ξ0

2is sufficient to imply

condition (3).

Using the given freedom, we choose b=1/nand a=(n−

1)/(n1T1), so that conditions (1)–(3) are satisfied with ξ0=

−1+1

n, in particular ξ(q)∈(−1,1]. Note that typically we

will have ξ(q)>0. We can express the Hamiltonian (6) in

terms of the GOP. Changing H(q,p)again by an additive

constant, we obtain

H2(q,p):=1

2pTM−1p+n

n−11T1[1 −ξ(q)](9)

2pTM−1p−1Tcosq+1T1, (10)

which has 0 as its minimum. We will show in Sec. Vthat the

GOP ξ(q)behaves better than the ordinary COP, especially for

sparsely coupled networks. Note that a similar GOP has been

formulated simultaneously to Ref. 23 in Ref. 22; however, it

has not been related to a pH formulation nor applied to power

grid networks.

IV. A NEW PHDAE FORMULATION FOR THE GKM

Similar to Ref. 10, we introduce (ρ,σ) :=(cosθ,sinθ),

viewed componentwise, and express (4) in x=(ω,ρ,σ).

Using gjk =0 when there is no edge connecting the j-th and

k-th node, and introducing again P=[Pj]inRn,D=

diag(dj),andG=[gjk]inRn,n, and the matrix functions Dρ=

diag(ρj),Dσ=diag(σj)with values in Rn,n, we can write the

GKM as

M˙ω=−Dω−DσGρ+DρGσ+P,

˙ρ=−Dσω,

˙σ=Dρω.

(11)

Since all terms with coefficients gjj for j=1,...,n

vanish in (11), we are free to choose gjj; our choice is

gjj =−k=jgjk,sothatG≤0. This can be interpreted

as a pHDAE formulation E˙x=(J−R)Qx +Bu with x=

[ωT,ρT,σT]T,u=P,E=diag(M,In,In),0≤R=RT=

diag(D,0n,0n),Q=diag(In,−G,−G),

J=−JT=⎡

⎣

0Dσ−Dρ

−Dσ00

Dρ00

⎤

⎦,B=⎡

⎣

0⎤

⎦,

and Hamiltonian function

H3(x)=1

2xTQTEx =1

2ωTMω−1

2ρTGρ−1

2σTGσ.

In particular, we have

2ρTGρ+1

2σTGσ=1

2

j,k

gjk cos(θk−θj)

=

j<k

gjk[cos(θk−θj)−1]

=1Tcosq−1T1,

101102-4 Mehrmann

et al

.Chaos 28, 101102 (2018)

so that H3(x)=H2(q,p). This model has a larger state space

than a modified version of (5) (3ninstead of 2n), but it has

several advantages.

• Using (ρ,σ)instead of θas variables helps removing some

redundancy from the state, without having to adjust θj∈

[−π,π] numerically at each step. One can also consider a

formulation with angles defined relative to a fixed oscilla-

tor; then, the zero state is the only fully synchronized state,

and the matrix Gbecomes invertible (see Ref. 23).

• Although (11) is still a nonlinear system, it is now poly-

nomial and has a quadratic Hamiltonian with respect to x;

an advantage of this is, e.g., that, for P=0andD=0, it

can be integrated preserving H2without discretization error

(see Ref. 24).

• While (11) is still an ODE, we are allowed to add algebraic

equations without losing the port-Hamiltonian structure.

These could be, e.g., the relations linking the mechan-

ical variables ω,ρ,σto the electrical variables (voltage

and currents), Kirchhoff’s laws, or the implicit relation

ρ2

j+σ2

j=1.

Since we omitted the dependence of ρand σon θ,the

variables ρ,σimplicitly need to satisfy ρ2

j+σ2

j=1. Unfor-

tunately, in numerical integrators for (11), such implicit rela-

tions are typically not preserved, due to discretization and

roundoff errors.25 To address this, we can add these conditions

explicitly to the system, by introducing Lagrange multipli-

ers μ=[μj]∈Rn, adding the equation 0 =Dρρ+Dσσ−μ

to the system, and modifying correspondingly the matrices

E,J,R,Q,B(see Ref. 23 for details). One must also make sure

that the initial conditions are consistent and include μ=1.

Employing geometric integrators helps one to preserve these

conditions, as we show in Sec. V.

V. NUMERICAL RESULTS

In this section, we illustrate the advantages of our

new pHDAE formulation by some numerical examples. We

use a model of the Italian high-voltage (380 kV) power

grid of n=127 nodes, divided in 34 sources and 93 con-

sumers, connected by 342 links, (see Ref. 26) and low aver-

age connectivity nc=2.865. See http://www.geni.org and

https://www.entsoe.eu/map/Pages/no-webgl.html for a map

and the data.

We use the GKM and assume for simplicity that each

node has the same inertia and damping constant (mj=m=6

and dj=d=1, j=1,...,n), and each transmission line has

the same coupling coefficient, i.e., either gjk =Kor gjk =0.

To have a stable fully locked state as a possible solution of (4)

(i.e., ω(t)≡0), it is necessary that n

j=1Pj=0, which we

get, e.g., by setting Pj=−PC=−1 for all consumers, and

Pj=PG=2.7353 for all generators. This corresponds to a

Kuramoto model with inertia with perfectly balanced bimodal

δ-distribution of the natural frequencies. Since in a real power

grid the power of the generators (or consumers) will not have

exactly the same values, we analyzed also the case where

the Pjare chosen as random variables with distribution given

by the superposition of two almost non-overlapping Gaussian

distributions. See Ref. 23 for further computational results.

FIG. 2. Conservation of Lagrange multiplier μfor different integration

schemes [Symplectic Euler (SE), Implicit Midpoint (IMP), Explicit Mid-

point (EMP)] and K=0.5. The orange and magenta curves differ in the

convergence threshold in the nonlinear solver (10−6and 10−12).

We have integrated (11) with different integration schemes

and compared the results for the symplectic Euler method,

the implicit and the explicit midpoint method, as well as a

4th order explicit Runge-Kutta scheme applied to the old

model (4). We used stepsize 0.002, transient time 100, and

simulation time 5000. As illustrated in Fig. 2, on short time

intervals, the different solvers deliver similar results for the

multiplier μ. However, if we want its exact conservation, a

geometric integrator like the implicit midpoint method should

be applied with a sufficient accuracy for the nonlinear solver at

every integration step. When the accuracy is not sufficient, the

multiplier drifts from the value 1, as for the explicit midpoint

method.

To understand the transition from the non-synchronized

state for small Kto synchronized state for large K,we

performed sequences of simulations by increasing the param-

eter Kwith random initial conditions for {θi}and {˙

θi},for

FIG. 3. Time-averaged order parameters vs coupling constant K, with (a) per-

fectly balanced bimodal δ-distribution or (b) Gaussian distributions of the

j.rKI is the average standard Kuramoto parameter, calculated via a 4th

order Runge-Kutta scheme in the original formulation [i.e., via integrating

Eq. (4)]. ¯rand ¯

ξ1/2are the average standard Kuramoto parameter and gener-

alized order parameter, respectively, obtained with the implicit midpoint rule

applied to the pHDAE formulation. The data have been obtained by changing

adiabatically the coupling constant, starting from K=0 and with K=0.5.

101102-5 Mehrmann

et al

.Chaos 28, 101102 (2018)

FIG. 4. Average order parameter ¯

ξ1/2vs the damping constant dfor different

coupling constants K.

both (4) and (11). In each case, the simulation was ini-

tialized by employing the last values of the previous sim-

ulation. The average COP ¯rand its GOP counterpart ¯

ξ1/2

(assuming ξ(q)≥0) both show a non-monotonic behav-

ior in K, especially the former (see Fig. 3). For the

bimodal δ-distributions [panel (a)] and small coupling val-

ues ¯r∝1/√N, we observe an abrupt jump for K=6.5.

Subsequently, ¯rdecreases, reaching a minimum at K=9.

For larger K,then¯rincreases monotonically toward full

synchronization. The generalized order parameter ¯

ξ1/2,on

the other hand, does not show such an irregular behav-

ior for small coupling values but has almost constant

value until K=6.5, where the transition to synchronization

takes place and from then on rapidly increases toward 1.

Since ξtakes into account the topology of the network, the

analysis of the synchronization level is more stable for small

coupling constants, while it is still able to identify the tran-

sition to synchronization at K=6.5. The correctness of the

critical value for the transition to synchronization can be

confirmed indirectly by calculating the maximal Lyapunov

exponent (see Ref. 23). The COP ¯rfails in identifying the cor-

rect level of synchronization because, for K<6.5, the system

splits in two clusters: one consisting of the sources, oscillat-

ing close to their proper frequency, and one consisting of the

consumers, which rotates with negative average velocity. This

non-trivial form of partial synchronization already present in

the network is not detected by ¯r, whose value around 1/√N

indicates that the system behaves asynchronously. By increas-

ing the coupling to K≥7.5, the two clusters merge into a

single cluster, which is reflected in a monotonic increase of

the average order parameters (see Ref. 23). For a bimodal

Gaussian frequency distribution, the average order param-

eters ¯r,¯

ξ1/2as a function of the coupling constant reveal

that, in this setup, it is more difficult to achieve synchroniza-

tion, due to the inhomogeneity of the natural frequencies [see

Fig. 3(b)]. Here, ¯ris still irregular and unstable, irrespective of

whether we use the original Kuramoto model with inertia or

the pHDAE formulation, while ¯

ξ1/2is stable and informative,

with a continuous transition to synchronization, instead of a

jump from partial to full synchronization. It should be noted

that, while in these examples the pHDAE and the original

FIG. 5. Time tsnecessary to reach good synchronization (ξ1/2>0.95), as

a function of the parameters K,m,d, which represent the values to which

the parameters are tuned initially in order to calculate the time needed to

reach synchronization. To efficiently control the system, when stability is

affected by emerging disturbances, it is more convenient to increase the cou-

pling constant. On the other hand, bigger inertia prevents the system from fast

synchronization recovery since it slows down the dynamics. tsis expressed in

nondimensional time units.

formulation present similar qualitative behaviour [with a few

noticeable differences in Fig. 3(b)], the discretization of the

quadratic pHDAE with the implicit midpoint is usually more

adequate in limiting situations.

We conclude this section by testing how adjusting the

inertia mand the damping constant dcould help with the syn-

chronicity control of the system. For larger inertia, the critical

coupling value for transition to synchronization increases, so

that larger coupling strength is necessary to achieve the syn-

chronized state; see also Ref. 23. On the other hand, for fixed

K, larger damping constants dhelp reaching synchronization,

as shown in Fig. 4. However, a consistent increase is possible

only for Kclose to the critical value K=6.5 as in Fig. 3, thus

illustrating that the coupling strength plays a more important

role in the transition than the damping constant. This effect

is also observed when varying the mass m, while keeping K

constant (see Ref. 23).

When a generator is diverging from the synchronized

regime, it is important to react to this perturbation as fast as

possible in order to avoid a shut-down of the generator. Start-

ing from a slightly asynchronous model (K=6), we can also

control the time to synchronization tsby changing K,d,or

m;seeFig.5. Adjusting mechanically some parameters of the

system can thus be useful both to not drift away from syn-

chronicity and to return to a synchronous state in a shorter

time.

VI. CONCLUSIONS

We have presented a new port-Hamiltonian differential-

algebraic formulation of the Kuramoto model of coupled

rotators, as well as a new definition of the order parameter.

The new pHDAE model formulation is structured in such a

way that it can be easily extended to analogous pHDAE for-

mulations of finer granularity as they are used in quantitative

101102-6 Mehrmann

et al

.Chaos 28, 101102 (2018)

stability and synchrony analysis of power systems. The new

order parameter is more robust in limiting situations. We have

also shown the advantage of the port-Hamiltonian formulation

in the preservation of conserved quantities. Future work will

include the analysis of the whole model hierarchy, whose low-

est level is represented by the Kuramoto model with inertia

described here.

ACKNOWLEDGMENTS

V.M. acknowledges support by Deutsche Forschungsge-

meinschaft (DFG) via Project A2 of SFB 910 and Einstein

Foundation Berlin via the Einstein Center ECMath. R.M. is

supported by Einstein Foundation Berlin via the Einstein Cen-

ter ECMath. S.O. and E.S. were supported by DFG via Project

A1 of SFB 910.

1G. Filatrella, A. H. Nielsen, and N. F. Pedersen, Eur. Phys. J. B 61, 485–491

(2008).

2P. Kundur, N. J. Balu, and M. G. Lauby, Power System Stability

and Control, EPRI Power System Engineering Series (McGraw-Hill,

1994).

3S. Olmi, A. Navas, S. Boccaletti, and A. Torcini, Phys.Rev.E90, 042905

(2014).

4F. Salam, J. Marsden, and P. Varaiya, IEEE Trans. Circuits Syst. 8, 673–688

(1984).

5T. Nishikawa and A. E. Motter, New J. Phys. 17, 015012 (2015).

6M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Phys. Rev. Lett. 6,

064101 (2012).

7F. Doerfler, M. Chertkov, and F. Bullo, Proc. Natl. Acad. Sci. 6, 2005–2010

(2013).

8N. Monshizadeh, C. De Persis, A. J. van der Schaft, and J. M. A. Scherpen,

IEEE Trans. Automat. Control 63, 1288–1299 (2018).

9A. J. van der Schaft and T. Stegink, “Perspectives in modeling for control

of power networks,” Annu. Rev. Control 41, 119–132 (2016).

10A. Sarlette and R. Sepulchre, “Synchronization on the circle,” in The Com-

plexity of Dynamical Systems: A Multi-disciplinary Perspective, edited

by J. Dubbeldam, K. Green, D. Lenstra, (Wiley, 2011), available at

https://arxiv.org/pdf/0901.2408.pdf.

11R. Ortega, A. J. van der Schaft, Y. Mareels, and B. M. Maschke, Control

Syst. Mag. 21, 18–33 (2001).

12G. Golo, A. J. van der Schaft, P. C. Breedveld, and B. M. Maschke, Nonlin-

ear and Hybrid Systems in Automotive Control, edited by A. Rantzer and

R. Johansson (Springer, Heidelberg, 2003), pp. 351–372.

13A. J. van der Schaft, Advanced Dynamics and Control of Structures and

Machines (Springer, New York, 2004), Vol. 444.

14P. C. Breedveld, Modeling and Simulation of Dynamic Systems Using Bond

Graphs (EOLSS Publishers Co. Ltd./UNESCO, Oxford, UK, 2008).

15A. J. van der Schaft, Surveys in Differential-Algebraic Equations I

(Springer, Berlin, 2013), pp. 173–226.

16C. Beattie, V. Mehrmann, H. Xu, and H. Zwart, preprint arXiv:1705.09081

(2017).

17C. I. Byrnes, A. Isidori, and J. C. Willems, IEEE Trans. Automat. Control

36, 1228–1240 (1991).

18J. Cervera, A. J. van der Schaft, and A. Baños, Automatica 43, 212–225

(2007).

19S. Gugercin, R. V. Polyuga, C. Beattie, and A. J. van der Schaft,

Automatica 48, 1963–1974 (2012).

20S. Fiaz, D. Zonetti, R. Ortega, J. M. A. Scherpen, and A. J. van der Schaft,

Eur. J. Control 19, 477–485 (2013).

21P. Monshizadeh, C. De Persis, N. van der Monshizadeh, and A. J. Schaft, in

2016 IEEE 55th Conference on Decision and Control (CDC) (IEEE, 2016),

pp. 4116–4121.

22M. Schröder, M. Timme, and D. Witthaut, Chaos: Interdiscip. J. Nonlinear

Sci. 27(7), 073119 (2017).

23V. Mehrmann, R. Morandin, S. Olmi, and E. Schöll, preprint

arXiv:1712.03160 (2017).

24L. Brugnano and F. Iavernaro, Line Integral Methods for Conservative

Problems (Chapman & Hall, 2015).

25P. Kunkel and V. Mehrmann, Differential-Algebraic Equations: Analysis

and Numerical Solution (EMS Publishing House, Zürich, 2006).

26L. Fortuna, M. Frasca, and A. Sarra Fiore, Int. J. Modern Phys. 26,25

(2012).