Response Based Emergency Control System for Power System Transient Stability

Abstract

A transient stability control system for the electric power system composed of a prediction method and a control method is proposed based on trajectory information. This system, which is independent of system parameters and models, can detect the transient stability of the electric power system quickly and provide the control law when the system is unstable. Firstly, system instability is detected by the characteristic concave or convex shape of the trajectory. Secondly, the control method is proposed based on the analysis of the slope of the state plane trajectory when the power system is unstable. Two control objectives are provided according to the methods of acquiring the far end point: one is the minimal cost to restore the system to a stable state; the other one is the minimal cost to limit the maximum swing angle. The simulation indicates that the mentioned transient stability control system is efficient.

1. Introduction

This is a follow-up paper of a series on closed loop control systems for power system transient stability. The previous work proposed a real-time approach to detect transient instability with high accuracy and wide applicability, while this paper mainly focuses on how and where to control the system after the instability is detected.

In China, transient stability control systems acquire some characteristic variables according to a large number of offline calculations and identify the stability by the combination of these characteristic variables. For unstable situations, the control quantity of different situations is obtained by repeated simulations to develop a cure table. When disturbances occur, protection devices operate on the basis of the cure table [1,2].

The samples consist of different situations, including power flow, grid topology and fault conditions making it hard to ensure an effective cure table. In order to decrease the number of samples and improve the utility of the cure table, a simulation which is helpful to detect transient instability and form control law adopts measurement results to update the online power flow with a refresh cycle of 3–5 min [3,4,5,6]. Corresponding approaches has been introduced in [5,6], including preplanned remedial action and system integrity protection schemes, software processes and hardware requirements. To speed up the simulation, the transient energy function method, the equal area criterion, topological energy function or quasi-real-time online transient analysis are employed in the time domain simulation, which can greatly improve the efficiency [7,8,9,10,11,12].

In any case these methods need repeated simulations to calculate control laws by taking into account the changes in operating conditions or parameter variations. The validity of the cure table depends on the similarity between the anticipated faults and actual faults and the accuracy of the model parameters. However, a practical operating model (especially a model of the load) or the system parameters are hard to acquire. Meanwhile, the rapid development and the operation variations lead to a large number of calculation samples.

With the development of computer science and communication technology, especially the application of phase measurement unit (PMU) based on global positioning system (GPS) in power systems, wider-area measurement systems (WAMS) offer an opportunity to develop real-time protection and control systems [13]. As a result, recently there has been a focus on online instability control, out-of-step protection and their corresponding theories.

It has been reported that the geometric characteristics of the system trajectory can be employed to detect the system instability [14]. Depending on the nature of the stability studied, Girgis has found that the characteristic shape (concave or convex) of a surface, on which a post-fault transient trajectory lies, can be used as an index for online instability detection [15]. The assumption is proved in a phase portrait, and an instability detection method is presented, which is independent of the network structure, system parameters and model because it only uses observation data. During the process of proof, a damping coefficient is considered in the system model to keep it more reasonable [16,17,18]. However, the discrete expression of the detection method in [18] requires differential calculations. The curve of the index looks ragged. For a stable situation after disturbances, it may cause some erroneous judgments.

This paper aims to introduce a real-time transient stability control system which employs the generator speed, power angle and unbalanced power to detect the system instability and develop a control strategy. In this paper, the system instability is detected by the characteristics concave or convex shape of the trajectory. Then, the control law is achieved based on the slope of the state plane trajectory, when the power system is unstable. Because the expected variables can be obtained from the trajectory information and the calculations are very simple, the control system can be realized in real time. The simulation indicates that the mentioned transient stability control system is fast, efficient, and realizable.

2. Instability Detection

2.1. Identification Method of Transient Instability for an Autonomous SMIB System

Previous research has found the geometric characteristics of the trajectory in a phase portrait can be employed to determine system instability through a great deal of simulation. It is found that stable trajectories are always concave with respect to the post-fault stable equilibrium point (SEP) and that unstable trajectories are convex with respect to the post-fault SEP immediately or a short time after the fault-clearing time as seen in Figure 1. The point at which the geometric characteristic of a system trajectory is convex, is defined as the no return point (NRP). All of these points (NRP) form of the interface that is defined as the no return point interface (NRPI).

Figure 1. System trajectories of a SMIB.

Figure 1. System trajectories of a SMIB.

It should be particularly pointed out that the definition in this paper is different from the usual definition. As a rule, the point where a trajectory shifts to accelerated progress ($\mathrm{\Delta }P=0$, $\mathrm{\delta }={\mathrm{\delta }}_{uep}$, $\text{Δ}P$ is the unbalanced power and ${\mathrm{\delta }}_{uep}$ is the unstable equilibrium point) is defined as the NRP. The difference between both definitions is shown in Figure 1. Although Figure 1 is obtained for a specific set of system parameters, the characteristics of the system trajectory are typical of a single machine infinite bus (SMIB) system. Here, we use a SMIB system as an example to explain the foregoing concept. For the SMIB system, the system dynamics can be expressed as:

$$\begin{array}{l}\stackrel{•}{\mathrm{\delta }}={\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega }\\ \stackrel{•}{\mathrm{\Delta }\mathrm{\omega }}=[P_m-P_{emi}\sin \mathrm{\delta }-D\mathrm{\Delta }\mathrm{\omega }]/M\end{array}$$

(1)

where $\mathrm{\delta }$ and $\mathrm{\Delta }\mathrm{\omega }$ are the generator angle and angular velocity with respect to a synchronous frame; $D$ is the generator damping coefficient; $M$ is the generator inertia; $P_m$ is the generator mechanical power input; and $P_{emi}$ is the generator maximum electric power output.

On the basis of the definition of NRP, the tangent slope at any point on the trajectory of system Equation (1) is given as follows:

$$\frac{\raisebox{1ex}{$d\mathrm{\Delta }\mathrm{\omega }$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.}{\raisebox{1ex}{$d\mathrm{\delta }$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.}=\frac{d\mathrm{\Delta }\mathrm{\omega }}{d\mathrm{\delta }}=k(\mathrm{\delta },\mathrm{\Delta }\mathrm{\omega })=\mathrm{\Delta }P/(M{\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega })$$

(2)

Consequently, the one-order derivative of the tangent slope at any point of the trajectory with respect to the generator angle (that is, the two-order derivative of the angular velocity with respect to the generator angle) can be expressed as:

$$l=\frac{dk}{d\mathrm{\delta }}=k^{\prime }(\mathrm{\delta },\mathrm{\Delta }\mathrm{\omega })=\frac{[-P_{emi}\cos (\mathrm{\delta })-D\frac{d\mathrm{\Delta }\mathrm{\omega }}{d\mathrm{\delta }}]M{\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega }-{[P_m-P_{emi}\sin \mathrm{\delta }-D\mathrm{\Delta }\mathrm{\omega }]}^2/\mathrm{\Delta }\mathrm{\omega }}{{(M{\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega })}^2}$$

(3)

On the ground of the definition of NRP, the point of system trajectory at which $l=0$ is actually NRP. The interface (NRPI) that all of NRP constitute can be described as:

$$-P_{emi}\cos \mathrm{\delta }M{\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega }-D[P_m-P_{emi}\sin \mathrm{\delta }-D\mathrm{\omega }]-{[P_m-P_{emi}\sin \mathrm{\delta }-D\mathrm{\Delta }\mathrm{\omega }]}^2/\mathrm{\Delta }\mathrm{\omega }=0$$

(4)

Figure 2 shows the NRPI divides the whole plane into two regions, i.e., a convex area and a concave area. In virtue of the mathematical definition, if $l>0$, we consider the geometrical characteristics of the post-fault trajectory is convex, while it is considered as concave if $l<0$.

We define the concavity and convexity of phase trajectory as follows:

(1)

The phase trajectory is convex if $l\cdot \mathrm{\Delta }\mathrm{\omega }>0$;

(2)

The phase trajectory is concave if $l\cdot \mathrm{\Delta }\mathrm{\omega }<0$;

(3)

The trajectory is on the inflexion point if $l\cdot \mathrm{\Delta }\mathrm{\omega }=0$.

In order to facilitate real-time computing, the unstable criterion can be expressed as:

$$\mathrm{\tau }=l\cdot \mathrm{\Delta }\mathrm{\omega }=\frac{dk}{d\mathrm{\delta }}\cdot \frac{d\mathrm{\delta }}{dt}\cdot \frac{1}{{\mathrm{\omega }}_0}=\frac{dk}{dt}\cdot \frac{1}{{\mathrm{\omega }}_0}\ge \text{0}$$

(5)

The discrete form of the index is given by:

$$\mathrm{\tau }(i)=\frac{\mathrm{\Delta }P(i)}{M\mathrm{\Delta }\mathrm{\omega }(i)}-\frac{\mathrm{\Delta }P(i-1)}{M\mathrm{\Delta }\mathrm{\omega }(i-1)}$$

(6)

Figure 2. Relationship of NRPI and the critical post-fault trajectory.

Figure 2. Relationship of NRPI and the critical post-fault trajectory.

If $\mathrm{\tau }$ is always less than zero, the system will be stable, while the system will lose synchronization immediately or a short time after the fault-clearing time if $\mathrm{\tau }$ is larger than zero. The proposition has proved that the proposed index is a sufficient and necessary condition for instability detection of SMIB system.

It should be pointed out that the moment when the trajectory passes through the NRI is just the detection time. As shown in Figure 1, it can be found that the unstable trajectory is sure to pass through the NRI before it arriving at dynamic saddle point (DSP). Thus, by employing the characteristic concave or convex shapes of the trajectory one can easily and quickly identify the transient instability.

2.2. Identification Theory of Transient Instability for a Non-Autonomous SMIB System

The equations of an electromechanical transient process for non-autonomous SMIB system can be described as follows:

$$\begin{array}{l}\frac{d\mathrm{\delta }}{dt}={\mathrm{\omega }}_0\mathrm{\Delta }\mathrm{\omega }\\ M\frac{d\mathrm{\Delta }\mathrm{\omega }}{dt}=P_m(t)-P_e(t,\mathrm{\delta })=\mathrm{\Delta }P(t,\mathrm{\delta })\end{array}$$

(7)

where the right function of the motion equation involves the variable of time, and thus $\mathrm{\Delta }P(t,\mathrm{\delta })$ doesn’t represent a strict sine function any more. As a consequence, the index $\mathrm{\tau }$ greater than zero in a certain time can only indicate the instability of the autonomous system under that parameter condition, but not the non-autonomous system:

$$\begin{array}{l}{P}_e(t_i,\mathrm{\delta })={\mathrm{\lambda }}_0(t_i)-{\mathrm{\lambda }}_1(t_i)\cos \mathrm{\delta }-{\mathrm{\lambda }}_2(t_i)\sin \mathrm{\delta }\\ \mathrm{\Delta }P(t_i,\mathrm{\delta })=P_c(t_i)-{\mathrm{\lambda }}_1(t_i)\cos \mathrm{\delta }-{\mathrm{\lambda }}_2(t_i)\sin \mathrm{\delta }\end{array}$$

(8)

To fortify the accuracy of the detection, the index is established as follows. The feature index $\mu$ of trajectory in state-plane of unbalanced power and power angle is provided.

The feature index of the trajectory in state plane of unbalanced power and power angle is defined as:

$${\frac{dl}{dt}|}_{l=0}=\frac{d^2(\mathrm{\Delta }P/M)}{d{\mathrm{\delta }}^2}$$

(9)

The discrete form is proposed as:

$$\mathrm{\mu }=\frac{\mathrm{\Delta }P(i)-\mathrm{\Delta }P(i-1)}{\mathrm{\delta }(i)-\mathrm{\delta }(i-1)}-\frac{\mathrm{\Delta }P(i-1)-\mathrm{\Delta }P(i-2)}{\mathrm{\delta }(i-1)-\mathrm{\delta }(i-2)}>0$$

(10)

The index means that the trajectory enters the convex region and does not return to the concave region in a short time. If the trajectory satisfies the index, the trajectory is unstable possibly.

Eventually, the instability criterion for a non-autonomous SMIB system can be described as:

$$\mathrm{\tau }>0\&\mathrm{\mu }>0$$

(11)

The detection method launches when the power system is disturbed. If the maximum swing angle is less than 50°, the detection method automatically ends.

3. Control Method

There are many transient stability control measures, including generator shedding and load shedding [19,20]. In North America, generator shedding has been proved to be one of the most effective discrete supplementary control means for maintaining stability [21]. In this paper, generator shedding is mainly discussed too.

3.1. The Slope of the State-Plane Trajectory

According to Equation (2), the expression of the slope is related to the unbalanced power which can be changed by generator shedding. It is supposed that the system returns to a stable state by generator shedding. If the value of the slope at control time is obtained, the corresponding generator shedding can be calculated by Equation (2). To obtain the value of the slope at control time, the characteristic of the slope is analyzed below. A single machine infinite bus system is shown in Figure 3.

Figure 3. The topology network of the SMIB system.

Figure 3. The topology network of the SMIB system.

Under normal conditions, the system operates at a stable equilibrium point (SEP). The initial state of the generator is 0.812 rad, 0 rad/s. A three-phase grounding fault occurs on one of the two transmission lines, and the fault is cleared by switching off the line. The clearing time is 0.18 s (critical clearing time) and 0.22 s. When the system is unstable, different control quantities are taken. The trajectory of the state-plane and the slope of the trajectory are shown in Figure 4 and Figure 5.

For the black curve, the system is stable without control. The trajectory swings back at the far end point (FEP) which exists only in a stable system. The slope of the trajectory keeps decreasing with the mutation at the FEP.

For the red curve, the system is unstable without control. The slope begins to increase after the inflexion point and grows to zero at the dynamic saddle point (DSP).

For the green curve, the system is unstable with insufficient control. A sudden change of the slope is caused by the control at the control time. The larger the control quantity is, the greater the change of the slope at the control time is. However, the slope begins to increase at another inflexion point after control.

For the blue curve, the system is stable with enough control. The control causes a sudden reduction of the slope at the control time, and the slope keeps decreasing after the control with the mutation at the FEP. Apparently, the trajectory of the state-plane turns back at the unstable equilibrium point (UEP) by the minimum control quantity.

Figure 4. State-plane representation of speed and angle.

Figure 4. State-plane representation of speed and angle.

Figure 5. The slope over power angle.

Figure 5. The slope over power angle.

The derivative of Equation (2) with respect to t is given by:

$$\frac{dk(\mathrm{\delta },\mathrm{\omega })}{dt}=\frac{dk(\mathrm{\delta },\mathrm{\omega })}{d\mathrm{\delta }}\frac{d\mathrm{\delta }}{dt}=l \cdot \mathrm{\Delta }\mathrm{\omega }=\mathrm{\tau }$$

(12)

A similar equation can be found in the detection algorithm mentioned in the previous article. Obviously the decrease of the slope indicates that the trajectory runs in the concave area, and the trajectory runs in the convex area after the slope begins to increase. If the control is appropriate, the system will return to a stable state and the slope of the trajectory will keep decreasing until it reaches FEP.

3.2. Calculation of Control Quantity for SMIB

To find the relationship between the slope and angular speed, integrating both sides of Equation (2) gives:

$${\displaystyle \underset{{\mathrm{\delta }}_a}{\overset{{\mathrm{\delta }}_b}{\int }}k(\mathrm{\delta },\mathrm{\omega })}d\mathrm{\delta }={\displaystyle \underset{{\mathrm{\delta }}_a}{\overset{{\mathrm{\delta }}_b}{\int }}\frac{d\mathrm{\Delta }\mathrm{\omega }}{d\mathrm{\delta }}}d\mathrm{\delta }=\mathrm{\Delta }{\mathrm{\omega }}_b-\mathrm{\Delta }{\mathrm{\omega }}_a$$

(13)

where ${\mathrm{\delta }}_a$ and ${\mathrm{\delta }}_b$ are the bounds of the integration; $\mathrm{\Delta }{\mathrm{\omega }}_a$ and $\mathrm{\Delta }{\mathrm{\omega }}_b$ are the angular speed corresponding to ${\mathrm{\delta }}_a$ and ${\mathrm{\delta }}_b$. Let ${\mathrm{\delta }}_b$ be the power angle at the control time $t_a$; and then $\mathrm{\Delta }{\mathrm{\omega }}_a$ is the angular speed at $t_a$. For stable state, $\mathrm{\Delta }{\mathrm{\omega }}_b$ is zero when ${\mathrm{\delta }}_b$ is the FEP. For unstable state, if $\mathrm{\Delta }{\mathrm{\omega }}_b$ can be zero by the control, the system will return from ${\mathrm{\delta }}_b$ to a stable state and ${\mathrm{\delta }}_b$ will be the FEP. Therefore Equation (13) can be written as:

$${\displaystyle \underset{{\mathrm{\delta }}_a}{\overset{{\mathrm{\delta }}_b}{\int }}k(\mathrm{\delta },\mathrm{\omega })}d\mathrm{\delta }=-\mathrm{\Delta }{\mathrm{\omega }}_a$$

(14)

Therefore, the value of the slope at control time can be calculated by Equation (14), and the corresponding control quantity can be calculated by Equation (2). There are two unknowns ${\mathrm{\delta }}_b$ and $k(\mathrm{\delta },\mathrm{\omega })$ in Equation (14). $k(\mathrm{\delta },\mathrm{\omega })$ is related to the unbalanced power on the basis of Equation (2). ${\mathrm{\delta }}_b$ can be preset as the FEP as needed, but it must be no more than the UEP, so when ${\mathrm{\delta }}_b$ is acquired, $k(\mathrm{\delta },\mathrm{\omega })$ can be calculated. If ${\mathrm{\delta }}_b$ is just the UEP, the corresponding control quantity is the minimum.

3.3. Approximation Method to Calculate Control Quantity

As can be seen from the trajectories in Figure 3, $k(\mathrm{\delta },\mathrm{\omega })$ is nonlinear and changes corresponding to different control quantities at $t_a$. Because of the nonlinear nature of $k(\mathrm{\delta },\mathrm{\omega })$, the value of $k(\mathrm{\delta },\mathrm{\omega })$ at $t_a$ could not be acquired exactly. Therefore, an approximation method is presented to acquire the value of $k(\mathrm{\delta },\mathrm{\omega })$ at $t_a$. Employing a constant $k^{\prime }$ instead of $k(\mathrm{\delta },\mathrm{\omega })$ in Equation (14):

$${\displaystyle \underset{{\mathrm{\delta }}_a}{\overset{{\mathrm{\delta }}_b}{\int }}{k}^{\prime }}d\mathrm{\delta }=k^{\prime }({\mathrm{\delta }}_b-{\mathrm{\delta }}_a)=-\mathrm{\Delta }{\mathrm{\omega }}_a$$

(15)

$$k^{\prime }=\frac{-\mathrm{\Delta }{\mathrm{\omega }}_a}{{\mathrm{\delta }}_b-{\mathrm{\delta }}_a}$$

(16)

Obviously, the angular speed and power angle at control time and the power angle at FEP are necessary to calculate $k^{\prime }$. Above all ${\mathrm{\delta }}_b$ must be preset less than or equal to UEP, then the system can be stable after control. Therefore, $k(\mathrm{\delta },\mathrm{\omega })$ will keep decreasing after the control and the trajectory will return at ${\mathrm{\delta }}_b$. Because that $k(\mathrm{\delta },\mathrm{\omega })$ will keep decreasing after the control, $k^{\prime }$ is surely less than the value of $k(\mathrm{\delta },\mathrm{\omega })$ at control time. Hence, the control quantity calculated by $k^{\prime }$ is a little bigger than actually needed.

For a SMIB system, suppose the ratio of generator shedding is $\mathrm{\lambda }$, the relation of $\mathrm{\lambda }$ and $k^{\prime }$ after control is as follows:

$$k^{\prime }=\frac{(1-\mathrm{\lambda })P_m-P_{ea}}{(1-\mathrm{\lambda })M\mathrm{\Delta }{\mathrm{\omega }}_a}$$

(17)

$$\mathrm{\lambda }=1-\frac{P_{ea}}{P_m-M\mathrm{\Delta }{\mathrm{\omega }}_a{k}^{\prime }}$$

(18)

where $P_{ea}$ is the generator output electrical power at $t_a$. When $k^{\prime }$ is obtained by Equation (16), the generator shedding ratio of SMIB system can be calculated by Equation (18).

3.4. Seeking for FEP

The parameters needed in calculation can be collected by WAMS except for FEP. Two methods to acquire the FEP are provided in this paper:

Method 1: In order to obtain the minimum control quantity, it is intended to slow down the angular speed to zero at UEP by the control. It can be supposed that the power equilibrium point of the system after control is the UEP. Generator electrical power can be written as Equation (19), and mechanical power is considered as a constant in a short time period. It is assumed that generator shedding has a corresponding change on mechanical power:

$$P_e=P_c(t)+{\mathrm{\lambda }}_1(t)\sin (\mathrm{\delta })+{\mathrm{\lambda }}_2(t)\cos (\mathrm{\delta })$$

(19)

where $P_c(t),{\mathrm{\lambda }}_1(t),{\mathrm{\lambda }}_2(t)$ are parameters to be identified at $t$ instant. $P_c(t),{\mathrm{\lambda }}_1(t),{\mathrm{\lambda }}_2(t)$ are considered as constant in a short time when no other operations occur. At $t_a$ instant, $P_c(t_a),{\mathrm{\lambda }}_1(t_a),{\mathrm{\lambda }}_2(t_a)$ are calculated by the least square method, and then the electrical power can be acquired.

In order to seek for the UEP after the control, the control should be known at first. Therefore an iterative method can be employed as follows:

(1)

Obtain the prediction curve of the electrical power;

(2)

Preset zero to the generator shedding ratio and UPE: ${\mathrm{\lambda }}^{(0)}=0$, ${{\mathrm{\delta }}_u}^{(0)}=0$;

(3)

Mechanical power decreases at the same ratio:

$${P_m}^{(k)}=(1-{\mathrm{\lambda }}^{(k-1)}){P_m}^{(k-1)}$$

(4)

Search for the power equilibrium point as ${{\mathrm{\delta }}_b}^{(n)}$;

(5)

Acquire the ${\mathrm{\lambda }}^{(k)}$ according to the ${{\mathrm{\delta }}_b}^{(n)}$;

(6)

If $\left|{{\mathrm{\delta }}_b}^{(n)}-{{\mathrm{\delta }}_b}^{(n-1)}\right|\le \mathrm{\epsilon }$, complete the iterator; or return to step 3.

where $\mathrm{\epsilon }$ is the convergence conditions; $n$ is the number of the iterator.

Method 2: Preset a FEP as needed. On the basis of the operating requirement, the FEP, which must be less than or equal to the power equilibrium point, can be preset within the limits as needed.

Two control objectives are provided according to the way of acquiring the FEP: scheme one is the minimal cost to help restore the system to a stable state; scheme two is the minimal cost to limit the maximum angle.

3.5. Control Method for Multi-Machine System

Using the assumption made by [22] that the disturbed multi-machine system separates into two groups, leading group S and lagging group A, the partial center of angles, angular speed, mechanical power and electrical power of the two-machine power system are shown as follows:

$${\mathrm{\delta }}_s=\frac{{\displaystyle \sum _{i\in S}{M}_i{\mathrm{\delta }}_i}}{{\displaystyle \sum _{i\in S}{M}_i}}{\mathrm{\delta }}_a=\frac{{\displaystyle \sum _{i\in A}{M}_i{\mathrm{\delta }}_i}}{{\displaystyle \sum _{i\in A}{M}_i}}$$

(20)

$$\mathrm{\Delta }{\mathrm{\omega }}_s=\frac{{\displaystyle \sum _{i\in S}{M}_i\mathrm{\Delta }{\mathrm{\omega }}_i}}{{\displaystyle \sum _{i\in S}{M}_i}}\mathrm{\Delta }{\mathrm{\omega }}_a=\frac{{\displaystyle \sum _{i\in A}{M}_i\mathrm{\Delta }{\mathrm{\omega }}_i}}{{\displaystyle \sum _{i\in A}{M}_i}}$$

(21)

$$\begin{array}{l}{P}_{ms}={\displaystyle \sum _{i\in S}{P}_{mi}}{P}_{ma}={\displaystyle \sum _{i\in A}{P}_{mi}}\\ P_{es}={\displaystyle \sum _{i\in S}{P}_{ei}}{P}_{ea}={\displaystyle \sum _{i\in A}{P}_{ei}}\end{array}$$

(22)

Further, the two-machine power system can be equivalent to the SMIB system:

$$\mathrm{\delta }={\mathrm{\delta }}_s-{\mathrm{\delta }}_a$$

(23)

$$\mathrm{\Delta }\mathrm{\omega }=\mathrm{\Delta }{\mathrm{\omega }}_s-\mathrm{\Delta }{\mathrm{\omega }}_a$$

(24)

$$\mathrm{\Delta }P=P_M-P_E$$

(25)

where

$$\begin{array}{l}{M}_T={\displaystyle \sum _{i=1}^n{M}_i},M=\frac{M_s{M}_a}{M_T},\\ P_M=\frac{M_a{P}_{ms}-M_s{P}_{ma}}{M_T},\\ P_E=\frac{M_a{P}_{es}-M_s{P}_{ea}}{M_T}.\end{array}$$

According to the calculation method applied to SMIB equivalent system, the expected slope $k^{\prime }$ of SMIB equivalent system can be acquired. Transform the expected slope to control quantity:

$$\mathrm{\Delta }{P}_{ms}=M_s(k-k^{\prime }){\mathrm{\omega }}_0(\mathrm{\Delta }{\mathrm{\omega }}_s-\mathrm{\Delta }{\mathrm{\omega }}_a)$$

(26)

where $k$ is the original slope before control.

4. Simulation Result

4.1. SMIB System

The SMIB system as shown in Figure 3 is employed to illustrate the algorithm. Under normal conditions, the system operates at a stable equilibrium point. A three-phase grounding fault occurs on line L2, and the fault is cleared by switching off the line. The clear time differ from 0.17 s (critical clearing time) to 0.22 s. The system parameters are shown in Table 1.

Table 1. The parameters of the SMIB system.

Component	Parameters
The initial state	${\mathrm{\omega }}_0=2\mathrm{\pi }f$, $f=50\text{Hz}$, ${\mathrm{\delta }}_0={34.49}^{\circ }$
The generator	$P_m=120\text{MW}$, $T_j=6\text{s}$, ${\text{x}}_{\text{d}}=1.83\text{p.u.}$, ${\text{x}}_{\text{q}}=1.83\text{p.u.}$, ${{\text{x}}_{\text{d}}}^{\prime }=0.3\text{p.u.}$, ${{\text{x}}_{\text{q}}}^{\prime }=1.83\text{p.u.}$, ${x_q}^{″}=0.25\text{p.u.}$, ${x_d}^{″}=0.25\text{p.u.}$
The transmission line	$x_1=x_2=0.486\text{Ω/km},x_0=4x_1$

Table 1. The parameters of the SMIB system.

Component	Parameters
The initial state	${\mathrm{\omega }}_0=2\mathrm{\pi }f$, $f=50\text{Hz}$, ${\mathrm{\delta }}_0={34.49}^{\circ }$
The generator	$P_m=120\text{MW}$, $T_j=6\text{s}$, ${\text{x}}_{\text{d}}=1.83\text{p.u.}$, ${\text{x}}_{\text{q}}=1.83\text{p.u.}$, ${{\text{x}}_{\text{d}}}^{\prime }=0.3\text{p.u.}$, ${{\text{x}}_{\text{q}}}^{\prime }=1.83\text{p.u.}$, ${x_q}^{″}=0.25\text{p.u.}$, ${x_d}^{″}=0.25\text{p.u.}$
The transmission line	$x_1=x_2=0.486\text{Ω/km},x_0=4x_1$

It is needed to point out that the simulation just acquires the power angle, angular speed, mechanical power, electrical power, inertia constant of the generator without other information, which can be acquired on-line. When the sampling time interval is 10 ms, the detection results based on the instability criterion mentioned above are shown in Table 2.

Table 2. The transient stability detection results by 100 Hz sample frequency.

The Moment of Fault-Clearing	Simulation Results	Detection Results	The Moment of Instability Detected	The Angle of Instability Detected
0.17 s	Stable	Stable	\	\
0.18 s	Stable	Stable	\	\
0.187 s	Stable	Stable	\	\
0.188 s	Unstable	Unstable	0.56 s	121.9°
0.19 s	Unstable	Unstable	0.46 s	114.9°
0.20 s	Unstable	Unstable	0.37 s	106.2°
0.21 s	Unstable	Unstable	0.34 s	102.9°
0.22 s	Unstable	Unstable	0.33 s	102.8°

Table 2. The transient stability detection results by 100 Hz sample frequency.

The Moment of Fault-Clearing	Simulation Results	Detection Results	The Moment of Instability Detected	The Angle of Instability Detected
0.17 s	Stable	Stable	\	\
0.18 s	Stable	Stable	\	\
0.187 s	Stable	Stable	\	\
0.188 s	Unstable	Unstable	0.56 s	121.9°
0.19 s	Unstable	Unstable	0.46 s	114.9°
0.20 s	Unstable	Unstable	0.37 s	106.2°
0.21 s	Unstable	Unstable	0.34 s	102.9°
0.22 s	Unstable	Unstable	0.33 s	102.8°

When the sampling time interval is 20 ms, the detection results based on the instability criterion mentioned above are shown in Table 3.

Table 3. The transient stability detection results by 50 Hz sample frequency.

The Moment of Fault-Clearing	Simulation Results	The Moment of Instability Detected	The Angle of Instability Detected	The Moment of Angle Reach the Threshold of 180°
0.17 s	Stable	Stable	-	-
0.18 s	Stable	Stable	-	-
0.187 s	Stable	Stable	-	-
0.188 s	Unstable	0.56 s	121.9°	1.44 s
0.19 s	Unstable	0.46 s	114.9°	1.14 s
0.20 s	Unstable	0.38 s	108.5°	0.86 s
0.21 s	Unstable	0.34 s	102.9°	0.76 s
0.22 s	Unstable	0.34 s	105.6°	0.68 s

Table 3. The transient stability detection results by 50 Hz sample frequency.

The Moment of Fault-Clearing	Simulation Results	The Moment of Instability Detected	The Angle of Instability Detected	The Moment of Angle Reach the Threshold of 180°
0.17 s	Stable	Stable	-	-
0.18 s	Stable	Stable	-	-
0.187 s	Stable	Stable	-	-
0.188 s	Unstable	0.56 s	121.9°	1.44 s
0.19 s	Unstable	0.46 s	114.9°	1.14 s
0.20 s	Unstable	0.38 s	108.5°	0.86 s
0.21 s	Unstable	0.34 s	102.9°	0.76 s
0.22 s	Unstable	0.34 s	105.6°	0.68 s

The results in Table 2 and Table 3 show that the instability criterion mentioned in this paper can distinguish unstable cases and stable cases correctly and rapidly. Even the critical cases of which the difference of fault-clearing time is only 1 ms can be rightly detected. The accuracy of the detection result is not influenced by the sample frequency. From the results in Table 3, the angles of instability detected vary from 102.9° to 121.9°, which are much less than 180°. Further, the detection moment of the method proposed in this paper is much earlier than it of the threshold of 180°, which can provide adequate time for implementation of control measures.

It is necessary to take some control actions when the system is detected as unstable. Two control objectives are provided by the way to acquire the FEP: scheme one is the minimal cost to help the system restore to stable state; scheme two is the minimal cost to limit the maximum angle.

When the system is unstable, the control quantity is calculated by scheme one. The results are shown in Table 4. The minimum control quantity calculated by scheme one is close to and a little bigger than the real minimum control quantity. The control quantity calculated by this paper is effective and conservative as expected. In addition, the conservative algorithm is helpful for the recovery of the power system.

The same condition introduced above, control quantity is calculated by scheme two, and the preset FEP is 130°.The result is shown in Table 5. The FEP after control is closed to and less than 130° as expected. It means that the scheme two intended to limit the FEP of the power system is effective. Transient stability margin is directly reflected in the angle swing range. The method can be applied to improve the transient stability margin when the power system suffers large disturbance.

Table 4. Calculated minimum control compare to the real minimum control.

Fault-Clearing Time	Distinction Time	Calculated Minimum Control Quantity (%)	Real Minimum Control Quantity (%)
0.18 s	None (stable)	None (stable)	None (stable)
0.19 s	0.46 s	5.61	4.55
0.20 s	0.37 s	14.1	13.5
0.21 s	0.34 s	22.4	20.6
0.22 s	0.33 s	30.1	28.4

Table 4. Calculated minimum control compare to the real minimum control.

Fault-Clearing Time	Distinction Time	Calculated Minimum Control Quantity (%)	Real Minimum Control Quantity (%)
0.18 s	None (stable)	None (stable)	None (stable)
0.19 s	0.46 s	5.61	4.55
0.20 s	0.37 s	14.1	13.5
0.21 s	0.34 s	22.4	20.6
0.22 s	0.33 s	30.1	28.4

Table 5. Calculated minimum control by limiting the maximum angle.

Fault-Clearing Time	Distinction Time	Calculated Controlled Quantity (%)	Real FEP after Control
0.18 s	None (stable)	None (stable)	None (stable)
0.19 s	0.46 s	8.91	129.1°
0.20 s	0.37 s	18.6	126.8°
0.21 s	0.34 s	22.5	128.1°
0.22 s	0.33 s	35.6	128.8°

Table 5. Calculated minimum control by limiting the maximum angle.

Fault-Clearing Time	Distinction Time	Calculated Controlled Quantity (%)	Real FEP after Control
0.18 s	None (stable)	None (stable)	None (stable)
0.19 s	0.46 s	8.91	129.1°
0.20 s	0.37 s	18.6	126.8°
0.21 s	0.34 s	22.5	128.1°
0.22 s	0.33 s	35.6	128.8°

4.2. IEEE 39-Bus System

The IEEE 39-bus system shown in Figure 6 is employed to illustrate the algorithm. Two control objectives are provided according to the way of acquiring the FEP: scheme one is the minimal cost to help the system restore to stable state; scheme two is the minimal cost to limit the maximum angle.

Figure 6. The diagram of IEEE 10-unit 39-bus power system.

Figure 6. The diagram of IEEE 10-unit 39-bus power system.

A three-phase grounding fault occurs on the transmission line between bus 4 and bus 14 and the fault duration is 0.4 s. The power angle curves are shown in Figure 7.

Figure 7. Multi-machine swing curves after fault.

Figure 7. Multi-machine swing curves after fault.

The results of the closed loop control system are shown in Table 6.

Table 6. Results of the control system.

Fault	A Three-Phase Grounding Fault Occurs on the Line between Bus 4 and Bus 14
Fault duration	0.23 s
Detection time	0.41 s
Detection angle	126.5904°
Control objective	scheme one	scheme two (145°)
Control law (MW)	G31 (439)	G31 (521)
FEP after control	148.0°	144.8°

Table 6. Results of the control system.

Fault	A Three-Phase Grounding Fault Occurs on the Line between Bus 4 and Bus 14
Fault duration	0.23 s
Detection time	0.41 s
Detection angle	126.5904°
Control objective	scheme one	scheme two (145°)
Control law (MW)	G31 (439)	G31 (521)
FEP after control	148.0°	144.8°

The power angle curves after scheme one are shown in Figure 8. It can be seen that both scheme one and scheme two quickly dampen the system oscillation and keep the system stable. The cost of scheme one is 439 MW, and the FEP after control is 148.0°. The cost of scheme two is 521 MW, and the FEP is 144.8° within the limit.

Figure 8. Power angle curves of scheme one.

Figure 8. Power angle curves of scheme one.

5. Conclusions

In conclusion, the detection method and the control method deduced in this paper just use real-time information rather than the information of system models, parameters and disturbances to deal with the transient instability problem. As a result, the control system could be applicable in real time. The detection technique can identify instability accurately and quickly. Based on the tangent slope of the state plane trajectory, a new control method proposed in this paper forms a transient stability control system with the detection method. Two control objectives which are successfully tested in the simulations are provided based on the method of acquiring FEP in this paper. The simulation results in SMIB system and IEEE 39-bus system indicate the effectiveness of the proposed transient stability control system.