Inverse Stability Problem and Applications to
Renewables Integration

arXiv:1703.04491v5 [] 14 Oct 2017

Thanh Long Vu∗ , Hung Dinh Nguyen, Alexandre Megretski, Jean-Jacques Slotine and Konstantin Turitsyn

Abstract—In modern power systems, the operating point, at
which the demand and supply are balanced, may take different
values due to changes in loads and renewable generation levels.
Understanding the dynamics of stressed power systems with
a range of operating points would be essential to assuring
their reliable operation, and possibly allow higher integration
of renewable resources. This letter introduces a non-traditional
way to think about the stability assessment problem of power
systems. Instead of estimating the set of initial states leading to
a given operating condition, we characterize the set of operating
conditions that a power grid converges to from a given initial
state under changes in power injections and lines. We term
this problem as “inverse stability”, a problem which is rarely
addressed in the control and systems literature, and hence, poorly
understood. Exploiting quadratic approximations of the system’s
energy function, we introduce an estimate of the inverse stability
region. Also, we briefly describe three important applications
of the inverse stability notion: (i) robust stability assessment
of power systems w.r.t. different renewable generation levels,
(ii) stability-constrained optimal power flow (sOPF), and (iii)
stability-guaranteed corrective action design.
Index Terms—Power grids, renewables integration, transient
stability, inverse stability, emergency control, energy function

I. I NTRODUCTION

R

ENEWABLE generations, e.g., wind and solar, are increasingly installed into electric power grids to reduce
CO2 emission from the electricity generation sector. Yet,
their natural intermittency presents a major challenge to the
delivery of consistent power that is necessary for today’s grid
operation, in which generation must instantly meet load. Also,
the inherently low inertia of renewable generators limits the
grid’s controllability and makes it easy for the grid to lose its
stability. The existing power grids and management tools were
not designed to deal with these new challenges. Therefore,
new stability assessment and control design tools are needed
to adapt to the changes in architecture and dynamic behavior
expected in the future power grids.
Transient stability assessment of power system certifies that
the system state converges to a stable operating condition after
the system experiences large disturbances. Traditionally, this
task is handled by using either time domain simulation (e.g.,
[1]), or by utilizing the energy method (e.g., [2], [3]) and the
Lyapunov function method (e.g., [4]) to estimate the stability
region of a given equilibrium point (EP), i.e., the set of initial
states from which the system state will converge to that EP.
In modern renewable power grids, the operating point may
take different values under the real-time clearing of electricity
∗ Corresponding author. Email: longvu@mit.edu. All the authors are with
the Massachusetts Institute of Technology, Cambridge, MA 02139, USA.

markets, intermittent renewable generations, changing loads,
and external disturbances. Dealing with the situation when the
EP can change over a wide range makes the transient stability
assessment even more technically difficult and computationally cumbersome.
In this letter, rather than considering the classical stability
assessment problem, we formulate the inverse stability assessment problem. This problem concerns with estimating the
region around a given initial state δ0 , called “inverse stability
region” A(δ0 ), so that whenever the power injections or power
lines change and lead to an EP in A(δ0 ), the system state will
converge from δ0 to that EP. Indeed, the convergence from δ0
to an EP is guaranteed when the system’s energy function is
bounded under some threshold [2], [3]. In [5], we observed
that if the EP is in the interior of the set P characterized by
phasor angular differences smaller than π/2, then the nonlinear
power flows can be strictly bounded by linear functions of
angular differences. Exploiting this observation, we show that
the energy function of power system can be approximated by
quadratic functions of the EP and the system state, and from
which we obtain an estimate of the inverse stability region.
The remarkable advantage of the inverse stability certificate
is making it possible to exploit the change in EP to achieve
useful dynamical properties. We will briefly discuss three
applications of this certificate, which are of importance to the
integration of large-scale renewable resources:
Robust stability assessment: For a typical power system
composed of several components and integrated with
different levels of renewable generations, there are many
contingencies that need to be reassessed on a regular
basis. Most of these contingencies correspond to failures of relatively small and insignificant components,
so the post-fault dynamics is probably transiently stable.
Therefore, most of the computational effort is spent on
the analysis of non-critical scenarios. This computational
burden could be greatly alleviated by a robust transient
stability assessment toolbox that could certify the system’s stability w.r.t. a broad range of uncertainties. In
this letter, we show that the inverse stability certificate
can be employed to assess the transient stability of power
systems for various levels of power injections.
Stability-constrained OPF: Under large disturbances, a
power system with an operating condition derived by
solving the conventional OPF problem may not survive.
It is therefore desirable to design operating conditions
so that the system can withstand large disturbances.
This can be carried out by incorporating the transient
stability constraint into OPF together with the normal

voltage and thermal constraints. Though this problem was
discussed in the literature (e.g., [6]), there is no way
to precisely formulate and solve the stability-constrained
OPF problem because transient stability is a dynamic
concept and differential equations are involved in the
stability constraint. Fortunately, the inverse stability certificate allows for a natural incorporation of the stability
constraint into the OPF problem as a static constraint of
placing the EP in a given set.
Stability-guaranteed corrective actions: Traditional
protection strategies focus on the safety of individual
components, and the level of coordination among component protection systems is far from perfect. Also,
they do not take full advantage of the new flexible
and fast electronics resources available in modern power
systems, and largely rely on customer-harmful actions
like load shedding. These considerations motivated us to
coordinate widespread flexible electronics resources as
a system-level customer-friendly corrective action with
guaranteed stability [7]. This letter presents a unconventional control way in which we relocate the operating
point, by appropriately redispatching power injections,
to attract the emergency state and stabilize the power
systems under emergency situations.
II. I NVERSE STABILITY PROBLEM OF POWER SYSTEMS
In this letter, we utilize the structure-preserving model to
describe the power system dynamics [8]. This model naturally
incorporates the dynamics of the generators’ rotor angle and
the response of load power output to frequency deviation.
Mathematically, the grid is described by an undirected graph
A(N , E), where N = {1, 2, . . . , |N |} is the set of buses and
E ⊆ N × N is the set of transmission lines {k, j}, k, j ∈ N .
Here, |A| denotes the number of elements of set A. The sets
of generator buses and load buses are denoted by G and L.
We assume that the grid is lossless with constant voltage
magnitudes Vk , k ∈ N , and the reactive powers are ignored.
Then, the grid’s dynamics is described by [8]:
X
mk δ¨k + dk δ˙k +
akj sin(δk − δj ) =Pk , k ∈ G, (1a)
{k,j}∈E

dk δ˙k +

X

akj sin(δk − δj ) =Pk , k ∈ L, (1b)

{k,j}∈E

where equation (1a) applies at the dynamics of generator
buses and equation (1b) applies at the dynamics of load
buses. Here akj = Vk Vj Bkj , where Bkj is the (normalized)
susceptance of the transmission line {k, j} connecting the
k th bus and j th bus. Nk is the set of neighboring buses
of the k th bus (see [9] for more details). Let δ(t) =
[δ1 (t) ... δ|N | (t) δ̇1 (t) ... δ̇|N | (t)]> be the state of the system
(1) at time t (for simplicity, we will denote the system state
by δ). Note that Eqs. (1) are invariant under any uniform
shift of the angles δk → δk + c. However, the state δ
can be unambiguously characterized by the angle differences
δkj = δk − δj and the frequencies δ̇k .
Normally, a power grid operates at an operating condition
of the pre-fault dynamics. Under the fault, the system evolves

according to the fault-on dynamics. After some time period,
the fault is cleared or self-clears, and the system is at the socalled fault-cleared state δ0 (the fault-cleared state is usually
estimated by simulating the fault-on dynamics, and hence, is
assumed to be known). Then, the power system experiences
the so-called post-fault dynamics. The transient stability assessment certifies whether the post-fault state converges from
δ0 to a stable EP δ ∗ . Mathematically, the operating condition
∗
>
δ ∗ = [δ1∗ ... δ|N
| 0 ... 0] is a solution of the power-flow like
equations:
X
akj sin δkj =Pk , k ∈ N .
(2)
{k,j}∈E

With renewable generations or under power redispatching,
the power injections Pk take different values. Also, the couplings akj can be changed by using the FACTS devices. Assume akj ≤ akj ≤ ākj . In those situations, the resulting EP δ ∗
also takes different values. Therefore, we want to characterize
the region of EPs so that the post-fault state always converges
from a given initial state δ0 to the EP whenever the EP is
in this region. Though the EP can take different values, it
is assumed to be fixed in each transient stability assessment
because the power injections and couplings can be assumed to
be unchanged in the very fast time scale of transient dynamics
(i.e., 1 to 10 seconds). We consider the following problem:
• Inverse (Asymptotic) Stability Problem: Consider a
given initial state δ0 . Assume that power injections and
the line susceptances can take different values. Estimate
the region of stable EPs so that the state of the system
(1) always converges from δ0 to the EP in this region.
This problem will be addressed with the inverse stability
certificate to be presented in the next section.
III. E NERGY FUNCTION AND INVERSE STABILITY
CERTIFICATE

A. Stability assessment by using energy function
Before introducing the inverse stability certificate addressing the inverse stability problem in the previous section, we
present a normal stability certificate for system with the fixed
power injections and line parameters. For the power system
described by Eqs. (1), consider the energy function:
X Z δkj
X mk δ̇ 2
∗
k
E(δ, δ ∗ ) =
+
akj (sin ξ − sin δkj
)dξ
∗
2
δkj
{k,j}∈E

k∈G

(3)
Then, along every trajectory of (1), we have
X
X
∗
Ė(δ, δ ∗ ) =
mk δ̇k δ̈k +
akj (sin δkj − sin δkj
)δ̇kj
{k,j}∈E

k∈G

=

X

δ̇k (Pk − dk δ˙k −

X

δ̇k (Pk − dk δ˙k −

X

X

akj sin(δk − δj ))

{k,j}∈E

k∈L

+

akj sin(δk − δj ))

{k,j}∈E

k∈G

+

X

∗
akj (sin δkj − sin δkj
)(δ̇k − δ˙j )

{k,j}∈E

=−

X
k∈N

dk (δ˙k )2 ≤ 0,

(4)

in which the last equation is obtained from (2). Hence,
E(δ, δ ∗ ) is always decreasing along every trajectory of (1).
Consider the set P defined by |δkj | ≤ π/2, {k, j} ∈ E,
and the set Φ = {δ ∈ P : E(δ, δ ∗ ) < Emin (δ ∗ )}, where
Emin (δ ∗ ) = minδ∈∂P E(δ, δ ∗ ) and ∂P is the boundary of
P. Φ is invariant w.r.t. (1), and bounded as the state δ is
characterized by the angle differences and the frequencies.
Though Φ is not closed, the decrease of E(δ, δ ∗ ) inside Φ
assures the limit set to be inside Φ. As such, we can apply
the LaSalle’s Invariance Principle and use a proof similar to
that of Theorem 1 in [5] to show that, if δ0 is inside Φ then
the system state will only evolve inside this set and eventually
converge to δ ∗ . So, to check if the system state converges from
δ0 ∈ P to δ ∗ , we only need to check if E(δ0 , δ ∗ ) < Emin (δ ∗ ).
B. Inverse stability certificate
1
0.8

R(δ0)

M

0.6
0.4

δ∗
y

0.2

δ0
B(δ0 )

0
−0.2
−0.4

Λ

−0.6
−0.8

P

−1
−1

−0.5

0

0.5

1

1.5

2

2.5

Fig. 1. For a power system with a xgiven initial state δ0 , if the EP δ ∗ is
inside the set A(δ0 ) = Λ ∩ B(δ0 ) surrounding δ0 then the system state will
converge from δ0 to the EP δ ∗ since E(δ0 , δ ∗ ) < Emin .

For a given initial state δ0 ∈ P, we will construct a region
surrounding it so that whenever the operating condition δ ∗ is
in this region then E(δ0 , δ ∗ ) < Emin (δ ∗ ). Hence, the grid state
will converge from δ0 to δ ∗ according to the stability certificate
in Section III-A. Indeed, we establish quadratic bounds of the
energy function for every δ ∗ in the set Λ defined by inequali1 − sin λ
ties |δkj | ≤ λ < π/2, ∀{k, j} ∈ E. Let g =
> 0. In
π/2 − λ
∗
[5], we observed that for δ ∈ Λ, ξ ∈ P,

∗ 2
X mk δ̇ 2
X
(δkj − δkj
)
k
E(δ, δ ) ≥
+g
akj
,
2
2
∗

(5)

{k,j}∈E

∗ 2
X mk δ̇ 2
X
(δkj − δkj
)
k
+
akj
.
2
2

(6)

{k,j}∈E

k∈G

Define the following functions
D(δ, δ ∗ ) = g

X
{k,j}∈E

F (δ, δ ∗ ) =

akj

∗ 2
(δkj − δkj
)
,
2

X
X mk δ̇ 2
(δkj −
k
+
ākj
2
2

k∈G

(9)

The following is our main result regarding inverse stability
of power system, as illustrated in Fig. 1.
Theorem 1: Consider a given initial state δ0 inside the set
P. Assume that the EP of the system takes different values in
the set A(δ0 ) = Λ ∩ B(δ0 ), where the set B(δ0 ) is defined as
in (9). Then, the system state always converges from the given
initial state δ0 to the EP.
Proof: See Appendix VI.

Remark 1: In this paper, we limit the grid to be described
by the simplified model (1) which captures the dynamics of
the generators’ rotor angle and the response of load power
output to frequency deviation. More realistic models should
take into account voltage variations, reactive powers, dynamics
of the rotor flux, and controllers (e.g., droop controls and
power system stabilizers). It should be noted that the results
in this paper is extendable to more realistic models. The
reason is that all the key results in this paper rely on the
analysis of the system’s energy function, in which we combine
the energy function-based transient stability analysis with the
quadratic bounds of the energy function. In high-order models,
the energy function is more complicated, yet it still can be
bounded by quadratic functions. Hence, we can combine the
energy function-based transient stability analysis of the more
realistic models in [2]–[4] with the quadratic bounds of the
energy function to extend the results in this paper to these
higher-order models. This approach may also extend to other
higher-order models in the port-Hamiltonian formulation (e.g.,
[10], [11]) by applying the appropriate approximation on the
Lyapunov functions established for these models.

A. Robust stability assessment

Hence, for all δ ∗ ∈ Λ, δ ∈ P, we have

E(δ, δ ∗ ) ≤

B(δ0 ) = {δ : F (δ0 , δ) ≤ R(δ0 )/4}.

IV. A PPLICATIONS OF I NVERSE S TABILITY C ERTIFICATE

∗ 2
∗
∗
∗ 2
g(ξ − δkj
) ≤ (ξ − δkj
)(sin ξ − sin δkj
) ≤ (ξ − δkj
) .

k∈G

For a given initial state δ0 inside the set P, we calculate the
“distance” from this initial state to the boundary of the set P :
R(δ0 ) = minδ∈∂P D(δ0 , δ). Let B(δ0 ) be the neighborhood
of δ0 defined by

The robust transient stability problem that we consider involves situations where there is uncertainty in power injections
Pk , e.g., due to intermittent renewable generations. Formally,
for a given fault-cleared state δ0 , we need to certify the
transient stability of the post-fault dynamics described by (1)
with respect to fluctuations of the power injections, which
consequently lead to different values of the post-fault EP δ ∗
as a solution of the power flow equations (2). Therefore, we
consider the following robust stability problem [5]:
•

∗ 2
δkj
)

(7)
.

{k,j}∈E

Using (5) and (6), we can bound the energy function as
D(δ, δ ∗ ) ≤ E(δ, δ ∗ ) ≤ F (δ, δ ∗ ), ∀δ ∈ P, δ ∗ ∈ Λ.

(8)

Robust stability assessment: Given a fault-cleared state
δ0 , certify the transient stability of (1) w.r.t. a set of stable
EPs δ ∗ resulted from different levels of power injections.

Utilizing the inverse stability certificate, we can assure
robust stability of renewable power systems whenever the
resulting EP is inside the set A(δ0 ) = Λ ∩ B(δ0 ). To check
that the EP is in the set Λ, we can apply the criterion from

[12], which states that the EP will be in the set Λ if the power
injections p = [P1 , ..., P|N | ]T satisfy
kL† pkE,∞ ≤ sin λ,

(10)

where L† is the pseudoinverse of the network Laplacian
matrix and the norm kxkE,∞ is defined by kxkE,∞ =
max{i,j}∈E |x(i) − x(j)|. On the other hand, some similar
sufficient condition could be developed so that we can verify
that the EP is in the set B(δ0 ) by checking the power
injections. This will help us certify robust stability of the
system by only checking the power injections.
B. Stability-constrained OPF
Stability-constrained OPF problem concerns with determining the optimal operating condition with respect to the voltage
and thermal constraints, as well as the stability constraint.
While the voltage and thermal constraints are well modeled
via algebraic equations or inequalities, it is still an open
question as to how to include the stability constraint into
OPF formulation since stability is a dynamic concept and
differential equations are involved [6].
Mathematically, a standard OPF problem is usually stated
as follows (refer to [6] for more detailed formulation):
s.t.

min c(P )

(11)

P (V, δ) = P

(12)

Q(V, δ) = Q

(13)

V ≤ V ≤ V̄

(14)

S ≤ |S(V, δ)| ≤ S̄

(15)

where c(P ) is a quadratic cost function, the decision variables P are typically the generator scheduled electrical power
outputs, the equality constraints (12)-(13) stand for the power
flow equations, and the inequality constraints (14)-(15) stand
for the voltage and thermal limits of branch flows through
transmission lines and transformers. Assume that the stability
constraint is to make sure that the system state will converge
from a given fault-cleared state δ0 to the designed operating
condition, and that the reactive power is negligible. With
the inverse stability certificate, the stability constraint can be
relaxed and formulated as δ ∈ A(δ0 ). Basically, the inverse
stability certificate transforms the dynamic problem of stability
into a static problem of placing the prospective EP into a set.
In summary, we obtain a relaxation of the stability-constrained
OPF problem as follows:
s.t.

C. Emergency control design

min c(P )

(16)

P (V, δ) = P

(17)

V ≤ V ≤ V̄

(18)

S ≤ |S(V, δ)| ≤ S̄

(19)

δ ∈ A(δ0 ).

(20)

Solution of this optimization problem in an optimal operating condition at which the cost function is minimized
and the voltage/thermal constraints are respected. Furthermore,
the stability constraint is guaranteed by the inverse stability
certificate, and the system state is ensured to converge from
the fault-cleared δ0 to the operating condition.

δ0

δ ∗1

δ ∗2 δ ∗N

δ ∗desired

Fig. 2. Power dispatching to relocate the stable EPs δi∗ so that the faultcleared state, which is possibly unstable if there is no controls, is driven
∗
through a sequence of EPs back to the desired EP δdesired
. The placement
of these EPs is determined by applying the inverse stability certificate.

Another application of the inverse stability certificate,
that will be detailed in this section, is designing stabilityguaranteed corrective actions that can drive the post-fault
dynamics to a desired stability regime. As illustrated in Fig.
2, for a given fault-cleared state δ0 , by applying the inverse
stability certificate, we can appropriately dispatch the power
injections Pk to relocate the EP of the system so that the postfault dynamics can be attracted from the fault-cleared state
∗
δ0 through a sequence of EPs δ1∗ , ..., δN
to the desired EP
∗
δdesired
. In other words, we subsequently redispatch the power
injections so that the system state converges from δ0 to δ1∗ ,
∗
∗
to δdesired
. This
and then, from δ1∗ to δ2∗ , and finally, from δN
type of corrective actions reduces the need for prolonged load
shedding and full state measurement. Also, this control method
is unconventional where the operating point is relocated as
desired, instead of being fixed as in the classical controls.
Mathematically, we consider the following problem:
• Emergency Control Design: Given a fault-cleared state
∗
δ0 and a desired stable EP δdesired
, determine the
feasible dispatching of power injections Pk to relocate
the EPs so that the post-fault dynamics is driven from
the fault-cleared state δ0 through the set of designed EPs
∗
to the desired EP δdesired
.
To solve this problem, we can design the first EP δ1∗ by
minimizing kL† pkE,∞ over all possible power injections. The
optimum power injection will result in an EP which is most
far away from the stability margin |δkj | = π/2, and hence,
the stability region of the first EP δ1∗ probably contains the
given fault-cleared state. To design the sequence of EPs, in
each step, we carry out the following tasks:
∗
∗
• Calculate the distance R(δi−1 ) from δi−1 to the boundary
∗
∗
of the set P, i.e., R(δi−1 ) = minδ∈∂P D(δ, δi−1
). Noting
∗
that minimization of D(δ, δi−1 ) over the boundary of the
set P is a convex problem with a quadratic objective
function and linear constraints. Hence, we can quickly
∗
obtain R(δi−1
).
∗
∗
• Determine the set B(δi−1 ) and the set A(δi−1 ).
∗
• The next EP δi will be chosen as the intersection of
∗
the boundary of the set A(δi−1
) and the line segment
∗
∗
connecting δi−1 and δdesired .
• The power injections Pk that we have to redispatch will
P
(i)
be determined by Pk = j∈Nk akj sin δi∗kj for all k.

0.05

G

δ1 (t)
δ2 (t)
δ3 (t)
δ4 (t)
δ5 (t)
δ6 (t)
δ7 (t)
δ8 (t)
δ9 (t)

G
0

2

8

7

9

3

Angles (rad)

−0.05

5

6

4

−0.1
−0.15
−0.2
−0.25

1
G

−0.3

Fig. 3. 3-generator 9-bus system with frequency-dependent dynamic loads.
2

0

1

2

3

4

5

6

time (s)

Fig. 5. Stable dynamics with power injection control: Convergence of buses
angles from the fault-cleared state to δ1∗ in the post-fault dynamics

1

δ 1(t)
δ 2(t)
δ 3(t)
δ 4(t)
δ 5(t)
δ 6(t)
δ 7(t)
δ 8(t)
δ 9(t)

−2
−3
−4
−5
−6
−7

0

0.5

1

1.5

2

2.5

3

3.5

7

D 2 (t)

6

Distance to δ2∗

Angles (rad)

0
−1

5
4
3
2
1
0
0

4

1

2

3

4

5

This power dispatch will place the new EP at δi∗ which
is in the inverse stability region of the previous EP
∗
δi−1
. Therefore, the controlled post-fault dynamics will
∗
converge from δi−1
to δi∗ .
This procedure strictly reduces the distance from EP to
∗
(it can be proved that there exists a constant d > 0 so
δdesired
that such distance reduces at least d in each step). Hence, after
∗
will be sufficiently near the desired EP
some steps, the EP δN
∗
δdesired so that the convergence of the system state to the
∗
desired EP δdesired
will be guaranteed.
Node
1
2
3
4
5
6
7
8
9

V (p.u.)
1.0284
1.0085
0.9522
1.0627
1.0707
1.0749
1.0490
1.0579
1.0521
TABLE I

Pk (p.u.)
3.6466
4.5735
3.8173
-3.4771
-3.5798
-3.3112
-0.5639
-0.5000
-0.6054

B US VOLTAGES , MECHANICAL INPUTS , AND STATIC LOADS .

To illustrate that this control works well in stabilizing some
possibly unstable fault-cleared state δ0 , we consider the 3machine 9-bus system with 3 generator buses and 6 frequencydependent load buses as in Fig. 3. The susceptances of the
transmission lines are as follows: B14 = 17.3611p.u., B27 =
16.0000p.u., B39 = 17.0648p.u., B45 = 11.7647p.u., B57 =
6.2112p.u., B64 = 10.8696p.u., B78 = 13.8889p.u., B89 =
9.9206p.u., B96 = 5.8824p.u. The parameters for generators are: m1 = 0.1254, m2 = 0.034, m3 = 0.016, d1 =
0.0627, d2 = 0.017, d3 = 0.008. For simplicity, we take dk =
0.05, k = 4 . . . , 9. Assume that the fault trips the line between

7

8

9

10

Fig. 6. Effect of power dispatching control: the convergence of the distance
D2 (t) to 0. Here, the Euclid distance
qP D2 (t) between a state δ and the second
9
∗
2
EP δ2∗ is defined as D2 (t) =
i=2 (δi1 (t) − δ2i1 ) .
0.5

Ddesired(t)
0.4
∗
Distance to δdesired

Fig. 4. Unstable post-fault dynamics when there is no controls: |δ45 | and
|δ57 | evolve to 2π, triggering the tripping of lines {4, 5} and {5, 7}.

6

time (s)

time (s)

0.3

0.2

0.1

0

0

1

2

3

4

5

6

7

8

9

10

time (s)

Fig. 7. Autonomous dynamics when we switch the power injections to the
desired values: the convergence of the distance Ddesired (t) to 0. Here,
∗
the distance Ddesired (t) between
qP a state δ and the desired EP δdesired
9
∗
2
is defined as Ddesired (t) =
i=2 (δi1 (t) − δdesired ) .
i1

buses 5 and 7, and make the power injections to fluctuate.
When the fault is cleared this line is re-closed. We also assume
the fluctuation of the generation (probably due to renewables)
and load so that the voltages Vk and power injections Pk of the
∗
post-fault dynamics are given in Tab. I. The stable EP δdesired
is calculated as [−0.1629 0.4416 0.3623 −0.3563 −0.3608 −
0.3651 0.1680 0.1362 0.1371]> . However, the fault-cleared
state, with angles [0.025 − 0.023 0.041 0.012 − 2.917 −
0.004 0.907 0.021 0.023]> and generators angular velocity
[−0.016 − 0.021 0.014]> , is outside the set P. It can be seen
from Fig. 4 that the uncontrolled post-fault dynamics is not
stable since |δ45 | and |δ57 | quickly evolve from initial values
to 2π, which will activate the protective devices to trip the
lines.
Using CVX software [13] to minimize kL† pkE,∞ , we
obtain the new power injections at buses 1-6 as follows: P1 = 0.5890, P2 = 0.5930, P3 = 0.5989, P4 =
−0.0333, P5 = −0.0617, and P6 = −0.0165. Accord-

ingly, the minimum value of kL† pkE,∞ = 0.0350 <
sin(π/89). Hence, the first EP obtained from equation (2)
will be in the set defined by the inequalities |δkj | ≤
π/89, ∀{k, j} ∈ E, and can be approximated by δ1∗ ≈
L† p = [0.0581 0.0042 0.0070 0.0271 0.0042 0.0070 −
0.0308 − 0.0486 − 0.0281]> . The simulation results confirm that the post-fault dynamics is made stable by applying
the optimum power injection control, as showed in Fig. 5.
Using the above procedure, after one step, we can find that
∗
δ2∗ = 0.9259δdesired
+ 0.0741δ1∗ is the intersection of the set
∗
A(δ1∗ ) and the line segment connecting δ1∗ and δdesired
. This
EP is inside the inverse stability region of δ1∗ , and hence the
system state will converge from δ1∗ to δ2∗ when we do the
(2)
power dispatching Pk corresponding to δ2∗ . On the other
∗
hand, δ2∗ is very near the desired EP δdesired
and it is easy to
∗
check that δdesired is in the inverse stability region of δ2∗ , and
thus the system state will converge from δ2∗ to the desired
∗
EP δdesired
. Such convergence of the controlled post-fault
dynamics is confirmed in Figs. 6-7.

Note that D(δ0 , M ) ≥ R(δ0 ) as R(δ0 ) is the distance from δ0
to the boundary of the set P. This, together with (21), leads
to E(M, δ ∗ ) + E(δ0 , δ ∗ ) ≥ R(δ0 )/2. From (8) and (9), we
have E(δ0 , δ ∗ ) ≤ F (δ0 , δ ∗ ) < R(δ0 )/4. Hence,

V. C ONCLUSIONS
Electric power grids possess rich dynamical behaviours,
e.g., nonlinear interaction, prohibition of global stability,
and exhibition of significant uncertainties, that challenge the
maintenance of their reliable operation and pose interesting
questions to control and power communities. This letter characterized a surprising property termed as “inverse stability”,
which was rarely investigated and poorly understood (though
some related inverse problems were addressed in [14]). This
new notion could change the way we think about the stability
assessment problem. Instead of estimating the set of initial
states leading to a given operating condition, we characterized
the set of operating conditions that a power grid converges to
from a given initial state under changes in power injections
and lines. In addition, we briefly described three applications
of the inverse stability certificate: (i) assessing the stability of
renewable power systems, (ii) solving the stability-constrained
OPF problem, and (iii) designing power dispatching remedial
actions to recover the transient stability of power systems.
Remarkably, we showed that robust stability due to the fluctuation of renewable generations can be effectively assessed,
and that the stability constraint can be incorporated as a static
constraint into the conventional OPF. We also illustrated a
unconventional control method, in which we appropriately
relocate the operating point to attract a given fault-cleared state
that originally leads to an unstable dynamics.

R EFERENCES

VI. A PPENDIX : P ROOF OF T HEOREM 1
For each EP δ ∗ ∈ A(δ0 ), let M be the point on the
boundary of the set P so that E(M, δ ∗ ) = Emin (δ ∗ ) =
minδ∈∂P E(δ, δ ∗ ), as showed in Fig. 1. From (8), we have
E(M, δ ∗ ) + E(δ0 , δ ∗ ) ≥ D(M, δ ∗ ) + D(δ0 , δ ∗ )
∗ 2
∗ 2
X
(δMkj − δkj
) + (δ0kj − δkj
)
=g
akj
2
{k,j}∈E

≥g

X
{k,j}∈E

akj

(δMkj − δ0kj )2
D(δ0 , M )
=
.
4
2

(21)

E(M, δ ∗ ) > R(δ0 )/4 > E(δ0 , δ ∗ ).

(22)

Therefore, for any δ ∗ ∈ A(δ0 ), we have E(δ0 , δ ∗ ) <
Emin (δ ∗ ). By applying the stability analysis in Section III-A,
we conclude that the initial state δ0 must be inside the stability
region of the EP δ ∗ and the system state will converge from
the initial state δ0 to the EP δ ∗ .

VII. ACKNOWLEDGEMENTS
This work was supported by the MIT/Skoltech, Ministry
of Education and Science of Russian Federation (Grant no.
14.615.21.0001.), and NSF under Contracts 1508666 and
1550015.

[1] I. Nagel, L. Fabre, M. Pastre, F. Krummenacher, R. Cherkaoui, and
M. Kayal, “High-Speed Power System Transient Stability Simulation
Using Highly Dedicated Hardware,” Power Systems, IEEE Transactions
on, vol. 28, no. 4, pp. 4218–4227, 2013.
[2] A. Pai, Energy Function Analysis for Power System Stability, ser.
Power Electronics and Power Systems. Springer US, 2012. [Online].
Available: https://books.google.com/books?id=1HDgBwAAQBAJ
[3] H.-D. Chiang, Direct Methods for Stability Analysis of Electric Power
Systems, ser. Theoretical Foundation, BCU Methodologies, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Mar. 2011.
[4] R. Davy and I. A. Hiskens, “Lyapunov functions for multi-machine
power systems with dynamic loads,” Circuits and Systems I: Fundamental Theory and Applications, IEEE Transactions on, vol. 44, no. 9,
pp. 796–812, 1997.
[5] T. L. Vu and K. Turitsyn, “A framework for robust assessment of power
grid stability and resiliency,” IEEE Transactions on Automatic Control,
vol. 62, no. 3, pp. 1165–1177, March 2017.
[6] A. Pizano-Martianez and C. R. Fuerte-Esquivel and D. Ruiz-Vega,
“Global Transient Stability-Constrained Optimal Power Flow Using an
OMIB Reference Trajectory,” IEEE Transactions on Power Systems, vol.
25, no. 1, pp. 392–403, Feb 2010.
[7] T. L. Vu, S. Chatzivasileiadis, H. D. Chiang, and K. Turitsyn, “Structural
emergency control paradigm,” IEEE Journal on Emerging and Selected
Topics in Circuits and Systems, vol. 7, no. 3, pp. 371–382, Sept 2017.
[8] A. R. Bergen and D. J. Hill, “A structure preserving model for power
system stability analysis,” Power Apparatus and Systems, IEEE Transactions on, no. 1, pp. 25–35, 1981.
[9] T. L. Vu, S. M. A. Araifi, M. S. E. Moursi, and K. Turitsyn, “Toward
simulation-free estimation of critical clearing time,” IEEE Trans. on
Power Systems, vol. 31, no. 6, pp. 4722–4731, Nov 2016.
[10] S. Fiaz and D. Zonetti and R. Ortega and J.M.A. Scherpen and A.J. van
der Schaft, “A port-Hamiltonian approach to power network modeling
and analysis,” European Journal of Control, vol. 19, no. 6, pp. 477 –
485, 2013.
[11] S. Y. Caliskan and P. Tabuada, “Compositional Transient Stability
Analysis of Multimachine Power Networks,” IEEE Transactions on
Control of Network Systems, vol. 1, no. 1, pp. 4–14, March 2014.
[12] F. Dorfler, M. Chertkov, and F. Bullo, “Synchronization in complex
oscillator networks and smart grids,” Proceedings of the National
Academy of Sciences, vol. 110, no. 6, pp. 2005–2010, 2013.
[13] Grant, Michael C. and Boyd, Stephen P. and Ye, Yinu, “CVX: Matlab
software for disciplined convex programming (web page and software),”
Available at http://cvxr.com/cvx.
[14] I. A. Hiskens, “Power system modeling for inverse problems,” IEEE
Transactions on Circuits and Systems I: Regular Papers, vol. 51, no. 3,
pp. 539–551, March 2004.

