Identification of coherent trajectories by modal characteristics and hierarchical agglomerative clustering

Electric Power Systems Research 158 (2018) 170–183

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Research Paper

Identification of coherent trajectories by modal characteristics and hierarchical agglomerative clustering Mario R. Arrieta Paternina a,∗ , Alejandro Zamora-Mendez b , José Ortiz-Bejar b,c , Joe H. Chow d , Juan M. Ramirez e a

National Autonomous University of Mexico (UNAM), Mexico City, Mex. 04510, Mexico Michoacan University of Saint Nicholas of Hidalgo (UMSNH), Morelia, Mich. 58030, Mexico c Centro de Investigacion e Innovacion en Tecnologias de la Informacion y Comunicacion, Aguascalientes, Ags. 20313, Mexico d Rensselaer Polytechnic Institute, Troy, NY 12180, USA e CINVESTAV, Guadalajara, Jal. 45019, Mexico b

a r t i c l e

i n f o

Article history: Received 7 August 2017 Received in revised form 30 November 2017 Accepted 29 December 2017 Keywords: Coherency Hierarchical agglomerative clustering Inter-area modes Modal identification Slow coherency

a b s t r a c t This paper introduces a novel method to identify coherent generators using the inter-area modal characteristics of power systems. The key idea is to extract the inter-area modes from the simulated data and then to apply a clustering strategy. Thus, the proposed method consists of extracting the phase of the oscillatory modes via a modal identification technique and applying a hierarchical agglomerative clustering technique together with the Elbow’s method to gather the phases of each mode, enabling to provide coherent trajectories of generators. The proposed method uses a Taylor-Fourier filtering strategy to remove noises and nonlinearity in the time evolution of coherent generators. Simulated signals with noise added are used for assessing the proposition. Results corroborate the proposed strategy for identifying coherent trajectories in large-scale power systems. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Coherent trajectories in power systems are useful for controlled islanding [1–3] and wide-area control [4,5]. They arise due to the synchronous rotating machines connected into the power network, because they behave as coupled swing dynamics. This concept is widely known as coherency. There are two main types of methods for identifying coherent generators [6,7], classified as: (i) model-based methods; and (ii) measurement-based methods. The former is generally based on the linearized state-space model, where methods such as slow coherency grouping (using the singular perturbations technique

Abbreviations: WAMS, wide-area monitoring systems; TFF, Taylor Fourier filters; FT, Fourier transform; HHT, Hilbert-Huang transform; PA, Prony analysis; ERA, eigensystem realization algorithm; MP, matrix pencil; TKEO, Teager-Kaiser energy operator; TFT, Taylor-Fourier transform; HACA, hierarchical agglomerative clustering algorithm; NE, New England; NPCC, northeast power coordinating council; SC, slow coherency; MB, measurements-based; FLOPS, floating point operations; SSE, sum of squared errors; TBM, tolerance-based method. ∗ Corresponding author. E-mail addresses: [email protected] (M.R. Arrieta Paternina), [email protected] (A. Zamora-Mendez), [email protected] (J. Ortiz-Bejar), [email protected] (J.H. Chow), [email protected] (J.M. Ramirez). https://doi.org/10.1016/j.epsr.2017.12.029 0378-7796/© 2018 Elsevier B.V. All rights reserved.

to display the time-scale separation of the inter-area modes and local modes) [8], tight slow coherency grouping [8], Zaborszky’s clustering technique [9], weak link, and Tolerance-Based Method (TBM) technique [10,11], are used to analize offline coherency studies. The TBM finds coherent generators using the right eigenvectors from a set of user-specified modes; where generators are considered coherent if the phase of their right eigenvector entries, relative to a common reference entry, are within a coherency threshold. Also, the chosen modes should be inter-area modes with large observability in the entire system. Meanwhile, the latter deals with the measured or simulated response of the power system taking advantage of Wide-Area Monitoring Systems (WAMS) composed of synchrophasors. Several methods have been proposed, including recent techniques based on Koopman modes [12,13], graph theory [1], hierarchical clustering [14], correlation coefficients [15], K-harmonic means clustering [16], independent component analysis (ICA) [17–19], support vector clustering [20,21], and frequency deviation signals (FDS) [22]. These techniques address the coherency problem using polluted signals, handling noise immunity between 20 and 50 dB for FDS and ICA [22]. This paper tackles the inter-area modes identification-based coherency from the simulated response of power systems, with the ability of identifying modal parameters even in presence of noise by the Taylor-Fourier

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

filters (TFF) introduced in power systems in [23,24] and in phasor estimation in [25]; which allow to extract the inter-area modes, as exhibits in [23,24]. Likewise, measurement-based methods allow coherency identification under changes in system operating condition and network configuration, because the grouping of coherent generators may vary. Thus, transient conditions imply that one coherent group may be split into smaller groups, or on the contrary, multiple groups may be joined into a bigger coherent group. The proposed approach is motivated by the coherency concept defined by Podmore in [26], where a clustering algorithm is used to process the approximate swing curves obtained by linear simulation and thereby determine the coherent generators, so that the coherency is based on similar dynamic behavior between two machines. The approach here is driven by the advance in modal extraction techniques based on the power system measured response, which are assumed reliable in [27–30]. Other approaches have been proposed for identifying oscillatory modes based on measurements. Some relevant techniques are discussed in [31,32]. Among these, the following are highlighted: Fourier transform (FT); Hilbert-Huang transform (HHT); Prony analysis (PA); eigensystem realization algorithm (ERA); matrix pencil (MP); and Teager-Kaiser energy operator (TKEO). In reference [33], comparisons are accomplished exhibiting the Taylor-Fourier transform (TFT) performance with respect to PA, ERA, and MP, unveiling a suitable performance on its applicability for identifying electromechanical modes. In order to get the coherent groups, the modal characteristics estimated by TFF and a clustering technique are used. As stated, previous works use the topological information of the grid for grouping the coherent trajectories. In this paper, the modal information is used for identifying the coherent groups, i.e., no information about the network structure is needed. The hypothesis is validated using an exhaustive clustering technique known as the hierarchical agglomerative clustering algorithm (HACA) [34,35], selecting the optimal number of clusters by the Elbow’s method presented in [36]. The applicability of this proposition has been assessed in two large-scale power systems: New England (NE) power system [37] and the Northeast Power Coordinating Council (NPCC) system [8]. Comparisons with the classical slow-coherency method are provided in order to corroborate the proposition. According to the above-mentioned, the major contributions of the paper are summarized as follows: 1. A novel strategy for determining inter-area-mode-based coherent trajectories including high noise tolerance is proposed by a filtering approach. 2. This paper introduces the coherent trajectory grouping based on modal characteristics from the Taylor-Fourier filters estimates. 3. A hierarchical agglomerative clustering technique and the Elbow’s method are used to identify coherency. The paper is organized as follows. Section 2 presents the technique of identifying oscillatory modes based on Taylor-Fourier filters. Then, the hierarchical agglomerative clustering algorithm and Elbow’s method are discussed in Section 3. In Section 4, the computational complexity is presented. Finally, the method is applied to two large-scale power systems in Section 5.

where p(t) = a(t)ejϕ(t) is known as a dynamic phasor and f1 is the fundamental frequency, so that the phasor can be approximated by the Kth Taylor polynomial centered at t0 as (t − t0 )K K!

˙ 0 )(t − t0 ) + · · · + pK (t0 ) p(t0 ) + p(t

pK (t) =

(2)

t0 − T/2 ≤ t ≤ t0 + T/2 To build the Taylor subspace, the signal model in (1) is discretized by t = ns Ts , where Ts = 1/(Nf1 ) is the sampling period, and N represents the samples per fundamental cycle. Thus, the Taylor subspace is constructed by T = [t0n

t2n /2! · · ·

tn

tKn /K!]

(3)

where tn is a diagonal matrix whose components are given by tn =− (K + 1)Ts (ns /2) to (K + 1)Ts (ns /2), ns corresponding to each sample of the Taylor’s interpolating polynomial at each sampling time (Ts ). The Taylor-Fourier (TF) subspace in [33] is shaped as

⎡

t0n1

⎢ 0 ⎢ tn2 ⎢ B=⎢ ⎢ .. ⎣. t0nC

t1n1

tKn1

···

⎤⎡

WN

.. .

⎥⎢ · · · tKn2 ⎥ ⎢ 0 ⎥⎢ ⎥ . .. ⎥ ⎢ .. . . ⎦ ⎣ ..

t1nC

···

t1n2

tKnC

0

0

···

0

WN

···

0

.. .

..

.. .

0

···

.

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(4)

WN

where C = K + 1 is the number of cycles, which corresponds to the slower swing modes into the oscillating signal, and WN is the h = ej2h/N in each Fourier matrix with harmonic phase factors ωN vector h = 0, . . ., N − 1. Thus, the Fourier matrix is

⎡1 1 ⎢ 1 ωN ⎢ ⎢ 2 ⎢ WN = ⎢ 1 ωN ⎢ ⎢ .. .. ⎣. .

(N−1)

1 ωN

1

···

2 ωN

···

4 ωN

.. . 2(N−1)

ωN

⎤

1

⎥ ⎥ 2(N−1) ⎥ ⎥ · · · ωN ⎥ ⎥ .. ⎥ .. . . ⎦ ···

(N−1)

ωN

(5)

(N−1)2

ωN

Note that in (4), the vectors of the Fourier matrix are harmonic modulators of the Taylor terms included in a Kth Taylor polynomial, K > 0. Once the Taylor-Fourier subspace is defined, the synthesis equation is established by ˆ sˆ = B,

(6)

where matrix B is known as the Taylor-Fourier matrix, sˆ is the estimated signal, and vector ˆ contains the estimated TF coefficients and their derivatives according to the signal model (1) proposed in [23,38]. By taking the least-squares technique as the Taylor-Fourier coefficients in [39,40]. The error between the input signal s and its approximated Kth Taylor interpolating polynomial Bˆ is defined as ˆ e = s − B.

(7)

Then, the best estimates of ˆ are those obtained by solving the normal equations in (8).

2. Identification of modal characteristics The Taylor-Fourier transform (TFT) for extracting modal information from oscillating signals in power systems is described here. The TFT is integrated by Taylor and Fourier subspaces (for more details, see [24,33]), allowing a choice on the number of non-zero Taylor terms used for the continuous-time signal defined by s(t) = Re{p(t)ej2f1 t } = a(t) cos(2f1 t + ϕ(t))

171

(1)

BH Bˆ = BH s.

(8)

By solving (8), the best parameters are attained, in the sense that they minimize the sum of the squared errors in (7). Thus, no iterative procedure is required. Then, the TF estimates are given by −1 ˆ = [BH B] BH s = B† s,

(9)

172

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

To analyze multiple signals using the TF filters, the synthesis and analysis equations, (6) and (9), respectively, are expanded to form a set of M signals as [24] [ˆs1

sˆ 2

···

sˆ m ] = B[ˆ 1

[ˆ 1

ˆ 2

···

ˆ m ] = B† [s1

ˆ 2 s2

··· ···

ˆ m ], sm ],

(13) (14)

where sm (n), m = 1, . . ., M, share a common set of modes. B stands for the Taylor-Fourier matrix in (4), i.e., the filter bank is assumed equal for all signals (a single filter bank is designed). The parameter estimation process for multiple channels is evaluated over the estimates provided by the analysis equation (14), that is, each element in {ˆ 1 ˆ 2 · · · ˆ m } will render amplitude ˆ estimates in (11) and (12). (ˆa) and phase () Notice that the signal model in (1) and matrix B in (4) allow the establishment of linear (K = 1) and nonlinear (K = 3) approaches for the Taylor-Fourier filter banks, depending on the Taylor’s interpolating polynomial used. Both linear and nonlinear analyses will be considered in this investigation. In the next section, the modal phase estimates obtained by TFF are used by a hierarchical agglomerative clustering strategy and Elbow’s method to form coherent groups. The reliability of the TF filters lies on the basis to shape the TF subspace and the small errors guaranteed by including Taylor’s interpolating polynomial terms greater than zero which in turn modulate the Fourier coefficients, and by applying the least-square solution to estimate the best Taylor-Fourier coefficients [38].

Fig. 1. Frequency response of Taylor-Fourier filters for K = 3.

3. Hierarchical agglomerative clustering algorithm and Elbow’s method The inputs to the clustering algorithm are the phase estimates in (12), obtained from the oscillatory modes identified by TFT at each time stamp. The proposed algorithm is implemented in two stages: an exhaustive centroid-based HACA and the determination of the optimal number of clusters by the Elbow’s method.

Fig. 2. Modal characteristics extraction by TF filters.

ˆ˙ ˆ¨ whose phasor’s estimate is ˆ and its time derivatives are , , . . ., ˆ (K) , containing their complex conjugates as well. ˆ = [ˆ

ˆ ˙

ˆ¨

···

T

ˆ (K) ] .

(10)

Thus, all the Taylor-Fourier estimates in (10) are analytically deduced and computed based on the least-square process. The columns of B contain the basis vectors for the TF subspace, and its pseudoinverse B† , whose rows correspond to the impulse responses of the TF filters in Fig. 1, where BH stands for the Hermitian matrix of B. Such filters possess maximally flat pass-bands at the frequencies of concern fi and zero flat stop-bands [41], Fig. 1. Once the TF subspace has been established, the modal parameters’ estimation is accomplished by the analysis equation (9), i.e., the amplitude and phase for the jth mode are estimated by aˆ j (t) = 2|ˆ j | ˆ j (t) = ∠ˆ j . 

(11) (12)

To assess the performance of the modal identification, an oscillating signal with 15 dB of noise is processed by TFT in Fig. 2a. The signal has two oscillatory modes in its Fourier spectrum, f1 and f2 , as exhibited in Fig. 2b. When the oscillating signal is processed by TFT, signal decompositions are estimated by (6) and (9), so that the time evolution per mode s1,2 (t) is achieved in Fig. 2c. Then, the modal ˆ 1,2 (t) in Fig. 2d and e are estimated by (11) parameters aˆ 1,2 (t) and  and (12), respectively. The Taylor-Fourier estimates are computed once the first analysis window passes by the filters, corresponding to the center of the window.

3.1. Hierarchical agglomerative clustering algorithm In the agglomerative strategy, each element in the data set is a cluster and the two nearest clusters are joined at each iteration, which is repeated until there is a single cluster with all the elements. Similarity measurements, as defined in [34,42], are used by the agglomerative approach to determine which pair of clusters are joined at each iteration. The hierarchical agglomerative algorithm requires an error measurement and a proximity function. Here, the sum of squared errors (SSE) is used for measuring the error. The SSE for a set of points X = {x1 , x2 , . . ., xn } with xi ∈ Rn , ∀i = 1, 2, . . ., n, and for any given set C of k clusters {C1 , C2 , . . ., Ck }, is defined as SSE(k) =

mi k  

d(xi,j , C¯ i )

2

(15)

i=1 j=1 2 where d(xi,j , C¯ i ) is the distance from the point xi,j ∈ Ci to the centroid of cluster C¯ i (computed as in [43] for equal masses), and mi is the number of elements in cluster Ci . The proximity function used in this paper is based on the centroid method. This criterion uses the minimum Euclidean distance among centroids in choosing the clusters to be joined. In the proposed HACA strategy, an initial set of clusters Cˆ n is given by

Cˆ n = {C1 , C2 , . . ., Cp , . . ., Cn }

(16)

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

173

Fig. 3. Finding the Elbow’s point.

where each Ci contains only one element, whose similarity measurement is quantified by the Euclidean distance d2 (C¯ 1 , C¯ 2 ) = (C¯ 1 − C¯ 2 ) · (C¯ 1 − C¯ 2 )

(17)

Thus, the hierarchical agglomerative iterative process is accomplished by starting at k = n, so that, {k ∈ N|n ≥ k ≥ 1}, and proceeds as follows:

i = 1, . . ., n,

The farthest point k is calculated using the straight line equation L : y = ax + b

1. Calculate and store SSE(k). 2. Build the distance matrix D as Di,j = d(C¯ i , C¯ j ),

Fig. 4. Flowchart of the proposal. TFT (gray box) and HACA (black box).

(18)

where j = 1, . . ., n

3. Find the pair of clusters with centroids (C¯ p , C¯ q ), such that Dp,q is the minimum distance in D. 4. Create a cluster Cpq = Cp ∪ Cq , where p, q are returned from (p, q) = mini,j Di,j . 5. Create a new set of clusters, such that: Cˆ k−1 = ((Cˆ k Cp )Cq ) ∪ Cpq . where stands for the set subtraction operator. 6. Update matrix D including the centroid C¯ pq for the new cluster Cpq . The algorithm generates n possible ways to cluster elements in X, starting with n clusters and repeating steps 1–6 iteratively until one set of clusters remains, where n «2n and 2n is the number of all possible ways to group n elements. 3.2. Elbow’s selection rule The Elbow’s method is a heuristic technique used to determine the optimal number of clusters in a data set. It is based on the interpretation of the SSE and is analyzed as a function of the number of clusters. Notice that the SSE(k) has a minimum value when the number of clusters is equal to the number of elements in C, and a maximum value when all the elements belong to a single cluster. SSE(k) is monotonically decreasing, as can be seen in Fig. 3, which is obtained by plotting k vs SSE(k). Such behavior results from the agglomerative strategy, because SSE(k) ≥ SSE(k + 1) [44]. 3.3. Finding the optimal number of clusters After calculating all the groups from n to 1, the optimal number of clusters is found by using a geometric criterion in two steps. First, it determines the straight line L between points (1, SSE(1)) and (n, SSE(n)), see Fig. 3. Second, it locates the farthest point k from the line to SSE(k). At completion, the Elbow’s selection rule returns k and Ck , where Ck is the Cp whose distance is the farthest from L to SSE(k).

a=

SSE(1) , 1−n

b=

nSSE(1) 1−n

(19)

The distance Di is calculated as Di = d(SSE(i), L) =

|a(i) − SSE(i) + b|



a2 + 1

, ∀i = 1, . . ., n.

(20)

Thus, the Elbow’s value will be given by k = maxi {Di }

(21)

4. Proposal and computational complexity 4.1. Proposed strategy based on TFT and HACA Fig. 4 shows the flowchart to identify the coherent trajectories by modal characteristics and HACA. The modal characteristics (dotted gray box in Fig. 4) are drawn employing expressions (3)–(12), where the input parameters for designing the filter bank are indicated. The filter bank is designed by (4) using the rows of its pseudoinverse. Then, the modal phases are computed by (10) through the Taylor-Fourier coefficients obtained by means of the analysis equation (7), and they are the entries to the HACA in order to find the optimal coherent trajectories. The clustering process (dotted black box in Fig. 4) starts building a set of clusters Cˆ n , where each single element Ci corresponds to the modal phases computed for the ith generator. Then, a distance matrix D is computed through the clusters centroids, and a vector SSE stores the error for each cluster configuration Cˆ k . Then, it computes a new cluster Cpq (starting from k = n − 1) by joining two clusters (Cp and Cq ), whose centroids are the nearest according to D. Then, it builds a new cluster of |k − 1| elements by subtracting Cp and Cq from Cˆ k and joining Cpq ; it also computes an error value for the new cluster, i.e. SSE(k − 1). Afterwards, it updates matrix D by removing the distances associated to Cp and Cq , and adding the distance to the new cluster Cpq . Finally, it gets the optimal number of clusters k by applying the Elbow’s rule. 4.2. Computational complexity The computational complexity and TFT burden have been discussed and compared with fast Fourier transform (FFT) in [24,30,33,38,45]. The number of operations in TFT can be reduced

174

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

white Gaussian noise (AWGN). To identify the low-frequency components associated with the electromechanical modes, TFT with a sliding window advancing one sample at a time is used. The performance of the strategy for identifying coherent groups is corroborated using two scenarios with the following disturbances: 1. Scenario 1. New England system: 3-phase short-circuit at Bus 32 lasting 3 cycles, cleared by tripping the line between Buses 32 and 33. 2. Scenario 2. NPCC system: 3-phase short-circuit at Bus 7 (Medway) lasting 3 cycles, which is cleared by tripping the line from Medway to Sherman Road (Bus 6).

Fig. 5. Test system: 16-machine and 68-bus New England Power System.

by a shorter window in order to provide good dynamic estimates, that is, with K = 1 and C = 2. Further computation reduction is achieved due to the fact that the Taylor-Fourier matrix only contains the frequencies of concern [24,30,45]. That is, once B is computed by (4), only the columns corresponding to the frequencies of concern (Cfi ) are taken into account, reducing the size of B from CNxCN to CNxCfi , without affect the feasibility of the TFT method. It is noteworthy that matrix B and its pseudoinverse are calculated just once. With respect to the hierarchical agglomerative clustering algorithm, the total computation burden is (3M2 ) + (2M2 logM) + (M logM) [35]. But, according to the theory of asymptotic analysis only the dominant terms are taken into consideration, since the lower-order terms are relatively insignificant for large inputs [46]. Thus, the complexity is given by 2M2 logM. The proposed methodology is based on a heuristic technique, finding the optimal number of clusters k and minimizing the SSE in the Euclidean space. It has been shown that when the decision relies on minimizing SSE and the parameter k is fixed, the clustering process is conceived as NP-hard problem. This optimization is even worst when k is not defined, such that the complexity may become NP-complete [47]. If the brute force approach is chosen, all the possible cluster configurations should be evaluated, that is, for

N  1 ≤ k ≤ n, i.e. Nk = O(2n ) combinations for each value of k=1 k. Since the implemented architecture uses an exhaustive strategy for handling the search, only a set of n possible clusters configurations is evaluated, resulting in a complexity order reduction n «2n for exploring the n configurations. On the other hand, the real-time implementation of the proposed methodology may be achievable, if the computing time required by the proposal (TFT+HACA) is executed within a sample time (Ts ) defined by the sampled frequency (60 Hz), i.e. Ts = 16.66 ms. It is noteworthy that most of the commercial PMUs have sampling rate at 20, 30, 50, and 60 Hz, that is, the time-stamps may still become higher than 16.66 ms, increasing the execution interval of a significant manner for the feasibility of the implementation. 5. Test cases and results The applicability of the proposal has been tested using two large-scale power systems: the New England (NE) system in Fig. 5 and the Northeast Power Coordinating Council (NPCC) system in Fig. 6. Transient stability simulations for disturbances are performed using the Power System Toolbox (PST) [48]. Then, time-frequency analysis is carried out on the generator rotor speed variables, sampled at a rate of 60 Hz, with 15 dB of added

Both scenarios will allow the illustration of the transient behavior when the network’s topology changes after a disturbance takes place, that is, transmission lines are tripped after a fault. Likewise, the oscillatory modes are excited with severe disturbances, so that their characteristics are time varying. The performance of Taylor-Fourier filters is investigated with respect to the Taylor’s interpolating polynomial with K = 1 and K = 3, that is, both linear and nonlinear approaches in the Taylor subspace in (4) are considered. Note the flat top of the filters with K = 3 in Fig. 1. Then, the signals from the test power systems are processed by means of (13) and (14), which are the impulse responses of the TF filters for identifying the oscillating modes. Similarly, the TF coefficients are obtained by the analysis equation. Subsequently, the modal trajectories with its corresponding phase estimate (t) in (12) are computed. The TF filters are able to capture and extract the coherent modes as the oscillations are evolving. Prior to applying the hierarchical agglomerative clustering technique at each time-stamp for the phase corresponding to the jth mode of the ith signal, the list of points X in Section 3.1 needs to be defined. For this purpose, phases with the same time-stamps provided by TFT make up the initial points for HACA, and initially s1 (t), s2 (t), . . ., si (t)}, ∀j = 1, 2, . . ., N ,∀i ∈ M. grouped as X ={m m mj mj j Nm stands for the number of modes. For the test systems, each signal is associated with one generator, so that the coherent grouping may be defined. All generators are assumed as synchronous machines in order to evaluate the proposed methodology. 5.1. Coherent trajectories in New England power system The identification of the coherent trajectories performed on the NE system is based on the dynamic simulation in Scenario 1. The generator rotor speeds with noise added are depicted in Fig. 7a, and their corresponding Fourier spectra are shown in Fig. 7b. The frequency components at 3 inter-area oscillatory modes of 0.39, 0.64, and 0.78 Hz are used to design a TF-filter bank with impulse responses (corresponding to the coefficients at each row in matrix B† , using Fb = 0.13 Hz) as displayed in Fig. 8 [24,33], so that the central frequencies are close to the frequency components of concern. Then, the TF estimates are obtained for modal trajectories in Fig. 9 by (6), amplitudes in Fig. 10 by (11), and phases in Fig. 11 by (12), using TF-filters with K = 1 and K = 3, respectively. Then, the hierarchical agglomerative clustering algorithm is applied on the phase estimates of the oscillatory modes exhibited in Fig. 11, for the estimates provided by TF-filters with K = 1 and K = 3, respectively. It takes an initial set of clusters into account, corresponding to the modal phase estimates for the generators rotor speed signals, that is, X in Section 3.1. The HACA identifies the coherent generators once the TF-filters deliver a new estimate, as seen in Figs. 12–14. Note that a coherent group can be identified at each time step despite location of the disturbance. Figs. 13 and 14 correspond to short-circuit applied at bus 32 and 24, in both cases adjacent lines to these buses are tripped in order

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

175

Fig. 6. Coherent trajectories in NPCC power grid. inter-area modes-based (MB) using TFT for K = 3 (dotted line), and slow-coherency (SC) (continuous line).

Fig. 7. Waveforms and their spectra for 16 signals. (a) Generator rotor speeds and (b) Fourier spectra.

Fig. 8. Frequency response of the TF filters for the 3 slowest inter-area modes using K = 1 (dotted line) and K = 3 (continuous line).

to clear the fault. This is corroborated with the time evolution of the clusters in Figs. 12–14, where the coherent trajectories are split into 5 areas consisting of generators: G1-G9, G10-G13, G14, G15, and G16. Notice that the same coherent areas are kept for all times in Figs. 12–14, meaning that the TF-filters are not influenced by the noisy condition for K = 1 and K = 3. Consequently the attained

grouping by both approaches is suitable, despite the length of the analysis window for providing the TF estimates. For this scenario, the units assigned to each area in Fig. 5 will be composed by equivalents, whose equivalent apparent power and equivalent inertia per area according to the proposed approach are illustrated in Table 1. It is noteworthy to remark that Areas 3 to

176

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

Fig. 9. Oscillatory modes. (a) Mode 1, (b) mode 2, and (c) mode 3.

Table 1 Relationship between equivalent inertia and apparent power per Area for Scenario 1.

Area 1 Area 2 Area 3 Area 4 Area 5

Rated MVA

Equivalent inertia [s]

Equivalent Apparent power [pu]

100 100 100 100 100

282.6 647.5 300.0 300.0 450.0

53.4476 65.6917 17.8608 10.0225 40.2637

5 are formed by equivalent generators in NE power grid. Meanwhile, Areas 1 and 2 will be represented by single generators, whose parameters can be determined by Podmore aggregation, inertial aggregation or the slow-coherency aggregation in [8]. In order to compare the estimated phases with the phases of the right eigenvectors components, the small signal analysis (SSA) is performed in New England power system in order to obtain the mode shapes corresponding to the swing modes at 0.39 Hz, 0.64 Hz, and 0.78 Hz, respectively, as displayed in Fig. 15. Mode 1 in Fig. 15a indicates that Areas 1 (G1–G9) and 2 (G10–G13) swing against Areas 3 (G14), 4 (G15), and 5 (G16) of Fig. 5. Mode 2 in Fig. 15b shows that generators in Areas 1 (G1–G9) and 2 (G10–G13) oscillate with the opposite phase. Meanwhile, Mode 3 at 0.78 Hz in Fig. 15c that Areas 3 (G14) and 4 (G15) swing against Area 5 (G16). Regarding the relationship between equivalent inertia and apparent power per Area for Scenario 1 (see Table 1), it is remarkable that the aggregated generators in the reduced-NE model (see Fig. 16) retain a similar coherent dynamic in contrast to the full

system in Fig. 15. The aggregation process is performed by using Podmore method in [8].

5.2. Low-frequency coherent trajectories in NPCC system The proposed method is also applied to the larger NPCC power system [8] in Fig. 6. The machine rotor speeds subject to the Scenario 2 disturbance and with a 15 dB noise added are displayed in Fig. 17a. The Fourier spectra of these signals are exhibited in Fig. 17b, showing the 4 slowest frequency components corresponding to the inter-area modes for the full model (Table 3.4 in reference [8]). The impulse response of the TF filters in Fig. 18 is used for identifying the trajectories for the oscillating modes. Likewise the estimate of their phases by (12) are shown in Fig. 19. It is interesting that the modal phase estimates in Fig. 19a contain more variations for case K = 1, advocating the use of the Taylor interpolating polynomial with K = 3. The results from K = 3 contains less variations under noisy conditions according to the impulse response obtained in Fig. 18, as observed in Fig. 19b. Once the phases estimates are available, the HACA in Section 3.1 starts the clustering process providing the coherent trajectories at each time interval. Several time instants are presented with the proposed generators grouping in Table 2, for the modal estimates using TF-filters. For K = 1, the generators are split into 12 areas, and for K = 3, they are split into 11 areas. Note that the main difference is Generator 38 corresponding to Area 7 in case K = 1, which is grouped in Area 6 for K = 3. Fig. 6 shows the clusters at t = 39.91 s, t = 41.48 s, and t = 49.23 s, and compares with the slow coherency approach.

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

177

Fig. 10. Amplitude estimates in Scenario 1 for: (a) mode 1; (b) mode 2; and (c) mode 3.

According to the results attained in Scenario 2 and assuming the clusters at t = 39.91 s, t = 41.48 s, and t = 49.23 s, the inertial aggregation will render equivalents, as exhibited in Table 3.

5.3. Discussion The Taylor-Fourier subspace allows the construction of bandpass filters in Figs. 8 and 18 for K = 1 and K = 3, which work with different window lengths. These filters introduce delays into the estimation process executed by TFT, such that the first phase estimates in Scenario 1 are provided around 7.7 s for K = 1 (Fig. 11a) and 16 sec for K = 3 (see Fig. 11b). For Scenario 1, both linear and nonlinear approaches yield identical results in Scenario 1, as corroborated in Figs. 12 and 13, respectively. Whereas, the filters in Scenario 2 have a delay corresponding to 15 s (Fig. 19a) and 31 s (Fig. 19b) for K = 1 and K = 3, respectively. On the other hand, the computing time for signal processing, parameters estimation, and generator clustering is evaluated for multi-channels strategy stated in [24] (16-oscillating signals in Scenario 1, and 48-oscillating signals in Scenario 2). Thus, the TFT’s computing time in both scenarios comprises the processing stage (analysis and synthesis equations) and post-processing stage (parameter estimation), regarding that each filter bank is only designed once, whereas, the computing time for the clustering stage corresponds to get M clusters and find the optimal using the HACA and Elbow’s rule, respectively. Table 4 exhibits the computing time and the number of Floating Point Operations (FLOPS)

for both scenarios, running one thousand times the TFT and hierarchical agglomerative process on a processor at 2.7 GHz with 4 GB RAM. As mentioned in Section IV.B, the implementation of the algorithm is achievable, if all operations required by itself are executed within a sample time defined by the sampling rate (60 Hz). Thus, the implementation of the proposed algorithm depends on the execution time of all its mathematical operations with respect to the sampled frequency. In this paper, the sampling time is 16.66 ms and the maximum number of FLOPS in scenarios 1 and 2 becomes 1,425,474 and 11,380,644, respectively, this can be noted in Table 4. There, exists digital processors able to work with 0.83 ns of instruction cycle time [49,50], which corresponds to a processing frequency of 1.2 GHz (PowerPC 970MP RISC Microprocessor) [51], this implies that the worst time in order to execute all the mathematical operations by the proposed algorithm is 9.48 ms, corresponding to M = 48, Nm = 4, and TFT using K = 3. Therefore, the time exhibited by the proposal becomes lower than that of the sampling time. Thus, a microprocessor executes the required calculations within the sampling interval. For the Slow-Coherency algorithm, the computing time corresponds to 133 ms, which is greater than 16.66 ms. In Scenario 2, it is noteworthy to remark that the modal phase estimates corresponding to K = 1 in Fig. 19a are more affected by noise than the estimates for K = 3 in Figs. 19b and 10b, highlighting the estimates with less variations for K = 3. Thus a different grouping is generated by K = 3 (see Table 2 and Fig. 6). Finally, it infers that

178

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

Fig. 11. Phase estimates in Scenario 1 for: (a) mode 1, (b) mode 2, and (c) mode 3.

Fig. 12. Evolution of coherent groups in Scenario 1, fault at Bus 32, for estimates using K = 1 at: (a) 7.7000 s, (b) 9.3667 s, (c) 11.0333 s, (d) 12.7000 s, (e) 16.0333 s, and (f) 19.3667 s.

the bandpass filters with flat top denoising more successfully to provide the better phase estimates for the clustering process. In both scenarios, the coherency obtained could be time varying. The nonlinear approach for K = 3 yields more consistent clusters as the result is less affected by noise. Note that in Table 2 (Scenario 2, NPCC system), the area’s participants could become different; for instance, at t = 15.38 s G31, G33–G35 are in Area 5 while at t = 30.76 s

G33–G35 are in Area 6. This fact is the result of energy exchange among generators (the capability of a system to perform work upon another depends on the energy state that both systems exhibit relatively to each other), which affects the physical properties due to the influence the transfer rate of energy among them. Likewise, the disturbance magnitude, location and proximity to the fault, also the system structure and parameters make generators to respond with

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

179

G15 G16 G14

G1 - G9 G10 - G13

Mode 3

G14 G10 - G13

G14

G10 - G13

5

0

Mode 2 -5 -5

0

5

Mode 1 G15

G1 - G9

0

2 0 -2 5

G16

G14

G10 - G13

5

0

Mode 1 G15 G16

G1 - G9

G14

G10 - G13

0

Mode 2 -5 -5

Mode 1

(e)

G14

G10 - G13

(d)

G16

0

0

Mode 2 -5 -5

G15 G1 - G9

2 0 -2 5

Mode 1

(c)

2 0 -2 5

5

0

Mode 2 -5 -5

G1 - G9

(b)

G1 - G9

0

G15 G16

Mode 2 -2 -5

G15

G16

2 0 -2 2

Mode 1

(a)

2 0 -2 5

5

0

Mode 3

0

Mode 2 -2 -5

Mode 3

Mode 3

5 0 -5 2

Mode 3

Mode 3

Fig. 13. Evolution of coherent groups in Scenario 1, fault at Bus 32, for estimates using K = 3 at: (a) t = 15.42 s, (b) t = 17.08 s, (c) t = 18.75 s, (d) t = 20.42 s, (e) t = 22.91 s, and (f) t = 24.56 s.

(f)

5

0

Mode 1

Fig. 14. Evolution of coherent groups in Scenario 1, fault at Bus 24, for estimates using K = 3 at: (a) t = 15.42 s, (b) t = 17.08 s, (c) t = 18.75 s, (d) t = 20.42 s, (e) t = 22.91 s, and (f) t = 24.56 s.

Fig. 15. Mode shapes in NE power grid for the inter-area modes.

120

90 1

60 0.5 30

150

Area 2

Area 4

180 Area 1

240

270 (a)

300

60 0.5 30

150

120 150

330

180

0 Area 1

210 240

270 (b)

300

60 0.5 30

180

0 Area 4

330

90 1

Area 1 Area 2

Area 2

0

210

90 1

120

210 240

Area 5 Area 3

270

330

300

(c)

Fig. 16. Mode shapes in Reduced-NE system to 5 areas for the inter-area modes at: (a) 0.39 Hz; (b) 0.64 Hz; (c) and 0.78 Hz.

180

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

Fig. 17. Waveforms and their spectra for 48 signals. (a) Generator rotor speeds and (b) Fourier spectra.

Fig. 18. Frequency response of the TF filters for the 4 slowest inter-area modes using K = 1 (dotted line) and K = 3 (continuous line).

Fig. 19. Time evolution of the modes’ phase estimates in Scenario 2. (a) Mode 1, (b) mode 2, (c) mode 3, and (d) mode 4.

a high or low intensity under disturbances. According to the equation of rotational motion of each generator, the damping factor may lead to rotor mechanical oscillations to become time varying. Thus, machines coherently oscillate during the first swing may not do so during the subsequent swings. From a macro perspective of energy transformation (energy conservation between injected mechanical kinetic energy and network electromagnetic energy), power systems transient stability is a reflection of unbalanced energy absorption within the power network. Therefore, analysis of the transient energy flow is critical for

determining stability, and is primarily composed of the exchange between kinetic and potential energy. Ref. [52] shows that when a large disturbance or fault takes place, the instability in a power system is generally between the critical machines cluster (CMC) and the non-critical machines cluster (NMC). System stability mainly depends on the relative energy between two generators in different clusters. More pragmatically, under some specific operational conditions, the generators in large-scale power systems might lose synchronism and as a result are clustered into more than two groups

– – – – – – G47 G47 G47 G47 G47 G47 G43–G46 G43–G46 G43–G46 G43–G46 G43–G46 G43–G46 G41–G42 G41–G42 G41,G42 G41,G42 G41,G42 G41,G42

G48 G48 G48 G48 G48 G48

G48 G48 G48 G48 G48 G48 G43–G47 G43–G46 G43–G46 G43–G46 G43–G46 G43–G47 G41,G42 G41,G42 G41,G42 G41,G42 G41,G42 G41,G42 G39 G39 G39 G39 G39 G39

– G47 G47 G47 G47 –

Area 10

– G38 G38 G38 G38 G38

G39 G39 G39 G39 G39 G39

G32,G37,G38,G40 G32–G35,G37,G40 G32–G35,G37,G40 G32–G35,G37,G40 G32–G35,G37,G40 G32–G35,G37,G40

G32–G38, G40 G37–G38, G40 G37,G38,G40 G32, G37,G38,G40 G32, G37,G38,G40 G32, G37,G38,G40

Area 8 Area 7 Area 6

Area 1 Area 2 Area 3 Area 4 Area 5 Area 6 Area 7 Area 8 Area 9 Area 10 Area 11

Rated MVA

Equivalent inertia [s]

Equivalent Apparent power [pu]

100 100 100 100 100 100 100 100 100 100 100

258.20 349.92 371.40 55.00 233.70 221.00 999.99 223.00 1182.99 3.50 999.99

54.7952 68.6225 61.1448 12.6461 41.8036 37.0060 7.0136 27.1348 55.2483 1.9929 1.1722

Table 4 Computing time and number of FLOPS. M

Nm

Elapsed time in ms Max.

16 48

3 4

16 48

3 4

16 48

3 4

HACA 3.97 5.99 TFT, K = 1 0.56 0.67 TFT, K = 3 0.65 2.90

G31,G33–G35 G31 G31 G31 G31 G31

– G32–G36 G32–G36 G33–G36 G33–G36 G33–G36

G27–G30 G27–G30 G27–G30 G27–G30 G27–G30 G27–G30

G31 G31 G31 G31 G31 G31

G11–G26,G36 G11–G26,G36 G11–G26,G36 G11–G26,G36 G11–G26,G36 G11–G26,G36

– – G27–G30 G27–G30 G27–G30 G27–G30

Coherent trajectories for K = 1 G1–G9 G10 G1–G9 G10 G1–G9 G10 G1–G9 G10 G1–G9 G10 G1–G9 G10

Coherent trajectories for K = 3 G1–G10,G27–G30 G11–G26, G36 G1–G10,G27–G30 G11–G26 G1–G10 G11–G26 G1–G10 G11–G26 G1–G10 G11–G26 G1–G10 G11–G26

Area 5 Area 4 Area 3

15.38 30.76 33.25 36.58 39.91 41.48

30.76 33.25 36.58 39.91 41.48 49.23

Area 2 Area 1 Time (s)

FLOPS

Min.

Avg.

1.94 3.19

2.25 3.62

3072 38,604

0.18 0.45

0.21 0.45

355,650 2,835,672

0.27 1.00

0.30 1.20

1,422,402 11,342,040

Table 5 Coherent trajectories by slow-coherency for scenarios 1–2. Scenario 1

Table 2 Time evolution of Coherent trajectories in Scenario 2 using TF-filters for K = 1 and K = 3.

181

Table 3 Relationship between equivalent inertia and apparent power per Area for Scenario 2.

Area 9

Area 11

Area 12

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

Scenario 2

Area

Slow-coherency for 5 slowest eigenvalues

Area

Slow-coherency for 9 slowest eigenvalues

1 2 3 4 5 – – – –

G1–G9 G10–G13 G14 G15 G16 – – – –

1 2 3 4 5 6 7 8 9

G1–G9 G11,G12,G36 G13-G26,G31,G33 G10,G27–G30 G34,G35 G32,G37,G38,G40,G41 G39,G42 G43–G47 G48

according to the similarity of disturbance response trajectory. However, a wide range of research [53], [54] shows that at the start of a disturbance, e.g. a small number of cycles, the power angle trajectory are normally clustered into two groups. Especially for realistic power grids, due to the presence of large inertia as well as synchronization effect, it is difficult to observe multi-clustering separation phenomenon right after the fault clearance [55]. Usually, three or more clusters of machines are gradually drawn from the generator group with the highest acceleration. But during the specific interval shortly after fault clearing, it is accepted that all the generators are considered to be contained in two clusters. For comparison, the coherency grouping from the hierarchical agglomerative clustering algorithm and the Slow-Coherency Grouping Algorithm presented in [8] are discussed. The latter is computed based on linearized models, attaining a set of generators for each scenario, which is shown in Table 5. Note that the areas for both algorithms completely match in Scenario 1 (see Figs. 5, 12 and 13). In Scenario 2, the areas for both algorithms are clustered according to Fig. 6, where MB refers to modal-trajectoryBased, and SC stands for the Slow-Coherency grouping algorithm. An existing system equivalent tool based on coherency methods is the DYNRED program which supports coherency identification by

182

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183

two methods [11]: weak links and Tolerance-Based Method. The former simplifies all synchronous machines as classical machines. Then, a coherency matrix is formed based on the angular response to current injections at each classical machine. The latter requires the eigenvector output file from OAM, part of the DSAToolsTM SSAT package. Thus, the main differences between the proposed approach and existing coherency methods are that the DYNRED uses the network topology’s information and works for one operating condition; moreover, a threshold and a small-signal analysis are required prior to applying the TBM. The proposed methodology does not need the number of desired groups or defining a threshold value. 6. Conclusions This paper has demonstrated the use of Taylor-Fourier filters for extracting modal components from disturbance trajectories, using the hierarchical agglomerative strategy together with the Elbow’s method, for clustering coherent groups. The properties of the frequency responses of the Taylor-Fourier filters with maximally flat pass-bands and zero flat stop-bands at the frequencies of concern fi , enable the spectral separation of the oscillatory modes. Both a linear and nonlinear filtering process have been illustrated. By extracting the inter-area modes from the disturbance trajectories, the proposed method is able to find slow coherent generators, which can then be aggregated to obtain reduced-order models. Furthermore, this method is able to find along time coherent groups, in contrast to almost all existing coherency grouping methods which yield fixed coherent groups. References [1] O. Gomez, M.A. Rios, Real time identification of coherent groups for controlled islanding based on graph theory, IET Gen. Trans. Distrib. 9 (8) (2015) 748–758. [2] H. You, V. Vittal, X. Wang, Slow coherency-based islanding, IEEE Trans. Power Syst. 19 (1) (2004) 483–491. [3] J. Li, C.C. Liu, K.P. Schneider, Controlled partitioning of a power network considering real and reactive power balance, IEEE Trans. Smart Grid 1 (3) (2010) 261–269. [4] B.P. Padhy, S.C. Srivastava, N.K. Verma, A coherency-based approach for signal selection for wide area stabilizing control in power systems, IEEE Syst. J. 7 (4) (2013) 807–816. [5] N. Kishor, L. Haarla, J. Turunen, M. Larsson, S.R. Mohanty, Controller design with model identification approach in wide area power system, IET Gen. Trans. Distrib. 8 (8) (2014) 1430–1443. [6] U.D. Annakkage, N.K.C. Nair, Y. Liang, A.M. Gole, V. Dinavahi, B. Gustavsen, T. Noda, H. Ghasemi, A. Monti, M. Matar, R. Iravani, J.A. Martinez, Dynamic system equivalents: a survey of available techniques, IEEE Trans. Power Deliv. 27 (1) (2012) 411–420. [7] R. Singh, M. Elizondo, S. Lu, A review of dynamic generator reduction methods for transient stability studies, 2011 IEEE Power and Energy Society General Meeting (2011) 1–8. [8] J.H. Chow, Power System Coherency and Model Reduction, Springer, NY, USA, 2013. [9] J. Zaborszky, K.-W. Whang, G. Huang, L.-J. Chiang, S.-Y. Lin, A clustered dynamic model for a class of linear autonomous systems using simple enumerative sorting, IEEE Trans. Circ. Syst. 29 (11) (1982) 747–758. [10] R. Nath, S.S. Lamba, K.S.P. Rao, Coherency based system decomposition into study and external areas using weak coupling, IEEE Power Eng. Rev. PER 5 (6) (1985) 52–53. [11] W. Price, G. Boukarim, J. Chow, R. Galarza, A. Hargrave, B. Hurysz, R. Tapia, Improved dynamic equivalencing software, final report, EPRI Project RP 2447-02. [12] Y. Susuki, I. Mezic, Nonlinear koopman modes and coherency identification of coupled swing dynamics, IEEE Trans. Power Syst. 26 (4) (2011) 1894–1904. [13] H.R.C. Mehrdad Ghandhari, R. Eriksson, Coherent groups identification under high penetration of non-synchronous generation, IEEE PES General Meeting (2016) 1–5. [14] D. Molina, G.K. Venayagamoorthy, J. Liang, R.G. Harley, Intelligent local area signals based damping of power system oscillations using virtual generators and approximate dynamic programming, IEEE Trans. Smart Grid 4 (1) (2013) 498–508. [15] M. Aghamohammadi, S. Tabandeh, A new approach for online coherency identification in power systems based on correlation characteristics of generators rotor oscillations, Int. J. Elect. Power Energy Syst. 83 (2016) 470–484.

[16] K. Tang, G.K. Venayagamoorthy, Online coherency analysis of synchronous generators in a power system, Innov. Smart Grid Tech. Conf., IEEE PES (2014) 1–5. [17] M.A.M. Ariff, B.C. Pal, Coherency identification in interconnected power system – an independent component analysis approach, IEEE Power Energy Society General Meeting (2013), pp. 1–1. [18] J. Wang, K. Guo, A on-line coherence identifying method based on independent component analysis, Image and Signal Processing, 7th Int. Cong. (2014) 1090–1094. [19] M.A.M. Ariff, B.C. Pal, Coherency identification in interconnected power system-an independent component analysis approach, IEEE Trans. Power Syst. 28 (2) (2013) 1747–1755. [20] R. Agrawal, D. Thukaram, Support vector clustering-based direct coherency identification of generators in a multi-machine power system, IET Gen. Trans. Distrib. 7 (12) (2013) 1357–1366. [21] R. Agrawal, D. Thukaram, Identification of coherent synchronous generators in a multi-machine power system using support vector clustering, Power and Energy Systems (ICPS), Int. Conf. (2011) 1–6. [22] A.M. Khalil, R. Iravani, A dynamic coherency identification method based on frequency deviation signals, IEEE Trans. Power Syst. 31 (3) (2016) 1779–1787. [23] J. de la, O. Serna, E. Vazquez Martinez, Smart grids part 2: synchrophasor measurement challenges, IEEE Instrum. Meas. Mag. 18 (1) (2015) 13–16. [24] A. Zamora, V.M. Venkatasubramanian, J.A. de la, O. Serna, J.M. Ramirez, M. Paternina, Multi-dimensional ringdown modal analysis by filtering, Electr. Power Syst. Res. 143 (2017) 748–759. [25] J. de la, O. Serna, Dynamic phasor estimates for power system oscillations, IEEE Trans. Instrum. Meas. 56 (5) (2007) 1648–1657. [26] R. Podmore, Power System Coherency and Model Reduction, Springer-Verlag, New York, 2013, pp. 15–38, Ch. Coherency in Power Systems. [27] V. Venkatasubramanian, J.R. Carroll, Oscillation Monitoring System. http:// www.naspi.org/. [28] D.J. Trudnowski, J.W. Pierre, Signal processing methods for estimating small-signal dynamic properties from measured responses, in: Analysis of Nonlinear and Non-stationary Inter-area Oscillations: A Time Frequency Perspective, Springer, 2009. [29] Z. Tashman, H. Khalilinia, V. Venkatasubramanian, Multi-dimensional Fourier ringdown analysis for power systems using synchrophasors, IEEE Trans. Power Syst. 29 (2) (2014) 731–741. [30] M.R.A. Paternina, J.M. Ramirez, A. Zamora-Mendez, Real-time implementation of the digital Taylor-Fourier transform for identifying low frequency oscillations, Electr. Power Syst. Res. 140 (2016) 846–853. [31] J. Sanchez-Gasca, D. Trudnowski, Identification of electromechanical modes in power system, IEEE Task Force on Identification of Electromechanical Modes of the Power System Stability, Power & Energy Society (2012, June), Tech. rep. [32] I. Kamwa, A.K. Pradhan, G. Joos, Robust detection and analysis of power system oscillations using the teager-kaiser energy operator, IEEE Trans. Power Syst. 26 (1) (2011) 323–333. [33] J. de la, O. Serna, J. Ramirez, A. Zamora Mendez, M. Paternina, Identification of electromechanical modes based on the digital Taylor-Fourier transform, IEEE Trans. Power Syst. 31 (1) (2016) 206–215. [34] P.-N. Tan, M. Steinbach, V. Kumar, Introduction to Data Mining, Addison-Wesley, 2013. [35] C.D. Manning, P. Raghavan, H. Schütze, Introduction to Information Retrieval, Cambridge University Press, New York, NY, USA, 2008. [36] T.S. Madhulatha, An Overview on Clustering Methods. http://arxiv.org/abs/ 1205.1117. [37] C. Canizares, T. Fernandes, E.G. Junior, G.L. Luc, M. Gibbard, I. Hiskens, J. Kersulis, R. Kuiava, L. Lima, F.D. Marco, N. Martins, B. Pal, A. Piardi, R. Ramos, J. dos Santos, D. Silva, A.K. Singh, B. Tamimi, D. Vowles, Benchmark models for the analysis and control of small-signal oscillatory dynamics in power systems, IEEE Trans. Power Syst. PP (99) (2016), 1–1. [38] J. de la, O. Serna, Taylor-Fourier analysis of blood pressure oscillometric waveforms, IEEE Trans. Instrum. Meas. 62 (9) (2013) 2511–2518. [39] J.A. de la, O. Serna, Dynamic phasor estimates for power system oscillations, IEEE Trans. Instrum. Meas. 56 (5) (2007) 1648–1657. [40] M.A. Platas-Garza, J.A. de la, O. Serna, Dynamic harmonic analysis through Taylor-Fourier transform, IEEE Trans. Instrum. Meas. 60 (3) (2011) 804–813. [41] J.A. de la, O. Serna, M.A. Platas-Garza, Maximally flat differentiators through WLS Taylor decomposition, Digit. Signal Process. 21 (2) (2011) 183–194. [42] D. Xu, Y. Tian, A comprehensive survey of clustering algorithms, Ann. Data Sci. 2 (2) (2015) 165–193. [43] P.J. Slater, Centers to centroids in graphs, J. Graph Theory 2 (3) (1978) 209–222. [44] R. Tibshirani, G. Walther, T. Hastie, Estimating the number of clusters in a dataset via the Gap statistic 63 (2000) 411–423. [45] M.R.A. Paternina, J.M. Ramirez-Arredondo, J.D. Lara-Jiménez, A. Zamora-Mendez, Dynamic equivalents by modal decomposition of tie-line active power flows, IEEE Trans. Power Syst. 32 (2) (2017) 1304–1314. [46] T.H. Cormen, C.E. Leiserson, R.L. Rivest, C. Stein, Introduction to Algorithms, edición: 3 Edition, MIT Press, Cambridge, MA, 2009. [47] A.V. Dolgushev, A.V. Kelmanov, On the algorithmic complexity of a problem in cluster analysis, J. Appl. Ind. Math. 5 (2) (2011) 191–194. [48] J.H. Chow, K.W. Cheung, A toolbox for power system dynamics and control engineering education and research, IEEE Trans. Power Syst. 7 (4) (1992) 1559–1564, http://dx.doi.org/10.1109/59.207380.

M.R. Arrieta Paternina et al. / Electric Power Systems Research 158 (2018) 170–183 [49] Texas Instruments Inc., Dallas, TX, TI. (2012, August). TMS320F28335 Digital Signal Controllers. Data Manual. http://www.ti.com/lit/ds/sprs439m/ sprs439m.pdf. [50] TI (1991, August). TMS320 Second-Generation Digital Signal Processors. Data Manual, Texas Instruments Inc., Dallas, TX. http://pdf.datasheetcatalog.net/ datasheets/700/315970 DS.pdf. [51] International Business Machines Corp., Hopewell Junction, NY, IBM (2008, January). PowerPC 970MP RISC Microprocessor. Version 1.3. [52] A. Pai, Energy Function Analysis for Power System Stability, Springer Science & Business Media, 2012.

183

[53] M. Pavella, D. Ernst, D. Ruiz-Vega, Transient Stability of Power Systems: A Unified Approach to Assessment and Control, Springer Science & Business Media, 2012. [54] Y. Xue, M. Pavella, Extended equal-area criterion: an analytical ultra-fast method for transient stability assessment and preventive control of power systems, Int. J. Electr. Power Energy Syst. 11 (2) (1989) 131–149. [55] L. Youbo, L. Yang, L. Junyong, L. Maozhen, M. Zhibo, G. Taylor, High-performance predictor for critical unstable generators based on scalable parallelized neural networks, J. Mod. Power Syst. Clean Energy 4 (3) (2016) 414–426.