Februari 2011
Influence of damping winding, controller settings and exciter on the damping of rotor angle oscillations in a hydroelectric generator The testing of a mathematical model
Jonathan Hanning
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
electrical feeders on the damping of rotor angle oscillations in a hydroelectric generator
Jonathan Hanning
This thesis has been performed in the university context for Master thesis 30 credits, which is a compulsory exercise in order to gain a degree in electrical engineering.
The thesis main objectives were to investigate how the damping and the stiffness of a hydroelectric generator changed depending on different parameter values, and to test a new mathematical model to calculate the damping and stiffness constants Kd and Ks.
The work has been performed at the request of VG Power, but has been performed at the division for electricity at Uppsala University. The reason for undertaking this thesis was to ensure that generators are robust. But also when building future models for generators, to have a system that can be used to compute robustness.
During this thesis a power cabinet has also been constructed to be able to test the simulated model on a real generator. Under the first five weeks a power cabinet was constructed in the laboratory at the division for electricity. The tests were then performed at a generator with a rated power of 75 kVA.
Handledare: Martin Ranlöf
Detta examensarbete har utförts vid avdelningen för elektricitetslära vid Uppsala universitet.
Uppgiften var att undersöka hur dämpning och styvhet påverkas av olika faktorer i en
generator. En del av arbetet bestod i att jämföra skillnaden mellan kontrollerad dämpning med hjälp av en automatisk spänningsregulator tillsammans med en PSS, mot ett system som använder kopparskenor för dämpning.
Den viktigaste slutsatsen som kan dras av detta examensarbete är att om man vill ange riktlinjer för tillverkare av generatorer, när systemet endast består av en generator och en regulator, bör riktlinjerna bestämmas utifrån massan. Eftersom denna faktor är den viktigaste för robustheten i systemet. Tanken med systemet skulle vara att för varje viktig variabel, så skulle ett värde erhållas och skulle sedan kunna kontrolleras mot en tabell för att säkerställa att inga farliga värden erhålls. Att konstruera denna tabell är ett annat examensarbete, som skulle kräva fler simuleringar på många fler maskiner, och därför bör utföras av någon med en bakgrund inom beräkningsvetenskap.
Den matematiska modellen som testats i detta arbeta behöver lite mer justering på grund av att den inte verkade matcha helt den nuvarande accepterade modellen. Det måste dock sägas, till den nyares försvar, att med vissa inställningar, så korrelerade den mycket bra med den äldre modellen. Men det kommer att behöva ändras och anpassas lite mer, särskilt i beaktande vid beräkningen av den synkrona vridmomentskoefficienten, som nästan alltid verkade vara 10 till 30 procent för låg.
D-Q-axis Direct and Quadrature axis DAE Differential-Algebraic Equation Et Terminal Voltage
H An inertia constant
Ka/Kp Gain constant in the feedback system
Kd (1) Damping constant in electric torque equation Kd (2) Derivative constant in the feedback system Ki Integrating constant in the feedback system Ks Synchronous constant in electric torque equation ODE Ordinary Differential Equation
Pf Power Factor
PhD “Philosophiæ Doctor” or Doctor of Philosophy PSS Power System Stabilizer
P.U. Per Unit
Re Resistance in the tie-line SMIB Single Machine, infinite bus
St Power output
Td Foresight of the time step Te Electrical torque
Xe Reactance in the tie-line UU Uppsala University
task was to conduct research about the damping and stiffness of a generator. One part was to compare the difference with controlled damping with the help of an automatic voltage regulator, together with a power system stabilizer. And also a system which used copper bars for damping.
The main conclusion that can be drawn from this thesis is that if you want to provide guidelines for manufacturers of generators, when the system contains only of the generator and a regulator, the guidelines should be determined by the mass. Since this factor is the most important one for the robustness of the system. The idea of the system would be that for each important variable, a number is acquired and could then be checked against a table to ensure that no dangerous values are obtained. To construct a table like this is another thesis, which would need a lot more machines to be simulated on, and therefore should be performed by someone with a background in scientific computing.
The mathematical model tested in this thesis need some more adjusting, due to the fact that it did not seem to match entirely to the current accepted model. It must be said, though the latter’s defense, that with some settings, the altered mathematical model matched very well.
But it will need to be modified and tuned some more, especially in regard to the calculation of the synchronous torque coefficient, which almost always seemed to be 10 to 30 percent to low.
Foreword
This thesis has been performed in the university context for Master thesis 30 credits, which is a compulsory exercise in order to gain a degree in electrical engineering.
The thesis was to investigate how the damping and the stiffness of a hydroelectric generator changed depending on different parameter values. And also to test an altered mathematical model to calculate the damping and stiffness constants Kd and Ks. The work has been performed at the division for electricity at Uppsala University, as a joint operation together with VG Power. The reason for undertaking this thesis was to ensure that generators are robust. But also when building future models for generators, to have a system that can be used to compute robustness.
I would like to especially thank my supervisor Martin Ranlöf for the time he has spent helping me over the threshold incurred during my work. My thanks are also directed to the people in the same working group, Johan Lidenholm, whose thesis has been very helpful, but who has also helped with understanding some problems in Matlab. Thanks also to Mattias Wallin for much practical instruction during the construction of the power cabinet. Also big thanks to Urban Lundin for the hydropower course and for making this thesis possible and I would also like to thank Kjartan Halvorsen for the help with the automatic control. Also thanks to Stefan Pålsson for this knowledge with Matlab. Last but absolutely not least, all the teachers who has put in much effort in my education so that the courses I have read has become much more interesting, thanks also to my examiner, Nora Masszi.
Jonathan Hanning January 2011 Uppsala
1 Introduction ... 4
1.1 Background ... 4
1.2 Method ... 5
1.3 Demarcation ... 5
1.4 Objectives ... 5
2 Theory ... 7
2.1 Synchronous generator ... 7
2.2 The damping and synchronous coefficient ... 7
2.3 System analysis ... 8
2.4 The mathematical model ... 9
2.4.1 Direct and Quadrature axis ... 9
2.4.2 Per unit representation ... 10
2.4.3 State-space representation ... 10
2.4.4 Automatic control ... 10
2.4.5 The nine basic equations ... 11
2.4.6 Ordinary differential equations solver ... 12
2.4.7 Standard parameters ... 12
2.5 Rotor angle oscillation ... 12
3 Method and construction ... 14
3.1 Mathematical model in matlab for simulation ... 14
3.1.1 The introducing of state space representation ... 14
3.2 The construction of the power cabinet ... 16
3.2.1 Modifying the generator with damper bars ... 16
3.3 Operating the generator ... 17
4 Results ... 18
4.1 Simulation results ... 18
4.1.1 Single machine without regulators ... 18
4.1.2 Single machine with an automatic voltage P-regulator ... 19
4.1.3 Single machine with an automatic voltage PD-regulator ... 20
4.1.4 Single machine with an automatic voltage PID-regulator ... 21
4.1.5 Single machine with both PID-regulator and PSS ... 22
4.1.6 The mathematical model ... 22
4.1.7 Unstable systems ... 23
4.2 Laboratory tests ... 23
4.2.1 Connecting the generator to the grid ... 23
5 Discussion ... 25
5.1 The simulation model ... 25
5.2 The constructed power cabinet ... 25
5.3 The results ... 25
5.4 Future work ... 26
5.5 Confounding ... 26
6 References ... 27
6.1 Literature ... 27
7 Appendix ... 28
1 Introduction
A different model, [1] compared with the accepted model, has been expanded and tested to calculate the synchronizing and damping components of electrical torque developed in a synchronous machine. The method is based on the numerical analysis of system response time, using least squares adjustment.
1.1 Background
Since the introduction of synchronous generators in the late 19th century, the way of operating a system with several generators has significantly improved. In the early years, it was not unusual to have power black-outs over a huge area of the grid. But when the regulation was modernized, it has become more and more unusual with power failure. Nowadays it is almost required a storm which destroys a cable to receive a power failure.
The stability of a power grid is depending both on the total grid, but also on its individual components. Usually in a grid there are power consumers, power producers, power transmission and power control. And since the producers are depending on the consumers, there has always been an interest of how the producing unit reacts to changes in consuming.
For example how the electrical torque changes when a huge load is connected to the grid. The electrical torque is built up by the synchronous and damping constants of the generator.
Therefore these constants have been of interest for some time.
The ability to calculate the damping and synchronizing constants has been an important problem since the expansion of power system interconnections. And since the improvement of digital computers and modern control theory, a better control of power systems has been gained. However, the method how to calculate these torque components has not improved at the same rate. This new approach is thus based on the time-domain analysis of system response. Precision depends on how good the accuracy was of the time response.
Due to imperfection in the system, a couple of oscillations will occur. The most interesting and important one is the rotor angle oscillation. This oscillation occurs when the power is raised or lowered, and the generator is trying to find its new equilibrium, the equilibrium between the torque from the turbine and the electrical torque. This oscillation gives a few other oscillations, which will be studied in this thesis. For example the oscillation in the power produced. The power produced is connected to the swing equation, equation 1b, which is connecting the rotor angle acceleration, the mechanical torque, and the electrical torque.
This is further described in chapter 2.2.
When measuring is performed on a generator, different variables are calculated in order to be able to compare the results. It is usual to use either damping time constant Td, which is the time required for the amplitude to decrease to a new value from its original value. Another value that is often use is the damping torque coefficient, Kd, which is used to identify the rotor angle stability of the system. Yet another constant that is interesting is the damping capacity b, which is the ability to absorb vibration by internal friction. In the current situation, there is
The purpose of this thesis is to show how the damping and synchronous constant is affected by the electrical design, selection of exciters, adding damping parameters and inertia. The equation relationship between attenuation and these factors are certainly well known, but are found in various places in the literature. The study of damping has been performed on a two- axle model of a generator.
1.2 Method
This thesis will analyze this problem in two different ways; first the expanded mathematic model which will be tested in Matlab. In this thesis are also included some testing on a real generator. To make this possible, a synchronization unit will be built to link up one of the division of electricity´s generator to the electrical grid. Since Uppsala does not have any great waterfall, the generator is driven by a motor which is connected to the rotor, instead of a turbine. Therefore there will be the possibility to perform quick torque changes. The idea is then that the natural damping of the generator, which occurs due to the copper in the rotor windings, together with the quick torque change, will give rise to oscillating revolutions per minute. This oscillation should continue until the generator has found its new steady-state level of operation. This should also provide an oscillating power curve. Another test will be to connect a network of copper bars to get a stronger damping, due to the current that will be induced in these which will counteract the change in torque.
1.3 Demarcation
The idea of this thesis is to develop a functional program for the mathematic model in Matlab, which can simulate different machines, with different parameter values. This program will then be modified so that you can connect an automatic voltage regulator in front of the generator, and the final version should also include a power system stabilizer. This is done to be able to compare the stability and robustness of a system, due to variation in settings on the control systems and various types of generators. A proposal will also be included of how the value of the included parameters in the regulator should be, to ensure stable systems. If the results of the simulations in Matlab give a distinct and unambiguous picture, a proposal for recommendation of generator design in terms of robustness and torque stability will be made.
Primary focus will be to investigate the influence of the automatic voltage regulator’s parameters on stability of the system, while the power system stabilizer will more or less, if not time permits, just be implemented. The simulations will only be on a single machine infinite bus. No simulation will be tested on island grids (weak bus). Primary focus will be on the synchronous coefficient and the damping coefficient, described more closely in chapter 2.2. Secondary focus will be on the changes in rotor angle velocity, and torque change and also how the electric angle between rotor and stator magnetic axes difference. Most units in this program will be measured in per unit, and the reason is that it gives more comparable data.
1.4 Objectives
The primary objective is to construct a program for analysis of the damping and synchronous coefficient. Secondary is to build a functional synchronizing unit to be able to connect a generator to the grid. This to make it possible to try the theory in reality, but also for further experiments being performed by other students and PhD under the division of electricity.
After completing the main objective, a series of add-ons is desirable. For example, the possibility to use an automatic voltage regulator and a power system stabilizer to increase the robustness and stability of the system. Another sub objectives there would be to investigate how the different parameters change the system. Another sub objective is to analyze how the new mathematical model works, compared to the old accepted model, with other words, if they correlate. Another objective is to try to make up a system so companies who design generators can have some kind of model for robustness and stability when designing the generators. It would also be desirable to look into the stability from an automatic control perspective, due to the fact that both the automatic voltage regulators as well as the power system stabilizers are feedback systems. A desirable and maybe final objective would be if the simulated results would correlate with those which can be measured in the laboratory, primary the change in rotor speed when a disturbance in torque is being done. This thesis is being done because there is a gap in the literature concerning how the stability is inflicted by an automatic voltage regulator and its parameters.
2 Theory
In this chapter we will familiarize ourselves with the theory underlying this thesis, if one wants to look further into this theory, reference [2] is recommended.
2.1 Synchronous generator
A large part of the power production in the world originates from power stations using generators directly connected to the grid, synchronous generators. A synchronous generator rotates with a speed that is proportional to the frequency of the current in the armature. The magnetic field that is created by the armature currents, rotates with the same speed as that created by the current on the rotor, the field current. If it is a strong grid, special preparations have to take place to meet the demands. In general there are five conditions that are required before synchronizing a generator to the grid, but phase sequence and waveform should be fixed by the construction of the generator and its connection to the system. But voltage, frequency and phase angle must be controlled each time a generator is to be connected to a grid.
1. The generator frequency is equal to the system frequency.
ω1 = ω2.
2. The generator voltage is equal to the system voltage.
E= V (Generator E = Grid V)
3. The generator voltage is in phase with the system voltage.
α = 0 (phase difference)
A voltage difference will result in a steady flow of reactive power and when this coincident with a frequency difference a substantial reactive and active power will flow back and forth to the grid under a short interval of time that could damage the generator.
2.2 The damping and synchronous coefficient
The equation that is the focus of this thesis is described below, equation 1a, chapter 2.2 in reference [2]. This describes the changes in the electrical torque ▲Te, depending on the synchronous coefficient, Ks, which is multiplied with the change in the electrical rotor angle
▲δ. Then it is the damping coefficient, Kd, which is multiplied with the change in rotor angle velocity ▲ω. The change in electrical torque is taken from its context where it usual belongs, the swing equation, equation 1b. Where J is the total moment of inertia of the rotor mass, Tm is the mechanical torque supplied by the prime mover. Te is the electrical torque output of the alternator, and θ is the angular position of the rotor in radians.
[Equation 1a]
[Equation 1b]
It is around this equation the thesis is built. But also about how to calculate Ks and Kd. The procedure to calculate the Ks and Kd works by using the time response of the torque, speed and angle. And then by applying the least squares adjustment to obtain the electric torque to the two signals [3], equation 1. Then the error can be determined by equation 2:
[Equation 2]
And to be able to see the summation of error over time, we are forced to integrate over the interval of oscillations, equation 3.
[Equation 3]
Damping and synchronizing torque coefficients Kd and Ks are calculated to minimize the integral of the least squares adjustment. They must fulfill these two equations, equation 4 and equation 5.
[Equation 4 (top)]
[Equation 5 (below)]
Kd and Ks are time-independent, hence also the differential equations and integration parameters are time-independent. This means that it is possible to change the differential equation and integrating the system. This makes it possible to rewrite equation 4 and equation 5, to equation 6 and equation 7:
[Equation 6]
[Equation 7]
By using the above equations, Kd and Ks could be estimated, the other values in the equations are calculated numerically in the simulation, which can be further studied in chapter 2.4.
2.3 System analysis
This thesis has been tested on several known systems. Which was given an electric or a torque disturbance, and then with help from equations 6 and 7, with the stated time integrations performed numerically, on the supplied data from a simulation and the resulting algebraic equations solved for Ks and Kd. The Equations can then be used with the three time responses to calculate the damping and synchronizing constants. The different settings on the machines can be found in appendix 1. The system was tested on a single machine, infinite bus.
2.4 The mathematical model
One of the main objectives of this thesis was to try an altered mathematical model, with the help of Matlab. The basic mathematical model was constructed in such a way that nine equations, see chapter 2.4.5, were needed to be met, in order to calculate Ks and Kd, both with the analytic and square alignment. Much more information can be found in appendix 2, which is the code for the Matlab-model. Here the mathematic equations are put into its context. Which may simplify the understanding on how they are used, therefore, they will be less described here in the text. These nine basic equations are divided into two different kinds of equations. The first five equations will return an actual value, with help from Matlab. This value will be obtained with help of numerical analysis. The other four equations are calculated to become zero. The reason why nine equations are needed is simply due to the fact of the numbers of unknown variables, later on in this thesis, more equations will be added to satisfy the need when new variables arose as a result of more complicated and automation systems.
2.4.1 Direct and Quadrature axis
Direct and quadrature axis is more known as d-q-axis. They are often used when calculating on a generator, instead of traditional x-y-axis. This due to the simplified equations received, when looking at an electrical point of view. They are simplified because the magnetic circuits and all rotor windings are symmetrical with respect to both polar axis and the inter-polar axis.
The direct (d) axis is defined by that it is centered magnetically in the centre of the north pole, and the quadrature (q) axis is defined by that it is 90 (electrical) degrees ahead of the d-axis.
Further on, the position of the rotor, relative to the stator, is measured by the angle theta, which is the angle between the d-axis and the magnetic axis of phase a winding, seen from above. The choice of the q-axis leading the d-axis is based on the IEEE standard definition [4]. Throughout this thesis the d-q-axis is used if nothing else is specified, for a graphical representation, see figure 1.
[Figure 1: shows the relationship between the rotor and stator-current together with the d-q- axis]
2.4.2 Per unit representation
It is very common in power systems calculations to use per unit system, due to the fact that it becomes much easier to compare results when you do not need to remember real values, the per unit system is functioning according to equation 8a. In this thesis, per unit system is standard for all units, such as voltage, current, resistance, and so on. One example is that since the simulations are performed on different machines, it would be hard to compare the results if you constantly would need to check-up the basic value of each machine, it is easier now to see the percentage difference instead and compare that way. So for each quantity that is used in the simulation, a base value is chosen, then use your actual quantity value and divide with your base value, to receive your quantity in the per unit system. An example is given in equation 8b, which will represent Machine I in appendix 1. This example will show the rated active power.
[Equation 8a]
[Equation 8b]
2.4.3 State-space representation
As mentioned earlier, in chapter 2.3, the need for more equations due to more variables will be discussed here. The way of solving the problem with more complex regulation system, machine 5 in Appendix 1, is to add more equations to the system. First, the state-space representation should be introduced, chapter 12.2.6 in reference [2], equation 9, which shows a state space representation for a time-domain solution. A, B, C and D are just constants, but with higher order systems, they will become rows and columns with constants. These constants can represent different matrices. In this thesis, A will be the state matrix. B is the input matrix. C is the output matrix, and D is the feedthrough, or feedforward, matrix. If the system only contains an automatic voltage regulator, then x is the automatic voltage regulator state-vector, and u is the difference between the wanted signal and the feedback connected signal.
[Equation 9]
This system is a time-domain representation of the transfer function, described further in chapter 3.1.1, obtained through inversed Laplace-transformation. The derivative dx/dt will be a column by changing constants when the system becomes more complex, these columns will provide the need for new equations as new variables appear.
2.4.4 Automatic control
In order not to totally rely on the strong grid, in maintaining the frequency, phase and voltage within permissible limits, various devices may be connected in front of the generator to help
can be of different complexity, from as simple as just a gain, that means that the error is amplified by a factor k, to quickly eliminate the difference between desired and actual value.
But the regulator might be of the degree to have gain, integrating step, and derivation step, together with limits on the signals, to ensure no transients disturbs the system. More information about the different regulators and power stabilizers can be found in appendix 3.
2.4.5 The nine basic equations
The number of equations that were needed was determined by the number of unknown variables, and without a regulator, there were nine unknown. That led to the need for nine equations that could determine nine unknowns. The equations represent state variables, which will receive a new value for every discrete time step that the simulation takes. Since many of the variables in equations 10-18 contains other variables which also need new values for every discrete time step, there will be more lines of code executed than these nine, see appendix 2 for more information. All equations are represented in the per unit scale, and also in d-q-axis. The equations were derived, chapter 13.3 in reference [2], and here they are presented in their final versions, equation 10-18, as they were used in the code. Equation 10 and 11 is part of the swing equation, which describes how the rotor speed is inflicted by the unbalance of the mechanical torque and the braking torque. Equation 10 shows how the rotor speed changes, and to get the right unit, it is multiplied with the base for angular velocity.
Equation 11 shows what will happen to the mass equation when the torque is changed, H is an inertia constant (described in equation 20), and Pm0 is the mechanical torque, and Pe is the electrical torque. For equation 12 to 18, mentioned should be that Psi stands for flux linkage, d and q is the axis, L is the inductance, R stands for resistance. Efd is the measured voltage going out from the regulator, and i stand for current. Ed and Eq is the terminal voltage of the d and q axis, while Ebd and Ebq is the d and the q components for the field voltage Eb. X represents the reactance in the tie line.
[Equation 10]
[Equation 11]
[Equation 12]
[Equation 13]
[Equation 14]
[Equation 15]
[Equation 16]
[Equation 17]
[Equation 18]
2.4.6 Ordinary differential equations solver
In order to function with sufficient precision in the calculations an ode-solver is used, which produces new values for each discrete time step the simulation takes. It is for the nine equations described in chapter 2.3.5, that the ODE-solver obtains new values in every time step. But since the variables in the nine equations difference as well, there will be more than nine lines of code which will be executed on every discrete time step. The ODE-solver used in these simulations were Matlabs ODE23t, since it is quite fast but still with high accuracy. It is also good to use when you got differential algebraic equations, which the last four equations of the nine basic ones are.
2.4.7 Standard parameters
Since this thesis has been engaged by a company, standard parameters may differ a bit from the manufactory to the university. The parameters you insert in this program are the standard parameters for companies, which are based on inductances. In the chosen Per Unit system, the inductances are equal to the corresponding reactances. To be able to present numbers which both can be satisfied with, a transfer-script was needed, to see these scripts, appendix 2 should be studied, StandardParam.m and StandardParamTieLine.m. An example can be found below, in example 1. Describing the saturated synchronous q-axis reactance, this is the standard parameter, while Laq and L1 are parameters taken from the generator manufacturer. These are the inductance in the phase a, on the q-axis, Laq and the leakage inductance, Ll.s
Xq = Laq + Ll;
[Example 1]
2.5 Rotor angle oscillation
In electric power engineering, there is a concept of rotor angle oscillation of a generating unit, consisting of generator, shaft and turbine, which slowly oscillate around the synchronous
taken in order to remove the fluctuations. Among the measures that are most common there is damper windings which provides natural attenuation, which add damping contribution from the current that is induced. The damping of rotor angle oscillations depends on the generator’s electrical design, with the reactance’s and time constants, inertias, the operating point,
controller settings and more.
Due to that synchronous generators are directly connected to the grid, change of load on the grid directly influence the frequency. For example, if a huge load connects to the grid, the frequency tends to lower a little bit from its stable operation point, until more power is produced. Either by letting more water through or connecting another hydropower station onto the grid. In this thesis the thought is to give the turbine, which actually is a motor, order to raise its revolutions per minute. And then hopefully get a rotor angle oscillation, which can be detected with the help of measuring equipment. The change in rotor angle should give rise to a few other oscillations, for example in power production, which is based on voltage and current.
3 Method and construction
In this chapter we will familiarize ourselves with the method and construction that has been performed during this thesis, this chapter is divided into two parts, the first part will address the theoretical part. While the second part will explain how the practical testing and construction of the power cabinet were built up.
3.1 Mathematical model in matlab for simulation
This part of the thesis is a continued study of Johan Lidenholm’s thesis [5], which included a Matlab-program, which several parts of this program has been aided by. The program sets up a generator with different parameters, which is connected to an infinite bus. To be able to compare results more easily, all parameters are converted to the per unit system. The generator is then running at nominal speed and nominal voltage level for a time t1. After that, either an electrical error occur, or a torque disturbance, this for a short period t2. At last, the voltage or torque disturbance is restored, and the generator is trying to get back to normal state of operation, during the time t3. And it is the behavior during the time t3 that is studied, and the behavior can be altered by changing parameters, or connect different kind of automatic control units in front of the generator. The interesting data is then saved from the simulation, and then treated in order to compare between different simulations to obtain a behavior that correlates with the changes in settings or parameters. Even without any automatic controller unit, the generator is supposed to have a moment of inertia in itself because of the mass, and the damping also gets a contribution from the resistance is the windings. Also the strong grid is supposed to try to get the generator back to normal steady state. Therefore a number of simulations were made on just plain generators connected to an infinite bus, to see how much the generators parameters would change the stability of the system. And then especially of interest would be the change in the electrical torque.
3.1.1 The introducing of state space representation
In order to use automatic voltage regulator, the feedback system needs to be reversed laplace- transformed back to the time-domain. Since its representation is in the frequency-domain. The complete systems can be found in appendix 3, but below in figure 2, one system is represented.
[Figure 2: shows the second automatic voltage used in this thesis, Kp is the gain-factor, Ki/s is the integrating part, and Kd*s/(1+sTd) is the derivative part]
As can be seen in equation 9, y is what comes out of the system, efd in figure 2. The equation 9 is a classical way of representing a state space equation. Therefore there is a function in Matlab which could be used. It works in the manner of that you insert the transfer function, from figure 2, which can be found in equation 19, and the function in matlab gives in return values for A, B, C and D. The system whose transfer function can be seen in equation 19 is a regulator with, gain, derivative and integrating part.
[Equation 19]
After retrieving these constants, they are used in the differential-algebraic equations. Used for calculating a new value on dx/dt, appendix 2, in simulation3.m. Since a new value must be obtained in every discrete time step, dx/dt is one of the time-dependent variables that the ordinary differential equations solver produces for each time step. But as for all automatic control units, the system is always using the last value, to control its next behavior. So the need for short discrete time step is inevitable. The change in e, is actually Vref (reference) minus the actual voltage value et. This is the furthest left summation in figure 2.
3.2 The construction of the power cabinet
In order to be able to try out the theory of this thesis, as well as trying to see if the results from the simulation could be seen in reality, a power cabinet had to be constructed in order to be able to connect one of the division of electricity’s generators onto the Swedish grid. This then could be seen as a single machine infinite bus. The construction work took an estimated three weeks of this thesis. All parts were ordered and the work of assembling them was done in the laboratory hall. The final result of the construction can be seen in appendix 4. The power cabinet contained a synchronization unit, and the interior can be studied below in figure 3a (to the left) and the exterior is shown in figure 3b (to the right).
[Figure 3a (to the left) shows the inside of the power cabinet, which housed the
synchronization unit, at the top one can see the cables going out to the connection to the grid, and at the bottom, is the cables that is connecting with the generator]
[Figure 3b (to the right) shows the outside of the power cabinet, with its voltmeter, ampere meter and inductive meter. The buttons is for connecting and disconnecting the unit]
3.2.1 Modifying the generator with damper bars
As one test out of many in this thesis, a series of damper bars were connected to the generator, these damper bars were made out of copper, and were connected to each others, by bridges between the bars. The reason for this modification was the idea of that the copper bars would increase the generators damping, due to the fact that the copper would induce a current which
3.3 Operating the generator
Once the generator was connected to the grid, a series of testes with different characters were performed. These tests were to try to show the oscillation in rotor angle, by changing the amount of power inserted to the motor. In normal case, if it would not have been connected to the grid, it would have raised its revolutions per minute. But now instead it would raise its torque. But before it reached its new steady state, an oscillation should be possible to measure. The measuring of revolutions per minute was done by using laser equipment and different color stripes on the rotor.
First of all, the generator had to be connected to the grid, this was done by using the earlier described synchronization unit. The speed of the motor, which acted as a turbine, was set a little bit higher then which were required to gain the 50 hertz. After that, the speed were set to a value, which would represent lower than 50 hertz, and hopefully in the transition in-between these values, a connection to the grid could be made. Otherwise one had to do the procedure again, the other way around, first to low speed, and then raise the power to the motor, and hope that the three requirements, described in chapter 2.1, would be met.
4 Results
This part of the thesis will not be presented in its entirety, since all the data that was obtained during the work is impossible to present. Therefore, just the part that was considered important to present is represented here, the rest of the data will be available in the appendix or in some sort of link attachment.
4.1 Simulation results
In this chapter the results of the simulation will be presented, due to the amount of data received during these simulations, most of the figures will be in appendix 6-12. The values chosen for these simulations were based upon recommendation, chapter 12.4 in reference [2]
and from manufacturers [6], but sometimes interesting results were followed up by simulations with values outside those boundaries.
4.1.1 Single machine without regulators
The first part of the thesis was to investigate different type of machines which should be considered to be equal real generators. The simulations were performed by first running with the standard settings, and then the parameters were changed. Just to try to get a grip of how much influence the different parameters had for the overall performance in concern of stability. The original settings for each machine can be found in appendix 1. Some parameters were changed to see the impact of stability, the parameters that were changed were:
Et: This is the terminal voltage
H: This is an inertia constant, which is based on equation 20. And it shows how H, which is used in equation 11, is depending on J, which is mass moment of inertia, wm , which is the rated mechanical angular frequency and S which is the apparent power in VA
Pf: This is the power factor
Re: This is the resistance in the tie-line Xe: This is the reactance in the tie-line St: This is the apparent Power output
[Equation 20]
The result of this simulation can be found in appendix 6-9, together with the results for the other machines. Conclusions that can be drawn from these simulations are that the most important factor for stability and robustness, which is shown in figure 4, of the system was when the parameter H was changed. This was done in such a manner that it was set to a value, and hence the equation 20 was overwritten, the values was in the per unit system. Important things to note from Figure 4 is that at lower H, the system becomes sensitive to disturbance, but also more rapid to return to stable operation after the disturbance. One can generally say that when the value of H is higher, the system is simply slower. And since H is dependent on the mass of the rotor, the mass should be the factor that sets the guidelines of other parameters.
[Figure 4: shows how the torque reacts when a disturbance is made, at time t1 = 1 sec, the red color is H=2, green 3, blue 4, magenta 5, yellow 6]
If mass is a parameter that could be changed during construction, it should be adapted in a way that if one want a robust system, the mass should be maximized, and if one want a system that can response to quick oscillations, the mass should be reduced.
4.1.2 Single machine with an automatic voltage P-regulator
This simulation was made with a P-regulator, which is a simple regulator with only a gain step which increases the error between the wanted signal and the actual signal. It should be said that this simulation were made without regards for stability criteria for the feedback system. A stable system is a system which has only negative poles, the poles are obtained by the solutions to the denominator roots. That means that some of the results might be on unstable systems, which will give some strange results, as can be seen in figure 5, Kd (the damping constant) is negative for higher value of Ka (gain constant in the feedback system).
This unstable feedback system may be the reason why suddenly the value 1000 on Ka seems to increase the damping constant, all data is consolidated in appendix 10. The machine used with the P-regulator is machine V, which can be closely studied in appendix 1, the main feature of machine V is that R1d and R1q is set to 1 p.u. which should represent a machine without damping.
[Figure 5: shows how the damping constant change with increasing Ka (gain value in the feedback system), negative value of Kd might represent an unstable feedback system]
4.1.3 Single machine with an automatic voltage PD-regulator
This simulation was made with a PD-regulator, which both has a gain and a derivative step.
The machine used with the PD-regulator is machine V, which can be closely studied in appendix 1, the main feature of machine V is that R1d and R1q is set to 1, which should represent a machine without damping. Except changing the gain and derivative constants, also the time constant Td was altered to see the impact on the synchronizing and damping constant. Td is the foresight of the time step. One thing that should be kept in mind when reading through the results, which can be found in appendix 11, is that some of the combinations, of the values of gain, the derivative and time step constants, might create unstable feedback systems. This might be the reason for why in figure 6, it is not entirely conclusive, but one should still be able to see the trends. One interesting and maybe alarming trend is that the result seems to vary very little with a low value on Ka (gain constant), blue line below. Every value which calculate Ks has been analyzed, and they all seems to be in the per unit system. But even so, a larger impact could be expected when you multiply the error between the wanted signal and the actual signal. But further analyzes are needed to distinguish if any errors has been done, and to investigate if the feedback systems are stable.
[Figure 6: shows how Ks change with Ka (gain) and Kd (derivative, not damping constant)]
4.1.4 Single machine with an automatic voltage PID-regulator
This simulation was made with a PID-regulator, which has a gain, a integrating and a derivative step. The machine used with the PID-regulator is machine V, which can be closely studied in appendix 1, the main feature of machine V is that R1d and R1q is set to 1, which should represent a machine without damping. Except changing the gain, integrating and derivative constants, also the time constant Td was altered. The possibilities to alter different settings made the amount of data enormous, for intense studies of specific cases, see appendix 12, below in figure 7, the synchronous constant can be studied, with different Kd and Ki.
[Figure 7: shows how Ks change with Kd (derivative part of regulator, not damping constant) and Ki (integrating part of regulator)]
The trend seems to say that with increasing integrating constant, the synchronous constant is lowered. The strange thing is that it does not seem to be linear, due to the fact of the order of the lines, which can be seen in figure 7. One could think that they should be in rising or in dropping order, but it seems more randomly then that, maybe due to unstable feedback systems.
4.1.5 Single machine with both PID-regulator and PSS
The results from this study should not be presented, since either the mathematical model or the computer power was not enough to calculate this model with accurate accuracy. Due to the many T-values (foresight in time step) in the PSS, the model in Matlab could never finish, even though it was given a 24 hour time span. Even so, the model is probably right, but need more computer power or other limits in the mathematical model, the PSS-code can be found in appendix 2 (DAE).
4.1.6 The mathematical model
The new mathematical model [1] was proven to work quite well, even thou it did not entirely match, it seems that it still need some folding to be a perfect match, as can be seen in figure 8, it is close to a perfect match. The black colored lines are the ones that is the new mathematical model, and the different colors are the old model. For a closer look how well they matched, appendix 13 is recommended which is a table of data from study I, machine I.
[Figure 8: shows how Ks and Kd changes when a disturbance is made, at time t1 = 0 sec, the red color is Xe 0.0 green 0.15 blue 0.30, magenta 0.45, yellow 0.6, the black is the new mathematic model]
4.1.7 Unstable systems
During some simulations, unexpected results occurred. And many of those might be due to unstable feedback systems. It is not always easy to see from the response on Ks and Kd, but it can be seen from the other graphs, see figure 9 for an example of an unstable feedback
system. And as can be easily spotted, when the integrating part of the feedback system is given a very high value, the whole system gets unstable.
[Figure 9: shows how the angle reacts when a disturbance is made, at time t1 = 1 sec, Ka is 100 and the red color is Ki 0.005, green 0.1 blue 0.4, magenta 1.0, yellow 2.0]
4.2 Laboratory tests
Due to the resources at the division of electricity at Uppsala University, this thesis could be tested with a generator. The test equipment is described in chapter 3.2. And pictures can be found in appendix 4 and 5. A set-up with different tests was performed with this equipment.
To try and establish if the simulation result could be transposed into the laboratory generator.
The generator used in these tests was a synchronous generator with rated power 75kVA.
4.2.1 Connecting the generator to the grid
The first test that were performed and recorded were the connection of the generator to the grid, to be able to perform this, the right voltage, the right frequency and the right phase was needed to be obtained. Instead of a river through the laboratory, a motor acted as the turbine, this motor could simulate different kind of disturbance. Mostly used was to try to change the revolutions per minute, but due to the connection to the strong grid, it was not possible to raise the RPM. Instead a power increase occurred, sadly due to the strong net, it was hard to
see the mechanical oscillations. They were so small that they were lost is the accuracy of the measurement. But the oscillations in power could be perfectly seen, see figure 10.
[Figure 10: Shows when a torque change is performed on the motor what happened to the total power output. The oscillations can be perfectly seen directly after the change, before stabilizing on a level, two changes are made in this figure]
Time (second) Power
(MVA)
5 Discussion
This thesis has been testing a new mathematical model for calculating the electric torque in a generator, with the help from the synchronous constant Ks and the damping constant Kd.
5.1 The simulation model
This thesis has provided a method in a program for using advanced regulators and PSS together with a generator. The PSS could not be implemented in a satisfying way, due to lack of computer power or more likely, the input to the ordinary differential equation solver in matlab is in such a way that the condition never could be meet. It was tested to let the program run for 24 hours, and still it had not produced a single value. That is if you think of the ten seconds it was supposed to run as a long chain of discreet data points. So the model probably needs new conditions that are put into the solver. Therefore it is recommended that someone with a background in scientific computing continue forward and investigates if the feedback system is stable, and also continue with the mathematical model and investigate if it needs folding, since it does not entirely match the common used model.
Also the results should perhaps be more consolidated with someone with scientific computing, since my own knowledge in the subject is limited. The program which was programmed to solve this task should be streamlined, since there is a need for many changes and clicks to produce one result. Perhaps even the layout of the program should be changed in order to make it more easily used by a third part person.
5.2 The constructed power cabinet
Except for some minor misses and flaws during construction, most of the work went without problems. The one thing I want to recommend to someone making their first power cabinet is to drag cable corridors. Even though it seems to work out in the beginning, there will be a lot of cables in the end. As do not be stingy when it comes to the use of cable, or you will regret it later. Now afterwards, I know how I should have constructed my first power cabinet.
As far as to this date, its functionality has worked flawlessly, the process of connecting it to the grid I would recommend is a two person job, since it is hard to both keep the control voltage level and the right revolutions per minute at the same time as trying to connect the generator to the grid.
5.3 The results
The results show clearly some trends which could be expected when tuning on the regulator.
For example, you get a faster system with a higher gain value on the regulator, but at the same time higher transients. But one should be alarmed when the tuning constants on the regulator are getting high values, since the probability of an unstable system is imminent. None the less this thesis could work as guidelines, both for manufacturers and for further research into this unexplored territory. Most of the results is quite expected, for example when the gain is raised, you get a faster system, but with higher transients and oscillations. A trivial but confusing matter is that both the damping constant and the derivative part of the feedback system has the abbreviation Kd, and maybe one should be altered, but they are both very standard to use. When using derivation and integration steps, one could think that the transients should get to a lower altitude. And sometimes that was the case, but far from
always, which might again be due to an unstable system. The reason for having integrating and derivation step is to try to optimize the return to steady state after a disturbance. And sometimes it helped with these steps, but not always, and for high values on these constants, Kd and Ki, it even had the reverse effect. The results concerning the automatic voltage regulator with a derivative step came up with some problems which also were encountered in reference [5]. Which was that sometimes it seemed like the system swung twice around a center point, see Appendix 11, full version for the figures.
5.4 Future work
I would recommend letting someone with a scientific computing background, maybe as a master thesis, look over the model and run some simulations with confirmed stable systems.
In order to make the results interesting for manufacturers of generators, the compared system all should be stable and tested on real models of generators. Also the model should be altered so that the PSS could be implemented, since it is an important part of the total system. Also it would be interesting to compare simulation on the same machine with and without different kind of regulators. To investigate how important they are, and how the generator changes its response due to different disturbances.
5.5 Sources of errors
Mainly I would suspect that many of the simulation had parameter values which gave an unstable system, which would produce results which is not trustworthy. Another thing that could be worth investigated is what limits the time step should be in, because with a value to big, the regulator had no influence. Probably because the changes are faster than the regulator could handle. The question is then what value is small enough to get a reliable result, because when you lower the time step, the amount of data for every simulation grows. And due to that, the simulation with a PSS could not be performed. For those tests performed with a real generator, the errors should strictly be bound to the inaccuracy of the measurement equipment.
6 References
Below are the references which have been used in this thesis.
6.1 Literature
[1] R.T.H. Alden and A.A. Shaltout, ”Analysis of Damping and Synchronizing Torques Part I – A General Calculation Method”, IEEE Transactions on Power Apparatus and Systems, vol.
PAS-98, Sept./Oct. 1979.
[2] P. Kundur “Power System Stability and Control”, McGraw-Hill inc. 1994.
[3] B.P. Lathi, “An introduction to Random Signals and Communication Theory”, International Textbook Company, Chapter 1, 1968.
[4] ANSI/IEEE Standard 100-1977, “IEEE Standard Dictionary of Electrical and Electronic Terms”.
[5] J Lidenholm “Power System Stabilizer Performance” ISSN 1401-5757, UPTEC F07 109 [6] M. Wahlén, “Transfer function for excitation system and automatic voltage regulator”, VG Power AB, 2004.
7 Appendix
These versions of the appendix concerning the results are a limited edition, the complete result appendixes can be given at a request.
Appendix 1 Settings on the the different machines Appendix 2 the matlab-code.
Appendix 3 example of AVR and PSS
Appendix 4 Pictures of the power cabinet that was built Appendix 5 Pictures of the generator and damping bars Appendix 6 Results from the simulations, study 1, machine 1 Appendix 7 Results from the simulations, study 1, machine 2 Appendix 8 Results from the simulations, study 1, machine 3 Appendix 9 Results from the simulations, study 1, machine 4
Appendix 10 Results from the simulations, study 2, AVR 1 (simple gain)
Appendix 11 Results from the simulations, study 2, AVR 2 (gain and derivative) Appendix 12 Results from the simulations, study 2, AVR 3 PID-regulator (gain
integrating and derivative)
Appendix 13 Results from study I, machine I. A table of Ks and Kd with the new mathematical model compared to the old model.
Settings on the the different machines, Appendix 1
% FUNDAMENTEL PARAMETRES IN [P.U] (weak damper) MACHINE I Lad = 0.65; Ra = 0.003;
Laq = 0.35;
Lfd = 0.1; Rfd = 0.0003;
L1d = 0.05; R1d = 0.005;
L1q = 0.2; R1q = 0.01;
Ll = 0.15;
Lad_u = 0.75;
% FUNDAMENTEL PARAMETRES IN [P.U] (strong damper) MACHINE II Lad = 0.65; Ra = 0.003;
Laq = 0.35;
Lfd = 0.1; Rfd = 0.0003;
L1d = 0.05; R1d = 0.005;
L1q = 0.05; R1q = 0.005;
Ll = 0.15;
Lad_u = 0.75;
% FUNDAMENTEL PARAMETRES IN [P.U] (weak damper, large synchronous reactance) MACHINE III
Lad = 0.95; Ra = 0.003;
Laq = 0.55;
Lfd = 0.1; Rfd = 0.0003;
L1d = 0.05; R1d = 0.005;
L1q = 0.2; R1q = 0.01;
Ll = 0.15;
Lad_u = 1.05;
% FUNDAMENTEL PARAMETRES IN [P.U] (strong damper, large synchronous reactance) MACHINE IV
Lad = 0.95; Ra = 0.003;
Laq = 0.55;
Lfd = 0.1; Rfd = 0.0003;
L1d = 0.05; R1d = 0.005;
L1q = 0.05; R1q = 0.005;
Ll = 0.15;
Lad_u = 1.05;
% FUNDAMENTEL PARAMETRES IN [P.U] (machine without damping) MACHINE V Lad = 0.65; Ra = 0.003;
Laq = 0.35;
Lfd = 0.1; Rfd = 0.0003;
L1d = 0.05; R1d = 1;
L1q = 0.2; R1q = 1;
Ll = 0.15;
Lad_u = 0.75;
Appendix 2
The Matlab-code
% File: DAE_3.m
% Formulation of the system of differential-algebraic equations
% of the SMIB-system before and after the fault.
function [dae] = DAE_3(T,X_Y,PARAM);
H = PARAM(1);
Pm0 = PARAM(2);
w_base = PARAM(3);
Rfd = PARAM(4);
Lfd = PARAM(5);
R1d = PARAM(6);
L1d = PARAM(7);
R1q = PARAM(8);
L1q = PARAM(9);
Xdb = PARAM(10);
Xqb = PARAM(11);
Ll = PARAM(12);
Ra = PARAM(13);
efd0 = PARAM(14);
EB = PARAM(15);
RE = PARAM(16);
XE = PARAM(17);
noStAE = PARAM(18);
Et = PARAM(19);
Lad = PARAM(20);
noStpss = PARAM(21);
AE = PARAM(22:30,1); %cheating code, find length of AE pss = PARAM(23+8:end,1); %cheating code, find length of pss
A_ae = zeros(noStAE,noStAE);
B_ae = zeros(noStAE,1);
C_ae = zeros(1,noStAE);
D_ae = 0;
for i=1:noStAE
A_ae(1,i) = AE(i);
A_ae(2,i) = AE(noStAE + i);
end
for i=1:noStAE
B_ae(i,1) = AE(2*noStAE + i);
end
end
D_ae = AE(4*noStAE + 1);
A_pss = zeros(noStpss,noStpss);
B_pss = zeros(noStpss,1);
C_pss = zeros(1,noStpss);
D_pss = 0;
for i=1:noStpss
A_pss(i,1) = pss(noStpss*(i-1)+1);
A_pss(i,2) = pss(noStpss*(i-1)+2);
A_pss(i,3) = pss(noStpss*(i-1)+3);
A_pss(i,4) = pss(noStpss*(i-1)+4);
A_pss(i,5) = pss(noStpss*(i-1)+5);
end
for i=1:noStpss
B_pss(i,1) = pss(5*noStpss + i);
end
for i=1:noStpss
C_pss(1,i) = pss(6*noStpss + i);
end
D_pss = pss(7*noStpss + 1);
% State variables DELTA = X_Y(1);
OMEGA = X_Y(2);
PSI_fd = X_Y(3);
PSI_1d = X_Y(4);
PSI_1q = X_Y(5);
xAE = X_Y(6:6+noStAE-1);
xpss=X_Y(7:7+noStpss-1);
PSI_ad = X_Y(8+noStAE+noStpss-2);
PSI_aq = X_Y(9+noStAE+noStpss-2);
ed = X_Y(10+noStAE+noStpss-2);
eq = X_Y(11+noStAE+noStpss-2);
Ladsb = Xdb - Ll;
Laqsb = Xqb - Ll;
% Infinite bus voltage in machine reference frame EBd = EB*sin(DELTA);
EBq = EB*cos(DELTA);
% Subtransient voltage sources Edb = Laqsb*(PSI_1q/L1q);
Eqb = Ladsb*(PSI_fd/Lfd + PSI_1d/L1d);
% CALCULATION OF STATOR CURRENTS id and iq (Kundur sid. 783) RT = Ra + RE;
XTd = XE + Xdb; XTq = XE + Xqb;
D = RT^2 + XTd*XTq;
EdN = Edb + EBd; EqN = Eqb - EBq;
id = (XTq*EqN - RT*EdN)/D;
iq = (RT*EqN + XTd*EdN)/D;
% Air-gap (breaking) power Pe = PSI_ad*iq - PSI_aq*id;
Vref = 1;
% Differential equations: xdot = f(x,y)
f_xy = [w_base*OMEGA;
(1/(2*H))*(Pm0 - Pe);
w_base*(efd0 + (Rfd/Lad)*(C_ae*xAE + ...
D_ae*(Et-sqrt(ed*ed+eq*eq))) + (PSI_ad - PSI_fd)*Rfd/Lfd);
%PSI_fd/dt (general solution)
w_base*(PSI_ad - PSI_1d)*(R1d/L1d);
%PSI_1d/dt
w_base*(PSI_aq - PSI_1q)*(R1q/L1q);
%PSI_1q/dt
A_ae*xAE + B_ae*(Et-sqrt(ed*ed+eq*eq)+C_pss*xpss+D_pss*OMEGA);
%dxAE/dt
A_pss*xpss+B_pss*OMEGA];
%dxpss/dt
g_xy = [PSI_ad + Ladsb*(id - PSI_fd/Lfd - PSI_1d/L1d); % PSI_ad PSI_aq + Laqsb*(iq - PSI_1q/L1q); % PSI_aq ed - EBd - RE*id + XE*iq; % ed eq - EBq - XE*id - RE*iq]; % eq dae=[f_xy;
g_xy;];
% FILE: ExtraktionKsKd.m
% Created by Martin Ranlöf 2010-08-27
% Extract Ks and Kd from the simulated responses of deltaTE, deltaOMEGA and
% deltaDELTA following a small disturbance.
% A least-square fitting approach is used.
% The script is designed to operate on data output from the MATLAB script
% "Simulering.m".
% Subscript P = "Post"
start = 200;
time = T_P(start:end,1) - T_P(start,1); % [s]
DELTA_res = X_Y_P(start:end,1); % [rad]
dOMEGA = w_base*X_Y_P(start:end,2); % [rad/s]
PSI_fd_resP = X_Y_P(start:end,3);
PSI_1d_resP = X_Y_P(start:end,4);
PSI_1q_resP = X_Y_P(start:end,5);
PSI_ad_resP = X_Y_P(start:end,13);
PSI_aq_resP = X_Y_P(start:end,14);
ifd_resP = (PSI_fd_resP-PSI_ad_resP)/Lfd; % [p.u] field current
i1d_resP = (PSI_1d_resP-PSI_ad_resP)/L1d; % [p.u] D-damper current
i1q_resP = (PSI_1q_resP-PSI_aq_resP)/L1q; % [p.u] D-damper current
id_resP = -(PSI_ad_resP - Lad*(ifd_resP + i1d_resP))/Lad; % [p.u] d-axis current
iq_resP = -(PSI_aq_resP - Laq*i1q_resP)/Laq; % [p.u] q-axis current
Te_anp = PSI_ad_resP.*iq_resP - PSI_aq_resP.*id_resP; % [p.u] moment som skall anpassas
% Calculate dTe and dDELTA Te_mean = Te0;
dTe = Te_anp - Te_mean;
DELTA_mean = delta0;
dDELTA = DELTA_res - DELTA_mean;
% Create time-interval vector time2 = time(2:end);
dT = time2 - time(1:end-1);
% Calculate least square integrals (see paper by Alden and Shaltout) B1 = sum(dTe(1:end-1).*dDELTA(1:end-1).*dT);
B2 = sum(dTe(1:end-1).*dOMEGA(1:end-1).*dT);
A1 = sum(dDELTA(1:end-1).*dDELTA(1:end-1).*dT);
A2 = sum(dDELTA(1:end-1).*dOMEGA(1:end-1).*dT);
A4 = sum(dOMEGA(1:end-1).*dOMEGA(1:end-1).*dT);
% Solve system of equations to find Ks and Kd
Ks_anp = (A4*B1 - A2*B2)/(A1*A4 - A2^2); % [p.u./rad]
Kd_anp = (B1 - A1*Ks_anp)/A2; % [p.u./rad/sec]
dTe_cntrl = Ks_anp*dDELTA + Kd_anp*dOMEGA;
figure(8)
plot(time,dTe,'y') hold on
plot(time,dTe_cntrl,'k') xlabel('time [s]');
ylabel('Ks+Kd') hold on
