PSO-Based Oscillatory Stability Assessment by Using the Torque Coefficients for SMIB

Abstract

This study discusses the evaluation of oscillatory stability based on the synchronizing $K_s$ and damping $K_d$ torque coefficients for a single-machine system connected to an infinite bus (SMIB). Particle swarm optimization (PSO) technique is used to determine $K_s$ and $K_d$ values and subsequently identify the oscillatory stability conditions in the SMIB. The ability of PSO is compared with those of evolutionary programming (EP) techniques and artificial immune system (AIS). The least square (LS) method is selected as the benchmark for $K_s$ and $K_d$ values determined by PSO, EP, and AIS. Simulation results show that PSO successfully estimated $K_s$ and $K_d$ values closest to LS compared with EP and AIS. PSO also uses lower computational time compared with those of the two other techniques. The proposed technique is suitable for solving oscillatory stability problems based on the determination of eigenvalues and minimum damping ratio.

1. Introduction

The increase in population worldwide has major implications to electricity consumption. Power generation companies always innovate to ensure that the energy supply is sufficient and stable. Accordingly, small signal stability is a frequently discussed topic, as reported in [1,2,3,4,5]. This system must be tracked online because its power system operating condition changes over time. Various indicators, such as damping ratio [6,7,8] and damping factor [9], have been proposed to determine the angular stability of a system. However, the eigenvalues obtained from the entire mathematical model of the system are needed to calculate these two indicators. In this study, the synchronizing $K_s$ and damping $K_d$ torque coefficients were introduced to verify system stability. $K_s$ and $K_d$ can be calculated on the basis of the information from three rotor responses, namely the changes in rotor angle, $\Delta \delta \left(t\right)$; rotor speed, $\Delta \omega \left(t\right)$; and electromechanical torque, $\Delta T_e\left(t\right)$. A system is considered stable when the values of $K_s$ and $K_d$ are positive [10,11,12]. However, the system is considered unstable when one of the values is negative.

Least square (LS) method is often used to calculate $K_s$ and $K_d$ values of a system by using a static parameter estimation approach [13,14,15]. However, the LS method needs a large data set to obtain the correct values. Extensive data collection also requires long computation time. Monitoring the oscillation period is also critical because the occurrence of slight data error leads to inaccurate $K_s$ and $K_d$ values. Thus, heuristic techniques are introduced to solve this problem. $K_s$ and $K_d$ values will be optimized from the beginning of the data until constant $K_s$ and $K_d$ values are obtained. Therefore, only a portion of the data is necessary for estimating $K_s$ and $K_d$. Minor data errors do not significantly affect the determination of $K_s$ and $K_d$ values.

In recent years, the use of artificial intelligence (AI) technology has become the preferred option in solving power system problems. The use of AI is introduced to solve the optimum values of a system or condition particularly in the fields of economic dispatch, capacitor placement and sizing, and assessment and improvement of voltage and oscillatory stability. Artificial neural networks [16,17], evolutionary programming (EP) [18,19,20], artificial immune systems (AIS) [21,22,23], and ant colony optimization (ACO) [24,25,26] are AI approaches that are commonly used in power systems. The EP algorithm is modeled on the biological evolution process of solving a complex problem. The main features of EP include the mutation process of the next generation and the selection of increasingly powerful genes. The AIS algorithm uses a concept similar to that of EP. Although both concepts are biologically based on living things, EP focuses on the evolution of living things, whereas AIS adopts the concept of the living immune system. The difference between AIS and EP algorithms is that AIS has an additional process of cloning called the clonal selection algorithm. However, the ACO approach is inspired by the true behavior of ants while searching for food and interacting with fellow ants. In ACO, artificial ants (the search agent) will communicate by using pheromones, which guide the searcher ants to solve the calculation problem by tracking the best route. Meanwhile, the particle swarm optimization (PSO) [27,28,29,30] concept mimics the movements of a herd, such as the behavior of schooling fish and swarming insects. This technique was originally founded based on the population of random particles, in which every particle is a potential solution. PSO can make adjustments to obtain balance between global and local explorations during the search process. This feature makes the PSO suitable for overcoming the problems caused by initial convergence and improving the ability to search.

This study presents the techniques for determining oscillatory stability based on the estimation of twin indicators called synchronizing $K_s$ and damping $K_d$ torque coefficients. A single machine linked to a large bus network or infinite bus (SMIB) is selected as the test system. The changes in rotor angle, $\Delta \delta \left(t\right)$; rotor speed, $\Delta \omega \left(t\right)$; and electromechanical torque, $\Delta T_e\left(t\right)$ are used to determine $K_s$ and $K_d$ values. This optimization process aims to minimize the error of both torque coefficients. In this study, the PSO technique was selected as the heuristic technique for solving this optimization problem. The results in PSO will be compared with those of EP and AIS. From the simulation using MATLAB, these three heuristic techniques will be compared based on the accuracy of the torque coefficient, the amount of iteration for the simulation process, and simulation time. The eigenvalues $\lambda$ and minimum damping ratio ${\xi }_{min}$, are also used to verify the system stability.

2. Materials and Methods

2.1. Oscillatory Stability Assessment

In oscillatory instability detection studies, early detection can provide the system time to improve its stability and prevent further system instability and eventual collapse. This study introduces $K_s$ and $K_d$ as the indicators for determining system stability. $K_s$ and $K_d$ are estimated periodically whenever new data are received from the system. To validate the result, $K_s$ and $K_d$ are estimated using the LS calculation method. Eigenvalue and damping ratio analysis were also conducted to verify the accuracy of the results.

2.1.1. $K_s$ and $K_d$

In the oscillatory stability analysis, electromechanical torque can be expressed in the form of synchronous and damping torque components. The relationship among the rates of change in the estimated electromechanical torque, $\Delta T_{es}\left(t\right)$; rotor angle, $\Delta \delta \left(t\right)$; in rotor speed, $\Delta \omega \left(t\right)$; synchronizing torque coefficients, $K_s$; and damping torque coefficients, $K_d$ can be expressed as follows:

$$\Delta T_{es}\left(t\right)=K_s·\Delta \delta \left(t\right)+K_d·\Delta \omega \left(t\right)$$

(1)

The stability evaluation of a linear system can be predicted based on $K_s$ and $K_d$ values. A stable system is guaranteed when $K_s$ and $K_d$ values are positive. Figure 1 illustrates a stable angle stability as resulted from the positive values of $K_s$ and $K_d$.

If the linear system has a positive $K_s$ and a negative $K_d$ value, then the system is under an unstable oscillatory condition, which is due to lack of adequate damping torque. The effect of the oscillatory instability condition can be detected from the increment of amplitude oscillations of the rotor. Figure 2 shows the unstable conditions of angle stability from the positive $K_s$ value and the negative $K_d$ value.

Non-oscillatory instability occurs when $K_s$ and $K_d$ show negative and positive values, respectively, because the absence of automatic voltage regulators results in the lack of sufficient synchronizing torque. This condition can be verified based on the steady increment of rotor angle response. Figure 3 presents the unstable conditions for angle stability from the negative $K_s$ and positive $K_d$ value.

2.1.2. LS Method

Several methods have been proposed to estimate the value of $K_S$ and $K_d$ which involved optimization approach. LS method can be one of the possible techniques in addressing this phenomenon, which has been used as static parameter estimation as reported in [18]. The LS method is a form of mathematical analysis that calculates the most appropriate solution based on a data set. In this study, the LS technique is used to minimize the sum of the squares of errors $e\left(t\right)$. Error $e\left(t\right)$ is the difference between the changes in the estimated electromechanical torque $\Delta T_{es}\left(t\right)$ and the changes in the electromechanical torque $\Delta T_e\left(t\right)$. From the previous equation, $\Delta T_{es}$ can be calculated based on the values of $\Delta \delta \left(t\right)$ and $\Delta \omega \left(t\right)$ from the data set, and also estimated values of $K_s$ and $K_d$. The equation of $e\left(t\right)$ is as follows:

$$e\left(t\right)=\Delta T_{es}\left(t\right)-\Delta T_e\left(t\right)=\mathsf{\Lambda }x-\Delta T_e\left(t\right)$$

(2)

$$\mathsf{\Lambda }=\left[\begin{array}{cc}\Delta \delta \left(t\right)& \Delta \omega \left(t\right)\end{array}\right]$$

(3)

$$x={\left[\begin{array}{cc}{K}_s& K_d\end{array}\right]}^T$$

(4)

For the overdetermined matrix system $\mathsf{\Lambda }x=\Delta T_{es}\left(t\right)$, the sum of the squared value of $e\left(t\right)$ can be defined as the following function $f\left(x\right)$:

$$f\left(x\right)=\sum _{t=1}^{N}e{\left(t\right)}^2=\sum _{t=1}^N{\Vert \mathsf{\Lambda }x-\Delta T_e\left(t\right)\Vert }^2$$

(5)

where N is the number of experimental samples.

To minimize $e\left(t\right)$, $K_s$ and $K_d$ need to be calculated throughout the entire oscillation period. The solution for this equation can be given as

$${\mathsf{\Lambda }}^T·\mathsf{\Lambda }·x={\mathsf{\Lambda }}^T·\Delta T_e\left(t\right)$$

(6)

Matrix ${\mathsf{\Lambda }}^T·\mathsf{\Lambda }$ is invertible. Thus, the solution of x is given by

$$x={\left[{\mathsf{\Lambda }}^T·\mathsf{\Lambda }\right]}^{-1}·{\mathsf{\Lambda }}^T·\Delta T_e\left(t\right)$$

(7)

The torque coefficients $K_s$ and $K_d$ can be calculated by solving Equation (7). Although the calculated values are accurate, the application of the LS method requires data for the entire swing [5,6]. The need for such large data collection also requires a long calculation time. Accordingly, a new indicator is needed.

2.1.3. Eigenvalue and Damping Ratio

There are various methods for determining the stability of a system, including the calculation of eigenvalues and damping ratios. Eigenvalues are derived from matrix arrays by linear system equations. However, the damping ratio is derived using a combination of real and imaginary values of eigenvalues. Eigenvalues and damping ratios are often used as indicators for measuring the stability of a system, as described in [18]. Eigenvalues $\lambda$ of a matrix representing a linear system are obtained as follows:

$$\mathsf{\Pi }\varphi =\lambda \varphi$$

(8)

where $\mathsf{\Pi }$ is a $(n\times n)$ matrix, $\varphi$ is a $(n\times 1)$ vector, and $\lambda$ is the eigenvalues matrix $(n\times n)$. To determine the eigenvalues, Equation (8) can be written as

$$(\mathsf{\Pi }-\lambda I)\varphi =0$$

(9)

The eigenvalues of $\mathsf{\Pi }$ are the n solutions of $\lambda ={\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$. The ith eigenvalue of a matrix representing a linear system are obtained as follows:

$${\lambda }_i={\sigma }_i\pm j{\omega }_i$$

(10)

where ${\sigma }_i$ and ${\omega }_i$ are the real and imaginary parts of the ith eigenvalue, respectively. The negative real part of all eigenvalues indicates that a linear system is stable. The damping ratio, ${\xi }_i$, for the ith eigenvalue is defined as [31,32]:

$${\xi }_i=-\frac{{\sigma }_i}{\sqrt{{\sigma }_i^2+{\omega }_i^2}}$$

(11)

The linear system is certainly stable when all the damping ratios have positive values. For simplification, only the minimum damping ratio, ${\xi }_{min}$, for the linear system is selected to verify the result.

2.2. Test System

In this study, an SMIB system is selected for evaluation. $\Delta \delta \left(t\right)$, $\Delta \omega \left(t\right)$, and $\Delta T_e\left(t\right)$ are generated in a MATLAB Simulink environment based on the block diagram of the SMIB model. According to Equation (1), the estimated electric torque, $\Delta T_{es}\left(t\right)$, is calculated by using the generated sample data for $K_s$ and $K_d$.

2.2.1. Philips–Heffron Model

Figure 4 shows the SMIB system with a classical generator model.

Here, $R_e$ and $X_e$ are resistance and reactance of the transmission line, respectively. $V_T$ and $V_{\infty }$ are terminal bus and infinite bus voltages, respectively. $\delta$ and $\alpha$ are the rotor angles at terminal bus and infinite bus, respectively. $E_q^{\prime }$ is the generator terminal voltage. The equations for rotor acceleration, rotor angle, and the field circuit are presented as follows [33,34]:

$$\frac{\Delta {\omega }_r}{\Delta t}=\frac{\Delta T_m-K_1\Delta \delta -D\Delta {\omega }_r-K_2\Delta {\psi }_{fd}}{2H}$$

(12)

$$\frac{\Delta \delta }{\Delta t}={\omega }_0\Delta {\omega }_r$$

(13)

$$\frac{\Delta {\psi }_{fd}}{\Delta t}=\frac{-(K_3{K}_4\Delta \delta +\Delta {\psi }_{fd}-K_3\Delta E_{fd})}{T_3}$$

(14)

Here, $\Delta {\omega }_r/\Delta t$ is the change of rotor acceleration, $\Delta \delta /\Delta t$ is the change of rotor angle, $\Delta {\psi }_{fd}/\Delta t$ is the change of the field circuit, $E_{fd}$ is a field voltage, $T_m$ is the rotor mechanical torque, H is the machine inertia constant, and D is the machine damping coefficients.

According to Equations (12)–(14), the Philips–Heffron block diagram model of SMIB is developed and shown in Figure 5 [1]. This figure illustrates that the constant K represents several variables, such as the electric torque, rotor speed, and rotor angle. From the Philips–Heffron block diagram model, the following are developed:

$$\Delta \dot{X}=A·\Delta X+B·\Delta U$$

(15)

$$A=\left[\begin{array}{ccc}-\frac{D}{2H}& -\frac{K_1}{2H}& -\frac{K_2}{2H}\\ {\omega }_0& 0& 0\\ 0& -\frac{{\omega }_0{R}_{fd}{m}_1{L}_{ads}^{\prime }}{L_{fd}}& -\frac{{\omega }_0{R}_{fd}}{L_{fd}}\left(1-\frac{L_{ads}^{\prime }}{L_{fd}}+m_2{L}_{ads}^{\prime }\right)\end{array}\right]$$

(16)

$$\Delta X=\left[\begin{array}{ccc}\Delta {\omega }_r& \Delta \delta & \Delta {\psi }_{fd}\end{array}\right]$$

(17)

$$B=\left[\begin{array}{cc}\frac{1}{2H}& 0\\ 0& 0\\ 0& \frac{{\omega }_0{R}_{fd}}{L_{adu}}\end{array}\right]$$

(18)

$$\Delta U=\left[\begin{array}{cc}\Delta T_m& \Delta E_{fd}\end{array}\right]$$

(19)

$$m_1=\frac{E_B(X_{Tq}sin{\delta }_0-R_{T}cos{\delta }_0)}{D}$$

(20)

$$m_2=\frac{X_{Tq}}{D}\frac{L_{ads}}{L_{ads}+L_{fd}}$$

(21)

Here, $R_{fd}$ and $L_{fd}$ are resistance and reactance of field circuit, respectively. $L_{ads}$ and $L_{aqs}$ are the generator d-axis and q-axis saturated value of mutual inductance, respectively. $X_{Tq}$ is total q-axis reactance of the system. $R_T$ is total system resistance. ${\delta }_0$ is the initial rotor angle. $L_{fd}$ is the field circuit reactance.

The function of the system parameters is separated into two system matrices, namely A and B. The details in the calculation of matrices A and B are presented in [18].

2.2.2. Objective Function

In this study, three optimization techniques, namely EP, AIS, and PSO, are used to determine the values of torque coefficients $K_s$ and $K_d$. The difference between $\Delta T_{es}\left(t\right)$ and $\Delta T_e\left(t\right)$ is defined as an error. This error is minimized at each iteration by using heuristic techniques when the new $K_s$ and $K_d$ values are calculated. In this study, the objective function is formulated based on the error as follows:

$$J=inv\left(1+\left|\frac{\Delta T_e\left(t\right)-\Delta T_{es}\left(t\right)}{\Delta T_e\left(t\right)}\right|\right)$$

(22)

By the implementation of this objective function will ensures that the difference between $\Delta T_{es}\left(t\right)$ and $\Delta T_e\left(t\right)$ versus $\Delta T_e\left(t\right)$ is close to or equal to 0. This objective function is developed such that the value of J will be in the range 0 and 1, and value 1 is the optimal value. The designed problem in this study can be represented by

$$Maximize\left(J\right)$$

(23)

These two torque coefficients are simultaneously optimized via the three optimization techniques for the different cases by using the proposed objective function.

2.2.3. Algorithm for $K_s$ and $K_d$ Estimation

The flowchart for $K_s$ and $K_d$ estimation process is shown in Figure 6.

3. Computational Intelligence Methods for Oscillatory Stability Assessment

In recent years, AI technology has been widely used in solving optimization problems in power systems. Evolutionary computation (EC) is a widely used AI technique. EC models evolutionary processes to develop strategies that seek optimal or almost optimal solutions for specific problems. Examples of EC techniques include EP, genetic algorithm, AIS, and PSO. In this study, AIS, EP, and PSO are selected as optimization techniques.

3.1. PSO

The PSO technique is inspired by the feeding process of certain animals, such as swarming birds and schooling fish. The PSO algorithm begins with initialization, followed by the update of velocity and position, fitness calculation, the best position update, and convergence test. The pseudocode that represents the PSO algorithm is illustrated in Algorithm 1. Detailed explanations of the PSO algorithm process can be found in [28].

In Algorithm 1, k is the number of iterations, i is the number of particles, $\omega$ is the inertia weight, $v_i$ and $x_i$ are the velocity and position for the ith particle, respectively. $c_1$ and $c_2$ are the acceleration coefficients, J is the objective function, as shown in Equation (22), $x_{b,i}$ is the personal best position for the ith particle, and $x_g$ is the global best position.

Algorithm 1 Pseudocode for the PSO algorithm [28].

initialize particle fork = 1 : maximum iteration do     for i = 1 : number of particle do         $v_i\left(k\right)=w·v_i(k-1)+c_1·x_{b,i}-x_i(k-1)+c_2·x_g-x_i(k-1)$         $x_i\left(k\right)=v_i\left(k\right)+x_i(k-1)$         calculate $J\left(x_i\left(k\right)\right)$     end for     if $J\left(x_i\left(k\right)\right)>J\left(x_{b,i}\right)$ then         $x_{b,i}=x_i\left(k\right)$     end if     rank $x\left(k\right)$ in ascending order of $J\left(x\right(k\left)\right)$     if $J{\left(x\left(k\right)\right)}_{max}>J\left(x_g\right)$ then         $J\left(x_g\right)=J{\left(x\left(k\right)\right)}_{max}$         define new $x_g$     end if     if ${|J\left(x\left(k\right)\right)}_{max}-J{\left(x\left(k\right)\right)}_{min}|<{10}^{-5}$ then         return     end if     $k=k+1$ end for

3.2. EP

The concept of EP is based on the theory of evolution of life through natural selection. EP is motivated by the process at the evolution stage (parents, mutation, offspring) but without the genetic evolution. The EP algorithm begins with the initialization, followed by the determination of fitness, mutations, combinations of parent and offspring, and ends with selection. Algorithm 2 demonstrates the pseudocode for the EP algorithm [18].

Algorithm 2 Pseudocode for the EP algorithm [18].

initialize population forl = 1 : maximum generation do     ///// Parents     for j = 1 : number of population do         define $x_{j,par}\left(l\right)$ and $J\left(x_{j,par}\left(l\right)\right)$     end for     ///// Offspring     for j = 1 : number of population do         $x_{j,off}\left(l\right)=\alpha ·(x_{j,par}{\left(l\right)}_{max}-x_{j,par}{\left(l\right)}_{min})·x_{j,par}\left(l\right)/J{\left(x_{j,par}\left(l\right)\right)}_{max}$         calculate $J\left(x_{j,off}\left(l\right)\right)$     end for     combine parents and offspring     rank $x\left(l\right)$ in ascending order of $J\left(x\right(l\left)\right)$     select top-half of $x\left(l\right)$ as new $x_{j,par}\left(l\right)$     if ${|J\left(x\left(l\right)\right)}_{max}-J{\left(x\left(l\right)\right)}_{min}|<{10}^{-5}$ then         return     end if     $l=l+1$ end for

In Algorithm 2, l is the number of generations, j is the number of populations, $\alpha$ is a mutation factor in EP, and $x_{i,par}$ and $x_{i,off}$ are the parents and offspring for the $jth$ population, respectively.

3.3. AIS

AIS is an optimization technique that attempts to biologically imitate the human immune system. This concept, which is practiced in the AIS technique, is similar to that of the EP technique. However, AIS has an additional criterion called cloning. The entire process is given in the form of a pseudocode in Algorithm 3 [23].

In Algorithm 3, m is the number of cycles, h is the number of cells, $\beta$ is a mutation factor in AIS, $x_h$ is the precloning of hth cells, $x_c$ is the post-cloning cells, and $x_{mut,h}$ is the mutated hth cells.

Table 1. List of parameters used in PSO, EP and AIS.

Techniques	Parameters
PSO	$c_1$ = 0.5, $c_2$ = 0.5, $\omega$ = 0.05
EP	$\alpha$ = 0.05
AIS	$\beta$ = 0.05

lists the parameters used in the PSO, EP, and AIS techniques. The value of these parameters are selected based on the values used in the reference [8].

Algorithm 3 Pseudocode for the AIS algorithm [23].

initialize cells form = 1 : maximum cycle do     for h = 1 : number of cells do         define $x_h\left(m\right)$ and $J\left(x_h\left(m\right)\right)$     end for     ///// Cloning     rank $x\left(m\right)$ according to $J\left(x\right(m\left)\right)$     select top-half of $x\left(m\right)$     clone $x\left(m\right)$ to become $x_c\left(m\right)$     clone $J\left(x\right(m\left)\right)$ to become $J\left(x_c\left(m\right)\right)$     ///// Mutate     for h = 1 : number of cells do         $x_{mut,h}\left(m\right)=\beta ·(x_c{\left(m\right)}_{max}-x_c{\left(m\right)}_{min})·x_c\left(m\right)/J{\left(x_c\left(m\right)\right)}_{max}$         calculate $J\left(x_{mut,h}\left(m\right)\right)$     end for     rank $x_{mut}\left(m\right)$ in ascending order of $J\left(x_{mut}\left(m\right)\right)$     select $x_{mut}\left(m\right)$ as new $x_h$     if then $|J{\left(x_{mut}\left(m\right)\right)}_{max}-J{\left(x_{mut}\left(m\right)\right)}_{min}|<{10}^{-5}$         return     end if     $m=m+1$ end for

4. Results and Discussion

This study estimates oscillatory stability by using $K_s$ and $K_d$ in various loading conditions of the SMIB system. The samples of $\Delta \delta \left(t\right)$, $\Delta \omega \left(t\right)$, and $\Delta T_e\left(t\right)$ are required in the calculated estimations generated in the MATLAB Simulink environment. EP, AIS, and PSO are used to estimate the values of $K_s$ and $K_d$. The results are compared with the benchmark values calculated via the LS method. The values of $\lambda$ and ${\xi }_{min}$ are also used to justify stability determination. For estimation with AI algorithms and LS, data size used is set to 20 s, while number of samples is set to 400 samples. The value of this 400 samples data proved to be effective for the AI algorithms and LS to make accurate calculations, based on the reference [18]. In this study, the simulation tests are conducted using Intel (R) Core (TM) i7-4770 CPU @ 3.40Ghz processor.

To determine the appropriate population size during the optimization process, the effect of population size for angle stability assessment of SMIB system has been studied. Three different population sizes, namely 5, 20 and 50 population sizes have been tested, using loading condition of P = 0.6 p.u. and Q = 0.7 p.u. The data of synchronizing torque coefficient, $K_s$, damping torque coefficient, $K_D$ and three different population sizes are tabulated in

Table 2. Results of $K_s$, $K_D$ and computation time for three different population sizes with P = 0.6 p.u. and Q = 0.7 p.u.

Population Sizes	${\mathit{K}}_{\mathit{s}}$, ${\mathit{K}}_{\mathit{d}}$	EP	AIS	PSO	LS
	$K_s$	0.6141	0.6153	0.6020	0.5928
5	$K_d$	1.1541	1.3063	0.3919	0.6282
	Computation Time (s)	5.7	4.3	5.1	-
	$K_s$	0.5942	0.5962	0.5926	0.5928
20	$K_d$	0.6564	0.7161	0.6014	0.6282
	Computation Time (s)	14.5	3.9	11.6	-
	$K_s$	0.5929	0.5954	0.5925	0.5928
50	$K_d$	0.6219	0.6958	0.6223	0.6282
	Computation Time (s)	24.5	8.2	17.9	-

.

From Table 2, it is clearly shown that the results obtained using PSO and LS are the same when optimized with 20 and 50 populations. The result for $K_S$ optimized using EP at population size of 50 is the same as that computed using LS. This indicates that 50 is the most suitable population size if closeness to LS technique is desired. This result is significant with the result of $K_D$. It is shown that PSO outperformed AIS and EP since it manages to achieve final solution close to the value obtained by LS. From the computation time point of view, it shows that computation time of optimization process of EP and PSO is proportional to the population size of the simulation. On the other hand, AIS method gives the shortest computation time when the population size is set to 20. This revealed that AIS managed to achieve optimal solution within population size of 20. It was also discovered that the value of fitness is always consistent even the size population is increased up to 50. According to this result, 20 was selected as population size during the optimization process.

In this study, two large events with different cases, namely Events A and B, are evaluated. In Event A, the value of active power, P, is randomly set to 0.5 p.u., whereas reactive power, Q, increases linearly by 0.1 p.u. from 0.1 p.u. to 1.0 p.u. In Event B, Q is set to 0.15 p.u., whereas P decreases by 0.1 p.u. from 1.0 p.u. to 0.1 p.u. The loading condition values are selected because they obtain significant results.

Table 3. $K_s$ and $K_d$ via EP, AIS, PSO, LS, eigenvalues and ${\xi }_{min}$ for Event A.

Cases	P, Q	${\mathit{K}}_{\mathit{s}}$, ${\mathit{K}}_{\mathit{d}}$	EP	AIS	PSO	LS	$\mathit{\lambda }$	${\mathit{\xi }}_{\mathit{min}}$
A-1	P = 0.5	$K_s$	0.8039	0.8243	0.8197	0.8197	−0.0698 ± j6.6437	0.0105
	Q = 0.1	$K_d$	1.2609	0.5720	0.9866	0.9791	−0.2610
A-2	P = 0.5	$K_s$	0.7959	0.7961	0.7931	0.7866	−0.0594 ± j6.5083	0.0091
	Q = 0.2	$K_d$	0.5350	0.5378	0.8483	0.8298	−0.2901
A-3	P = 0.5	$K_s$	0.7266	0.7266	0.7580	0.7552	−0.0510 ± j6.3772	0.0080
	Q = 0.3	$K_d$	0.5686	0.5686	0.7172	0.7162	−0.3160
A-4	P = 0.5	$K_s$	0.7266	0.7266	0.7226	0.7226	−0.0444 ± j6.2381	0.0071
	Q = 0.4	$K_d$	0.5672	0.5686	0.6231	0.6199	−0.3393
A-5	P = 0.5	$K_s$	0.6976	0.6988	0.6871	0.6866	−0.0392 ± j6.0809	0.0064
	Q = 0.5	$K_d$	0.5636	0.3757	0.5428	0.5506	−0.3606
A-6	P = 0.5	$K_s$	0.6209	0.6209	0.6458	0.6458	−0.0353 ± j5.8974	0.0060
	Q = 0.6	$K_d$	0.2278	0.2278	0.4962	0.4922	−0.3804
A-7	P = 0.5	$K_s$	0.6061	0.5926	0.5992	0.5992	−0.0324 ± j5.6805	0.0057
	Q = 0.7	$K_d$	0.6167	0.1780	0.4547	0.4546	−0.3993
A-8	P = 0.5	$K_s$	0.5654	0.5667	0.5445	0.5462	−0.0304 ± j5.4234	0.0056
	Q = 0.8	$K_d$	0.4351	0.4727	0.4330	0.4274	−0.4174
A-9	P = 0.5	$K_s$	0.4952	0.4989	0.4861	0.4865	−0.0294 ± j5.1185	0.0057
	Q = 0.9	$K_d$	0.4345	0.5546	0.4104	0.4095	−0.4348
A-10	P = 0.5	$K_s$	0.4248	0.4248	0.4199	0.4200	−0.0294 ± j4.7561	0.0062
	Q = 1.0	$K_d$	0.2889	0.2893	0.4170	0.4135	−0.4515

presents the values of $K_s$ and $K_d$ estimated by using EP, AIS, PSO, LS, ${\xi }_{min}$, and $\lambda$ for Event A. All the cases show the negative and positive values of $\lambda$ and ${\xi }_{min}$, indicating that all the cases are stable. In each case, a total of 10 replications for the simulation process were performed. The results obtained are consistent and within the range of not more than 1%. The LS method is selected as the standard value for EP, AIS, and PSO estimation methods.

The results in Table 3 show that the estimated values of $K_s$ and $K_d$ using the PSO method obtains more accurate values than EP and AIS. Particularly for the data of $K_s$ for Cases A-1, A-4, and A-7, the PSO estimation technique obtains exactly the same values as those calculated via the LS method. The AIS technique achieves the worst results, such that most of the estimated values are different from those calculated with the LS method.

Fitness of the EP, AIS, and PSO methods for Event A is illustrated in graph form in Figure 7. Figure 7a shows that the PSO technique obtains the highest fitness values for all cases. Fitness for EP and AIS are in the range of 0.92 to 0.98, whereas that of PSO is in the range of 0.99 to 1.00. This finding implies that PSO calculates better fitness results than EP and AIS.

The results of the computation time for Event A for all three estimation techniques are demonstrated in Figure 7b. The graph indicates that EP takes the longest computation time with an average time of approximately 60 s. The AIS method is the fastest among the three techniques because most of the simulation cases achieved an optimal solution below 30 s.

For Event A, PSO achieves better results in the objective function compared with EP and AIS. However, AIS simulates the estimation process in a short computation time and small iteration number. Overall, PSO is the best method in most of the discussed situations.

Table 4. $K_s$ and $K_d$ via EP, AIS, PSO, LS, eigenvalues and ${\xi }_{min}$ for Event B.

Cases	P, Q	${\mathit{K}}_{\mathit{s}}$, ${\mathit{K}}_{\mathit{d}}$	EP	AIS	PSO	LS	$\mathit{\lambda }$	${\mathit{\xi }}_{\mathit{min}}$
B-1	P = 1.0	$K_s$	0.7239	0.7247	0.7104	0.7104	−0.0249 ± j6.1853	0.0040
	Q = 0.15	$K_d$	0.2389	0.2448	0.3471	0.3493	−0.3593
B-2	P = 0.9	$K_s$	0.8063	0.8109	0.8065	0.8060	−0.0703 ± j6.5879	0.0107
	Q = 0.15	$K_d$	1.1691	1.3715	0.9681	0.9814	−0.2486
B-3	P = 0.8	$K_s$	0.8781	0.8640	0.8910	0.8912	−0.1068 ± j6.9272	0.0154
	Q = 0.15	$K_d$	1.5674	1.4611	1.4751	1.4975	−0.1604
B-4	P = 0.7	$K_s$	0.9864	0.9884	0.9876	0.9879	−0.1246 ± j7.2930	0.0171
	Q = 0.15	$K_d$	1.7431	1.8278	1.7439	1.7474	−0.1220
B-5	P = 0.6	$K_s$	1.0167	1.0223	1.0498	1.0514	−0.1351 ± j7.5236	0.0179
	Q = 0.15	$K_d$	1.8866	2.0740	1.8947	1.8925	−0.1016
B-6	P = 0.5	$K_s$	1.0165	1.0223	1.0760	1.0751	−0.1441 ± j7.6078	0.0189
	Q = 0.15	$K_d$	1.9198	2.0740	1.9984	2.0153	−0.0842
B-7	P = 0.4	$K_s$	1.1035	1.1068	1.0697	1.0699	−0.1531 ± j7.5894	0.0202
	Q = 0.15	$K_d$	2.0607	2.0699	2.1582	2.1455	−0.0667
B-8	P = 0.3	$K_s$	−1.0601	−1.0601	−1.0510	−1.0511	7.3361, −7.6625	−1.0000
	Q = 0.15	$K_d$	2.9546	2.9546	2.6562	2.6560	−0.0473
B-9	P = 0.2	$K_s$	−1.0136	−0.9449	−0.9875	−0.9991	7.1650, −7.5150	−1.0000
	Q = 0.15	$K_d$	3.1763	3.5296	1.8054	2.0112	−0.0245
B-10	P = 0.1	$K_s$	−0.8589	−0.8642	−0.9437	−0.9437	−7.3137, 6.9347	−1.0000
	Q = 0.15	$K_d$	−1.7892	−1.9245	2.7579	2.7580	0.0034

tabulates the values of $K_s$ and $K_d$ estimated via EP, AIS, PSO, LS, ${\xi }_{min}$, and $\lambda$ for Event B. From the 10 cases, the first 7 cases show positive ${\xi }_{min}$ and negative $\lambda$, indicating that the first 7 cases are in a stable condition. The three final cases are considered unstable because of the negative value of the minimum damping ratio, ${\xi }_{min}$, and positive eigenvalues, $\lambda$. The estimated values of $K_s$ and $K_d$ via the LS method supports the result of ${\xi }_{min}$ and $\lambda$. From the estimation process for the first seven cases, which used the LS method, the values of $K_s$ and $K_d$ are positive. These findings indicate that Cases B-1 to B-7 are stable. For Cases B-8, B-9, and B-10, the LS results give negative values for $K_s$ and positive values for $K_d$, thereby proving that all three final cases are considered non-oscillatory instability cases. In each case, a total of 10 replications for the simulation process were performed. The results obtained are within the range of not more than 1%.

From Table 4, PSO method exhibits more comparable results with respect to AIS and EP, particularly for the two final cases, namely Cases B-9 and B-10. In these two cases, the estimated value of $K_d$ using EP and AIS are far deviated from PSO and LS method.

The profile for the fitness of EP, AIS, and PSO method for Event B is presented in Figure 8a. From the graph, the PSO technique acquires the highest fitness values as compared with EP and AIS. The fitness values for EP and AIS are in the range of 0.92 to 0.98, whereas PSO is in the range of 0.99 to 1.00. This finding implies that PSO outperformed EP and AIS in determining the values of $K_s$ and $K_d$ with high accuracy.

Figure 8b lists the results of computation time for Event B for the EP, PSO, and AIS estimation techniques. EP is the slowest among the three techniques with the computation time in a range of 40–80 s. For the AIS method, 5 out of 10 cases finished the computation process in 20 s.

The results of Event B indicate that PSO can optimize the best value of torque coefficients $K_s$ and $K_d$. PSO reaches higher fitness values than EP and AIS. AIS obtains better computation time, but EP is computationally burdensome.

5. Conclusions

In this study, an efficient real-time estimation technique of torque coefficients $K_s$ and $K_d$ for solving angle stability problems is proposed. $K_s$ and $K_d$ are estimated on the basis of $\Delta \delta \left(t\right)$, $\Delta \omega \left(t\right)$, and $\Delta T_e\left(t\right)$. The values of $K_s$ and $K_d$ based on the LS technique are selected as the benchmark for evaluating the performance of the proposed optimization method. Overall, PSO is excellent for calculating the values of $K_s$ and $K_d$ because it obtains almost the same values as those calculated via the LS method. Although EP and AIS can determine the stability condition of all the cases, the value differences of $K_s$ and $K_d$ when compared with LS for EO and AIS are larger than those obtained by using PSO. From the perspective of computation time and iterations, AIS obtains the fastest simulation time in this event, followed by PSO and EP. Despite the large iteration number, the time consumed for the PSO simulation process is still minimal and acceptable. From the point of view of population size, the performance of calculation accuracy for all the three techniques significantly increased when simulations are conducted in 20 and 50 population sizes compared with that in 5 population sizes. The computation time consumed and number of iterations increase significantly regardless of the calculation accuracy performance in each optimization technique as the population sizes increase. The results of the last three events suggest that 20 is the recommended population size for calculating the accurate value of $K_s$ and $K_d$ via the PSO technique.