Parameter Estimation of Electromechanical Oscillation Based on a Constrained EKF with C&I-PSO

Abstract

By combining together the extended Kalman filter with a newly developed C&I particle swarm optimization algorithm (C&I-PSO), a novel estimation method is proposed for parameter estimation of electromechanical oscillation, in which critical physical constraints on the parameters are taken into account. Based on the extended Kalman filtering algorithm, the constrained parameter estimation problem is formulated via the projection method. Then, by utilizing the penalty function method, the obtained constrained optimization problem could be converted into an equivalent unconstrained optimization problem; finally, the C&I-PSO algorithm is developed to address the unconstrained optimization problem. Therefore, the parameters of electromechanical oscillation with physical constraints can be successfully estimated and better performed. Finally, the effectiveness of the obtained results has been illustrated by several test systems.

1. Introduction

Over the past few decades, there have been a lot of concerns about estimation of electromechanical modes since it can offer a substantial amount of important information about power system stability [1,2,3,4,5]. Therefore, it is necessary to track these parameters in real time to monitor the system stability and prevent blackouts [6,7,8]. In general, electromechanical oscillations are divided into two groups: local and inter-area, which refer to the oscillations between nearby generators and distant generators, respectively [3].

Generally speaking, there are mainly two different methods for the parameters estimation of power system oscillation: one of them is the model-based method, and the other is the measurement-based method. In the model-based method, the governing equations of the studied model are linearized near the present operating point. In contrast, a linear model can be estimated directly from measurement data in the measurement-based method. Therefore, in recent years, for complex modern power systems, the measurement-based method is generally considered easier than the model-based method and is widely used [3,9,10,11,12,13,14,15,16,17]. Normally, two types of measurement data can be obtained by using phasor measurement units (PMUs). One of them is the ambient data, which can be sampled from a power system at a stable operating point [10]; the other one is designed as the ringdown data, which is generated from a power system with a major disturbance. The ringdown data contains key oscillatory information and is the main focus of this paper. In recent years, for ringdown data detection, several methods have been developed. In [11,12], some methods are investigated that are mainly based on frequency estimation and harmonic detection. Furthermore, the matrix pencil method [13], the Kung′s algorithm [14], the rotational invariance technique [15], and the linear prediction methods such as the Prony method [16,17] were proposed and widely utilized for ringdown detection. However, most of the aforementioned approaches cannot be used for real-time application since they are sliding window-based methods, in which the calculations are performed only when a window of data is received [3]. Recently, a novel method for ringdown detection was proposed based on the traditional EKF in [3], where the damping ratio and frequency of ringdown data were modeled in the state vector; then these parameters could be easily estimated by using EKF. However, it should be pointed out that, in the above mentioned results, some critical physical constraints on parameters were ignored; thus, the practical value of the proposed method could be limited.

The Kalman filter is known as an optimal estimator and has been widely used in aerospace, the state estimation of power systems and other fields [18,19,20,21,22,23]. EKF is a famous nonlinear estimator, which is developed by combining the conventional Kalman filter and the linearized technology-based Taylor series expansion. Due to its advantages and convenience, it has been used in many nonlinear state estimation problems. However, in the application of EKF, some known information of the signal is sometimes either ignored or addressed heuristically. For example, the damping ratio and the frequency of oscillation signal are usually meet with the positive constraint condition. In general, for a constrained estimation problem, it is radically difficult within the recursive framework for the traditional EKF to incorporate the constraints into system states. To deal with this problem, many approaches have been developed, such as the projection method, the mean square method [24,25], etc., by which the constrained estimation problem can be successfully converted to the constrained optimization problem after each iteration of EKF [26,27].

The constrained optimization problems are common in many practical applications, such as optimal power flow calculations, allocation problems, and structural optimizations [28]. Until now, various solutions have been developed to address these optimization problems. The most popular method of them is particle swarm optimization (PSO) [29,30,31,32], which has been widely used in many optimization problems. It has the ability to solve difficult optimization problems and converge quickly to a solution [33,34,35,36]. Although the PSO was developed primarily to address unconstrained optimization problems; nevertheless, it performs well when applied to constrained optimization problems. In this case, by adding a penalty term into the objective function, an equivalent unconstrained optimization problem can be obtained; then PSO can be used to address the unconstrained optimization problem. However, in practice, although the traditional PSO algorithm can converge quickly, it is always easy to fall into local optimum [37,38,39,40,41], so the global optimization solutions cannot be obtained easily. Therefore, in this paper, a modified PSO algorithm (namely, C&I-PSO) is developed to address the equivalent unconstrained optimization problem, in which a constrictive factor and a linear decreasing inertia weight (LDIW) are introduced to improve the performance of global searches and local searches [42]. By this method, in the early search stage, a large inertia weight is utilized to extend the search region and avoid the occurrence of premature problems. In the later stage, a small inertia weight is used to obtain a more accurate local search such that a more accurate global optimization solution can be achieved. To the best of the authors’ knowledge, in existing work, there have been little literature that considers parameter estimations of electromechanical oscillation with physical constraints. Based on the above discussion, with the consideration of the physical constraints on parameters, a constrained EKF in combination with C&I-PSO algorithm is proposed in this paper for parameter estimation of the ringdown signal.

The rest of this paper is organized as follows: In Section 2, the state space model for the parameters estimation of ringdown signal is provided. In Section 3, by incorporating the projection method, the traditional EKF with inequality constraints on system states is presented. In Section 4, the C&I-PSO algorithm for the equivalent unconstrained optimization problem is proposed, and the new parameter estimation algorithm is presented. In Section 5, some simulation results are provided to verify the usefulness and effectiveness of the proposed results. Finally, the conclusions are presented in Section 6.

2. State-Space Model of the Ringdown Signal

As stated in [3,7], the ringdown signal can be expressed by the sequence of exponentially damped sinusoids (EDS):

$$y(t)={\displaystyle \sum _{i=1}^N{A}_i}{e}^{-{\delta }_{i}t}\cos ({\omega }_{i}t+{\varphi }_i)+n(t),$$

(1)

where $N$ represents the number of EDSs, $A_i$ is a known constant coefficient, ${\delta }_i$ is the damping ratio, ${\omega }_i$ represents the frequency, ${\varphi }_i$ is the initial phase value, the subscript $i$ indicates the $i$th component of signal, and $n(t)$ is considered a white Gaussian noise sequence with zero mean.

Assume that the signal sampling frequency is $f_s$. Then the ringdown signal containing multiple EDS signals can be deduced:

$$S_k={\displaystyle \sum _{i=1}^N{A}_i}{e}^{-{\delta }_i\frac{k}{f_s}}\cos ({\omega }_i\frac{k}{f_s}+{\varphi }_i)+n_k.$$

(2)

Define 4 $N$ state variables as [3]:

$$\{\begin{array}{l}{x}_{4i-3,k}=A_i{e}^{-{\delta }_i\frac{k}{f_s}}\cos ({\omega }_i\frac{k}{f_s})\\ x_{4i-2,k}=A_i{e}^{-{\delta }_i\frac{k}{f_s}}\sin ({\omega }_i\frac{k}{f_s})\\ x_{4i-1,k}={\omega }_i\\ x_{4i,k}={\delta }_i\end{array},$$

(3)

where the subscript $i$ indicates the $i$th component of the signal. Based on Equation (3), the $k+1$ instant state equation of the ringdown signal can be formulated as:

$$\{\begin{array}{l}{x}_{4i-3,k+1}=e^{-\frac{x_{4i,k}}{f_s}}\left(x_{4i-3,k}\cos (\frac{x_{4i-1,k}}{f_s})-x_{4i-2,k}\sin (\frac{x_{4i-1,k}}{f_s})\right)+{\omega }_{4i-3,k}\\ x_{4i-2,k}=e^{-\frac{x_{4i,k}}{f_s}}\left(x_{4i-3,k}\sin (\frac{x_{4i-1,k}}{f_s})+x_{4i-2,k}\cos (\frac{x_{4i-1,k}}{f_s})\right)+{\omega }_{4i-2,k}\\ x_{4i-1,k}=x_{4i-1,k}+{\omega }_{4i-1,k}\\ x_{4i,k}=x_{4i,k}+{\omega }_{4i,k}\end{array},$$

(4)

where ${\omega }_{4i-l,k}(l=0,1,2,3)$ denotes the system noise and modeling variations, which are usually taken as the zero-mean white noises. The in-phase signal is given by the state variables $x_{4i-3,k}$, $x_{4i-2,k}$ is the quadrature signal, $x_{4i-1,k}$ represent the frequency ${\omega }_i$ to be identified, $x_{4i,k}$ is the damping factor to be estimated, and the subscript $i$ indicates the $i$th component of ringdown signal.

The observation can be obtained by:

$$y_k={\displaystyle \sum _{i=1}^N{k}_{2i-1}{x}_{4i-3,k}+k_{2i}{x}_{4i-2,k}+n_k},$$

(5)

where $N$ indicates the total number of EDSs, the subscript $i$ indicates the $i$th component of the ringdown signal, $k_{2i-1}=\cos ({\varphi }_i),k_{2i}=-\sin ({\varphi }_i)$, and $n_k$ is a white Gaussian noise sequence with zero mean.

Remark 1.

It should be pointed out that most of the monitoring data of electromechanical oscillation in real power systems are stable or positively damped, which means the damping factor is positive, even in some poorly damped cases. Therefore, it is reasonable to assume the damping ratio and frequency are positive under this case. Based on the above discussion, following the idea proposed in [3], the next step of this research is to explore new methods to estimate parameters of electromechanical oscillation by taking into account the physical constraints.

3. Extended Kalman Filter with Inequality Constraints

In this section, the traditional EKF for parameter identification will be introduced in advance. Then, the projection method is utilized to turn the constrained estimation problem into an equivalent constrained optimization problem.

3.1. Traditional Extended Kalman Filter

EKF is usually used for states or the parameters estimation of nonlinear systems, which is developed by combining the conventional Kalman filter and the linearized technology based on Taylor series expansion.

Consider the following nonlinear stochastic system:

$$\{\begin{array}{l}{x}_{k+1}=f(x_k,u_k)+w_k\\ y_k=h(x_k)+v_k\end{array},$$

(6)

where $f(\cdot )$ is the system function, $h(\cdot )$ is the output function, both of them are nonlinear functions that could be linearized by Taylor series expansion, $x_k$ is the n-dimensional system state, $u_k$ is the input vector, $w_k$ is the state noise process vector denoting disturbances and modeling errors, $y_k$ is the m-dimensional observation vector, and $v_k$ is the measurement noise. Here $w_k$ and $v_k$ are usually considered as the zero-mean Gaussian white noises.

For convenience, the following formulas are defined in advance:

$${\widehat{x}}_k=E[(x_k|y_k^{\ast })],{\widehat{P}}_k=E[(x_k-{\widehat{x}}_k){(x_k-{\widehat{x}}_k)}^T]$$

(7)

$${\tilde{x}}_k=E[(x_k|y_{k-1}^{\ast })],{\tilde{P}}_k=E[(x_k-{\tilde{x}}_k){(x_k-{\tilde{x}}_k)}^T]$$

(8)

where ${\widehat{x}}_k$ is the state estimation, ${\widehat{P}}_k$ is the estimation error covariance matrix, ${\tilde{x}}_k$ is the state prediction, and ${\tilde{P}}_k$ is the prediction error covariance matrix, the subscript $k$ indicating time instant.

The prediction step:

$$\{\begin{array}{l}{\tilde{x}}_{k+1}=f({\widehat{x}}_k,u_k)\\ {\tilde{P}}_{k+1}=F_k{\widehat{P}}_k{F}_k^T+Q_k\end{array},$$

(9)

The filtering step:

$$\{\begin{array}{l}{\widehat{x}}_{k+1}={\tilde{x}}_{k+1}+G_{k+1}[y_{k+1}-h({\tilde{x}}_{k+1})]\\ G_{k+1}={\tilde{P}}_{k+1}{H}_{k+1}^T{(H_{k+1}{\tilde{P}}_{k+1}{H}_{k+1}^T+R_{k+1})}^{-1}\\ {\widehat{P}}_{k+1}=(I-G_{k+1}{H}_{k+1}^{}){\tilde{P}}_{k+1}\end{array},$$

(10)

where $G_{k+1}$ is the Kalman filter gain vector at $k+1$ time instant, and the Jacobian matrices $F_k$ and $H_{k+1}$ can be derived by:

$$F_k=\frac{\partial f(x_k,u_k)}{\partial x_k}{|}_{x_k={\widehat{x}}_k},\hspace{1em}{H}_{k+1}=\frac{\partial h(x_{k+1},u_{k+1})}{\partial x_{k+1}}{|}_{x_{k+1}={\tilde{x}}_{k+1}}$$

(11)

$Q_k$ and $R_k$ are the covariance matrices of state noise and measurement noise, respectively, and they are usually defined as:

$$Q_k=E[w_k{w}_k^T],\hspace{1em}{R}_k=E[v_k{v}_k^T]$$

(12)

Based on the model of ringdown signal (4) and (5), by using the EKF algorithm, some related equations for estimating parameters of the ringdown signal are obtained:

$$f(x_k)=\left[\begin{array}{c}{M}_1\\ M_2\\ ⋮\\ M_i\\ ⋮\\ M_N\end{array}\right],\hspace{1em}{M}_i=\left[\begin{array}{c}{x}_{4i-3,k}\\ x_{4i-2,k}\\ x_{4i-1,k}\\ x_{4i,k}\end{array}\right]$$

(13)

The function of $h(x_k)$, and the covariance matrix $Q_k$, $R_k$ can be also obtained:

$$H=(\begin{array}{ccccc}{k}_1{k}_{2}00& \cdots & k_{2i-1}{k}_{2i}00& \cdots & k_{2N-1}{k}_{2N}00\end{array}),h(x_k)=H{x}_k$$

(14)

$$Q_k=\left(\begin{array}{cccc}{q}_1& 0& \cdots & 0\\ 0& q_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & q_{4N}\end{array}\right),R_k=E[v_k{v}_k^T]=[r]$$

(15)

Remark 2.

In [3], a new method was proposed for ringdown detection, where the traditional EKF method was used to estimate the parameters. However, it should be pointed out that both the frequency and damping ratio of the ringdown signal are usually positive. Thus, the practical value of the method will be inevitably reduced if such state constraints are ignored. In order to obtain a better and more practical estimation algorithm, a new constrained method is proposed to estimate the parameters of the ringdown signal by accounting for the physical constraints.

3.2. The Projection Method

In most cases, it is necessary to specify constraints on parameter values. Therefore, in order to solve the constrained parameter identification problem, the constrained estimation problem can be successfully converted by using the projection method to the equivalent minimum constrained optimization problem at each iteration.

Consider the nonlinear dynamical system (6) with the constraint as follows:

$$U_{c}x\le u_c,$$

(16)

where $U_c$ indicates a known matrix, which has s lines and n columns. s represents the number of constraints, n indicates the number of state variables, usually the rank of $U_c$ is assumed to be s.

By using the projection method, the unconstrained state estimate $\tilde{x}$ can be projected onto the constraint surface directly, then the constrained EKF problem is expressed by:

$$\underset{\tilde{x}}{\min }{\left(\tilde{x}-\widehat{x}\right)}^T{W}_p\left(\tilde{x}-\widehat{x}\right),\hspace{1em}s.t.\hspace{1em}{U}_c\tilde{x}\le u_c$$

(17)

where $W_p>0$ is a weighting matrix. For simplicity, $W_p$ is chosen as the identify matrix in this paper.

The constrained EKF could be inferred by seeking out an estimation $\tilde{x}$, which satisfies the constraint condition (16) and can be obtained by solving the minimum constrained optimization problem expressed in (17). Therefore, with the projection method, the constrained state estimation problem can be addressed successfully by solving the constrained optimization after each iteration.

4. The C&I Particle Swarm Optimization

In this section, the penalty function method is utilized to change the minimum constrained optimization problem into an equivalent unconstrained optimization problem. Then, the C&I-PSO algorithm is developed to solve the equivalent unconstrained optimization problem.

4.1. The Penalty Function Approach

Due to its simple principle, the penalty function method has been used in many constrained optimization problems; after adding a penalty term, it can be applied to the constrained optimization problem, and an equivalent unconstrained optimization problem can be obtained.

Motivated by the idea proposed in [34], a non-stationary multistage assignment penalty function will be utilized to address the minimum constrained optimization, which can be expressed by:

$$f_a\left(x\right)=f\left(x\right)+h_d\left(k\right)H_p\left(x\right),x\in R^n$$

(18)

where $f\left(x\right)$ represents the minimum constrained objective function expressed in (17). $h_d\left(k\right)$ indicates a dynamical modified penalty value, k represents the iteration number of the algorithm, and $H_p\left(x\right)$ denotes a penalty factor given by:

$$H_p(x)={\displaystyle \sum _{i=1}^m{\theta }_f({\beta }_i(x))}{\beta }_i{(x)}^{\alpha ({\beta }_i(x))},$$

(19)

where ${\beta }_i(x)=\max \left\{0,g_i(x)\right\},(i=1,\dots s)$, which is a relative violated function of the constraints, and $g(x)=U_{c}x-u_c$ are the constraints expressed in (17); ${\theta }_f({\beta }_i(x))$ represents a multi-stage assignment function; $\alpha ({\beta }_i(x))$ is the power of the penalty function.

Based on the above method, the minimum constrained optimization problem expressed in (17) can be turned into an equivalent minimum unconstrained optimization problem, and the unconstrained objective function can be expressed by (18). In order to address this minimum unconstrained optimization problem, the C&I-PSO approach will be introduced in the next part of this paper.

4.2. The C&I-PSO Algorithm

The PSO is one of the most popular swarm intelligence optimization algorithms [29]; its main idea is based simulating simplified social models, such as bird flocking and fish schooling. It begins with an initial population, which is generated in random positions in the search space. Without loss of generality, the $D$-dimensional search space is assumed to be $S\subseteq R^D$; the swarm has $N$ particles that each having neither weight nor volume and each holding its own position and velocity. These particles update their positions over finite iterations. For the $i$th particle at iteration $L$, the position and velocity vectors are defined as $X_i^L=(x_{i1}^L,x_{i2}^L,,\cdots x_{iD}^L)$ and $V_i^L=({\nu }_{i1}^L,{\nu }_{i2}^L,\cdots {\nu }_{iD}^L)$, respectively. The best historical position visited by particle $i$ is taken as a point in $S$, which is denoted as $Pbes{t}_i^L=(pbes{t}_{i1}^L,pbes{t}_{i2}^L,\cdots ,pbes{t}_{iD}^L)$. The best $Gbes{t}^L=\left(gbes{t}_1^L,gbes{t}_2^L,\cdots gbes{t}_D^L\right)$. In the traditional PSO, each particle moves toward the Particle Best and Global Best positions by the specified velocity:

$$V_i^{L+1}=\psi \left(L\right)\times V_i^L+c_1{r}_{1,i}\left(Pbes{t}_i^L-X_i^L\right)+c_2{r}_{2,i}\left(Gbes{t}^L-X_i^L\right),$$

(20)

$$X_i^{L+1}=X_i^L+\gamma \times V_i^{L+1},$$

(21)

where $\gamma$ is a constrictive factor introduced to constrict and control velocities, and $\psi \left(L\right)$ is a linear decreasing inertia weight, which can be obtained by:

$$\psi (L)={\psi }_{start}\times ({\psi }_{start}-{\psi }_{end})\times (T_{\max }-L)/T_{\max },$$

(22)

where ${\psi }_{start}$ represents the inertia weight of initial iteration point, wend indicates the inertia weight of final iteration point, usually chosen as ${\psi }_{start}=0.9$ and ${\psi }_{end}=0.4$ respectively, $L$ is the number of iteration and $T_{\max }$ represents the maximum iteration number.

Remark 3.

It is worth pointing out that, for solving the optimization problem in parameter estimation of electromechanical oscillation, a modified PSO named the C&I-PSO is proposed, in which a constrictive factor and a linearly decreasing inertia weight are introduced to improve the performance of local and global search. By using the proposed C&I-PSO, the premature problem of the traditional PSO can be successfully addressed, and much better solutions to the global optimization problem can be achieved if compared with the traditional PSO. Now we can present the proposed approach for parameter estimation, which is formulated as the following Algorithm 1:

Algorithm 1: Constrained EKF with C&I-PSO

1: Initialization: Set appropriate values for ${\widehat{P}}_0$, ${\widehat{x}}_0$, $Q_0$, $R_0$, $S_t$, $M_t$;

2: for $k=0$ to $S_t$ do

3: Calculate the value of state prediction ${\tilde{x}}_{k+1}$ and prediction error covariance ${\tilde{P}}_{k+1}$:

4:            ${\tilde{x}}_{k+1}\leftarrow f\left({\widehat{x}}_k,u_k\right)$, ${\tilde{P}}_{k+1}\leftarrow F_k{\widehat{P}}_k{F}_k^T+Q_k$;

5: Compute the Kalman gain matrix at time instant $k+1$:

6:             $G_{k+1}\leftarrow {\tilde{P}}_{k+1}{H}_{k+1}^T{\left(H_{k+1}{\tilde{P}}_{k+1}{H}_{k+1}^T+R_{k+1}\right)}^{-1}$;

7: Update the state estimation ${\widehat{x}}_{k+1}$ and estimation error covariance:

8:       ${\widehat{x}}_{k+1}\leftarrow {\tilde{x}}_{k+1}+G_{k+1}\left[y_{k+1}-h\left({\tilde{x}}_{k+1}\right)\right]$, ${\widehat{P}}_{k+1}\leftarrow \left(I-G_{k+1}{H}_{k+1}\right){\tilde{P}}_{k+1}$;

9:  if $\min \left({\omega }_1,{\delta }_1,{\omega }_2,{\delta }_2\cdots {\omega }_N,{\delta }_N\right)<0$ then

10:               $\underset{\tilde{x}}{\min }{\left(\tilde{x}-\widehat{x}\right)}^T{W}_p\left(\tilde{x}-\widehat{x}\right)$, s.t. $U_c\tilde{x}\le u_c$;

11:  while $L\le M_t$ do

12:             $f_a\left(\tilde{x}\right)=f\left(\tilde{x}\right)+h_d\left(L\right)H_p\left(\tilde{x}\right)$, $\tilde{x}\in R^n$;

13:   for $i=1$ to $D$ do

14:     Update $Gbest$, $Pbes{t}_i$ and calculate particle new velocity and position:

15:              $\psi (L)\leftarrow {\psi }_{start}\times ({\psi }_{start}-{\psi }_{end})\times (M_t-L)/M_t$,

16:         $V_i^{L+1}\leftarrow \psi \left(L\right)\times V_i^L+c_1{r}_{1,i}\left(Pbes{t}_i^{iter}-X_i^L\right)+c_2{r}_{2,i}\left(Gbes{t}^L-X_i^L\right)$,

17:              $X_i^{L+1}=X_i^L+\gamma \times V_i^{L+1}$;

18:   end for

19:                 $L=L+1$;

20:  end while

21:                 ${\widehat{x}}_{k+1}=Gbest$;

22:  end if

23:                  $k=k+1$;

24: end for

Remark 4.

By using the proposed algorithm, the problem of parameter estimation for the ringdown signal with physical constraints can be solved successfully. It should be noted that not only can the parameters of the ringdown signal be truly identified with physical constraints, but a shorter convergence time can also be expected if compared with the existing results in [3]. Therefore, the proposed algorithm is expected to be much more realistic than the previous work.

5. Simulation Results and Discussions

In this section, three different test systems are provided to demonstrate the performance of the proposed method. The parameters in the penalty function are chosen using the same values as those in [34]. Specifically, the function $h_d\left(k\right)$ is set as $h_d\left(k\right)=k\sqrt{k}$, other related function values are given by:

$$\alpha \left({\beta }_i\left(x\right)\right)=\{\begin{array}{c}1 \mathrm{when} {\beta }_i\left(x\right)<1\\ 2 \mathrm{when} {\beta }_i\left(x\right)\ge 1\end{array},$$

(23)

$${\theta }_f\left({\beta }_i\left(x\right)\right)=\{\begin{array}{ll}10& \mathrm{when} {\beta }_i\left(x\right)<0.001\\ 20& \mathrm{when} 0.001\le {\beta }_i\left(x\right)<0.1\\ 100& \mathrm{when} 0.1\le {\beta }_i\left(x\right)\le 1\\ 300& \mathrm{when} {\beta }_i\left(x\right)>1\end{array}.$$

(24)

In addition, in order to evaluate the whole performance of the proposed approach and conventional EKF method effectively, the root-mean-square deviation (RMSD) used in [43] is adopted:

$$RMSD=\sqrt{\frac{{\displaystyle \sum _{k=1}^{k=n}{\displaystyle \sum _{i=1}^{i=m}{\left({\widehat{x}}_{k,i}-x_{k,i}\right)}^2}}}{n}},$$

(25)

where ${\widehat{x}}_{k,i}$ is the estimation result of the $i$th corresponding parameter, $x_{k,i}$ is the true value of the $i$th corresponding parameter, $m$ is the number of parameters, and $n$ is the number of time steps.

5.1. Test System 1: Ringdown Signal Composed of One EDS Signal

At first, the ringdown signal in [7] is utilized to evaluate performance of the proposed algorithm. The damping ratio and frequency of this ringdown signal are 0.01 and 1 rad/s, respectively, which is provided by:

$$y_1(t)=e^{-0.01t}\cos (t)+n(t).$$

(26)

The input ringdown signal is shown in Figure 1. By utilizing the Matrix Pencil [13], Prony [17], EKF [3] methods and the proposed algorithm, the damping factor and frequency could be estimated. Figure 2 shows the estimation results, and the RMSD comparison of different approaches is also provided in Figure 3. It can be found that the proposed approach is able to estimate parameters correctly in a short time, while EKF and other methods need a long time to converge to the real values. It follows from these simulation results that, a better convergence and steady state performance could be obtained if compared with those by using EKF, Matrix Pencil and Prony methods.

In order further to demonstrate the performance of the proposed algorithm, the final identification results of different methods are also shown in

Table 1. Identified results of different methods.

Method	$\mathbf{Damping}\mathbf{Factor}{\mathit{\delta }}_{\mathbf{1}}$	$\mathbf{Frequency}{\mathit{\omega }}_{\mathbf{1}}(\mathbf{rad}/\mathbf{s})$
Matrix Pencil [13]	0.0093	1.0104
Prony [17]	0.0087	1.0223
EKF [3]	0.0091	1.0109
Proposed Method	0.0099	1.0001

. It is observed from the results that the proposed method with smallest estimation errors than other approaches in [3,13,17]; which proves that the proposed method can estimate the parameters of ringdown signal accurately.

5.2. Test System 2: Ringdown Signal Composed of Two EDS Signals

In this case, a ringdown signal consisting of two EDS signals is considered:

$$y(t)=y_2\left(t\right)+y_3\left(t\right)+n(t),$$

(27)

where:

$$y_2\left(t\right)=e^{-0.005t}\cos \left(0.2t\right), y_3\left(t\right)=e^{-0.01t}\cos \left(0.6t\right).$$

(28)

and $n(t)$ denotes white noise. The input signal is presented in Figure 4. The estimation results of frequency and damping factor by using different methods are shown in Figure 5 and Figure 6. Figure 7 shows the RMSD comparison of the methods. It can be seen from these simulation results that a better convergence and steady state performance could be obtained if compared with those when using other methods in [3,13,17].

In addition, the final identification results of different methods are also provided in

Table 2. Identified results of different methods.

Method	$\mathbf{Damping}\mathbf{factor}{\mathit{\delta }}_{\mathbf{2}}/{\mathit{\delta }}_{\mathbf{3}}$	$\mathbf{Frequency}{\mathit{\omega }}_{\mathbf{2}}/{\mathit{\omega }}_{\mathbf{3}}(\mathbf{rad}/\mathbf{s})$
Matrix Pencil [13]	0.0049/0.0094	0.1983/0.5887
Prony [17]	0.0045/0.0089	0.1864/0.5843
EKF [3]	0.0048/0.0092	0.1967/0.5876
Proposed Method	0.0050/0.0101	0.2001/0.6001

, which further proves superior performance of the proposed method.

5.3. Test System 3: WSCC Model

In this case, the WSCC model is taken as the test system [42], which has three generators and nine buses. The network is depicted in Figure 8a. The normal frequency of G2 is 60 HZ, at t = 5.1 s, the load at bus 5 is removed. This disturbance causes an inter-area sustained oscillation in G2 with the frequency of ${\omega }_4=2.4rad/s$, and the damping factor is ${\delta }_4=0.3$. Figure 8b presents the event recorded in the bus frequency signal.

Figure 9 shows the estimated frequency and damping factor when using the proposed algorithm and other methods. Figure 10 illustrates the RMSD comparison of different methods. It could be found from these simulation results that, in this case, the proposed algorithm can estimate the parameters effectively with better convergence and steady state performance than EKF and other methods. Finally, the final identification results of other methods are also provided in

Table 3. Identified results of different methods.

Method	$\mathbf{Damping}\mathbf{Factor}{\mathit{\delta }}_{\mathbf{4}}$	$\mathbf{Frequency}{\mathit{\omega }}_{\mathbf{4}}(\mathbf{rad}/\mathbf{s})$
Matrix Pencil [13]	0.2808	2.3875
Prony [17]	0.2732	2.3873
EKF [3]	−0.0488	2.3666
Proposed Method	0.2998	2.4000

, which demonstrate the superior performance of the proposed mehthod. Therefore, it can be concluded that the proposed algorithm possesses much more practical value than those shown in previous work.

6. Conclusions

In this paper, with the consideration of the physical constraints on parameters, the constrained parameter estimation problem of ringdown signals has been discussed in detail. Based on the EKF algorithm and the projection method, the constrained estimation problem was converted into a minimum constrained optimization issue. Thus, the constrained optimization issue can be solved successfully by the proposed C&I-PSO with the penalty function method. Finally, the proposed algorithm was successfully applied to estimate the constrained parameters of the ringdown signal. Simulation results have demonstrated in most cases that the proposed algorithm can achieve much better performance and makes more practical sense than using the Matrix Pencil, Prony, and traditional EKF approaches.