ECEN 667 Power System Stability: Lecture 13: Governors, PID Controllers
ECEN 667 Power System Stability: Lecture 13: Governors, PID Controllers
ECEN 667 Power System Stability: Lecture 13: Governors, PID Controllers
2
PID Controllers
• Governors and exciters often use proportional-integral-
derivative (PID) controllers
– Developed in 1890’s for automatic ship steering by observing
the behavior of experienced helmsman
• PIDs combine
– Proportional gain, which produces an output value that is
proportional to the current error
– Integral gain, which produces an output value that varies with
the integral of the error, eventually driving the error to zero
– Derivative gain, which acts to predict the system behavior.
This can enhance system stability, but it can be quite
susceptible to noise
3
PID Controller Characteristics
• Four key characteristics
of control response are
1) rise time, 2) overshoot,
3) settling time and
4) steady-state errors
Bus 4
90 MW
slack
Bus 3
11.59 Deg 10 MW
1.0971 pu
Further details on
tuning are covered in
IEEE Std. 1207-2011
9
Tuning PID Controller Example
• Figure shows the Ziegler-Nichols for a P, PI and PID
controls. Note, this is for stand-alone, not
interconnected operation
10
Example KI and KP Values
• Figure shows example KI and KP values from an
actual system case
About 60%
of the models
also had a
derivative term
with an average
value of 2.8,
and an average
TD of 0.04 sec
11
Non-windup Limits
• An important open question is whether the governor PI
controllers should be modeled with non-windup limits
– Currently models show no limit, but transient stability
verification seems to indicate limits are being enforced
• This could be an issue if frequency goes low, causing
governor PI to "windup" and then goes high (such as in
an islanding situation)
– How fast governor backs down depends on whether the limit
winds up
12
PI Non-windup Limits
• There is not a unique way to handle PI non-windup
limits; the below shows two approaches from IEEE Std
421.5
Another
common
approach
is to cap the
output and
put a non-
windup limit
on the
integrator
13
PI Limit Problems with Actual
Hydro Models
• A recent research project comparing transient
stability packages found there were significant
differences between hydro model implementations
with respect to how PI limits were modeled
– One package modeled limits but did not document them,
another did not model them; limits were recommended
at WECC MVWG in May 2014
14
PIDGOV Model Results
• Below graph compares the Pmech response for a two
bus system for a frequency change, between three
transient stability packages
Packages
A and B
both say they
have no
governor
limits, but
B seems to;
PowerWorld
can do
either
approach
15
GGOV1
• GGOV1 is a relatively newer governor model
introduced in early 2000's by WECC for modeling
thermal plants
– Existing models greatly under-estimated the frequency drop
– GGOV1 is now the most common WECC governor, used with
about 40% of the units
• A useful reference is L. Pereira, J. Undrill, D. Kosterev,
D. Davies, and S. Patterson, "A New Thermal Governor
Modeling Approach in the WECC," IEEE Transactions
on Power Systems, May 2003, pp. 819-829
16
GGOV1: Selected Figures from
2003 Paper
Diablo Canyon is California’s last nuclear plant, with Unit 1 now scheduled to
shutdown in 2024 and Unit 2 in 2025. 17
GGOV1 Block Diagram
GGOV1 and
the related
GGOV3 are
the most
common
governors in
WECC, with
more than
40% in 2019
18
Transient Stability
Multimachine Simulations
• Next, we'll be putting the models we've covered so far
together
• Later we'll add in new model types such as stabilizers,
loads and wind turbines
• By way of history, prior to digital computers, network
analyzers were used for system stability studies as far
back as the 1930's (perhaps earlier)
– For example see, J.D. Holm, "Stability Study of A-C Power
Transmission Systems," AIEE Transactions, vol. 61, 1942, pp.
893-905
• Digital approaches started appearing in the 1960's
19
Transient Stability
Multimachine Simulations
• The general structure is as a set of differential-
algebraic equations
– Differential equations describe the behavior of the machines
(and the loads and other dynamic devices)
– Algebraic equations representing the network constraints
In EMTP applications the
transmission line delays
decouple the machines; in
transient stability they are
assumed to be coupled
together by the algebraic
network equations
20
General Form
• The general form of the problem is solving
x f (x, y , u)
0 g ( x, y )
where x is the vector of the state variables (such
as the generator 's), y is the vector of the algebraic
variables (primarily the bus complex voltages), and
u is the vector of controls (such as the exciter voltage
setpoints)
21
Transient Stability
General Solution
• General solution approach is
– Solve power flow to determine initial conditions
– Back solve to get initial states, starting with machine
models, then exciters, governors, stabilizers, loads, etc
– Set t = tstart, time step = Dt, abort = false
– While (t <= tend) and (not abort) Do Begin
• Apply any contingency event
• Solve differential and algebraic equations
• If desired store time step results and check other conditions
(that might cause the simulation to abort)
• t = t + Dt
– End while
22
Algebraic Constraints
• The g vector of algebraic constraints is similar to the
power flow equations, but usually rather than
formulating the problem like in the power flow as real
and reactive power balance equations, it is formulated
in the current balance form Simplest
cases
Ι (x, V ) Y V or Y V Ι (x, V ) 0 can have
where Y is the n n bus admittance matrix I independent
(Y G jB), V is the complex vector of of x and V,
allowing a
the bus voltages, and I is the complex direct solution;
vector of the bus current injections otherwise we
need to iterate
23
Why Not Use the
Power Flow Equations?
• The
power flow equations were ultimately derived
from
I=YV
• However, the power form was used in the power flow
primarily because
– For the generators the real power output is known and either
the voltage setpoint (i.e., if a PV bus) or the reactive power
output
– In the quasi-steady state power flow time frame the loads can
often be well approximated as constant power
– The constant frequency assumption requires a slack bus
• These assumptions do not hold for transient stability
24
Algebraic Equations for
Classical Model
• To introduce the coupling between the machine models
and the network constraints, consider a system
modeled with just classical generators and impedance
loads
In this example
because we are
using the classical
model all values
are on the network
reference frame
Image Source: Fig 11.15, Glover, Sarma, Overbye, Power System Analysis and Design, 5th Edition, Cengage Learning, 2011
25
Algebraic Equations for
Classical Model
• Replace the internal voltages and their impedances by
their Norton Equivalent
Ei i 1
Ii , Yi
Rs ,i jX d ,i Rs ,i jX d ,i
26
Swing Equation
• The first two differential equations for any
synchronous machine correspond to the swing
equation
d i
i s i
dt
2H i di 2H i d i
TMi TEi Di i
s dt s dt
with TEi de,i i qi qe,i i di
27
Swing Equation Speed Effects
• There is often confusion about these equations because
of whether speed effects are included
– Recognizing that often w ws (which is one per unit), some
transient stability books have neglected speed effects
• For a rotating machine with a radial torque,
power = torque times speed
• For a subtransient model
E V ( Rs jX ) I , Ed jEq q j d
TE dI q qI and
PE TE Ed jEq I d jI q EdI d EqI q
28
Classical Swing Equation
• Often in an introductory coverage of transient stability
with the classical model the assumption is w w ws
so the swing equation for the classical model is given
as d i
i s i
dt
2H i d i
PMi PEi Di i
s dt
with PEi Ei i Ei i Vi Yi
30
Outline for Next Several Slides
• The next several slides will provide basic coverage
of the solution process, partitioned explicit, then
the simultaneous-implicit approach
• We'll start out with a classical model supplying an
infinite bus, which can be solved by embedded the
algebraic constraint into the differential equations
We'll start out just solving x f ( x)
and then will extend to solving the full problem of
x f (x, y , u)
0 g(x, y )
31
Classical Swing Equation with
Embedded Power Balance
• With a classical generator at bus i supplying an infinite
bus with voltage magnitude Vinf, we can write the
problem without algebraic constraints as
d i
i s i i. pus
dt
d i , pu 1 Ei Vinf
PMi sin i Di i , pu
dt 2H i X th
Ei Vinf
with PEi sin i Note we are using
X th the per unit speed
approach
32
Explicit Integration Methods
• As covered on the first day of class, there are a
wide variety of explicit integration methods
– We considered Forward Euler, Runge-Kutta, Adams-
Bashforth
• Here we will just consider Euler's, which is easy to
explain but not too useful, and a second order
Runge-Kutta, which is commonly used
33
Forward Euler
• Recall the Forward Euler approach is approximate
dx x
x f (x(t )) as
dt t
Then
x(t t ) x(t ) t f ( x(t ))
34
Infinite Bus GENCLS Example
using the Forward Euler's Method
• Use the four bus system from before, except now gen 4
is modeled with a classical model with Xd'=0.3, H=3
and D=0; also we'll reduce to two buses with
equivalent
Bus 2
line reactance, moving the gen from bus 4 to 1
GENCLS Bus 1
Infinite Bus
X=0.22
slack
In this example Xth = (0.22 + 0.3), with the internal voltage
giving E'1=1.281 and d1=
35
Infinite Bus GENCLS Example
• The
associated differential equations for the bus 1
generator are
d 1
1, pus
dt
d 1, pu 1 1.281
1 sin 1
dt 23 0.52
• The value of PM1 = 1 is determined from the initial
conditions, and would stay constant in this simple
example without a governor
• The value d 1= is readily verified as an equilibrium
point (which is 0.418 radians)
36
Infinite Bus GENCLS Example
• Assume a solid three phase fault is applied at the
generator terminal, reducing PE1 to zero during the
fault, and then the fault is self-cleared at time Tclear,
resulting in the post-fault system being identical to
the pre-fault system
– During the fault-on time the equations reduce to
d 1
1, pus That is, with a solid fault
dt on the terminal of the
d 1, pu 1 generator, during
1 0 the fault PE1 = 0
dt 23
37
Euler's Solution
• The initial value of x is
(0) 0.418
x(0)
(0
pu ) 0
38
Euler's Solution
• At t = Tclear the fault is self-cleared, with the
equations changing to
d
pus
dt
d pu 1 1.281
1 sin
dt 6 0.52
• The integration continues using the new equations
39
Euler's Solution Results (Dt=0.02)
• The below table gives the results using Dt = 0.02 for
the beginning time steps
Time
0
Gen 1 Rotor Angle, Degrees Gen 1 Speed (Hz)
23.9462 60
This is saved as
0.02 23.9462 60.2 PowerWorld case
0.04 25.3862 60.4
0.06 28.2662 60.6
B2_CLS_Infinite.
0.08 32.5862 60.8 The integration
0.1 38.3462 61
0.1 38.3462 61 method is set to
0.12 45.5462 60.8943 Euler's on the
0.14 51.9851 60.7425
0.16 57.3314 60.5543 Transient Stability,
0.18 61.3226 60.3395
0.2 63.7672 60.1072
Options, Power
0.22 64.5391 59.8652 System Model page
0.24 63.5686 59.6203
0.26 60.8348 59.3791
0.28 56.3641 59.1488
40
Generator 1 Delta: Euler's
• The below graph shows the generator angle for varying
values of Dt; numerical instability is clearly seen
41
Second Order Runge-Kutta
• Runge-Kutta methods improve on Euler's method by
evaluating f(x) at selected points over the time step
• One approach is a second order method (RK2) in which
1
x t t x t k 1 k 2 This is also known
2
where as Heun's method
or as the Improved
k1 t f x t Euler's or Modified
k 2 t f x t k 1 Euler's Method
• That is, k1 is what we get from Euler's; k2 improves on
this by reevaluating at the estimated end of the time step
• Error varies with the cubic of the time ste
42
Second Order Runge-Kutta (RK2)
• Again assuming a time step Dt = 0.02 seconds, and a
Tclear of 0.1 seconds, then using Heun's approach
(0 ) 0.418
x(0)
pu (0 ) 0
0 0 0.418
k 1 0.02 , x(0 ) k 1
0 .1667 0 . 00333 0 .00333
1.257 0.0251
k 2 0.02
0 .1667 0 . 00333
0.418 1 0.431
x(0.020 ) k1 k 2
0 2 0 .00333
43
RK2 Solution Results (Dt=0.02)
• The below table gives the results using Dt = 0.02 for
the beginning time steps
This is saved as
Time Gen 1 Rotor Angle, Degrees Gen 1 Speed (Hz)
0 23.9462 60 PowerWorld case
0.02 24.6662 60.2
0.04 26.8262 60.4 B2_CLS_Infinite.
0.06 30.4262 60.6 The integration
0.08 35.4662 60.8
0.1 41.9462 61 method should be
0.1 41.9462 61 changed to Second
0.12 48.6805 60.849
0.14 54.1807 60.6626 Order Runge-Kutta on
0.16 58.233 60.4517 the Transient Stability,
0.18 60.6974 60.2258
0.2 61.4961 59.9927 Options, Power System
0.22 60.605 59.7598 Model page
0.24 58.0502 59.5343
0.26 53.9116 59.3241
0.28 48.3318 59.139
44
Generator 1 Delta: RK2
• The below graph shows the generator angle for varying
values of Dt; much better than Euler's but still the
beginning of numerical instability with larger values of
Dt
45