Coherency and Online Signal Selection Based Wide Area Control of Wind Integrated Power Grid

TL;DR

This paper presents a real-time wide area control method for wind-integrated power grids that uses online coherency grouping and signal selection to improve damping of interarea oscillations.

Contribution

It introduces a novel WAC design combining Kalman filtering and residue-based signal selection for real-time adaptation in wind power systems.

Findings

Effective damping of interarea oscillations demonstrated

Online coherency grouping improves control accuracy

Method outperforms traditional local and offline WAC approaches

Abstract

This paper introduces a novel method of designing wide area control (WAC) based on a discrete linear quadratic regulator and Kalman filtering based state-estimation that can be applied for real-time damping of interarea oscillations of wind integrated power grid. The main advantages of the proposed method are that the architecture provides online coherency grouping that properly characterizes real-time changes in the power grid and online wide-area signal selection based on residue method for proper selection of the WAC signals. The proposed architecture can, thus, accurately monitors changes in the power grid and select the appropriate control signal for more effectively damping the interarea oscillation when compared to the conventional local signal based power system stabilizers or offline based WAC designs. The architecture is tested on a wind integrated two-area system and the IEEE 39 bus system in order to show the capability of the proposed method.

Tables

Table 1. TABLE I : COMPARISON OF COHERENCY GROUPING METHODS

Test System	Slow Coherency	Spectral Clustering
Two Area	Group-1:1,2	Group-1:1,2
	Group-2:3,4	Group-2:3,4
	$t$ = 0.4826s	$t$ = 0.021s
IEEE 39-BUS	Group-1:4,5,6,7	Group-1:4,5,6,7,9
	Group-2:1,8,9	Group-2:1,8
	Group-3:2,3	Group-3:2,3
	Group-4:10	Group-4:10
	$t$ = 0.8138s	$t$ = 0.0423s
*$t$ = Computational Time

Table 2. TABLE II : COHERENCY GROUPING OF IEEE-39 BUS SYSTEM (SPECTRAL CLUSTERING)

Fault Bus
No of Clusters (k)	Bus-14	Bus-19	Bus-6
4	Computational Time = 0.0423s	Computational Time = 0.0303s	Computational Time = 0.0696s
	Group-1:4,5,6,7,9	Group-1:6,7,9	Group-1:5,9
	Group-2:1,8	Group-2:1,2,3,8	Group-2:1,4,6,7,8
	Group-3:2,3	Group-3:4,5	Group-3:2,3
	Group-4:10	Group-4:10	Group-4:10

Table 3. TABLE III : CONTROL LOOP BASED ON RESIDUE APPROACH

Type	Residue analysis using MIMO identification				Residue analysis of state space matrices
Coherent groups	Group 1		Group 2		Group 1		Group 2
Speed/Field voltage	$u_1$	$u_2$	$u_3$	$u_4$	$u_1$	$u_2$	$u_3$	$u_4$
$\Delta {\omega }_1$	0.5640	0.5280	0.6398	0.6323	0.4973	0.6081	0.5263	0.6537
$\Delta {\omega }_2$	0.4323	0.6571	1.0000	0.6385	0.3551	0.4342	0.3758	0.4668
$\Delta {\omega }_3$	0.8758	0.7143	0.7366	0.9602	0.7608	0.9303	0.8051	1.0000
$\Delta {\omega }_4$	0.7093	0.6609	0.8584	0.7826	0.6733	0.8232	0.7125	0.8850
	Computational Time (Signal Selection) = 1.28s				Computational Time (Signal Selection) = 0.078s
	Computational Time (WAC Output Signal) = 0.0116s

Table 4. TABLE IV : IEEE 39-BUS CONTROL LOOP (FAULT ON BUS-14)

	Control Loop	Residue
Group-1	$\Delta {\omega }_{10}$ $\to$ $u_7$	0.678
Group-2	$\Delta {\omega }_7$ $\to$ $u_8$	0.49
Group-3	$\Delta {\omega }_6$ $\to$ $u_2$	1
Group-4	$\Delta {\omega }_1$ $\to$ $u_{10}$	0.027
Computational Time (Signal Selection) = 3.16s
Computational Time (WAC Output Signal) = 0.0192s

Equations

$$S_{ij}=e^{(-\frac{\Vert x_i-x_j\Vert }{2{\sigma }^2})}$$

$$\displaystyle\textbf{S}=\left[{\begin{array}[]{cc}\textbf{A}&\textbf{B}\\ \bf{B^{T}}&\textbf{C}\\ \end{array}}\right],\textbf{W}=\left[{\begin{array}[]{c}\textbf{A}\\ \bf{B^{T}}\\ \end{array}}\right]$$

$$\displaystyle\textbf{S}\approx\bar{\textbf{S}}=\left[{\begin{array}[]{cc}\textbf{A}&\textbf{B}\\ \bf{B^{T}}&\bf{B^{T}A^{-1}B}\\ \end{array}}\right]$$

$$\bar{\text{L}}=\text{I}-{\text{D}}^{-\frac{1}{2}}\bar{\text{S}}{\text{D}}^{-\frac{1}{2}}$$

$$\bar{\text{L}}=\bar{\text{V}}\bar{\mathbf{\Sigma }}{\bar{\text{V}}}^T$$

$$\bar{\mathbf{V}}=[\overrightarrow{v_1},\overrightarrow{v_2},....,\overrightarrow{v_n}]$$

$${\mathbf{U}}_{\mathbf{im}}=\frac{{\bar{\mathbf{V}}}_{\mathbf{im}}}{\sqrt{{\sum }_{\mathbf{r}=1}^{\mathbf{k}}{{\bar{\mathbf{V}}}_{\mathbf{ir}}}^2}},i=1,...,n\text{and}m=1,...,j$$

$$\begin{array}{rl}\mathbf{\Delta }\omega (\mathbf{z})& =\mathbf{G}(\mathbf{z})\mathbf{U}(\mathbf{z})\end{array}$$

$$\displaystyle\begin{split}\bf{\Delta\omega}(z)&=\left[{\begin{array}[]{c}\Delta\omega_{1}(z)\\ .\\ .\\ \Delta\omega_{m}(z)\end{array}}\right],\bf{U}(z)=\left[{\begin{array}[]{c}u_{1}(z)\\ .\\ .\\ u_{p}(z)\\ \end{array}}\right]\\ \bf{G}(z)&=\left[{\begin{array}[]{cccc}G_{11}(z)&.&.&G_{1p}(z)\\ .&.&.&.\\ .&.&.&.\\ G_{m1}(z)&.&.&G_{mp}(z)\\ \end{array}}\right]\end{split}$$

$$\begin{array}{rl}{\mathbf{G}}_{\mathbf{mp}}(\mathbf{z})& =\frac{\Delta {\omega }_m(z)}{u_p(z)}\\ & =\frac{b_0^h+b_1^h{z}^{-1}+...+b_k^h{z}^{-k}}{1+a_1{z}^{-1}+a_2{z}^{-2}+...+a_k{z}^{-k}}\end{array}$$

$${\mathbf{X}}_{\mathbf{His}}^{\mathbf{h}}={\mathbf{X}}_{\mathbf{Num}}^{\mathbf{h}}+{\mathbf{X}}_{\mathbf{Den}}^{\mathbf{h}}$$

$$\displaystyle\bf{X_{His}^{h}}=\left[{\begin{array}[]{c}\Delta\omega_{1}(z)\\ .\\ .\\ \Delta\omega_{m}(z)\end{array}}\right]$$

$$\displaystyle\bf{X_{Num}^{h}}=\left[{\begin{array}[]{ccc}u_{p}(j-1)&.&u_{p}(j-k)\\ .&.&.\\ .&.&.\\ u_{p}(j-N)&.&u_{p}(j-N+1-k)\\ \end{array}}\right]\left[{\begin{array}[]{c}b_{0}^{h}\\ .\\ .\\ b_{k}^{h}\\ \end{array}}\right]$$

$$\displaystyle\bf{X_{Den}^{h}}=\left[{\begin{array}[]{ccc}\Delta\omega_{m}(j-1)&.&\Delta\omega_{m}(j-k)\\ .&.&.\\ .&.&.\\ \Delta\omega_{m}(j-N)&.&\Delta\omega_{m}(j-N+1-k)\\ \end{array}}\right]\left[{\begin{array}[]{c}a_{1}\\ .\\ .\\ a_{k}\\ \end{array}}\right]$$

$$\displaystyle\left[{\begin{array}[]{c}\bf{X_{His}^{1}}\\ \bf{X_{His}^{2}}\\ .\\ \bf{X_{His}^{h}}\\ \end{array}}\right]=\left[{\begin{array}[]{c}\bf{X_{Num}^{1}}\\ \bf{X_{Num}^{2}}\\ .\\ \bf{X_{Num}^{h}}\\ \end{array}}\right]+\left[{\begin{array}[]{c}\bf{X_{\Delta\omega}^{1}}\\ \bf{X_{\Delta\omega}^{2}}\\ .\\ \bf{X_{\Delta\omega}^{h}}\\ \end{array}}\right]\left[{\begin{array}[]{c}a_{1}\\ .\\ .\\ a_{k}\\ \end{array}}\right]$$

$$\displaystyle\left[{\begin{array}[]{c}\bf{X_{\Delta\omega}^{1}}\\ \bf{X_{\Delta\omega}^{2}}\\ .\\ \bf{X_{\Delta\omega}^{h}}\\ \end{array}}\right]\left[{\begin{array}[]{c}a_{1}\\ .\\ .\\ a_{k}\\ \end{array}}\right]=\left[{\begin{array}[]{c}\bf{X_{His}^{1}}\\ \bf{X_{His}^{2}}\\ .\\ \bf{X_{His}^{h}}\\ \end{array}}\right]-\left[{\begin{array}[]{c}\bf{X_{Num}^{1}}\\ \bf{X_{Num}^{2}}\\ .\\ \bf{X_{Num}^{h}}\\ \end{array}}\right]$$

$$\displaystyle\left[{\begin{array}[]{c}\bf{X_{up}^{h}}\end{array}}\right]\left[{\begin{array}[]{c}b_{0}^{h}\\ .\\ .\\ b_{k}^{h}\\ \end{array}}\right]=\left[{\begin{array}[]{c}\bf{X_{His}^{h}}\end{array}}\right]-\left[{\begin{array}[]{c}\bf{X_{\Delta\omega}^{h}}\end{array}}\right]\left[{\begin{array}[]{c}a_{1}\\ .\\ .\\ a_{k}\\ \end{array}}\right]$$

$$\begin{array}{rl}{\mathbf{G}}_{\mathbf{mp}}(\mathbf{s})& =\frac{\Delta {\omega }_m(s)}{u_p(s)}\\ & =\frac{r_{mp}(1)}{s-p_1}+\frac{r_{mp}(2)}{s-p_2}+..+\frac{r_{mp}(j)}{s-p_j}+k_{mp}(s)\end{array}$$

$$\begin{array}{rl}\mathbf{x}(\mathbf{p}+1)& =\mathbf{Ax}(\mathbf{p})+\mathbf{Bu}(\mathbf{p})\\ \mathbf{y}(\mathbf{p})& =\mathbf{Cx}(\mathbf{p})+\mathbf{Du}(\mathbf{p})\end{array}$$

$$\mathbf{J}$$

$$u_p$$

$${\mathbf{K}}_{\mathbf{p}}$$

$$\begin{array}{rl}{\mathbf{P}}_{\mathbf{p}-1}& =\mathbf{Q}+{\mathbf{A}}^{\mathbf{T}}{\mathbf{P}}_{\mathbf{p}}\mathbf{A}\\ & -{\mathbf{A}}^{\mathbf{T}}{\mathbf{P}}_{\mathbf{p}}\mathbf{B}(\mathbf{R}+{\mathbf{B}}^{\mathbf{T}}{\mathbf{P}}_{\mathbf{p}}\mathbf{B}{)}^{-1}{\mathbf{B}}^{\mathbf{T}}{\mathbf{P}}_{\mathbf{p}}\mathbf{A}\end{array}$$

$$\bar{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}$$

$$\bar{\mathbf{L}}=\mathbf{A}\mathbf{L}{\mathbf{A}}^{\mathbf{T}}+\mathbf{Q}$$

$$\mathbf{G}=\bar{\mathbf{L}}{\mathbf{H}}^{\mathbf{T}}{(\mathbf{H}\bar{\mathbf{L}}{\mathbf{H}}^{\mathbf{T}}+\mathbf{R})}^{-1}$$

$$\mathbf{x}=\bar{\mathbf{x}}+\mathbf{G}(\mathbf{z}-\mathbf{H}\bar{\mathbf{x}})$$

$$\mathbf{L}=\bar{\mathbf{L}}-\mathbf{KH}\bar{\mathbf{L}}$$

Full text

I Introduction

Interconnected power systems exhibit dominant inter-area oscillations with the frequency of oscillation between 0.1-0.8 Hz [1]. During inter-area oscillations, coherent generators tend to swing together in groups against other groups based on changes in system conditions [2]. Inter-area oscillations poses significant problems in the operation of power system as it limits power transfer capability of tie-lines and also deteriorates power system security [3]. To suppress the inter-area oscillations, power system stabilizers (PSS) have been designed. PSS provides supplementary damping through synchronous generator exciters [4], high-voltage direct current (HVDC) links [5], and flexible alternating current transmission systems (FACTS) devices [6].

Conventional damping controllers are designed considering specific operating conditions that are typical for the power system [7]. Effectiveness of such controllers decreases when actual operating conditions of the power system deviate from the specified conditions used for the design of such damping controllers. With the increase in integration of renewable energy to the grid, the operating conditions of power grid changes more frequently especially due to the variability in the power output of these renewable resources. Under such variable operating conditions, local robust damping controllers fails to perform well. For example, in [8, 9, 10] a robust damping controller is designed considering a dominant operating condition with bounded uncertainty to make the controller effective to use under varying operating conditions. As the design is considering a dominant frequency at one operating condition, performance of these controllers may not be optimal as the operating condition deviates from the dominant one.

Due to the recent developments in wide-area measurement system (WAMS), adaptive damping controller design has drawn a lot of attention. Work performed in [11],[12] and development of algorithms for applying system identification techniques in power system [13] are critical in these aspect. However, most of the research till now is focused on tuning of local adaptive controllers while the issue of coordination among all the controller is not addressed well yet. It is worth noting that both adaptive design of local individual controller and the coordination between different controllers can be achieved simultaneously through a properly designed WAC.

There are two major challenges in designing an effective real-time WAC. First, there should be a robust method for obtaining a fast online model of the power system that can capture all low-frequency oscillation modes, and second, there should be an online method for selecting most observable and controllable signals for the control loop. For successfully managing these challenges, an online coherency grouping method is required, as such grouping reduces the number of signals to be monitored for the design of WAC [14]. The reduction in the number of signals decreases computation burden and at the same time facilitates online design and implementation of the WAC. Two types of coherency grouping methods are reported in the literature. The first method is a modal based approach in which system eigenvalues and eigenvectors are analyzed [2]. Modal based methods are further classified into time domain based slow coherency, and weak link methods. These methods are groups the generators at a specific operating point independent of fault locations. The second method is measurement based [15], which generally uses dynamic response captured from the measured data that characterizes underlying phenomena of oscillations. In this approach coherency grouping is performed by processing measured data using mathematical functions like Fast Fourier transform [16], Hilbert-Huang transform [17], energy function based [18], artificial neural networks [19], self-organizing maps [20], graph theory [21], principle component analysis (PCA) [20, 22], and k-means clustering [23]. Since these methods are based on measurements, it automatically considers and adapts to changes in system operating conditions. However excessive computational burden and bandwidth requirements limit the success of existing coherency grouping based on measured data. In general, WAC designs are based on the selection of control loops. In [24, 25] many techniques for selection of control loop for WAC design are illustrated. State-of-the-art methods include residue, and relative gain array (RGA) based approaches. However, they are offline algorithms evaluated based on linearizing the system at a specific operating point and analyzing the eigenvalues and eigenvectors at that instant, making the approach to fail when the operating point changes.

Significant advances have been made on applications of system identification on power system [14, 26, 27, 28, 29]. References [26, 27] shows the application of system identification for estimating electromechanical modes and mode shapes. Subspace state model of the system has also been designed for different operating conditions of the power system by using subspace state-space system identification (N4SID) method [29], or recursive adaptive stochastic subspace identification (RASSI) method [28]. However, obtaining subspace state-space model is computationally challenging. Reference [14] shows that identifying the power system using auto-regressive moving average exogenous (ARMAX) model structure based multiple input multiple outputs (MIMO) can accurately capture dynamics of the actual power system. Further in [14] an adaptive and coordinated damping controller is designed based on the ARMAX-MIMO model.

In this paper, a novel method of designing online WAC is introduced. This paper is an extended version of [30]. In this extended work, three major novel contributions are proposed. First, an online computationally robust coherency grouping of machines is performed based on spectral clustering which uses measurement signals for grouping. This method considers current system operating condition. In order to overcome the computational burden and to reduce the execution time, only the slow eigenvalues of the Laplacian matrix are used for clustering. Second, the control/loop for WAC is selected based on residue analysis of gain matrix formulated online using MIMO identification which considers the changes in operating conditions of the system. Third, a novel controller is designed based on discrete LQR and state estimation using Kalman filtering. The efficacy of the proposed method is verified by evaluating performance of the architecture to successfully damp the inter-area oscillations on a WTG integrated two-area and IEEE 39 bus system models in real-time using RSCAD/RTDS real-time digital simulator. To summarize, the major contributions of this work are:

$\bullet$

A new coherency grouping method based on online spectral clustering.

$\bullet$

A new method for online selection of wide-area control loop using MIMO identification.

$\bullet$

A novel method for WAC based on discrete LQR and state estimation using Kalman filtering.

The rest of the paper is organized as follows: In section II the proposed framework for the online WAC design is discussed. In section III, the results on a power grid using two-area test system and IEEE 39-bus system are discussed and Section IV concludes the paper.

II Proposed Framework for the Online Wide Area Controller

The proposed framework for the online WAC design involves following steps: a) coherency grouping algorithm, b) signal selection algorithm for WAC based on residue approach and, c) WAC design algorithm.

II-A Proposed coherency grouping algorithm

In this approach coherency grouping of generators is performed based on spectral clustering using measured data points which is the speed of the generators. Data points $x_{1},x_{2},...x_{n}$ for a window length of $n$ are considered for clustering. Further, using the data points a similarity matrix $S\subseteq{R}^{mn}$ is formulated, where $S_{ij}$ gives the relation between $x_{i}$ and $x_{j}$. The similarity matrix information is used to group $x_{1},x_{2},...x_{n}$ into $k$ clusters. The similarity matrix is a based on a Gaussian function represented as in (1)

[TABLE]

where $\sigma$ is a scaling factor. Here $S$ is dense and is of the order $n\times n$. The size of $S$ increases with increase in the number of data points under consideration, but this slows the simulation speed. To increase the online clustering speed Nystrom method is used which uses sub-matrix of the dense matrix S. Let a sub-matrix A of the dense matrix S can be represented as an $l\times l$ matrix (where $l<<n$), sub-matrix B as $l\times(n-l)$ matrix, and $C$ as $(n-l)\times(n-l)$ matrix. Upon rearranging the columns and rows S can be represented as

[TABLE]

where $\textbf{A}\subseteq R^{l\times l}$, $\textbf{B}\subseteq R^{l\times(n-l)}$ , and $\textbf{C}\subseteq R^{(n-l)\times(n-l)}$. $S$ is approximated by Nystrom method using $A$ and $B$. The approximation $S$ is given in (9)

[TABLE]

The normalized laplacian matrix using the approximated similarity matrix S can then be represented as (10)

[TABLE]

where $\textbf{D}=\sum_{j=1}^{n}\bar{S_{ij}}$ is a diagonal matrix. The decomposition of $\bar{\textbf{L}}$ gives the eigenvalues and corresponding eigenvectors according to (11).

[TABLE]

where $\bar{\bf{\Sigma}}$ contains eigenvalues, and $\bar{\bf{V}}$ are the corresponding eigenvectors. If $j$ slow eigenvalues are considered, then the eigenvectors corresponding to $j$ eigenvalues, written as an $R^{n\times j}$ matrix can be formulated as in (12)

[TABLE]

where $\vec{v_{i}}\subseteq R^{n},i=1.....j$ and $j$ is the number of eigenvectors. The normalized eigenvector matrix can then be written as

[TABLE]

The $n$ rows of $\bf{U}$ can easily be clustered into groups using k-means method (see Algorithm 1). In Algorithm-1 the inner-loop repeatedly assigns each set of entities $u_{i}$ to the closest cluster center $\mu_{j}$, and recalculate the new cluster center ($\mu_{j}$) based on the statistical mean of the data points assigned to it. More details of k-means method can be found in [31, 32].

Fig. 1 shows the flow chart of the proposed algorithm. The efficacy of the proposed algorithm is verified by implementing the architecture on Kundur two area (Fig. 2) and IEEE-39 bus power system models (Fig. 3), and then comparing the method with state-of-the-art offline slow coherency based grouping. For two-area system a 3-ph fault is created on bus-9 for a duration of 0.1 sec. The coherent groups obtained based on the proposed algorithm and comparisons are shown in Table I. It can be seen that the proposed online method provides same grouping as offline method but with lower computational time. Also, for IEEE 39 bus system, group 1 is different with generator 9 included in group 1. This is wrongly interpreted in the slow coherency method. Fig. 4 to Fig. 7 shows the comparison of coherency grouping for a fault on bus-14.

Further to validate the algorithm for changing operating conditions in real-life, a sequence of disturbances are studied at various locations of IEEE-39 bus system. For this, 3-ph faults are created for a duration of 0.1 sec at bus-14, bus-19, and bus-6 at 5, 31 and 61 sec respectively. Table II shows the coherency grouping comparison of 39-bus system for various operating conditions. Fig. 8 shows the coherency grouping for various operating conditions for IEEE-39 bus system. It can be observed that with the proposed coherency grouping method the generators can be grouped based on current operating condition whereas in offline methods like slow coherency the generator grouping does not change irrespective of changing operating conditions. The dynamics observed at different operating conditions in Fig. 8 supports this statement. It is also worth noting that the computational time required for the proposed online coherency grouping method is much less confirming that the algorithm is feasible to implement on a real-life system.

For MIMO identification of the power system, input and output signals need to be selected. The output signal is taken as the speed of synchronous machines and the input signal is taken as the voltage input to the automatic voltage regulator (AVR) and exciter of the synchronous machines. The measurement points for input and output is shown in Fig. 9. If there are $m$ output and $p$ input signals selected, then (14) represents the MIMO system representation of the power system in discrete time domain. For signal selection, the MIMO identification is performed first to identify the transfer functions in (14). A Recursive least square identification (RLS) algorithm [33, 34, 35, 36, 37] is proposed for the MIMO system identification.

Further, based on the MIMO identification, transfer function of the power system can be developed as shown in equation (15) [33, 34, 35, 36, 37]. In (15), $\Delta\omega_{m}$ is the difference between rated and actual speed of the generator, $u_{p}$ is the input signal (Fig. 9), $k$ the order of transfer function and $h$ is the element number in matrix. From (15), it can be seen that all transfer functions in (14) have same denominator coefficients but unique numerator coefficients. This ensures that coupling between different loops are considered while identifying the MIMO system. Also, this indicates that the power systems modes as seen from any input are same however their corresponding residues differ. The identification of MIMO system is divided into three steps which are discussed in the following subsections.

[TABLE]

where

[TABLE]

[TABLE]

II-B *Signal selection algorithm for WAC design *

II-B1 Step-1

For $N$ observation window length, and $j$ samples, rewrite (15) as shown in (16).

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

II-B2 Step-2

Concatenate $\bf{X_{His}^{h}}$ for all input-output combinations as shown in (53)

[TABLE]

II-B3 Step-3

Calculate denominator coefficients and numerator coefficients iteratively. For this, in the first iteration numerator coefficients are initialized. Then the denominator coefficients $(a_{1},a_{2}...a_{k})$ are calculated by applying least squares technique to (70). Further the numerator coefficients $(b_{0}^{h},b_{1}^{h}...b_{k}^{h})$ are calculated again as shown in (82)

[TABLE]

[TABLE]

The numerator and denominator are calculated iteratively until the desired tolerance is achieved such that $||\left([\bf{X_{His}^{h}}-(\bf{X_{Num}^{h}}+\bf{X_{Den}^{h}})])\right||\leq 0.0001$. The MIMO identification flow is shown in Algorithm 2.

After identifying the entire system with different inputs and outputs for $G_{mp}(z)$, (15) is formulated. Then $G_{mp}(z)$ is converted from discrete to continuous time domain ($G_{mp}(s)$) and transformed into a partial fraction form, such that (15) takes the form of (83).

[TABLE]

where $r_{mp}(j)$ is the residue of $G_{mp}(s)$ corresponding to the pole $p_{j}$. The residue $r_{mp}(j)$ provides information about how the mode $p_{j}$ is affected by input $u_{p}$ and how observable is it from $\Delta\omega_{m}$. This indicates that residue is a measure of joint controllability and observability index, where larger the value of residue, the stronger is the control loop.

The validity of the proposed algorithm for signal selection is verified by implementing this approach on a wind integrated two-area power system model. For this, a dominant mode of 0.6038 Hz is observed. Table III shows the comparison of proposed signal selection algorithm with the residue analysis [24] of state space matrices corresponding to the dominant mode. Since, there are two coherent groups in two-area system, two WAC are required; one in each group. Generators 1 and 2 are in one coherent group and the remaining generators are in second coherent group. It is worth noting that, as the value of residue gets larger, the stronger is the control loop to damp oscillations.

The control loop required to damp the observed mode can be derived from Table III as follows. In coherent group-1 the residue value is largest between speed of generator-3 ($\Delta\omega_{3}$) and input signal of generator-1 ($u_{1}$ ), so generator-1 is the most controllable machine. Likewise, in coherent group-2 the most observable signal is $\Delta\omega_{2}$, and generator 4 is most controllable machine. Table IV shows the WAC control loop for the case of IEEE-39 bus system with a fault on bus-14 and with four clusters. In group-1, generator 8 is the most controllable machine and the most observable signal is $\Delta\omega_{7}$, in group-2, generator 10 is the most controllable machine and the most observable signal is $\Delta\omega_{1}$, in group-3, generator 2 is the most controllable machine and the most observable signal is $\Delta\omega_{6}$, and in group-4, generator 7 is the most controllable machine and the most observable signal is $\Delta\omega_{10}$. However the residue of control loop for group-2 is very low, which can be verified from the fact that the inertia of generator-10 is very high compared to other generators and so it is the less controllable machine. Hence the controller for group-2 is ignored.

II-C WAC design architecture

The proposed WAC algorithm has three parts 1) identification of control loop transfer function, 2) Discrete-time Linear Quadratic Regulator (DLQR), and 3) state estimation using Kalman filtering. The details are discussed below. Fig. 10 shows the architecture of the proposed WAC design.

II-C1 Identification of control loop transfer function

In the proposed WAC, the state space matrices of the control loop between process input $(u_{p})$ and the process output $\Delta\omega_{m}$ is calculated using the transfer function coefficients obtained from (14) to (82). In general, the order $k$ is not known a priori, so $k$ is assumed to be a large number which is limited by number of samples $l$ and computational burden. Then, modes with negligible residues are discarded [38] and new order $p$ is identified (i.e the order of (15) is reduced from $k$ to $p$). The discrete state space matrices are formulated using (15) which can be written as in (84).

[TABLE]

where $\bf{A}$ is the state matrix, $\bf{B}$ is the input matrix, $\bf{C}$ is the output matrix, $\bf{D}$ is the feedforward matrix, $u$ represents state vector, and $y$ represents input and output respectively.

II-C2 Discrete-time linear quadratic regulator

In the proposed algorithm, the discrete-time linear quadratic regulator is used to minimize the cost function for the discrete-time linear system model as described in (84). The cost function[39] is represented as shown in (85a) where $N$ is the sample horizon, $\bf{Q}=C^{T}C$ and $R=\rho I(\rho>0)$ are the weight matrices. The cost function is minimized by calculating the optimal gain matrix $\bf{K}$, such that the optimal control sequence is given by (85b) and the optimal control gain $\bf{K}$ is given by (85c). $\bf{P}$ is the solution of the discrete time algebraic Riccati equation, which is represented as (85d) by iterating backward in time until the solution converges, initially during start of the iteration $\bf{P_{p}}=0$.

[TABLE]

II-C3 State estimation using Kalman filtering

The optimal control sequence as given by (85b) requires states of the system to be estimated, which is done using Kalman filtering technique. The state space representation in (84) obtained from RLS identification is used for this estimation. The predictor step for state vector and co-variance is given by (86) and (87) respectively.

[TABLE]

where $L$ is the covariance of state vector estimate and $Q$ is the process noise covariance. The Kalman gain factor ($G$) is calculated as shown in (88)

[TABLE]

where $H$ is the observation matrix and $R$ is the measurement noise covariance. The corrector step is given by (89) and (90).

[TABLE]

Additional details of state estimation using Kalman filtering are provided in [40].

III Simulation Results

The proposed controller is implemented on the two-area (Fig. 1) and IEEE 39-bus (Fig. 2) system interconnected with a 150MW WTG [41] with each WTG is rated at 2MW. The experimental test-bed is as shown in Fig. 11.

III-A Implementation test results using two-area system

Based on coherency grouping as shown in Table I, it can be seen that generators 1 and 2 are in one group and generators 3 and 4 are in the other group, so two WACs are required which are to be placed one in each group. The signal selection for WAC control loop is discussed in section II. The simulation results with the proposed controller are compared with the results from a system with both Exciter and PSS, and with Exciter only. Fig. 12 shows the wind speed profile used for simulations. Fig. 13 shows the active power. It is seen that (see Fig. 13) injection of the variable active power of WTG into the grid changes operating condition of the power system, which initiates synchronous generators to oscillate. At this time, a disturbance is created by initiating a fault on bus-8 at 13 sec for a duration of 0.1 sec and another disturbance is created at 33 sec by dropping load connected to bus-9 for a duration of 0.1 sec. Fig. 14 shows the relative speed of generator 2 w.r.t generator 3, Fig. 15 show the speed of generator 1 w.r.t generator 4, Fig. 16 shows the relative speed of generator 2 w.r.t generator 4, and Fig. 17 shows the WAC controlling signal to generator 1 and generator 3.

From the above results, with the WAC controlling signal (Fig. 17) to generator 1 and generators 3, it is observed that the oscillations are effectively damped (Fig. 14 to Fig. 16). From Fig. 14, it is seen that the frequency of inter-area oscillation is 0.6135Hz (approx.). Since the control loop is selected for inter-area oscillation (0.6038Hz), the inter-area oscillations are damped out effectively. From these results, it can be concluded that the proposed algorithm based on the online coherency grouping performs very effectively when compared to PSS working alone, and can be implemented online.

III-B Implementation test results on IEEE 39-bus system

Based on coherency grouping as shown in Table I, generators 4, 5, 6, 7, and 9 are in one group (group-4), generators 1, and 8 are in group-1, generators 2, 3 are in group-3, so three WAC controllers are required which are to be placed one in each group. The signal selection for WAC control loop is shown in Table IV. The simulation results with the proposed controller are compared with the results from a system with both Exciter and PSS, and with Exciter only. Fig. 18 shows the wind speed profile used for simulations. Fig. 19 shows the active power. It is seen that (see Fig. 19), injection of the variable active power of WTG into the grid changes the operating condition of the power system initiating oscillations in synchronous generators. At this time, a disturbance is created by initiating a fault on bus-14 at 5 sec for a duration of 0.1 sec and another disturbance is created at 41 sec by dropping load connected to bus 26 for a duration of 0.1 sec. Fig. 20 shows the relative speed of generator 5 w.r.t generator 2, Fig. 21 show the speed of generator 7 w.r.t generator 2, Fig. 22 shows the relative speed of generator 8 w.r.t generator 2, and Fig. 23 shows the WAC controlling signal to generator 2, 7, and 8. From Fig. 20 to Fig. 22, it can be seen that after 17sec the oscillations are damped effectively hence the WAC output is dropped.

From the above results, with the WAC controlling signal (Fig. 23) to generator 2, 7, and 8, it is observed that the oscillations are effectively damped (Fig. 20 to Fig. 22). From Fig. 20 it is seen that the frequency of inter-area oscillation is 0.683Hz (approx.). Since the control loop is selected for inter-area oscillation (0.6038Hz), the inter-area oscillations are damped out effectively. Hence it can be concluded that the proposed algorithm based on the online coherency grouping performs very effectively when compared to PSS working alone, and can be implemented online.

IV Conclusion

In this paper an online coherency based WAC is proposed. The proposed WAC algorithm is an efficient way of damping inter-area modes. The location of WAC is decided based on online coherency grouping and signal selection is based on online residue analysis. The spectral clustering and residue using MIMO system identification are computed during the simulation, so any changes in system operating conditions especially due to renewable energy sources are taken into consideration. The WAC gains calculation and state estimation are performed online, hence the controller parameters are updated based on system conditions. The efficacy of the proposed method is verified by testing it on WTG integrated two-area and IEEE 39-bus power system, models. The proof of concept illustrates that with the proposed WAC, inter-area modes can be damped much effectively.

Appendix A Parameters of Wind Turbine Generator [41]

[TABLE]

Appendix B Parameters of Wound Rotor Induction Machine

[TABLE]