Data-Driven Transient Stability Boundary Generation for Online Security Monitoring

Stability boundary can be formulated as a set of constraints shown as follows:

where $\bm{P}$ are the power flow equations, $\bm{J}$ are the static operation constraints, $\bm{Z}$ are the initial condition equations, and $\bm{\Phi}^{k}$ are the transient constraints for the $k$-th contingency. $\overline{\bullet}$ and $\underline{\bullet}$ are the upper and lower bound of the variable. While, as for variables in Eq. (1), $u$ are static variables of the given operation point, including active and reactive power of generators, and amplitude and phase angle of voltages. $x_{0}$ are the initial values of state variables. $\bm{x}^{k}$ and $\bm{y}^{k}$ are the time-variant state variables and operating variables respectively that describe dynamic response of the system after the $k$-th contingency by solving the differential equations $\bm{F}^{k}$ and algebraic equations $\bm{G}^{k}$ shown in the following equations:

As a matter of fact, only the last equation in Eq.(1) is time-consuming while formulating the boundary because of the simulation part in Eq.(2). Therefore, in this paper, we only focus on the transient constraint or boundary of this type.

As for the transient stability inequality constraint, there are several ways to formulate. In this paper, we choose the absolute deviation of rotor angle with respect to the center of inertia (COI). It is shown in Eq.(3),

where $i=1,2,...,N_{G}$, $N_{G}$ stands for the number of generators in the given power system, $T$ is the length of simulation time window while solving the DAEs defined in Eq. (2), $\delta_{i}^{k}(t)$ denotes the absolute rotor angle of the $i$-th generator under the $k$-th contingency at time $t$, and $\delta_{COI}^{k}(t)$ indicates COI which is defined as:

where $T_{Ji}$ denotes inertia time constant of the $i$-th generator.

Since it is difficult to handle Eq. (3) directly, the constraint transformation technique is applied here to define a equivalent transient stability index as follows:

where $\theta^{k}[\delta^{k}(t|u)]$ is defined as:

In Eq.(6), $\varLambda$ indicates a constant, in order to distinguish the stable scenario from unstable one, and it is given by

Therefore, $\bm{\Phi}^{k}(u)$ can be seen as an index to measure the dynamic performance of a given OP.

To guide the generated OPs near the stability boundary, it is of significance to obtain the first-order derivative information of the stability index $\bm{\Phi}$. To do so, forward and adjoint sensitivity analysis are two major approaches. The latter one is preferred due to its lower computational burden. The first-order derivative information can be calculated as Eq. (8).

where $\bm{H}^{k}[u,\bm{x}^{k}(t),\bm{y}^{k}(t),\bm{\lambda}^{k}(t),\bm{\beta}^{k}(t)]$ is the Hamilton function corresponding to the transient constraint for the $k$-th contingency, and can be defined as follows:

where $\bm{\lambda}^{k}$ and $\bm{\beta}^{k}$ are the time solution to co-state equation:

The detailed proof of Eq.(8)-(10) can be referred to [13].

In order to obtain more data samples close to the boundary, three main indices are required: one measures the distance between current sampled OP and boundary, the other two are the direction and step size guiding the current OP moving close to the boundary. In terms of distance measurement, we choose $\bm{\Phi}^{k}(u)$ reflecting how stable/unstable the system is. In terms of direction, the first-order derivative information $\triangledown_{u}\bm{\Phi}^{k}(u)$ is employed. While, for step size, we design a variable step algorithm, in order to move the sampling OPs near the stability boundary as fast as possible. In that case, more critical OPs near the TSB can be sampled and simulated. Considering that larger step size can be set when the current OP is far from the boundary, while the opposite is smaller. Therefore, the step size $\zeta_{m}$ is can be calculated as:

where $u_{max}$ indicates a vector consisting with the maximum value of the controllable variables in the OPs, while $\nu_{i,j}$ is the coefficient with respect to the stability index $\bm{\Phi}^{k}(u)$. Thus, the next sampling OPs can be calculated as follows:

After repeating this step until the new OP is found to be a critical one (close to the boundary), just as the route shown in Fig 1. However, some search processes cross critical area as the route considering highly nonlinear character of transient stability problem, a small trick like binary search is employed to handle such cases. Additionally, the new transient samples can also be positioned perpendicular to the gradient direction with small step size based on the existing samples near the boundary, as the route , to further enrich transient dataset close the boundary. Throughout the data collection process of these three routes, data samples with high information entropy are obtained for data-driven security assessment.

After obtaining sufficient data samples using the strategy in previous subsection, TSB can be generated using data-driven algorithms, like decision trees, support vector machine, k-nearest neighbors algorithm, convolutional neural network, and etc. No matter what exact algorithm we choose, the essence of the boundary can be summarized as a complicated function with respect to the OP $u$ and contingency $k$ given as follows:

However, it is an undeniable fact that there might exist some sampling gap close to the security bound, although the sampling strategy improve the information entropy. Some parts of boundary may be inaccurate due to such problem. Thus, a re-sampling mechanism is proposed to find out new critical scenarios ignored in the previous search process. Essentially, this issue can be re-stated as an optimization problem that maximizes the distant to the existing sampling data closest to the stability boundary. Assume that $N_{N}$ data points $\{(u_{j1},...,u_{j\iota},...u_{j(N_{G}-1)})\}_{j=1,2,...,N_{N}}$ have been generated, and $N_{N(\text{new})}$ new points $\{(\widetilde{u}_{\varrho 1},...,\widetilde{u}_{\varrho\iota},...\widetilde{u}_{\varrho(N_{G}-1)})\}_{\varrho=1,2,...,N_{N(\text{new})}}$ are required in the sampling gap area. Thus, the optimization problem is:

where $j\in[1,2,...,N_{N}]$, $\iota\in[1,2,...,N_{G}-1]$, and $\varrho\in[1,2,...,N_{N(\text{new})}]$.

As we can see from the equation shown above, it is a min-max optimization problem which is hard to solve. To simplify, an auxiliary variable $\gamma$ is introduced into this optimization problem, realizing the chordal decomposition of minimum and maximum problem. Thus, the object function of Eq. (14) can be replaced by Eq.(15) shown as follows:

In fact, it is not necessary to find out the precise global optimal solution, considering that it is only a step for re-sampling. To reduce the required computational time, the duality gap can be preset a larger value. Meanwhile, more initialization points are selected before solving this optimization problem, considering the limitations of nonlinear optimization algorithm.

After focusing on data samples generation strategy, we should evaluate the termination criterion of the generation process. The auxiliary variable $\gamma$ defined in Eq. (15) measures the distance vector between new sampling data and the existing data. Thus, the minimum element of $\gamma$ can be utilized for data generation termination criterion and is defined as follows:

where $\gamma_{\text{cri}}$ is a preset termination threshold. A reasonable value is less than 1% of maximum controllable variables in OPs. Lower threshold improves the accuracy but sacrifices efficiency.

Normally, the proposed method aims at dealing with the changes of possible OPs and contingencies. However, the overall TSB is far too complicated to cover all possible scenarios, and thus we are trying to refresh TBS periodically by tracking current OP in time. In that case, the search space can be reduced significantly, and satisfied the limitation of online computational burden.

As shown in Fig. 2, the search space can be divided into two main dimensional categories: uncontrollable variables and controllable ones. In terms of the former ones, it mainly includes all kinds of loads, uncontrollable generators (like photovoltaic and wind generators). The exact value can be forecast day (hour or minute) ahead with the forecasting error less than 3% using state-of-the-art technique. As a result, variables in this part can be sampled randomly within a given small range. In terms of the latter category, it includes most normal generators. Considering the ramping constraint of each generator, the controllable variables are also be limited in a small range (although it is significantly larger than uncontrollable ones). Noted that the search area has a direct relationship with model periodic refreshing frequency. Overall, this process reduces the search space according to the latest scheduling OP.

After determining the preliminary search space, the data relevance of all possible OPs are required to be analyzed. First of all, the contingency is assumed to be the same. Under such circumstances, we define a matrix $\bm{\Psi}^{k}$ with respect to the impact on different OPs $u$ under the same given contingency $k$ using stability index, that is, $\bm{\Psi}^{k}=\triangledown_{\Upsilon}\bm{\Phi}^{k}(u)$ equals to:

where $\Upsilon$ is the controllable operation variable (e.g., controllable generators), $\epsilon$ and $N_{\Upsilon}$ are respectively the index and the dimension of variable $\Upsilon$. $u$ is the OPs, $\varsigma$ and $N_{u}$ are the index and the dimension of variable $u$. $k$ denotes the index of the given contingency. $\bm{\Psi}^{k}$ is a $N_{u}\times N_{\Upsilon}$ matrix, and the $\varsigma$-th element $\triangledown_{\Upsilon}\bm{\Phi}^{k}(u_{\varsigma})$ has $N_{\Upsilon}$ column, and can be described as:

In order to measure the similarity of each OP, Spearman Correlation (SC) Algorithm is employed to classify all possible OPs into several clusters, based on the response of generator rotor angles for a given OP. To calculate the SC value, the matrix defined in Eq. (17) have to be converted into rank vector. For example, the $r$-th row of $\bm{\Psi}^{k}$ indicates the derivative information of $r$-th OP, and the controllable variable corresponding to the $\rho$-th lowest value is assigned rank $\rho$. Based on this definition, we can calculate the value of SC between the $r$-th and $s$-th OP by the following equation.

where $\overline{r_{u_{r}}}$ and $\overline{r_{u_{s}}}$ are the average value of rank vector with respect to the $r$-th and $s$-th OP, and are defined as follows.

After we got the SC value, spectral clustering technique is employed to classify the contingencies into several clusters. The most severe OP in each cluster can be regarded as the representative of others in the corresponding cluster.

Considering that the power system may undergo all kinds of contingencies, all possible disturbances have to be taken into account. As discussed in the previous subsection, grouping information has been obtained under the same contingency. In this subsection, therefore, we are trying to identify the representative contingencies under various OPs. If the grouping information is similar for two contingencies, only one needs to be analyzed in the next stage. Based on this idea, we employed adjusted rand index (ARI) to represent the relationship between contingencies and then cluster them into several sets. Given an OPs set with $N_{u}$ elements, and two contingency clusters of these elements, namely $p=\{p_{1},...,p_{f},...,p_{N_{p}}\}$ and $q=\{q_{1},...,q_{g},...,q_{N_{q}}\}$. The overlap between these two contingency scenarios $p$ and $q$ are summarized in Eq.(21).

where, each entry $n_{fg}$ (the $f$-th row, $g$-th column) denotes the number of OPs in common between $p_{f}$ and $q_{g}$. Thus, the ARI between contingency $f$ and $g$ can be calculated as follows:

where $C$ denotes the combinatorial operator.

Similar contingency clustering process is carried out as critical OPs selection. So far, several critical contingencies are selected to reduce the computational burden.

Although the aforementioned approaches have declined the number of OPs significantly, the real OPs vary from the ones in the day-ahead scheduling period, and the OPs matching mechanism is required. It is noted that the most critical generators (MCGs), in other words, the controllable generators that have the most significant impact on system transient stability, are almost the same in each OPs cluster. Therefore, MCGs can be determined once the new OP is matched into one certain cluster. We introduce the multivariate Gaussian model to do so. Suppose $N_{u}$ OPs under each contingency, the parameter $u_{\mu}$ and $u_{\Sigma}$ of a given cluster can be calculated:

Then, the possibility of a new OP $u_{\text{new}}$ in cluster $c$ is:

Lastly, the new OP is regarded to be in the cluster with largest possibility calculated in Eq.(24), and MCGs can be chosen accordingly.

In the first case study, the proposed algorithm is applied on a small system: IEEE 9-bus test system. The steady-state and transient parameters are referred in MATPOWER [14] and PSAT manual [15], respectively.

Considering that there are only 3 generators in this grid, it is still possible to visualize the generative process of TSB without dimension reduction and to verify all possible OPs that violate the security constraints of the given system. More specifically, it is highlighted in this subsection that the proposed algorithm is able to generate more data samples close to the TSB, in order to improve efficiency without compromising on accuracy of the boundary. Therefore, we measure not only efficiency but accuracy improvement visually compared with several existing methods.

In this study, all loads are assumed to change randomly and independently within $\pm$10$\%$ of its reference level. The contingency preset is initialed by a three-phase to ground fault at bus 5, and cleared after 0.2 seconds by tripping a line between bus 5 and 7. The stability performance is assessed using TDS and determined by maximum rotor angle difference. For comparison, accurate TSB is generated by brute force, that is, performing TDS at all possible OPs with discretized interval of 1MW.

As shown in Fig. 5(a), 20 initial OPs are randomly selected using Latin Hypercube Sampling method within the given output range of generators. Among them, 11 OPs with the ’$\otimes$’ mark are found in the infeasible area by static security check after solving power flow. The remaining 9 OPs are regarded as the initial seeds to generate the rest samples. Fig. 5(b) shows the process of getting close to the boundary using specially designed critical data sampling strategy introduced in Section II. It can be easily found that the step size varies according to the distance to the boundary. If the sample is getting closer to the boundary, the next step size becomes smaller. This ensures that the sampling near the boundary is sufficient, while traversing the rest less-informative part with fewer samples. Occasionally, some search processes cross the boundary due to the highly nonlinear nature of TSB, binary search is then employed to recover from such cases.

Additionally, the new transient samples can also be searched perpendicular to the gradient direction based on the existing samples near the boundary, as shown in Fig. 5(c). However, it can be observed that there exists the sampling gap on part of TSB. Therefore, more data samples in that area are required to generate a more accurate boundary. A rough boundary is generated as the red line and new sampled points are gotten in the gap area as green triangle markers in Fig. 5(d).

So far, new sampled points can be regarded as the new initial seeds and repeat the above procedures, so as to generate a more accurate boundary with the increase of data samples.

Fig. 6 shows gradient direction information for all possible OPs in static stable area, with discretization interval of 5MW. It is observed that almost all gradient direction arrows point to transient stability boundaries using the proposed transient index and algorithm. Therefore, it proves effective to generate more data samples close to the stability boundary using this specially designed index and method.

Throughout the whole process of the data and boundary generation, the proposed method shows significant effect on efficiency without compromising accuracy. We compare the proposed method with four existing ones, summarized in Table I in details. For comparison more conveniently, we set a minimum accuracy requirement for all methods. Enabled by the proposed method, the data size required to reach 99.9% accuracy declines by nearly 80% and the time required reduces by 32%, compared with the most state-of-the-art method (importance sampling based method). It is highlighted that the proposed method shows better performance.

Meanwhile, the scatter plot of the data samples using different methods is illustrated in Fig. 7. We can easily distinguish the stability boundary through the last two plots (importance sampling based and the proposed methods) rather than the first three. Among them, the data samples generated by the proposed method are much closer to the boundary, and thus shows the superior performance.

In the second case study, the proposed method is applied to a larger and more complex power system named NESTA 162-bus system. Considering the large number of controllable generators and possible ”N-1” contingencies in this grid, computational burden increases geometrically in order to generate TSB. Therefore, we focus on the most critical scenarios selection under different circumstances, to ensure that the online computational burden is under control.

1,000 OPs are selected randomly according to the load prediction and dispatching plan within the given scheduling period. Meanwhile, 512 contingency, which is initiated by a three-phase-to-ground fault at any line close to bus of one end and cleared after 0.2 seconds by tripping the line, is also selected on all these 1,000 OPs. So far, 512,000 samples are selected based on day-ahead scheduling, to find out the most representative scenarios in the next scheduling period.

Various OPs under a single contingency are, firstly, to be clustered into several categories. The matrix $\triangledown_{\Upsilon}\bm{\Phi}\in\Re^{1000\times 12}$ as Eq.(17) is constructed, and employed for clustering. Fig. 8(a) illustrates the eigenvalue of normalized Laplacian matrix in spectral clustering of all 1,000 OPs under contingency #1. It is observed that the first 6 eigenvalues are relatively small, while the others are large. Therefore, cluster number is set to 6. Note that, the number of cluster varies from 2 to 8 for different contingencies. Fig. 8(b) shows the OPs clustering result under contingency #1. The colored matrix shows the correlation relationship between all different OPs under the same contingency. Darker matrix elements indicate weak correlation between the two OPs, and vice versa. Among them, most OPs belong to the first two clusters, and only small amount of OPs fall into other four clusters. Therefore, only 6 critical OPs are taken into account, since they represents all possible OPs within a given scheduling period. In other words, this reduces the number of OPs dramatically from 1,000 to 6.

After obtaining all grouping and critical OPs information for these 512 contingency scenarios, critical contingencies are also required to be identified. Similar algorithm is applied in this task with result shown in Fig. 8(c). As observed, it can be divided into 25 categories in total. So, only 25 contingencies are required to represent 512 preset contingencies.

Additionally, it is also necessary for us to analyze the efficiency of the clustering process, although it is carried out off-line. The time consuming of OPs clustering is only 2.8-3.0 seconds per each contingency. Considering this step is of natural parallel characteristic, asynchronous parallel algorithm can be employed here if necessary. While the time consuming for contingencies clustering is 24.8 seconds. In a word, it is time-effective to find out the most critical OPs and contingencies, compared to doing TDS for a huge amount of scenarios.

In on-line operation stage, OP varies from time to time with different circumstances. Although huge amount of data samples are employed off-line, it is still impossible for almost all real OPs to match with the existing samples exactly. Considering that all clusters are difficult to be distinguished on a 2-D plane, we select three of them, shown in Fig. 9, as an example to clarify this issue. In this figure, three selected OPs clusters are marked with different colors, together with the probability distribution contour plot. As seen in this figure, the dots with magenta asterisk mark belong to cluster #2, and MCGs in these scenarios are #1 and #12. In other words, these two generators are the key for operator to monitor and control in and near the current OP to prevent possible contingencies. Similar results are found for the other clusters except the critical generators. While encountering a new OP in real-world scheduling, multivariate Gaussian distribution probability result is utilized to evaluate the most possible cluster, and to determine MCGs. 1,000 new scenarios are evaluated and 96.1% of them obtain the same index of MCGs. Although the remaining 3.9% scenarios are not the most 2 critical generators, they still rank top 3 or 4. Considering the continuously updating in the scheduling period of the proposed method, it has little impact on the security assessment.

Besides, in order to evaluate the effect of periodic refreshing mechanism using the proposed algorithm in higher dimensional power system, we generate the dynamically updated TSB in a one-hour scheduling period with its load rate ranging from 0.8 p.u. to 1.1 p.u.(see Fig. 10). As shown in Fig. 11(a), all OPs in the search area (in the area surrounded by dotted line) are stable, when the load level is relatively low. At this stage, it only takes less than 10 seconds to generate the boundary within the search area, and make a conclusion that the state of the current OP is safe. It can be observed that, however, with the increase of load level, some OPs in the search area fall into transient unstable area shown in Fig.11(b). As a result, the operator are encouraged to reduce the outputs of generator #1 and #12 (see Fig.11(c)), while increasing others to some extent. By periodic refreshing the boundary, operators can adjust the current OP continuously, maintaining sufficient stability margin all the time. Even if the load level reaches its maximum as Fig.11(d), the OPs still in the stable area. More details can be referred in GIF animation¹¹1Results for more time sections in GIF animation format can be downloaded at http://genggc.org/files/YanTSB2020.gif.

It is noted that the time consuming to generate the TSB is less than 50 seconds (see Fig. 10 for more details) according to the length of boundary in search area. Thus, TSB can be refreshed every minute. Moreover, parallel technique can be employed in the future to generate data samples, because different search path (see route , and in Fig.1) is with the character of naturally parallel. Additionally, it helps in further reducing the refreshing interval to improve the hardware conditions in practical applications.