- Email: [email protected]
Synchronous generator modelling and parameters estimation using least squares method
Simulation Modelling Practice and Theory 16 (2008) 678–689
Contents lists available at ScienceDirect
Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat
Synchronous generator modelling and parameters estimation using least squares method Emile Mouni *, Slim Tnani, Gérard Champenois University of Poitiers, Laboratoire d’Automatique et d’Informatique Industrielle, Bâtiment mécanique, 40, avenue du recteur Pineau 86022, Poitiers, France
a r t i c l e
i n f o
Article history: Received 19 April 2006 Received in revised form 17 March 2008 Accepted 8 April 2008 Available online 27 May 2008
Keywords: Synchronous generator Parameters estimation Short circuit test Park’s transformation State space modelling
a b s t r a c t In this paper, a technique for estimating the synchronous machine’s parameters using sudden short circuit test, is proposed. Before implementing estimation algorithms, a special method of the machine modelling is given. This last one allows to perform tests such as short-circuit, load impact and shedding test, in an easier way than the models usually developed in the literature. Thanks to the well known electrical equivalent circuit of the generator, the relationships between parameters generally used in the industry (i.e., reactances and time constants) and those used in researcher’s domains will be given. Finally, simulation results of the proposed method, allows to show that the algorithm is capable of providing very good estimated parameters fitting with the actual parameters. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction A number of modelling of synchronous machine methods have been already developed. With the increasing cost of detailed prototyping of electrical machine, it is becoming necessary to replace or supplement it with mathematical methods and computer simulation. Early works, see [1–6] have shown the crucial importance of a good model of synchronous machine taking into account dampers and other elements which are sometimes ignored or neglected in simplified modelling. In this paper, a new method of synchronous generator modelling taking into account an inside infinite resistance will be presented [7]. The particularity of such a modelling method is to make the performing of tests, usually used to validate or identify the machine, easy: short circuit test, load impact test, shedding test, etc. Once the synchronous machine has been constructed, manufacturers use programs based on various parameters (e.g., reactances and times constants) which are been graphically estimated to check the finest structural details of this one. To be in accordance with manufacturers methods, relationships between transient and sub transient quantities on one hand and the modelling parameters such as mutual and main inductances on the other hand,will be presented thank to the classical electrical circuits of the synchronous generator. In the last part of this work a numerical algorithm of synchronous machine parameters estimation based on least squares method, will be presented. A discussion on the simulation results will be done at the end of this work to validate the algorithm and the proposed modelling.
* Corresponding author. Tel.: +33 549453506. E-mail addresses: [email protected] (E. Mouni), [email protected] (S. Tnani), [email protected] (G. Champenois). 1569-190X/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.simpat.2008.04.005
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
679
d-axis q-axis
b
main eld D
Q
θe
c
a
Fig. 1. Synchronous generator windings with dampers.
2. Synchronous generator modelling Considering a synchronous generator with dampers at the rotor, the simplified scheme of the machine is given in Fig. 1: D and Q represent, respectively, d-axis and q-axis dampers. a, b and c are the three phases of the synchronous generator and he its electrical angle depending on the poles number. 2.1. Synchronous generator electrical equivalent circuit As we can see on the figure above, dampers are synthesized by short circuited inductances. From this figure and adopting generator convention, we can write machine equations in three axes frame as follows:
vabc ¼ r s iabc þ
d Wabc dt
d Wf dt d 0 ¼ rD iD þ WD dt d 0 ¼ rQ iQ þ WQ dt vf ¼ r f if þ
ð1Þ
where iD and iD are the direct and transverse dampers’ currents, WD and WQ are the direct and transverse dampers’ total flux, Wabc is stator total flux, Wf is the main field total flux. r s is the stator resistance, rf is the main field resistance, r D and r Q are the dampers resistances. The study will be done in Park’s frame thanks to Park’s matrix defined below with the electrical angle of the machine he :
Pðhe Þ ¼
cosðhe Þ
cosðhe 23pÞ
cosðhe þ 23pÞ
!
sinðhe Þ sinðhe 23pÞ sinðhe þ 23pÞ
such as:
Pðhe Þ vabc ¼ vdq
ð2Þ
Then, the global equation of the machine becomes
d Wd xe Wq dt d vq ¼ rs iq þ Wq þ xe Wd dt d vf ¼ r f if þ Wf dt d 0 ¼ rD iD þ WD dt d 0 ¼ rQ iQ þ WQ dt d J xe ¼ T e T r dt vd ¼ rs id þ
ð3Þ
T e is the electromechanical torque depending on machine current and given by [8]
Te ¼
3 pðWd iq Wq id Þ 2
ð4Þ
680
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
d-axis Ψσ sd
Ψσ dD
Damper
Ψσ f D
q-axis
{ΨΨσσDQ
Ψσ f
Ψσ sq 1 2
Legend 1 two poles rotor
2 stator armature Fig. 2. Synchronous generator leakage and linkage flux.
T r is resistant torque depending on the external load 2 p is the machine’s poles number and J is the machine inertia. In order to make the simulations in accordance with our future experimental conditions, the electrical speed is supposed to be constant. Indeed, in the test bench which is being achieved for validating the algorithms and control laws that we developed, the synchronous generator is involved by a DC motor. This last one is controlled by a motor drive: VNTC4075 from Alstom Company. Then, the state space equation of the machine will be given by using this assumption. In order determine synchronous generator equivalent circuit, a two salient poles machine will be considered. The scheme of this one is given in Fig. 2: where WrD and WrQ are the direct and transverse dampers leakage flux, Wrsd and Wrsq are the direct and transverse stator leakage flux, WrdD is the linkage flux between direct axis and direct dampers, WrfD linkage flux between main field and direct dampers, Wad and Waq are the direct and transverse main flux but implicitly omitted in Fig. 2. The effect of main field on the stator ðWrfs Þ is neglected, then the following relationships can be written
8 Wd ¼ Wad þ Wrsd þ WrdD ¼ lad ðid þ iD þ if Þ lrsd id þ lrdD ðiD id Þ > > > > > Wq ¼ Waq þ Wrsq ¼ laq ðiq þ iQ Þ lrsq iq > > > < Wf ¼ Wad þ Wrf þ WrfD ¼ lad ðid þ iD þ if Þ þ lrfD ðif þ iD Þ þ lrf if > > WD ¼ Wad þ WrfD þ WrdD þ WrD > > > > ¼ lad ðid þ iD þ if Þ þ lrfD ðif þ iD Þ þ lrD iD þ lrdD ðid þ iD Þ > > : WQ ¼ Waq þ WrQ ¼ laq ðiq þ iQ Þ þ lrQ iQ
ð5Þ
From these equations, we can deduce the synchronous generator electrical scheme (Fig. 3) [9,10]: where lrsd and lrsq are the direct and transverse stator leakage inductances, lrf is the main field leakage inductance, lad and laq are the direct and transverse stator main inductances, lrdD is the linkage inductance between stator d-axis and the direct damper, lrfD is the linkage inductance between rotor and the direct damper, lrD and lrQ are dampers leakage inductances.
rs
lσ sd lσ f D
rf
lσ f i f
lad
vd
rD
vf lσ D iD
lσ dD id ωe Ψq
rs
vq
lσ sq laq
iQ rQ lσ Q
ωe Ψd
iq
Fig. 3. d-axis and q-axis electrical equivalent circuits.
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
681
It is well known that reactance and inductance are linked by x ¼ lx; so thanks to the Eq. (5) we can deduce the following relationships between main and mutual inductances on one hand and reactances on the other hand:
8 x l ¼ xad ; laq ¼ xaqe ; ld ¼ xxde ¼ xad þxrxsdeþxrdD > > < ad xe x x þx rfD þxrdD lq ¼ xqe ¼ aqxe rsq ; lD ¼ xad þxrD þx xe > > x þx : l ¼ aq rQ ; l ¼ xad þxrf þxrfD ; m ¼ xad ; m ¼ xaq Q sQ f sf xe xe xe xe
ð6Þ
The relations between main reactances and machines parameters are:
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > < xad ¼ T 0d0 r f xe ðxd x0d Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > : xaq ¼ xq r Q xe ðT 00q0 T 00q Þ
ð7Þ
where xd is steady state reactance, x0d is the direct transient reactance, x00d is the direct sub transient reactance, x00q is the transverse sub transient q-reactance, T 0d is the direct transient time constant, T 00d is the direct sub transient time constant, T 0d0 is the open direct transient time constant, T 00q0 is the open transverse sub transient time constant, and xe is the machine electrical speed corresponding to the time derivative of he . Thus, from the expressions in (7), the machine can be represented by reactances and time constants and fluent relationships between parameters usually used in industry and those form academic domains can be deduced. 2.2. Synchronous machine modelling by state equations The model used in this paper is given by Fig. 4. In this modelling, a star-connected ‘‘infinite” resistance rin ð106 XÞ is incorporated. This allows ones to generate threephase voltage and then to create terminals A, B and C on which a three-phase load can be connected. On Fig. 4, vf and the output currents are used in the input vector. For this, inside currents ia , ib and ic are transformed into id and iq on one hand and load currents isa , isb and isc into idl and iql on the other hand. Consequently, the output voltage in Park’s framework can be expressed as
vd ¼ r in ðid idl Þ
ð8Þ
vq ¼ r in ðiq iql Þ with
idl iql
0
1 isa B C ¼ Pðhe Þ@ isb A isc
Finally the global equation is
1 0 1 0 1 id id rin idl C B B C B C i i r i B in ql C B qC B qC C B B C B C B vf C ¼ R B if C þ M a d B if C C B B C B C dt C B B C B C @ 0 A @ iD A @ iD A 0
0
iQ
ð9Þ
iQ
isa
ia rin vf
rin
isb
A
B
rin isc
ic x = Ax + Bu y = Cx + Du Fig. 4. Synchronous generator with infinite inside load.
C
682
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
where
0
r s r in
B B ld xe B R¼B 0 B B 0 @ 0
l q xe r s r in
0
0
xe msf xe msD
xe msQ 0
0
rf
0
0
0
0
rD
0
0
0
0
rQ
1
0
C C C C; C C A
B B 0 B Ma ¼ B B msf B @ msD 0
ld
0
msf
msD
lq
0
0
0
lf
mfD
0
mfD
lD
msQ
0
0
0
1
C msQ C C 0 C C C 0 A lQ
Then, the synchronous generator state-space equation is given by
x_ ¼ Ax þ Bu
ð10Þ
y ¼ Cx þ Du
rin 0 0 0 0 where A ¼ M 1 is the state matrix, B ¼ M1 C¼ is the observation matrix,D ¼ a R a , 0 rin 0 0 0 T rin 0 0 0 T is the state vector, y ¼ ð vd vq Þ is the output vector and , x ¼ id iq if iD iQ 0 rin 0 0 T u ¼ idl iql vf 0 0 is the exogenous inputs vector containing the excitation vector vf . 3. Sudden short circuit principle parameters determination The aim of the modelling above is to perform some validation tests. The one that we will perform in this paper is the sudden short circuit. The principle is described below: The machine is involved at the rated speed without load until the system reaches the steady state. During this steady state, a short circuit is performed on its three phases and then, currents and voltage are measured. This test allows to determine synchronous machine parameters and then to validate or not the achieved model. Below are IEEE’s recommendations according to direct and transverse axes of Park’s framework. 3.1. Direct axis parameters measurement After the performing of sudden short circuit, the current on each phase can be described as following:
1 1 1 t 1 1 t : exp 0 þ 00 0 exp 00 cosðx t þ h0 Þ þ V m þ 0 xd xd xd xd xd Td Td ! ! " # 1 1 t 1 1 t : exp : exp cosðh cosð2 þ Þ þ x t þ h Þ 0 0 x00d x00q Ta x00d x00q Ta
i ¼ Vm
ð11Þ
where V m is the voltage maximum value prior the short circuit applying. Notice that the above expression can be simplified in considering the current which aperiodic component is null (h0 ¼ p2 . If the three currents contain aperiodic part, a little manipulation can be used to eliminate it. This manipulation consists in subtracting the exponential curve included in the short circuit current and representing the aperiodic component’s contribution. Among the parameters that we used in the simulation, transverse sub transient q-reactance x00q and direct sub transient d-reactance x00d are the same. Definitely the Eq. (11) becomes:
1 1 1 t 1 1 t exp 0 þ 00 0 exp 00 : cosðx t þ h0 Þ þ 0 xd xd xd xd xd Td Td
envelope phase current
1000 Current (Amperes)
i ¼ Vm
500 0 –500 –1000 1
2 time (sec.)
Fig. 5. Current and envelope.
3
ð12Þ
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
683
The IEEE standards, [12,13], recommend to draw an envelope which fits the best with output current. Then calculations are made with this last one. The figure below shows this draw applied to simulated generator and corresponds to a current without aperiodic component (see Fig. 5). The envelope equation is obtained in considering the current peaks. Thus the Eq. (12) becomes:
ienv ¼ V m
1 1 1 t 1 1 t : exp 0 þ 00 0 : exp 00 þ 0 xd xd xd xd xd Td Td
ð13Þ
3.1.1. Direct steady state reactance This reactance is easy to calculate. It corresponds to the reactance of the machine when it works at steady state. In the Eq. (13), as the exponential functions are decreasing, direct steady state reactance can be deduced as following:
xd ¼
Vm isteady
ð14Þ
3.1.2. Direct transient reactance and time constant Once the steady current value is found, it is subtract to ienv as presented below:
1 1 t 1 1 t : exp þ : exp ienv is ¼ V m : x0d xd x00d x0d T 0d T 00d
ð15Þ
The IEEE standards recommendations are about the use of semi logarithmic frame to determine reactances and time constants. These recommendations allow to go from an exponential curve to a sum of real straight curves. This leads to have easier calculations. Nevertheless an assumption is made on the time constants: 00 Assumption: Transient time constant is very high besides sub transient time constant ðT 0d Td Þ. Consequently The sub transient component decreases quickly in relation to the transient one. So its effects can be neglected from a certain time. This assumption leads to make an approximation of the above current difference ienv is . That means:
1 1 t : exp ienv is V m : 0 x0d xd Td
ð16Þ
Using semi logarithmic method, we can say that it exists two quantities A and B such as
lnðienv is Þ A t þ B
ð17Þ
A is the slope and B the value at the frame origin. The transient parameters can be then obtained by solving the following equations system:
8 0 < T d ¼ A1
: ln V m : x10 x1 ¼B
ð18Þ
d
d
3.1.3. Direct sub-transient reactance and time constant The last part of this section is about sub transient parameters calculation. The method is the same as described above. Therefore, the approximation is made on the quantity below:
lnðienv is itrans Þ A0 t þ B0
ð19Þ
where itrans is the transient current calculated above with steady and transient parameters. This leads to the following new system of equations:
8 00 < T d ¼ A10
: ln V m x100 x10 ¼ B0 d
ð20Þ
d
Note that, several other methods are used to determine direct axis parameters, see [14–16]. The one presented in this paper is easy to implement and the results we got are satisfactory. 3.2. Transverse axis parameters measurement To determine q-axis parameters, a Park transformation on currents is used. In this part, only q-axis current is used. The figure below shows this current after sudden short circuit application: (see Fig. 6) According to IEEE standards the shape above can be represented by the following equation.
iq ¼
Vm t sinðx t þ h0 Þ exp Ta x00q
ð21Þ
684
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
1500 transverse current 1000
q
i (A)
500 0 –500 –1000 20
20.05
20.1
20.15
20.2
Time (sec) Fig. 6. Current on q-axis.
We can use the same method as described above, but it is also simple to consider some peaks and to solve an equation system as following
8
> < iq1 ¼ Vxm00 exp Tt1a q
> : iq2 ¼ Vxm00 exp Tt2 a
ð22Þ
q
Because of the easiness of the calculation of transverse parameters, we will only apply the least squares method to the direct parameters. 3.3. Time constants in open circuit These time constants are deduced from the above calculations. According to the IEEE standards, to calculate open circuit time constants, the following relationships are used: Open circuit direct transient.
xd T 0d0 ¼ 0 x0d Td
ð23Þ
Open circuit direct sub transient.
x0d T 00d0 ¼ 00 x00d Td
ð24Þ
Open circuit transverse sub-transient. 00 xq T q0 ¼ x00q T 00q
ð25Þ
4. Parameters estimation using least squares method The Fig. 7 above shows the strategy used to implement the algorithm. Initially, a sudden short circuit test is performed on the synchronous machine and the measured outputs are obtained. These data are used to make an estimation of synchronous machine parameters. The estimation method is based on the IEEE standards described above. Each set of parameters is used to simulate a new system. From this last one an error between the actual model and the simulated model is calculated such as (for the kth set of parameters for instance):
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX JðkÞ ¼ t ðiactual ðiÞ isim ði; kÞÞ2
ð26Þ
i¼1
where k is the iteration order or set of parameters order, iactual ðiÞ is the actual current value for ith sampling point, isim ði; kÞ is the simulated current value for ith sampling point with the kth set of parameters, N is the number of sampling points, JðkÞ the criterion value with the kth set of parameters.
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
685
Test data loading Parameters Estimation
Simulated System
no Least squares Algorithm
Recording
end? yes
Minimum estimation Fig. 7. Parameters estimation framework.
For each iteration the criterion value, extended to the whole range of variation is stored in a vector. The operation is done again until the instruction end is reached. From there, the vector J allows to get the optimal set of parameters. The following flow chart given by Fig. 8 explains how the algorithm is performed: First of all, initialization is done (imax , jmax , i ¼ 0; . . .), then steady reactance is calculated and some ranges are defined to determine transient set of parameters. These ranges are chosen far enough from short circuit start point and their width depends on the machine and the simulation step. For instance, the simulation we performed uses 1000 samples per transient range and the simulation step is 0:25 ms. As regards to sub transient ranges they are chosen near the start point. Fig. 9 shows the subdivisions made on the current envelope. From the Fig. 8, we can notice that for each transient range, jmax sub transient parameters are calculated. At the end of the process, the criterion J is a vector which length is imax jmax . Each value of this vector corresponds to a particular set of parameters. Minimizing this criterion allows to get the parameters which fit the best with accurate values. When this criterion is plotted, this leads to Fig. 10. The Fig. 10 shows the different parts corresponding to transient ranges. For each transient range, a minimum can be found as specified below by Fig. 11. The algorithm, we elaborated, works in order to find a particular value kopt corresponding to the general minimum such as:
8vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi9
ð27Þ
k2K
where K is the iterations group which length is defined above ðimax jmax Þ, N the envelope length. When this optimum iteration ðkopt Þ is found, the different validation curves can be plotted to check whether actual quantities and estimated ones match each other. 5. Simulation results and discussion Once the state space modelling with inside high enough load is done, the estimation algorithm based on the least squares method is implemented. The data used to perform this algorithm are shown in the Table 1. The machine involved in the modelling presents the following characteristics: Rated phase to phase voltage U n ¼ 530 V, Rated current In ¼ 243 A. The rated impedance of the synchronous machine zn can then be deduced by
Un zn ¼ pffiffiffiffiffiffiffi ð3Þ In
ð28Þ
This last value allows to convert reactances in per units (p.u). 5.1. Validation of machine modelling The implementation of the synchronous machine is done by SimulinkTM/S-function. Before performing a short circuit on the simulated machine, a load test is achieved. The results, corresponding to the steady state, are presented below (see Figs. 12 and 13. The currents and the voltages shown by the figures above, are sinusoidal and well balanced.
686
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
Initialization jmax imax k 0 i=0 xd
i
No
imax ?
end of the process
Yes j=0 xd
i Tdo
i 1 Td
j=j+1 xd
k
j
Tdo
Td ”
k 1 J(k) j
jmax
jmax
Fig. 8. Flow chart for parameters estimation.
envelope (Amperes) j
1500 j
jmax
1 sub transient ranges
i i
1 Transient imax ranges
isteady 2.4 0.4 Time sec Fig. 9. Transient and sub transient ranges subdivision.
5.2. Parameters estimation results As mentioned above, the simulated machine is involve until the steady state is reached. Then a sudden short circuit is applied on its three phases. After this short circuit currents and voltages are recorded and the developed algorithm allows to obtain reactances and time constants given in Table 2: The error we calculated is referred to the accurate parameters values in Table 1 by using the following relation
error ¼
ðaccurate valueÞ ðsimulated valueÞ accurate value
ð29Þ
687
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
6
x 10 4
5
criterion J
4 3
zoom area
2 1 0 0
500 1000 iterations (k)
jmax
1500
Fig. 10. Criterion curve.
2500
criterion J
2000
1500
1000 optimum
500 0
200
400
600 800 1000 1200 1400 iterations (k)
Fig. 11. Zoom on criterion curve.
Table 1 Table of Synchronous machine parameters Parameters
xd (%)
x0d (%)
x00d (%)
xq (%)
x00q (%)
T 0d0 (s)
T 0d (s)
T 00d0 (ms)
T 00d (ms)
T 00q0 (ms)
T 00q (ms)
Accurate values
176
36.1
26.5
150
26.5
3.12
0.64
41
30
41
7
s = second; ms = millisecond.
400
va vb vc
300 Voltage (V)
200 100 0 –100 –200 –300 –400 9.985
9.99
9.995 10 10.005 10.01 Time (Sec.)
Fig. 12. Steady state output voltage in abc frame.
688
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
300
ia ib ic
Current (A)
200 100 0 –100 –200 –300 9.985
9.99
9.995 10 Time(sec.)
10.005
10.01
Fig. 13. Steady state output current in abc frame.
Table 2 Table of synchronous machine estimated parameters Parameters
xd (%)
x0d (%)
x00d (%)
x00q (%)
T 0d0 (s)
T 0d (s)
T 00d0 (ms)
T 00d (ms)
Simulated values Error%
175.7 0.09
34.5 4.1
23 12.9
26.8 1.5
3.28 5.1
0.645 0.8
39 4.9
26 13.3
accurate envelope estimated envelope
2000
ienv (A)
1500
1000
500 0
2
4 6 time (sec.)
8
10
Fig. 14. Actual and simulated current envelopes.
5.3. Discussion on results The Table 2 shows that the estimation of reactances by least squares method provide good results. The precision is such as the error is less than 0.1% in the case of transient reactance. Neverthless, an error of 12.9% is found in the direct subtransient reactance estimation. Even if this last value is high, it is less than those obtained with the traditional method which is widely graphical and which errors sometimes reach 20% or 25% of accurate values, see [17]. As regards the estimated time constants, they match very well with the accurate parameters. Even the sub transient time constant which is difficult to determine with the traditional method, is satisfying. For this last value we got 26 ms instead of 30 ms. Globally, the obtained results are in accordance with what we are supposed to get. Below is the reconstructed envelope from the estimated quantities (i.e., reactances and time constants). In this figure, a comparison with the actual short circuit envelope is done: As we can see on the Fig. 14, there is an accordance between actual signal and the estimated one. This figure reinforces the conclusions we made from Table 2 on the effectiveness of the algorithm we proposed. 6. Conclusion In this paper, a complete modelling of synchronous machine taking in account of the existence of dampers and using only reactances and time constants as parameters is presented. The Park’s framework is used and the modelling is performed with
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
689
the state space modelling. The particularity of this one is the use of an inside infinite resistance which allows to create terminals on the model. Therefore a series of external loads can be connected and disconnected without modifying the proper structure of the model. This model is handier than those found in literature which usually include load [7,11]. Indeed, thanks to such a model, tests often performed on the machine to validate it [12] or to identify it [1], become easier to achieve.In the second part of this paper, a statistical technique for determining synchronous machine parameters is proposed. This method is based on only electrical quantities (currents and voltages), so it could be easily achieved by using standard equipments. The method is very efficient and the obtained results suit strongly with those that we were supposed to get. Thus, the method presented in this paper allows not only to determine synchronous machine parameters but also to validate the model built on Matlab/SimulinkTM ). References [1] A. Launay-Querré, J. Gabano, G. Champenois, Parameter estimation of a synchronous generator, Electrimacs, MontrTal, Canada. [2] F. Bruck, F. Himmelstoss, Modelling and simulation of synchronous machine, in: IEEE Computers in Power Electronics, 6th Workshop on 19–22 July, 1998, pp. 81–86. [3] P. Pillay, R. Krishnan, Modeling, simulation, and analysis of permanent-magnet motor drives. Part 1: The permanent-magnet synchronous motor drive, IEEE Transactions on Industry Application 25 (1989). [4] T. Sebastien, G.R. Slemon, M. Rahman, Modelling of permanent synchronous motors, IEEE Transactions on Magnetics Mag 22 (5) (1986) 1069–1071. [5] K. Srivasta, B. Berggren, Simulation of synchronous machine in phase coordinates including magnetic saturation, Elsevier Science Electric Power Systems Research 56 (2000) 177–283. [6] I. Kamwa, P. Viarouge, E. Dickinson, Identification of generalised models of synchronous machines from time domain tests, IEEE Proceedings-c 138 (6) (1991) 485–498. [7] E. Mouni, S. Tnani, G. Champenois, Comparative study of three modelling methods of synchronous generator, in: Conference of the IEEE Industrial Electronics Society, Paris, France, 2006. [8] P.C. Krause, oleg Wasynczuk, S.D. Sudhoff, Analysis of Electric Machinery, IEEE Press, 1994. [9] J. Verbeeck, Standstill frequency response measurement and identification methods for synchronous machines, Thesis, Vrije Unversiteit Brussel, 2000. [10] H. Guesbaoui, C. Iung, Reduced models and characteristic parameters of the synchronous machine obtained by a multi-time scale simplification, in: IEEE Industry Applications Society Annual Meeting, 1994, Conference Record of the 1994 IEEE 2–6 October 1994, vol. 1, pp. 262–271. [11] J.D. Gabano, G. Champenois, Synchronous generator modeling using a non-steady state park model, in: 6th International Symposium on Advanced Electromechanical Motion Systems, Switzerland, 2005. [12] IEEE: Std 1158-1995, IEEE guide: Test procedures for synchronous machines, IEEE Standards, 1995. [13] J.D. Ghali, Détermination des paramétres de la machine synchrone, Ecole Polytechnique Fédérale de Lausanne, 1994. [14] S. Moreau, R. Kahoul, J. Louis, Parameters estimation of permanent magnet synchronous machine without adding extra signal as input excitation, ISIE’04 Ajaccio. [15] F.L. Alvarado, C. Cañizares, Synchronous machine parameters from sudden short tests by back solving, IEEE Transactions on Energy Conversion 4 (2) (1989) 4. June. [16] A. Tumageanian, A. Keyhani, Identification of synchronous machine linear parameters from standstill step voltage input data, IEEE Transactions on Energy Conversion 10 (2) (1995) 10. June. [17] C. Xingang, J. Luming, W. Xusheng, Y. Kaisheng, W. Zhifei, A new approach to determine parameters of synchronous machine using wavelet transform and Prony algorithm, in: 2004 International Conference on Power System Technology – POWERCON 2004 Singapore.
Contents lists available at ScienceDirect
Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat
Synchronous generator modelling and parameters estimation using least squares method Emile Mouni *, Slim Tnani, Gérard Champenois University of Poitiers, Laboratoire d’Automatique et d’Informatique Industrielle, Bâtiment mécanique, 40, avenue du recteur Pineau 86022, Poitiers, France
a r t i c l e
i n f o
Article history: Received 19 April 2006 Received in revised form 17 March 2008 Accepted 8 April 2008 Available online 27 May 2008
Keywords: Synchronous generator Parameters estimation Short circuit test Park’s transformation State space modelling
a b s t r a c t In this paper, a technique for estimating the synchronous machine’s parameters using sudden short circuit test, is proposed. Before implementing estimation algorithms, a special method of the machine modelling is given. This last one allows to perform tests such as short-circuit, load impact and shedding test, in an easier way than the models usually developed in the literature. Thanks to the well known electrical equivalent circuit of the generator, the relationships between parameters generally used in the industry (i.e., reactances and time constants) and those used in researcher’s domains will be given. Finally, simulation results of the proposed method, allows to show that the algorithm is capable of providing very good estimated parameters fitting with the actual parameters. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction A number of modelling of synchronous machine methods have been already developed. With the increasing cost of detailed prototyping of electrical machine, it is becoming necessary to replace or supplement it with mathematical methods and computer simulation. Early works, see [1–6] have shown the crucial importance of a good model of synchronous machine taking into account dampers and other elements which are sometimes ignored or neglected in simplified modelling. In this paper, a new method of synchronous generator modelling taking into account an inside infinite resistance will be presented [7]. The particularity of such a modelling method is to make the performing of tests, usually used to validate or identify the machine, easy: short circuit test, load impact test, shedding test, etc. Once the synchronous machine has been constructed, manufacturers use programs based on various parameters (e.g., reactances and times constants) which are been graphically estimated to check the finest structural details of this one. To be in accordance with manufacturers methods, relationships between transient and sub transient quantities on one hand and the modelling parameters such as mutual and main inductances on the other hand,will be presented thank to the classical electrical circuits of the synchronous generator. In the last part of this work a numerical algorithm of synchronous machine parameters estimation based on least squares method, will be presented. A discussion on the simulation results will be done at the end of this work to validate the algorithm and the proposed modelling.
* Corresponding author. Tel.: +33 549453506. E-mail addresses: [email protected] (E. Mouni), [email protected] (S. Tnani), [email protected] (G. Champenois). 1569-190X/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.simpat.2008.04.005
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
679
d-axis q-axis
b
main eld D
Q
θe
c
a
Fig. 1. Synchronous generator windings with dampers.
2. Synchronous generator modelling Considering a synchronous generator with dampers at the rotor, the simplified scheme of the machine is given in Fig. 1: D and Q represent, respectively, d-axis and q-axis dampers. a, b and c are the three phases of the synchronous generator and he its electrical angle depending on the poles number. 2.1. Synchronous generator electrical equivalent circuit As we can see on the figure above, dampers are synthesized by short circuited inductances. From this figure and adopting generator convention, we can write machine equations in three axes frame as follows:
vabc ¼ r s iabc þ
d Wabc dt
d Wf dt d 0 ¼ rD iD þ WD dt d 0 ¼ rQ iQ þ WQ dt vf ¼ r f if þ
ð1Þ
where iD and iD are the direct and transverse dampers’ currents, WD and WQ are the direct and transverse dampers’ total flux, Wabc is stator total flux, Wf is the main field total flux. r s is the stator resistance, rf is the main field resistance, r D and r Q are the dampers resistances. The study will be done in Park’s frame thanks to Park’s matrix defined below with the electrical angle of the machine he :
Pðhe Þ ¼
cosðhe Þ
cosðhe 23pÞ
cosðhe þ 23pÞ
!
sinðhe Þ sinðhe 23pÞ sinðhe þ 23pÞ
such as:
Pðhe Þ vabc ¼ vdq
ð2Þ
Then, the global equation of the machine becomes
d Wd xe Wq dt d vq ¼ rs iq þ Wq þ xe Wd dt d vf ¼ r f if þ Wf dt d 0 ¼ rD iD þ WD dt d 0 ¼ rQ iQ þ WQ dt d J xe ¼ T e T r dt vd ¼ rs id þ
ð3Þ
T e is the electromechanical torque depending on machine current and given by [8]
Te ¼
3 pðWd iq Wq id Þ 2
ð4Þ
680
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
d-axis Ψσ sd
Ψσ dD
Damper
Ψσ f D
q-axis
{ΨΨσσDQ
Ψσ f
Ψσ sq 1 2
Legend 1 two poles rotor
2 stator armature Fig. 2. Synchronous generator leakage and linkage flux.
T r is resistant torque depending on the external load 2 p is the machine’s poles number and J is the machine inertia. In order to make the simulations in accordance with our future experimental conditions, the electrical speed is supposed to be constant. Indeed, in the test bench which is being achieved for validating the algorithms and control laws that we developed, the synchronous generator is involved by a DC motor. This last one is controlled by a motor drive: VNTC4075 from Alstom Company. Then, the state space equation of the machine will be given by using this assumption. In order determine synchronous generator equivalent circuit, a two salient poles machine will be considered. The scheme of this one is given in Fig. 2: where WrD and WrQ are the direct and transverse dampers leakage flux, Wrsd and Wrsq are the direct and transverse stator leakage flux, WrdD is the linkage flux between direct axis and direct dampers, WrfD linkage flux between main field and direct dampers, Wad and Waq are the direct and transverse main flux but implicitly omitted in Fig. 2. The effect of main field on the stator ðWrfs Þ is neglected, then the following relationships can be written
8 Wd ¼ Wad þ Wrsd þ WrdD ¼ lad ðid þ iD þ if Þ lrsd id þ lrdD ðiD id Þ > > > > > Wq ¼ Waq þ Wrsq ¼ laq ðiq þ iQ Þ lrsq iq > > > < Wf ¼ Wad þ Wrf þ WrfD ¼ lad ðid þ iD þ if Þ þ lrfD ðif þ iD Þ þ lrf if > > WD ¼ Wad þ WrfD þ WrdD þ WrD > > > > ¼ lad ðid þ iD þ if Þ þ lrfD ðif þ iD Þ þ lrD iD þ lrdD ðid þ iD Þ > > : WQ ¼ Waq þ WrQ ¼ laq ðiq þ iQ Þ þ lrQ iQ
ð5Þ
From these equations, we can deduce the synchronous generator electrical scheme (Fig. 3) [9,10]: where lrsd and lrsq are the direct and transverse stator leakage inductances, lrf is the main field leakage inductance, lad and laq are the direct and transverse stator main inductances, lrdD is the linkage inductance between stator d-axis and the direct damper, lrfD is the linkage inductance between rotor and the direct damper, lrD and lrQ are dampers leakage inductances.
rs
lσ sd lσ f D
rf
lσ f i f
lad
vd
rD
vf lσ D iD
lσ dD id ωe Ψq
rs
vq
lσ sq laq
iQ rQ lσ Q
ωe Ψd
iq
Fig. 3. d-axis and q-axis electrical equivalent circuits.
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
681
It is well known that reactance and inductance are linked by x ¼ lx; so thanks to the Eq. (5) we can deduce the following relationships between main and mutual inductances on one hand and reactances on the other hand:
8 x l ¼ xad ; laq ¼ xaqe ; ld ¼ xxde ¼ xad þxrxsdeþxrdD > > < ad xe x x þx rfD þxrdD lq ¼ xqe ¼ aqxe rsq ; lD ¼ xad þxrD þx xe > > x þx : l ¼ aq rQ ; l ¼ xad þxrf þxrfD ; m ¼ xad ; m ¼ xaq Q sQ f sf xe xe xe xe
ð6Þ
The relations between main reactances and machines parameters are:
8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > < xad ¼ T 0d0 r f xe ðxd x0d Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > : xaq ¼ xq r Q xe ðT 00q0 T 00q Þ
ð7Þ
where xd is steady state reactance, x0d is the direct transient reactance, x00d is the direct sub transient reactance, x00q is the transverse sub transient q-reactance, T 0d is the direct transient time constant, T 00d is the direct sub transient time constant, T 0d0 is the open direct transient time constant, T 00q0 is the open transverse sub transient time constant, and xe is the machine electrical speed corresponding to the time derivative of he . Thus, from the expressions in (7), the machine can be represented by reactances and time constants and fluent relationships between parameters usually used in industry and those form academic domains can be deduced. 2.2. Synchronous machine modelling by state equations The model used in this paper is given by Fig. 4. In this modelling, a star-connected ‘‘infinite” resistance rin ð106 XÞ is incorporated. This allows ones to generate threephase voltage and then to create terminals A, B and C on which a three-phase load can be connected. On Fig. 4, vf and the output currents are used in the input vector. For this, inside currents ia , ib and ic are transformed into id and iq on one hand and load currents isa , isb and isc into idl and iql on the other hand. Consequently, the output voltage in Park’s framework can be expressed as
vd ¼ r in ðid idl Þ
ð8Þ
vq ¼ r in ðiq iql Þ with
idl iql
0
1 isa B C ¼ Pðhe Þ@ isb A isc
Finally the global equation is
1 0 1 0 1 id id rin idl C B B C B C i i r i B in ql C B qC B qC C B B C B C B vf C ¼ R B if C þ M a d B if C C B B C B C dt C B B C B C @ 0 A @ iD A @ iD A 0
0
iQ
ð9Þ
iQ
isa
ia rin vf
rin
isb
A
B
rin isc
ic x = Ax + Bu y = Cx + Du Fig. 4. Synchronous generator with infinite inside load.
C
682
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
where
0
r s r in
B B ld xe B R¼B 0 B B 0 @ 0
l q xe r s r in
0
0
xe msf xe msD
xe msQ 0
0
rf
0
0
0
0
rD
0
0
0
0
rQ
1
0
C C C C; C C A
B B 0 B Ma ¼ B B msf B @ msD 0
ld
0
msf
msD
lq
0
0
0
lf
mfD
0
mfD
lD
msQ
0
0
0
1
C msQ C C 0 C C C 0 A lQ
Then, the synchronous generator state-space equation is given by
x_ ¼ Ax þ Bu
ð10Þ
y ¼ Cx þ Du
rin 0 0 0 0 where A ¼ M 1 is the state matrix, B ¼ M1 C¼ is the observation matrix,D ¼ a R a , 0 rin 0 0 0 T rin 0 0 0 T is the state vector, y ¼ ð vd vq Þ is the output vector and , x ¼ id iq if iD iQ 0 rin 0 0 T u ¼ idl iql vf 0 0 is the exogenous inputs vector containing the excitation vector vf . 3. Sudden short circuit principle parameters determination The aim of the modelling above is to perform some validation tests. The one that we will perform in this paper is the sudden short circuit. The principle is described below: The machine is involved at the rated speed without load until the system reaches the steady state. During this steady state, a short circuit is performed on its three phases and then, currents and voltage are measured. This test allows to determine synchronous machine parameters and then to validate or not the achieved model. Below are IEEE’s recommendations according to direct and transverse axes of Park’s framework. 3.1. Direct axis parameters measurement After the performing of sudden short circuit, the current on each phase can be described as following:
1 1 1 t 1 1 t : exp 0 þ 00 0 exp 00 cosðx t þ h0 Þ þ V m þ 0 xd xd xd xd xd Td Td ! ! " # 1 1 t 1 1 t : exp : exp cosðh cosð2 þ Þ þ x t þ h Þ 0 0 x00d x00q Ta x00d x00q Ta
i ¼ Vm
ð11Þ
where V m is the voltage maximum value prior the short circuit applying. Notice that the above expression can be simplified in considering the current which aperiodic component is null (h0 ¼ p2 . If the three currents contain aperiodic part, a little manipulation can be used to eliminate it. This manipulation consists in subtracting the exponential curve included in the short circuit current and representing the aperiodic component’s contribution. Among the parameters that we used in the simulation, transverse sub transient q-reactance x00q and direct sub transient d-reactance x00d are the same. Definitely the Eq. (11) becomes:
1 1 1 t 1 1 t exp 0 þ 00 0 exp 00 : cosðx t þ h0 Þ þ 0 xd xd xd xd xd Td Td
envelope phase current
1000 Current (Amperes)
i ¼ Vm
500 0 –500 –1000 1
2 time (sec.)
Fig. 5. Current and envelope.
3
ð12Þ
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
683
The IEEE standards, [12,13], recommend to draw an envelope which fits the best with output current. Then calculations are made with this last one. The figure below shows this draw applied to simulated generator and corresponds to a current without aperiodic component (see Fig. 5). The envelope equation is obtained in considering the current peaks. Thus the Eq. (12) becomes:
ienv ¼ V m
1 1 1 t 1 1 t : exp 0 þ 00 0 : exp 00 þ 0 xd xd xd xd xd Td Td
ð13Þ
3.1.1. Direct steady state reactance This reactance is easy to calculate. It corresponds to the reactance of the machine when it works at steady state. In the Eq. (13), as the exponential functions are decreasing, direct steady state reactance can be deduced as following:
xd ¼
Vm isteady
ð14Þ
3.1.2. Direct transient reactance and time constant Once the steady current value is found, it is subtract to ienv as presented below:
1 1 t 1 1 t : exp þ : exp ienv is ¼ V m : x0d xd x00d x0d T 0d T 00d
ð15Þ
The IEEE standards recommendations are about the use of semi logarithmic frame to determine reactances and time constants. These recommendations allow to go from an exponential curve to a sum of real straight curves. This leads to have easier calculations. Nevertheless an assumption is made on the time constants: 00 Assumption: Transient time constant is very high besides sub transient time constant ðT 0d Td Þ. Consequently The sub transient component decreases quickly in relation to the transient one. So its effects can be neglected from a certain time. This assumption leads to make an approximation of the above current difference ienv is . That means:
1 1 t : exp ienv is V m : 0 x0d xd Td
ð16Þ
Using semi logarithmic method, we can say that it exists two quantities A and B such as
lnðienv is Þ A t þ B
ð17Þ
A is the slope and B the value at the frame origin. The transient parameters can be then obtained by solving the following equations system:
8 0 < T d ¼ A1
: ln V m : x10 x1 ¼B
ð18Þ
d
d
3.1.3. Direct sub-transient reactance and time constant The last part of this section is about sub transient parameters calculation. The method is the same as described above. Therefore, the approximation is made on the quantity below:
lnðienv is itrans Þ A0 t þ B0
ð19Þ
where itrans is the transient current calculated above with steady and transient parameters. This leads to the following new system of equations:
8 00 < T d ¼ A10
: ln V m x100 x10 ¼ B0 d
ð20Þ
d
Note that, several other methods are used to determine direct axis parameters, see [14–16]. The one presented in this paper is easy to implement and the results we got are satisfactory. 3.2. Transverse axis parameters measurement To determine q-axis parameters, a Park transformation on currents is used. In this part, only q-axis current is used. The figure below shows this current after sudden short circuit application: (see Fig. 6) According to IEEE standards the shape above can be represented by the following equation.
iq ¼
Vm t sinðx t þ h0 Þ exp Ta x00q
ð21Þ
684
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
1500 transverse current 1000
q
i (A)
500 0 –500 –1000 20
20.05
20.1
20.15
20.2
Time (sec) Fig. 6. Current on q-axis.
We can use the same method as described above, but it is also simple to consider some peaks and to solve an equation system as following
8
> < iq1 ¼ Vxm00 exp Tt1a q
> : iq2 ¼ Vxm00 exp Tt2 a
ð22Þ
q
Because of the easiness of the calculation of transverse parameters, we will only apply the least squares method to the direct parameters. 3.3. Time constants in open circuit These time constants are deduced from the above calculations. According to the IEEE standards, to calculate open circuit time constants, the following relationships are used: Open circuit direct transient.
xd T 0d0 ¼ 0 x0d Td
ð23Þ
Open circuit direct sub transient.
x0d T 00d0 ¼ 00 x00d Td
ð24Þ
Open circuit transverse sub-transient. 00 xq T q0 ¼ x00q T 00q
ð25Þ
4. Parameters estimation using least squares method The Fig. 7 above shows the strategy used to implement the algorithm. Initially, a sudden short circuit test is performed on the synchronous machine and the measured outputs are obtained. These data are used to make an estimation of synchronous machine parameters. The estimation method is based on the IEEE standards described above. Each set of parameters is used to simulate a new system. From this last one an error between the actual model and the simulated model is calculated such as (for the kth set of parameters for instance):
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX JðkÞ ¼ t ðiactual ðiÞ isim ði; kÞÞ2
ð26Þ
i¼1
where k is the iteration order or set of parameters order, iactual ðiÞ is the actual current value for ith sampling point, isim ði; kÞ is the simulated current value for ith sampling point with the kth set of parameters, N is the number of sampling points, JðkÞ the criterion value with the kth set of parameters.
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
685
Test data loading Parameters Estimation
Simulated System
no Least squares Algorithm
Recording
end? yes
Minimum estimation Fig. 7. Parameters estimation framework.
For each iteration the criterion value, extended to the whole range of variation is stored in a vector. The operation is done again until the instruction end is reached. From there, the vector J allows to get the optimal set of parameters. The following flow chart given by Fig. 8 explains how the algorithm is performed: First of all, initialization is done (imax , jmax , i ¼ 0; . . .), then steady reactance is calculated and some ranges are defined to determine transient set of parameters. These ranges are chosen far enough from short circuit start point and their width depends on the machine and the simulation step. For instance, the simulation we performed uses 1000 samples per transient range and the simulation step is 0:25 ms. As regards to sub transient ranges they are chosen near the start point. Fig. 9 shows the subdivisions made on the current envelope. From the Fig. 8, we can notice that for each transient range, jmax sub transient parameters are calculated. At the end of the process, the criterion J is a vector which length is imax jmax . Each value of this vector corresponds to a particular set of parameters. Minimizing this criterion allows to get the parameters which fit the best with accurate values. When this criterion is plotted, this leads to Fig. 10. The Fig. 10 shows the different parts corresponding to transient ranges. For each transient range, a minimum can be found as specified below by Fig. 11. The algorithm, we elaborated, works in order to find a particular value kopt corresponding to the general minimum such as:
8vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi9
ð27Þ
k2K
where K is the iterations group which length is defined above ðimax jmax Þ, N the envelope length. When this optimum iteration ðkopt Þ is found, the different validation curves can be plotted to check whether actual quantities and estimated ones match each other. 5. Simulation results and discussion Once the state space modelling with inside high enough load is done, the estimation algorithm based on the least squares method is implemented. The data used to perform this algorithm are shown in the Table 1. The machine involved in the modelling presents the following characteristics: Rated phase to phase voltage U n ¼ 530 V, Rated current In ¼ 243 A. The rated impedance of the synchronous machine zn can then be deduced by
Un zn ¼ pffiffiffiffiffiffiffi ð3Þ In
ð28Þ
This last value allows to convert reactances in per units (p.u). 5.1. Validation of machine modelling The implementation of the synchronous machine is done by SimulinkTM/S-function. Before performing a short circuit on the simulated machine, a load test is achieved. The results, corresponding to the steady state, are presented below (see Figs. 12 and 13. The currents and the voltages shown by the figures above, are sinusoidal and well balanced.
686
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
Initialization jmax imax k 0 i=0 xd
i
No
imax ?
end of the process
Yes j=0 xd
i Tdo
i 1 Td
j=j+1 xd
k
j
Tdo
Td ”
k 1 J(k) j
jmax
jmax
Fig. 8. Flow chart for parameters estimation.
envelope (Amperes) j
1500 j
jmax
1 sub transient ranges
i i
1 Transient imax ranges
isteady 2.4 0.4 Time sec Fig. 9. Transient and sub transient ranges subdivision.
5.2. Parameters estimation results As mentioned above, the simulated machine is involve until the steady state is reached. Then a sudden short circuit is applied on its three phases. After this short circuit currents and voltages are recorded and the developed algorithm allows to obtain reactances and time constants given in Table 2: The error we calculated is referred to the accurate parameters values in Table 1 by using the following relation
error ¼
ðaccurate valueÞ ðsimulated valueÞ accurate value
ð29Þ
687
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
6
x 10 4
5
criterion J
4 3
zoom area
2 1 0 0
500 1000 iterations (k)
jmax
1500
Fig. 10. Criterion curve.
2500
criterion J
2000
1500
1000 optimum
500 0
200
400
600 800 1000 1200 1400 iterations (k)
Fig. 11. Zoom on criterion curve.
Table 1 Table of Synchronous machine parameters Parameters
xd (%)
x0d (%)
x00d (%)
xq (%)
x00q (%)
T 0d0 (s)
T 0d (s)
T 00d0 (ms)
T 00d (ms)
T 00q0 (ms)
T 00q (ms)
Accurate values
176
36.1
26.5
150
26.5
3.12
0.64
41
30
41
7
s = second; ms = millisecond.
400
va vb vc
300 Voltage (V)
200 100 0 –100 –200 –300 –400 9.985
9.99
9.995 10 10.005 10.01 Time (Sec.)
Fig. 12. Steady state output voltage in abc frame.
688
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
300
ia ib ic
Current (A)
200 100 0 –100 –200 –300 9.985
9.99
9.995 10 Time(sec.)
10.005
10.01
Fig. 13. Steady state output current in abc frame.
Table 2 Table of synchronous machine estimated parameters Parameters
xd (%)
x0d (%)
x00d (%)
x00q (%)
T 0d0 (s)
T 0d (s)
T 00d0 (ms)
T 00d (ms)
Simulated values Error%
175.7 0.09
34.5 4.1
23 12.9
26.8 1.5
3.28 5.1
0.645 0.8
39 4.9
26 13.3
accurate envelope estimated envelope
2000
ienv (A)
1500
1000
500 0
2
4 6 time (sec.)
8
10
Fig. 14. Actual and simulated current envelopes.
5.3. Discussion on results The Table 2 shows that the estimation of reactances by least squares method provide good results. The precision is such as the error is less than 0.1% in the case of transient reactance. Neverthless, an error of 12.9% is found in the direct subtransient reactance estimation. Even if this last value is high, it is less than those obtained with the traditional method which is widely graphical and which errors sometimes reach 20% or 25% of accurate values, see [17]. As regards the estimated time constants, they match very well with the accurate parameters. Even the sub transient time constant which is difficult to determine with the traditional method, is satisfying. For this last value we got 26 ms instead of 30 ms. Globally, the obtained results are in accordance with what we are supposed to get. Below is the reconstructed envelope from the estimated quantities (i.e., reactances and time constants). In this figure, a comparison with the actual short circuit envelope is done: As we can see on the Fig. 14, there is an accordance between actual signal and the estimated one. This figure reinforces the conclusions we made from Table 2 on the effectiveness of the algorithm we proposed. 6. Conclusion In this paper, a complete modelling of synchronous machine taking in account of the existence of dampers and using only reactances and time constants as parameters is presented. The Park’s framework is used and the modelling is performed with
E. Mouni et al. / Simulation Modelling Practice and Theory 16 (2008) 678–689
689
the state space modelling. The particularity of this one is the use of an inside infinite resistance which allows to create terminals on the model. Therefore a series of external loads can be connected and disconnected without modifying the proper structure of the model. This model is handier than those found in literature which usually include load [7,11]. Indeed, thanks to such a model, tests often performed on the machine to validate it [12] or to identify it [1], become easier to achieve.In the second part of this paper, a statistical technique for determining synchronous machine parameters is proposed. This method is based on only electrical quantities (currents and voltages), so it could be easily achieved by using standard equipments. The method is very efficient and the obtained results suit strongly with those that we were supposed to get. Thus, the method presented in this paper allows not only to determine synchronous machine parameters but also to validate the model built on Matlab/SimulinkTM ). References [1] A. Launay-Querré, J. Gabano, G. Champenois, Parameter estimation of a synchronous generator, Electrimacs, MontrTal, Canada. [2] F. Bruck, F. Himmelstoss, Modelling and simulation of synchronous machine, in: IEEE Computers in Power Electronics, 6th Workshop on 19–22 July, 1998, pp. 81–86. [3] P. Pillay, R. Krishnan, Modeling, simulation, and analysis of permanent-magnet motor drives. Part 1: The permanent-magnet synchronous motor drive, IEEE Transactions on Industry Application 25 (1989). [4] T. Sebastien, G.R. Slemon, M. Rahman, Modelling of permanent synchronous motors, IEEE Transactions on Magnetics Mag 22 (5) (1986) 1069–1071. [5] K. Srivasta, B. Berggren, Simulation of synchronous machine in phase coordinates including magnetic saturation, Elsevier Science Electric Power Systems Research 56 (2000) 177–283. [6] I. Kamwa, P. Viarouge, E. Dickinson, Identification of generalised models of synchronous machines from time domain tests, IEEE Proceedings-c 138 (6) (1991) 485–498. [7] E. Mouni, S. Tnani, G. Champenois, Comparative study of three modelling methods of synchronous generator, in: Conference of the IEEE Industrial Electronics Society, Paris, France, 2006. [8] P.C. Krause, oleg Wasynczuk, S.D. Sudhoff, Analysis of Electric Machinery, IEEE Press, 1994. [9] J. Verbeeck, Standstill frequency response measurement and identification methods for synchronous machines, Thesis, Vrije Unversiteit Brussel, 2000. [10] H. Guesbaoui, C. Iung, Reduced models and characteristic parameters of the synchronous machine obtained by a multi-time scale simplification, in: IEEE Industry Applications Society Annual Meeting, 1994, Conference Record of the 1994 IEEE 2–6 October 1994, vol. 1, pp. 262–271. [11] J.D. Gabano, G. Champenois, Synchronous generator modeling using a non-steady state park model, in: 6th International Symposium on Advanced Electromechanical Motion Systems, Switzerland, 2005. [12] IEEE: Std 1158-1995, IEEE guide: Test procedures for synchronous machines, IEEE Standards, 1995. [13] J.D. Ghali, Détermination des paramétres de la machine synchrone, Ecole Polytechnique Fédérale de Lausanne, 1994. [14] S. Moreau, R. Kahoul, J. Louis, Parameters estimation of permanent magnet synchronous machine without adding extra signal as input excitation, ISIE’04 Ajaccio. [15] F.L. Alvarado, C. Cañizares, Synchronous machine parameters from sudden short tests by back solving, IEEE Transactions on Energy Conversion 4 (2) (1989) 4. June. [16] A. Tumageanian, A. Keyhani, Identification of synchronous machine linear parameters from standstill step voltage input data, IEEE Transactions on Energy Conversion 10 (2) (1995) 10. June. [17] C. Xingang, J. Luming, W. Xusheng, Y. Kaisheng, W. Zhifei, A new approach to determine parameters of synchronous machine using wavelet transform and Prony algorithm, in: 2004 International Conference on Power System Technology – POWERCON 2004 Singapore.