Sync Generator Rotor Position Estimation

DEGREE PROJECT IN ELECTRICAL ENGINEERING,
SECOND CYCLE, 30 CREDITS

STOCKHOLM, SWEDEN 2016
Synchronous generator rotor

position estimation from stator
and rotor terminal currents and
voltages
WILLE KIURU
KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ELECTRICAL ENGINEERING
Abstract
This thesis presents an alternative method for measuring the rotor position of a synchronous
generator, without the use of sensors to measure the rotor position. The method uses the
current and voltage measurements on the stator and rotor terminals from the protection
relays. The first step is to investigate how the rotor position can be estimated based
on the available measurement data. It was found out that the rotor position couldn’t
be estimated during dynamic state operation, when directly solving it from measurement
data. This is due to the presence of the field and d-axis damper winding current dynamics,
including the q-axis damper winding current. In order to achieve this, multiple steps
consisting of different methods for estimating the rotor position were utilized. The first
step is to estimate the rotor position during steady state. The second step estimates the
rotor position iteratively based on initial values. This method is useful during a dynamic
state operation. The initial values are obtained from the first step at the start of a dynamic
state operation. The final step predicts the rotor position based on previous data. It is
used during cases when the second step fails to estimate the rotor position accurately.
Due to insufficient data for estimating the rotor position during saturation, the effects of
saturation for the simulated synchronous machine were not considered.
Sammanfattning
Detta projekt presenterar en alternativ metod för estimering av rotorvinkeln för en synkro-
nmaskin med utpräglade poler, utan användning av sensorer. Metoden som används är
att estimera rotorvinkeln genom uppmätt data av strömmar och spänningar från stator
och rotor terminaler inom skyddsreläer. Första steget är att undersöka ifall rotorvinkeln
kan bestämmas direkt ur mätdata. Det visade sig att rotorvinkeln kunde inte estimeras
vid en dynamisk fas, genom att lösa rotorvinkeln direkt från mätdata. Detta på grund
av strömmderivatan i fält och dämplindningen i d-led, samt strömmen i dämplindningen i
q-led. För att uppnå detta, flera olika steg baserad på olika metoder användes för att es-
timera rotorvinkeln. Det första steget går ut på att estimera rotorvinkeln i ett stabilt läge.
Det andra steget estimerar rotorvinkeln iterativt baserad på initialvärden. Denna metod
är användbar vid ett icke-stabilt läge. Initialvärdena erhålls från steg ett vid början av ett
icke-stabilt läge. Det slutliga steget förutsäger rotorvinkeln baserad på tidigare data. Den
används då steg två misslyckas att estimera rotorvinkeln noggrant. På grund av otillräcklig
data för att estimera rotorvinkeln under mättnad så togs inte mättnadsfenomen till hänsyn
för den simulerade synkronmaskinen.
Acknowledgement
I would like to thank my supervisor Nathaniel Taylor, researcher at KTH for the Electro-
magnetic Engineering departement and Jianping Wang, Senior Principal Scientist at ABB
for guiding me during the course of the master thesis.
List of symbols
ia Phase a current [A]

ib Phase b current [A]
ic Phase c current [A]
id d-axis current [A]
iq q-axis current [A]
i1d d-axis damper winding current [A]
i1q q-axis damper winding current [A]
if d Field current [A]
ea Instantaneous phase a voltage [V]
eb Instantaneous phase b voltage [V]
ec Instantaneous phase c voltage [V]
ed d-axis voltage [V]
eq q-axis voltage [V]
ef d Field voltage [V]
Ψa Phase a flux linkage [Wb]
Ψb Phase b flux linkage [Wb]
Ψc Phase c flux linkage [Wb]
Ψd d-axis flux linkage [Wb]
Ψq q-axis flux linkage [Wb]
Ψ1d d-axis damper winding flux linkage [Wb]
Ψ1q q-axis damper winding flux linkage [Wb]
Ψf d Field winding flux linkage [Wb]
Ra Stator resistance [Ω]
R1d d-axis damper winding resistance [Ω]
R1q q-axis damper winding resistance [Ω]
Rf d Field winding resistance [Ω]
laa Self-inductance of stator winding a [H]
lbb Self-inductance of stator winding b [H]
lcc Self-inductance of stator winding c [H]
lab Mutual inductance between stator winding a and b [H]
lbc Mutual inductance between stator winding b and c [H]
lca Mutual inductance between stator winding c and a [H]
laf d Mutual inductance between the field winding and stator [H]
la1d Mutual inductance between the d-axis damper winding and stator winding a [H]
la1q Mutual inductance between the q-axis damper winding and stator winding a [H]
Lad d-axis mutual inductance [H]
Laq q-axis mutual inductance [H]
Ll Stator leakage inductance [H]
Lf f d Self-inductance of field winding [H]
L11d Self-inductance of d-axis damper winding [H]
L11q Self-inductance of q-axis damper winding [H]
Lf 1d Mutual inductance between field winding and d-axis damper winding [H]
Ld d-axis synchronous inductance [H]
Lq q-axis synchronous inductance [H]
L0d d-axis transient inductance [H]
L0q q-axis transient inductance [H]
L00d d-axis subtransient inductance [H]
L00q q-axis subtransient inductance [H]
0
τd0 Open-circuit d-axis transient time constant [s]
00
τd0 Open-circuit d-axis subtransient time constant [s]
τd0 Short-circuit d-axis transient time constant [s]
τd00 Short-circuit d-axis subtransient time constant [s]
00
τq0 Open-circuit q-axis subtransient time constant [s]
τq00 Short-circuit q-axis subtransient time constant [s]
θ, θre Electrical rotor angle/position [rad]
ωr , ωre Electrical rotor angular velocity [rad/s]
Contents
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Purpose, aim and objective of the work . . . . . . . . . . . . . . . . . . . . 2
1.3 Outline of the report . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Theoretical background for the testing environment 3

2.1 Mathematical description of a synchronous machine . . . . . . . . . . . . . 3
2.1.1 Stator circuit equations . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.2 Stator self-inductances . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.3 Stator mutual inductances . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.4 Mutual inductance between stator and rotor windings . . . . . . . . 6
2.1.5 Rotor circuit equations . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Park’s transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Stator flux linkages in dq0 components . . . . . . . . . . . . . . . . . 8
2.2.2 Rotor flux linkages in dq0 components . . . . . . . . . . . . . . . . . 8
2.2.3 Stator voltage equations in dq0 components . . . . . . . . . . . . . . 9
2.3 Per unit representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.1 Per unit stator voltage and flux linkage equations . . . . . . . . . . . 9
2.3.2 Per unit rotor voltage and flux linkage equations . . . . . . . . . . . 10
2.3.3 Complete set of electrical equations in per unit . . . . . . . . . . . . 10
2.4 Derivation of machine parameters . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4.1 Classical expression for time constants . . . . . . . . . . . . . . . . . 13
2.4.2 Classical expression for transient and subtransient inductances . . . 13
2.5 Power system components . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5.1 ∆−Y step up transformer . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5.2 Transmission line model . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Rotor position estimation 17

3.1 Extended Electromotive Force based on the voltage concept . . . . . . . . . 17
3.1.1 Dynamic model of a salient-pole synchronous machine . . . . . . . . 17
3.1.2 Model Reconstruction Based on Voltage Concept . . . . . . . . . . . 18
3.1.3 Reconstructed salient-pole synchronous machine model . . . . . . . . 19
3.2 Dynamic state rotor position estimator . . . . . . . . . . . . . . . . . . . . . 20
3.2.1 Estimating d-axis current . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.2 Rotor position estimation . . . . . . . . . . . . . . . . . . . . . . . . 23
4 Designing the simulation environment and the estimator 25

4.1 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.1 Power system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.1.2 Synchronous generator . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.1.3 Output vectors from the synchronous generator . . . . . . . . . . . . 27
4.2 Estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2.1 Designing the EEMF subsystem . . . . . . . . . . . . . . . . . . . . 29
4.2.2 Designing the dynamic state rotor position estimator . . . . . . . . . 30
4.2.3 Calculating the rotor position estimation error . . . . . . . . . . . . 35
5 Results and discussion 36
5.1 Input data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.2 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2.1 Single-phase-to-ground fault . . . . . . . . . . . . . . . . . . . . . . . 38
5.2.2 Phase-to-phase fault . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.2.3 Three-phase-fault . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.3.1 Choice of power system model . . . . . . . . . . . . . . . . . . . . . 45
5.3.2 Choice of estimator model . . . . . . . . . . . . . . . . . . . . . . . . 46
6 Conclusion and future work 47

6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Appendix A 48
A.1 Per Unit System for Stator and Rotor Quantities . . . . . . . . . . . . . . . 48
A.2 Derived Standard parameters from time constants . . . . . . . . . . . . . . 48
Appendix B 49
B.1 Complementary results for P-G fault . . . . . . . . . . . . . . . . . . . . . . 49
B.2 Complementary results for P-P fault . . . . . . . . . . . . . . . . . . . . . . 50
B.3 Complementary results for 3φ fault . . . . . . . . . . . . . . . . . . . . . . . 51
Appendix C 54
C.1 Calculating id and i1d in Simulink . . . . . . . . . . . . . . . . . . . . . . . 54
C.2 Other Simulink systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
Appendix D 57
D.1 Magnetic saturation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
D.2 Measuring the field current for a static excitation system . . . . . . . . . . . 58
1 Introduction
1.1 Background
Sensor-less estimation of the rotor position has been an important application for different
types of electrical machines, including synchronous [1]–[7], asynchronous [8] and DC elec-
trical machines [9], [10]. The advantages for a sensor-less rotor position estimation is to
limit the amount of external components. This in turn will lead to reduced costs and better
system reliability [7]. The disadvantages are the higher requirement of control algorithms
and complicated electronics [11]. For this thesis, a synchronous generator is considered,
where the rotor position is to be estimated from current and voltage measurement data,
already available from the generator protection relay. The rotor field current and voltage
are assumed to be known as well, which provides additional information in order to esti-
mate the rotor position.
From a vast amount of available sources and documents, it is clear that sensor-less rotor
position estimation technology has been in development as early as the 1990s, where a
lot of research has been done to improve vector controlled synchronous PM (permanent
magnet) machines [1]–[7]. Due to the simplicity of PM magnet synchronous machines,
the methods presented in [1]–[7], [12], [13] are only practical for certain PM synchronous
machines, where the effects of damping are negligible. The focus in this master thesis is
to implement a sensor-less rotor position estimator for a non-PM wound-rotor salient pole
synchronous generator, providing power to a large power system. Therefore conventional
methods used in vector controlled PM machines to estimate the rotor position will not
be valid for certain cases, ranging from power swings to symmetrical and unsymmetrical
faults. A more complicated scheme is needed to account for these types of faults.
Before designing the rotor position estimator, the first step is to make sure there is a
proper environment for the estimator to be tested. The choice was to use Simulink (Mat-
lab version R2016a) to simulate a synchronous machine connected to a large power system.
To simplify such a complicated system, the power system is represented as an infinite bus,
i.e. an ideal voltage source with fixed frequency. From this assumption, the system is
represented as a single-machine-infinite-bus (SMIB) system.
Figure 1.1: Layout of the SMIB system.
The SMIB layout as shown in figure 1.1 is connected to a ∆−Y step-up transformer
(13.8/200 kV) which is connected to a transmission line consisting of two parallel trans-
mission lines connected to an infinite bus.
1
1.2 Purpose, aim and objective of the work
As mentioned previously, the purpose for a sensor-less rotor position estimation is to reduce
the system cost and to improve the reliability. Both factors are achieved due to one extra
component being omitted. The aim for this project was to design an accurate and practical
method for estimating the rotor position. This however did not fulfill all of the expectations
due to the inability to estimate the rotor position when the rotor core is saturated. The
objective was to use the available current and voltage measurements on the stator and
rotor terminals, to estimate the rotor position.
1.3 Outline of the report

The main content of this report starts in Section 2 with a detailed description of the theory
for the simulation environment, including the theoretical background for an synchronous
machine and the power system. The background theory for estimating the rotor position
is explained in Section 3, starting from a simpler EEMF (Extended Electromotive Force)
model, to include a system to account for higher dynamic cases. The estimator design is
explained in Section 4. Testing and the results are presented in Section 5, where the system
is subjected to different types of faults, ranging from symmetric to unsymmetrical faults.
The conclusion and suggested future work is presented in Section 6.
2
2 Theoretical background for the testing environment
2.1 Mathematical description of a synchronous machine
This chapter explains the step by step derivations of the synchronous machine parameters
and equations. A cross section of a salient pole synchronous machine is shown in figure
2.1, where the outer layer is referred as the stator, consisting of the three phase windings
connected to the power grid. The innermost part is referred to as the rotor, consisting of
a field winding and two damper windings on each side or axis. The field winding is used
to magnetize the rotor with the help of a separately added field current through a DC
voltage source or a static excitation system (which is provided in Appendix D). The main
concern is to express stator and rotor self- and mutual inductances. The problem is that
the mutual inductances between the rotor and stator are not constant, but vary depending
on the rotor position θ. For a salient pole the permeance (which is the inverse of magnetic
reluctance) is not uniform along the periphery as shown in figure 2.2. The term α in figure
2.2 is the angular distance from the d-axis along the periphery [14].
Figure 2.1: Three phase synchronous machine [14].
3
Figure 2.2: Permanence as a function of rotor position [14].
It will be shown that the mathematical expression of a synchronous machine can be rep-
resented with constant machine parameters.
All mathematical models described in this chapter are derived from [14].
2.1.1 Stator circuit equations

The stator circuit voltage equations are expressed as follows
dΨa
ea = − Ra ia = pΨa − Ra ia (2.1)
dt
eb =pΨb − Ra ib (2.2)
ec =pΨc − Ra ic (2.3)
Where the above equations represent the three phase voltage equations. The term p is
the differential operator ”d/dt”, while Ra is the armature resistance per phase. The flux
linkage for each phase is expressed as follows
Ψa = −laa ia − lab ib − lac ic + laf d if d + la1d i1d + la1q i1q (2.4)
Further on, the terms laa , lbb , lcc denotes the self-inductances of stator windings. The mu-
tual inductances between stator windings are denoted as lab , lbc , lca , while the mutual induc-
tances between stator and rotor windings are laf d , la1d , la1q (phase a). The self-inductances
of the rotor circuits are denoted as Lf f d , L11d , L11q . The terms laa , lbb , lcc , lab , lbc , lca , laf d , la1d ,
la1q all vary depending on the rotor position. The self-inductances of the rotor circuits how-
ever do not vary with rotor position. This is due to the constant permeance of the rotor
circuits, due to the cylindrical structure of the stator.
2.1.2 Stator self-inductances

The stator self-inductance laa is equal to the ratio of flux linking to phase a winding to the
current ia , which is proportional to the permeance, considering figure 2.2 the inductance
4
laa varies from having a maximum for θ = 0◦ and a minimum for θ = 90◦ , and so on.
The magnetomotive force mmf has its peak at the center of the phase a axis. The peak
amplitude of the mmf is equal to Na ia , where Na is the effective turns per phase. The peak
mmf components for d and q-axis are expressed as follows.
MMFad(peak) =Na ia cos θ (2.5)

MMFaq(peak) = − Na ia sin θ (2.6)
From the above equations, the air-gap fluxes per pole along the d and q-axis are as follows.
Φgad = (Na ia cos θ) Pd (2.7)

Φgaq = (−Na ia sin θ) Pq (2.8)
Where Pd and Pq are the d, respectively q-axis permeance coefficients.
The total air-gap flux linking phase a is
Φgaa =Φgad cos θ − Φgaq sin θ = Na ia (Pd cos2 θ + Pq sin2 θ)

(2.9)

Pd + Pq Pd − Pq
=Na ia + cos 2θ
2 2
The self inductance due to air-gap flux is

Na Φgaa Pd + Pq Pd − Pq
lgaa = = Na2 + cos 2θ = Lg0 + Laa2 cos 2θ (2.10)
ia 2 2
The leakage flux not crossing the air-gap is represented as Lal . The total self inductance
laa is obtained by adding the self-inductance due to the air-gap flux and Lal .
laa = Lal + lgaa = Lal + Lg0 + Laa2 cos 2θ = Laa0 + Laa2 cos 2θ (2.11)
Similarly the b and c phases are expressed as follows

4π
lbb = Laa0 + Laa2 cos 2θ − (2.12)
3

4π
lcc = Laa0 + Laa2 cos 2θ + (2.13)
3
5
2.1.3 Stator mutual inductances
The mutual inductance lab between the stator windings (any two stator winding) is calcu-
lated first by evaluating the air-gap flux Φgba , linking phase b when only phase a is excited.

2π 2π
Φgba =Φgad cos θ − − Φgaq sin θ −
3 3

2π 2π
=Na ia Pd cos θ cos θ − + Pq sin θ sin θ − (2.14)
3 3

Pd + Pq Pd − Pq 2π
=Na ia − + cos 2θ −
4 2 3
The mutual inductance due to the air-gap flux between phases a and b is

Na Φgba 1 2π
lgba = = − Lg0 + Lab2 cos 2θ − (2.15)
ia 2 3
Therefore the mutual inductance between phase a and b can be written as

2π
lab = lba = − Lab0 + Lab2 cos 2θ −
3
(2.16)
2π
= − Lab0 − Lab2 cos 2θ +
3
The mutual inductance between phase b-c and c-a are
lbc = lcb = − Lab0 − Lab2 cos(2θ − π) (2.17)

π
lca = lac = − Lab0 − Lab2 cos 2θ − (2.18)
3
2.1.4 Mutual inductance between stator and rotor windings

The mutual inductances between stator and rotor windings are expressed as follows
laf d =Laf d cos θ (2.19)

la1d =La1d cos θ (2.20)
π
la1q =La1q cos θ + = −La1q sin θ (2.21)
2
Equation (2.4) can now be expressed for all three phases as follows

2π
Ψa = − ia [Laa0 + Laa2 cos 2θ] + ib Lab0 + Lab2 cos 2θ +
3
h π i (2.22)
+ ic Lab0 + Lab2 cos 2θ − + if d Laf d cos θ
3
+ i1d La1d cos θ − i1q La1q sin θ
6

h π i 2π
Ψb = − ia Laa0 + Laa2 cos 2θ + − ib Lab0 + Lab2 cos 2 θ −
3 3

2π
+ ic [Lab0 + Lab2 cos (2θ − π)] + if d Laf d cos θ − (2.23)
3

2π 2π
+ i1d La1d cos θ − − i1q La1q sin θ −
3 3
h π i
Ψc = − ia Laa0 + Laa2 cos 2θ − − ib [Lab0 + Lab2 cos (2θ − π)]
3
4π 2π
+ ic Lab0 + Lab2 cos 2θ + + if d Laf d cos θ + (2.24)
3 3

2π 2π
+ i1d La1d cos θ + − i1q La1q sin θ +
3 3
2.1.5 Rotor circuit equations

The rotor circuit voltage equations are expressed as follows
ef d =pΨf d − Rf d if d (2.25)
e1d =pΨ1d − R1d i1d = 0 (2.26)
e1q =pΨ1q − R1q i1q = 0 (2.27)
Since the damper windings are short circuited, both damper winding voltages e1d and e1q
are zero.
The rotor circuit flux linkages are expressed as follows

2π 2π
Ψf d =Lf f d if d + Lf 1d i1d − Laf d ia cos θ + ib cos θ − + ic cos θ + (2.28)
3 3

2π 2π
Ψ1d =Lf 1d if d + L11d i1d − La1d ia cos θ + ib cos θ − + ic cos θ + (2.29)
3 3

2π 2π
Ψ1q =L11q i1q − La1q ia sin θ + ib sin θ − + ic sin θ + (2.30)
3 3
2.2 Park’s transformation

Equations (2.2)-(2.3), (2.22)-(2.24) and (2.25)-(2.30) contains inductance terms that are
dependent on the rotor angle, which is dependent on time. In order to express the men-
tioned equations in terms of constant inductance parameters, all three phase components
7
are expressed as their equivalent rotor reference frame (dq0-axis) components. This is done
bu using the Park’s transformation, as shown in equation (2.31).
fd cos θ cos(θ − 2π
3 ) cos(θ + 2π
3 ) fa
    
f  2 − sin θ − sin(θ − 2π ) − sin(θ + 2π )  f 
 q =  3 3   b (2.31)
3 1 1 1
f0 2 2 2 fc
Where the term f denotes the current, voltage or flux linkage. The inverse Park’s transfor-
mation is given by
fa cos θ − sin θ 1 fd
    
 f  2 − cos(θ − 2π ) − sin(θ − 2π ) 1 fq 
 b =  (2.32)
 
3 3
3

fc cos(θ + 2π3 ) − sin(θ + 2π
3 ) 1 f0
2.2.1 Stator flux linkages in dq0 components

Expressing the equations (2.22)-(2.24) in terms of dq0 components yields
Ψd = − Ld id + Laf d if d + La1d i1d (2.33)

Ψq = − Lq iq + La1q i1q (2.34)
Ψ0 =L0 i0 (2.35)
Where the newly defined inductances are are
3
Ld =Laa0 + Lab0 + Laa2 (2.36)
2
3
Lq =Laa0 + Lab0 − Laa2 (2.37)
2
L0 =Laa0 − Lab0 (2.38)
2.2.2 Rotor flux linkages in dq0 components

Similarly, expressing equations (2.28)-(2.30) in terms of dq0 components yields
3
Ψf d =Lf f d if d + Lf 1d i1d − Laf d id (2.39)
2
3
Ψ1d =Lf 1d if d + L11d i1d − La1d id (2.40)
2
3
Ψ1q =L11q i1q − La1q iq (2.41)
2
8
2.2.3 Stator voltage equations in dq0 components
By applying the Park’s transformation on the stator voltage equations (2.1)-(2.3) gives
ed =pΨd − Ψq pθ − Ra id (2.42)
eq =pΨq + Ψd pθ − Ra iq (2.43)
e0 =pΨ0 − Ra i0 (2.44)
2.3 Per unit representation

In many cases it is convenient to normalize the system variables in the form of per unit
quantities. The purpose is to express the equations (2.33)-(2.45) in per unit quantities,
which in chapter 3.2 is shown to be a better alternative than to use physical units.
A quantity in per unit is represented as the physical quantity divided with its base value.
The base values differ with regards to unit and whether they are stator or rotor quantities.
The derivation of the stator and rotor base quantities is provided in Appendix A.
2.3.1 Per unit stator voltage and flux linkage equations

Expressing equations (2.42)-(2.44) in per unit quantities, noting that es,base = is,base Zs,base =
ωbase Ψs,base gives

ed 1 Ψd Ψq
ωr Ra id
=p − − (2.45)
es,base ωbase Ψs,base Ψs,base ωbase Zs,base is,base

eq 1 Ψq Ψd ωr Ra iq
=p + − (2.46)
es,base ωbase Ψs,base Ψs,base ωbase Zs,base is,base

e0 1 Ψ0 Ra i0
=p − (2.47)
es,base ωbase Ψs,base Zs,base is,base
In per unit notation, equations (2.45)-(2.47) may be written as
1
ēd = pΨ̄d − Ψ̄q ω̄r − R̄a īd (2.48)
ωbase
1
ēq = pΨ̄q + Ψ̄d ω̄r − R̄a īq (2.49)
ωbase
1
ē0 = pΨ̄0 − R̄a ī0 (2.50)
ωbase
The stator flux linkage equations using the relationship Ψs,base = Ls,base is,base gives
Ψ̄d = − L̄d īd + L̄af d īf d + L̄a1d ī1d (2.51)

Ψ̄q = − L̄q īq + L̄a1q ī1q (2.52)
Ψ̄0 =L̄0 ī0 (2.53)
9
2.3.2 Per unit rotor voltage and flux linkage equations
Similarly the per unit rotor equations can be written as follows
1
ēf d = pΨ̄f d − R̄f d īf d (2.54)
ωbase
1
ē1d = pΨ̄1d − R̄1d ī1d = 0 (2.55)
ωbase
1
ē1q = pΨ̄1q − R̄1q ī1q = 0 (2.56)
ωbase
Ψ̄f d =L̄f f d īf d + L̄f 1d ī1d − L̄af d īd (2.57)

Ψ̄1d =L̄f 1d īf d + L̄11d ī1d − L̄a1d īd (2.58)
Ψ̄1q =L̄11q ī1q − L̄a1q īq (2.59)
2.3.3 Complete set of electrical equations in per unit

An alternative representation of the equations (2.48)-(2.59) are as follows
1
ed = pΨd − Ψq ωr − Ra id (2.60)
ωbase
1
eq = pΨq + Ψd ωr − Ra iq (2.61)
ωbase
1
e0 = pΨ0 − Ra i0 (2.62)
ωbase
Ψd = − (Lad + Ll )id + Lad if d + Lad i1d (2.63)

Ψq = − (Laq + Ll )iq + Laq i1q (2.64)
Ψ0 =L0 i0 (2.65)
1
ef d = pΨf d − Rf d if d (2.66)
ωbase
1
e1d = pΨ1d − R1d i1d = 0 (2.67)
ωbase
1
e1q = pΨ1q − R1q i1q = 0 (2.68)
ωbase
10
Ψf d =Lf f d if d + Lf 1d i1d − Lad id (2.69)
Ψ1d =Lf 1d if d + L11d i1d − Lad id (2.70)
3
Ψ1q =L11q i1q − Laq iq (2.71)
2
2.4 Derivation of machine parameters

So far the mathematical expressions for a synchronous machine and its parameters have
been presented. However the problem regarding the operational parameters is that they
cannot be directly determined from measured responses of the synchronous machine [14].
Instead the operational parameters are derived from observed behaviour of the synchronous
machine, in terms of time constants.
The first step is to express the relationship between the d- and q-axis flux linkages, currents
and voltages in the following incremental form
∆Ψd (s) =G(s)∆ef d (s) − Ld (s)∆id (s) (2.72)
∆Ψq (s) = − Lq (s)∆iq (s) (2.73)
Where G(s) is the stator to field transfer function, Ld (s) and Lq (s) are the d- and q-axis
operational inductances. The s is the Laplace operator and the prefix ∆ denotes the incre-
mental values. The d-axis transfer functions are calculated by expressing equations (2.63),
(2.66), (2.67), (2.69) and (2.70) in incremental form.
The first step is to express the relationship between the terminal flux linkage, current
and voltages as follows
Ψd (s) = − Ld id (s) + Lad if d (s) + Lad i1d (s) (2.74)

Ψf d (s) = − Lad id (s) + Lf f d if d (s) + Lf 1d i1d (s) (2.75)
Ψ1d (s) = − Lad id (s) + Lf 1d if d (s) + L11d i1d (s) (2.76)
The rotor voltage equations are
ef d (s) =sΨf d (s) − Ψf d (0) + Rf d if d (s) (2.77)

0 =sΨ1d (s) − Ψ1d (0) + R1d i1d (s) (2.78)
Expressing the d-axis rotor voltage equations in incremental form yields
∆ef d (s) =sΨf d (s) + Rf d if d (s)

(2.79)
= − sLad ∆id (s) + (Rf d + sLf f d )∆if d (s) + sLad ∆i1d (s)
11
0 =sΨ1d (s) + R1d i1d (s)
(2.80)
= − sLad ∆id (s) + sLad ∆if d (s)(R1d + sL11d )∆i1d (s)
The currents if d and i1d are eliminated by expressing them as a function of ef d and id
quantities. From equation (2.79) and (2.80) the following equations are obtained
(R1d + sL11d )∆ef d (s) + sLad (R1d + sL1d ∆id (s))

∆if d (s) = (2.81)
s (L11d Lf f d − L2ad ) + s(L11d Rf d + Lf f d R1d ) + R1d Rf d
2
−sLad ∆ef d (s) + sLad (Rf d + sLf d )∆id (s)

∆i1d (s) = (2.82)
s (L11d Lf f d − L2ad ) + s(L11d Rf d + Lf f d R1d ) + R1d Rf d
2
By substituting equation (2.77) and (2.78) in the incremental form of equation (2.74), it
is possible to express the d-axis quantities in the form of equation (2.72). The transfer
functions are
1 + (T4 + T5 )s + T4 T6 s2
Ld (s) = Ld (2.83)
1 + (T1 + T2 )s + T1 T3 s2
(1 + sT1d )
G(s) = G0 (2.84)
1 + (T1 + T2 )s + T1 T3 s2
Where
Lad L1d
G0 = T1d =
Rf d R1d
Lad + Lf d Lad + L1d
T1 = T2 =
Rf d R1d
(2.85)
1 Lad Lf d 1 Lad Ll
T3 = L1d + T4 = Lf d +
R1d Lad + Lf d Rf d Lad + Ll

1 Lad Ll 1 Lad Lf d Ll
T5 = L1d + T6 = L1d +
R1d Lad + Ll R1d Lad Ll + Lad Lf d + Lf d Ll
In factored form, equation (2.83) and (2.84) can be written as follows
(1 + sTd0 )(1 + sTd00 )

Ld (s) = Ld 0 )(1 + sT 00 ) (2.86)
(1 + sTd0 d0
(1 + sT1d )
G(s) = G0 0 )(1 + sT 00 ) (2.87)
(1 + sTd0 d0
Similarly with the q-axis, the operational inductance is
12
(1 + sTq0 )(1 + sTq00 )
Lq (s) = Lq 0 )(1 + sT 00 ) (2.88)
(1 + sTq0 q0
(2.89)
For the d-axis, the following relationships apply for the time constants
0 00
(1 + sTd0 )(1 + sTd0 ) = 1 + s(T1 + T2 ) + s2 (T1 T3 ) (2.90)
(1 + sTd0 )(1 + sTd00 ) 2
= 1 + s(T4 + T5 ) + s (T4 T6 ) (2.91)
2.4.1 Classical expression for time constants

If it often enough for a synchronous machine model to behave approximately as close to
reality as possible, for the sake of simplicity. One approximate method is to use the clas-
sical definition for the time constants, as explained in [14]. Considering that T2 and T3
are much smaller than T1 . T5 and T6 are much smaller than T4 . Thus equation (2.90) and
(2.91) can be simplified to
0 00
(1 + sTd0 )(1 + sTd0 ) ≈ (1 + sT1 )(1 + sT3 ) (2.92)
(1 + sTd0 )(1 + sTd00 ) ≈ (1 + sT4 )(1 + sT6 ) (2.93)
where
0 00
Td0 ≈ T1 Td0 ≈ T3 Td0 ≈ T4 Td00 ≈ T6 (2.94)
The q-axis open circuit time constant for a salient pole synchronous machine is given by
00 Laq + L1q
Tq0 = (2.95)
R1q
Note that there is no open circuit q-axis transient time constant for a salient pole machine.
2.4.2 Classical expression for transient and subtransient inductances

Considering equation (2.88), during steady state when s = 0, the d-axis operational induc-
tance is equal to the d-axis inductance. During a rapid transient, where s tends to infinity,
equation (2.86) is
Td0 Td00 00
Ld (∞) = Ld 0 T 00 = Ld (2.96)
Td0 d0
13
With no damper winding in the d-axis, equation (2.96) is
Td0 0
Ld (∞) = Ld 0 = Ld (2.97)
Td0
thus
Lad Lf d L1q
L00d =Ll + (2.98)
Lad Lf d + Lad L1q + Lf d L1q
Lad Lf d
L0d =Ll + (2.99)
Lad + Lf d
For a salient pole machine, the subtransient q-axis inductance is
Laq L1q
L00q = Ll + (2.100)
Laq + L1q
2.5 Power system components

The external system consists of transformers, transmission lines, substations and other
components within a power system. For this master thesis the only necessary compo-
nents to include is a step up delta-wye three phase transformer and the transmission line
components.
2.5.1 ∆−Y step up transformer

There are a wide variety of transformers ranging from (Y−Y), (∆−∆), (∆−Y) and (Y−∆).
The (∆−Y) type of transformer is the most commonly used step up transformer from the
low voltage side at the generator to the high voltage side at the first end of the transmission
line [15]. The reason for this is that a (∆−Y) connected transformer has the advantage
of reduced costs and more stable with respect to unbalanced loads [15]. The layout for a
three phase (∆−Y) transformer is as follows.
14
Figure 2.3: ∆ − Y transformer
Based on observations in Simulink, an unfortunate property of the three phase transformer

is that it adds a harmonic component equal to the net frequency to the field current. This
means that further harmonics are present during steady state operation, which will cause
further inaccuracies for the estimated rotor position. The per phase equivalent circuit for
a three phase transformer is shown in figure 2.4.
Figure 2.4: Per phase equivalent circuit of a transformer [15].
2.5.2 Transmission line model

The transmission line is modeled after the lumped parameter pi-line model, consisting of a
series resistance and inductance between two line to ground capacitances per transmission
line segment. The single segment three phase equivalent circuit for the pi-line model is
shown in figure 2.5.
15
Figure 2.5: Single segment pi-model of a three phase transmission line, where C denotes
the Clarke transform.
To simplify the equations describing the behavior of the transmission line, the following
Clarke’s transform is used
 √ 
1 2 0
1  −1
p 
T =√ 1 √ 2
3/2 (2.101)
3 −1
p 
1 √ 2
3/2
Where the simplified equations for the transmission line are expressed as follows
R + 2Rm L + 2M
   
 dI10
V10 − V20 =  R − Rm  0  L−M
 I1 + 

dt

R − Rm L−M
(2.102)
Cg
 
0 0 Cg + 3Cl  dV20
I1 − I2 =  (2.103)

dt

Cg + 3Cl
I10 =T I1 (2.104)
I20 =T I2 (2.105)
V10 =T V1 (2.106)
V10 =T V1 (2.107)
16
3 Rotor position estimation
3.1 Extended Electromotive Force based on the voltage concept
The method of estimating the rotor position for a salient pole synchronous generator, only
by considering stator and rotor measurement values is a very complicated task. Therefore
the methodology to achieve this is to start with simpler methods, and gradually taking ad-
vantage of more complicated methods to improve the rotor position estimation. The first
task is to use a conventional rotor position estimation, which later on will be evident that
it only works during steady state cases for a wound-rotor synchronous generator. Steady
state in this sense implies that there are no current transients in the field winding and the
d-axis damper winging, along with currents induced in the damper windings. The chosen
method is based on the voltage concept of deriving the rotor position, from the available
stator voltage and current measurement data. The method is used in [1], [2] to estimate
the rotor position for PM-based synchronous machines.
For simplicity the three phase notation is represented in a two phase (αβ) equivalent
reference frame. The relationship is expressed as
# fa
 
− 12 − 21
" # "
fα 2 1 √ √ f 
= 3 3  b (3.1)
fβ 3 0 − 2 − 2
fc
It should be noted that for this case the synchronous machine parameters and variables are
in physical quantities, not per unit quantities. The reason is because the following sources
[1], [2] use physical quantities to derive the rotor position. It is also more convenient since
only stator d- and q-axis parameters are needed in order to estimate the rotor position.
3.1.1 Dynamic model of a salient-pole synchronous machine

The dynamic model of an synchronous machine, relating the d- and q-axis voltages to d-
and q-axis stator currents rotor and damper winding currents is expressed as follows
" # " #" # " #
ed −Ra − pLd ωre Lq id Lad pif d + Lad pi1d − Laq ωre i1q
= + (3.2)
eq −ωre Ld −Ra − pLq iq Laq pi1q + Lad ωre if d + Lad ωre i1d
The transformation from d-q axis components to αβ components is

" # " #" #
fα cos(θre ) − sin(θre ) fd
= (3.3)
fβ sin(θre ) cos(θre ) fq
Applying equation (3.3) in equation (3.1) gives
17
" # " # " # " #
eα iα −L − ∆L cos(2θre ) −∆L sin(2θre ) iα
= Ra + ·p
eβ iβ −∆L sin(2θre ) −L + ∆L cos(2θre ) iβ
" #
cos(θre )
+ (Lad pif d + Lad pi1d − Laq ωre i1q ) (3.4)
sin(θre )
" #
− sin(θre )
+ (Laq pi1q + Lad ωre if d + Lad ωre i1d )
cos(θre )
where
Ld + Lq Ld − Lq
L= ∆L =
2 2
The problem with equation (3.4) is that both arguments θre and 2θre are present. It is
therefore difficult to solve the rotor position θre directly from equation (3.4). What is
necessary is to reconstruct equation (3.4) as suggested in [1], [2], where the first step is
express equation (3.2) in a voltage/flux model as follows
eα = − Ra iα − Ld p(id cos(θre ) + Lq p(iq sin(θre )) +vd (cos(θre )) + vq (sin(θre ))

| {z }
pλα
(3.5)
eβ = − Ra iβ − Ld p(id sin(θre ) − Lq p(iq cos(θre )) +vd (− sin(θre )) + vq (cos(θre ))
| {z }
pλβ
where
vd =Lad pif d + Lad pi1d − Laq ωre i1q

vq =Laq pi1q + Lad ωre if d + Lad ωre i1d
Equation (3.5) contains both voltage terms and derivatives of flux terms denoted as pλα
and pλβ . Rearranging equations (3.5) yields
" # " # " #" #
eα iα L 0 iα
= − Ra −p
eβ iβ 0 L iβ
" # " #! " # " #
cos(θre ) sin(θre ) id cos(θre ) − sin(θre ) (3.6)
− ∆Lp + vd + vq
sin(θre ) − cos(θre ) iq sin(θre ) cos(θre )
| {z }
Position related terms V (θre )
3.1.2 Model Reconstruction Based on Voltage Concept

The position related terms of the set of equations (3.6) contains both sine and cosine re-
lated terms. For the PMSM in [1], [2] the vd term is zero, because there is no field current
or damper windings. Thus it is possible to eliminate all cosine terms from the first row
18
and sine terms from the second row, by eliminating the d- and q-axis currents from the
position related terms V (θre ). This is done by using the dq - αβ transform (3.3). To ease
further calculations, the position related terms are rearranged as follows
" #
−∆Lp(id cos(θre + iq sin(θre )) + vd cos(θre ) − vq sin(θre )
V (θre ) = (3.7)
−∆Lp(id sin(θre − iq cos(θre )) + vd sin(θre ) + vq cos(θre )
By applying the dq - αβ transform to the currents, the following set of equations is obtained
" #
−∆Lp(iα + 2iq sin(θre )) + vd cos(θre ) − vq sin(θre )
V (θre ) = (3.8)
−∆Lp(iβ − 2iq cos(θre )) + vd sin(θre ) + vq cos(θre )
By applying (3.3) two more times, the following set of equations are obtained
" #
−∆Lp(iα ) − 2ωre ∆Liq cos(θre ) − (vq + 2∆Lp(iq )) sin(θre ) + vd cos(θre )
V (θre ) =
−∆Lp(iβ ) − 2ωre ∆Liq sin(θre ) + (vq + 2∆Lp(iq )) cos(θre ) + vd sin(θre )
" #
−∆Lp(iα ) − 2ωre ∆Liβ − (vq + 2ωre ∆Lid + 2∆Lp(iq )) sin(θre ) + vd cos(θre )
=
−∆Lp(iβ ) − 2ωre ∆Liα + (vq + 2ωre ∆Lid + 2∆Lp(iq )) cos(θre ) + vd sin(θre )
(3.9)
3.1.3 Reconstructed salient-pole synchronous machine model

With the new expression of the rotor position terms, the set of equations (3.6) can be
represented as
" # " #" # " #" # " #
eα −Ra −ωre (Ld − Lq ) iα Ld 0 iα − sin(θ re )
= −p + vq0
eβ ωre (Ld − Lq ) −Ra iβ 0 Ld iβ cos(θre )
" #
cos(θre )
+ vd
sin(θre )
(3.10)
where
vq0 = vq + 2ωre ∆Lid + 2∆Lp(iq )
Equations (3.10) can be expressed in a simpler form as

" 0# " # " #
eα
0 − sin(θre ) cos(θre )
= vq + vd (3.11)
e0β cos(θre ) sin(θre )
where
e0α
" # " # " #" # " #" #
eα −Ra −ωre (Ld − Lq ) iα Ld 0 iα
= − +p
e0β eβ ωre (Ld − Lq ) −Ra iβ 0 Ld iβ
19
It is evident that when vd is equal to zero, which is a valid assumption when the derivative
of the field and d-axis current is zero, along with the q-axis current being zero. The rotor
position can be then calculated as
θ̂re = − arctan 2(e0α , e0β ) (3.12)
Here the term arctan2 refers to the arctan function spanning from −π to π, compared to
the conventional arctan function that is limited to the first and the fourth quadrant (i.e
− 12 π to 12 π). The term θ̂re denotes the estimated rotor position.
For a wound-rotor synchronous machine, equation (3.12) is not valid when vd is non-zero.
By eliminating vq0 from equations (3.11) and solving for the rotor position, the following
equation is obtained
 
vd  − arctan 2(e0α , e0β )
θ̂re = arcsin  q
02
eα + eβ 02
  (3.13)
Lad pif d + Lad pi1d − Laq ωre i1q 
= arcsin  q − arctan 2(e0α , e0β )
02
eα + eβ02
Due to the arcsine term being limited to the first and fourth quadrant, the rotor position
cannot be directly solved from equation (3.13), even with the knowledge of field and damper
winding currents. This means that the equivalent EEMF rotor position estimation based on
the voltage concept is valid close to steady state, for a wound-rotor synchronous machine.
In Section 3.2, a method to compensate for the error induced by high d-axis field/damper
currents along with the q-axis damper winding current is presented.
3.2 Dynamic state rotor position estimator

The equivalent EEMF based method was shown to only function within the assumption
that the system is close to steady state. The dynamic state rotor position estimator
works by computing initial state dq-axis currents and voltages, when the field current
dynamics reaches a certain level. Then based of the initial values the further dq-axis
values are calculated numerically by solving a set of first order differential equations, until
the system is back in steady state. The idea is to use the equivalent EEMF based on
the voltage concept to determine initial values necessary for estimating dq-axis or rotor
reference frame values. The latest sample values are used as initial conditions for the
dynamic state rotor position estimator. The estimated d and q axis components contains
information about the rotor position. Due to the lack of information in the q-axis, it is
only possible to estimate the d-axis stator current, without the use of the rotor position
output as a feedback loop. Predicting the q-axis current from previous rotor position data
has the unfortunate attribute for making the estimator unstable. An another aspect is the
accumulation of large errors for each iteration. This will only worsen the results and is not
recommended.
20
3.2.1 Estimating d-axis current
Before explaining the algorithm for estimating id , a list and explanation of the simplifica-
tions and assumptions made for estimating id are explained first.
• The effects of saturation are neglected, thus Lad is assumed to be constant.
• The backwards Euler’s method is used for differentiation.
• The initial value for the d-axis damper winding current is assumed to be zero.
• The field voltage ef d is known, which for the static excitation system has not been
explained how it is derived.
The measured stator components include the stator currents and voltages. The measured
rotor components include the field current and voltage. The machine parameters are as-
sumed to be known as well and the effects of saturation are neglected. Saturation will only
affect Lad , hence if the effects of saturation were to be considered, an external method of
estimating the mutual inductance Lad for each sample would be necessary.
The flux linkage and flux linkage dynamic equations are formulated as follows for the
”Synchronous Machine Salient Pole (standard)” block in the ”SimPowerSystems compo-
nents” library in Simulink.
Ψd = −(Lad + Ll )id + Lad if d + Lad i1d (3.14)

Ψf d = −Lad id + Lf f d if d + L11d i1d (3.15)
Ψ1d = −Lad id + Lf 1d if d + L11d i1d (3.16)
1
Ψ̇f d = ef d − Rf d if d (3.17)
ωbase
1
Ψ̇1d = −R1d i1d (3.18)
ωbase
The indexes fd and 1d denote the field and d-axis damper winding values. All equations
are expressed in per unit quantities.
The algorithm is based on the equations (3.14)-(3.18) listed above. The procedure is ex-
plained as follows.
First, the initial values are determined, which corresponds to the last sample values ob-
tained from the steady state estimator (denoted as the sample number k), before the
dynamic state rotor position estimator is switched on. From the initial values, the initial
value for the field flux linkage is obtained from equation (3.15). It is assumed that the
damper winding current i1d [k] is close to zero. So equation (3.15) reduces to the following
equation.
21
Ψf d [k] = −Lad id [k] + Lf f d if d [k]
Now that all necessary initial conditions have been determined, the following method is
used to estimate id [k + i], where i stands for the i’th iteration in the algorithm. All deriva-
tives are calculated by using the backward Euler’s method.
1 Ψf d [k + i] − Ψf d [k + i − 1]
Ψ̇f d [k + i] = ef d [k + i] − Rf d if d [k + i] = =>
ωbase Ts · ωbase
Ψf d [k + i] = (Ts · ωbase ) · (ef d [k + i] − Rf d if d [k + i]) + Ψf d [k + i − 1]
The currents id [k + i] and i1d [k + i] can now be obtained from equation (3.15), (3.18) and
by differentiating (3.16).
Ψf d = − Lad id + Lf f d if d + L11d i1d

1
Ψ̇1d = − R1d i1d
ωbase
1 id [k + i] − id [k + i − 1] if d [k + i] − if d [k + i − 1]
Ψ̇1d = − Lad + Lf 1d
ωbase Ts ωbase Ts ωbase
i1d [k + i] − i1d [k + i − 1]
+ L11d
Ts ωbase
From the above equation, id [k + i] is solved by the following equation.
Lf f d if d [k + i] + Lad i1d [k + i] − Ψf d [k + i]
id [k + i] = (3.19)
Lad
Where i1d [k + i] is obtained from the following equation.
−Lf f d if d [k + i] + Ψf d [k + i] − L11d i1d [k + i − 1]

i1d [k + i] =
Lad − L11d − R1d Ts ωbase
(3.20)
Lad id [k + i − 1] + Lad (if d [k + i] − if d [k + i − 1])
+
Lad − L11d − R1d Ts ωbase
22
3.2.2 Rotor position estimation
The idea for estimating the rotor position, based on the knowledge of d-axis current is to
express id as a function of measured abc current phase components and the rotor posi-
tion. Expressing the abc components as their two-phase equivalent alpha-beta components
yields the following equation.
id [k + i] = iα [k + i] cos(θre [k + i]) + iβ [k + i] sin(θre [k + i]) (3.21)
Representing the right hand side of the above equation as a single sinusoidal term gives.
q
id [k + i] = i2α [k + i] + i2β [k + i] sin(θre [k + i] + ϕ[k + i]) (3.22)
Where ϕ[k + i] = arctan2(ialpha [k + i], ibeta [k + i]). Solving the rotor position from equation
(3.22) yields the following solutions.
 
id [k + i]  − arctan 2(ialpha [k + i], ibeta [k + i]) (3.23)
θre1 [k + i] = arcsin  q
i2α [k + i] + i2β [k + i]
 
id [k + i]
θre2 [k + i] = π − arcsin  q  − arctan 2(ialpha [k + i], ibeta [k + i])
2 2
iα [k + i] + iβ [k + i]
(3.24)
Here the problem is to determine which solution is the correct one. A method to analyze
which solution is the correct one is to check the behavior of the rotor speed and its deriva-
tive for each solution. A reasonable range for the rotor speed is around its rated speed, for
this case approximately 60 Hz, or 377 rad/s. The derivative should be close to zero. The
following equation is used to determine which solution to be the correct one.

D = | ωre1 [k + i] − 377 | + | ω̇re1 [k + i] | − | ωre2 [k + i] − 377 |
(3.25)
+ | ω̇re2 [k + i] | +
In equation (3.25), if D < 0, then θre1 [k + i] is chosen as the correct solution, respectively
for D > 0, then θre2 [k + i]. The term is used as a weighting term, since the behavior of
the rotor speed for both cases will be almost identical during steady state.
Equation (3.25) will not be useful during cases when the q-axis stator current is close
to zero. The reason is that the estimated d-axis current is not one hundred percent ac-
curate. The reason for the large errors might be due to the high relative error when id is
close to zero. The results give rise to pulse like errors within short time intervals. It is
therefore possible to use an integrating part in the estimator to eliminate them.
The integrator setting functions in parallel to the dynamic state rotor position estima-
tor. There are different methods to make sure the rotor position estimator switches to the
23
integrating setting, when the dynamic state rotor position estimator detects an abnormal
change or rotor position level, similar to that of equation (3.25). An another method is
to predict the current rotor position from previous rotor position estimations. A predictor
term ωbase · Ts is added to the previous rotor position θre [k + i − 1]. The q-axis stator
current can then be estimated, where any boundary close to iq = 0 switches the estimator
from the dynamic state to the integrating setting. A simplified equation for the integrating
part is as follows.
Z
θre [k + i] = θre [k 0 ] + ωbase dt (3.26)
Where sample k’ is the sample value when the integrating part is switched on. The integral
part of (3.26) resets to zero when the integrating setting is switched off.
At the moment, a mix of the first and the second method is used, as shown in equa-
tion (3.27) below.
c
D0 =| ωre [k + i − 1] − ωbase | + | ωre [k + i − 2] − ωbase | + (3.27)
1+ | îq [k + i] | /b
Where îq [k +i] is the predicted q-axis current. For D0 > a, the dynamic state rotor position
estimator is switched to the integrating part. The constant parameters in equation (3.27)
can be adjusted in order to improve the accuracy for the estimated rotor position.
24
4 Designing the simulation environment and the estimator
4.1 Simulation
The first step in the designing process is to chose an appropriate testing or simulation
environment, in order to test the functionality of the rotor position estimator. The re-
quirement for the estimator is to test it for a salient pole synchronous generator connected
to a large power system. Such a system will be difficult and time consuming to design
in a software program, depending on the size of the power system. However the entire
power system can be represented as an infinite bus, where the generator is connected to a
step up transformer, which is connected through one or multiple transmission lines to the
infinite bus. For this case only two parallel transmission lines were used, where the second
transmission line will be subjected to three different types of line faults.
The chosen test environment was Simulink, where the SimPowerSystems application was
used. The single-machine-infinite-bus (SMIB) system as designed in Simulink is as follows
Figure 4.1: SMIB system in Simulink.
As illustrated in figure 4.1, the power system configuration consists of the power system
blocks, mechanical blocks and the solver configuration block. There is also a fault block
which generates a particular type of fault for a certain amount of time. This will be used
to verify the functionality of the rotor position estimator. During simulation the output
values are connected to a multiplexer port (the black box), where the values are saved in
a matrix containing each sample values corresponding to a specific time instance. This
matrix is then saved into a separate file. The convenience for this is to make the estimator
simulation run faster. Because for each simulation Simulink has to recalculate all the
25
output parameters, which is unnecessary considering that the estimator runs parallel to
the power system simulation. This is not much of an importance for obtaining the final
results, but reducing the simulation time helps to run more extensive simulations. Due to
the limited computational performance the simulations were limited from 0 to 15 seconds,
which implies from start that the synchronous generator has enough time to reach steady
state, being exposed to a fault and reaching steady state again.
4.1.1 Power system
Figure 4.2: Power system in Simulink.
The external part or the power system is is simulated according to figure 4.2. The syn-
chronous machine port ”∼” is connected to the low voltage side ”∼2” port of the step
up three phase transformer. The step up transformer block is an ideal ∆−Y transformer,
where the effects of core saturation are neglected. The high voltage side ”∼1” of the trans-
former is connected to the transmission network consisting of two 1 km long parallel lines.
The point of fault which is located in the middle the second line is modeled as two 0.5 km
long transmission lines with the a time based fault generator connected in between. The
circuit breakers are set to operate at t = 13.05 s, hence the duration of fault is 0.05 s or
50 ms. In reality both lines would have circuit breakers installed, since the tests are run
by having the fault occur at line 2, it is unnecessary to add the additional circuit breakers
at both ends of line 1. The infinite bus is represented as a constant AC voltage source.
4.1.2 Synchronous generator

The input systems for the synchronous generator consists of the excitation system and the
mechanical system with a constant input torque. To make the system more simple, the
input signals consists of a constant torque and field voltage. Since the estimator should
work regardless of the input signals, it is not necessary to include any further control
systems if the generator still manages to maintain its synchronism.
26
Figure 4.3: Synchronous generator together with the excitation system and constant torque
input.
The input to the excitation system is a constant DC voltage, which corresponds to the field
voltage. Due to time constraint a static excitation system were not considered, instead the
simple constant DC field voltage input is considered only. The output is the field current
If dpu . The electrical conserving ports fd+ and fd− act as signals between the field winding
and other parts of the machine, to determine the output field current.
The lines and blocks marked in green is the mechanical system, which consists of an ideal
torque source, machine inertia and a mechanical reference. The torque source act as an
ideal source of mechanical energy to generate a torque proportional to the input torque
signal, which corresponds to the step torque signal. The input torque signal is the physical
torque value and not its per unit value.
4.1.3 Output vectors from the synchronous generator

The output vectors should represent the three phase stator voltages and currents, along
with the field current. Since the SimPowerSystems synchronous machine block only out-
puts dq0 quantities, the actual rotor position is needed in order to obtain the abc or α-β
components. The rotor position is not an output vector and thus it has to be obtained
from the rotor speed. The conversion from dq0 to α-β components is done in the same
Simulink file where the estimator is located at. This choice was done arbitrarily to ex-
clude further configurations for the simulation. As seen in figure 4.4 the chosen output
vectors are converted from per unit quantities to their physical quantities. This was also
an arbitrary choice due to the EEMF estimator using physical quantities.
27
Figure 4.4: The chosen output vectors from the synchronous machine block.
4.2 Estimator
The EEMF and the dynamic state rotor position estimator (as explained in section 3) are
the most significant parts of the entire system. But they rely on the output values from
the simulated SMIB system, which are the inputs of the estimator. The input values are
sampled with a fixed sampling frequency fS , and then connected to the EEMF and the
dynamic state rotor position estimator. This system is shown in figure 4.5 (input values
connected to the estimator) and inside the estimator in figure 4.6.
Figure 4.5: Input values feed into the estimator.
28
Figure 4.6: The estimator function blocks, where the dq-αβ transform is marked with red,
the sampling blocks in blue, the estimator blocks in green and computing the actual rotor
position in magenta.
Considering figure 4.5, the input values are loaded from a file and connected to a demulti-
plexer block, which splits the different input values to be feed into the estimator. In figure
4.6, the d- and q-axis currents and voltages are converted into α-β components before
sampling (marked with a red square) with the help of the actual rotor speed, which rep-
resent the two-phase equivalent values which can be measured. As mentioned previously,
the Simulink application SimPowerSystems only outputs dq components.
The actual rotor speed is the integral of the output rotor speed. To limit the rotor position
within the range [−π, π] the rotor speed after it is integrated is separated into its sine
and cosine values, where the rotor position is then calculated with the help of the arctan2
function. This will limit the rotor position within this interval and is necessary in order to
calculate the estimated rotor position error.
All input values are sampled, except the actual rotor speed where the actual rotor po-
sition is obtained by the integrator. The rotor position is used for obtaining the rotor
position error by subtracting the actual rotor position from the estimated. The EEMF and
the dynamic estimators are marked with green.
4.2.1 Designing the EEMF subsystem

The EEMF system is based on equation (3.11) for solving the EEMF terms.
29
Figure 4.7: EEMF based rotor position estimator.
The current derivative is solved by using the backward Euler’s method. As shown in figure
4.7 the [k − 1] sample is obtained from the delay block. The current as sample k − 1
is subtracted from the current at sample k and then divided by the sampling time Ts .
However the multiplication of Ld and the division of Ts are done simultaneously at the
gain blocks Ld /Ts . The rotor angular speed is represented as a constant input value of the
base angular speed, which is a valid assumption in steady state. The final step consists of
solving equation (3.12) from the equivalent EEMF terms.
4.2.2 Designing the dynamic state rotor position estimator

The dynamic state rotor position estimator is the most complex part, hence it will be
necessary to explain the design in different parts, rather than explaining the function as
a whole. The first part is where all the initial values are defined. The initial values are
defined when the difference of the field current at sample k and k−1 reaches a certain value.
To compensate for the loss of accuracy during this initial dynamic state, the input values
are calculated by using a previous rotor position value and then predicted by assuming the
rotor speed maintains its constant steady state value during the beginning of the dynamic
30
state.
Figure 4.8: Initial dynamic state observer (green), Initial state estimation (blue) and the
estimated rotor speed from the EEMF (red).
In the blue region in figure 4.8, the initial d-axis current is estimated with the help of
equation (3.15), assuming the d-axis damper winding current is zero. The initial dynamic
state is detected by the configuration shown in the green region. Since the difference
between the field winding current at the current and previous sample is not exactly zero at
steady state, thus the dynamic state boundary for the current difference is chosen arbitrarily
at a certain level. As seen in figure 4.8, the dynamic state rotor position estimator is
activated for current difference values higher than 0.0001 pu. During this case the switch
connects the output to the uppermost input and thus integrates a positive value. The
integrator is floored at zero, so that the output is limited to a digital output signal of 0
and 1. Any positive value at the output of the integrator will generate an output signal
as 1. When the current difference is below the boundary value, the switch will connect
the output to the second port, which is connected to a constant negative value. The
output value from the integrator will decrease until it reaches zero, thus the dynamic state
rotor position estimator will be switched off after a certain amount of time. The current
difference will reach values under the boundary limit for short periods of time, therefore it
is necessary to have a lower constant negative input value to make sure the output value
doesn’t drop to zero during a period of dynamic operation.
31
Figure 4.9: Estimating the d-axis current.
The d-axis current is estimated with the help of equation (3.19) and (3.20), which cor-
responds to the subsystem blocks as seen in figure 4.9. The estimated d-axis current is
converted into its physical non per unit value in the gain block, within the marked region
it is possible to compare both estimated and actual d-axis current to check the accuracy
of the d-axis estimator. To avoid errors accumulating in the feedback loop within the d-
axis damper winding current solver, the feedback loop starts to operate when the dynamic
state rotor position estimator is switched on. In other cases the previous sample value or
i1d [k − 1] is set to zero.
32
Figure 4.10: Detecting the correct solution from the two calculated rotor positions from
the d-axis current.
The red marked region in figure 4.10 is the entire system of detecting the right solution
between equation (3.23) and (3.24). The green arrow points out where the output value
of equation (3.23) is located. Similarly the blue arrow points out where the output of
equation (3.24) is. The system which solves the inequality (3.25) is in the magenta region.
As mentioned previously in section 3.2.2 the estimated rotor position will experience inac-
curacies when the q-axis current is close to zero. The integrating part which functions by
integrating from the initial rotor position with an angular rotor speed assumed to be the
rated angular speed, when the q-axis current gets close to zero, along with a high derivative
of the rotor position. Figure 4.11 shows the input ports to the integrating part. Figure
4.12 is the inside of the subsystem of the integrating part.
33
Figure 4.11: Input values for the integrating part.
Figure 4.12: Integrating part.
The output from the corrected rotor position estimation is feed into a block that switches
the dynamic state rotor position and the integrator. The other is feed into an speed observer
(marked in red) which purpose is to detect high or low rotor speeds below rated speed.
The green region includes the predicted q-axis current and the inequality (3.27) within
the magenta colored region. Since the estimated rotor position at sample k is unknown, a
predicted rotor position at sample k is used instead. The integrator is shown in the blue
region and consists of an initial rotor position detector and the integrator itself. To avoid
34
large initial value errors, the initial rotor position is predicted from previous sample values.
For this case from the previous 5 samples or k − 5 samples, which gives a predictor term of
5 · ωbase . When the inequality is less than ”a” the predicted rotor position is continuously
updated in the sample and hold block. At the moment when the inequality is higher than
”a” the switching signal for the sample and hold block is zero, hence the rotor position
maintains a constant value which is equal to the last sample value when the input signal
was one. This serves as the constant term θre [k 0 ] in equation (3.26). When the integrating
part is switched off, the input for the integrator is switched to a constant negative value,
high enough to drive the output value quickly to zero (where the minimum output value
from the integrator is set to zero).
4.2.3 Calculating the rotor position estimation error
Figure 4.13: Calculating the rotor position error.
To determine if the estimated rotor position is accurate it has to be compared to the actual
rotor position. Since the estimated rotor position is a discrete quantity, for convenience
the actual rotor position were sampled as well, in order to reduce the in between sample
errors. If necessary the estimated rotor position could be reconstructed to generate a non
discrete output value. In figure 4.13 the subsystem in the green area is where the actual
rotor position is sampled, the red area is the EEMF based rotor position estimation, the
blue is the dynamic state rotor position estimator. To compare the performance between
the different steps of how to estimate the rotor position, from top to bottom in the right of
figure 4.13, only the EEMF estimator is used, then with the included dynamic state rotor
position estimator and the final step of estimating the rotor position error.
35
5 Results and discussion
5.1 Input data
The tests were performed with a hydro power generator, where the system parameters are
shown in table 1. The parameters were obtained from [16] for the H6 type of hydro power
generator.
Table 1: Machine parameters
Parameter Parameter description Value

Ld d-axis inductance 1.000 pu
Lq q-axis inductance 0.620 pu
Ll Leakage inductance 0.170 pu
L0d d-axis transient inductance 0.260 pu
L00d d-axis subtransient inductance 0.235 pu
L00q q-axis subtransient inductance 0.264 pu
Ra Stator resistance 0.004 pu
0
τd0 Open-circuit d-axis transient time constant 7.1 s
τd00 Short-circuit d-axis subtransient time constant 0.035 s
τq00 Short-circuit q-axis subtransient time constant 0.035 s
p Number of pole pairs 10
The input parameters for the step-up transformer and the transmission lines are shown in
table 2 respectively table 3. Due to the difficulty of finding reliable transformer and trans-
mission line data, the default values from the step-up transformer and the transmission line
blocks in SimPowerSystem were used. The default values for the line-ground capacitance
and the mutual resistance are both zero, while the number of segments is equal to one.
This suggests that the default values are based of a simplified model.
Table 2: Transformer parameters

V2 /V1 Turns ratio 200 kV / 13.8 kV
Rlp Primary resistance 0.01 pu
Xlp Primary leakage reactance 0.001 pu
Rls Secondary resistance 0.01 pu
Xls Secondary leakage reactance 0.001 pu
Rsh Shunt magnetizing resistance 500 pu
Xsh Shunt magnetizing reactance 500 pu
36
Table 3: Line parameters

R Line resistance 0.02 Ω/km
L Line inductance 0.5 mH/km
M Mutual inductance 0.1 mH/km
Cl Line-line capacitance 0.3 µF/km
Cg Line-ground capacitance 0 µF/km
Rm Mutual resistance 0 Ω/km
Rp Parasitic series resistance 1 µΩ/km
Gp Parasitic parallel conductance 1 µS/km
Ns Number of segments 1
The estimator input data are the same as the values presented in table 1, except the other
circuit parameters have been determined by the time constants, transient and subtransient
inductances. The machine parameters derived from the given time constants, transient and
subtransient inductances are refered to appendix A.2. As for the constants for equation
(3.25) and the integrating part in equation (3.27), the constant values were chosen as follows
= 10, a = 6, b = 60, c = 10
Other input data included the initial conditions for the synchronous generator, step up
transformer and the infinite bus. The choice of initial can be chosen so that the syn-
chronous generator starts from a steady state condition. As long as the system converges
to a steady state relatively fast, it is not necessary to make the work more complicated in
order to stabilize the system at start. The initial values were chosen by trial an error to
make the system stable enough for it to converge fast. The initial conditions are shown in
table 4.
Table 4: Initial conditions
Type Parameter description Value

Generator Terminal voltage magnitude 13.80 kV
Generator Terminal voltage angle 0 deg
Generator Terminal active power 0.5 pu
Generator Terminal reactive power 0 pu
Transformer Initial primary currents [53.8558, −53.8558, 0] A
Transformer Initial secondary currents [107.7116, −53.8558, −53.8558] A
Transformer Initial magnetizing currents [0, 0, 0] A
Transformer Initial fluxes [0, 0, 0] Wb
5.2 Simulation results

The simulations were carried out for a maximum of 15 seconds. During start, the generator
starts from a non steady state operation, over time it will settle down to steady state
operation. Three different cases of faults consisting of line-to-ground, phase-to-phase and
37
three phase fault were considered when running the simulation. The faults were generated
at t = 13 s where after 50 ms from the beginning of the fault line 2 is disconnected by the
circuit breakers. The purpose is to check how the rotor position estimator will react for
different types of faults and how accurate the estimated rotor position is.
5.2.1 Single-phase-to-ground fault

The phase-to-ground fault occurs from phase a to ground as shown in figure 5.1.
Figure 5.1: Illustration of a single-phase-to-ground (P-G) fault.
From the EEMF based method, the resulting rotor position error is shown in figure 5.2.
During start the stator and rotor quantities are not at their steady state levels, hence
why the estimation error is very high. Subsequently when the system has converged to
steady state, the error maintains a very low level, which is shown more in detail in figure
5.4 with the included dynamic state rotor position estimator before fault. As observed
at t = 13 s, the fault will cause the system to go into non steady state, hence why the
error during the fault is very high, which after clearing the fault quickly reaches steady
state. The dynamic state rotor position estimator is set to operate at t > 12 s, in order to
avoid it to operate at start. The resulting rotor position estimation error with the included
dynamic state rotor position estimator is shown in figure 5.3. For t = 13 s when the fault
occurs there is a very small change in the rotor position error. A closeup between the
fault intervall figure 5.4 shows clearly how the error behaves during pre-fault, during fault
and post fault. During pre fault at steady state operation, the estimated rotor position
error oscillates with an amplitude of approximately 0.001 rad/s, observed between 12 and 13
seconds. This is observed to be caused by the step up transformer, as mentioned previously.
During fault there are positive and negative error spikes present, these correspond to the
inaccuracies caused by the q-axis passing the zero boundary. At the peak values of the
spikes the integrator is switched on, otherwise they would continue up to very high values.
During post fault (after the fault is cleared) the system is stabilizing and the estimated
rotor position error maintains a steady value. When the dynamics of the field current has
maintained low values, the dynamic state rotor position estimator is switched off, as seen
around t = 13.65 s.
38
Figure 5.2: Rotor position estimation error for a P-G fault without the dynamic state rotor
position estimator.
Figure 5.3: Rotor position estimation error for a P-G fault without the dynamic state rotor
position estimator.
39
Figure 5.4: Rotor position estimation error for a P-G fault with the dynamic state rotor
position estimator.
Figure 5.5: Rotor position estimation error for a P-G fault with the dynamic state rotor
position estimator (close to fault).
5.2.2 Phase-to-phase fault

For this case the phase-to-phase (or line-to-line) fault occurs from phase a to b as shown
in figure 5.4.
40
Figure 5.6: Illustration of a phase-to-phase (P-P) fault.
Considering the case for the single-phase-to-ground fault, during fault the shape of the
estimation error as function of time displayed an high frequency oscillation, combined with
spikes or pulses. The same form can be observed in the phase-to-phase fault as shown in
figure 5.8. The difference is the larger spikes during a fault, which suggests that the q-axis
current not only passes through the zero boundary for the peak value to reach a relatively
high negative value, but the peak value to reach relatively close to the zero boundary. The
d- and q-axis currents as a function of time is referred to the appendix.
During post fault the offset error is larger than that of the single-phase-ground fault.
The reason might be due to a larger error when estimating the initial value for the d-
axis damper winding current. The initial d-axis current should not make that much of a
difference since it is based on a predicted value from a previous sample data.
Figure 5.7: Rotor position estimation error for a P-P fault without the dynamic state rotor
position estimator.
41
Figure 5.8: Rotor position estimation error for a P-P fault without the dynamic state rotor
position estimator.
Figure 5.9: Rotor position estimation error for a P-P fault with the dynamic state rotor
position estimator.
42
Figure 5.10: Rotor position estimation error for a P-P fault with the dynamic state rotor
5.2.3 Three-phase-fault
A three phase fault as illustrated in figure 5.7 is known as a symmetrical fault, due to the
simultaneous short circuit across all three phases [15].
Figure 5.11: Illustration of a three phase (3φ) fault.
The three phase fault for this particular case does not display very large error spikes or
pulses during a fault, as compared to the two previous cases. The offset error during post
fault is not as large compared to the line-to-line fault case. To make sure if this is true for
most three phase faults, other results obtained from different power system layouts and
fault locations should be included. Since the purpose is to test the performance of the rotor
position estimator, as a start it is enough to check its behaviour during different types of
faults.
43
Figure 5.12: Rotor position estimation error for a 3φ fault without the dynamic state rotor
position estimator.
Figure 5.13: Rotor position estimation error for a 3φ fault without the dynamic state rotor
position estimator.
44
Figure 5.14: Rotor position estimation error for a 3φ fault with the dynamic state rotor
position estimator.
Figure 5.15: Rotor position estimation error for a 3φ fault with the dynamic state rotor
5.3 Discussion
5.3.1 Choice of power system model
The choice of simulation environment for the power system was done to make it behave
as realistically as possible. Such a system would be incredibly complicated and futile if
accurate results can be achieved from a much simpler system. The simulation environment
was designed in Simulink with the help of the built in package ”simpowersystems”. As
compared to a real system, the salient pole synchronous generator was modeled on the
classical system and a simple field saturation model. Hence, if further tests were to be
performed considering the effects of saturation, an analysis should be made to determine
45
if a more detailed saturation model will be necessary in order to obtain more realistic
results. The type of step up transformer commonly used between a generator and the
transmission network is of the ∆−Y type, which was the most important consideration
when choosing the right transformer model. Although the transformer block from the
SimPowerSystems was an ideal transformer model in the sense that the effects of iron core
saturation were neglected. However the model did consider non-ideal properties such as
leakage inductance, eddy currents etc. The notion that a large power system network can
be represented as a constant ac voltage source, hence an infinite bus saves a lot of time
and reduces computational effort to run such a system.
5.3.2 Choice of estimator model

The most difficult part of this thesis was the choice of the right type of estimator design,
best adapted to a salient pole synchronous generator. A lot of studies have been made to
estimate the rotor position for a permanent magnet synchronous machine, with no presence
of damper windings. This is the reason why the choice of estimation model was difficult
for a more complicated synchronous generator. A lot of effort went to directly estimating
the rotor position from the known machine parameters and the measured rotor and stator
currents and voltages. The final choice was to split the estimator into two parts, the first
one considering steady state operation and the second one using the initial values from
the first estimator to solve a set of differential equations to calculate the d-axis current,
which contains the information of the rotor position. This was not enough considering that
by calculating the d-axis current by itself would yield two different solutions to the rotor
position. The method was to analyze the behavior of the rotor speed, for both solutions
and chose the one that gave the most reasonable behavior for a synchronous generator,
regarding its rate of acceleration and angular speed. For a generator with high inertia,
rapid accelerations should not be expected. For a machine connected to the power system,
a speed close to the rated speed is expected, unless during a case when the synchronous
machine loses synchronism, then higher rotor angular speeds are expected.
46
6 Conclusion and future work
6.1 Conclusion
The results from the three cases of line faults shows that the rotor position is improved
tremendously with the added dynamic state rotor position estimator for all three different
types of faults, with the effects of saturation disabled. The EEMF based rotor position
estimator worked well during steady state operation, however the intended purpose from
the beginning was to estimate the rotor position for both steady-state and dynamic cases,
which turned out to be impossible due to the present d-axis field and damper winding
current dynamics, along with the q-axis damper winding current. A much simpler method
for estimating the rotor position could have been utilized, but due to the versatility of the
EEMF method it was used instead, presumably since it is independent of the d and q axis
current dynamics. The dynamic state rotor position estimator could have been made even
more accurate by adjusting the integrating setting constants to more optimal values. As
seen for the phase to phase fault there are still instances when the initial rotor position
during the integrating phase has a relatively high error. An another method to improve the
estimation of the initial values for the dynamic state rotor position estimator is to make it
operate before a fault occurs. An another possible method is to start estimating the d-axis
current from previous sample data to the current sample data, within the sampling time
interval from the time the dynamic state rotor position estimator starts to operate. Due
to the time constraint, any further improved estimator designs/models were not consid-
ered. More details on how to make the rotor position estimation more suited for practical
applications will be discussed in Future work.
6.2 Future work

As an continuation with the current estimator design, the most important addition is to
include a method for estimating the d-axis mutual inductance to account for the effects of
saturation. If it is necessary, other methods could be chosen to estimate the rotor position
which are more suitable for taking the effects of saturation into account. Further it could
be useful to only limit the use for a sensor-less rotor position estimator for only a certain
cases. As example the rotor position could be estimated during normal operation only,
hence assuming that the rotor position should only operate when no faults occurs close to
the synchronous generator. What this implies is that the rotor position should not only be
estimated during steady state operation, but during power swings from other generators
or machines close to the generator from which the rotor position is estimated.
There is also a possibility for estimating the load angle with the swing center voltage
method. The load angle could be proven useful for estimating the rotor position.
47
Appendix A
A.1 Per Unit System for Stator and Rotor Quantities
es,base =11267.7 V
is,base =2070.82 A
fbase =60 Hz
ωbase =2πfbase
es,base
Zs,base =
is,base
Zs,base
Ls,base =
ωbase
es,base
Ψs,base =Ls,base is,base =
ωbase
if d,base =1300 A
ef d,base =13.5926 V
A.2 Derived Standard parameters from time constants

The following parameters are derived from (2.87) with the classical expression
Lf d =0.1009 pu
L1d =0.2342 pu
Rf d =3.4779 · 10−4 pu
R1d =0.0222 pu
L11d =Lad + L1d pu
48
Appendix B
B.1 Complementary results for P-G fault
Figure B.1: Measured field current if d in per-unit.
Figure B.2: Estimation error for the d-axis current [A] rms.
49
Figure B.3: Predicted q-axis current [A] rms.
B.2 Complementary results for P-P fault
50
B.3 Complementary results for 3φ fault
51
52
53
Appendix C
C.1 Calculating id and i1d in Simulink
Figure C.1: Calculating the d-axis damper winding current i1d .
Figure C.2: Calculating the d-axis stator current id .
C.2 Other Simulink systems
54
Figure C.3: dq to α − β transform.
Figure C.4: α − β to dq transform.
55
Figure C.5: Sampling.
56
Appendix D
D.1 Magnetic saturation
The model used in [14] to account for the effects of saturation is based on the following
assumptions.
1. Leakage inductances are independent of saturation
2. Leakage fluxes do not contribute to the iron saturation.
3. The saturation relationship under loaded conditions is the same as under no-load
conditions. Hence the open-circuit saturation curve can be used to simulate the
effects of saturation.
4. There is no magnetic cross coupling between the d- and q-axes.
5. Saturation effects for the quadrature axes are neglected, due to a salient pole rotor
having a larger air gap at the q-axis.
From the above assumptions, only the d axis inductance is affected by saturation. Hence
the effects of saturation is modeled as
Lad = Ksd Ladu (D.1)
where Ksd is the saturation factor and Ladu is the unsaturated value of Ladu . The saturation
factor is a function of the d- and q-axis flux linkage Ψat , where the saturation factor curve
is obtained from the open-circuit characteristics (assumption 3). Under no load conditions,
the open circuit characteristics curve is obtained by expressing the air-gap voltage Ea as a
function of field current if d , in per unit the air gap voltage is equal to the flux linkage Ψat .
During OCC the air gap voltage is equal to the terminal voltage, which makes it possible
to express Lad as a function of the flux linkage.
57
Figure D.1: Open-circuit characteristic saturation curve and the non-saturation curve [14].
The inductance Lad obtained by differentiating the OCC curve. The air-gap voltage is
expressed as
Ea = Et + (Ra + jXl )It (D.2)
During no load condition, the following conditions apply
id = iq = Ψq = ed = 0
Et = eq = Ψd = Lad if d
Therefore during no load condition the air-gap voltage equals the terminal voltage (equation
(D.1)).
D.2 Measuring the field current for a static excitation system

This part explains the method for measuring the rotor DC current, by the use of current
transformers. The idea is based on a patent owned by the company ABB [17].
58
Figure D.2: Structure of a static excitation system, together with a current measuring
device [17].
The idea is illustrated in figure D.2, where a three phase signal port is connected to an
external IED (intelligent electronic device) is connected between the secondary side of the
rotor excitation transformer and the thyristor bridge. Alternatively between the stator and
the primary side of the excitation transformer.
The dc side current denoted as IDC is calculated from the following equation
| Ia | + | Ib | + | Ic |
IDC = (D.3)
2
59
References
[1] Y. Zhao, Z. Zhang, W. Qiao, et al. “An Extended Flux Model-Based Rotor Position
Estimator for Sensorless Control of Salient-Pole Permanent-Magnet Synchronous Ma-
chines”. In: IEEE Transactions on Power Electronics 30.8 (Aug. 2015), pp. 4412–
4422. issn: 0885-8993. doi: 10.1109/TPEL.2014.2358621.
[2] Zhiqian Chen, M. Tomita, S. Doki, et al. “An extended electromotive force model
for sensorless control of interior permanent-magnet synchronous motors”. In: IEEE
Transactions on Industrial Electronics 50.2 (Apr. 2003), pp. 288–295. issn: 0278-
0046. doi: 10.1109/TIE.2003.809391.
[3] Y. Lee, Y. C. Kwon, and S. K. Sul. “Comparison of rotor position estimation per-
formance in fundamental-model-based sensorless control of PMSM”. In: 2015 IEEE
Energy Conversion Congress and Exposition (ECCE). Sept. 2015, pp. 5624–5633.
doi: 10.1109/ECCE.2015.7310451.
[4] Y. Zhao, Z. Zhang, W. Qiao, et al. “An extended flux model-based rotor position
estimator for sensorless control of interior permanent magnet synchronous machines”.
In: 2013 IEEE Energy Conversion Congress and Exposition. Sept. 2013, pp. 3799–
3806. doi: 10.1109/ECCE.2013.6647204.
[5] J. Liu, T. A. Nondahl, P. B. Schmidt, et al. “Rotor Position Estimation for Syn-
chronous Machines Based on Equivalent EMF”. In: IEEE Transactions on Industry
Applications 47.3 (May 2011), pp. 1310–1318. issn: 0093-9994. doi: 10.1109/TIA.
2011.2125935.
[6] Chuanyang Wang and Longya Xu. “Investigation of a practical approach for rotor
position estimation of PM machines without rotor saliency”. In: Power Electronics
and Motion Control Conference, 2000. Proceedings. IPEMC 2000. The Third Inter-
national. Vol. 1. Aug. 2000, 329–335 vol.1. doi: 10.1109/IPEMC.2000.885424.
[7] M. J. Corley and R. D. Lorenz. “Rotor position and velocity estimation for a salient-
pole permanent magnet synchronous machine at standstill and high speeds”. In:
IEEE Transactions on Industry Applications 34.4 (July 1998), pp. 784–789. issn:
0093-9994. doi: 10.1109/28.703973.
[8] R. V. Jacomini and E. Bim. “Sensorless rotor position based on MRAS observer for
doubly fed induction generator”. In: 2015 IEEE 13th Brazilian Power Electronics
Conference and 1st Southern Power Electronics Conference (COBEP/SPEC). Nov.
2015, pp. 1–6. doi: 10.1109/COBEP.2015.7420139.
[9] M. Tomita, T. Senjyu, S. Doki, et al. “New sensorless control for brushless DC motors
using disturbance observers and adaptive velocity estimations”. In: IEEE Transac-
tions on Industrial Electronics 45.2 (Apr. 1998), pp. 274–282. issn: 0278-0046. doi:
10.1109/41.681226.
[10] J. P. Johnson, M. Ehsani, and Y. Guzelgunler. “Review of sensorless methods for
brushless DC”. In: Industry Applications Conference, 1999. Thirty-Fourth IAS An-
nual Meeting. Conference Record of the 1999 IEEE. Vol. 1. 1999, 143–150 vol.1. doi:
10.1109/IAS.1999.799944.
60
[11] V. Hubik, M. Sveda, and V. Singule. “On the development of BLDC motor con-
trol run-up algorithms for aerospace application”. In: Power Electronics and Motion
Control Conference, 2008. EPE-PEMC 2008. 13th. Sept. 2008, pp. 1620–1624. doi:
10.1109/EPEPEMC.2008.4635499.
[12] S. Morimoto, K. Kawamoto, M. Sanada, et al. “Sensorless control strategy for salient-
pole PMSM based on extended EMF in rotating reference frame”. In: IEEE Trans-
actions on Industry Applications 38.4 (July 2002), pp. 1054–1061. issn: 0093-9994.
doi: 10.1109/TIA.2002.800777.
[13] S. Koonlaboon and S. Sangwongwanich. “Sensorless control of interior permanent-
magnet synchronous motors based on a fictitious permanent-magnet flux model”. In:
Fourtieth IAS Annual Meeting. Conference Record of the 2005 Industry Applications
Conference, 2005. Vol. 1. Oct. 2005, 311–318 Vol. 1. doi: 10 . 1109 / IAS . 2005 .
1518326.
[14] P. Kundur. Power System Stability and Control. McGraw Hill, 1994.
[15] H. Saadat. Power System Analysis. PSA Publishing, 2010.
[16] A. A. Fouad Paul M. Anderson. Power System Control and Stability. Wiley-IEEE
Press, 2002.
[17] Z. Gajić, D. Trišić, and S. Roxenborg. “Rotor DC current measurement by utilizing
current transformers”. In: Developments in Power System Protection (DPSP 2014),
12th IET International Conference on. Mar. 2014, pp. 1–6. doi: 10.1049/cp.2014.
0023.
61
TRITA 2016:149
ISSN 1653-5146
www.kth.se

Sync Generator Rotor Position Estimation

Uploaded by

Copyright:

Available Formats

Sync Generator Rotor Position Estimation

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Sync Generator Rotor Position Estimation

Uploaded by

Copyright:

Available Formats

DEGREE PROJECT IN ELECTRICAL ENGINEERING,

SECOND CYCLE, 30 CREDITS

Synchronous generator rotor

KTH ROYAL INSTITUTE OF TECHNOLOGY

ia Phase a current [A]

2 Theoretical background for the testing environment 3

3 Rotor position estimation 17

4 Designing the simulation environment and the estimator 25

6 Conclusion and future work 47

Figure 1.1: Layout of the SMIB system.

1.3 Outline of the report

Figure 2.1: Three phase synchronous machine [14].

2.1.1 Stator circuit equations

Ψa = −laa ia − lab ib − lac ic + laf d if d + la1d i1d + la1q i1q (2.4)

2.1.2 Stator self-inductances

MMFad(peak) =Na ia cos θ (2.5)

Φgad = (Na ia cos θ) Pd (2.7)

Where Pd and Pq are the d, respectively q-axis permeance coefficients.

The total air-gap flux linking phase a is

Φgaa =Φgad cos θ − Φgaq sin θ = Na ia (Pd cos2 θ + Pq sin2 θ)

The self inductance due to air-gap flux is

Similarly the b and c phases are expressed as follows

Therefore the mutual inductance between phase a and b can be written as

The mutual inductance between phase b-c and c-a are

lbc = lcb = − Lab0 − Lab2 cos(2θ − π) (2.17)

2.1.4 Mutual inductance between stator and rotor windings

laf d =Laf d cos θ (2.19)

2.1.5 Rotor circuit equations

The rotor circuit flux linkages are expressed as follows

2.2 Park’s transformation

2.2.1 Stator flux linkages in dq0 components

Ψd = − Ld id + Laf d if d + La1d i1d (2.33)

Where the newly defined inductances are are

2.2.2 Rotor flux linkages in dq0 components

2.3 Per unit representation

2.3.1 Per unit stator voltage and flux linkage equations

In per unit notation, equations (2.45)-(2.47) may be written as

Ψ̄d = − L̄d īd + L̄af d īf d + L̄a1d ī1d (2.51)

Ψ̄f d =L̄f f d īf d + L̄f 1d ī1d − L̄af d īd (2.57)

2.3.3 Complete set of electrical equations in per unit

Ψd = − (Lad + Ll )id + Lad if d + Lad i1d (2.63)

2.4 Derivation of machine parameters

Ψd (s) = − Ld id (s) + Lad if d (s) + Lad i1d (s) (2.74)

The rotor voltage equations are

ef d (s) =sΨf d (s) − Ψf d (0) + Rf d if d (s) (2.77)

Expressing the d-axis rotor voltage equations in incremental form yields

∆ef d (s) =sΨf d (s) + Rf d if d (s)

(R1d + sL11d )∆ef d (s) + sLad (R1d + sL1d ∆id (s))

−sLad ∆ef d (s) + sLad (Rf d + sLf d )∆id (s)

In factored form, equation (2.83) and (2.84) can be written as follows

(1 + sTd0 )(1 + sTd00 )

Similarly with the q-axis, the operational inductance is

2.4.1 Classical expression for time constants

2.4.2 Classical expression for transient and subtransient inductances

For a salient pole machine, the subtransient q-axis inductance is

2.5 Power system components

2.5.1 ∆−Y step up transformer

Based on observations in Simulink, an unfortunate property of the three phase transformer

Figure 2.4: Per phase equivalent circuit of a transformer [15].