1

Adaptive Nonlinear Model Reduction for Fast
Power System Simulation
Denis Osipov, Student Member, IEEE, and Kai Sun, Senior Member, IEEE


Abstract--The paper proposes a new adaptive approach to
power system model reduction for fast and accurate time-domain
simulation. This new approach is a compromise between linear
model reduction for faster simulation and nonlinear model
reduction for better accuracy. During the simulation period, the
approach adaptively switches among detailed and linearly or
nonlinearly reduced models based on variations of the system
state: it employs unreduced models for the fault-on period, uses
weighted column norms of the admittance matrix to decide which
functions to be linearized in power system differential-algebraic
equations for large changes of the state, and adopts a linearly
reduced model for small changes of the state. Two versions of the
adaptive model reduction approach are introduced. The first
version uses traditional power system partitioning where the
model reduction is applied to a defined large external area in a
power system and the other area defined as the study area keeps
full detailed models. The second version applies the adaptive
model reduction to the whole system. The paper also conducts
comprehensive case studies comparing simulation results using
the proposed adaptively reduced models with the linearly
reduced model on the Northeast Power Coordinating Council
140-bus 48-machine system.
Index Terms--Linear model reduction, nonlinear model
reduction, power system partitioning, power system simulation.

P

I. INTRODUCTION

OWER system simulation is very important for grid
operations and planning at electricity utilities. It can assess
dynamic security under a certain operating condition of the
studied power system following a given contingency such as
loss of a transmission line or generator unit. Essentially,
power system simulation is to obtain a time-series trajectory
of the system state for a specified simulation window by
solving the initial value problem of a set of nonlinear
differential-algebraic equations co-determined by the
mathematical model of the whole system, the operating
condition and the contingency. Nowadays, the fast growth in
electricity demand and a relatively slow construction of new
transmission infrastructure are pushing power transmission
systems to be operated closer to their stability limits, and
motivating the transition of power system simulation from the
offline planning stage to the real-time operation environment.
This work was supported by NSF grant ECCS-1610025 and the CURENT
Engineering Research Center.
D. Osipov and K. Sun are with the Department of Electrical Engineering
and Computer Science, University of Tennessee, Knoxville, TN 37996 USA
(email: dosipov@vols.utk.edu, kaisun@utk.edu).

One way to increase the speed of simulation of a complex
power grid is to conduct network partitioning and then model
reduction. For example, a traditional approach defines a study
area, which is an important small part of the system for
dynamic security assessment, considering all the details, and
reduces the model of the rest of the system, i.e. the external
area. In practice, an additional buffer zone is often defined in
between with moderate reduction to connect the study and the
external areas [1]. The methods for model reduction on the
external area can be divided into two main groups: the ones
that preserve the structure of a power system and the ones that
use mathematical transformations from original states to
nonphysical states that are subsequently reduced.
The most widely used methods from the first group are
coherency-based methods [2]-[4], which were originally
developed for power system model reduction and conduct the
following steps: coherency identification, aggregation of
coherent generators, and network reduction. After the first
step, the generators that oscillate together following a
disturbance are included into one group. The groups of
coherent generators are then aggregated into individual
equivalent generators connected with each other by equivalent
branches and with the study area by a reduced system
network. This creates a unique boundary between the external
area and the study area and does not allow arbitrary division
between areas. In addition, if the topology of the original
system changes, it can affect the coherency and consequently
the split between the study area and the generator grouping of
the external area. This can change the boundary between the
study area and the external area. Thus, the grouping of
generators based on coherency has an inherent limitation on
the way a system can be partitioned.
The second group of methods does not have that limitation
as the states are transformed into a new state space. Thus, the
system can be partitioned in any way. These methods came
from the control field of engineering. In the most used
methods from this category, the external area model is first
linearized and then reduced using balanced truncation [5]-[7]
or Krylov subspace methods [8]-[10]. The linearization of the
model gives acceptably accurate results when concerned
disturbances are small but once the size of the disturbance
increases the linearized model cannot guarantee an accurate
representation of the original part of the system. To improve
the accuracy of large-disturbance simulation, nonlinear model
reduction methods can be used [11]-[15]. However, as shown
in [16], application of nonlinear model reduction cannot give

2

substantial computational time decrease as compared to the
original model. In addition, some of the methods require
training simulation data to create a reduced model, which
cannot guarantee adequate performance during all possible
disturbances. If the disturbance is very different from that with
the training set, the model reduction error can substantially
increase [17].
Compared to the existing work this paper proposes a new
adaptive model reduction approach, which is a compromise
between linear model reduction and nonlinear model
reduction in terms of accuracy and speed of time-domain
simulation using the reduced model. A comprehensive study is
also presented to compare this adaptive model reduction
approach with the linear model reduction approach. During
the simulation period, the approach adaptively switches
among detailed and linearly or nonlinearly reduced models
based on variations of the system state: it employs unreduced
models for the fault-on period, uses weighted column norms
of the admittance matrix to decide which functions to be
linearized in the power system model for large changes of the
state, and adopts a linearly reduced model for small changes
of the state.
The version of the adaptive approach described above uses
traditional topological power system partitioning with the
study area and the external area. This partitioning creates an
additional error that can affect the performance of the model
reduction. To address the partitioning error the second version
of the adaptive approach is proposed where the model
reduction is performed to the whole system.
Thus, the main contributions of this paper are in the
following aspects: 1) linearization of only certain functions
and the selection criterion of the functions to be linearized; 2)
adaptive switching among detailed and linearly or nonlinearly
reduced models and the switching criterion; 3) application of
power system model reduction to the system without
partitioning.
The remainder of this paper is organized as follows. Section
II describes the model and the approach used in power system
partitioning, presents the linear model reduction approach and
proposes the adaptive model reduction approach including the
model reduction method as a hybrid of nonlinear and linear
model reduction techniques, the algorithm enabling adaptive
switching among models of different levels of details and the
second version of the adaptive approach that can be applied to
the whole power system. In Section III the proposed
approaches are tested together with the linear model reduction
approach on the Northeast Power Coordinating Council
(NPCC) 48-machine 140-bus system. Finally, conclusions are
drawn in Section IV.
II. PROPOSED ADAPTIVE MODEL REDUCTION
A. Power System Partitioning
As it was mentioned above, in power system model
reduction the system is divided into two areas: 1) the study
area, which is the main interest of an investigator, where all

details are preserved and all disturbances are originated from;
2) the external area, which can be simplified and reduced. The
partitioned power system is shown in Fig. 1.

External
Area

V1e1e

V1s 1s

V jeej

V js sj

VNeb eNb

Study
Area

VNs b sNb

Fig. 1. Partitioned power system.

Each area of the original system is connected to other area
by several tie-lines. For every area, each tie-line is treated as a
simple fictitious generator with the internal voltage phasor
equal to the voltage phasor of the corresponding boundary bus
in the opposite area and with the armature resistance and
transient reactance equal to the resistance and reactance of the
corresponding tie-line. These fictitious generators are treated
as constant voltage sources during each iteration and represent
the electrical power injections from one area to the other area.
Therefore, voltage magnitudes and voltage angles of boundary
buses in one area are the inputs to the model of the other area.
At every iteration of the system simulation, each area is
calculated separately, then boundary bus voltages of all areas
are recalculated and their values are sent as inputs to the
corresponding area to perform the next iteration.
In this paper, power system generators are modeled based
on the detailed two-axis machine model, the non-reheat steam
turbine model, the first-order governor model and the IEEE
type 1 exciter model [18]. In each area of the partitioned
system every generator is described by the following nine
differential equations:
 i  base  i  1

Tchi Pmi   Pmi  Pgvi


 i  1

Tgvi Pgvi   Pgvi  Prefi 
Ri


K K
TA VR  VR  K A R f  Ai Fi E fd  K A Vref  Vt
i
i i
i
i
i
i
TFi
 i i

K Fi

(1)
E fdi
TFi R fi   R fi 
TFi


B E
TE E fd    K E  AE e Ei fdi  E fd  VR
i
i
i
i
i
i




 E d   Ed  X q  X q I q
Tqo
i
i
i
i
i
i

 E q   Eq  X d  X d I d  E fd
Tdo
i
i
i
i
i
 i i
2 H i 
 i  Pm  Ed I d  Eq I q  Di  i  1
i
i i
i i













3

n  9 N g and m  2 Nb respectively denote the number of

where
I di 

Ng

 Ed j Gij cos(i   j )  Bij sin(i   j ) 
j 1

Ng

  Eq j Gij sin(i   j )  Bij cos(i   j ) 
j 1

Nb

where x  δ Pm Pgv VR Rf Efd Ed Eq ω

j 1

x  R n is the state vector; u  R m is the input vector; y  R n
is the output vector.

  V j Gij sin(i   j )  Bij cos(i   j ) ,
Ng

 Ed j  Gij sin(i   j )  Bij cos(i   j ) 

I qi 

j 1

Ng

  Eq j Gij cos(i   j )  Bij sin(i   j ) 
j 1

Nb

  V j Gij cos(i   j )  Bij sin(i   j ) ,
j 1

Here, i and i are the rotor angle and speed of generator i
in rad and rad/s, respectively; base  120 rad/s is the base
speed; Pm ,

states and the number of inputs of the external area and let the
outputs of the system be the states of the system.
Then the system (1) can be described as the nonlinear
system:
x  f  x, u 
(2)

y  x

Pgvi

and Prefi

the mechanical power, the

governor output power and the reference power, respectively;
Ri is the speed regulation factor; Ed i , Eq i , X di , X qi ,
are respectively the d- and q-axis
X d i , X q i , I di and I qi
internal voltages, synchronous reactances, transient reactances
and currents all in p.u.; VRi , K Ai , R fi , K Fi , E fdi , Vrefi ,
Vti , K Ei , AEi and BEi are the voltage regulator input,
amplifier gain, rate feedback, feedback gain, field voltage,
reference voltage, terminal bus voltage, exciter gain and
exciter saturation coefficients; Time constants Hi, Tchi , Tgvi ,
 and Tdo
 are respectively the generator
TAi , TFi , TEi , Tqo
i
i
inertia, turbine charging time, governor time constant,
amplifier time constant, feedback time constant, the exciter
time constant, the q-axis open circuit time constant and the daxis open circuit time constant all in s; Di is the inertia and
damping coefficient; N g is the number of generators; Nb is

the number of boundary buses; Gij and Bij are conductance
and susceptance between generator i and generator j in p.u.;
V j is the voltage magnitude at boundary bus j in the
opposite area in p.u.;  j is the voltage angle at boundary bus

j in the opposite area in rad.
B. Model Reduction
If model reduction is applied to the external area it is
necessary to define states and inputs of the system.
Considering that every generator is described by nine
differential equations and every boundary bus has the voltage
magnitude and the voltage angle as its parameters, let

T

and u   θ V ;
T

1) Linear Model Reduction:
The system (2) can be linearized around an equilibrium
point as:
 x  Ax  Bu
(3)

 y  Cx
where x , u and y are the deviation variables of

respectively the original states, inputs and outputs; A  R nn
is the matrix of partial derivatives of the functions in (1) with
respect to each state evaluated at the equilibrium point;

B  R nm is the matrix of partial derivatives of the functions
in (1) with respect to each input evaluated at the equilibrium
point; C  R nn is the identity matrix.
The system (3) can be reduced using a linear reduction
method, for example, the balanced truncation method [5]. To
apply this method Lyapunov equations are solved to get
controllability Gramian Wc and observability Gramian Wo :
T
T
 AWc  Wc A  BB  0
 T
T
 A Wo  Wo A  C C  0

(4)

The Gramians are then used to calculate transformation
.
matrix T and its inverse T
Matrix T transforms the states from the original state
space to a new balanced state space: x  Tx.
In the resulted new balanced system, the states are arranged
in a such way that the first state is the most controllable and
the most observable and the last state is the least controllable
and the least observable. Henkel singular values show this
relationship:
1  2    i    n 1  n  0,

(5)

where
i   i  Wc Wo    ii .

Considering the above-mentioned fact only the first r states
can be kept and the rest can be truncated. H  norm of the
error of balanced truncation is bounded by the following
expression:

  2

n



i r 1

i .

(6)

4

The transformation matrix and its inverse are recalculated as
follows:
  TP
 T,
(7)
T
T  PT ,
where P   I 0  is the identity matrix, the last  n  r  rows
of which are deleted.
Thus, the balanced truncated system is represented as
follows:
 x  TBu
 x  TAT
(8)

 x
 y  CT
The system in (8) can be written in a more compact form:
 x  B u
 x  A

 x
 y  C

where

(9)

2) Proposed Hybrid Model Reduction:
In this paper, a model reduction method is proposed as a
hybrid of nonlinear and linear model reduction techniques. As
 can
shown in [19]-[20], the transformation matrices T and T
be used to reduce the nonlinear system as well. In this case,
the system can be represented as follows:





is

the

(10)

(12)

admittance

matrix

representing
N N

T
R b g
connections between fictitious generators; Y21  Y12
is the admittance matrix representing connections between the
generators of the external area and the fictitious generators.
Thus, column norms of absolute values of elements in
matrix Y21 can be used to determine which function to be
linearized as the norms describe how close electrically each
generator is to the boundary between the external area and a
study area.
Column norms are calculated by:
Nb

 Y21 ji

2

.

(13)

j 1

The nonlinear functions that correspond to the generators
with large column norms are kept nonlinear and the nonlinear
generator functions with small column norms are linearized.
Thus, the hybrid reduced system can be represented as
follows:

 ,u 
 fˆ Tx


x  T 

(14)

f
l




y  Tx



where

The system in (10) has fewer states than the original
system but it is still necessary to compute all nonlinear
functions in f. Thus, there is basically no reduction in
computation time.
To address that problem, reference [21] suggests
eliminating some of the functions. However, as it is shown in
[22] it can create large errors due to the model reduction.
In the proposed hybrid model reduction approach, the
functions that have the least contributions to the dynamics
between the external area and a study area are not eliminated
but linearized. To evaluate contributions of the functions, let
us consider the expressions for the d-axis current and the qaxis current of generators in the external area as these
expressions have most nonlinearities and are used in 33% of
all differential equations in (1).
Nonlinearities in the expressions are cosine and sine
functions and coefficients of these functions are conductances
and susceptances between generators including fictitious
generators representing the boundary between the external
area and the study area. These values are real and imaginary
parts of elements of the admittance matrix:
(11)
Yij  Gij  jBij .
The matrix can be divided into four submatrices:
Y12 
Y
Y   11
,
 Y21 Y22 

Y22  R Nb  Nb

i 

  R r r ,
  TAT ,
A
A
B  R rm ,
B  TB ,
  CT ,
  R nr .
C
C

 ,u
x  Tf Tx


y  Tx

N N

where Y11  R g g is the admittance matrix representing
connections between generators inside the external area;



ˆ x  Bˆ u  xˆ 0 , xˆ 0  Px
ˆ 0,
fl  A
ˆ  PAT
ˆ ,
A
ˆ ,
Bˆ  PB

ˆ  R (nq)r ,
A
Bˆ  R ( nq)m ,

Vector f̂ comprises the functions that are kept nonlinear; fl
has the linearized functions; x0 is the initial state vector; P̂ is
the identity matrix with deleted rows that correspond to the
functions in fˆ ; q is the number of nonlinear functions in fˆ.
The system in (14) can be rewritten as:

 ,u


fˆ Tx
x  T 


A
ˆ x  Bˆ u  xˆ 0 





y  Tx





(15)

C. Adaptive Switching Algorithm
Considering that the linearly reduced system gives
satisfactory performance during small disturbances, the
duration of a large disturbance is short, and majority of the
time a system is under small or no disturbance it is reasonable
to change the type of the model reduction of the external
system to increase the accuracy and speed of the system
simulation. Proposed adaptive algorithm changes the
complexity of model reduction of the external area based on
the current condition. The adaptive algorithm is shown in
Fig. 2.
During the fault-on period, the original, fully detailed
system model is used as the maximum accuracy of the system
model is required and the duration of a fault is limited to tens

5

of milliseconds, which does not increase much the simulation
time. In the post-fault period, when the angle deviation i of
any generator in a study area exceeds a preset threshold
max , the external area is reduced using the hybrid model
reduction method, which keeps the balance between accuracy
and speed of simulation when the disturbance is large. In prefault and post-fault periods when all angle deviations are
within the threshold, the external area is reduced using a linear
model reduction method. This guarantees that most of the time
when there is no disturbance or variation of the state is very
small, the fastest model reduction method is applied.
To calculate the rotor angle deviation the generator with
the smallest column norm is used as a reference; i.e., the
reference generator is the electrically furthest generator from a
study area and reacts to disturbances in study areas the least. If
the time for a rotor angle deviation i exceeding threshold
max , denoted by tth , is longer than a preset limit tthmax a

iteration. To eliminate the partitioning error the second
version of the adaptive approach is proposed. This version is
applied to the whole system that is treated as just one area.
Without the partitioning, there is no need in the concept of
fictitious generators representing boundary buses of the study
area and the external area and the expressions in (1) for the daxis current and the q-axis current are simplified as:

 and B of
new operating condition is obtained and matrices A
the linearly reduced system are recalculated. This action
corrects the adaptive algorithm after a large change of the
system state.

As the only area of the system contains all generators
including the generators from the study area whose dynamics
are of the main interests, the transformation and truncation of
the states are not performed and the performance
improvement comes only from the linearization of nonlinear
functions. In the absence of inputs from the boundary between
the study area and the external area, the control matrix B is
eliminated from (15) and the system used in the second
version of the adaptive approach is simplified as:

 fˆ  x, u  
x   

 Ax  xˆ 0 
(16)




y  x

Start
Calculate
reduced system
matrices, select
nonlinear
functions

Yes

Fault‐on
condition?

Get the system
operating
condition

No

i  max

Simulate the
system reduced
by the linear
model reduction

 Ed j Gij cos(i   j )  Bij sin(i   j ) 
j 1

Ng

  E q j  Gij sin( i   j )  Bij cos( i   j )  ,
j 1

I qi 

Ng



j 1

E d j   G ij sin(  i   j )  B ij co s(  i   j ) 

Ng

  E q j  G ij co s(  i   j )  B ij sin (  i   j )  .
j 1



ˆ ,
A  PA
A  R ( n  q )n .
All nonlinear functions representing generators of the study

No

Simulate the
system reduced
by the hybrid
model reduction

Yes

No

Ng

where
Yes

Simulate the
original system

I di 

tth  tthmax

No

Is simulation
over?
Yes
End

Fig. 2. Adaptive switching algorithm.

D. Adaptive Model Reduction of the System without
Partitioning
Partitioning the system into the study area and the external
area creates a specific error. The error is caused by the fact
that the inputs (boundary bus voltage magnitudes and
boundary bus voltage angles) are calculated at the previous
iteration of the simulation, i.e. the inputs are lagging by one

area are contained in f̂ . The list of linearized functions
corresponding to generators of the external area is the same as
in the adaptive approach applied to the partitioned system
described above. The adaptive switching is performed
between the system in (16) and the simplified version of the
linearized system in (3):
x  Ax
(17)

y  Cx
As all generators of the whole system including the
generators of the study area are linearized in (17) the angle
deviation threshold max of the adaptive switching
algorithm is set to small value to enable switching after the
oscillations are damped enough.
III. CASE STUDIES
Comprehensive case studies are conducted to compare the
proposed adaptive model reduction approaches with the
traditional linear model reduction approach. A realistic power
system model is tested. For the system, the study area is
defined and retained with original, detailed models and the
rest of the system is defined as the external area to be reduced
respectively by different approaches. Then, time-domain

6

contingency simulation using each reduced system model is
conducted and compared with the simulation using the
original system model. In addition, the approach is validated
during different post-fault operating conditions.
A. Temporary bus fault tests
The linear model reduction approach and the adaptive
approaches described above are applied to the NPCC 140-bus
48-generator power system. The study area is set to be the
New England region of the NPCC system, which has 9
generators. The external area is set to be the rest of the system,
which has 39 generators. The external area has 39×9=351
state variables and 39×4=156 nonlinear functions as the first 5
states of each generator is described by linear differential
equations. The external area is connected to the study area by
two tie-lines. The partitioned NPCC system is shown in Fig.
3.
92

External Area
(39 generators)

97

Hydro Quebec

96
95

IESO
99

100

98

48
47
49

93
46

ISO‐NE

94

100 MVA. Setting the base voltage to the common generator
terminal voltage of 20 kV the new thresholds is calculated:
S
100
Yth  1 base  1
 0.25S
(19)
2
Vbase
202
When the proposed approach is applied to a different
system, the threshold is converted to a per unit value using the
new system base power value and the base voltage of 20 kV.
Another case study is performed to select the angle
deviation threshold max for the adaptive switching
algorithm applied to the partitioned system. The threshold is
decreased with 1-degree increments from 180 degrees until all
rotor angle errors for all generators following any of the
contingencies are below 6 degrees. After this case study,
max is set to 67 degrees. The threshold for the adaptive
approach applied to the unpartitioned system is set to 6
degrees to ensure rotor angle errors of generators representing
study area are also below 6 degrees.
The simulations are performed in MATLAB R2015a on the
computer with 4-GHz AMD FX-8350 processor and 8 GB of
memory. The duration of simulation is set to 16 seconds,
integration time step is set to 0.01 seconds.

26
20

TABLE I
COMPARISON OF ROOT MEAN SQUARE ERROR OF ROTOR ANGLE

NYISO

91
108

109

111
12 0

1 04

119

1 02

54

11 0
106

10 3

56
55

1 14

1 07

89
90

10 1

57

45

50

19
44

52
51
43

53

27

62

88

42
65

40

72

29

85

87

121

60

28

37

38
63

86
1 22

16

59

1 15
1 17

84

83

18

64

61

116
1 13

1 18

17

41

58
10 5

8

30
1 12

71

68

137

7
15

39
9

13 6

70

67

69

6

66

12

1 23

13
138

31
5
35

73

32

34

3

33

3

2
4

77

36

74

1

14
10

24

11

75

76

80
126

13 5

82
78

79

81
134

125

1 32

12 4
1 27

1 33

12 8

131

130
129

Study Area
(9 generators)

1 39

Fig. 3. Partitioned NPCC system.

The H  norm of the balanced truncation error is set to be
equal to 105 :

105  2

n



i r 1

i .

Gen.

CCT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

26
26
23
26
23
23
26
23
26
26
23
23
23
23
23
23
26
23
26
26
23
23
23
23
23
26
23
26
23
23
23
26
26
23
26
23

0.13
0.13
0.39
0.12
0.09
0.08
0.05
0.05
0.06
0.13
0.19
0.1
0.12
0.13
0.07
0.05
0.03
0.04
0.03
0.03
0.22
0.19
0.20
0.18
0.22
0.03
0.10
0.12
0.06
0.06
0.08
0.11
0.11
0.15
0.14
0.21

22

23

21

PJM

MISO

Bus

25

23

(18)

This corresponds to truncation of 200 out of 351 state
variables. To find the threshold value of the column norm of
the admittance matrix that decides if generator is electrically
close to the boundary a case study is performed. The system in
(15) is used in this case study. The threshold is decreased with
0.1-p.u. increments from 10 until all rotor angle errors for all
generators following any of the contingencies are below 6
degrees. After the case study, the threshold is set to 1 p.u. The
nonlinear functions corresponding to generators with column
norms of admittance matrix less than 1.0 p.u. are linearized.
This corresponds to linearization of 136 out of 156 nonlinear
functions of the external area in the case of partitioned system
and 136 out of 192 nonlinear functions of the unpartitioned
system.
The column norm threshold per unit value in not useful if
the proposed approach is applied to a different system with a
different base power. To make it more practical it is converted
to value in siemens. The NPCC system has the base power of

Partitioned
Linear
Adapt.
8.85
1.33
3.10
0.86
25.94
5.51
1.50
0.75
0.59
0.48
0.72
0.60
0.28
0.28
4.87
4.41
3.26
2.83
10.49
2.90
8.81
1.40
1.89
0.98
2.57
0.97
2.83
1.00
0.38
0.38
0.87
0.78
0.30
0.30
0.33
0.33
0.11
0.11
0.27
0.27
12.97
1.65
6.55
1.32
3.83
1.07
5.66
1.46
20.45
2.22
0.35
0.34
0.29
0.27
0.28
0.30
3.15
3.48
1.21
1.35
0.40
0.43
1.25
0.79
1.05
0.65
3.79
0.89
3.14
0.82
4.03
0.92

Unpartitioned
Adapt.
1.06
0.57
5.12
0.41
0.29
0.30
0.47
0.29
0.29
1.26
1.46
0.44
0.49
0.51
0.33
0.31
0.59
0.66
0.95
1.05
1.89
1.06
0.71
0.91
3.87
0.86
0.46
0.42
0.43
0.38
0.25
0.42
0.39
0.63
0.52
0.61

7

To compare the approach performance during a large
disturbance, a three-phase short circuit fault lasting for the
critical clearing time (CCT) is created separately at every bus
of the study area. The errors of all outputs (state variables) of
all generators in the study area are analyzed and the rotor
angle state variable has the largest error for each of the
generators. For every fault, the generator with the largest error
of the rotor angle is found and used to compare the
approaches. The results of comparison of rotor angle root
mean square errors of the linear model reduction approach and
the adaptive approach are shown in Table I. As it can be seen
the linear model reduction approach cannot guarantee
satisfactory performance during large disturbances generating
the errors of more than 25 degrees while the proposed
adaptive approach keeps the error for all disturbances within 6
degrees both for the case with the system partitioned into the
study and the external area and the case with only one area
representing the whole system.
170

Rotor angle (degrees)

140

l2
l4

110
80
50
20

Original system

 10

Linearly reduced
Partitioned
Adaptively reduced

 40
 70
 70

0

2

4

6

8

10

12

14

16

System
States

, degrees

Unpartitioned
adaptively reduced
5.1  100

Pm , p.u.

1.8  103

4.1  104

3.4  104

Pgv , p.u.

2.3  102

5.1  103

4.5  103

140

VR , p.u.

1.7  101

3.7  102

4.1  102

110

R f , p.u.

1.3  102

3.4  103

3.2  103

E fd , p.u

9.8  102

2.2  102

2.3  102

Ed , p.u.

6.9  102

1.4  102

1.5  102

Eq , p.u.

1.1  102

2.8  103

2.6  103

, p.u.

4.7  103

1.0  103

9.7  104

16

170

Rotor angle (degrees)

(19)
,
N
where N is the number of simulation steps; xij and xˆij are

Partitioned
adaptively reduced
5.5  100

170

l4

j 1

TABLE II
COMPARISON OF ROOT MEAN SQUARE ERROR OF STATES OF GENERATOR 23

Fig. 4. Rotor angle of generator 23 following the fault at bus 3.

l3

i 

2

Linearly
reduced
2.24  101

Timet (s)

0

l1

N

  xij  xˆij 

respectively the values of i-th state variable of generator 23 of
the original system and the reduced system at time step j.
The results of calculation are shown in Table II. The
proposed adaptive approaches reduce the error by 74% to 81%
for every state variable compared to the linear model
reduction approach.
A fault with the duration equal to CCT is the worst-case
scenario that can cause the largest error. The worst-case errors
of the proposed adaptive model reduction approach are small
enough to justify its applicability to power system stability
studies.

170

l1

following this fault, and thus its rotor angle is used for the
comparison. The results of simulation are shown in Fig. 4
when the adaptive approach is applied to partitioned system
and in Fig. 5 when the adaptive approach is applied to
unpartitioned system. As it can be seen from Fig. 4 and Fig. 5,
if the external system is reduced using linear model reduction,
the rotor angle trajectory differs significantly from the rotor
angle trajectory of the original system whereas the rotor angle
trajectory of the system reduced by the adaptive approaches
follows accurately the original trajectory.
To present in more detail the difference between
approaches, the root mean square errors of all states are
calculated using the following expression:

80
50
20

Original system

 10

Linearly reduced
Unpartitioned
Adaptively reduced

 40
 70
 70

0
0

2

4

6

8

10

12

Timet (s)

14

16
16

Fig. 5. Rotor angle of generator 23 following the fault at bus 3.

The fault at bus 3 is the largest disturbance of the NPCC
system and generator 23 has the largest rotor angle error

In addition to the accuracy the approaches are compared in
terms of simulation time as shown in Table III. In Table III the
large disturbance corresponds to the fault at bus 3 and the
small disturbance corresponds to the fault at bus 17. The
results show that the proposed adaptive approach applied to
the partitioned system reduces the simulation time by 57%
during the large disturbance and by 59% during the small
disturbance compared to the original system. If the adaptive
approach is applied to the unpartitioned system the simulation
time is reduced respectively by 73% and 84%. Better speed

8

performance of the second version of the adaptive approach is
caused by the fact that the whole system is switched to the
linearized model which is the fastest model. Especially it is
clear in the case of the small disturbance as the switching to
the linear model happens earlier.
Thus, the proposed adaptive approach provides both high
accuracy and high simulation speed.
TABLE III
COMPARISON OF SIMULATION TIME

System
Original unpartitioned
Partitioned and linearly reduced
Partitioned and adaptively reduced
Unpartitioned and adaptively reduced

Simulation time, seconds
Large
Small
disturbance
disturbance
3.7
0.9
1.6
1.5
1.0
0.6

[4]

[5]

[6]

[7]

[8]

[9]
[10]

B. Test of operating condition change
To test how robustly the adaptive approach performs
against a change of the operating condition, the temporary
fault at bus 3 representing the largest disturbance is changed
to a permanent fault cleared by tripping one of the lines
connected to the bus or by complete isolation of the bus from
the system by tripping all lines connected to the bus.
Following the bus fault, the operating condition of the system
is changed. The results of simulations are shown and
compared in Table IV, from which the proposed adaptive
approach can reduce the system after the change of the
operating condition maintaining the accuracy of simulation.

[13]

TABLE IV
COMPARISON OF FAULTS LEADING TO A NEW OPERATING CONDITION

[16]

Fault

Temporary fault at bus 3
Fault at line 3-2 followed by the line trip
Fault at line 3-4 followed by the line trip
Fault at bus 3 followed by the bus trip

Rotor angle RMS error,
degrees
Partitioned Unpartitioned
adaptively
adaptively
reduced
reduced
5.51
5.12
5.48
5.09
5.28
4.80
2.74
2.04

IV. CONCLUSIONS
This paper has proposed two versions of a new adaptive
model reduction approach that can be applied to traditional
partitioned system or to unpartitioned system for fast power
system simulation. The approach is capable of accurate
representation of the original power system model with
significant reduction in computational time.
V. REFERENCES
[1]
[2]

[3]

J. H. Chow, Power system coherency and model reduciton. Springer,
2013, pp. 1-3.
J. H. Chow, R. J. Galarza, P. Accari, W. W. Price, “Inertial and slow
coherency aggregation algorithms for power system dynamic model
reduction,” IEEE Trans. Power Systems, vol. 10, no. 2, pp. 680–685,
May 1995.
R. J. Galarza, J. H. Chow, W. W. Price, A. W. Hargrave, P. M. Hirsch,
“Aggregation of exciter models for constructing power system dynamic

[11]

[12]

[14]

[15]

[17]

[18]
[19]

[20]

[21]

[22]

equivalents,” IEEE Trans. Power Systems, vol. 13, no. 3, pp. 782–788,
Aug. 1998.
F. Ma, V. Vittal, “Right-sized power system dynamic equivalents for
power system operation,” IEEE Trans. Power Systems, vol. 26, no. 4, pp.
1998–2005, Nov. 2011.
B. C. Moore, “Principal component analysis in linear systems:
Controllability, observability and model reduction,” IEEE Trans.
Automatic Control, vol. 26, no. 1, pp. 17–32, Feb. 1981.
S. Liu, “Dynamic-data driven real-time identification for electric power
systems,” Ph.D. dissertation, Univ. of Illinois at Urbana-Champaign,
Urbana, IL, 2009.
C. Sturk, L. Vanfretti, Y. Chompoobutrgool, H. Sandberg, “Coherencyindependent structured model reduction of power systems,” IEEE Trans.
on Power Systems, vol. 29, no. 5, pp. 2418–2426, Sep. 2014.
L. T. Pillage, R. A. Rohrer, “Asymptotic waveform evaluation for timing
analysis,” IEEE Trans. on Computer-Aided Design of Integrated
Circuits and Systems, vol. 9, no. 4, pp. 352–366, Apr. 1990.
E. J. Grimme, “Krylov projection methods for model reduction,” Ph.D.
dissertation, Univ. of Illinois at Urbana-Champaign, Urbana, IL, 1997.
D. Chaniotis, “Krylov subspace methods in power system studies,”
Ph.D. dissertation, Univ. of Illinois at Urbana-Champaign, Urbana, IL,
2001.
J. Hahn, T. F. Edgar, “An improved method for nonlinear model
reduction using balancing of empirical Gramians,” Computers and
Chemical Engineering, vol. 26, no. 10, pp. 1379–1397, Oct. 2002.
M. Rathinam, L. R. Petzold, “A new look at proper orthogonal
decomposition,” SIAM Journal on Numerical Analysis, vol. 41, no. 5,
pp. 1893–1925, 2003.
K. Fujimoto, J. M. A. Scherpen,”Singular value analysis and balanced
realizations for nonlinear systems,” in Model Order Reduction: Theory,
Research Aspects and Applications, Springer, 2008, pp. 251–272.
E. I. Verriest, “Time variant balancing and nonlinear balanced
realizations”, in Model Order Reduction: Theory, Research Aspects and
Applications, Springer, 2008, pp. 213–250.
J. Qi, J. Wang, H. Liu, A. D. Dimitrovsky, “Nonlinear model reduction
in power systems by balancing of empirical controllability and
observability covariances,” IEEE Trans. on Power Systems, in press.
M. Striebel, J. Rommes, “Model order reduction of nonlinear systems in
circuit simulation: status and applications,” in Model Reduction for
Circuit Simulation, Springer, 2011, pp. 289–301.
S. Wang, S. Lu, N. Zhou, G. Lin, M. Elizondo, M. A. Pai, “Dynamicfeature extraction, attribution, and reconstruction (DEAR) method for
power system model reduction,” IEEE Trans. on Power Systems, vol. 29,
no. 5, pp. 2049–2059, Sep. 2014.
P. W. Sauer, M. A. Pai, Power System Dynamics and Stability.
Englewood Cliffs, NJ: Prentice Hall, 1998, pp. 183-188.
A. A. Desrochers, R. Y. Al-Jaar, “A method for high order linear system
reduction and nonlinear system simplification,” Automatica, vol. 21, no.
1, pp. 93–100, Jan. 1985.
X. Ma, J. A. De Abreu-Garcia, “On the computation of reduced order
models of nonlinear systems using balancing technique,” in Proc.
Decision and Control, 1988, pp. 1165–1166.
A. Verhoeven, J. ter Maten, M. Striebel, R. Mattheij, “Model order
reduction for nonlinear IC models,” in System Modeling and
Optimization, Springer, 2007, pp. 476–491.
M. Striebel, J. Rommes, “Model order reduction of nonlinear systems:
status, open issues, and applications,” Chemnitz University of
Technology, Chemnitz, Germany, CSC/08-07, Nov. 2008.

