5.3.1 N-1 Contingency Monte Carlo Simulation

The DVRs are proposed in order to take advantage of the great potential of PMUs to allow performing post-contingency DVA that can be used to trigger corrective control actions in real time. Conventionally, these corrective actions are set to occur when specific pre-established operational conditions are reached, and they are unable to work under unconsidered contingencies that could begin cascading events.

Monte Carlo simulation allows obtaining the dynamic PMU post-contingency variables, which are used to determine the DVRs for the bulk power system under analysis. Figure 5.2 depicts the Monte Carlo simulation procedure in detail.

In order to adequately apply Monte Carlo simulation and obtain realistic results of the post-contingency dynamic response, some considerations and modeling requirements are taken into account:

• The DVRs can be defined for either short-term or long-term phenomena. This chapter focus on the three phenomena defined as short-term stability, which comprises transient stability, short-term voltage stability, and short-term frequency stability (TVFS) [9].
• Since the accuracy of the DVR boundaries depends on the accuracy of the models, the dynamic components (generators, motors, loads) and relevant control systems (such as excitation control system and speed governor systems) must be modeled with enough detail [8].
• The simulated events are based on N-1 contingencies, and the vulnerability criterion consists in the possibility of this kind of disturbance driving the system to further undesirable events (i.e., N-2 contingencies), which are considered as the beginning of a cascading event. For this purpose, three types of local protection devices are modeled in this chapter: out-of-step relays (OSR), under and over voltage relays (VR), and under and over frequency relays (FR). The tripping of one or more of these protection relays is used as the indicator of an imminent cascading event, and so it becomes a sensor of system vulnerability status. Thus, after each MC iteration, a vulnerability status flag gets the value of “1” (i.e., status of vulnerable) if one or more of the above mentioned local protections were tripped during the time domain simulation. Otherwise, in the case that none of the local protections was tripped, this vulnerability status flag acquires the value of “0” (i.e., status of non-vulnerable). This vulnerability status indicator is stored together with the post-contingency dynamic data. At the end of the MC simulation, a vulnerability status vector (vs) will contain all the vulnerability status flags corresponding to each iteration.

5.3.2 Post-contingency Pattern Recognition Method

The required mathematical tool for achieving the tasks involved in the post-contingency DVA has to be capable of predicting the post-contingency system security status and of specifying the actual dynamic state relative position with respect to the DVR boundaries, with reduced computational effort and quick-time response. Pattern recognition-based methods have the potential to effectively achieve this goal.

Pattern recognition is concerned with the automatic discovery of similar characteristics in data by applying computer algorithms and using them to take some action such as classifying the data into different categories [10], called classes [11].

Dynamic electrical signals can exhibit certain regularities (patterns) signalling a possibly vulnerable condition. However, these patterns are not necessarily directly evident in the electrical signal, though it may contain hidden information which can be uncovered by a proper pattern recognition tool [3].

This chapter applies a novel methodology for recognizing patterns in post-contingency bus features in order to numerically determine the DVRs of a specific Electric Power System. The method uses a time series data mining technique based on EOFs. The data consist of measurements of some dynamic post-disturbance PMU variables (voltage phasors or frequencies) at m spatial locations (buses where PMUs are installed) at r different times of n Monte Carlo cases. The measurements are arranged so that p time measures are treated as variables () and every Monte Carlo simulation plays the role of observation, constituting an () time series data matrix (F).

In the case of voltage angles, and since transient stability depends on angle separation, it is recommended not to use the voltage angles directly but to employ normalized angles obtained via the previous application of (5.2):

5.2

where x_i(t) represents the voltage angle on bus i at time t; x_i0 is its initial value prior to the perturbation; and is the average value over the number of buses N_b at time t:

5.3

Once matrix F has been structured, the method for obtaining EOFs is applied to this matrix. EOFs are the result of applying singular value decomposition (SVD) to time series data [12, 13]. Some authors consider that Principal Component Analysis (PCA) and EOF are synonymous [13]; nevertheless, this chapter uses a different interpretation of these two techniques. Whereas PCA is a data mining method that allows reducing the dimensionality of the data, EOF is a time series data mining technique that allows decomposing a discrete function of time f(t) (such as voltage angle, voltage magnitude, or frequency) into a sum of a set of discrete pattern functions, namely EOFs. Thus, EOF transformation is used in order to extract the most predominant individual components of a compound signal waveform (similar to Fourier analysis), which allows revealing the main patterns buried in the signal.

The main approaches related to EOFs have been developed for use in the analysis of spatio-temporal atmospheric science data, whereas their application in other scientific fields continues to be scarce. The data concerned consist of measurements of specific variables (such as sea level pressure, temperature, etc.) at n spatial locations at p different times [14]. As mentioned before, this chapter employs a variation of this definition, where the n spatial locations are replaced by n different post-contingency power system states (obtained from MC simulation), and the p different times consist of PMU instant values of post-disturbance dynamic variables (voltage phasors or frequencies) measured on m buses, at r different instants that depend on the selected time window (i.e., ).

Therefore, an () data matrix of discrete functions (F) is structured, where the post-contingency measurements at different power system states (n) are treated as observations, and the PMU samples belonging to a pre-specified time window (p time points) play the role of variables. Since the different power system states result from the application of MC-based simulations, n is conceptually greater than , and so F is a rectangular matrix:

5.4

where f_k is the k-th discrete function of time obtained in the k-th MC-repetition that consists of p samples.

Formally, the SVD of the real rectangular matrix F of dimensions () is a factorization of the form [15]

5.5

where U is an orthogonal matrix whose columns are the orthonormal eigenvectors of FF′, V′ is the transpose of an orthogonal matrix whose columns are the orthonormal eigenvectors of F′F, and Λ^1/2 is a diagonal matrix containing the square roots of eigenvalues from U or V in descending order, which are known as the singular values of F.

Taking into account that , this matrix decomposition can be written, as a finite summation, as follows:

5.6

where u_i and v_i are the i-th column eigenvectors belonging to U and V respectively, and λ_i^1/2 is the i-th singular value of F.

From (5.6), and after some computations, each element of F (each discrete function) can be written as follows:

5.7

It is worth mentioning that the expression shown by (5.7) actually represents the decomposition of the discrete function of time f_k into a sum of a set of discrete functions (v_j) which are orthogonal in nature (since they are the orthonormal eigenvectors of F′F), weighted by real coefficients resulting from the product of the j-th singular value of F by the j-th element of the eigenvector u_k. Thus, v_j represents the j-th EOF and its coefficient is called the EOF score.

Based on a generalization of (5.7), it is possible to reconstruct the complete matrix F (i.e., the original data) using the EOFs and their corresponding EOF scores as given in (5.8):

5.8

where a_i is the i-th vector whose elements are all the a_ij EOF scores.

Comparing (5.6) and (5.8), it is easily concluded that . Then, all a_ij EOF scores can be calculated by using its matrix form as follows:

5.9

where A is the EOF score matrix.

From (5.5) and (5.9), it can be determined that

5.10

From the last equation, and based on the fact that V is an orthogonal matrix (whose main feature is that its transpose is equal to its inverse), the EOF score matrix A can be computed by (5.11):

5.11

where matrix V contains the corresponding EOFs of F (i.e., the eigenvectors of F′F).

The sum of the singular values of F (λ_i^1/2) is equivalent to the total variance of the data matrix, and each i-th singular value offers a measurement of the explained variability (EV_i) given by EOF_i as defined by (5.12). Thus, the number of the chosen EOFs depends on the desired explained variability:

5.12

It is worth mentioning that the main advantage of EOFs is their ability to determine the orthogonal functions that better adapt to the set of dynamic functions. This feature enables the mining of the patterns in the signal, and allows EOF to beat other signal processing tools such as Fourier analysis, which always employs the same pre-defined trigonometric pattern functions (i.e., sine and cosine), which are not always suited to representing specific dynamic behavior [16].

Once the EOFs (coefficients and scores) are computed, the corresponding EOF scores represent the projection of the data in the direction of the EOFs and form vectors of real numbers that represent the post-contingency dynamic behavior of the system. These pattern vectors permit mapping of the DVRs represented in the coordinate system formed by the EOFs.

Additionally, since some features in greater numeric ranges can dominate those of smaller numeric ranges, a numerical normalization is recommended before mapping the DVRs. Then, the pattern vectors have to be scaled in order to avoid this drawback [17]. In this chapter, a linear scaling in the range of [0, 1] is proposed.

Each pattern vector has a specific associated “class label” depending on the resulting vulnerability status vector (vs). These class labels might correspond to a non-vulnerable case (label 0) or a vulnerable case (label 1). Using the resulting pattern vectors and their associated vulnerability status class labels, the DVRs can be numerically mapped in the coordinate system formed by the main EOFs. Figure 5.3 presents the scheme of the proposed post-contingency pattern recognition method.

5.3.3 Definition of Data-Time Windows

In order to adequately show the system response for the different stability phenomena, several time windows (TWs) have to be defined. These time windows are established depending on the Monte Carlo statistics of the relay tripping times, influenced also by the WAMPAC communication time delay (t_delay), which represent the delays caused by measurement, communication, processing, and tripping. It is noteworthy that TW corresponds to the length of each time window, whose beginning (t₀) always accords with the instant of the fault clearing (t_cl).

Since three different types of relays are considered, three different tripping time statistical analyses will be carried out. First, the minimum MC tripping time (t_min) has to be determined. This time represents the maximum admissible delay for the actuation of any corrective control action, which has to be also affected by the WAMPAC communication time delay:

5.13

where t_OSRi, t_VRi, and t_FRi are the tripping times of OSRs, VRs, and FRs of the n MC repetitions, respectively. Typically, OSRs present the fastest tripping time due to the fast time-frame evolution of TS.

Since the post-contingency data comprise the samples taken immediately after the fault is cleared, the first time window (TW₁) is defined by the difference between t_min and the clearing time (t_cl):

5.14

The rest of the time windows are defined based on the statistical concept of confidence interval related to Chebyshev's inequality [18], which requires that at least 89% of the data lie within three standard deviations (3σ):

5.15

where std{·} represents the standard deviation (σ) of the relay tripping time that most intersects the corresponding time window TW_k.

5.3.4 Identification of Post-contingency DVRs—Case Study

For illustrative purposes, this section shows the recognition of the post-contingency DVRs of the IEEE New England 39-bus, 60 Hz [19] and 345 kV test system, slightly modified in order to satisfy the N-1 security criterion and to include dynamic load models, automatic voltage regulators (AVR), governors (GOV), and the three mentioned types of relays: OSR, VR, and FR. Additionally, six PMUs have already been located at buses 2, 9, 11, 17, 22, and 28, based on the results presented in [20]. The test system single-line diagram is depicted in Figure 5.4, which also includes the PMU location.

The placement PMUs are considered to have a one-cycle (16.67 ms) updating period. Also, it is assumed that the WAMPAC communication time delay is (i.e., a conservative value derived from the reference data presented in Chapter 1).

Simulations were carried out in DIgSILENT PowerFactory and consisted of several contingencies where the causes of vulnerability could be transient instability, short-term voltage instability, or short-term frequency instability (TVFS). These contingencies are generated by applying the Monte Carlo method.

Two types of events are considered: three phase short circuits and generator outage. The short circuits are applied at different locations of the transmission lines, depending on the Monte Carlo simulation. The disturbances are applied at 0.12 s, followed by the opening of the corresponding transmission line at 0.2 s (i.e., 80 ms fault clearing time t_cl). Likewise, the generator to be tripped is also chosen by the Monte Carlo method, and this type of contingency is applied at 0.2 s. Several operating states have been considered by varying the loads of PQ buses, depending on three different daily load curves. Then, optimal power flow (OPF) is performed in order to establish each operating state, using the MATPOWER package. After that, time domain simulations have been performed using the DIgSILENT PowerFactory software, so that the post-contingency PMU dynamic data could be obtained.

Both voltage components (magnitude and angle), and bus frequencies are considered as potential input variables. A total number of 10 000 cases have been simulated, from which 7600 are stable or “non-vulnerable” and 2400 are unstable or “vulnerable”: 1308 are transient unstable, 682 are frequency unstable, and 410 are voltage unstable. It is worth to mention that 84 cases belonging to the 7600 TVFS non-vulnerable instances correspond to oscillatory unstable cases. Nevertheless, since oscillatory phenomena are treated in a different way in this book, these affairs will be separately analyzed in Chapter 6.

For illustrative purposes, the post-contingency dynamic responses of two different MC iterations are presented in Figures 5.5 and 5.6. These figures depict the simulation results of one non-vulnerable case and one voltage unstable case (vulnerable case). It is possible to observe in the figures the particular signal dynamic behavior for each case, which expose the existence of certain patterns revealing the actual system security level as well as the future system vulnerability status tendency.

By using the MC simulation results, consisting of the post-contingency PMU dynamic data, the relay tripping times, and the vulnerability status vector, the DVRs are identified for this specific power system. First, several time windows are defined based on the Monte Carlo statistics of the relay tripping times. Figure 5.7 presents the 39-bus-system relay tripping time histograms, where the reference frame corresponds to the time simulation.

The out-of-step relay time has a mean of 1.2252 s, a standard deviation of 0.3746 s, and a minimum value of 0.8342 s. Thus, vulnerability assessment has to be done in less than 0.5842 s (). For this reason, an adequate data window (TW₁) for TS phenomenon would be 300 ms () starting from fault clearing.

In this test system, the tripping of voltage relays presents a mean of 4.1275 s, a standard deviation of 1.6872 s, and a minimum value of 3.22 s, whereas the frequency relay tripping time has a mean of 10.6829 s, a standard deviation of 2.4921 s, and a minimum value of 5.987 s. Using these values, and (5.15), the rest of the time windows are determined. This time window definition is summarized in Table 5.1.

TimeWindow	std{tOSR/VR/FR} (s)	(s)	TW(s)
TW1	–	–	0.30
TW2	0.3746	1.4238	1.50
TW3	0.3746	2.6238	2.70
TW4	0.3746	3.8238	3.90
TW5	1.6872	8.9616	9.00

After TW definition, the corresponding EOF coefficients and scores (i.e., pattern vectors) are determined using the resulting MC dynamic data. After this, the EOF scores are computed. For instance, the first four TW₃ voltage-magnitude-based EOFs are presented in Figure 5.8. Similarly, EOFs are determined for each electric variable at every TW.

Once the EOF scores are computed, the corresponding DVRs can be mapped. For instance, Figure 5.9 presents the three dimensional distribution of the pattern vectors obtained from the voltage magnitudes corresponding to TW₃. In the figure, the ellipsoidal area (enclosing the “vulnerable” pattern vectors represented by circles) represents the vulnerable region; whereas the white area (relating to the “non-vulnerable” diamond pattern vectors) corresponds to the non-vulnerable region.

These areas have been empirically delimited, bordering the obtained pattern shapes which depend on the pattern vector spatial locations. The plot presents the DVRs corresponding to the first three EOFs depicted in Figure 5.8.

5.4.1 Support Vector Classifier (SVC) Training

In dynamic vulnerability status classification, the target values resemble the MC vulnerability status “class label.” These class labels might correspond to non-vulnerable cases (label 0) or vulnerable cases (label 1), depending on the triggering status of the different relays. Hence, a two-class classifier has to be used. The attributes, instead, are represented by the post-contingency pattern vectors that better represent the DVRs belonging to each time window. In this connection, and based on the fact that some signals better show the evolution of specific phenomena than others, a procedure for improving the feature extraction and selection has to be followed in order to choose the appropriate pattern vectors before training the SVC. Then, it is necessary to choose an adequate number of EOFs that allows maintaining as much of the variation presented in the original variables as possible. This analysis constitutes the proposed feature extraction stage. For this purpose, the i-th EOF desired explained variability (EV_i), represented by (5.12), is used as a weight measure. After several tests, it has been established that the number of EOFs chosen will be that which permits obtaining an EVi of more than 97%. In the case where the EVi requirement is satisfied already with the first EOF, two EOFs will be used in order to ensure the mapping of at least two-dimensional DVRs. On the other hand, based on the premise that some electrical variables better show the evolution of specific phenomena than others, it is necessary to select those signals that allow the best DVR representation for each TW. The selection will permit obtaining the best classification accuracy. This analysis can be seen as the feature selection stage. Therefore, a decision tree (DT) method for feature selection, originally introduced in [23], has been adapted to this problem, based on the classification accuracy obtained considering a combinatorial analysis of the potential electrical variables. Additional details of this feature extraction and feature selection method can be found in [22].

Once the optimal SCV parameters have been identified and the best features have been determined, the SVC is trained for each data-time window. The objective is to obtain a robust enough SVC that allows predicting the vulnerability status by using as input the best pattern vectors obtained from the EOFs. For this purpose, 10-fold cross validation (CV) is used to analyze the robustness of the classifier training as explained in [22]. For illustrative purposes, the performance of the trained SVC for each of the five TW determined in the previously presented case study is presented in Table 5.2. This table also includes the performance of other classifiers, such as: decision tree classifier (DTC), pattern recognition network (PRN: a type of feed-forward network), discriminant analysis (DA), and probabilistic neural networks (PNN: a type of radial basis network). The performance of the classification is evaluated, in this comparison, by means of the mean of each CV-fold classification accuracy (CA_i) for each TW [22]. The Library for Support Vector Machines (LIBSVM) [21] is used for designing, training, and testing the SVC models.

From the results, it is possible to observe that the SVC outperforms all other classifiers in terms of classification accuracy. In addition, the SVC is the only classifier that permits obtaining more than 99% of CA for all TWs. These results allow us to form the conclusion that the SVC presents an excellent performance when adequate parameters and features are used. Then, suitable parameters permit exploiting the SVC advantage of being more robust to over-fitting problems, which is reflected in a more robust 10-fold CV accuracy. This pre-requisite has been optimally satisfied via the proposed MVMO^S-based parameter identification method.

5.4.2 SVC Real-Time Implementation

For real-time implementation, the off-line trained SVCs will be in charge of classifying the post-contingency dynamic vulnerability status of the power system, using the actual post-contingency PMU voltage phasors and frequencies as data. First, the dynamic signals must be transformed to the corresponding EOF scores. For this purpose, the measured data have to be multiplied by the EOFs determined in the off-line training and stored in the control center processor. Then, the obtained EOF scores will be the input data to the trained SVCs, which will automatically classify the system status into “vulnerable” or “non-vulnerable.” This vulnerability status prediction will be then used by the “performance-indicator-based real-time vulnerability assessment” module (for more details, see Chapter 6) that is responsible for computing several vulnerability indices, which are the final indicators of the system vulnerability condition. Since this classification uses patterns hidden in the input signal as input data, it is capable of predicting the post-contingency dynamic vulnerability status before the system reaches a critical state. Then, the methodology presented in this chapter allows the assessment of the second aspect regarding the vulnerability concept, whose causes concern TVFS issues. Figure 5.12 illustrates how the proposed TVFS vulnerability status prediction method might be implemented in a control center.

In order to validate the complete TVFS vulnerability status prediction presented in this chapter, the power system protection concepts of dependability (Dep (5.19): ability for tripping when there exists the necessity to trip) and security (Sec (5.20): ability not to trip when there is no necessity to trip) are used. It is worth mentioning that security is more critical than dependability because a miss-classification might suggest wrong corrective control actions:

5.19

5.20

Table 5.3 shows a summary of the number of vulnerable and non-vulnerable cases corresponding to the case study of the 39-bus system, where the percentages of complete classification accuracy, security, and dependability are also presented. Results highlight the excellent performance of the methodology regarding security and dependability, offering more than 99% confidence in both measures. The small percentage of errors is due to the existence of overlapping zones between vulnerable and non-vulnerable regions.

Considering the necessity of ensuring very quick response for WAMPAC applications, the methodology presented in this chapter must deliver in a very short elapsed time. Thus, in order to verify the speed of the real-time computations of the proposed method, a complete classification test has been made using only one of the cases analyzed in the previous simulations, which would represent a real-time event. Simulations have been run in Windows 7 Ultimate-Intel (R) Core (TM) i3-2350M-2.30GHz-6GB RAM. In this instance, the resulted SVC elapsed time was 0.68 ms, and the entire process presented in Figure 5.12 takes only 1.47 ms. This very quick response of the proposed methodology matches the very short time delay requirement of this type of real-time application.