Data-Driven Dynamic Security Partition Assessment of Power Systems Based on Symmetric Electrical Distance Matrix and Chebyshev Distance

Abstract

A rapid dynamic security assessment (DSA) is crucial for online preventive and restoration decision-making. The deep learning-based DSA models have high efficiency and accuracy. However, the complex model structure and high training cost make them hard to update quickly. This paper proposes a dynamic security partition assessment method, aiming to develop accurate and incrementally updated DSA models with simple structures. Firstly, the power grid is self-adaptively partitioned into several local regions based on the mean shift algorithm. The input of the mean shift algorithm is a symmetric electrical distance matrix, and the distance metric is the Chebyshev distance. Secondly, high-level features of operating conditions are extracted based on the stacked denoising autoencoder. The symmetric electrical distance matrix is modified to represent fault locations in local regions. Finally, DSA models are constructed for fault locations in each region based on the radial basis function neural network (RBFNN) and Chebyshev distance. An online incremental updating strategy is designed to enhance the model adaptability. With the simulation software PSS/E 33.4.0, the proposed dynamic security partition assessment method is verified in a simplified provincial system and a large-scale practical system in China. Test results demonstrate that the Chebyshev distance can improve the partition quality of the mean shift algorithm by approximately 50%. The RBFNN-based partition assessment model achieves an accuracy of 98.96%, which is higher than the unified assessment with complex models. The proposed incremental updating strategy achieves an accuracy of over 98% and shortens the updating time to 30 s, which can meet the efficiency of online application.

1. Introduction

Dynamic security assessment (DSA) evaluates the security and stability of the power system under specified operating conditions (OCs) and anticipated faults [1]. It can provide valuable information for preventive and restoration control decision-making [2,3]. With the integration of renewable energy and frequent extreme weather events, the uncertainty of the OCs has increased significantly. Preventive and restoration control decision-making needs to be performed online to adapt to fast-varying OCs. A rapid DSA is required. However, the traditional DSA method based on time-domain simulation is computationally complex [4]. It cannot meet the high assessment efficiency for massive candidate schemes generated during the control decision-making process.

The data-driven DSA method utilizes machine learning to directly learn the map between the input and the dynamic security of power systems [5,6]. It can realize an online assessment due to its high efficiency. The input features are fundamental to data-driven DSA. For a specific fault type, the dynamic security of the power system is closely related to the pre-fault OCs and fault locations. Therefore, the DSA inputs need to incorporate both OCs and the fault location information. To ensure high evaluation efficiency, the DSA typically employs steady-state features as inputs [7,8,9]. Pre-fault OCs can be represented by steady-state power flow features, such as the node injection power, voltage, and line transmission power [10,11]. For fault locations, the traditional DSA method constructs a dedicated assessment model for each fault. It has high assessment accuracy [12,13]. However, an increase in fault locations can lead to a large scale of DSA model sets and high maintenance costs.

Existing DSA methods with a unified assessment model adapted to fault locations can be divided into three categories. The first category leverages machine learning techniques to handle multiple fault locations. Assuming that each OC corresponds to a single anticipated fault, Refs. [14,15] employ transfer learning to build DSA models for multiple faults, of which the accuracy is over 97%. However, they do not include fault location information in the inputs and cannot assess different fault locations under the same OCs. Ref. [16] employs multi-label learning to construct a single DSA model for multiple faults, and improves the DSA accuracy to 99%. The number of considered fault locations is determined by training sample labels. Therefore, the DSA model has no generalization performance in unseen fault locations. The second category employs discrete encoding methods such as binary [17], integer [18,19,20,21], and one-hot encoding [22,23] to represent fault location features. DSA models distinguish different fault locations through these discrete encodings and are generally able to achieve an assessment accuracy of more than 97%. However, discrete encoding fails to reflect the physical differences between fault locations. The discrete encoding-based DSA model has poor generalization performance in unseen fault locations. The third category constructs electrical coordinates to represent fault locations [24,25]. The electrical coordinate is composed of the electrical distance to preselected reference nodes. It can reflect the electrical differences among fault locations. The electrical coordinate-based DSA model has high generalization performance and obtains an assessment accuracy of 98%.

The mapping relationship among OCs, fault locations, and the system dynamic security is strongly nonlinear. To ensure high assessment accuracy, the training set needs to cover a wide variation range of OCs and anticipated fault locations. The size of the training set is thereby large. The learning for massive samples containing a nonlinear map requires a machine learning model with strongly nonlinear fitting capabilities. Typically, increasing the structural complexity of machine learning models can enhance nonlinear fitting capabilities, e.g., long short-term memory (LSTM) neural networks [26], convolutional neural networks (CNNs) [27], and recurrent neural networks (RNNs) [28]. However, these highly complex models face issues such as high training costs, overfitting risks, poor interpretability, and difficulties in model updating. They are hard to update timely to support online preventive and restoration control decision-making.

According to the principle of Occam’s Razor in machine learning, the models with simple structures are preferentially adopted within acceptable training errors. Therefore, this paper proposes a dynamic security partition assessment method of the power system. Instead of a single high-complexity unified assessment model, multiple low-complexity partition assessment models are employed without an accuracy loss. The aim is to achieve an accurate and incrementally updated DSA using low-complexity models. The contributions of this paper are as follows:

(1)

A self-adaptive power grid partition method that strictly accounts for electrical connections is proposed based on the mean shift algorithm. The node locations are represented by a symmetric electrical distance matrix, and the Chebyshev distance is introduced to measure the differences between node location features.

(2)

To build the input feature set for DSA, high-level steady-state features of OCs are extracted based on a stacked denoising autoencoder (SDAE). Meanwhile, the symmetric electrical distance matrix is modified to represent fault locations in local regions, which reduces the redundancy of fault location features.

(3)

Multiple DSA models are constructed for fault locations in each local region. The model structure is the radial basis function neural networks (RBFNN) integrated with the Chebyshev distance, which can achieve a close assessment accuracy to complex models. Moreover, an incremental updating strategy is designed to enhance the adaptability of the DSA model.

The rest of this paper is organized as follows: Section 2 introduces a self-adaptive power grid partition method that strictly accounts for electrical connections. The input feature set for DSA is constructed in Section 3. Section 4 constructs DSA models for fault locations in each local region and designs an incremental updating strategy. Test results for a simplified provincial system and a larger-scale practical system in China are presented in Section 5. Conclusions are drawn in Section 6.

2. Power Grid Partition with a Strict Consideration of Electrical Connections

2.1. Description of the Power Grid Partition Problem for DSA

The impacts of different fault locations on the power system dynamic security have local similarities. That is, fault locations with tight electrical connections have similar effects on the dynamic security and vice versa. The tightness of electrical connections is typically quantified as the impedance-based electrical distance [29,30]. Therefore, the nodes with close electrical distances should be within the same region, while nodes with distant distances are in separate regions. Constructing dedicated models for each region can reduce the difficulty of training DSA models.

The impedance-based electrical distance between nodes in a power grid is typically defined as the equivalent impedance [30]:

$$d_{ij}=Z_{ii}+Z_{jj}-2Z_{ij}$$

(1)

where d_ij represents the electrical distance between nodes i and j; Z_ii and Z_jj denote the self-impedances of nodes i and j, respectively; and Z_ij denotes the mutual impedance between nodes i and j.

To quantify the quality of the power grid partition, an evaluation index of electrical modularity Q_e is defined as Equation (2) based on complex network theory and the concept of weighted modularity [31].

$$Q_{\mathrm{e}}={\displaystyle \frac{1}{2M}}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n\left(w_{ij}-\frac{W_i{W}_j}{2M}\right)}}·\delta \left(i,j\right)$$

(2)

$$M={\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{w}_{ij}}}$$

(3)

where n denotes the number of nodes in the power grid; W_i denotes the sum of edge weights of node i; δ(i, j)∈{0, 1} denotes whether nodes i and j belong to the same region; δ(i, j) = 1 indicates the same partition, while 0 indicates different regions; and M denotes the sum of the edge weights among all nodes. The edge weight w_ij between nodes i and j is represented by electrical distance. A higher Q_e indicates a better partition quality.

Clustering algorithms are common methods for partitioning power grids [32,33], which can be categorized into prototype, density, and hierarchical clustering. Prototype clustering, e.g., the K-means algorithm, requires a predetermination of the partition number [34], while hierarchical clustering has high computational complexity in large-scale power grids. Without specifying the partition number, the mean shift algorithm in density clustering can achieve a self-adaptive power grid partition with low computational costs [35,36]. Therefore, the mean shift algorithm is adopted in this paper.

Considering the characteristics of the DSA, the symmetric electrical distance matrix and Chebyshev distance are integrated into the mean shift algorithm for the power grid partition. The mean shift-based power grid partition is introduced from four aspects: feature representation, algorithm structure, partition steps, and distance metric.

2.2. Feature Representation of the Node Location Based on a Symmetric Electrical Distance Matrix

Clustering algorithms rely on input features to distinguish electrical differences among node locations. Therefore, the feature representation of a node location is supposed to reflect these electrical differences.

The symmetric electrical distance matrix, denoted as D, is composed of the electrical distance between any two nodes. It can be represented as follows:

$$D=\left(\begin{array}{ccc}{d}_{11}& \dots & d_{1n}\\ ⋮& \ddots & ⋮\\ d_{n1}& \cdots & d_{nn}\end{array}\right)$$

(4)

D is a n-order and symmetric matrix due to the symmetry of the electrical distance (i.e., d_ij = d_ji). Each row of the matrix D is a feature vector of the corresponding node location. For example, the location of node i can be represented as an n-dimensional feature vector x_i = (d_i1, d_i2, …, d_in), which is called the “electrical coordinate” of node i. Ref. [24] pointed out that the electrical coordinate can effectively reflect the electrical differences among nodes. Therefore, the electrical coordinates are selected as the node location features for the mean shift algorithm.

The row vector can be further called the full-dimensional electrical coordinate since the matrix D in Equation (4) is a complete symmetry matrix of the order n × n. If D is an incomplete and unsymmetric matrix of the order n × k (k < n), the feature vector of node i is a k-dimensional electrical coordinate, denoted as x_i = (d_i1, d_i2, …, d_ik). The electrical coordinate is composed of the electrical distances to preselected reference nodes. Its dimension is equal to the number of reference nodes. The reference nodes of the full-dimensional electrical coordinate include all nodes in the whole power grid.

The symmetric electrical distance matrix is sensitive to changes in the power grid topology. It needs to be recalculated when the outage of a branch occurs. The common recalculation method is to modify the node admittance matrix and obtain the new node impedance matrix by inversion. The new symmetric electrical distance matrix is obtained by Equation (1). The computation time is high due to matrix inversion. To improve the recalculation efficiency, the branch addition method is adopted.

As shown in Figure 1, Z₀ is the line impedance of line l_ij. Disconnecting line l_ij is equivalent to adding a new line of −Z₀ to line l_ij. The new node impedance matrix Z′ can be represented as follows:

$$Z^{\prime }=Z-{\displaystyle \frac{Z{A}_{\mathrm{M}}{\left(Z{A}_{\mathrm{M}}\right)}^T}{A_{\mathrm{M}}^{\mathrm{T}}Z{A}_{\mathrm{M}}-Z_0}}$$

(5)

where A_M is a n-order column vector, which can be represented as [0, …, 1, …, −1, …,0]. Its i-th element is 1, and the j-th element is −1. Without the complex matrix inversion, the new node impedance matrix can be obtained rapidly. The recalculation efficiency of the symmetric electrical distance matrix is thereby improved significantly.

2.3. Mean Shift Algorithm

The underlying assumption of the mean shift algorithm is that sample points in different clusters follow different probability density distributions [37,38]. Each sample point performs iterations along the direction of the steepest increase in density and converges to the region with the highest local density. Sample points converging to the same region belong to the same cluster. With the two-dimensional feature plane in Figure 2 as an example, the iteration process of density centers for sample points is elaborated as follows. The colored circle represents the distribution of the sample point in the feature plane. The arrow represents the direction of the shift vector.

A sample point x is randomly selected as the initial density center. Based on all sample points within a circle centered at x with radius r, the shift vector xx⁽¹⁾ is calculated. The direction of the shift vector is toward the steepest increase in the density of sample distribution. The endpoint x⁽¹⁾ of the vector xx⁽¹⁾ is selected as the new density center, and a new shift vector x⁽¹⁾x⁽²⁾ is calculated. The above process is repeated until the final shift vector ||x^(m−1)x^(m)|| is less than a predefined threshold. The iteration process converges, and the endpoint x^(m) is considered the converged density center for sample point x. The shift vector M(x) for sample point x is calculated as follows [35]:

$$M(x)=\left[{\displaystyle \sum _{x_i\in S_x}G\left(\frac{x_i-x}{\sigma }\right)\left(x_i-x\right)}\right]/{\displaystyle \sum _{x_i\in S_x}G\left(\frac{x_i-x}{\sigma }\right)}$$

(6)

$$G\left({\displaystyle \frac{x_i-x}{\sigma }}\right)={\displaystyle \frac{1}{\sqrt{2\mathsf{\pi }}\sigma }}{e}^{-{\displaystyle \frac{{‖x_i-x‖}^2}{2{\sigma }^2}}}$$

(7)

where x represents the feature vector of the sample point x, i.e., electrical coordinate; S_x denotes the high-dimensional spherical region around the sample point x; G((x_i-x)/σ) represents the Gaussian weight [39]; and σ is the bandwidth of the Gaussian weight. After one shift, the feature vector of the density center point x⁽¹⁾ can be represented as follows:

$$x^{(1)}=x+M(x)={\displaystyle \sum _{x_i\in S_x}\left[G\left(\frac{x_i-x}{\sigma }\right)x_i\right]}/{\displaystyle \sum _{x_i\in S_x}G\left(\frac{x_i-x}{\sigma }\right)}$$

(8)

The density center iteration is performed for all sample points until the iteration process of all sample points converges. The clusters are determined according to the converged density centers of sample points. Sample points converging to the same region are divided into the same cluster, while those in different regions are divided into different clusters.

2.4. Self-Adaptive Power Grid Partition Based on the Mean Shift Algorithm

The power grid partition based on the mean shift algorithm is self-adaptive without specifying the partition number beforehand. The detailed partition steps are as follows:

(1)

Compute symmetric electrical distance matrix, extract the electrical coordinate of each node, and set the iteration termination threshold d_s and aggregation threshold d_a;

(2)

Calculate the shift vector M(x) of each node;

(3)

Move each node according to its shift vector M(x);

(4)

Repeat Steps (2) and (3) until ||M(x)||<d_s and the density center iteration process for each node has converged;

(5)

If the distance between the converged density centers of two nodes is less than d_a, these two nodes are aggregated to the same cluster, i.e., the same region; otherwise, these two nodes are assigned to different regions.

Step (1) involves feature input and parameter setting. Steps (2) to (4) constitute the density center iteration process for each node. Step (5) is the node aggregation process. To ensure the topological connectivity of each local region after partitioning, a hierarchical aggregation method is adopted in Step (5).

The node aggregation process considering topological connectivity is shown in Figure 3. The point denotes the node; the blue point denotes the aggregation center node. lines represent topological connections between nodes, numbers on the points indicate the connectivity level of each node to the aggregation center node, green denotes division into the same region, and white denotes separation into different regions. One-level connectivity indicates a direct connection, while two-level connectivity indicates connection through one intermediate node. Firstly, as shown in Figure 3a, a node is randomly selected as an aggregation center node, i.e., the blue point. Nodes that are 1-level connected to the blue point are identified based on the adjacency matrix. If the distance between the converged density centers of the 1-level connected node and the blue point is less than d_a, these two nodes are partitioned into the same region; otherwise, these two nodes belong to different regions. Then, as shown in Figure 3b, 2-level connected nodes that are directly connected to 1-level connected nodes in the current region are identified. The distance between the converged density centers of 2-level connected nodes and the blue point is calculated, and the nodes within the same region are screened. The node aggregation process continues until no more higher-level connected nodes within the same region can be found. The above process is performed repeatedly for another unpartitioned node until all nodes are partitioned. The aggregation process of the power grid is terminated.

In large-scale power systems, the computation cost of the mean shift algorithm is high, and the power grid partition efficiency is decreased. To boost the partition efficiency, this paper makes the following two improvements: (1) Reduce the number of nodes to be clustered. The generator bus belongs to a one-degree node, which is only connected to one terminal bus. To enhance partition efficiency, generator buses do not participate in clustering. They are directly partitioned into the same region as the connected terminal buses. (2) Adopt parallel computation. In the mean shift algorithm, the density center iteration for each node is independent of others and only depends on the initial node distribution. It can be accelerated by parallel computation and is not affected by each other.

During the clustering process, small regions consisting of a few nodes may be formed. Constructing a DSA model for each small region can increase the number of models, occupying many computation and storage resources. To reduce the number of DSA models, the partition results of the power grid need to be fine-tuned. Small regions are aggregated into the adjacent regions.

2.5. Chebyshev Distance between Node Location Features

In clustering algorithms, the distance metric between features, a.k.a. feature distance, is a crucial factor affecting the quality of power grid partition. To achieve a power grid partition that strictly considers the electrical distance, the feature distance is supposed to reflect the electrical distance between nodes. Typically, clustering algorithms use the Euclidean distance to measure distances between sample points. However, the Euclidean distance between electrical coordinates is not equivalent to the electrical distance between nodes. In this paper, the Chebyshev distance is selected as the distance metric for further improving the quality of the power grid partition.

With electrical coordinates, the feature distance dist(x_i, x_j) between nodes i and j can be represented in the following norm form:

$$dist\left(x_i,x_j\right)={‖x_i-x_j‖}_p,p>0$$

(9)

where p represents the norm, p = 2 is the Euclidean distance, and p = ∞ is the Chebyshev distance. As shown in Equation (9), the feature distance dist(x_i, x_j) is determined by the electrical coordinate x and p-norm. Ref. [40] has proven that the Chebyshev distance between full-dimensional electrical coordinates is equivalent to the electrical distance. That is, Equation (10) holds.

$$dist\left(x_i,x_j\right)={‖x_i-x_j‖}_{\infty }=\max \left(d_{i1}-d_{j1},d_{i2}-d_{j2},\cdots ,d_{in}-d_{jn}\right)=d_{ij}$$

(10)

Therefore, full-dimensional electrical coordinates and the Chebyshev distance are employed in the mean shift algorithm to achieve the power grid partition by strictly considering electrical connections. The application of the Chebyshev distance in the mean shift algorithm includes the determination of the high-dimensional spherical region in Equation (6), the distance metric in Equation (7), the length calculation of ||M(x)|| in Step 4 of Section 2.4, and the distance metric between converged density centers in Step 5 of Section 2.4.

3. Construction of the Input Feature Set for DSA

Multiple DSA models are constructed for the anticipated fault locations in each local region. For a specific fault type and a power grid with a fixed topology, the dynamic security of the power system is determined by pre-fault OCs and fault locations. Input data must include OCs and fault location information. In this paper, pre-fault OCs are represented by steady-state power flow features, which are widely used by existing research. Anticipated fault locations are represented by electrical coordinates.

3.1. Feature Extraction of OCs Based on SDAE

The steady-state power flow features are used to represent pre-fault OCs. These features include the active and reactive power of each generator, the active and reactive power of each load, the voltage magnitude and phase angle of each bus, and the active and reactive transmission power on each line. It is noted that the OCs are intended to describe the operating status of the whole power system rather than a local region. Therefore, steady-state power flow features across the whole power system are adopted. These features are generally high-dimensional and exhibit feature redundancy, deteriorating the training accuracy and efficiency of the data-driven DSA model. Thus, high-level features must be extracted from high-dimensional steady-state features.

SDAE is an unsupervised multilayer neural network that extracts abstract high-level features by stacking multiple denoising autoencoders (DAEs) [41,42]. It has a robust feature extraction capability. In this section, SDAE is adopted as an extractor for steady-state power flow features, and the extracted high-level features are used as inputs of each DSA model for different local regions.

Each DAE in the SDAE consists of an encoder and a decoder. The training objective is to minimize the error between the input of the encoder and the reconstruction output of the decoder. The loss function L_D of DAE can be represented as follows [43]:

$$\min L_D={\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^N{‖F_i-z_i‖}^2}$$

(11)

where N represents the number of training samples; the i-th training sample feature vector, denoted as F_i, serves as the encoder input; and z_i is the reconstructed feature vector, which serves as the decoder output. The encoder output is the high-level features extracted by DAE. In the SDAE, the DAE in the first layer uses steady-state power flow features as input for training. The outputs of the first layer serve as inputs of the next layer. The SDAE is trained layer by layer. The outputs of the final layer serve as the high-level features, which are used as one part of the DSA model inputs.

3.2. Redundancy Analysis of Electrical Coordinates for Fault Locations in a Local Region

The full-dimensional electrical coordinate of the order n is used to represent the node location in the mean shift-based power grid partition. However, it is redundant for DSA to represent fault locations in the local region. As shown in Figure 4, Regions A and B are two interconnected local regions with l tie lines. A_i and B_i are the nodes at both ends of the tie lines. The electrical distance between nodes P₁ and Q₁ is calculated as follows:

$$d_{{\mathrm{P}}_1{\mathrm{Q}}_1}=Z_{{\mathrm{P}}_1{\mathrm{P}}_1}+Z_{{\mathrm{Q}}_1{\mathrm{Q}}_1}-2Z_{{\mathrm{P}}_1{\mathrm{Q}}_1}$$

(12)

Based on the physical meanings of the mutual impedance, Z_P1Q1 is numerically equal to the voltage at node P with a unit current injected into node Q₁ and no current injected into other nodes, denoted as Z_P1Q1 = ${\dot{U}}_{{\mathrm{P}}_1}{|}_{{\dot{I}}_{Q1}=1}$. With a unit current injected into node Q₁, the currents on tie lines are denoted as k₁, k₂, …, k_l. According to the superposition theorem, the voltage at node P₁ is equal to the sum of the voltages at node P₁ with region B removed and currents k₁, k₂, …, k_l, respectively, injected at nodes A₁, …, A_l. It can be represented as follows:

$$Z_{{\mathrm{P}}_1{\mathrm{Q}}_1}={\dot{U}}_{{\mathrm{P}}_1}{|}_{{\dot{I}}_{Q1}=1}={\displaystyle \sum _{i=1}^l{U_{{\mathrm{P}}_1}}^{\prime }{|}_{{\dot{I}}_{\mathrm{A}i}=k_i}}={\displaystyle \sum _{i=1}^l{k}_i}{Z_{{\mathrm{P}}_1{\mathrm{A}}_i}}^{\prime }$$

(13)

where Z_P1Ai′ represents the mutual impedance between node P₁ and A_i after region B and the tie lines are removed. Substituting Equation (13) into Equation (12) yields the following equation:

$$d_{{\mathrm{P}}_1{\mathrm{Q}}_1}=-2{\displaystyle \sum _{i=1}^l{k}_i}{Z_{{\mathrm{P}}_1{\mathrm{A}}_i}}^{\prime }+Z_{{\mathrm{P}}_1{\mathrm{P}}_1}+Z_{{\mathrm{Q}}_1{\mathrm{Q}}_1}$$

(14)

Assume there are m reference nodes located in region B, denoted as Q₁, …, Q_m. The currents on tie lines are denoted as k_i1, k_i2, …, k_il when a unit current is injected into node Q_i. The m-dimensional electrical coordinate composed of the electrical distances between node P and m node Q_i, denoted as x_p, can be represented as follows:

$$x_{\mathrm{P}}={\left[\begin{array}{cccc}{d}_{{\mathrm{PQ}}_1}& d_{{\mathrm{PQ}}_2}& \cdots & d_{{\mathrm{PQ}}_m}\end{array}\right]}^{\mathrm{T}}=\left(\begin{array}{cccc}-2k_{11}& \cdots & -2k_{1l}& 1\\ ⋮& \ddots & ⋮& ⋮\\ -2k_{m1}& \cdots & -2k_{ml}& 1\end{array}\right)\left(\begin{array}{c}{Z_{{\mathrm{PA}}_1}}^{\prime }\\ {Z_{{\mathrm{PA}}_2}}^{\prime }\\ ⋮\\ {Z_{{\mathrm{PA}}_l}}^{\prime }\\ Z_{\mathrm{PP}}\end{array}\right)+\left(\begin{array}{c}{Z}_{{\mathrm{Q}}_1{\mathrm{Q}}_1}\\ Z_{{\mathrm{Q}}_2{\mathrm{Q}}_2}\\ ⋮\\ Z_{{\mathrm{Q}}_m{\mathrm{Q}}_m}\end{array}\right)$$

(15)

Equation (15) can be rewritten into the following form:

$$x_{\mathrm{P}}=K_{(m\times \left(l+1\right))}·Z_{(l+1)}+C_m$$

(16)

where K_(m×(l+1)) is a constant matrix of order m × (l + 1); C_m is an m-dimensional constant vector; and Z_(l+1) is an (l + 1)-dimensional vector related to nodes P, A₁, …, A_l. The m-dimensional electrical coordinate of the node in region A is only related to Z_(l+1). Due to the sparse electrical connections between regions A and B, the number of tie lines is small. The feature dimension m is greater than l + 1, resulting in feature redundancy.

To reduce the feature redundancy, the electrical coordinate reference nodes of the fault location in the local region should be selected within that region. The symmetric electrical distance matrix is modified as shown in Equation (17), where D_A and D_B are used to represent the fault locations occurring at regions A and B.

$$D=\left(\begin{array}{cc}{D}_{\mathrm{A}}& 0\\ 0& D_{\mathrm{B}}\end{array}\right)=\left(\begin{array}{cc}\begin{array}{ccc}{d}_{{\mathrm{A}}_1{\mathrm{A}}_1}& \dots & d_{{\mathrm{A}}_1{\mathrm{P}}_s}\\ ⋮& \ddots & ⋮\\ d_{{\mathrm{P}}_s{\mathrm{A}}_1}& \cdots & d_{{\mathrm{P}}_s{\mathrm{P}}_s}\end{array}& \\ & \begin{array}{ccc}{d}_{{\mathrm{B}}_1{\mathrm{B}}_1}& \dots & d_{{\mathrm{B}}_1{\mathrm{Q}}_m}\\ ⋮& \ddots & ⋮\\ d_{{\mathrm{Q}}_m{\mathrm{B}}_1}& \cdots & d_{{\mathrm{Q}}_m{\mathrm{Q}}_m}\end{array}\end{array}\right)$$

(17)

The dimension of steady-state power flow features can be significantly reduced by the feature extraction of SDAE. The dimension of electrical coordinates can also be reduced using block matrices from the symmetric electrical distance matrix to represent fault locations in local regions. Consequently, the dimension of input features for DSA does not drastically increase with the growth of the system scale. The proposed input feature representation method for DSA is applicable in large-scale power systems.

4. Dynamic Security Partition Assessment

The DSA model is constructed based on the RBFNN due to its flexible distance metric. It can leverage the electrical coordinates and Chebyshev distance to distinguish similar samples. To improve the adaptability of DSA model, an incremental updating strategy is designed to ensure its high assessment accuracy for different non-linear transient changes.

4.1. Construction of the RBFNN-Based DSA Model for a Local Region

The anticipated fault locations in the same region have similar impacts on the system dynamic security. The nonlinearity of the input–output map within the training set for each local region is reduced. The map is suitable to be learned by shallow models. RBFNN is a three-layer feedforward neural network with a single hidden layer. It has a simple structure, high training efficiency, and an interpretable decision-making process [44]. It can approximate any continuous function with arbitrary precision. Therefore, multiple DSA models are constructed based on the RBFNN for the anticipated fault locations in each local region.

As shown in Figure 5, the input to the RBFNN consists of OC and anticipated fault location features. Each neuron in the RBFNN hidden layer corresponds to a center vector obtained from the training set by the clustering algorithm. The radial basis function Φ() measures the difference between the sample and the center vector, which is the activation function of the hidden layer neurons. The Gaussian function is a widely used radial basis function, which can be represented as follows [39]:

$$\varphi \left(F_i,F_{\mathrm{c}}\right)={\displaystyle \frac{1}{\sqrt{2\mathsf{\pi }}\sigma }}{e}^{-{\displaystyle \frac{dist{(F_i,F_{\mathrm{c}})}^2}{2{\sigma }^2}}}$$

(18)

where F_i is the feature vector of the i-th sample and F_c is the center vector of the hidden layer neuron.

The feature vector is defined as F = (o, x), where o are high-level steady-state features of OCs, and x is the electrical coordinate of the anticipated fault location. The OC features represented by steady-state features are typically treated as Euclidean data, which are measured by the Euclidean distance. The anticipated fault location features, represented by electrical coordinates, are non-Euclidean data, which are measured by the Chebyshev distance. Considering these two types of feature distance metrics, dist(F_i, F_c) is defined as follows:

$$dist(F_i,F_{\mathrm{c}})=\alpha {‖o_i-o_{\mathrm{c}}‖}^2+\beta {‖x_i-x_{\mathrm{c}}‖}^{\infty }$$

(19)

where α and β are weight coefficients and α + β = 1.

Both the mean shift-based grid partition and RBFNN employ the Chebyshev distance, which focuses on the maximum value and disregards other values. The mean shift algorithm and RBFNN expect to measure differences between locations according to electrical distance, and the Chebyshev distance between full-dimensional electrical coordinates is equivalent to the electrical distance. Therefore, the Chebyshev distance is sufficient in this paper.

The Chebyshev distance between electrical coordinates is impacted by the power grid topology. If there is an island region in the power grid, the electrical coordinate and Chebyshev distance will be infinite. Both the mean shift algorithm and RBFNN cannot handle infinite feature values. To ensure the feasibility of the methods proposed in this paper, island regions need to be identified based on the adjacency matrix and removed.

4.2. Online Incremental Update of the DSA Model

The typical mode of current data-driven DSA methods is “offline training, online assessment”. However, the OCs of the power grid are constantly changing, and the training set generated offline cannot fully cover actual OCs and anticipated fault locations. During the online assessment, significant differences between the samples to be evaluated and the training set may lead to low generalization accuracy of the DSA model. To enhance the assessment accuracy of the DSA model, it is significant to generate new training samples to update the DSA model online. This section designs an online incremental updating method of the DSA model based on distance metrics and RBFNN. The new sample set for model updating is denoted as X_new. The detailed model updating steps are as follows:

(1)

The new OCs and anticipated fault locations are obtained based on real-time and future prediction information. High-level steady-state features and electrical coordinates are extracted, and an unlabeled new sample set is generated.

(2)

The minimum feature distance between the new sample and the RBFNN hidden layer center vectors is calculated according to Equation (19). If the minimum exceeds the predefined threshold D_t, a dynamic security label of the new sample is obtained based on a time-domain simulation and stored in X_new. If the minimum is lower than D_t, the new sample is abandoned.

(3)

Central samples of X_new are obtained based on the mean shift algorithm and Equation (19).

(4)

New neurons in the RBFNN hidden layer are added with central samples from X_new as the center vectors. The weights between the new neurons and the output layer are trained based on X_new, while the other structures and weights remain unchanged.

With the above incremental updating method, the DSA model can efficiently acquire knowledge from the new sample set. Moreover, it can avoid catastrophic forgetting during the model updating process. By repeatedly performing incremental updates, the rolling update of the DSA model can be achieved, ensuring that the model consistently maintains high generalization accuracy.

4.3. Framework of the Dynamic Security Partition Assessment

The mode of “offline training, rolling update, online assessment” for the dynamic security partition assessment is adopted. The assessment framework is illustrated in Figure 6. It mainly includes three blocks: offline training, rolling updating, and online assessment. Each block is elaborated as follows:

(1)

Offline training provides initial training data and a model for DSA. The power grid is firstly partitioned into several regions based on the mean shift algorithm. Multiple training sets for each region are generated, which include high-level steady-state features extracted by SDAE and a block symmetric electrical distance matrix. Multiple DSA models are constructed based on the RBFNN and each training set. The OC features of the RBFNN models in each region are the same.

(2)

Rolling updating periodically performs an update of DSA models to guarantee high DSA performance all the time. During the rolling updating stage, new training samples are generated periodically. DSA models are updated online and incrementally based on the proposed updating method. The number of DSA models is not increased, while catastrophic forgetting can be avoided. In this block, T represents the model updating cycle.

(3)

Online assessment gives the final assessment results. Combinations of OCs and anticipated faults are obtained, and the region where the fault is located is first determined. Next, the steady-state features and electrical coordinates of the combination are extracted. Finally, the DSA model corresponding to the fault location is adopted to obtain the assessment results.

The quality of input data is crucial for the data-driven DSA. Poor data may arise in two stages. (1) In the offline training stage, insufficient training samples and class imbalance may occur. To generate sufficient samples for model training, a wider range of OCs and more anticipated fault locations can be simulated by a time-domain simulation. The oversampling technique can be employed to increase the number of samples in minority classes. (2) In the online assessment stage, the acquired OC data may contain noise. This paper employs SDAE to process OC features, which exhibit strong robustness to noise. The high-level features extracted by SDAE are used as inputs of the RBFNN, reducing the impact of noise on the assessment results of the RBFNN.

The comprehensive flow chart of the proposed data-driven dynamic security partition assessment is shown in Figure 7. It includes the power grid partition, sample generation, DSA model construction, model updating, and online assessment. The DSA of the power system depends on communication and big data technology. It is vulnerable to false data injection attack (FDIA), which injects false information into the cyber system of power grid and can result in the misjudgments of DSA model. To reduce the risk of FDIA, false data detection is performed for new training samples before the model update. The false data are identified and eliminated, and then the DSA model is updated online for new preprocessed training samples.

The power system DSA includes transient rotor angle stability, voltage stability, and frequency stability assessments. The proposed partition assessment method is applicable to the above three scenarios. In this paper, the transient rotor angle stability assessment is adopted as the application scenario. The output of each regional DSA model provides the transient angle stability results of the whole power system under a combination of OCs and anticipated fault locations. The power system is defined as unstable when the maximum rotor angle difference between synchronous generators exceeds 180° during the transients.

The power grid is partitioned into multiple regions and dedicated DSA models are constructed for fault locations in each region. Each DSA model requires the input of global OC information, and the output is the global dynamic security results for the entire power system. These global data are typically consolidated at the provincial power dispatch center and not shared among the different grid regions. Therefore, the construction, storage, calling, and updating of DSA models are usually carried out at the provincial power dispatch center. The central computing resources are limited. When central computing resources are insufficient, distributed computing technology can be employed. The power dispatch center assigns updating tasks for DSA models to the computing resources in each region, and then the updated model parameters are sent back to the power dispatch center.

5. Experimental Results and Discussions

The effectiveness of the proposed dynamic security partition assessment method is verified in a simplified provincial power system in China. The example system includes 112 buses, 32 generators, and three HVDC links. The transmission power of three HVDC links are 8 GW, 8 GW, and 4 GW. Nine wind farms are integrated into the system. The simulation software is PSSE 34. Numerical simulation and machine learning are performed on a server with 64 GB of RAM and a 24-core Intel Xeon Gold 6148 CPU (2.40 GHz).

5.1. Test of Power Grid Partition

The power grid is partitioned based on the mean shift algorithm integrated with the symmetric electrical distance matrix and Chebyshev distance. The symmetric electrical distance matrix is a 112-order matrix. The radius of the high-dimensional spherical region is set to the average value of the symmetric electrical distance matrix. The bandwidth σ in the Gaussian weight is set to 20. The iteration termination threshold d_s is set to 0.01. The aggregation threshold d_a is a key factor in the partition quality. The electrical modularity Q_e under different aggregation thresholds d_a is shown in Figure 8. The Q_e reaches the maximum when d_a = 1.6. Therefore, d_a is determined to be 1.6.

The number of buses in different regions before fine-tuning is shown in Figure 9a. Different colors represent different regions. “32.1% (36)” denotes that the region includes 36 buses and accounts for 32.1% of the total bus count. During fine-tuning process, three small regions where the node number is less than five are aggregated into adjacent regions. The number of buses in different regions after fine-tuning is shown in Figure 9b. Three small regions are aggregated into adjacent regions. The distributions of the remaining five regions are shown in Figure 10. Different colors represent different regions. The number represents the region number.

To verify the effectiveness of the Chebyshev distance, mean shift algorithms with different distance metrics are compared in Figure 11. The distance metrics include Chebyshev, Euclidean, Manhattan, and cosine distances. The results show that the power grid partition quality with the Chebyshev distance metric is the best. It demonstrates that the Chebyshev distance is applicable to the mean shift-based power grid partition.

5.2. Test of Dynamic Security Partition Assessment

Multiple DSA models are constructed for anticipated fault locations in each local region. To generate the dataset for DSA, the time-domain simulations are performed as follows. The power of each load in the system fluctuates randomly between 80% and 120%. The outputs of the generators are adjusted based on load levels and spinning reserve capacity. A total of 1000 OCs are obtained. The anticipated faults are N-1 three-phase short circuits on 10%, 20%, …, 90% of each line. The fault duration is 0.25 s. For each fault location, 100 OCs are randomly selected. Following the above simulation settings, a dataset including 92,700 samples is generated, which includes 50,514 stable samples and 42,186 unstable samples. The number of anticipated fault locations and samples in each region is shown in

Table 1. The number of anticipated fault locations and samples in each region.

Region No.	The Number of Anticipated Fault Locations	The Number of Samples	The Number of Stable Samples	The Number of Unstable Samples
1	207	20,700	6709	13,991
2	144	14,400	5628	8772
3	27	2700	0	2700
4	261	26,100	17,687	8413
5	288	28,800	20,490	8310
Total	927	92,700	50,514	42,186

. It is worth noting that the proposed partition method discovers a special region 3, where no stable samples exist. It indicates that the instability risk of transient rotor angle caused by faults in region 3 is high. Due to the lack of stable samples, the faults in region 3 are considered to cause instability without building the specified DSA model.

The steady-state power flow features to represent OCs are extracted by SDAE. The structure of SDAE is (50, 40, 30). The dimension of the extracted high-level features is 30. The symmetric electrical distance matrix of each region is used to represent the regional fault locations. Except for region 3, the samples in other regions are randomly divided into training and test samples at 3:1, respectively. Then, four RBFNN-based DSA models are constructed on the training samples in regions 1, 2, 4, and 5. The hyper-parameters of RBFNN are shown in

Table 2. Hyperparameters of the RBFNN.

Parameter Type	Value	Parameter Type	Value
Hidden-layer neurons	100	Optimizer	Adam
−1/2σ2 of RBF kernel	1	Learning rate	0.001
Distance metric for OC features	Euclidean distance	Training epochs	100
Distance metric for fault location features	Chebyshev distance	Training strategy	Mini-batch
α in Equation (19)	0.5	Batch size	300
β in Equation (19)	0.5

. Center vectors of hidden-layer neurons are obtained from training samples based on the default K-means algorithm.

Figure 12 presents the comparison results of partition assessment and unified assessment with different machine learning models. The unified assessment constructs one model for fault locations in all regions. The machine learning models for comparison include RBFNN, CNN, deep belief network (DBN) [45], gated recurrent unit (GRU) [13], RNN [28], random forest (RF) [46], and extreme gradient boosting (XGBoost) [47]. The structures of the above machine learning models are shown in

Table 3. Structures of machine learning models.

Models	Structure
CNN	Conv (10, 3) + MaxPooling (2) + Conv (10, 3) + GlobalMaxPooling () + Dense (50)
DBN	(100, 100, 50)
GRU	GRU (50, 20) + Dense (50)
RNN	RNN (50, 20) + Dense (50)
RF	n_estimators = 100
XGBoost	n_estimators = 100

. As shown in Figure 12, the accuracy of the partition assessment is higher than the unified assessment, whatever the machine learning model is. It demonstrates that the proposed partition assessment method is effective.

Table 4. Comparison of different machine learning models for the partition assessment.

Region	RBFNN	CNN	DBN	GRU	RNN	RF	XGBoost
1	99.38%	99.04%	98.98	97.82	98.86	98.86	99.32
2	98.44%	98.44%	95.91	96.92	97.36	95.69	97.94
3	\	\	\	\	\	\	\
4	99.25%	99.13%	98.71	98.50	98.81	98.71	99.25
5	98.64%	98.54%	97.51	95.88	98.21	98.68	98.86
Total	98.96%	98.91%	97.94	97.25	98.40	98.25	98.93

further shows the comparison of the results from different machine learning models for the partition assessment. The accuracy of RBFNN is close to CNN, RNN, and XGBoost. It demonstrates that the constructed RBFNN-based DSA model can achieve a close assessment performance to complex models.

The model for fault locations in region 1 is denoted as “Model 1”. With Model 1 as an example, the effectiveness of the proposed incremental updating method of the DSA model is verified. A total of 4600 new samples are generated. The OCs and fault locations in new samples are not included in old samples. Samples are divided as shown in Figure 13. Old and new samples are divided into training and test samples at 3:1. The proposed incremental updating method is compared with retraining and continued training. Model 1 with incremental updating is trained on new training samples, and 26 new neurons in the hidden layer are added. The incremental updating is achieved by the freeze of a layer in Keras. The model with retraining is trained on old and new training samples. The training epochs are set to 100. The continued training is designed to train the initial model with new training samples. An assessment of the accuracy of the updated Model 1 on old and new test samples and updating time are shown in

Table 5. Comparisons of the retraining and incremental updating of Model 1.

Updating Strategy	Accuracy on Old Test Samples	Accuracy on New Test Samples	Updating Time
Retraining	99.04%	99.17%	500 s
Continued training	92.12%	99.22%	97 s
Incremental updating	99.46%	98.48%	30 s

.

The assessment accuracy with continued training on old test samples is significantly lower than the other two. It demonstrates that continued training leads to catastrophic forgetting of old knowledge. The accuracy with incremental updating on old test samples is the highest, while the assessment accuracy on new test samples is close to the other two. Moreover, the incremental updating process takes only 30 s. Test results demonstrate that the proposed incremental updating method has high updating performance and efficiency while avoiding catastrophic forgetting of old knowledge.

5.3. Tests on a Large-Scale Practical Power System in China

The practical power system is a regional system in China. It includes 2624 buses, 527 generator buses, 919 load buses, 2349 AC lines, two DC lines with 800 kV, eight DC lines with 500 kV, and 1634 transformers.

The power grid is partitioned based on the mean shift algorithm. The number of buses in different regions before fine-tuning is shown in Figure 14a. During the fine-tuning process, two small regions are aggregated into adjacent regions. The power grid is finally partitioned into three regions. The clustering results of buses in the three-dimensional feature space are shown in Figure 15, where d₁₀ represents the electrical coordinate feature of the tenth dimension. Different colors represent different clusters.

To generate the dataset, the time-domain simulations are performed as follows. The level of each load fluctuates randomly between 80% and 120%. The outputs of the generators are adjusted based on the load levels and spinning reserve capacity. A set of 7000 OCs is obtained. The anticipated faults are N-1 three-phase short circuits on each bus except for generator buses. The fault duration is set to 0.25 s. For each anticipated fault location, 50 OCs are randomly selected. A total of 104,850 samples are generated. The sample distribution in each region is shown in

Table 6. Distribution of samples in each region.

Region	Samples	Stable Samples	Unstable Samples
1	33,400	22,121	11,279
2	34,650	22,777	11,873
3	36,800	22,721	14,079
Total	104,850	67,619	37,231

.

Figure 16 presents a comparison of the results from partition and unified assessments with RBFNN. The accuracy of the partition assessment is significantly higher than that of the unified assessment. It demonstrates that the proposed dynamic security partition method is applicable to large-scale power systems.

6. Conclusions

This paper proposes a data-driven dynamic security partition assessment method for power systems, aiming to develop accurate and incrementally updated DSA models with simple structures. By the mean shift algorithm integrated with a symmetric electrical distance matrix and Chebyshev distance, the power grid is partitioned strictly considering electrical connections. Multiple DSA models are constructed based on the RBFNN and Chebyshev distance, achieving high assessment accuracy with a simple model structure. DSA models are updated online based on the incremental updating strategy, which achieves high updating efficiency and avoids catastrophic forgetting during the updating process. The case study demonstrates that the partition assessment has higher accuracy than the unified assessment. In the partition assessment, the RBFNN-based assessment models have close accuracy to deep learning and ensemble learning models. The proposed incremental updating method can greatly reduce the updating time while maintaining high accuracy.

In addition to pre-fault operating conditions and fault locations, the topology is also a significant factor affecting the dynamic security of the power system. In this paper, the input feature set of the DSA model does not include power grid topology information. The DSA model must be reconstructed to adapt to the change of the power grid topology. The topological change can be reflected by the symmetric electrical distance matrix. A feature representation method of the topology will be studied based on the symmetric electrical distance matrix. Furthermore, the construction of a low-complexity DSA model considering OCs, fault locations, and topology is also the focus of future research.