The graph neural network uses the G ( x , A ) matrix to describe the topological relationship of the power system network, where x represents the collected power grid information feature vector and A represents the adjacency matrix of the topology. The application of the graph neural network in the power system is to aggregate the information of G ( x , A ) having a topological relationship. It can aggregate the characteristics of nodes themselves and neighbors and generate a new feature G ( x ˜ , A ) . G ( x ˜ , A ) , containing original information and topological information, has a higher correlation with the output results. x ˜   represents a new feature vector with higher correlation with the evaluation results, as shown in Figure 1.

The attention mechanism is introduced to weighted summation of the features of adjacent nodes. The weight is completely determined by the features of nodes, which is not affected by the dynamic topology. Each node in the model represents a monitoring node in the grid topology. The input of the attention layer of the graph is the node feature vector set x, as given below: x = ( x 1 → , x 2 → , … , x N → ) T ,     x i → ∈ R F . (1)

Here, x i → is the eigenvector of node i, N is the number of system nodes, and F is the characteristic number of nodes.

The output of each layer is a new node feature vector set x ˜ , as given below: x ˜ = ( x 1 ˜ → , x 2 ˜ → , … , x N ˜ → ) T ,   x i ˜ → ∈ R F . (2)

Here, F ′ is the feature number of the new node, similar to the feature extractor.

For N nodes, input node features predict output new node features. The operation of G ( x ˜ , A ) obtained from G ( x , A ) first needs to calculate the attention coefficient a i j , as given below: a i j = s o f t m a x ( e i j ) = exp ( η ( β ⁡ cos ( W x i → , W x j → ) ) ) ∑ k ∈ N i exp ( η ( β ⁡ cos ( W x i → , W x k → ) ) ) . (3)

Here, the matrix W is initialized, β is the training parameter, cos is the cosine similarity, and η is the LeakyReLU nonlinear activation function. a i j ˜ is determined in Eq. 4 through the adjacency matrix A a i j ˜ = { a i j       A i j = 1 0         A i j = 0 . (4)

Different attention coefficients are generated during each training. Therefore, this study proposes a network structure using a new attention head mechanism. It sets the weight for the important attention head containing topology information, that is, 0 to 1. The model pays attention to the node information of the important attention head and improves the interpretability of the model.

For a GAT layer with k attention heads, each attention head contains a different set of parameters W and a i j . m i represents a vector composed of k attention head weights. The working process is shown in Figure 2.

Among them, node 1 has three attention heads in its neighborhood. Different arrow styles and colors represent independent attention heads. Different soft gates aggregate and control the features of each head to get more relevant feature vectors. The weight formula of attention head is shown below: m i = [ m i ( 1 ) , m i ( 2 ) , … , m i ( k ) ] , = F C θ m α ( x i ⊕ m a x j ∈ N i ( { F C θ g ( x j ) } ) ⊕ ∑ j ∈ N i x j | N i | ) , (5)

Here, the combined maximum pool and average pool are used to construct the network. x j is the adjacent node of node i, θ g means to map the features of neighbor nodes to the vector of g dimensions as small as possible, N i is the number of neighbor nodes of node i, F C θ α = α ( w x + b ) represents the single and full connection layers, α is the sigmoid activation function, ⊕ is the linker, k is the number of attention heads, θ m maps connected features to the k-dimensional space, and m i ( k ) is the weight of the kth attention head of the ith node.

Each node i aggregates the new output features of the attention head topology information, as follows: x i ˜ = ∥ k = 1 K m i ( k ) σ ( ∑ j ϵ N j a i j k ˜ W k x j ) . (6)

Here, after x i ˜ going by the GSTGNN layer, node i contains the output new features of the feature information of different spatial adjacent nodes. The attention head k satisfies inequality k ≥ 2 . σ represents the GELU activation function (Hendrycks and Gimpel., 2016).

In order to prevent overfitting of the model, the random regularization idea is introduced to control the performance on the training set, as follows: GELU s ( x ) = x ϕ ( x ) , = 0.5 x ( 1 + t a n ⁡ h [ 2 π ( x   + 0.044715 x 3 ) ] ) . (7)

Here, φ ( x ) uses N (0,1) normal distribution.

When x decreases, the output value will depend on the input value randomly according to the probability. It improves the generalization ability of the model.

As mentioned before, GSTGNN is introduced to extract power grid topology information to improve the evaluation accuracy. However, the evaluation time is not considered. The system studied in this study is a high-dimensional dynamic system. Not only each node has power grid topology information but also its physical quantity has time series characteristics. The GRU network has the characteristics of forgetting and selective memory. It learns and retains the timing characteristics of input data for later use. In transient assessment, the stability after fault removal can be deduced from the fault state. Therefore, the GRU layer is added to capture the historical sequence and current information of nodes and predict the future time stability. It does not need to input all time data. Thus, it can shorten the evaluation time.

In conclusion, in order to process the sequence information with a complex topological structure and time correlation, GSTGNN and GRU are combined to form a gated spatiotemporal GNN framework. At each time, the input x t and hidden state h t − 1 get new features x ˜ t , h ˜ t − 1 with topology information through the GSTGNN layer, as follows:   x ˜ t = G ( x t , A ) , h ˜ t − 1 = G ( h t − 1 , A ) . (8)

Here, two gating states z t and r t are obtained by the state transmitted from the previous node h t − 1 ˜ and the input of the current node x t ˜ ; z t is the gating of control update, which decides the unit to update its active content to reduce the risk of gradient disappearance; r t is the gating of control reset, which determines the degree of combining the new input information with the previous memory state features. The expression is shown as follows: z t = α ( W z x t [ i , : ] ˜ + b z ) ,   r = α ( W r x t [ i , : ] ˜ + U r h t − 1 [ i , : ] ˜ + b r ) . (9)

Here, α is the sigmoid activation function with a range of (0,1).

Through r t , the new memory will store the information related to the past in the following: h t [ i , : ] ^ = tanh ( W h x t [ i , : ] ˜ + U h ( r t ⊙ U h h t − 1 [ i , : ] ˜ ) + b h ) . (10)

Here, ⊙ means the Hadamard product and multiplies the corresponding elements in the matrix to identify retained and forgotten previous information. h t [ i , : ] ^ remembers the state of previous time by resetting.

The memory state of current time step h t [ i , : ] needs z t to forget and select memory at the same time as follows: h t [ i , : ] = ( 1 − z t ) ⊙ h t − 1 [ i , : ] ˜ + z t ⊙ h t [ i , : ] ^ . (11) Here, ( 1 − z t ) ⊙ h t − 1 [ i , : ] ˜ means selective “forgetting” of the original hidden state and removing unimportant information. z t ⊙ h t [ i , : ] ^ selects some information in h t [ i , : ] ^ , which means forgetting the state information of the past moment and adding some state information input by the current node. In order to improve the generalization ability and evaluation accuracy of the model, topology information is included in the transmission of the abovementioned state information.

Finally, new features aggregated are input into the softmax separator for classification. The output prediction value y ˜ at the last time is obtained in Eq. 12. The value range is (0,1) y ˜ = s o f t m a x ( W o h l a s t [ i , : ] + b o ) . (12)

Here, W o and b o are the parameters learned by gradient back propagation.

In order to start the emergency control immediately after the fault is removed, TSA needs to be performed quickly and accurately. In Sun et al. (2020) and Tian et al. (2020), most of the existing TSA methods use a fixed length observation window; that is, the evaluation time is constant. However, this static evaluation time may not be able to cope with fast transient instability. The system models with different fault degrees need different observation window lengths. After the fault is removed, the dynamic data of the system are observed in the dynamic time window. The stability of the system is evaluated in the future time window, and the observation time window is gradually adjusted. As long as the evaluation system loses stability in the future, the emergency control will be started immediately.

The stability threshold is used to evaluate the evaluation results. Because the GRU layer is used in this study, y ˜ i only focuses on the stability index generated at time T. When i = 1, 2,..., T − 1 is ignored, the rule is shown as follows: S t a b i l i t y = { S t a b l e     1 ≥ y ˜ i > δ U n s t a b l e     y ˜ i < 1 − δ           i = 1,2 , ... , T U n k n o w                 O t h e r s . (17)

Here, δ ∈ ( 0.5 , 1 ) is the stability threshold.

It is necessary to search and set the appropriate threshold δ to balance the accuracy of TSA with the average evaluation time. An adaptive TSA process is proposed in Figure 3.

In this study, x i is input to the model according to the observation time window after fault removal. The evaluation results are predicted moment by moment. When it is within the judgment range, the results are directly output. Otherwise, the hidden state h i at this time and the next time x i + 1 are input to the next neuron for TSA. If the evaluation time exceeds the observation time window, the sliding time window method is adopted until the reliable evaluation result. y i ˜ is obtained, or the maximum allowable evaluation time ( T m a x ) is reached. If T m a x has not been determined, it is regarded as unstable. In this process, h i and x i + 1 are needed to go through the GSTGNN layer to get h i and x i + 1 with topological relationship and inputed to the next evaluation time. It is found that T m a x is generally set to 10 cycles.

With the continuous development of machine learning, for the classification problem, the index based on the confusion matrix in Table 1 can better evaluate the applicability of the classification model than only using the accuracy.

In the power system transient stability analysis, researchers pay more attention to whether the instability is classified correctly. As a result, the following evaluation indexes are mainly based on whether the instability can be correctly judged. In Table 1, TP is the number of stable samples correctly predicted, FP is the number of unstable samples incorrectly predicted, FN is the number of mispredicted stable samples, and TN is the number of unstable samples correctly predicted. In the confusion matrix, the more the number of TP and TN is, the better the number of FP and FN is. However, there are a large number of sample data. It is difficult to measure the model only by the number. Therefore, in order to comprehensively evaluate the TSA performance, four indexes are obtained as follows to calculate the accuracy (ACC), miscalculation rate (recall), precision, and comprehensive evaluation index ( F 1 ) of the model instability.

ACC is the proportion of the number of correct assessments to the total number of assessments, defined as follows: A C C = T P + T N T P + T N + F P + F N × 100 %   . (18)

Once an unstable situation is misjudged as a stable situation, it may lead to a large area power outage of the system. Thus, the unstability recall is used to express the misjudgment rate. The closer the recall is to 1, the lower the possibility of misjudgment of the unstability situation. The formula of recall is shown as follows: r e c a l l = T N T N + F P × 100 %   . (19)

Unstability precision is the proportion of correctly predicted unstable samples, defined as follows: p r e c i s i o n = T N T N + F N . (20)

The higher the precision is, the more accurate it is to predict the unstable situation. However, recall and precision are inversely proportional and cannot reach a high ratio at the same time. Therefore, F 1 , the weighted harmonic average of precision and recall, is introduced in Eq. 21. The two performance indicators are comprehensively considered. F 1 is distributed between [0,1]. The closer it is to 1, the stronger the feature extraction ability of the model and the better the evaluation performance F 1 = 2 × p r e c i s i o n × r e c a l l p r e c i s i o n × r e c a l l . (21)

The accuracy mainly includes ACC and F 1 . Although the deep learning algorithm studied in this study can reach 1 (or 100%) in theory, it will lead to overfitting of the model. The data model does not have enough prediction ability. Therefore, this study introduces the idea of stochastic regularization, which is as close as 100% when the prediction requirements are met.

Besides the accuracy, the average response time (ART) in Eq. 22 is also an important index to evaluate the performance of TSA A R T = ∑ i = 1 T m a x [ T i × C ( T i ) ] / ∑ i = 1 T m a x C ( T i ) . (22)

Here, T m a x is the maximum allowable evaluation time (10 cycles). C ( T i ) is the total number of classified instances in the current evaluation cycle. ART refers to the average TSA time after fault removal.

The total TSA process is shown in Figure 4:

1) Off-line training: In this study, the three-phase permanent short-circuit fault is set under different fault points of different lines. It can realize the diversity of sample data. The time domain simulation method generates data, including voltage amplitude, phase angle, injected active power, injected reactive power, and the maximum power angle difference δ m a x at the end of the simulation. δ m a x obtains the stability index t; when t > 0, it indicates the system transient stability, and the label is “1"; on the contrary, the system is unstable and the label is “0". This study selects the dynamic data of each node on the bus as the input data to determine the system stability in order to be limited to a certain range to reduce the difference of the data. The training data and test data are generally normalized. The weights under instability and stability are obtained from the training data so as to form a weighted cross entropy loss function. The model parameters are continuously adjusted. The Adam method is used to minimize the loss function until the optimal evaluation model is obtained.

2) Online evaluation: Test data and actual data are added, and performance is evaluated by using the best trained model.

This study uses a representative New England 10-machine 39-node power system in Figure 5 to evaluate the performance of the method. The reference frequency is 50 Hz. The time domain simulation data of PSA-BPA software is used, similar to the data generated by PMU in real time. Python is used to program.

The sampling step is set as 0.02 s. Various faults are simulated in various scenarios to obtain complete data sets. A total of 11 load levels of 75–125% (step size 5%) are set. Generator output is changed accordingly to ensure power flow convergence. The Three-phase permanent short-circuit fault is set for the fault type. Fault points are set for 0, 25, 50, 75, and 100% of lines. The fault clearing time is set for 0.2 s. The simulation time is 5 s. Fault samples are selected considering the whole wiring system and N-1 accident. The fault samples are 8855. A total of 6000 training set samples are randomly selected, and test set samples are selected according to 4:1, including 5211 stable samples and 2289 unstable samples.

In order to train the evaluation model with the best performance, these control parameters are defined later. In terms of the model structure, this study sets up two GSTGNN input layers to avoid excessive smoothing caused by too many layers. One GSTGNN middle layer and four attention heads are set to extract and aggregate the features with the topological structure of the power grid. The GRU layer is set to two layers, with 128 batch sizes in each layer, which further improves the evaluation speed. The last layer is the dense layer, which uses the nonlinear activation function to fit the nonlinear problem to improve the evaluation accuracy. In terms of training parameters, the dropout is set to 0.05. The learning rate is set to 0.001. The training iterations are set to 300. The aforementioned settings can improve the performance of the evaluation model.

In the test, each evaluation cycle of the adaptive TSA method based on the sampling step is 0.02 s. The maximum evaluation time is set to 0.2 s ( T max = 10). The observation window T and δ are set to 8 and 0.62, respectively. The specific setting experiment is described later.

In this study, the idea of KNN is introduced to calculate the weight of each sample. The weight is added to the training parameters in the loss function to obtain the best model. A total of 6000 samples with or without imbalance are taken as the training data. The number of iterations is 300. When T = 8, δ = 62; the training process is shown in Figure 6.

After the fault occurs, with the increase of iteration times, the accuracy (ACC) increases rapidly until it reaches 50 times. The accuracy tends to increase gently. Finally, the evaluation accuracy of balanced processing can reach 99.237%. Therefore, data quality has an important influence on the training high-performance model.

For the same set of data, the performance of this method is improved compared with other deep learning models, as shown in Table 2.

According to Table 2, the recall of SVM which belongs to the shallow model is less than 0.9, so it cannot extract features accurately and is prone to misjudgment. GAT, LSTM, and GRU models all have higher ACC. However, when the scale of system topology expands, the F 1 value of the GAT model is higher than that of LSTM and GRU models. It indicates the importance of topological relationship between nodes for TSA. Because the CNN model does not take into account the topological structure and time series information of the power grid, ACC is relatively low. Compared with the previous model, the ACC of the framework established in this study is closer to 100%, and the F 1 value is relatively close to 1. It still has strong generalization performance for high-dimensional data. The framework is relatively more stable. Its performance is better than that of GAT-GRU, GAT, SVM, GRU, LSTM, and CNN models. This indicates that this framework can better mine the essential characteristics of data. Its training parameter sharing overcomes the complex problems of the traditional adaptive evaluation system.

The following table shows the comparison of art and training time in Table 3.

The ART of the GRU model and LSTM model are basically the same (Yu et al., 2018). Because the parameters of the GRU model are less than those of the LSTM model, the training time is relatively short (Chen and Wang. 2021). The extreme learning machine (ELM) network (Zhang et al., 2015) needs to train 10 classifiers whose parameters are not shared. A large number of parameters lead to the longest training time. A simple SVM network structure cannot effectively extract features, resulting in the longest ART. This method extracts the topological relationship and time series information of important attention heads. This method feature extraction ability is higher than that of the GAT-GRU model, which only extracts the topological relationship of attention heads on average. Thus, the output results are easier to reach the stability threshold, and the TSA evaluation speed is faster. In the case of ensuring a high accuracy, the ART of this study is 1.01 cycles, which meets the requirements of emergency control. At the same time, although the model layers of this study and GAT-GRU are more and the offline training time is longer than those of the GRU model and LSTM model, it does not affect the effect of online evaluation. This model is still practical.

In this study, these two parameters T and δ affect the evaluation performance. The following figures in Figure 7 and Figure 8 are the histograms of ACC and ART, respectively.

When T = 8, it means that after the fault is cleared, the input data of the first eight cycles are used for training and continuous time series data are used for testing. The maximum accuracy is about 99.27%, which shows that the input data affect the system performance. With the increase of T, the overall trend of ACC increases (there is a small fluctuation), and the TSA time is longer. The shorter T may damage the integrity of the input data, reducing the accuracy. Therefore, in the case of T as small as possible to ensure a high accuracy, T = 8 is the best choice.

When T = 8, δ for the histogram of ACC and ART is obtained in Figure 9 and Figure 10.

As with the observation time window, the size of δ also affects ACC and ART. The smaller δ may lead to earlier evaluation. However, the accuracy rate may be reduced, while the larger δ produces more accurate results at the expense of evaluation speed. As can be seen from Figures 9, 10, the increase of δ leads to the reduction of ACC and the increase of ART. Small ART and high ACC are the expected results, so the stability threshold of 0.62 is the best.

To sum up, T is set to 8. δ is set to 0.62 to get the fastest evaluation time without reducing the accuracy.

After the fault is removed, the adaptability of the model to topology changes is tested. Four changes are listed in the following table, and 1500 test sets are selected for the performance test in Table 4.

It is found from Table 4 that the accuracy of TSA caused by generator shutdown is the lowest, but it is still more than 99% when the load is cut off. In general, the TSA accuracy of this model is still about 99% in the case of topology changes. The results show that the features selected in this study are enough to reveal the system state of the power grid in the case of large disturbance. The selected model can fully mine the transient variation law and extract data features. The model with strong robustness is suitable for complex topology.

t-Distributed stochastic neighbor embedding (t-SNE) for the extracted feature data is conducive to visually verify the effectiveness of the algorithm. The following is the feature map of different layers of the GSTGNN framework in Figure 11, Figure 12, and Figure 13.

t-SNE is used to establish linear projection. The mapping relationship between approximate high-dimensional data space and low-dimensional embedded space is obtained. The high-dimensional input data are mapped to the two-dimensional space. The stability classification is explained intuitively. Figure 11, Figure 12, and Figure 13 show the characteristic diagrams of the input layer, GRU layer, and full dense layer of the model. With the deepening of the layer, the boundary between the two categories is more and more obvious, and the spatial overlapping data are less and less. It shows that the GSTGNN framework plays an important role in deep structure level extraction. Finally, it achieves an obvious classification effect.