Transient Stability Analysis of the IEEE-9 Bus System under Multiple Contingencies

The determination of the transient stability of an electric power system is a crucial step in power system analysis. This paper investigates the transient stability of an IEEE-9 bus system consisting of three generators and nine buses. At first, a load flow analysis is conducted in order to determine the prefault conditions. Secondly, fault analysis is performed to analyze post fault conditions like the fast fault clearing time and load switching in order to determine the system stability. For transient stability analysis, Euler and Runga methods are compared and applied on the frequency and rotor angle of the system to analyze the system variations under different fault conditions. The simulations were done on the Power World Simulator (PWS) software. It is concluded that Critical Fault Clearing Time (CFCT) is a very important factor in keeping the power system within the stability bounds. A slight increase in Clearing Time (CT) from the critical value causes un-synchronism. Keywords-transient stability; IEEE-9 bus system; critical time; rotor angle; power world simulator

INTRODUCTION A power system is designed to supply continuous power with good quality by maintaining voltage stability even in the presence of lightning, short-circuit, or ground faults [1]. Due to these faults, one or many generators may act abnormally causing a wide gap between demand and supply [2,3]. Globally, power demand is escalating, which results in the installation of bigger systems raising, among others, the issue of synchronization [4]. Therefore, an increase in electricity demand has made the forecasting of upcoming development very complex [5]. This increase in electricity demand also causes power line overloading [6]. The overloading causes the voltage at each bus to reduce and the efficiency of the generators to supply power to the system under a power system abnormal operation representing a contingency situation, decreases [7]. Operation and system operators need to take steps to bring back the power system in safe operation. To resolve these glitches of overloaded transmission lines and increasing load demand, two resolutions can be useful: either increasing the generated power or constructing new transmission lines [8]. The analysis of transient stability is a vital step in designing a power system and its components. Transient stability assesses the electrical power system's capability to sustain disturbances and to keep the system in normal operating conditions [9]. The power system instabilities can be faults such as transmission line short circuits, losses of generator, and losses of load [10]. They all result in a large deviation of the generator rotor angle and also effect power flow, bus voltage, and other system variables causing a partial or total loss of the transmission network [11]. Steady state stability only deals with operating environments, whereas transient stability deals with operating environments and disturbances [12]. A frequent investigation is required for various types of disturbances [13]. Classification of electrical power system stability is critical for the good understanding of power reliability. Power system stability is consequently classified according to the nature of system instability (voltage and frequency stability), the magnitude of disturbances (large and small disturbances), and long term and short-term stability (time-frame stability) [14]. The instability of the power system is a complex issue since it can take many forms and can be related to many factors. A detailed classification defines three main categories of power system stability. The considerations that are taken to counter system instability include the physical behavior of the instability, loss of synchronization, low bus voltages, high frequency deviation, size of the disturbance, and the measures taken to improve power system stability [14].
Rotor angle stability is the competency of the Power System to sustain synchronism in a system after a disturbance has occurred [16]. Being endangered to instabilities, the synchronous machines start to fluctuate w.r.t. each other. If one machine starts running faster, the deviation in angular position will increase and the fast machine will supply more load. System speed will then decrease allowing the other machines' speed to increase until an equilibrium is reached [17]. Instability occurs when the equilibrium is not reached and the speed of some machines increases until these machines are tripped [15,18]. Rotor angle stability is a highly nonlinear, multi-dimensional problem [3]. Instability can be categorized into two different parts: aperiodic instability and oscillatory instability [19]. Aperiodic instability is caused by a lack of synchronizing torque that is associated with rotor angle deviation, while oscillatory instability is affected by the lack of damping torque associated with speed deviation. Instability problems are aperiodic and are mainly due to insufficient synchronizing torque [20][21]. Rotor angle stability consists only of short-term phenomena from a few cycles or a few seconds up to a few tens of seconds [18,15]. Large disturbance angular stability (transient stability) is a subsection of rotor angle stability and contains the function of the power system to keep angular stability after major disturbances [22]. The maximum time during which a system holds disturbances without losing stability is called Critical Clearing Time (CCT) [23]. Due to the highly nonlinear nature of transient stability, time domain simulation is used to solve the algebraic/differential equations using a step-by-step calculation procedure. The direct method can also be used to determine stability using Newton Raphson's second method [24]. This method is not often used due to the complication of finding an adequate function and its inability of defining a practical stability domain [3,15]. Pattern recognition can also be a method for transient stability analysis by referring to past experience and applying it to current stability properties [3]. Small disturbance angular stability is another subsection of rotor angle stability and it comprises of functions of the power system that sustain stability after experiencing small disturbances [25]. A small disturbance is considered if the system still exhibits the original system kinematics under disturbances. Modal analysis method can be applied when a system is linearized. The modal analysis calculates the characteristics of means (mode shapes, transfer functions and coupling coefficients) that are useful for improving damping. It uses the time domain simulation to visualize the results [26]. The branching method is another way to evaluate the effect of large changes in small signal stability, unlike modal analysis, where all system parameters are fixed [27][28]. As mentioned above, two types of rotor angular instability can occur. The presence of a complex conjugate eigenvalue having positive real part causes vibration instability. If the eigenvalues are real and positive, instability is an acyclic type.
Frequency stability is the capability of the Power System to uphold stable frequencies following disturbances due to the discrepancy between generation and demand [30]. The system is said to have frequency stability, if it can balance or restore power and load with minimal load trips. Time domain simulations are mainly used for frequency stability analysis as an appropriate representation of dynamic devices [15]. Numerical integration techniques can be applied to obtain approximate solutions of nonlinear differential equations. Many algorithms are available for numerical integration such as Euler's method of integration, modified Euler's method of integration, and Runga Kutta method. Euler's method is the simplest and the least accurate. Euler's method is a first order method. It is a straight-forward method that estimates the next point based on the rate of change at the current point and it is easy to code. Notably, it is unconditionally unstable for undamped oscillating systems in space discretization. It can be used for basic numerical analysis, but may fail in complex problems and/or boundary conditions. This method is not commonly used for spatial discretization but is some times used in time discretization and is not recommended for hyperbolic differential equation. The order of convergence of this scheme with grid refinement is very poor. Extending Euler's method to higher order method is easy and straight forward.
The time-domain numerical integration is not suitable for on-line security analysis due to the long CPU run times for simulation. A typical time-domain numerical integration of 2s takes more than 120s depending on the step size of the integration. Larger step size that reduces time causes inaccurate and less reliable results. The Euler method is actually the simplest Runga Kutta (RK) method (1 stage, first order). The RK methods approximate the Taylor series solution, however unlike the formal Taylor series solution, the RK methods do not require explicit evaluation of derivatives higher than the first. The effects of higher derivatives are included by several evaluations of the first derivative. Depending on the number of terms effectively retained in the Taylor series, we have RK methods of different orders. Higher order accurate RK methods are multi-stage because they involve slope calculations at multiple steps at or between the current and next discrete time values. The next value of the dependent variable is calculated by taking a weighted average of these multiple stages based on a Taylor series approximation of the solution. The weights in this weighted average are derived by solving non-linear algebraic equations which are formed by requiring cancellation of error terms in the Taylor series. Developing higher order RK methods is tedious and difficult without using symbolic tools for computation. The most popular RK method is RK4 since it offers a good balance between the accuracy order and computation cost. Usually the error in Euler's method is higher than higher order RK method's. If the exact solution to the differential equation is a polynomial of order n, it will be solved exactly by an n-th RK method. For example, Euler will be exact if the solution is a line. RK4 will be exact if the solution is a polynomial of degree 4 or less.
In this paper, transient stability analysis is studied under three-phase balance fault. RK and Euler methods are used for analysing the swing behavior of the system. A large number of simulations have been performed during the planning stage of the system to obtain its accurate result. The simulation of the IEEE-9 bus system model was carried out in Power World Simulator. The proposed work is tested on the IEEE-9 bus power systems under normal and abnormal operating conditions [31].

II. POWER FLOW STUDY
Transient stability analysis requires the knowledge of the pre-fault voltage magnitude. Key information in the power flow study includes the voltage magnitude of buses, the phase angle of bus voltage, the transmission line power flow and losses, the generator bus power flow, etc. [32]. From load flow analysis, pre-fault conditions can be obtained by using Newton Raphson method or fast de-coupled method. In Newton-Raphson method, convergence is not affected by the reserve bus selection. This method starts with a preliminary guessing of all unknown parameters [33]. Firstly, a study is performed for the existence of a load bus and a distant PQ bus. For the ith bus: III. METHODOLOGY Transient stability addresses the effects of accidental outages such as sudden outages, fault incidences, abstraction or application of sudden loads [34]. The standard IEEE-9 bus system given in [7] was used. Load flow is completed on the system model. Load flow analysis was performed using the adaptive Newton Raphson method as shown in Table I [36]. After load flow analysis, the voltage at a bus, its frequency, voltage angle and the generator rotor angle are shown in Table  II. Multiple contingencies cleared before and after CCT are applied and the output results are obtained by Power World Simulator using the Runga and Euler methods. The results of both approaches are analyzed and compared. Important observations can be made about Jacobian matrix elements. If there is no linking between the ith bus and the kth bus, ܻ݅݇ =0. This procedure remains until the stop condition is met.

IV. RESULTS AND DISCUSSION
The current section at first considers five cases of fault clearness after CCT. Simulation is also performed with respect to blackouts in which, a system can no longer supply load [37] and two cases of fault clearness before CCT. The simulation results of the rotor angle and frequency for the different contingencies of the multi-machine power system of [7] are shown in Figures 1-13 A three phase short circuit balanced fault occurred on the transmission line between bus 5 and bus 7 at 1.0s and was cleared at 1.1s. The rotor angle of all three generators is shown in Figure 1.  Rotor angle deviation is maximum in generators 2 and 3 due to their smaller distance from the fault point in comparison with generator 1. After some time, the system becomes stable. Increasing the CCT also increases the angle difference. This means that the system is in an unstable mode. The system will be out of synchronization, if more time is required to clear the fault, so the CCT must be very small to retain the system in synchronization [38]. The frequency of all buses is shown in Figure 2. The frequency of the connected system also changes due to the fault on the transmission line between buses 5 and 7. The frequency of the system increases because of less loading on the system as the load at bus 5 is removed due to fault.

2) Case II
In this case, the load is increased by 50% on bus 5 at 1.0s and it is immediately removed at 1.1s. The rotor angles of the three affected generators are shown in Figure 3. After some time the system becomes stable. The frequency response of the system at all buses is shown in Figure 4. The frequency of the connected system also changes due to sudden load changes at bus 5. It decreases because the system has a heavy load at bus 5 and increases due to sudden loading. Rotor angles of the three generators for a sudden 50% increase in load on bus 5 at 1s. Bus frequencies for a sudden 50% increase in load by 50% on bus 5 at 1s.

3) Case III
In this case, the load is decreased by 50% on bus 5 at 1.0s and it is removed at 1.1s. The rotor angle of the three generators is shown in Figure 5. Rotor angles of the three generators for a sudden 50% decrease in load on bus 5 at 1s. Fig. 6.
Bus frequencies for a sudden 50% decrease in load on bus 5 at 1s.
There is a rotor angle deviation in all three generators, however after some time the system becomes stable. But increasing the CCT time also increases the angle difference. This means that the system is in an unstable mode. The system will be out of synchronization if more time is required to clear the fault, so the CCT must be very short to retain the system in synchronization [38]. The frequency response of the system buses is shown in Figure 6. The frequency of the connected system also changes due to the sudden load change at bus 5.

4) Case IV
A three-phase balanced fault occurs on bus 5 at 1.0s and the fault is cleared at 1.1s. The rotor angle of the three generators is shown in Figure 7. Rotor angle deviation is maximum in generators 2 and 3 due their smaller distance from the fault point in comparison with generator 1. After some time the system is stable. The frequency of all buses is shown in Figure  8. The frequency of the connected system changes due to the three-phase balanced fault at bus 5. The frequency of the system increases because the system has less load at bus 5 due to the three-phase balanced fault. Rotor angles of the three generators for a three-phase balanced fault on bus 5 at 1s.

5) Case V
In this case, we consider an outage of generator 2 at 1.0s. The rotor angles of the generators are shown in Figure 9. There is a rotor angle deviation in generators 1 and 3. Generator 2 has no rotor angle deviation due to the outage from the system. The frequency of all buses is shown in Figure 10. The frequency of the connected system also changes due to the outage of generator 2. The frequency of the system decreases because the system has less generation due to the outage.

1) Case VI
Power outages in the power system can occur by several causes, such as an overload of the transmission line, ice coating of the track, protection or malfunction of the control system, etc. Blackouts as shown in Figure 11 can be resolved through appropriate control schemes in a power system.

C. Fault Clearness before Critical Clearing Time
When the load on bus 5 increases, a transmission line trips between bus 4 and bus 5. As the line trips, the rest of the line needs to deliver power, which consumes more reactive power and reduces the load center voltage without affecting the frequency. If the reactive power is not enough, the voltage drop increases under the load line load condition. Therefore, the voltage is the key parameter of the power system and not frequency. In addition, the voltage drop at the load center shows that the system will experience low frequency after it crashed [39]. To analyze the effect of fault clearing on transient stability, fault is being removed after 3.05s. The maximum rotor angle changes by 116.874 o by the Runga method and 118.623 o by the Euler method. Comparing the two methods, the Runga method provides a faster response as shown in Figure 14. The maximum change in rotor angle is 118.623 o , the system is cleared after 3.05s and after some time the system is stable. Increasing the CCT increases the rotor angle difference. This means that the power system will be then in an unstable mode and it will be out of synchronization and more time will be required to clear the fault. Because of this reason the CCT must be kept short in order to keep the system in synchronization [38,40]. The simulation results of the rotor angle and V_Angle for different contingencies on the multimachine power system are shown in Figures 12-13.

1) Case VII
A balanced three phase fault occurs from bus 2 to bus 7. Figure 12(a) shows that the voltage angle stays constant during the transient at -2.835025 o (V_ Angle) throughout the time interval of 3.05s. This undesired change is due to the occurrence of the transient. Figure 12(b) shows that the per-unit value is not changed from its initial value and remains the same at 1p.u. during the time interval of the transient. Figure 13(a) shows that the speed of the generator remains constant during the mismatch between the mechanical and electrical power. In the initial state, it changes from 0 to 1 giving a new constant response across time which means that the machine has selected a new value for synchronous operation after the fault clearance while the generator rotor angle continues to remain in synchronism after the disturbance. Due to the increase of the current in the rotor field, the graph steadily shows an MVar rise blink on the x-axis, which marks the y-axis somewhere above 2400. This causes the magnetic coupling of the rotor to the stator and increases the MVar output during this instant. Figure 13(b) shows the received MW of the generator, a step blink at 1s marks somewhere around 8300 the y-axis and finally chooses the new increased and damped MW value. During the mentioned transient, the p.u. voltage of the generator remains almost constant. The generator rotor angle will swing at zero during the 3.05s transient and the speed will remain constant. As generator's reactive power rises abruptly above 24000 up to 25000 for 1s and remains at the constant during the CCT.  V. CONCLUSION In order to analyze the power system behavior, different types of fault have been simulated on the IEEE-9 bus system for transient stability analysis. For this purpose, load flow analysis and fault analysis were performed. The simulation results show that the balanced three fault is being removed after 3.05s by Runga Kutta and Euler methods, however the maximum rotor angle changes to 116.874 o by the Runga Kutta method and 118.623 o by the Euler method. Comparing the two methods, the Runga method provides a faster response than the Euler method. During transient stability analysis it was shown that increasing CCT increases the rotor angle difference. Moreover, the results of different case studies of fault occurrence on power system show that fault clearness after CCT makes the system unstable while fault clearness before CCT keeps the system in synchronism.