Power System Transient Stability Analysis Using Sum

Of Squares Programming
M Tacchi, B. Marinescu, M Anghel, S Kundu, S. Benahmed, C Cardozo

To cite this version:

M Tacchi, B. Marinescu, M Anghel, S Kundu, S. Benahmed, et al.. Power System Transient Stability
Analysis Using Sum Of Squares Programming. Power Systems Computation Conference - PSCC, Jun
2018, Dublin, Ireland. �10.23919/PSCC.2018.8442484�. �hal-02516468�

HAL Id: hal-02516468

https://hal.archives-ouvertes.fr/hal-02516468
Submitted on 23 Mar 2020

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Power System Transient Stability Analysis Using
Sum Of Squares Programming
M. Tacchi∗ , B. Marinescu∗ , M. Anghel† , S. Kundu‡ , S. Benahmed∗ and C. Cardozo§
∗
Ecole Centrale Nantes-LS2N-CNRS, France
[email protected], [email protected], [email protected]
† Los Alamos National Laboratory, USA

[email protected]
‡ Pacific Northwest National Laboratory, USA

[email protected]
§ RTE-R&D, France

[email protected]

Abstract—Transient stability is an important issue in power do not exist for grids with losses [1].
systems but difficult to quantify analytically. Most of the ap-
proaches in this case lie on a very simplified model of the Recent advances in semidefinite programming relaxation
generators, usually reduced to the swing equation. In this work,
a more detailed model which includes voltage dynamics and both and the use of the Sum Of Squares (SOS) decomposition
voltage and frequency regulators is considered to get more real- to check non-negativity have opened the way for efficient
istic results. The used methodology for algorithmic construction algorithmic analysis of systems with polynomial vector fields
of Lyapunov functions is based on recent advances in the field of [3]. This approach was used in [2] for transient stability
positive polynomials. This analysis framework uses an algebraic analysis of a power system. However, the generators were
reformulation technique that recasts the systems dynamics into a
set of polynomial differential algebraic equations in conjunction modeled only by their swing equations. As the voltage and
with a sum of Squares method to search for a Lyapunov function. frequency regulators have an important impact on stability
Next, linear matrix inequalities are used to expand the region analysis, in the present paper this approach is extended by
of attraction of the considered equilibrium point. The results considering a generator with voltage dynamics and both
are checked against an experimental (by extensive numerical voltage and frequency regulations. This also provides a way
simulations) evaluation of the region of attraction.
to take into account high penetration of power electronics
Index Terms—Lyapunov theory, nonlinear systems, power elements in the grid due to the integration of renewable
systems analysis, region of attraction, sum of squares. energies and HVDC transmission lines, thus having an
important impact on the transient stability of the system.
I. I NTRODUCTION SOS method gives an analytical solution for the construction
of Lypunov functions in order to estimate the Region
Transient stability is one of the most important power Of Attraction (ROA) for a locally asymptotically stable
systems analysis problems. From a physical viewpoint, equilibrium point of the system.
transient stability can be defined as the ability of a power
system to maintain its synchronism when subjected to large The paper is organised as follows: the problem is formulated
and transient disturbances [1]. Widespread assessment tools in Section II on a Single Machine Infinite Bus (SMIB) system.
are generally based on indirect methods which rely on the SOS formalism is brieffly recalled in Section III. The change
numerical integration of nonlinear differential equations of variables needed to reach a SOS form (i.e., recasting the
describing system dynamics. However, this approach is not model of the system with trigonometric non linearities into a
suited to synthesize controllers with a direct quantification set of polynomial differential algebraic equations) is given in
of stability margin since it does not provide an analytical Section IV. The recasting procedure is proved to be a Lie-
characterization of stability [2]. Alternatively, direct methods Bäcklund transformation, which means that the transformed
are based on the estimation of the stability domain of the system has equivalent trajectories and stability properties [4].
equilibrium point; they ensure that all the trajectories initiated Next, in section V we relax Lyapunov's conditions for stability
in this domain converge to the equilibrium point [2]. The and model constraint equations to suitable SOS conditions
main drawback of direct methods is that they rely on the using theorems from real algebraic geometry in order to
identification of Lyapunov functions which are hard to formulate the problem as an optimization one. Hence, a
determine. Indeed, there is not a systematic method for the Lyapunov function for the asymptotically stable equilibrium
construction of Lyapunov functions. Moreover, it has been point is constructed using the expanding interior algorithm
shown that quadratic (i.e., energy type) Lyapunov functions developed in [3]. An estimate of the ROA is given by a
level set of the Lyapunov function. We finally test the ROA where Ta , Ka are the parameters of the voltage regulator
estimation error in Section VI by numerically computing the (AVR), Vref is the reference, and
real one in all state directions via full nonlinear simulations. q
The codes are implemented in MATLAB using SOSTOOLS Vt = vd2 + vq2 (3a)
[5], [6] which is a free, third-party MATLAB toolbox that dPm
solves SOS problems. The analysis is carried out in the case Tg = −Pm + Pref + Kg (ωref − ω) (3b)
dt
of a single machine-infinite bus system.
where Tg , Kg are the parameters of the turbine regulator
(governor) and Pref , ωref are the references.
II. P ROBLEM FORMULATION We then introduce a temporary short-circuit at node S (see
Figure 1), according to the following protocol:
We consider a synchronous machine G connected to an
• At a time tc a short-circuit occurs and we switch from the
infinite bus N∞ through two lines in parallel (see Figure
1). The infinite bus imposes a nominal voltage of amplitude nominal system (Figure 1) to a new short-circuited system
Vs = 1pu and frequency ωs = 1pu. The synchronous machine (Figure 2), and thus we are no longer at an equilibrium
has two rotating axes d (direct) and q (quadratic), which point
• The system leaves the equilibrium point of the nominal
support the currents id and iq , generating a voltage (vd , vq ).
The synchronous machine is characterized by a resistance model and follows the short-circuit equations (see below)
r and equivalent reactances xd , x0d and xq . It receives a for a while, until the short-circuit is eliminated
• At a time tcl = tc + ∆t, we switch back to the nominal
mechanical power Pm which makes it rotate at frequency ω,
generating e.m.f e0q with voltage field Ef d . Let δ be the phase topology and equations; the problem is to know whether
between the machine and the infinite bus and H the machine’s the system reaches an equilibrium point or not, i.e.
mechanic inertia constant. whether it is still in the ROA of one equilibrium points
of the nominal model or not.
2(R + jX)
2(R + jX)

G N∞
∞ G N∞
R + jX R + jX ∞
S R + jX R + jX
VG ∼ (vd , vq ) V∞ ∼ (Vs , ωs ) S
VG ∼ (vd , vq ) V∞ ∼ (Vs , ωs )

Figure 1: Synchronous machine connected to an infinite bus

Figure 2: The short-circuited system

The equations describing the dynamics of this system are During the short-circuit the system’s equations are the same
given in [7] (p.105). They can be written as follows as (1), but with (R, X) replaced by 32 (R, X) and Vs replaced
by 13 Vs :
de0q
Td0 0 = −e0q − (xd − x0d )id + Ef d (1a) de0q
dt Td0 0 = −e0q − (xd − x0d )id + Ef d (4a)
dω dt
2H = Pm − (vd id + vq iq + ri2d + ri2q ) (1b) dω
dt 2H = Pm − (vd id + vq iq + ri2d + ri2q ) (4b)
dδ dt
= ω − ωs (1c) dδ
dt = ω − ωs (4c)
(X + x0d )Vs sin δ − (R + r)(Vs cos δ − e0q ) dt
iq = (1d) 3x0d )Vs
sin δ − (2R + 3r)(Vs cos δ − 3e0q )
(R + r)2 + (X + x0d )(X + xq ) iq =
(2X +
(4d)
X + xq 1 (2R + 3r)2 + (2X + 3x0d )(2X + 3xq )
id = iq − Vs sin δ (1e) 2X + 3xq 1
R+r R+r id = iq − Vs sin δ (4e)
vd = xq iq − rid (1f) 2R + 3r 2R + 3r
vd = xq iq − rid (4f)
vq = Riq + Xid + Vs cos(δ) (1g)
3vq = 2Riq + 2Xid + Vs cos(δ). (4g)
where Td0 0 is a characteristic time.
Finally, the regulation equations (2), (3a) and (3b) are not
The machine is governed by two regulators whose equations affected by the short-circuit.
are the following The aim of this work is to estimate the ROA of a given
dEf d equilibrium point so as to be able to give an analytical stability
Ta = −Ef d + Ka (Vref − Vt ) (2) assessment of the power system.
dt
 III. T HE SOS APPROACH where p ∈ Σn and V is a yet unknown Lyapunov function
SOS are real polynomials that can be written as sums of of which level set Ω1 = {x ∈ Ω|V (x) ≤ 1} approximates the
k ROA. In order for ∆ to satisfy the conditions of Lyapunov’s
p2j where p1 , . . . , pk ∈ Rn .
P
squares of polynomials: s = theorem, we must have 0 < V ≤ 1 =⇒ V̇ < 01 . Since we
j=1
Here, Rn is the set of real polynomials in n variables and we want V to be positive definite over ∆ while V and ∆ are
denote by Σn the set of SOS in n variables. initially unknown, the best way to ensure that is to look for
The SOS approach is detailed in [3]. It is based on the a V which is positive definite on Rn . This problem can be
Positivstellensatz (P-satz, [8]): written with set emptiness constraints as

Let g1 , . . . , gβ , f1 , . . . , fα , h1 , . . . , hγ ∈ Rn . Then, the

following statements are equivalent: maximize β (6a)
V ∈Rn ,V (0)=0
• The set
  s.t., {x ∈ Rn |V (x) ≤ 0, x 6= 0} = ∅ (6b)
f1 (x) ≥ 0, . . . , fα (x) ≥ 0  n
{x ∈ R |p(x) ≤ β, V (x) ≥ 1, V (x) 6= 1} = ∅ (6c)

x ∈ Rn g1 (x) 6= 0, . . . , gβ (x) 6= 0

h1 (x) = 0, . . . , hγ (x) = 0
 {x ∈ Rn |V (x) ≤ 1, V̇ (x) ≥ 0, x 6= 0} = ∅ (6d)
is empty.
which can be reformulated as
• ∃F ∈ P(f ), G ∈ M(g), H ∈ I(h) ;
F + G2 + H = 0 (5) maximize β
V ∈Rn ,V (0)=0
(such F, G, H are called P-satz certificates, or P-satz s.t., {x ∈ Rn | − V (x) ≥ 0, `1 (x) 6= 0} = ∅
refutations), where n
{x ∈ R |β − p(x) ≥ 0, V (x) − 1 ≥ 0, V (x) − 1 6= 0} = ∅

Yβ

 {x ∈ Rn |1 − V (x) ≥ 0, V̇ (x) ≥ 0, `2 (x) 6= 0} = ∅
k
M(g) := gj j |k1 , . . . , kβ ∈ N

j=1
 with `1 , `2 ∈ Σn are positive definite ensuring the non
polynomial constraint x 6= 0. Then, we can use the P-satz
with the convention M(∅) = {1}. and a simplification given in [3] to write the problem in a
form which is suitable for the algorithm presented above:
( P
)
X
P(f ) := s0 + si bi |P ∈ N, s ∈ ΣP
n
+1
,b ∈ M(f ) P

i=1
maximize β (7a)
V ∈Rn ,V (0)=0, s1 ,s2 ,s3 ∈Σn
with the convention P(∅) = Σn .
( γ ) subject to, V − `1 = s4 ∈ Σn (7b)
− ((β − p)s1 + (V − 1)) = s5 ∈ Σn
X
I(h) := hk pk |p1 , . . . , pγ ∈ Rn (7c)
 
i=1 − (1 − V )s2 + V̇ s3 + `2 = s6 ∈ Σn . (7d)
with the convention I(∅) = {0}.
Since Σn and the set of positive semi-definite matrices are
The idea is simply to estimate the ROA of the system using isomorphic, this is reduced to an LMI problem which can be
a set whose complement can be described like in the P-satz solved within an iterative algorithm (e.g., [9], [10]) (see Figure
formulation. Indeed, we know that for a given vector field 3).
F : Rn → Rn , an equilibrium point x (without any loss of
Some features of SOSTOOLS include the setting of poly-
generality, we assume x = 0) of the equation ẋ = F (x) is
nomial optimization problems and the search for a polynomial
stable iff there exists a function V defined on a neighbourhood
Lyapunov function after expressing the SOS problem as a LMI
of 0 (let us call it Ω), and a domain D ⊂ Ω on which V is
feasibility problem. In [2], SOSTOOLS is used to implement
an LF. In such case, any level set Ωβ = {x ∈ Ω|V (x) ≤ β}
the expanding interior algorithm for a dynamic model without
such that Ωβ ⊂ D is a positively invariant subset of the
regulation.
system’s ROA. The SOS approach consists in finding the
largest possible Ωβ ⊂ D with varying V and D, which
is possible as soon as F ∈ Rn . Two algorithms allowing IV. R ECASTING THE POWER SYSTEMS MODEL
to perform this research are discussed in [3] and [2]: the
expanding D algorithm and the expanding interior algorithm. The system (1), (2), (3a), (3b) can be reformulated as an
Since the latter is more efficient than the former, we only autonomous ODE:
implemented the expanding interior algorithm.
ẏ = F (y) (8)
The aim of the SOS program is to find the largest
Pβ := {x ∈ Rn |p(x) ≤ β} ⊂ ∆ := {x ∈ Rn |V (x) ≤ 1} , 1 With the well known convention V̇ = ∇V · F
 choose β (0) > 0 and V (0) where Y eq is the value of the variable Y at the considered
Lyapunov function on equilibrium point. One can then easily compute the dynamics
∆ = {x ∈ Rn |V (0) (x) ≤ 1}
such that of the recasted system:
(0)
Pβ := {p(x) ≤ β (0) } ⊂ ∆
ż = H(z) (11)
(i+1)
search for s6,8,9 ∈ Σn and a with H : R8 −→ R8 defined by:
maximal β (i+1) > β (i) satisfying 
(7) for V (i+1) = V (i) 
 H1 (z) = (1 − z2 )z3

linear search for β



 H2 (z) = z1 z3
(i+2)
∈ Σn ,

 1
search for s6  H3 (y) =
 (z6 + Pref − (vd (z)id (z)+
V (i+2) ∈ Rn and a 2H



maximal β (i+2) > β (i+1) satisfying
(z) + rid (z)2 + riq (z)2 )


vq (z)i

 q
(i+2) (i+1)

(7) for s8,9 = s8,9

 
eq
 1 0eq 0
H4 (y) = Td00  −(z4 + eq ) − (xd − xd )id (z) + z5 + Ef d



H5 (y) = T1a −(z5 + Efeqd ) + Ka (Vref − (z7 + V eq )

(12)
test β (i+2) − β (i+1) < βtol 


 H6 (y) = T1g (−z6 − Kg z3 )
H7 (z) = z8 + V1eq

P


 (vd (z)∂zi vd (z) +
i∈{1,2,4}


The resulting 

YES ∆ = {x ∈ Rn |V (i) (x) ≤ 1} NO v (z)∂ v (z)) · Hi (z)



3 qP zi q

is an estimate of 0’s R.O.A. 
1




 H8 (z) = − z8 + V eq (vd (z)∂zi vd (z)+

 i∈{1,2,4}
Figure 3: The expanding interior algorithm 

vq (z)∂zi vq (z)) · Hi (z)

  and iq (z) = iq (y), id (z) = id (y), vq (z) = vq (y) and

δ vd (z) = vd (y).
 ω 
 0  5
where y =   eq  ∈ R is the state vector and the vector The size of the recasted system has a direct impact on the

 Ef d 
computation burden of the ROA estimation. Let us notice that,
Pm 2
if (Vref − V (y)) in F4 (y) is replaced by (Vref − V (y)2 ), i.e.,
field F : R5 −→ R5 is defined by: if
1
 F4 (y) = (−y4 + Ka (Vref − Vt (y))) (13)

 F1 (y) = y2 − ωs Ta
1




 F2 (y) = (y5 − (vd (y)id (y) + vq (y)iq (y)+ is used in (9), we no longer need to introduce z7 and z8 in the


 2H recasted system which consists in this case in only the first 6
rid (y)2 + riq (y)2 )


(9) equations of (12):
1 0
F3 (y) = Td00 (−y3 − (xd − xd )id (y) + y4 )

  
 H1 (z)
F4 (y) = T1 (−y4 + Ka (Vref − Vt (y)))



 a
H2 (z)
F5 (y) = 1 (−y5 + Pref + Kg (ωref − y2 ))
  
H3 (z)
Tg
ż = H 0 (z) = 
H4 (z) .
 (14)
 
with iq (y), id (y), vq (y) and vd (y) introduced in (1d), (1e), H5 (z)
(1g) and (1f). H6 (z)
Here one can see that F is not a polynomial vector field, From a physical point of view, instead of comparing the
as requested in the SOS approach. However, the following magnitude (or, modulus) of the voltage (phasor) over the
variables allowed us to recast the system as a polynomial ODE: synchronous machine to a reference, we are comparing its
squared magnitude to the square of the reference. As the model
z1 = sin(δ − δ eq ) (10a) is written in per-unit variables, voltages are around 1 and this
eq approximation has little impact (this has also been checked in
z2 = 1 − cos(δ − δ ) (10b)
z3 = ω − ωs (10c) simulation).
Thus, H 0 is a polynomial vector field: H 0 ∈ R6 . However,
z4 = e0q − e0eq
q (10d)
for the new equations to model the same system as the
z5 = Ef d − Efeqd (10e) former, it is necessary to make sure that the change of
z6 = Pm − Pref (10f) variable is a Lie-Bäcklund transformation (see [4]). Indeed,
z7 = Vt − V eq (10g) this is necessary to ensure that the new equations and the
1 1 former ones have the same trajectories.
z8 = − eq (10h)
Vt V
 Definition : Let M be a smooth manifold, possibly of From this definition, it is obvious that an endogenous
infinite dimension, and F : M −→ T M a smooth vector field transformation (which is a particular case of Lie-Bäcklund
on M. isomorphism) preserves the trajectories (and so the stability
The pair (M, F ) is a system iff there exists a smooth properties) of a system, which is exactly what we need for
fiber bundle π : M −→ (Rm )N , for a certain m ∈ N∗ , our transformations not to modify the ROA of the considered
such that every fiber is finite-dimensional with locally constant equilibrium point.
dimension, and for all ξ ∈ M (
R5 −→ R6
∇π(ξ) · F (ξ) = Fm (π(ξ)) (15) In the present case Φ : is not a Lie-
y 7−→ z
This allows us to define: Bäcklund isomorphism : indeed, F and H are Φ-related, but
• local coordinates ξ = (x, u) , where u = π(ξ) and x ∈ Φ obviously does not have a smooth inverse Ψ : R6 −→ R5 .
Rn , in which In fact, it is intuitive that one has to add some algebraic
m X constraints on z so as for Φ to be a licit change of variables:
∂ X (k+1) ∂
F (ξ) = f (x, u) + ui (16)
∂x i=1 ∂u(k)
k≥0 i
G(z) = 0. (20)
with f depending on a finite number of coordinates.
˙
• trajectories t 7−→ ξ(t) := (x(t), u(t)) such that ξ(t) = Here, the algebraic constraint is:
F (ξ(t)) i.e. ẋ(t) = f (x(t), u(t))
This way, one obtains a controlled differential system with
finite dimension. However, the definition of the state variables G(z) = z12 + z22 − 2z2 . (21)
and the control variables entirely depends on the choice of
π, which makes this definition fit also the non-controlled One can verify that Φ : R5 −→ {z ∈ R6 |G(z) = 0} is an
systems. In fact, the presence of a control does not depend endogenous transformation of the system. It is then possible to
on the system, but on the projection π one uses. apply the SOS approach to the differential algebraic equation
(DAE):
Then we can introduce the notion of Lie-B”acklund (
transformation. ż = H 0 (z)
(22)
G(z) = 0
Definition : Let (M, F ), (N , H) be two systems, Φ :
C∞
M −→ N , p ∈ M and q := Φ(p) ∈ N . Then, V. SOS ESTIMATION OF THE ROA
• if ξ is a trajectory of (M, F ) in a neighbourhood of p,
then ζ := Φ ◦ ξ stays in a neighbourhood of q and we In the expanding interior algorithm, the addition of the
have equality constraints G(z) = 0 only influences the definition
ζ̇(t) = ∇Φ(ξ(t)) · F (ξ(t)) (17) of the domain
which holds even in infinite dimension: everything de- ∆ = {z ∈ R6 |G(z) = 0 and V (z) ≤ 1}
pends only on a finite number of coordinates. and of
• Φ is an endogenous transformation iff, for any ξ in a
neighbourhood of p Pβ = {z ∈ R6 |p(z) ≤ β ; G(z) = 0}. (23)
∇Φ(ξ) · F (ξ) = H(Φ(ξ)) (18)
This leads, according to the P-satz, to the expanding interior
(we say that F and H are Φ-related at (p, q); then ζ̇ =
problem (7) enriched as
H(ζ)) and Φ has a smooth inverse Ψ (then, H and F are
automatically Ψ-related).
• Φ is a Lie-Bäcklund isomorphism iff we locally have
V − `1 − q1 G = s4 ∈ Σn (24a)
T Φ(span(F )) = span(H) − ((β − p)s1 + (V − 1)) − q2 G = s5 ∈ Σn (24b)
 
and Φ has a smooth inverse Ψ such that T Ψ(span(H)) = − (1 − V )s2 + V̇ s3 + `2 − q3 G = s6 ∈ Σn (24c)
span(F )2 . In other words, for ξ in a neighbourhood of p
and ζ in a neighbourhood of q, we should have
where q1,...,3 are free polynomial variables (qi G is therefore
(∇Φ(ξ))(R · F (ξ)) = R · H(Φ(ξ)) (19a) the definition of an element of the ideal I(G)). The expanding
(∇Ψ(ζ))(R · H(ζ)) = R · F (Ψ(ζ)) (19b) interior algorithm then returned the following estimation of the
equilibrium point’s ROA:
2 A distribution D :=span(V , . . . , V ) being intuitively defined on a man-
1 k
ifold M by D(x) =span(V1 (x), . . . , Vk (x)) ⊂ T M, where V1 , . . . , Vk are
Pβ = {z ∈ R6 |p(z) := kzk2 + z32 ≤ β = 0.08544} ⊂ ∆
vector fields on M. T Φ is a notation to denote the application ξ 7→ ∇Φ(ξ)·•. with
 V (z) = 2.010 z12 + 0.07823 z1 z2 + 3.1961 z1 z3
−2.244 z1 z4 − 0.02231 z1 z5 + 0.2172 z1 z6
+0.9483 z22 + 3.422 z2 z3 − 2.246 z2 z4
−0.003099 z2 z5 + 0.1913 z2 z6 + 22.92 z32 (25)
−0.07196 z3 z4 − 0.07616 z3 z5 + 2.998 z3 z6
+4.058 z42 − 0.0003899 z4 z5 − 0.1467 z4 z6
+0.004611 z52 + 0.008518 z5 z6 + 0.2425 z62 .
Two-dimensional projections of the resulting ROA in origi-
nal coordinates are plotted in Fig. 4 and 5 in red lines. These
estimations are quite large, especially when compared to the
(a)
exact ROA (blue lines in 5) numerically computed as explained
in the next section.
VI. ROA ESTIMATION VIA NUMERICAL TIME - DOMAIN
SIMULATIONS
We aim to validate the estimate of the ROA found by the
SOS approach, by testing, in simulation, the limits of stability
of the system in all the state-space directions. For this, the
system is systematically initialized at a starting point far from
the considered equilibrium point, and we check by simulation
if it goes back to equilibrium or not.
Since the system has 5 state variables which means a huge
number of combinations and because δ and ω are the most
important state variables for transient stability analysis, we
decide to make the test in a projection of the state space on
the plane (δ, ω), (Pm , ω) and (e0q , Ef d ). (b)
Figure 5 shows that the estimated ROA is inside the real
one (computed numerically by simulation) and this validates
the previous results. The arrows in the plots show that the real
ROA is larger in the direction of the arrows.
The same figure shows that the trajectories of the system
initialized in several points in the ROA converge to the
considered equilibrium point asymptotically, and thus the
computed ROA is compliant with the Lyapunov conditions
for asymptotic stability.
VII. C ONCLUSIONS
SOS approach and tools have been sucessfully used to
quantify transient stability of a SMIB system for which the
generator has been modeled in more detail as in previous
studies. Indeed, voltage dynamics and voltage and frequency (c)
regulations were taken into account in this formalism. First, Figure 4: (a) estimated ROA and the LF, together in a projec-
this provides more accuracy in estimation of the stability tion of the state space on the plane (e0q , Ef d ). (b) estimated
margin in terms of the ROA. Indeed, the estimated ROA ROA and the LF, together in a projection of the state space
is large enough compared to the exact ROA computed by on the plane (Pm , ω). (c) estimated ROA and the LF, together
simulation. Next, the Lyapunov approach is well suited for in a projection of the state-pace on the plane (δ, ω).
the control synthesis and this quantification can be further
exploited to build/tune regulators in order to maximize ROA.
As a matter of fact, in the SOS optimization one can next
include the regulators’ parameters as decision variables. Future • inclusion of the non linearities of the machine related to
work will focus on saturations of the actuating variables
• estimation with the full model (without the approximation • application to larger grids
(13)) • extension to the tuning of regulators’ parameters
 Td0 0 = 9.67 xd = 2.38 x0d = 0.336 xq = 1.21
100 H=3 r = 0.002 ωs = ωref = 1 R = 0.01
X = 1.185 Vs = 1 Ta = 1 Ka = 70
Vref = 1 Tg = 0.4 Kg = 0.5 Pref = 0.7
50 The considered equilibrium point (of (9) approximated with
(13)):
Efeqd = 2.459 ; e0eq
q = 1.070 ; δ
eq
= 1.539
eq eq
0 ω = ωref = 1 ; Pm = Pref = 0.7
R EFERENCES
[1] H. D. Chiang, Direct Methods for Stability Analysis of Electric Power
-50
-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5
Systems, Wiley 2011.
[2] M. Anghel, F. Milano, and A. Papachristodoulou, Algorithmic construc-
tion of Lyapunov functions for power system stability analysis Circuits
(a) and Systems I: Regular Papers, IEEE Transactions on, vol. 60, no. 9, pp.
25332546, Sept 2013.
[3] Z. W. Jarvis-Wloszek, Lyapunov based analysis and controller synthesis
for polynomial systems using sum-of-squares optimization, Ph.D. disser-
0.4
tation, University of California, Berkeley, CA, 2003.
[4] M. Fliess, J. Levine, P. Martin, and P. Rouchon. A Lie-Bäcklund
0.2
approach to equivalence and flatness of nonlinear systems, IEEE Trans-
actions on Automatic Control,44(5): 922-937, may 1999.
0 [5] S. Prajna, A. Papachristodoulou, and P.A. Parrilo. SOS-
TOOLS: Sum of squares optimization toolbox for Matlab,
-0.2 http://www.cds.caltech.edu/sostools, 2002.
[6] S. Prajna, A. Papachristodoulou, and P.A. Parrilo. Introducing sostools:
-0.4 A general purpose sum of squares programming solver, in Proceedings
of the Conference on Decision and Control, pp. 741746, 2002.
-0.6 [7] P.W. Sauer and M.A. Pai. Power System Dynamics and Stability. Prentice
Hall, 1998.
-0.8 [8] J. Bochnak, M. Coste, and M.-F. Roy. Geometrie Algebrique Reelle.
-8 -6 -4 -2 0 2 4 6 Springer, 1986.
[9] J.F. Sturm, Using sedumi 1.02, a matlab toolbox for optimizaton over
(b) symmetric cones, Optimization Methods and Software, 1112:pp. 625 653,
1999.
[10] Y. Nesterov, and A. Nemirovskii. Interior-Point Polynomial Algorithms
in Convex Programming, Studies in Applied and Numerical Mathematics,
0.3 SIAM, 1994.

0.2

0.1

-0.1

-0.2

-0.3

-1 -0.5 0 0.5 1

(c)
Figure 5: (a) real ROA in blue (the largest one) and the SOS
estimated one in red, together in a projection of the state-space
on the plane (e0q , Ef d ). (b) real ROA in blue (the largest one)
and the SOS estimated one in red, together in a projection of
the state space on the plane (Pm , ω). (c) real ROA in blue
(the largest one) and the SOS estimated one in red, together
in a projection of the state space on the plane (δ, ω).

VIII. A PPENDIX

Parameters of the test system