IEEE PES GM, July 26 ‐ 30, 2015
2015 PES Prabha Kundur Award Winner Talk
PSDP Committee Meeting
PSS Phase Shaping Design & Equivalents for Transmission Networks Containing Distributed Parameter Lines
Nelson Martins
2
Agenda
• Part I: PSS Phase Shaping Design
• Part II: Equivalents for Transmission Networks
Containing Distributed Parameter Lines
3
Part I:
PSS Phase Shaping Design 4
Outline of Part I
• Synthetic system method:
– Basic assumptions
“A robust application of the residues method"
– PSS phase shaping for generator & exciter set
– PSS tuning verification
– PSS tuning validation in multimachine test system
• Application to a large practical multimachine system
• Conclusions for Part I
5
Synthetic System Method (1/5)
The synthetic system method for PSS phase shaping:
• An extension of the classical GEP method.
• Uses a set of single‐machine infinite‐bus (SMIB)
systems (low‐order models).
• Produces a set of electromechanical poles whose
frequencies range from intra‐plant to inter‐area, as
well as their sensitivities to the parameters of the PSS
under design.
• Designs the PSS phase compensation to
simultaneously contribute to the damping of
electromechanical oscillations.
6
Synthetic System Method (2/5)
• Generator dispatch and line terminal voltages are kept constant
in family of 2‐bus power flows that abide to the constraint:
Pg = P12 + PL = Constant. Flows are solved analytically.
• The shunt load (PL) size increases with the line reactance X12
while the line loading (P12) and frequency of the
electromechanical mode are reduced.
7
Synthetic System Method (3/5)
The voltage phase angle across the line is varied such that for small
(large) values of X12 the shunt load is near to its minimum
(maximum) value: θ 12 = A B X 12
8
Synthetic System Method (4/5)
• Ideal phase compensation (IPC):
– Calculated from the pole sensitivity to the stabilizer's gain.
– Ideally, the compensated sensitivities should have an angle
of 180° at all frequencies.
• The PSS phase characteristic is required to lie within the PCB.
The preferred PCB in US
is +/‐ 30º with respect
to the IPC curve
9
W(s)/Vref(s) Residues for electromechanical poles
10
Variation in the ZIP load characteristic only impacts the residue magnitude
11
Synthetic System Method (5/5)
The value for the PSS gain is chosen to push the poles with
frequency above 1 Hz (in the electromechanical locus line,
below) to a minimum damping target of 10%.
12
Basic PSS Tuning for Multimachine Systems
• Each generator & exciter set is connected to the Synthetic
System test bench to conduct the electromechanical frequency
scan test to determine its PSS phase shaping to best damp
oscillations over the inter‐area to intra‐plant frequency range;
• The PSS gain for each generator is initially set to adequately
damp the modes in the upper half of the electromechanical
locus line.
• The same procedure is repeated to the generators that are best
candidates for oscillation damping control in power systems of
any topology.
• The gains for all PSSs in multimachine systems need be
coordinated, considering multiple generation/loading scenarios
and relying on nonlinear simulation studies.
13
Method validation on 3‐generator test system
Reactances given on 100 MVA base
14
Multi‐machine system description
ƒ Generators #1 and #2:
ƒ Typical hydro‐generator parameters.
ƒ Each with half the capacity of the generator
utilized in the 2‐bus synthetic system studies .
ƒ Generator #3: typical turbo‐generator parameters.
ƒ Load types:
Bus
#4
#5
Active power
80% constant‐I
20% constant‐Z
100% constant‐P
Reactive Power
100% constant‐Z
100% constant‐Z
15
Rotor speed mode‐shapes for
3‐generator test system
Intra-plant
Intra-area
Inter-area
16
Comparing the root‐loci results
Blue x: results for 3‐generator system;
Red + and continuous lines: synthetic system results.
17
3‐generator system with no PSSs
Residues for (Δω/ΔVref) transfer function of generator #1
calculated with no PSS (blue x: 3‐generator system results;
red +: synthetic system results).
18
3‐generator system with one PSS
Residues of the (ΔPSSoutput/ΔVref) transfer function of
generator #1 calculated with incremental gain PSS (blue x:
3‐generator syst. results; red +: synthetic system results).
19
System with two PSS (gen. #1 and #2)
Blue x: results for 3‐generator system;
Continuous line: synthetic system results.
20
Basic PSS Tuning for Multimachine Systems
• Step 1: Modal analysis of chosen critical
scenarios
• Step 2: Characterization of critical modes
• Step 3: Selection of candidate generators
• Step 4: PSS Phase shaping
• Step 5: Coordinated PSS gain tuning (linear estimate)
• Step 6: Transient stability studies (updates step 5)
21
Application Example (1/6)
2012 Central Interconnected System of Chile (CISC)
(2x165MVA)(4x175MVA)
•
•
•
•
System is 2,400 km long. Peak (Minimum) load is 6,400 (4,650)MW.
Main generation is located in the Southern and Central areas.
The load is concentrated in the Central area.
Full nonlinear simulations in DIgSILENT’s PowerFactory v14: – 110 generators equipped AVR and turbine governor.
– 55 PSSs from different manufacturers and technologies.
22
Application Example (2/6)
• 60 power flow scenarios:
– Base case and N‐1 scenarios, various loading levels and maximized
power transfers.
– Typical electromechanical range of CISC: 0.5–3.5 Hz.
• Rotor speed mode shapes:
– The CISC has six interarea modes ranging from 0.50 to 1.15 Hz.
– Five modes are highly observable in the Guacolda plant: three
intraplant, one regional (Guacolda oscillating against Taltal) and one
interarea (the North versus the rest of the system).
Speed mode shape
Center - South mode
(0.87 Hz)
23
Application Example (3/6)
• Speed participation factors:
– Pointed the need to install PSSs at 7 additional generators.
– Guacolda’s PCB frequency range:
• Lower freq. end: set at 0.8 × 0.55 Hz = 0.44 Hz
• Higher freq. end: set by highest intra‐plant frequency (1.52 Hz).
Center - South mode (0.87Hz)
24
Application Example (4/6)
• From the synthetic system tests:
– change 6 PSSs from ΔP to ΔP‐ω;
– retune 35 PSSs;
– 14 existing PSSs had acceptable
performances and their tunings were
left unchanged.
• Electromechanical pole loci of the SMIB
family with (blue) and without (black) PSS.
• Root locus branches produced by varying
the PSS gain between 0 and 20 pu (red) for
9 synthetic system scenarios.
• PSS phase compensation example:
– nearly ideal at low PCB frequencies;
– some undercompensation at higher
frequencies.
Guacolda U3
25
Application Example (5/6)
• PSSs’ gain coordination: linear coordination achieved in a few iterations
in engineer’s interactive process (analysing critical eigenvalue shifts
produced by simultaneous changes to PSS gains).
• All the inter‐area modes
were
shifted
almost
horizontally to the left
when PSSs are enabled.
• Minimum damping ratio
levels kept above 10 and
5 % (grid code values) for
base case and N–1
network
scenarios,
respectively, through PSS
stabilization.
26
Application Example (6/6)
• Transient stability runs considering various disturbances:
– PSS performance extensively verified for large disturbances;
– PSS phase compensation left unchanged, but PSS gain, filters and
non‐linear blocks were added or modified when necessary;
– Sample result: fault simulation in the presence and absence of PSSs
at Guacolda.
North-System mode is notably damped out
27
Conclusions for Part I
• The synthetic system PSS phase tuning method was systematically
applied to 62 generators of a 110‐generator power system:
– Grid‐code compliant small and large‐signal performances were
achieved;
– The technical report for this study, with the proposed PSS
settings and nonlinear simulation results, was delivered to the
Chilean system operator;
• Results confirmed the effectiveness of applying the synthetic
system method to selected generators in a multimachine power
system, for the phase shaping of their PSSs followed by the
coordinated tuning of their gains, to damp all the local and inter‐
area modes.
• The methodology appears to be applicable to generic power
systems.
28
Part II:
Equivalents for Transmission Networks Containing Distributed Parameter Lines
Outline of Part II
• Introduction;
• Modeling electrical network components in the formulations
Descriptor System (DS) and Y(s) matrix;
• Distributed parameter transmission line model for Y(s) matrix
• The Sequential MIMO Dominant Pole Algorithm (SMDPA) for
computing the dominant poles and residue matrices associated
with MIMO TFs of infinite systems;
• Performance of a multi‐bus equivalent (MIMO ROM) for a
transmission network with distributed parameter lines (poles
computed by SMDPA).
Introduction to Part II (1/2)
• Modal Analysis
– Involves the calculation of the system matrix, its poles & zeros
and their sensitivities to system parameters;
– Provides system structural information: mode shapes,
participation factors, TF dominant poles, reduced order models;
– Matrix models are used for the study of different power system
phenomena:
• Eletromechanical transients
(Algebraic network modeling, R+jX);
• Subsynchronous resonance
(Lumped R‐L‐C
modeling);
• Harmonic performance;
• Electromagnetic Transients.
(High‐frequency network modeling,
all
transmission
lines
having
distributed parameters)
dynamic
network
Introduction to Part II (2/2)
• High Frequency Modeling of Electrical Networks
– 3 formulations: State Space (SS), Descriptor Systems (DS) and Y(s) matrix;
– The distributed parameter nature of transmission lines (TL) can be
modeled by transcendental functions having infinite poles – Infinite
systems;
– Infinite systems are neatly modeled by the Y(s) matrix formulation;
– Finite approximations of infinite systems can be modeled in the SS and
DS formulations, where TLs are represented by cascaded RLC circuits;
– Various NLA methods exist to efficiently compute ROMs for large scale DS
models;
– A main disadvantage of the Y(s) matrix formulation is the inexistance of
robust and efficient algorithms for the computation of the poles and
residue matrices of multivariable TFs.
Descriptor Systems (1/5)
•
Basic equations: T x& (t ) = A x(t ) + B u(t )
y (t ) = CT x(t ) + D u(t )
•
•
•
The components of the system are described by first‐order
ordinary differential equations and algebraic equations as well;
The Kirchhoff Law of Currents for each individual node of the
network is then added to these equations, to define the
connection among the various existing system components;
The DS model is a generalization of the SS model and leads to a
simpler and more efficient computer implementation.
Descriptor Systems (2/5)
™
Transfer Functions
TF MIMO →
TF SISO →
™
H (s ) = CT (s T − A ) B + D
−1
H (s ) =
y (s ) T
= c (s T − A )−1 b + d
u (s )
Frequency Response
H ( jω) = cT ( jω T − A )−1 b + d
™
Time Response (trapezoidal rule of integration)
⎛ 2
⎞
⎛ 2
⎞
⎜ T − A ⎟ x(t + Δt ) = ⎜ T + A ⎟ x(t ) + B[ u(t ) + u(t + Δt ) ]
⎝ Δt
⎠
⎝ Δt
⎠
y (t + Δt ) = CT x(t + Δt ) + D u(t + Δt )
Descriptor Systems (3/5)
RLC Series
L
dikj
dt
Voltage Source
RLC Parallel
Lf
= −vC − R ikj + vk − v j
C
dvC
= ikj
dt
dv
1
C C = −iL − vC + ikj
dt
R
di
L L = vC
dt
vk − v j − vC = 0
di f
dt
= − R f i f − vk + v j + v f
Kirchhoff Current Law
Node k
→
∑i
mk
=0
m∈Ω
Ω → Nodes connected to k
Descriptor Systems (4/5)
Matlab script vs PSCAD Validation
Parameters for 3‐bus system
v1
Input → i2 = 1 pu
i2
Output → v1 (pu)
Nominal Frequency: 50 Hz
Nominal Voltage: 20 kV
MVA base: 10 MVA
Descriptor Systems (5/5)
Matlab script vs PSCAD Validation
Voltage (pu)
Voltage at Bus # 1 following a step in the current injected in Bus 2
Time (ms)
⎡C1
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
(7) → ⎢⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
(13) → ⎢
⎢
⎣
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥d
⎥ dt
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
L2
C2
1
⎡
⎢− 1
⎢
1
⎢
⎢
− 1 − 1 R2
⎢
⎢
−1
⎢
⎢
(7) → ⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
−1
⎢
(13) → ⎢
⎢
⎣
L3
C3
L12
L13
Lf
⎡ vC1 ⎤
⎢ ⎥
⎢ i10 ⎥
⎢ iL ⎥
⎢ 2⎥
⎢vC 2 ⎥
⎢i ⎥
⎢ 20 ⎥
⎢ iL3 ⎥
⎢v ⎥
⎢ C3 ⎥ =
⎢ i30 ⎥
⎢ ⎥
⎢ i12 ⎥
⎢ i13 ⎥
⎢ ⎥
⎢ if ⎥
⎢v ⎥
⎢ 1 ⎥
⎢ v2 ⎥
⎢v ⎥
⎣ 3 ⎦
1
1
− 1 − 1 R3
1
−1
− R12
− R13
−1
1
−1
−1
−1
1
− Rf
1
C3
dvC3
1
= − iL3 − vC3 + i30 (7 )
dt
R3
0 = −i20 + i12 + i2 (13)
⎤
⎤ ⎡ vC1 ⎤ ⎡
⎥
⎥
⎢
⎢
⎥
i
1
⎥
⎥ ⎢ 10 ⎥ ⎢
⎥
⎥ ⎢ iL2 ⎥ ⎢
⎥
⎥⎢ ⎥ ⎢
⎥
⎥ ⎢vC2 ⎥ ⎢
⎥
⎥
⎢
⎢
⎥
i20
1
⎥
⎥⎢ ⎥ ⎢
⎥ v
⎥ ⎢ iL3 ⎥ ⎢
⎥⎡ f⎤
⎥ ⎢vC ⎥ ⎢
⎥ ⎢ i2 ⎥
⎥ ⎢ 3⎥+⎢
1 ⎥ ⎢ i30 ⎥ ⎢
⎥⎢ ⎥
⎥ ⎢⎣ i3 ⎥⎦
⎥ ⎢i ⎥ ⎢
1 −1
12
⎥
⎥⎢ ⎥ ⎢
⎥
1
− 1⎥ ⎢ i13 ⎥ ⎢
⎥
⎥⎢i ⎥ ⎢
−1
⎥
⎥ ⎢ f ⎥ ⎢1
⎥
⎥ ⎢ v1 ⎥ ⎢
⎥
⎥⎢ ⎥ ⎢
⎥ ⎢ v2 ⎥ ⎢ 1 ⎥
⎥ ⎢ v3 ⎥ ⎢⎣
1⎥⎦
⎦⎣ ⎦
Descriptor System Matrices (2/2)
•
Let us consider the nodal voltages as the output variables:
⎡ v1 ⎤ ⎡
⎢v ⎥ = ⎢
⎢ 2⎥ ⎢
⎢⎣v3 ⎥⎦ ⎢⎣
⎡ vC1 ⎤
⎢i ⎥
⎢ 10 ⎥
⎢ i L2 ⎥
⎢ ⎥
⎢v C 2 ⎥
⎢ i20 ⎥
⎢ ⎥
⎢ iL3 ⎥
1
⎤⎢ ⎥ ⎡
vC 3
⎥
⎥+⎢
⎢
1 ⎥
⎢ i30 ⎥ ⎢
1⎥⎦ ⎢ ⎥ ⎢⎣
⎢ i12 ⎥
⎢ i13 ⎥
⎢ ⎥
⎢ if ⎥
⎢v ⎥
⎢ 1⎥
⎢ v2 ⎥
⎢v ⎥
⎣ 3⎦
⎤ ⎡v f ⎤
⎥ ⎢i ⎥
⎥⎢ 2⎥
⎥⎦ ⎢⎣ i3 ⎥⎦
Y(s) Matrix Formulation
• Basic equations: Y ( s ) x (s ) = B u (s )
y (s ) = CT x(s ) + D u(s )
• Elements
• Diagonal yii: Summation of the all elementary admittances
connected to node i;
• Off‐diagonal yij: negative value of summation of all elementary
admittances connected between nodes i and j;
• Voltage sources are modeled by additional equations;
• The derivative of Y(s) with respect to s, for the computation of the
system poles, is automatically built by coding simple rules that are
similar to those used for building Y(s).
Y(s) Matrix Formulation
™
Transfer Function
MIMO TF →
SISO TF →
™
H(s ) = CT Y(s ) B + D
−1
H (s ) =
y (s ) T
= c Y(s )−1 b + d
u (s )
Frequency Response
H ( jω) = cT Y( jω)−1 b + d
Y(s) Matrix Formulation – Basic Elements
Series RLC
yseries =
dyseries
ds
Parallel RLC
1
1
R+s L+
sC
1
2
s
C
=
2
⎛
1 ⎞
⎜⎜ R + s L +
⎟⎟
sC⎠
⎝
−L+
Voltage Source
n
y parallel
1 1
= + + sC
R sL
dy parallel
ds
1
=C− 2
s L
∑y
j =1
kj
v j − i fk = 0
vk + z f k i f k = v f k
where:
z fk = R fk + s L fk
dz f k
ds
= L fk
Y(s) Matrix – 3‐bus System Equations
⎡ y11
⎢y
⎢ 21
⎢ y31
⎢
⎣1
y12
y13
y22
0
0
y33
0
0
− 1⎤ ⎡ v1 ⎤ ⎡0
0 ⎥ ⎢v2 ⎥ ⎢1
⎥⎢ ⎥=⎢
0 ⎥ ⎢v3 ⎥ ⎢0
⎥⎢ ⎥
z f ⎦ ⎣i f ⎦ ⎢⎣0
0 0⎤
⎡ i2 ⎤
⎥
0 0 ⎢ ⎥
⎥ i3
1 0⎥ ⎢ ⎥
⎥ ⎢⎣v f ⎥⎦
0 1⎦
Y (s )
⎞
⎞ ⎛
⎛
1
1
⎟⎟
⎟⎟ + ⎜⎜
y11 = s C1 + ⎜⎜
R
s
L
R
s
L
+
+
⎝ 12
12 ⎠ ⎝ 13
13 ⎠
⎛
⎞
1
⎟⎟
y13 = y31 = −⎜⎜
+
R
s
L
⎝ 13
13 ⎠
Voltage Source
LCK
y11 v1 + y12 v2 + y13 v3 − i f = 0
LTK
v1 + z f i f = v f
Y(s) Matrix – 3‐bus System Equations
0 ⎤
⎡ dy11 ds dy12 ds dy13 ds
⎢
0
0 ⎥
dY(s ) ⎢dy21 ds dy22 ds
⎥
=
0
0 ⎥
dy33 ds
⎢ dy31 ds
ds
⎢
⎥
0
0
0
dz
ds
f
⎣
⎦
⎤
⎡
⎤ ⎡
dy11
− L12
− L13
= C1 + ⎢
+
2⎥ ⎢
2⎥
ds
(
)
(
)
R
+
s
L
R
+
s
L
⎣ 12
⎦ ⎣ 13
⎦
12
13
dy13 dy31
L13
=
=
ds
ds (R13 + s L13 ) 2
Fonte de Tensão
dz f
ds
=
d (R f + s L f )
ds
= Lf
Y(s) Matrix – 3‐bus System Equations
⎡ y11
⎢y
⎢ 21
⎢ y31
⎢
⎣⎢ 1
y12
y22
0
y13
0
y33
0
0
− 1⎤ ⎡ v1 ⎤ ⎡0
0 ⎥⎥ ⎢⎢v2 ⎥⎥ ⎢1
=⎢
0 ⎥ ⎢ v3 ⎥ ⎢0
⎥⎢ ⎥ ⎢
z f ⎦⎥ ⎣⎢i f ⎦⎥ ⎣0
0 0⎤
⎡ i2 ⎤
⎥
0 0⎥ ⎢ ⎥
i3 ⎥
⎢
⎥
1 0
⎥ ⎢⎣v f ⎥⎦
0 1⎦
⎡ v1 ⎤
⎡ v1 ⎤ ⎡1 0 0 0⎤ ⎢ ⎥ ⎡0 0 0⎤ ⎡ i2 ⎤
⎢v ⎥ = ⎢0 1 0 0 ⎥ ⎢v 2 ⎥ + ⎢0 0 0 ⎥ ⎢ i ⎥
⎥⎢ 3⎥
⎢ 2⎥ ⎢
⎥ ⎢ v3 ⎥ ⎢
⎢⎣v3 ⎥⎦ ⎢⎣0 0 1 0⎥⎦ ⎢ ⎥ ⎢⎣0 0 0⎥⎦ ⎢v f ⎥
⎣ ⎦
⎢⎣i f ⎥⎦
A compact description
Y (s ) x(s ) = B u (s )
y (s ) = C x(s ) + D u (s )
MIMO System
Y(s) matrix‐ Freq. response results for 3‐bus system
DS compared to Y(s) Matrix:
Transfer Impedance
buses 1 & 2
between
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
1
= yc csch (γ l )
z
z = zc sinh (γ l )
⎛γl⎞
y = yc tanh⎜ ⎟ = yc [coth (γ l ) − csch (γ l )]
⎝ 2 ⎠
yc =
1
⎡ik ⎤ ⎡⎢ y +
z
⎢ ⎥=⎢
⎢ ⎥ ⎢ 1
−
⎣i j ⎦
⎣ z
ys = y +
Y (s )
Z (s )
1 ⎤ v
⎡ k ⎤ ⎡ ys
⎥
z ⎢ ⎥=⎢
⎥
1 ⎥ ⎢v ⎥ ⎢
y+
j⎦
⎣ − ym
⎣
z⎦
−
zc =
1
yc
γ = Z (s ) Y (s )
− ym ⎤ ⎡vk ⎤
⎥⎢ ⎥
⎥⎢ ⎥
y s ⎦ ⎣v j ⎦
1
= yc [coth (γ l ) − csch (γ l )] + yc csch (γ l )
z
ys = yc coth (γ l )
ym = yc csc h (γ l )
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
ys = yc coth (γ l )
dy s dyc
dγ
=
coth (γ l ) − yc
l csch (γ l )
ds
ds
ds
ym = yc csc h (γ l )
dym dyc
dγ
=
csch (γ l ) − yc
l csch (γ l ) coth (γ l )
ds
ds
ds
yc =
γ =
Y (s )
Z (s )
Z (s ) Y (s )
dyc
1
=
ds 2 γ
⎡ dY
2 dZ ⎤
⎢⎣ ds − yc ds ⎥⎦
dZ ⎤
dγ 1 ⎡ dY
=
+
Y
Z
ds ⎥⎦
ds 2 γ ⎢⎣ ds
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
dY
= C1
ds
Y = s C1
Z =Z
(e )
+Z
Z (e ) = s L1(e )
(i )
+Z
(g )
dZ dZ (e ) dZ (i ) dZ ( g )
=
+
+
ds
ds
ds
ds
dZ (e )
= L1(e )
ds
C1 , L1(e ) → Positive sequence capacitance & indutance, computed by
matrix reduction considering ideal conductor and soil.
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
Z (i ) =
k (s ) =
k (s ) n ( s )
ns m(s )
sμ 1
σ 2 π re
n(s ) = I 0 (ρ1 ) K1 (ρ0 ) + I1 (ρ0 ) K 0 (ρ1 )
m(s ) = I1 (ρ1 ) K1 (ρ0 ) + I1 (ρ0 ) K1 (ρ1 )
ρ0 = ri
sμσ
ρ1 = re
sμσ
I0, I1 → Modified Bessel functions of first kind for integer orders 0 & 1,
respectively.
K0, K1 → Modified Bessel functions of first kind for integer orders 0 & 1,
respectively.
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
Z (i ) =
k (s ) n(s )
n s m (s )
dn(s )
dm(s ) ⎤
⎡
(
)
(
)
−
m
s
n
s
dZ
1 ⎢ dk (s ) n(s )
ds
ds ⎥
= ⎢
+ k (s )
⎥
2
ds
ns ⎢ ds m(s )
m (s )
⎥
⎣
⎦
(i )
dk (s ) k (s )
=
ds
2s
dn(s ) dI 0 (ρ1 ) dρ1
dK (ρ ) dρ0 dI1 (ρ0 ) dρ0
dK (ρ ) dρ1
K1 (ρ0 ) + I 0 (ρ1 ) 1 0
+
K 0 (ρ1 ) + I1 (ρ0 ) 0 1
=
ds
dρ1 ds
dρ 0
ds
dρ0 ds
dρ1
ds
dm(s ) dI1 (ρ1 ) dρ1
dK (ρ ) dρ0 dI1 (ρ0 ) dρ0
dK (ρ ) dρ1
K1 (ρ0 ) + I1 (ρ1 ) 1 0
−
K1 (ρ1 ) − I1 (ρ0 ) 1 1
=
ds
d ρ1 d s
dρ 0 d s
dρ0 ds
dρ1 d s
dρ 0
ds
=
ρ0
2s
dρ1 ρ1
=
ds 2 s
dI 0 (ρ)
= I 1 (ρ)
dρ
dI1 (ρ) I 0 (ρ) + I 2 (ρ)
=
2
dρ
K (ρ) + K 2 (ρ )
dK1 (ρ )
=− 0
dρ
2
dK 0 (ρ )
= − K 1 (ρ )
dρ
Matriz
Y(s) – Modelagem
de Redes
Y(s) Matrix –
Distrib. Param. Transmission Line
⎡
μ0 ⎢
(g )
Z =s
⎢
6π⎢
⎢⎣
(g )
μ0 ⎡
dZ
=
⎢
ds
6π⎢
⎣
p=
1
3
∑
i =1
3
∑
i =1
⎛ Hi + p ⎞
⎟⎟ −
ln⎜⎜
⎝ Hi ⎠
3
3
i =1
j =1
j ≠i
∑∑
⎛H + p⎞
s
dp
⎟⎟ +
−
ln⎜⎜ i
⎝ H i ⎠ H i + p ds
s μ (σ + s ε )
⎛
⎜
ln⎜
⎜
⎝
3
3
i =1
j =1
j ≠i
∑∑
(H i + H j + 2 p ) +
(H i + H j )2 + Dij2
2
⎛
⎜
ln⎜
⎜
⎝
μ (σ + 2 s ε ) p 3
dp
=−
2
ds
(H i + H j + 2 p )2 + Dij2
dHˆ ij ( p )
2 (H i + H j + 2 p )
=
ds
(H i + H j + 2 p )2 + Dij2
Hˆ ij =
dp
ds
Dij2
⎤
⎞⎥
⎟
⎟⎥
⎟⎥
⎠⎥
⎦
⎤
⎞
ˆ
s dH ij ( p )⎥
⎟
+
ˆ ( p ) ds ⎥
2
2 ⎟
H
⎟
H i + H j + Dij ⎠
ij
⎥
⎥⎦
(
Hˆ ij
)
Y(s) Matrix – Distrib. Param. Transmission Line
⎡ y11
⎢y
21
Y(s ) = ⎢
⎢ 0
⎢
⎢⎣ 1
y12
y22
y32
0
0
y23
y33
0
− 1⎤
0 ⎥⎥
0⎥
⎥
z f ⎥⎦
y33 (s ) =
52
1
1
+ s C3 + Y +
R3 + s L3
Z
yc (s ) coth[γ (s ) l ]
0
0 ⎤
⎡ dy11 ds dy12 ds
⎢
0 ⎥⎥
dY(s ) ⎢dy21 ds dy22 ds dy23 ds
=
− L3
d
⎢ 0
dy32 ds dy33 ds
0 ⎥ dy33 =
ds
{ yc (s ) coth[γ (s ) l ] }
+
C
+
3
2
⎢
⎥ ds
ds
(R3 + s L12 )
dz f ds ⎥⎦
0
0
⎢⎣ 0
Y(s) Matrix – Distrib. Param. Transmission Line
⎡ y11
⎢y
⎢ 21
⎢ 0
⎢
⎣⎢ 1
y12
0
y22
y32
0
y23
y33
0
− 1⎤ ⎡ v1 ⎤ ⎡ 0 ⎤
0 ⎥⎥ ⎢⎢v2 ⎥⎥ ⎢ 0 ⎥
= ⎢ ⎥ vf
0 ⎥ ⎢v3 ⎥ ⎢ 0 ⎥
⎥⎢ ⎥ ⎢ ⎥
z f ⎦⎥ ⎣⎢i f ⎦⎥ ⎣ 1 ⎦
⎡ v1 ⎤
⎢v ⎥
2
v3 = [0 0 1 0] ⎢ ⎥
⎢v3 ⎥
⎢ ⎥
⎢⎣i f ⎥⎦
Y (s ) x = b v f
v3 = ct Y(s ) b v f
−1
b t = [0 0 0 1]
v3 = ct x
G31 (s ) =
(
v3
v1 = v f
)
= ct Y(s ) b
ct = [0 0 1 0]
−1
53
Sequential MIMO Dominant Pole Algorithm (SMDPA)
SMDPA – Fundamental Concepts (1/3)
™
Direct Matrix
H (s ) = CT Y(s )−1 B + De
−1
Di
D = Di + D e
™
Partial Fraction Expansion
H(s ) =
∞
∑
i =1
Ri
+ s K + Di + D e
s − λi
dH (s )
s → ∞ ds
K = lim
Di = lim H(s ) − s K − De
s →∞
SMDPA – Fundamental Concepts (2/3)
™
Strictly Proper part of H(s)
∞
H (s ) = ∑
i =1
Ri
+ sK + D
s − λi
ˆ (s ) = H (s ) − s K − D
H
™
Dominant Pole
Let λk = αk + j βk be a pole with an associated residue matrix Rk, then:
ˆ ( j β ) = − Rk +
H
k
αk
∞
∑
i =1
i≠k
Ri
j βk − λ i
The pole λk will be dominant in H(s) if the magnitude of ( ||Rk||2 / |αk| )
is sufficiently large so as to cause a peak in the plot of σmax[Ĥ(jω)] in
the close neighborhood of the frequency βk.
SMDPA – Fundamental Concepts (3/3)
™
Reduced Order Model (ROM)
H (s ) ≅ H N ( s ) =
™
N
∑
i =1
λi ∈ Ω
Ri
+ sK +D
s − λi
Ω → Set of N dominant poles and
associated residue matrices.
MIMO ROM Deviation TF
H (s ) = H (s ) − H N (s ) = H (s ) − s K − D −
™
Norm of MIMO ROM Deviation TF
ε MOR ( jω) = σ max [H ( jω)]
N
∑
i =1
λi∈ Ω
Ri
s − λi
SMDPA – Newton Method (1/3)
™
™
The set of dominant poles of H(s) may be efficiently
computed only when eliminating from H(s):
• The N previously computed poles (deflation);
• Matrices K and D.
The Newton method should therefore be applied to the
MIMO ROM deviation TF:
H (s ) = H (s ) − s K − D −
N
∑
i =1
λi∈ Ω
Ri
s − λi
SMDPA – Newton Method (2/3)
™
Pole of H (s )
[
lim μ min H (s )
s →λ
™
−1
] =0
μmim → minimum eigenvalue of H (s )−1
vmin, wmin → eigenvectors associated with μmin.
Newton equationing
s (k +1) = s (k ) + Δs
f (s ) = μ min (s ) = 0
Δs (k ) =
1
(k )
( )
( )
* dH ( s )
v min s (k )
μ min ⎛⎜⎝ s (k ) ⎞⎟⎠ w min s (k )
ds
dH (s ) dH (s ) dH N (s )
=
−
ds
ds
ds
d H (s )
dY(s )
= − XTC (s )
X B (s )
ds
ds
dH N (s )
=−
ds
N
Ri
∑
(s − λ )
i
=1
λi∈ Ω
2
i
+K
(μmim, vmin, wmin)
function eig of Matlab.
Y (s ) X B (s ) = B
Y(s ) X C (s ) = C
T
The sparse solution of these
two
matrix
equations
sistemas require a single LU
factorization and various
solves.
SMDPA – Newton Method (3/3)
™
Determining the initial pole estimates
s (0 ) = 2 π j [ f1 L
f1
f2 f3
f4 f5
f6
f8 ]
f7 f8
SMDPA – Pole Residue Matrix (1/2)
ƒ Computation of the Residue Matrix for a Pole
™
Definition of the Integration Curve
jω
jω
λ
λ
C
P3
1
2πj
∫ H(s ) ds
C
P4
Real(s)
Real(s)
R=
P1 ≡ P5
P2
4 P( k +1)
∫ H(s ) ds = ∑ ∫ H(s ) ds
C
k =1 Pk
SMDPA – Pole Residue Matrix (2/2)
ƒ Legendre‐Gauss Method with Error Control
Pk + ΔPk
Pk
Pk + (l − 1) ΔPk
Pk + 2 ΔPk
Pk + (mk − 1) ΔPk
Pk + l ΔPk
Pk +1P=
k +1mk ΔPk
s = Pk + [ l + 0.5 (ξ − 1) ] ΔPk
4
mk
k =1
4
ΔPk
H (s ) ds ≅
k =1 2
l =1 P + (l −1) ΔP
∫ H(s ) ds = ∑∑
C
Pk +l ΔPk
∑
∫
k
k
mk
M
wi J (ξi , l , Pk )
∑∑
l =1 i =1
J (ξ, l , Pk ) = H ( Pk + [ l + 0.5 (ξ − 1) ] ΔPk
1
R≅
2πj
4
ΔPk
k =1 2
∑
mk
M
∑∑ w J(ξ , l, P )
i
i
k
l =1 i =1
1 ΔPk
Rk ≅
2πj 2
mk
)
M
∑∑ w J(ξ , l , P ) , k = 1, K, 4
i
i
k
l =1 i =1
• Determining the Error
q
mk = 2 , q = 0, 1, 2, K
ε = R (q +1) − R (q )
k
k
k
(q +1)
2
Rk
2
ε=
4
ε k ≤ 4 ε max
∑
k =1
SMDPA – 34‐bus Test System Results (1/4)
ƒ The 34-bus test system with 25 distributed parameter TLs
Element
Quant.
Buses
34
TLs
25
Branches
12
Trafos
16
Loads
16
Generators
10
SMDPA – 34‐bus Test System Results (2/4)
Most dominant pole spectrum for 34‐bus system (2x2) TF. Color code identifies
iterations required for SMDPA convergence from the initial set of 37 estimates
SMDPA – 34‐bus Test System Results (3/4)
Infinito
MOR − 151
SMDPA – 34‐bus Test System Results (4/4)
Conclusions for Part II (1/2)
• Electrical network modeling with its RLC series and paralell components,
current and voltage sources, in the DS and Y(s) matrix formulations;
• Development of the first reliable Newton algorithm for computing the
dominant poles of SISO and MIMO TFs of infinite systems (SMDPA). The
method’s reliability comes from the very effective pole deflation
procedure and the accurate computation of the pole residue matrices;
– The residue is numerically computed as the path integral around the
pole which was here obtained by the Legendre‐Gauss quadrature
method.
• The modeling accuracy of the DS and Y(s) formulations was verified by
the close matching between their simulation results and those obtained
with ATP or PSCAD for various test systems.
Conclusions for Part II (2/2)
• SMDPA yields high fidelity ROMs over a specified frequency window,
for use as equivalents in transmission network electromagnetic
transient studies.
• Multi‐bus ROMs produced by SMDPA directly from infinite system
models are a more practical option than LMA (Linear Matrix
Approximation) models;
• In the attempt to obtain accuracy over a wider frequency range, LMA
models may soon reach very large dimensions and present severe
numerical stiffness.
• Use of SMDPA in other areas of engineering, physics and
mathematics is yet to be exploited.
References
•
•
•
•
•
•
DE MARCO, F.J., MARTINS, N., FERRAZ, J.C., An Automatic Method for Power System Stabilizers Phase Compensation Design. IEEE Transactions on Power Systems, USA, Vol. 28, Issue: 2, p. 997‐1007, May 2013.
VARRICCHIO, S.L., FREITAS, F.D., MARTINS, N., VELIZ, F.C., Computation of Dominant Poles and Residue Matrices for Multivariable Transfer Functions of Infinite Power System Models. IEEE Transactions on Power Systems, USA, Vol. 30, Issue:3, p. 1131 ‐ 1142, May 2015.
GOMES Jr, S.; MARTINS, N.; PORTELA, C. Sequential Computation of Transfer Function Dominant Poles of s‐Domain System Models. IEEE Transactions on Power Systems, USA, Vol. 24, No. 2, p. 776‐784, May 2009.
GOMES Jr, S.; PORTELA, C.; MARTINS, N., Detailed Model of Long Transmission Lines for Modal Analysis of ac Networks. Proc. International Conference on Power System Transients, Rio de Janeiro, Brazil, June 2001. ROMMES, J. and MARTINS, N., Efficient Computation of Multivariable Transfer Function Dominant Poles Using Subspace Acceleration. IEEE Transactions on Power Systems , USA, Vol. 21, No. 4, p. 1471‐1483, November 2006.
MARTINS, N., BOSSA, T.H. S. A Modal Stabilizer for the Independent Damping Control of Aggregate Generator and Intraplant Modes in Multigenerator Power Plants. IEEE Transactions on Power Systems, USA, Vol. 29, Issue: 6, p. 2646 ‐ 2661, November 2014. 70
Thank you!
Nelson Martins