Power grid fault propagation path identification method and system based on graph neural network

Abstract

The application provides a power grid fault propagation path identification method and system based on a graph neural network, relates to the technical field of power grid fault identification, and comprises the following steps: obtaining power grid topology and fault observation data, constructing a direction perception graph neural network, realizing one-way message passing based on potential energy difference, obtaining node potential energy values and embedding representation, screening edges to generate a directed graph, and generating a candidate path along the potential energy gradient direction; the invariance variable is used to verify the graph neural network to evaluate the satisfaction degree of the path to physical constraints, and finally the path with the minimum topological entropy is selected. The application improves the accuracy and reliability of fault propagation path identification.

Description

Technical Field

[0001] This invention relates to power grid fault identification technology, and more particularly to a method and system for identifying power grid fault propagation paths based on graph neural networks. Background Technology

[0002] Traditional methods for identifying power grid fault propagation paths primarily rely on physical model analysis and expert judgment. Physical model analysis involves establishing a mathematical model of the power system and analyzing changes in system parameters before and after a fault; however, this method is computationally complex and requires high model accuracy. Expert judgment relies on human experience, analyzing fault records and alarm information to infer fault propagation paths; however, this method is inefficient and highly subjective. With the increasing scale and complexity of power grids, these traditional methods face growing challenges. In recent years, with the development of data acquisition technologies and artificial intelligence algorithms, data-driven power grid fault analysis methods have gradually emerged. Some researchers have attempted to use machine learning algorithms to analyze fault data, identify fault modes, and predict fault propagation. However, existing technologies still have the following shortcomings:

[0003] Existing methods are ineffective in representing the complex topology of power grids and its impact on fault propagation, and they neglect the interactions between nodes and network effects, resulting in low identification accuracy in complex network environments.

[0004] Existing methods lack effective modeling of fault propagation directionality and cannot accurately characterize the asymmetric propagation characteristics of faults in power systems, thus affecting the accuracy of fault propagation path identification.

[0005] Existing technologies often neglect the physical constraints of power systems, such as power conservation and time-series causality, when identifying fault propagation paths. This leads to discrepancies between the identification results and the actual operating laws of the physical system, reducing the reliability and practicality of the method. Summary of the Invention

[0006] This invention provides a method and system for identifying power grid fault propagation paths based on graph neural networks, which can solve the problems in the prior art.

[0007] A first aspect of this invention provides a method for identifying power grid fault propagation paths based on graph neural networks, comprising:

[0008] Acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node;

[0009] A direction-aware graph neural network is constructed. The fault feature vector is input into the direction-aware graph neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential value and node embedding representation of each node.

[0010] Based on the propagation potential energy value, edge filtering is performed on the graph structure to obtain a directed graph structure. Starting from the fault source node, a set of candidate propagation paths is generated by traversing the directed graph structure along the potential energy gradient direction.

[0011] An invariant verification graph neural network is constructed. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network. The output is a score of the degree of satisfaction of the physical invariant constraint set, which includes topological connectivity constraints, power conservation constraints, and temporal causality constraints.

[0012] Paths with satisfaction scores higher than a preset score threshold are selected, and the path with the minimum topological entropy is chosen as the fault propagation path.

[0013] undefined.

[0014] The steps for obtaining the propagation potential value and node embedding representation of each node through unidirectional message passing based on directional attention weights according to potential energy difference include:

[0015] The fault feature vector is decomposed into state components and time components, which are encoded separately and then fused through cross-attention to obtain the initial node embedding; the initial propagation potential value of the node is calculated based on the fault time.

[0016] For each edge in the graph structure, the direct potential difference, embedding space potential difference, and temporal potential difference are calculated based on the initial propagation potential energy value of the connected nodes, the initial node embedding, and the fault time component, respectively. These are then concatenated into a multi-view potential difference feature and subjected to nonlinear transformation to obtain the encoded potential difference representation.

[0017] The encoded potential difference is combined with the initial propagation potential value of the connected node to calculate the potential gradient attention score of the edge. The potential gradient attention score is then discretized and top-k sparsified using the Gumbel-Softmax method to obtain the directional attention weight.

[0018] The directional attention weights are applied to the initial node embedding of the source node to generate a weighted message vector. For each target node, the weighted message vectors of all incoming edges are aggregated and fused with the initial node embedding of the target node itself through a gated recurrent unit to obtain a node embedding representation. Based on the node embedding representation, the propagation potential value is calculated through a potential energy prediction network.

[0019] The steps of combining the encoded potential difference representation with the initial propagation potential value of the connected node to calculate the potential gradient attention score of the edge, and then differentiating and discretizing the potential gradient attention score using the Gumbel-Softmax method and performing top-k sparsification to obtain the directional attention weight include:

[0020] For each edge in the graph structure, the encoded potential difference is used to calculate the basic potential gradient score by bilinear interaction with the initial propagation potential value of the connected node. The potential energy directionality bias term is calculated based on the difference between the initial propagation potential values ​​of the source node and the target node. The potential energy gradient basic score and the potential energy directionality bias term are added together to obtain the potential energy gradient attention score of the edge.

[0021] The temperature decay coefficient is calculated based on the ratio of the current training round to the total training rounds. The initial temperature value is multiplied by the temperature decay coefficient to obtain the Gumbel temperature parameter. Random noise sampled from the Gumbel distribution is added to the potential energy gradient attention score to obtain the perturbed attention score. The perturbed attention score is divided by the Gumbel temperature parameter and normalized by Softmax to obtain a differentiable discretizable attention weight distribution.

[0022] For all incoming edges of each node, calculate the entropy value of the differentiable discretizable attention weight distribution. Based on the comparison result of the entropy value and the preset entropy threshold, determine the sparsity parameter k value. Select the k incoming edges with the largest attention weights and renormalize them to obtain the directional attention weights.

[0023] The steps of performing edge filtering on the graph structure based on the propagation potential energy value to obtain a directed graph structure, and generating a set of candidate propagation paths by traversing the directed graph structure along the potential energy gradient direction starting from the fault source node, include:

[0024] For each edge in the graph structure, calculate the propagation potential energy difference between the connected nodes, retain the edges with a propagation potential energy difference greater than zero, and assign the edge a direction attribute from the high potential energy node to the low potential energy node according to the potential energy gradient direction, thus forming a directed graph structure.

[0025] Starting from the fault source node as the current traversal node, select the node with the smallest propagation potential value among the outgoing edge neighbors of the current traversal node as the next hop traversal node, and add the directed edge between the current traversal node and the next hop traversal node to the path sequence.

[0026] When the propagation potential energy of the next hop traversed node is a local minimum or the path sequence length reaches the preset depth limit, the traversal stops and the path sequence is recorded as a candidate propagation path. For branch nodes with multiple outgoing neighbor nodes during the traversal, independent traversal is performed on each branch to generate the corresponding candidate propagation path, and the traversal results of all branches are summarized to form a set of candidate propagation paths.

[0027] The steps of inputting the subgraph structure of each path in the candidate propagation path set and the node embedding representation into the invariant verification graph neural network, and outputting a score on the degree of satisfaction of the physical invariant constraint set, wherein the physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints, include:

[0028] For each path in the candidate propagation path set, extract the subgraph structure formed by all nodes and edges in the path, and concatenate the node embedding representation of each node in the subgraph structure with the position code of the node in the path to obtain the position-aware node representation. Calculate the directional edge weight for the edges in the subgraph structure based on the difference between the initial propagation potential energy values ​​of the connected nodes and the electrical impedance parameters of the edges.

[0029] The location-aware node representation and the directional edge weights are input into the invariant verification graph neural network. The updated node representation is obtained by message passing through the graph convolutional layer. The updated node representation is then pooled along the path sequence direction to obtain the path-level representation vector.

[0030] Based on the path-level representation vector, the topological connectivity score, power conservation score, and temporal causality score are calculated respectively. The temporal causality score is calculated by verifying the temporal increasing relationship of the node failure time in the path. The topological connectivity score, power conservation score, and temporal causality score are weighted and fused to obtain the degree of satisfaction of the physical invariant constraint set.

[0031] The steps for calculating the satisfaction score of the set of physical invariant constraints include:

[0032] The path-level representation vector is used to predict the connection probability between each pair of adjacent nodes in the path. The connection probability is compared with the actual connection relationship of the corresponding node pair in the power grid topology adjacency matrix and aggregated to obtain the topology connectivity score.

[0033] The path-level representation vector is mapped to the predicted power inflow and outflow vectors of each node in the path. Based on the predicted power inflow and outflow vectors and the line impedance parameters of the edges between adjacent nodes in the path, the cumulative power loss value along the path is calculated. The cumulative power loss value is compared with the initial propagation potential energy value of the starting node of the path to obtain the power conservation score.

[0034] The topological connectivity score is used as a gating signal to filter the power conservation score to obtain the topologically constrained power score. The temporal causality score is used as a directional weight to weight and correct the topologically constrained power score to obtain the temporally corrected power score. The degree of satisfaction of the physical invariant constraint set is obtained by weighted fusion of the temporally corrected power score, the topological connectivity score and the temporal causality score.

[0035] A second aspect of the present invention provides a power grid fault propagation path identification system based on graph neural networks, comprising:

[0036] The first unit is used to acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node;

[0037] The second unit is used to construct a direction-aware graph neural network. The fault feature vector is input into the direction-aware graph neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential value and node embedding representation of each node.

[0038] The third unit is used to perform edge filtering on the graph structure based on the propagation potential energy value to obtain a directed graph structure, and to generate a set of candidate propagation paths by traversing the directed graph structure along the potential energy gradient direction starting from the fault source node;

[0039] The fourth unit is used to construct an invariant verification graph neural network. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network, and the output is a score of the degree of satisfaction of the physical invariant constraint set. The physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints.

[0040] The fifth unit is used to filter paths whose satisfaction rating is higher than a preset rating threshold, and select the path with the minimum topological entropy as the fault propagation path.

[0041] This invention utilizes a direction-aware graph neural network to perform unidirectional message passing based on the directional attention weights of potential energy difference. This effectively extracts the directional features of fault propagation along the network topology, significantly improving the shortcomings of traditional undirected graph models in accurately characterizing the directional propagation characteristics of power grid faults and enhancing the accuracy of fault propagation path identification.

[0042] This invention introduces an invariant verification graph neural network to evaluate candidate propagation paths. By simultaneously considering physical invariant constraints such as topological connectivity, power conservation, and temporal causality, it achieves a deep fusion of physical laws, effectively avoiding misjudgment problems caused by single constraints and enhancing the reliability and interpretability of the recognition results.

[0043] This invention uses the topological entropy minimization criterion to select the optimal propagation path, which conforms to the physical characteristic that electrical faults propagate along the shortest impedance path. While ensuring that the identification results meet physical constraints, it improves computational efficiency and can adapt to the real-time monitoring needs of large-scale power grids, providing strong support for power grid fault diagnosis and preventive maintenance. Attached Figure Description

[0044] Figure 1 This is a flowchart illustrating the power grid fault propagation path identification method based on graph neural networks according to an embodiment of the present invention.

[0045] Figure 2 The flowchart shows the verification and scoring mechanism for the multidimensional constraints of physical invariants in power grid fault propagation based on graph neural networks. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0047] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0048] Figure 1 This is a flowchart illustrating the power grid fault propagation path identification method based on graph neural networks according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0049] Acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node;

[0050] A direction-aware graph neural network is constructed. The fault feature vector is input into the direction-aware graph neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential value and node embedding representation of each node.

[0051] Based on the propagation potential energy value, edge filtering is performed on the graph structure to obtain a directed graph structure. Starting from the fault source node, a set of candidate propagation paths is generated by traversing the directed graph structure along the potential energy gradient direction.

[0052] An invariant verification graph neural network is constructed. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network. The output is a score of the degree of satisfaction of the physical invariant constraint set, which includes topological connectivity constraints, power conservation constraints, and temporal causality constraints.

[0053] Paths with satisfaction scores higher than a preset score threshold are selected, and the path with the minimum topological entropy is chosen as the fault propagation path.

[0054] For example, power grid topology data is obtained through a data interface and includes a node table and an edge table. The node table fields include node identifier, node type, and rated voltage level; the edge table fields include edge identifier, source node identifier, target node identifier, resistance value, and reactance value. The data is represented as a graph structure and stored using an adjacency list, which is a dictionary data structure where the key is the node identifier and the value is a list of adjacent edges. Fault observation data is obtained through a monitoring system message queue, and its fields include node identifier, fault status, fault time, fault type, and voltage and current before and after the fault. The fault time is a millisecond-level timestamp, and after receiving the data, time alignment is performed; the earliest fault time is subtracted from all fault timestamps to obtain the relative time.

[0055] Based on fault observation data, each node is assigned a fault feature vector, which is a 16-dimensional floating-point vector. The first four elements of the vector are the one-hot encoding of the fault type, the fifth element is the fault state encoding, the sixth element is the normalized fault time, the seventh and eighth elements are the normalized voltage and current difference, and the ninth to sixteenth elements are reserved extended fields filled with 0s. The fault source node is determined by traversing the nodes whose fault state is true, and selecting the node with the smallest normalized fault time.

[0056] A direction-aware graph neural network is constructed, which adopts a three-layer graph attention network architecture. The input layer receives a 16-dimensional fault feature vector, the input embedding layer is a fully connected layer with a weight matrix of dimension 16×128 and a bias vector of dimension 128. The activation function is a modified linear unit function, and the output is a 128-dimensional initial embedding representation of the nodes.

[0057] The core of the direction-aware graph neural network consists of a potential energy difference calculation module and a directional attention mechanism. The potential energy difference calculation module extracts the normalized fault times of the source and target nodes for each edge, calculates the time difference by subtracting the target node's time from the source node's time, and the potential energy difference is equal to the time difference multiplied by a potential energy scaling factor, which defaults to 10.0. A positive time difference indicates that the potential energy decreases along the edge direction, while a negative time difference indicates that the potential energy decreases in the opposite direction.

[0058] The directional attention mechanism processes all incoming edges for each node, defined as edges whose target node is the current node. For each incoming edge, a 128-dimensional node embedding representation of the source node and a 128-dimensional node embedding representation of the current node are extracted, and these two embedding vectors are concatenated to form a 256-dimensional vector. This concatenated vector is input to the attention computation layer, a single-layer fully connected network with a weight matrix of 256 x 1 dimensions, a default bias scalar of 0, a Leaky ReLU activation function, and a negative slope parameter of 0.2. The attention computation layer outputs a scalar value as the unnormalized attention score, which is multiplied by the potential difference value corresponding to the incoming edge to achieve directional modulation. Edges with negative or zero potential differences are masked, setting the directionally modulated attention score to negative infinity. A Softmax normalization function is applied to the directionally modulated attention scores of all incoming edges of the current node; the normalized values ​​are used as directional attention weights. Edges with a mask of negative infinity have their weights set to 0 after normalization, ensuring that messages are only transmitted in the direction of decreasing potential energy.

[0059] The one-way message passing mechanism uses directional attention weights multiplied by the 128-dimensional node embedding representation of the corresponding source node to obtain a weighted embedding vector. The weighted embedding vectors of all incoming edges of the current node are summed to obtain a 128-dimensional aggregated message vector. This aggregated message vector is concatenated with the node's own embedding representation to form a 256-dimensional vector, which is input to the node update layer. The node update layer is a fully connected layer with a weight matrix of 256 x 128 dimensions, a bias vector of 128 dimensions, and a modified linear unit activation function. The output is a 128-dimensional updated node embedding representation. The graph attention layer is executed three times. After each execution, a residual connection is applied, which adds the node embedding representation before and after the update and then divides the sum by 2.

[0060] The propagation potential prediction module receives the node embedding representation output from the third-layer graph attention layer. It includes a fully connected layer for dimensionality reduction mapping, with a weight matrix dimension of 128 x 1, a default bias scalar value of 0, and a linear activation function. The output scalar values ​​are truncated to the range of 0 to 10 using a truncation function, and the propagation potential value of the fault source node is forcibly set to the maximum value of 10.

[0061] Edge selection is performed based on propagation potential energy values. For each edge, the propagation potential energy value of the source node is subtracted from the propagation potential energy value of the target node to obtain the edge's potential energy difference. Edges with a potential energy difference less than or equal to a preset threshold are deleted (default threshold is 0.1). The retained edges are marked as directed edges. Starting from the fault source node, a depth-first search is used to traverse the directed graph structure along the potential energy gradient direction. Outgoing target nodes are sorted in ascending order of propagation potential energy values, prioritizing nodes with lower potential energy values. The maximum path length limit is set to 10 by default, generating a set of candidate propagation paths.

[0062] An invariant verification graph neural network is constructed, employing a three-layer graph convolutional network architecture. The input includes the subgraph structure and node embedding representations of the candidate propagation paths. The subgraph node feature matrix extracts the 128-dimensional node embedding representations from the output of the third layer of the orientation-aware graph neural network; the number of rows in the matrix equals the number of nodes in the path. The subgraph edge feature matrix extracts the potential energy difference, resistance, and reactance values ​​of the edges; the number of rows in this matrix equals the number of edges in the path, and the number of columns is 3.

[0063] The graph convolutional layer of the invariant verification graph neural network achieves neighborhood feature aggregation. In the neighborhood aggregation stage, for each node, its neighboring nodes in the subgraph are collected, and a 128-dimensional node embedding representation is extracted for each neighboring node. A 3-dimensional edge feature vector is extracted from the edge from the current node to its neighboring nodes. The 128-dimensional neighbor embedding representation and the 3-dimensional edge feature vector are concatenated to form a 131-dimensional vector, which is then input to the edge feature encoding layer. This edge feature encoding layer is a fully connected layer with a weight matrix of dimension 131 x 64, a bias vector dimension of 64, and a modified linear unit activation function. The output is a 64-dimensional edge encoding vector.

[0064] The aggregated feature vector is obtained by summing the 64-dimensional edge encoding vectors of all neighbors of the current node. This aggregated feature vector is then concatenated with the current node's own 128-dimensional node embedding representation to form a 192-dimensional vector. This concatenated vector is input to the node update layer, a fully connected layer with a weight matrix of 192 x 128 dimensions, a bias vector of 128 dimensions, and a modified linear unit activation function. The output is a 128-dimensional updated node embedding representation. The graph convolutional layer is executed three times. After each execution, a layer normalization operation is applied. Layer normalization calculates the mean and standard deviation for each of the 128 dimensions of the node embedding representation. The value of each dimension is subtracted from its mean and then divided by its standard deviation. If the standard deviation is less than 0.001, it is replaced with 0.001.

[0065] The output layer of the invariant verification graph neural network generates a score on the satisfaction of physical invariant constraints. The output layer contains three parallel score calculation branches, calculating the topological connectivity score, power conservation score, and temporal causality score, respectively. The topological connectivity score calculation branch traverses adjacent node pairs in the path, queries the connection state of the corresponding position in the adjacency matrix. A connection state of 1 indicates a connection consistency of 1 for that node pair, and 0 indicates a consistency of 0. The arithmetic mean of the connection consistency of all adjacent node pairs is used to obtain the topological connectivity score. The power conservation score calculation branch calculates the power imbalance by subtracting the total outflow power from the total inflow power for each intermediate node in the path. The absolute value of the power imbalance is divided by the total inflow power to obtain the power imbalance rate. The arithmetic mean of the power imbalance rates of all intermediate nodes is then subtracted from 1 to obtain the power conservation score. The temporal causality score calculation branch traverses adjacent node pairs in the path, calculates the time difference by subtracting the normalized fault time of the previous node from the normalized fault time of the subsequent node, and divides the number of node pairs with a time difference greater than 0 by the total number of adjacent node pairs to obtain the temporal causality score.

[0066] The satisfaction score fusion module receives scores from three branches: topological connectivity score (weighted by 0.3), power conservation score (weighted by 0.4), and temporal causality score (weighted by 0.3). These scores are then weighted and summed to calculate the overall score. Paths with satisfaction scores higher than a preset threshold (0.7 by default) are selected. The path with the lowest topological entropy is chosen as the fault propagation path. Topological entropy is equal to the standard deviation of the degrees of all nodes. The output is an ordered list of node identifiers.

[0067] For example, the power grid topology contains 8 nodes, N1 to N8, with edges from N1 to N2, N2 to N3, N2 to N4, N3 to N5, N4 to N6, N5 to N7, N6 to N7, and N7 to N8. Fault observation data shows that the fault time at node N1 is 0 seconds, N3 is 1.2 seconds, N5 is 2.5 seconds, and N7 is 3.8 seconds, thus identifying N1 as the fault source node. The direction-aware graph neural network calculates the propagation potential energy values ​​for node N1 as 10, N2 as 8.5, N3 as 7.2, N5 as 5.6, N7 as 3.9, and N8 as 2.1. Depth-first search generates candidate propagation paths including [N1, N2, N3, N5, N7, N8] and [N1, N2, N4, N6, N7, N8]. The invariant verification graph neural network calculates the satisfaction score of the first path as 0.91 and the satisfaction score of the second path as 0.79. The topological entropy of the first path is 0.82, and the topological entropy of the second path is 0.94. The first path is selected as the fault propagation path, and the output result is [N1, N2, N3, N5, N7, N8].

[0068] In one optional implementation, the step of assigning a fault feature vector containing the fault state and the fault time to a node based on the fault observation data and determining the fault source node includes:

[0069] Extract the protection device action information and circuit breaker tripping information of each node from the fault observation data, encode the protection device action information as a fault state characterizing the node, and mark the time corresponding to the circuit breaker tripping information as the fault time of the node; combine the fault state and the fault time to form a fault feature vector and assign it to the corresponding node.

[0070] Based on the fault time, a time-series sorting sequence is constructed, and the node with the earliest fault time in the time-series sorting sequence is identified as the fault source node.

[0071] For example, after acquiring power grid topology data and fault observation data, the protection device operation information and circuit breaker tripping information for each node are extracted from the fault observation data. The protection device operation information includes the relay protection device's operation type identifier, operating element number, operation criterion flag, and operation duration. This information is collected in real-time by a distributed data acquisition unit from the protection devices at each node of the power grid and uploaded to a centralized processing platform via a dedicated communication channel. The circuit breaker tripping information includes the circuit breaker status change timestamp, tripping reason code, and current amplitude before and after tripping. The timestamp is synchronously recorded using a unified clock source with millisecond-level precision to ensure the comparability of timing information between different nodes.

[0072] The action information of protection devices is encoded into fault states representing nodes. The encoding process establishes a mapping relationship between protection action types and fault state labels: overcurrent protection actions are encoded as state value 1, distance protection actions as state value 2, differential protection actions as state value 3, failure protection actions as state value 4, and no action as state value 0. When multiple protection devices operate simultaneously at a single node, a priority encoding rule is used: differential protection has the highest priority, followed by distance protection, then overcurrent protection, and finally failure protection has the lowest priority. The action type with the highest priority is selected as the fault state encoding value for that node. The encoding results are stored as an integer state vector, with each node corresponding to one state value.

[0073] The earliest tripping time is extracted from the circuit breaker state change timestamp sequence as the fault time of the node. When a node is connected to multiple circuit breakers and these circuit breakers trip at different times, the values ​​of all tripping timestamps are compared, and the time with the smallest timestamp value, i.e., the earliest tripping time, is selected as the fault time of the node. Timestamps are stored as long integer values ​​in milliseconds, using the UNIX timestamp format for easy cross-platform calculation and comparison. For cases where timestamps are abnormal due to communication delays or clock drift, a timestamp validity check threshold is set. If the difference between a node's tripping timestamp and the tripping timestamps of adjacent nodes exceeds a preset physical propagation delay upper limit, such as 500 milliseconds, a timestamp correction mechanism is triggered, and the fault time of the node is re-estimated using the weighted average of the timestamps of adjacent nodes.

[0074] The fault feature vector is constructed as a two-dimensional numerical array. The first dimension stores the fault state code value, and the second dimension stores the normalized value of the fault time. Normalization is achieved by subtracting the global minimum timestamp from the original millisecond-level timestamp and then dividing by the global timestamp range to obtain a relative time value between 0 and 1. This avoids instability in subsequent neural network calculations caused by excessively large absolute timestamp values. The feature vector assignment process uses the node identifier as the index key. The fault feature vector of each node is stored in a hash table data structure, with the key being a unique node identifier string and the value being a floating-point array containing the fault state and the normalized fault time. For nodes that have not experienced a fault or have not observed a protection action, their fault state code is assigned a value of 0, and their fault time is assigned the normalized value of 1 corresponding to the global maximum timestamp, indicating that the node did not participate in fault propagation during the observation period.

[0075] The time-series sorting process iterates through the fault time values ​​of all nodes, extracts nodes with non-zero fault state code values ​​as the set of valid fault nodes, and arranges the valid fault nodes in ascending order of their fault time values. The sorting algorithm is implemented using quicksort, with a time complexity of log-linear. The sorted results are stored as an ordered array of node identifiers. When multiple nodes have the same fault time, a secondary sorting rule is introduced, arranging them according to the lexicographical order of the nodes' unique identifiers in the power grid topology, ensuring the determinism and repeatability of the sorting results.

[0076] The identification process directly reads the first element of the time-series sorted sequence, and the node identifier corresponding to this element is the fault source node. The identifier of the fault source node is set as a global state variable for use in subsequent fault propagation path generation and verification processes. When the time-series sorted sequence is empty, meaning all nodes are inactive, it is determined that the current observation data does not contain valid fault information, the subsequent processing process is terminated, and an empty path result is returned. When the time-series sorted sequence contains multiple nodes, and the difference between the fault time of the first node and the fault time of the second node is less than a preset source node discrimination threshold, such as 50 milliseconds, all nodes with time differences within the threshold are marked as a set of candidate fault source nodes. Topological proximity analysis is used to select the node with the closest electrical distance as the primary fault source node, and the remaining nodes are marked as secondary fault source nodes for analysis of multi-source fault scenarios.

[0077] This invention makes full use of existing power grid monitoring data, requiring no additional equipment investment, thus improving the practicality and adaptability of the method. The mechanism for determining fault source nodes based on time-series sequencing conforms to the physical laws of power system fault propagation, making source identification more accurate and reliable.

[0078] In one optional implementation, the step of obtaining the propagation potential value and node embedding representation of each node through unidirectional message passing based on directional attention weights according to potential energy difference includes:

[0079] The fault feature vector is decomposed into state components and time components, which are encoded separately and then fused through cross-attention to obtain the initial node embedding; the initial propagation potential value of the node is calculated based on the fault time.

[0080] For each edge in the graph structure, the direct potential difference, embedding space potential difference, and temporal potential difference are calculated based on the initial propagation potential energy value of the connected nodes, the initial node embedding, and the fault time component, respectively. These are then concatenated into a multi-view potential difference feature and subjected to nonlinear transformation to obtain the encoded potential difference representation.

[0081] The encoded potential difference is combined with the initial propagation potential value of the connected node to calculate the potential gradient attention score of the edge. The potential gradient attention score is then discretized and top-k sparsified using the Gumbel-Softmax method to obtain the directional attention weight.

[0082] The directional attention weights are applied to the initial node embedding of the source node to generate a weighted message vector. For each target node, the weighted message vectors of all incoming edges are aggregated and fused with the initial node embedding of the target node itself through a gated recurrent unit to obtain a node embedding representation. Based on the node embedding representation, the propagation potential value is calculated through a potential energy prediction network.

[0083] For example, after obtaining the fault feature vector of each node, the fault feature vector is decomposed into state components and time components. The first dimension of the fault feature vector, i.e., the fault state code value, is extracted as the state component, and the second dimension, i.e., the normalized fault time, is extracted as the time component. The state component is encoded through an embedding layer, which maps the discrete fault state code values ​​into a continuous high-dimensional vector representation. The embedding dimension is set to 64, and the weight matrix of the embedding layer is initialized with a uniform distribution, with values ​​ranging from -0.1 to +0.1. The time component is encoded through a position encoder, which uses a sine-cosine encoding scheme to map the normalized fault time value into a periodic feature vector. The encoding dimension is also set to 64, and the period parameters of the sine and cosine functions are set to a negative power sequence of 10000 to ensure that time patterns of different frequencies are captured.

[0084] The encoded state and time components are fused through cross-attention to obtain the initial node embedding. The cross-attention mechanism uses the state component as the query vector and the time component as the key and value vectors, respectively. Attention weights are calculated by the dot product of the query and key vectors. The dot product result is scaled by dividing by the square root of the embedding dimension (8). The scaled value is then normalized using the Softmax function to obtain the attention distribution. The weighted sum of the attention distribution and the value vector produces the fused representation. This representation is added to the original state component through a residual connection and then normalized by a layer. The epsilon parameter for normalization is set to 0.000001 to avoid division by zero. The fused representation undergoes a nonlinear transformation through two fully connected layers. The first fully connected layer has an output dimension of 256 and uses GELU as the activation function. The second fully connected layer has an output dimension of 128, resulting in the initial node embedding. The initial node embedding is stored as a floating-point array, with each node corresponding to a 128-dimensional vector.

[0085] The initial propagation potential is calculated by taking the negative of the normalized fault time of the node and adding 1. That is, the initial propagation potential equals 1 minus the normalized fault time. This calculation method results in nodes with earlier fault times having higher initial propagation potential values; the initial propagation potential value of the fault source node is 1, and the initial propagation potential value of the latest fault node is close to 0. For nodes that have not experienced a fault, their normalized fault time is 1, and the calculated initial propagation potential value is 0, indicating that the node does not have the ability to propagate the fault outward. The initial propagation potential value is stored as a single-precision floating-point number, ranging from 0 to 1.

[0086] For each edge in the graph structure, the direct potential difference is calculated by subtracting the initial propagation potential value of the target node from the initial propagation potential value of the source node. A positive value indicates that the potential energy decreases from the source node to the target node, while a negative or zero value indicates that the potential energy does not conform to the propagation direction. The embedding space potential difference is calculated by the cosine similarity between the initial node embeddings of the source and target nodes. The cosine similarity is calculated by dividing the dot product of the two embedding vectors by the product of their respective magnitudes, and the value ranges from -1 to +1. The higher the similarity, the closer the potential energy distribution of the two nodes in the embedding space. The temporal potential difference is calculated by subtracting the normalized failure time of the target node from the normalized failure time of the source node. A negative value indicates that the failure time of the source node is earlier than that of the target node, which conforms to the causal propagation direction, while a positive or zero value indicates that the temporal relationship does not conform to the propagation logic.

[0087] The direct potential difference, embedded spatial potential difference, and temporal potential difference are concatenated to form a multi-view potential difference feature, which is then subjected to a nonlinear transformation to obtain the encoded potential difference representation. The three potential difference values ​​are concatenated to form a three-dimensional vector, which is then input into a fully connected network for encoding. The fully connected network consists of two layers: the first layer has an output dimension of 32 and uses ReLU activation, while the second layer has an output dimension of 16 and uses Tanh activation. The encoded potential difference is represented as a 16-dimensional floating-point vector, stored in an edge attribute hash table. The key is a unique identifier string for the edge, formed by concatenating the source node identifier and the target node identifier.

[0088] The potential energy modulation representation is obtained by element-wise multiplying the 16 dimensions of the encoded potential energy difference representation with the initial propagation potential energy value of the source node. This potential energy modulation representation is a 16-dimensional floating-point vector. The potential energy modulation representation is input to a single-layer fully connected network. The weight vector of the fully connected network is a 16-dimensional floating-point array, with the bias term set to 0. The network performs a dot product operation between the potential energy modulation representation and the weight vector to obtain a scalar output. The scalar output undergoes a nonlinear transformation using the Sigmoid function, compressing the result to the interval between 0 and 1. This result is the edge's potential energy gradient attention score. The potential energy gradient attention score reflects the activation strength of the edge during fault propagation; a higher score indicates that the edge is more likely to participate in the actual propagation path.

[0089] The Gumbel-Softmax method introduces Gumbel noise into the potential gradient attention scores of all outgoing edges at each node. This noise is generated by sampling from a uniform distribution, taking the negative logarithm twice, and then taking the negative logarithm again. The attention score, after adding the Gumbel noise, is divided by the temperature parameter `tau`, which is set to 0.5 to control the smoothness of the discretization. The division result is normalized using the Softmax function to obtain a softened discrete distribution of all outgoing edges. For each node's outgoing edge set, the top k edges with the largest attention weights are selected and retained. The value of `k` is set to 3, indicating that each node propagates the fault to a maximum of 3 target nodes. The directional attention weights of unselected edges are reset to 0, while the directional attention weights of selected edges retain their Softmax-normalized values. The directional attention weights are stored in a sparse matrix format, recording only the edge indices and weight values ​​with non-zero weights.

[0090] For each edge in the graph structure, the 128-dimensional initial node embedding of the source node is multiplied element-wise with the directional attention weights of the edge. The result of the multiplication is the weighted message vector passed by the edge. The weighted message vector has the same 128-dimensional dimension as the initial node embedding and is stored in the edge's message cache.

[0091] For each target node, the weighted message vector of all incoming edges is aggregated and fused with the target node's own initial node embedding through a gated recurrent unit to obtain the node embedding representation. The aggregation process traverses all incoming edges of the target node, extracts the weighted message vector of each incoming edge, and performs a summation operation. The summation result is the aggregated message received by the target node. The aggregated message and the target node's own initial node embedding are input into the gated recurrent unit, which contains three components: a reset gate, an update gate, and a candidate hidden state. The reset gate is calculated by concatenating the aggregated message and the initial node embedding, then passing it through a fully connected layer and a sigmoid activation function. The output dimension of the fully connected layer is 128. The update gate is calculated in the same way as the reset gate, also outputting a 128-dimensional vector. The candidate hidden state is calculated by concatenating the aggregated message and the initial node embedding modulated by the reset gate, then passing it through a fully connected layer and a Tanh activation function. The node embedding representation is obtained by weighted fusion of the initial node embedding and the candidate hidden state through the update gate. Specifically, the update gate is multiplied by the candidate hidden state, then 1 is added, and the result is subtracted from the update gate and multiplied by the initial node embedding. The fused node embedding is represented as a 128-dimensional floating-point vector.

[0092] The node embedding represents the input potential energy prediction network, which is a three-layer fully connected network. The first layer has an input dimension of 128 and an output dimension of 64, using LeakyReLU activation with a negative slope of 0.2. The second layer has an input dimension of 64 and an output dimension of 32, also using LeakyReLU activation. The third layer has an input dimension of 32 and an output dimension of 1, using Sigmoid activation. The output is the propagation potential energy value, ranging from 0 to 1. The propagation potential energy value reflects the strength of a node's fault propagation capability after one-way message transmission; a higher value indicates that the node is more likely to become a relay node in subsequent propagation paths.

[0093] This invention proposes a directional message passing mechanism based on potential energy difference. By decomposing the fault feature vector and fusing multi-dimensional information, it achieves accurate modeling of the fault propagation direction. In particular, the construction of the multi-view potential energy difference representation simultaneously considers direct potential energy difference, embedded spatial potential energy difference, and temporal potential energy difference, comprehensively capturing the physical characteristics and temporal dependencies of fault propagation. The introduction of the Gumbel-Softmax method enables the model to achieve differentiable discretization and sparsity, ensuring both the continuity of gradient backpropagation and efficient sparsity in message passing, thus improving computational efficiency.

[0094] In an optional implementation, the step of combining the encoded potential difference representation with the initial propagation potential value of the connected node to calculate the potential gradient attention score of the edge, and then differentiating and discretizing the potential gradient attention score using the Gumbel-Softmax method and performing top-k sparsification to obtain the directional attention weight includes:

[0095] For each edge in the graph structure, the encoded potential difference is used to calculate the basic potential gradient score by bilinear interaction with the initial propagation potential value of the connected node. The potential energy directionality bias term is calculated based on the difference between the initial propagation potential values ​​of the source node and the target node. The potential energy gradient basic score and the potential energy directionality bias term are added together to obtain the potential energy gradient attention score of the edge.

[0096] The temperature decay coefficient is calculated based on the ratio of the current training round to the total training rounds. The initial temperature value is multiplied by the temperature decay coefficient to obtain the Gumbel temperature parameter. Random noise sampled from the Gumbel distribution is added to the potential energy gradient attention score to obtain the perturbed attention score. The perturbed attention score is divided by the Gumbel temperature parameter and normalized by Softmax to obtain a differentiable discretizable attention weight distribution.

[0097] For all incoming edges of each node, calculate the entropy value of the differentiable discretizable attention weight distribution. Based on the comparison result of the entropy value and the preset entropy threshold, determine the sparsity parameter k value. Select the k incoming edges with the largest attention weights and renormalize them to obtain the directional attention weights.

[0098] For example, for each edge in the graph structure, the bilinear interaction process constructs a bilinear transformation matrix with dimensions 16 x 1, where 16 corresponds to the dimension of the encoded potential difference representation and 1 corresponds to the scalar dimension of the initial propagation potential value. The encoded potential difference is represented as a 16-dimensional floating-point vector, and the initial propagation potential value of the source node is a single floating-point scalar. The bilinear interaction obtains an intermediate representation vector by multiplying the encoded potential difference representation with the bilinear transformation matrix. This intermediate representation vector is a 1-dimensional numerical value, which is then multiplied by the initial propagation potential value of the source node to obtain the basic potential gradient score. The weight parameters of the bilinear transformation matrix are set through random initialization using a Xavier uniform distribution, and their value range is adaptively determined based on the input and output dimensions. The basic potential gradient score is a scalar floating-point number with an unrestricted value range, which can be positive, negative, or zero.

[0099] The potential energy directionality bias term is calculated by subtracting the initial propagation potential energy value of the target node from the initial propagation potential energy value of the source node. This difference is then mapped using a nonlinear transformation function. The nonlinear transformation function is the hyperbolic tangent function Tanh, with the potential energy difference as input and the output compressed to the range of -1 to +1. The transformed value is multiplied by a learnable scaling factor to obtain the potential energy directionality bias term. The scaling factor is initialized to 1.0 and ranges from 0.1 to 10.0, adaptively adjusted during training. The potential energy directionality bias term reflects the strength of the potential energy flow directionality between edge-connected nodes. A positive value indicates that the potential energy decreases from the source node to the target node, conforming to the propagation direction; a negative value indicates that the potential energy flows in the opposite direction, violating the propagation logic. The absolute value represents the strength of the directional constraint.

[0100] The potential gradient attention score of an edge is obtained by adding the basic potential energy gradient score to the potential energy directionality bias term. The potential energy gradient attention score comprehensively considers the multi-view features in the encoded potential energy difference representation and the potential energy directionality constraints between nodes; a larger value indicates a higher importance of the edge in the fault propagation path. The potential energy gradient attention score is stored in the edge attribute dictionary, with the key being the edge's unique identifier and the value being the calculated score.

[0101] The current training epoch is obtained through a global counter maintained by the training control module. The total number of training epochs is preset as a hyperparameter, with a default value of 500 epochs. The progress value is calculated by dividing the current training epoch by the total number of training epochs, resulting in a progress value between 0 and 1. The temperature decay coefficient is calculated using an exponential decay strategy, specifically by multiplying the progress value by the decay rate parameter and then taking the negative power of the exponent. The decay rate parameter is set to 5.0 to ensure that the temperature is higher in the early stages of training and decreases rapidly in the later stages. The temperature decay coefficient ranges from 0 to 1, approaching 1 in the early stages of training and approaching 0 in the later stages.

[0102] The initial temperature value is set as a hyperparameter to 1.0, with an adjustable range from 0.1 to 10.0. The Gumbel temperature parameter is obtained by multiplying the initial temperature value by the temperature decay coefficient, and its value changes dynamically with the training progress, approaching the initial temperature value in the early stages of training and approaching 0 in the later stages. The temperature parameter controls the smoothness of the discretization of the Gumbel-Softmax method; the higher the temperature, the closer the output is to a uniform distribution, and the lower the temperature, the closer the output is to a discrete one-hot distribution.

[0103] Gumbel distribution sampling is achieved by generating standard uniformly distributed random numbers and then performing two negative logarithmic transformations. Specifically, a value u is sampled from a uniform distribution of 0 to 1, the negative logarithm is calculated to obtain negative ln(u), and then the negative logarithm of the result is calculated to obtain negative ln(negative ln(u)). This value follows a Gumbel distribution. Gumbel noise is added to the potential gradient attention score of each edge, and the addition operation is a direct summation. The perturbed attention score contains the deterministic information of the original attention score and the random perturbation of Gumbel noise. The random perturbation makes the attention allocation process exploratory during the training phase, avoiding premature convergence to a local optimum.

[0104] The perturbed attention score is divided by the Gumbel temperature parameter and then normalized using Softmax to obtain a differentiable discretizable attention weight distribution. The division operation performs element-wise division between the perturbed attention score and the Gumbel temperature parameter, and the result is input into the Softmax function. Softmax normalization is performed across all outgoing edges of each node. For each outgoing edge, the perturbed attention score is calculated by dividing the result by the temperature parameter. The result is then exponentially applied to all results, and the exponent value of each edge is divided by the sum of the exponent values ​​of all outgoing edges to obtain the normalized attention weight. The normalized attention weights satisfy the probability distribution property: the sum of the weights of all outgoing edges is 1, and the weight of each edge ranges from 0 to 1. The differentiable discretizable attention weight distribution maintains continuous differentiability during the training phase, supporting gradient backpropagation, and gradually approaches a discrete distribution during the inference phase as the temperature parameter approaches 0.

[0105] Entropy calculation iterates through all incoming edges of a node, obtaining the attention weight for each edge. For each weight, the product of its weight and its own logarithm is calculated. The sum of all products is then taken as the negative number to obtain the entropy value. The entropy value reflects the degree of uncertainty in the distribution of attention weights; a larger entropy value indicates a more even distribution of attention weights among the incoming edges, while a smaller entropy value indicates that attention weights are concentrated on a few incoming edges. During entropy calculation, the logarithmic product of incoming edges with a weight of 0 is defined as 0 to avoid numerical calculation errors. The entropy value is stored as a single-precision floating-point number, ranging from 0 to the logarithm of the number of incoming edges.

[0106] The preset entropy threshold is dynamically set based on the number of incoming edges to a node. It is calculated by multiplying the logarithm of the number of incoming edges by a sparsity factor, which is set to 0.5 and ranges from 0.1 to 0.9. The entropy value is compared to the preset entropy threshold. If the entropy value is greater than the threshold, it indicates that the attention weight distribution is too dispersed, and the sparsity parameter k is set to half the number of incoming edges, rounded down. If the entropy value is less than or equal to the threshold, it indicates that the attention weight distribution is relatively concentrated, and the sparsity parameter k is set to one-third of the number of incoming edges, rounded down. The minimum value of the sparsity parameter k is limited to 1, and the maximum value is limited to the number of incoming edges, ensuring that at least one incoming edge is retained and the total number of incoming edges is not exceeded.

[0107] The selection process sorts the attention weights of all incoming edges of a node in descending order and extracts the top k incoming edges as the retained edge set. The attention weights in the retained edge set are then re-normalized. The sum of the attention weights of the retained edges is calculated, and the original attention weight of each retained edge is divided by the sum of the weights to obtain the normalized directional attention weight. The sum of the normalized directional attention weights is 1, and the weight of each retained edge ranges from 0 to 1, reflecting its relative importance in the retained edge set. The directional attention weight of unselected incoming edges is set to 0, indicating that these edges do not participate in message passing in the current propagation path. The directional attention weights are stored in a sparse tensor format, recording only the edge index and weight value for non-zero weights, saving storage space and improving subsequent computational efficiency.

[0108] In one optional implementation, the step of performing edge filtering on the graph structure based on the propagation potential energy value to obtain a directed graph structure, and generating a set of candidate propagation paths by traversing the directed graph structure along the potential energy gradient direction starting from the fault source node includes:

[0109] For each edge in the graph structure, calculate the propagation potential energy difference between the connected nodes, retain the edges with a propagation potential energy difference greater than zero, and assign the edge a direction attribute from the high potential energy node to the low potential energy node according to the potential energy gradient direction, thus forming a directed graph structure.

[0110] Starting from the fault source node as the current traversal node, select the node with the smallest propagation potential value among the outgoing edge neighbors of the current traversal node as the next hop traversal node, and add the directed edge between the current traversal node and the next hop traversal node to the path sequence.

[0111] When the propagation potential energy of the next hop traversed node is a local minimum or the path sequence length reaches the preset depth limit, the traversal stops and the path sequence is recorded as a candidate propagation path. For branch nodes with multiple outgoing neighbor nodes during the traversal, independent traversal is performed on each branch to generate the corresponding candidate propagation path, and the traversal results of all branches are summarized to form a set of candidate propagation paths.

[0112] For example, for each edge in the graph structure, the propagation potential energy difference is calculated by subtracting the propagation potential energy value of the target node from the propagation potential energy value of the source node. The source and target nodes are determined according to the initial defined direction of the edge. The propagation potential energy value is extracted from the output of the potential energy prediction network. Each node corresponds to a scalar floating-point number representing the propagation potential energy value, ranging from 0 to positive infinity. The propagation potential energy difference is a scalar floating-point number, which can be positive, negative, or zero. A positive value indicates that the potential energy decreases from the source node to the target node, a negative value indicates that the potential energy increases in the opposite direction, and zero indicates that the potential energy remains constant. The propagation potential energy difference is stored in the edge attribute dictionary, where the key is the unique identifier of the edge, and the value is the calculated difference value.

[0113] The edge selection process traverses all edges in the graph structure, reads the propagation potential energy difference of each edge, and compares the difference with zero. If the propagation potential energy difference is greater than zero, the edge is retained in the directed graph structure, and a direction attribute is assigned to it, pointing from the node with the higher propagation potential energy value to the node with the lower propagation potential energy value. The direction attribute is implemented by recording the source node identifier and the target node identifier in the edge attributes; the source node identifier corresponds to the high potential energy node, and the target node identifier corresponds to the low potential energy node. If the propagation potential energy difference is less than or equal to zero, the edge is not retained and is removed from the directed graph structure. The directed graph structure is stored using an adjacency list data structure. Each node maintains a list of outgoing edges, which contains key-value pairs of the target node identifier and the edge attributes. The node set of the directed graph structure is the same as that of the original graph structure, and the edge set is a subset of the edges in the original graph structure whose propagation potential energy difference is greater than zero.

[0114] The identifier of the fault source node is obtained through input parameters, and the currently traversed node is initialized with the identifier of the fault source node. The traversal state maintains a stack structure to record traversal path and branch point information; each element of the stack is a tuple of node identifiers and path sequences. The path sequence is initialized as an empty list to store the directed edges traversed during the traversal. The state information of the currently traversed node includes the node identifier, the node's propagation potential value, and the list of its outgoing neighbor nodes, which are retrieved from the directed graph structure. The node with the smallest propagation potential value among the outgoing neighbor nodes of the currently traversed node is selected as the next-hop traversal node. The list of outgoing neighbor nodes is extracted from the adjacency list of the directed graph structure, containing the identifiers of all target nodes whose outgoing edges are from the currently traversed node. The list of outgoing neighbor nodes is traversed, and the propagation potential value of each neighbor node is retrieved from the node attribute dictionary. The propagation potential values ​​of all neighbor nodes are compared, and the node with the smallest propagation potential value is selected as the next-hop traversal node. If multiple neighboring nodes have the same minimum propagation potential value, the node corresponding to the smallest identifier is selected according to the lexicographical order of the node identifiers, ensuring the determinism of the selection process. The identifier and propagation potential value of the next hop traversal node are recorded in the traversal state.

[0115] Directed edges are represented by ordered pairs of source and target node identifiers, where the source node identifier is the identifier of the currently traversed node, and the target node identifier is the identifier of the next hop traversed node. Directed edges are added to the end of the path sequence, which is maintained using a list data structure where each element is an ordered pair of directed edges. The length of the path sequence is counted by the number of list elements; the initial length is zero, and the length increases by one for each added directed edge.

[0116] Local minima are determined by checking if the propagation potential value of the next-hop traversed node is less than or equal to the propagation potential values ​​of all its outgoing neighbor nodes. The list of outgoing neighbor nodes for the next-hop traversed node is extracted from the directed graph structure. If the list is empty, it indicates that the node has no outgoing edges and is considered a local minimum. If the list is not empty, all outgoing neighbor nodes are traversed, and the propagation potential value of each neighbor node is queried. The propagation potential value of the next-hop traversed node is compared with the propagation potential values ​​of all neighbor nodes. If the propagation potential value of the next-hop traversed node is less than or equal to the propagation potential values ​​of all neighbor nodes, it is determined to be a local minimum. A preset depth limit is set as a hyperparameter, with a default value of 10 and a range of 5 to 50. The path sequence length is obtained by counting the number of directed edges in the path sequence and comparing the length with the preset depth limit. If the length reaches or exceeds the preset depth limit, traversal stops. When traversal stops, the current path sequence is copied as a candidate propagation path, and the candidate propagation path is stored in the candidate propagation path set. The candidate propagation path set is maintained using a list data structure, where each element is a copy of the path sequence.

[0117] If the next-hop traversal node still has outgoing neighbor nodes, then the next-hop traversal node is set as the new current traversal node, and the operation of selecting the outgoing neighbor node with the smallest propagation potential value is repeated. The list of outgoing neighbor nodes of the next-hop traversal node is queried from the directed graph structure. If the list is not empty and no local minimum or depth limit stopping conditions are triggered, the next-hop traversal node is updated as the new current traversal node. The identifier and propagation potential value of the current traversal node are updated to the corresponding values ​​of the next-hop traversal node, and the traversal state transitions to the next-hop node. The operation of selecting the node with the smallest propagation potential value from the outgoing neighbor nodes of the current traversal node is repeated, forming a loop traversal process. The loop traversal process maintains a set of visited nodes, which records the identifiers of all nodes traversed during the traversal process to avoid forming loops. When selecting the next-hop traversal node, nodes in the set of visited nodes are excluded, and only the node with the smallest propagation potential value is selected from the unvisited outgoing neighbor nodes. If all outgoing neighbor nodes of the currently traversed node have been visited, the stopping condition is triggered, and the current path sequence is recorded as a candidate propagation path.

[0118] Branch node identification is achieved by checking the number of outgoing neighbor nodes of the currently traversed node. If the number of outgoing neighbor nodes is greater than one, the node is a branch node. Branch node traversal employs a depth-first search strategy, with a traversal stack structure used to manage branch points and traversal states. When a branch node is encountered, the traversal stack records the identifier of the currently traversed node, a copy of the current path sequence, a copy of the currently visited node set, and a list of untraversed outgoing neighbor nodes. The node with the smallest propagation potential value in the outgoing neighbor list is selected as the next hop traversal node for the current branch, and depth-first traversal of that branch continues until a stopping condition is triggered. After the stopping condition is triggered, the branch point state is popped from the traversal stack, restoring the branch node's traversal state. The next node with the second smallest propagation potential value in the untraversed outgoing neighbor list is selected, and traversal of the new branch begins. Each branch's traversal process independently maintains its path sequence and visited node set, ensuring state isolation between different branches. After all branches have been traversed, the traversal stack is empty, ending the overall traversal process.

[0119] The candidate propagation path set accumulates during the traversal. Each time a branch traversal triggers the stopping condition, the path sequence is added to the candidate propagation path set. The candidate propagation path set is stored using a list data structure, where each list element is a path sequence object, and each path sequence object contains an ordered list of directed edges. The candidate propagation path set supports deduplication; two paths are determined to be duplicates by comparing the ordered sets of directed edges in the path sequences, and duplicate paths are retained only once. The length of the candidate propagation path set reflects the number of different propagation paths generated during the traversal, ranging from one to an exponential order of the number of branch nodes. The candidate propagation path set is stored in memory or persistent storage for use by subsequent path scoring and filtering modules.

[0120] This invention transforms an undirected power grid topology into a directed graph structure reflecting the direction of fault propagation by filtering based on potential energy differences and defining directions. The traversal strategy along the potential energy gradient simulates the physical propagation process of the fault in the power grid, conforming to the principle of energy minimization. The branching mechanism can explore multiple possible propagation paths, enhancing the algorithm's exploration breadth, while the preset depth limit effectively controls the search space, balancing identification accuracy and computational efficiency, making it suitable for complex power grid topologies.

[0121] In one optional implementation, the step of inputting the subgraph structure of each path in the candidate propagation path set and the node embedding representation into the invariant verification graph neural network, and outputting a score on the degree of satisfaction of the physical invariant constraint set, wherein the physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints, includes:

[0122] For each path in the candidate propagation path set, extract the subgraph structure formed by all nodes and edges in the path, and concatenate the node embedding representation of each node in the subgraph structure with the position code of the node in the path to obtain the position-aware node representation. Calculate the directional edge weight for the edges in the subgraph structure based on the difference between the initial propagation potential energy values ​​of the connected nodes and the electrical impedance parameters of the edges.

[0123] The location-aware node representation and the directional edge weights are input into the invariant verification graph neural network. The updated node representation is obtained by message passing through the graph convolutional layer. The updated node representation is then pooled along the path sequence direction to obtain the path-level representation vector.

[0124] Based on the path-level representation vector, the topological connectivity score, power conservation score, and temporal causality score are calculated respectively. The temporal causality score is calculated by verifying the temporal increasing relationship of the node failure time in the path. The topological connectivity score, power conservation score, and temporal causality score are weighted and fused to obtain the degree of satisfaction of the physical invariant constraint set.

[0125] Combination Figure 2 The flowchart illustrates the multi-dimensional constraint verification and scoring mechanism for physical invariants of power grid fault propagation based on graph neural networks. For example, candidate propagation paths are stored as an ordered list of directed edges, each represented by an ordered pair of source and target node identifiers. The directed edge list of candidate propagation paths is traversed, collecting all occurrences of node identifiers to form a subgraph node set. This node set is maintained using a set data structure for automatic deduplication. The subgraph edge set is directly constructed from the directed edge list of candidate propagation paths, with each edge associated with both a source and target node identifier. The subgraph structure is stored using an adjacency list data structure, where the key is the node identifier and the value is a list of edges connected to that node. The edge list contains the target node identifier and an edge attribute dictionary. The number of nodes in the subgraph structure is the path length plus one, and the number of edges equals the path length, which is counted by the number of elements in the directed edge list. The subgraph structure extracts attribute data for corresponding nodes and edges from the complete graph structure. Node attributes include fields such as node embedding representation, propagation potential value, and fault time; edge attributes include fields such as electrical impedance parameters and the difference between initial propagation potential values.

[0126] Node embeddings are extracted from the node attribute dictionary and are fixed-dimensional floating-point vectors with a default dimension of 128, ranging from -1 to +1. The position encoding of a node within a path is obtained by calculating its index in the directed edge list, starting from 0. The starting node's index is 0, and the ending node's index equals the path length. Position encoding uses a sine-cosine coding method, applying sine and cosine functions of different frequencies to the position index values ​​to generate a fixed-dimensional encoding vector. The dimension of the position encoding vector is the same as the node embedding dimension, both being 128-dimensional. The frequency parameters of the sine-cosine coding are calculated by dividing the position index by a power of 10000, with the exponent increasing from 0 to half the encoding dimension. The position encoding vector and the node embedding vector are concatenated by appending elements of the position encoding vector to the end of the node embedding vector, forming a 256-dimensional position-aware node representation vector. Location-aware nodes are stored in a node attribute dictionary, with the key being the node identifier and the value being a 256-dimensional floating-point vector.

[0127] The initial propagation potential energy difference is extracted from the edge attribute dictionary. This difference is calculated during the edge selection stage and represents the result of subtracting the propagation potential energy of the target node from the source node's propagation potential energy. The electrical impedance parameter is also extracted from the edge attribute dictionary. This parameter is a scalar floating-point number representing the resistance or impedance of the physical line corresponding to the edge, ranging from 0.01 to 100 ohms. The directional edge weight is calculated by dividing the initial propagation potential energy difference by the electrical impedance parameter. Before the division operation, the electrical impedance parameter is truncated to a lower bound of 0.01 to avoid division by zero. The directional edge weight is a floating-point number ranging from negative infinity to positive infinity. A positive value indicates that the potential energy decreases along the edge direction and the impedance is small, while a negative value indicates that the potential energy is reversed or the impedance is extremely large. The directional edge weights are stored in the edge attribute dictionary, with the key being the edge's unique identifier and the value being the calculated weight value. After the directional edge weights of all edges in the subgraph structure have been calculated, the edge attribute dictionary contains complete weight information for subsequent use by the graph neural network.

[0128] The invariant verification graph neural network adopts a graph convolutional network architecture. The network input layer receives position-aware node representations and directional edge weights. The input interface is defined as a node feature matrix and an edge weight matrix. The number of rows in the node feature matrix equals the number of nodes in the subgraph, and the number of columns is 256. Each row corresponds to a position-aware node representation vector of a node. The edge weight matrix is ​​a sparse matrix, with row and column indices corresponding to node identifiers. The matrix elements are directional edge weights, and elements at positions where no edge exists have a value of 0. The node feature matrix and edge weight matrix undergo standardization preprocessing. For each column of the node feature matrix, the column mean is subtracted and then divided by the column standard deviation. For each row of the edge weight matrix, the row's L2 norm is divided. Rows with a norm less than 0.01 retain their original values ​​to avoid numerical instability.

[0129] The graph convolutional layer contains three convolutional modules, each performing neighborhood aggregation and feature transformation operations. In the neighborhood aggregation stage of the convolutional module, the position-aware node representations of the incoming neighbor nodes are collected for each node. These neighbor node representations are then weighted and summed with the directional edge weights of their corresponding edges. The weights are the normalized values ​​of the directional edge weights. The normalization process calculates the sum of the absolute values ​​of the directional edge weights of all incoming edges of a node. The weight of each edge is divided by this sum to obtain the normalized weight; if the sum is 0, the normalized weight is set to 0. The weighted sum is the aggregated feature vector, with the same 256 dimensions as the position-aware node representation. In the feature transformation stage of the convolutional module, the node's own position-aware node representation is concatenated with the aggregated feature vector to form a 512-dimensional vector. This vector is then input into the fully connected layer for linear transformation. The fully connected layer's weight matrix has a dimension of 512 x 256, the bias vector has a dimension of 256, and the activation function is a modified linear unit function. The three convolutional layers are executed sequentially. The first layer has an output dimension of 256, the second layer has an output dimension of 256, and the third layer has an output dimension of 128. Residual connections and layer normalization operations are inserted between each convolutional layer. The residual connections add the input features to the output features, and the layer normalization subtracts the mean from each dimension of the feature vector and divides it by the standard deviation. The 128-dimensional vector output by the third convolutional layer is used as the updated node representation. The updated node representations of all nodes form the updated node feature matrix, with the number of rows equal to the number of nodes in the subgraph and 128 columns.

[0130] The path sequence direction is defined by a list of directed edges for candidate propagation paths. The order of the directed edges indicates the traversal order of the path from the start point to the end point. Pooling operations extract the updated node representation of each node in the path, arranging them according to their index in the path to form a sequence matrix. The number of rows in the sequence matrix equals the path length plus one, and the number of columns is 128. Average pooling is applied to the sequence matrix, calculating the arithmetic mean of all row elements in each column to obtain a 128-dimensional path-level representation vector. This path-level representation vector reflects the global feature information of the entire path, and the range of values ​​for each vector element is determined by the range of values ​​for the updated node representations. The path-level representation vector is stored in a path attribute dictionary, with the key being the path identifier and the value being a 128-dimensional floating-point vector.

[0131] The topological connectivity score is calculated by evaluating the continuity of nodes and edges in the path. The topological connectivity score module receives the path-level representation vector and subgraph structure, inputs them into a fully connected layer for feature mapping, and the weight matrix of the fully connected layer has a dimension of 128 x 64, outputting a 64-dimensional hidden vector. The hidden vector is input into the second fully connected layer, where the weight matrix has a dimension of 64 x 1, and the output scalar value is mapped to the interval of 0 to 1 using the Sigmoid function. The mapped value is used as the topological connectivity score. A topological connectivity score close to 1 indicates strong path topological connectivity, while a score close to 0 indicates weak connectivity.

[0132] The power conservation score is calculated by verifying the conservation relationship of power flow in the path. The power conservation score module receives the path-level representation vector and the power injection and outflow data of the nodes in the path. Node power data is extracted from the node attribute dictionary, including the active and reactive power of the node, in kilowatts and kilovars respectively. The path-level representation vector is input to the fully connected layer and mapped to a 64-dimensional hidden vector. The hidden vector is concatenated with the statistical features of the power data of all nodes in the path, including the total power, power variance, maximum power, and minimum power. The concatenated vector is input to the second fully connected layer, with a weight matrix of 68 x 1. The output scalar value is mapped to the interval between 0 and 1 using the Sigmoid function. The mapped value is used as the power conservation score. A power conservation score close to 1 indicates good power conservation of the path, while a score close to 0 indicates poor conservation.

[0133] The temporal causality score is calculated by verifying the increasing temporal relationship of node failure times in the path. Node failure times are extracted from the node attribute dictionary, presented as timestamps representing the absolute time or the number of seconds relative to the start time of the failure, with millisecond precision. The calculation process iterates through adjacent node pairs in the path, extracting the failure time for each pair. The time difference is obtained by subtracting the failure time of the preceding node from the current node's failure time. A positive time difference indicates an increasing temporal relationship consistent with causality, while a negative or zero time difference indicates a temporal anomaly. The number of pairs of adjacent node pairs with positive time differences is counted and divided by the total number of adjacent node pairs to obtain the temporal causality compliance rate. The temporal causality compliance rate is a floating-point number between 0 and 1. This compliance rate is concatenated with the path-level representation vector and input into a fully connected layer. The weight matrix of the fully connected layer has a dimension of 129 x 64, outputting a 64-dimensional hidden vector. The hidden vector is input into the second fully connected layer. The weight matrix has a dimension of 64 x 1. The output scalar value is mapped to the interval between 0 and 1 using the Sigmoid function. The mapped value is used as the temporal causality score. A temporal causality score close to 1 indicates strong temporal causality of the path, while a score close to 0 indicates weak causality.

[0134] The three scores are multiplied by their respective weighting coefficients: the topological connectivity score has a default weighting coefficient of 0.3, the power conservation score has a default weighting coefficient of 0.4, and the temporal causality score has a default weighting coefficient of 0.3. The sum of these three weighting coefficients is 1. The three weighted scores are then summed to obtain the satisfaction score, which is a floating-point number between 0 and 1. The weighting coefficients can be adjusted according to the application scenario, ranging from 0 to 1 with a precision of 0.01. The satisfaction score is stored in the path attribute dictionary, with the path identifier as the key and the score value as the path name. The satisfaction score is used for subsequent path sorting and filtering; a higher score indicates that the path better conforms to the physical invariant constraints.

[0135] The training process of the invariant validation graph neural network adopts a supervised learning approach. The training data includes candidate propagation paths and manually labeled physical invariant satisfaction labels. The training dataset is collected from historical failure cases, with each case containing real failure propagation paths and spurious interference paths. The satisfaction label for real paths is 1, and the satisfaction label for spurious paths is 0. The training dataset is divided into training, validation, and test sets in an 8:1:1 ratio. The training process uses the cross-entropy loss function, which calculates the difference between the predicted satisfaction score and the actual satisfaction label. The optimizer is an adaptive moment estimation optimizer with an initial learning rate of 0.001 and a batch size of 32. The training iterations consist of 100 epochs, traversing all samples in the training set in each epoch. The validation set is used to monitor overfitting, and the test set is used to evaluate the final model performance. Model parameters are updated through gradient backpropagation, with a gradient clipping threshold of 1.0 to avoid gradient explosion. The trained model parameters are saved to persistent storage. During the inference phase, the model parameters are loaded to calculate the satisfaction score for new candidate propagation paths.

[0136] This invention transforms abstract physical laws into computable scoring metrics. The design of position-aware node representations and directional edge weights enables the model to accurately capture the structural information and electrical characteristics of the path. Parallel verification of three types of physical constraints ensures the consistency of the identification results with the physical laws of the power system, significantly improving the reliability and interpretability of the results. This helps power system experts understand and verify the algorithm output and promotes human-machine collaborative analysis.

[0137] In one alternative implementation, the step of calculating the satisfaction score of the set of physical invariant constraints includes:

[0138] The path-level representation vector is used to predict the connection probability between each pair of adjacent nodes in the path. The connection probability is compared with the actual connection relationship of the corresponding node pair in the power grid topology adjacency matrix and aggregated to obtain the topology connectivity score.

[0139] The path-level representation vector is mapped to the predicted power inflow and outflow vectors of each node in the path. Based on the predicted power inflow and outflow vectors and the line impedance parameters of the edges between adjacent nodes in the path, the cumulative power loss value along the path is calculated. The cumulative power loss value is compared with the initial propagation potential energy value of the starting node of the path to obtain the power conservation score.

[0140] The topological connectivity score is used as a gating signal to filter the power conservation score to obtain the topologically constrained power score. The temporal causality score is used as a directional weight to weight and correct the topologically constrained power score to obtain the temporally corrected power score. The degree of satisfaction of the physical invariant constraint set is obtained by weighted fusion of the temporally corrected power score, the topological connectivity score and the temporal causality score.

[0141] For example, when calculating the satisfaction score of the set of physical invariant constraints, the path-level representation vector is first input into the connection probability prediction module. The connection probability prediction module employs a bilinear mapping structure, receiving a 128-dimensional path-level representation vector as input. It maps the path-level representation vector to a 64-dimensional intermediate feature vector through a fully connected layer. The weight matrix of the fully connected layer has a dimension of 128 x 64, the bias vector has a dimension of 64, and the activation function is a modified linear unit function. The intermediate feature vector is then input into the node pair embedding generation layer, which generates an independent node pair embedding representation for each pair of adjacent nodes in the path. During the node pair embedding generation process, the corresponding row vector in the updated node feature matrix is ​​extracted for each pair of adjacent nodes in the path. The 128-dimensional vectors of the two nodes are concatenated to form a 256-dimensional original node pair feature. The original features of each node pair are concatenated with the 64-dimensional intermediate feature vector to obtain a 320-dimensional joint feature vector. This joint feature vector is then input into a bilinear layer for interactive modeling. The bilinear layer contains two weight matrices: the first weight matrix has a dimension of 320 x 32, and the second weight matrix has a dimension of 32 x 1. The product of these two matrices outputs a scalar value. This scalar value is mapped to the interval 0 to 1 using the sigmoid function. The mapped value serves as the connection probability of the node pair. A connection probability close to 1 indicates a high confidence that there is a physical connection between the two nodes, while a probability close to 0 indicates a low confidence. After calculating the connection probabilities of all adjacent node pairs in the path, a connection probability vector is formed. The length of this vector is equal to the path length, and its elements are the connection probability values ​​of each node pair.

[0142] The power grid topology adjacency matrix is ​​extracted from the power grid topology database. The adjacency matrix is ​​a two-dimensional square matrix, where the row and column indices correspond to the identifiers of all nodes in the power grid. Matrix elements take values ​​of 0 or 1; a value of 1 indicates a physical connection between nodes in the corresponding row and column, while a value of 0 indicates no connection. The adjacency matrix is ​​stored in a sparse matrix format, using a coordinate list format, recording only the row and column indices of elements with a value of 1 to reduce storage space. For each pair of adjacent nodes in the path, the actual connection relationship of the corresponding node pair is queried from the adjacency matrix. The query operation accesses the corresponding element in the adjacency matrix using the node identifier as the index; an element value of 1 indicates an actual connection, and an element value of 0 indicates no connection. The predicted connection probability is compared with the actual connection relationship, and a consistency metric is calculated during the comparison process. When the actual connection relationship is 1, the consistency metric equals the connection probability itself; when the actual connection relationship is 0, the consistency metric equals 1 minus the connection probability. The topology connectivity score is obtained by arithmetically averaging the consistency metrics of all node pairs in the path. The topology connectivity score is a floating-point number between 0 and 1, with a value closer to 1 indicating a higher consistency between the path's topology connectivity and the actual power grid topology. The topology connectivity score is stored in the path attribute dictionary, with the path identifier as the key and the score value as the value.

[0143] The path-level representation vector is then input into the power flow prediction module, which employs a sequence-to-sequence mapping structure. The path-level representation vector is first expanded into an unfolded matrix by multiplying the number of path nodes by the power flow dimension (default 2), corresponding to active and reactive power respectively. The number of rows in the unfolded matrix equals the path length plus one, and the number of columns is two. The dimension of the fully connected layer weight matrix is ​​128 multiplied by the number of nodes multiplied by 2. Each row of the unfolded matrix corresponds to the predicted power inflow / outflow vector for a node in the path. The first element of the vector represents the difference between the predicted active power inflow and outflow values ​​for that node, and the second element represents the difference between the predicted reactive power inflow and outflow values, in kilowatts and kilovars respectively, ranging from -1000 to +1000. The predicted power inflow / outflow vector undergoes a nonlinear transformation using a residual activation function. The residual activation function adds the input value to the output value of the modified linear unit function and divides by 2 to ensure the smoothness of the output value.

[0144] The cumulative power loss along the path is calculated based on the predicted power inflow and outflow vectors, traversing all adjacent node pairs in the path. For each pair of adjacent nodes, the line impedance parameters connecting the two nodes are extracted from the edge attribute dictionary. These parameters include resistance and reactance values; resistance is in ohms and ranges from 0.01 to 100, while reactance is in ohms and ranges from 0.01 to 50. The line loss power of an edge is calculated using the difference between the predicted power inflow and outflow vectors of adjacent nodes and the line impedance parameters. Active power loss equals the square of the difference in active power between the two nodes multiplied by the resistance value, and reactive power loss equals the square of the difference in reactive power between the two nodes multiplied by the reactance value. The cumulative active power loss is obtained by summing the active power losses of all edges in the path, and the cumulative reactive power loss is obtained by summing the reactive power losses of all edges. The cumulative power loss is a vector combination of the cumulative active power loss and cumulative reactive power loss values, with a vector dimension of 2 and units of kilowatts and kilovars. The validity of the cumulative power loss value is verified by comparing the residual with the initial propagation potential energy value of the path starting node. The active power loss ratio is calculated by dividing the cumulative active power loss value by the initial propagation potential energy value of the path starting node, and the reactive power loss ratio is calculated by dividing the cumulative reactive power loss value by the initial propagation potential energy value of the path starting node. The reasonable range for the loss ratio is 0.05 to 0.5; values ​​outside this range indicate abnormal power conservation. The power conservation score is calculated by mapping the loss ratio to the interval between 0 and 1. The mapping function uses a piecewise linear function. When the loss ratio is between 0.05 and 0.3, the score decreases linearly from 1 to 0.5; when the loss ratio is between 0.3 and 0.5, the score decreases linearly from 0.5 to 0. A score of 1 is given when the loss ratio is less than 0.05, and a score of 0 is given when the loss ratio is greater than 0.5. The active power conservation score and the reactive power conservation score are weighted and averaged, with an active power weight of 0.6 and a reactive power weight of 0.4. The weighted average result is used as the final power conservation score, which is stored in the path attribute dictionary.

[0145] When using topology connectivity score as a gating signal to filter power conservation score, a soft gating mechanism is employed. This mechanism converts the topology connectivity score into a gating strength coefficient using a threshold function, which is a smooth step function. The gating strength coefficient is 0 when the topology connectivity score is less than 0.5 and 1 when it is greater than 0.7. Within the range of 0.5 to 0.7, the gating strength coefficient increases linearly from 0 to 1. The power conservation score is multiplied by the gating strength coefficient to obtain the topology-constrained power score. The topology-constrained power score ranges from 0 to the value of the original power conservation score. When the topology connectivity score is low, the topology-constrained power score is strongly suppressed and approaches 0; when the topology connectivity score is high, the topology-constrained power score approaches the original power conservation score. The topology-constrained power score reflects the degree of power conservation while satisfying the topology connectivity constraints.

[0146] The temporal causality score is used as a directional weight to weight and correct the power score after topological constraints. The weighting correction operation employs an adaptive weight adjustment strategy. The temporal causality score generates directional weight coefficients through a nonlinear transformation. The nonlinear transformation function is an exponential weighting function. Subtracting 0.5 from the temporal causality score, multiplying by 2, and then inputting it into the exponential function, the output value of the exponential function is divided by the output value of the exponential function plus 1 to obtain the directional weight coefficient. The value of the directional weight coefficient ranges from 0 to 1. When the temporal causality score is 0.5, the directional weight coefficient is 0.5; when the temporal causality score is greater than 0.5, the directional weight coefficient is greater than 0.5; and when the temporal causality score is less than 0.5, the directional weight coefficient is less than 0.5. Multiplying the power score after topological constraints by the directional weight coefficient yields the temporally corrected power score. The value of the temporally corrected power score is modulated by the temporal causality score; a strong temporal causality enhances the score, while a weak temporal causality weakens the score. The temporally corrected power score serves as a core indicator for verifying integrated physical constraints.

[0147] The satisfaction score of the physical invariant constraint set is obtained by weighted fusion of the time-corrected power score, topological connectivity score, and time-series causality score. The weighted fusion operation adopts a three-component linear combination method. The default weight coefficient of the time-corrected power score is 0.5, the default weight coefficient of the topological connectivity score is 0.3, and the default weight coefficient of the time-series causality score is 0.2. The sum of the three weight coefficients is 1.0, and the precision of the weight coefficients is 0.01. The adjustment range is 0.1 to 0.6 for each item, which can be adjusted according to the application scenario to ensure the weight and constraint. The satisfaction score is obtained by multiplying the three scores by their corresponding weight coefficients and then summing them. The satisfaction score is a floating-point number between 0 and 1 with a numerical precision of 0.001. A satisfaction score greater than 0.8 indicates that the path highly conforms to the physical invariant constraints, 0.5 to 0.8 indicates moderate conformity, and less than 0.5 indicates non-conformity or weak conformity. The satisfaction score is stored in the path attribute dictionary, with the key being the path identifier and the value being the score value. The score is used for subsequent path sorting and filtering, and the path with the higher score is given priority as the inference result of the fault propagation path.

[0148] This invention designs a multi-level fusion mechanism for constraint scoring. It uses topological connectivity as a gating signal to filter power conservation scores, and then uses time-series causality scores for directional weighting correction, achieving hierarchical interaction between constraints. This design not only considers the independent influence of various constraints but also models their interdependencies, forming a comprehensive evaluation system that better reflects the physical characteristics of power systems. The optimal path selection strategy based on topological entropy follows Occam's razor, selecting the simplest interpretation while satisfying physical constraints, thus improving the reliability and interpretability of the identification results.

[0149] A second aspect of the present invention provides a power grid fault propagation path identification system based on graph neural networks, comprising:

[0150] The first unit is used to acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node;

[0151] The second unit is used to construct a direction-aware graph neural network. The fault feature vector is input into the direction-aware graph neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential value and node embedding representation of each node.

[0152] The third unit is used to perform edge filtering on the graph structure based on the propagation potential energy value to obtain a directed graph structure, and to generate a set of candidate propagation paths by traversing the directed graph structure along the potential energy gradient direction starting from the fault source node;

[0153] The fourth unit is used to construct an invariant verification graph neural network. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network, and the output is a score of the degree of satisfaction of the physical invariant constraint set. The physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints.

[0154] The fifth unit is used to filter paths whose satisfaction rating is higher than a preset rating threshold, and select the path with the minimum topological entropy as the fault propagation path.

Claims

1. A method for identifying power grid fault propagation paths based on graph neural networks, characterized in that, include: Acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node; A direction-aware graph neural network is constructed. The fault feature vector is input into the neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential energy value and node embedding representation of each node. This includes: decomposing the fault feature vector into state components and time components, encoding them separately, and then fusing them through cross-attention to obtain the initial node embedding; calculating the initial propagation potential energy value of the node based on the fault time; for each edge in the graph structure, calculating the direct potential energy difference, embedding space potential energy difference, and temporal potential energy difference based on the initial propagation potential energy value of the connected node, the initial node embedding, and the fault time component, respectively, concatenating them into a multi-view potential energy difference feature, and then... A nonlinear transformation yields the encoded potential difference representation. This encoded potential difference representation is then combined with the initial propagation potential energy value of the connected nodes to calculate the edge's potential gradient attention score. The potential gradient attention score is then discretized and top-k sparsified using the Gumbel-Softmax method to obtain directional attention weights. These directional attention weights are applied to the initial node embedding of the source node to generate a weighted message vector. For each target node, the weighted message vectors of all incoming edges are aggregated and fused with the target node's own initial node embedding using a gated recurrent unit to obtain the node embedding representation. Based on this node embedding representation, the propagation potential energy value is calculated using a potential energy prediction network. Based on the propagation potential energy value, the graph structure is edge-filtered to obtain a directed graph structure. Starting from the fault source node, the graph structure is traversed along the potential energy gradient direction to generate a set of candidate propagation paths. This includes: for each edge in the graph structure, calculating the propagation potential energy difference between connecting nodes, retaining edges with a propagation potential energy difference greater than zero, and assigning a direction attribute from a high potential energy node to a low potential energy node according to the potential energy gradient direction, thus forming a directed graph structure; starting from the fault source node as the current traversal node, selecting the node with the smallest propagation potential energy value among the outgoing edge neighbors of the current traversal node as the next hop traversal node, and adding the directed edge between the current traversal node and the next hop traversal node to the path sequence; stopping the traversal when the propagation potential energy value of the next hop traversal node is a local minimum or the path sequence length reaches a preset depth limit, and recording the path sequence as a candidate propagation path; for branch nodes with multiple outgoing edge neighbors during the traversal process, each branch is independently traversed to generate a corresponding candidate propagation path, and all branch traversal results are summarized to form a set of candidate propagation paths. An invariant verification graph neural network is constructed. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network. The output is a score of the degree of satisfaction of the physical invariant constraint set, which includes topological connectivity constraints, power conservation constraints, and temporal causality constraints. Paths with satisfaction scores higher than a preset score threshold are selected, and the path with the minimum topological entropy is chosen as the fault propagation path.

2. The method according to claim 1, characterized in that... The steps of assigning a fault feature vector containing the fault state and fault time to a node based on the fault observation data and determining the fault source node include: Extract the protection device action information and circuit breaker tripping information of each node from the fault observation data, encode the protection device action information as a fault state characterizing the node, and mark the time corresponding to the circuit breaker tripping information as the fault time of the node; combine the fault state and the fault time to form a fault feature vector and assign it to the corresponding node. Based on the fault time, a time-series sorting sequence is constructed, and the node with the earliest fault time in the time-series sorting sequence is identified as the fault source node.

3. The method according to claim 1, characterized in that, The steps of combining the encoded potential difference representation with the initial propagation potential value of the connected node to calculate the potential gradient attention score of the edge, and then differentiating and discretizing the potential gradient attention score using the Gumbel-Softmax method and performing top-k sparsification to obtain the directional attention weight include: For each edge in the graph structure, the encoded potential difference is used to calculate the basic potential gradient score by bilinear interaction with the initial propagation potential value of the connected node. The potential energy directionality bias term is calculated based on the difference between the initial propagation potential values ​​of the source node and the target node. The potential energy gradient basic score and the potential energy directionality bias term are added together to obtain the potential energy gradient attention score of the edge. The temperature decay coefficient is calculated based on the ratio of the current training round to the total training rounds. The initial temperature value is multiplied by the temperature decay coefficient to obtain the Gumbel temperature parameter. Random noise sampled from the Gumbel distribution is added to the potential energy gradient attention score to obtain the perturbed attention score. The perturbed attention score is divided by the Gumbel temperature parameter and normalized by Softmax to obtain a differentiable discretizable attention weight distribution. For all incoming edges of each node, calculate the entropy value of the differentiable discretizable attention weight distribution. Based on the comparison result of the entropy value and the preset entropy threshold, determine the sparsity parameter k value. Select the k incoming edges with the largest attention weights and renormalize them to obtain the directional attention weights.

4. The method according to claim 1, characterized in that, The steps of inputting the subgraph structure of each path in the candidate propagation path set and the node embedding representation into the invariant verification graph neural network, and outputting a score on the degree of satisfaction of the physical invariant constraint set, wherein the physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints, include: For each path in the candidate propagation path set, extract the subgraph structure formed by all nodes and edges in the path, and concatenate the node embedding representation of each node in the subgraph structure with the position code of the node in the path to obtain the position-aware node representation. Calculate the directional edge weight for the edges in the subgraph structure based on the difference between the initial propagation potential energy values ​​of the connected nodes and the electrical impedance parameters of the edges. The location-aware node representation and the directional edge weights are input into the invariant verification graph neural network. The updated node representation is obtained by message passing through the graph convolutional layer. The updated node representation is then pooled along the path sequence direction to obtain the path-level representation vector. Based on the path-level representation vector, the topological connectivity score, power conservation score, and temporal causality score are calculated respectively. The temporal causality score is calculated by verifying the temporal increasing relationship of the node failure time in the path. The topological connectivity score, power conservation score, and temporal causality score are weighted and fused to obtain the degree of satisfaction of the physical invariant constraint set.

5. The method according to claim 4, characterized in that, The steps for calculating the satisfaction score of the set of physical invariant constraints include: The path-level representation vector is used to predict the connection probability between each pair of adjacent nodes in the path. The connection probability is compared with the actual connection relationship of the corresponding node pair in the power grid topology adjacency matrix and aggregated to obtain the topology connectivity score. The path-level representation vector is mapped to the predicted power inflow and outflow vectors of each node in the path. Based on the predicted power inflow and outflow vectors and the line impedance parameters of the edges between adjacent nodes in the path, the cumulative power loss value along the path is calculated. The cumulative power loss value is compared with the initial propagation potential energy value of the starting node of the path to obtain the power conservation score. The topological connectivity score is used as a gating signal to filter the power conservation score to obtain the topologically constrained power score. The temporal causality score is used as a directional weight to weight and correct the topologically constrained power score to obtain the temporally corrected power score. The degree of satisfaction of the physical invariant constraint set is obtained by weighted fusion of the temporally corrected power score, the topological connectivity score and the temporal causality score.

6. A power grid fault propagation path identification system based on graph neural networks, used to implement the method of any one of claims 1-5, characterized in that, include: The first unit is used to acquire power grid topology data and fault observation data, represent the power grid topology data as a graph structure, assign fault feature vectors containing fault state and fault time to nodes based on the fault observation data, and determine the fault source node; The second unit is used to construct a direction-aware graph neural network. The fault feature vector is input into the direction-aware graph neural network, and unidirectional message passing is performed through directional attention weights based on potential energy difference to obtain the propagation potential value and node embedding representation of each node. The third unit is used to perform edge filtering on the graph structure based on the propagation potential energy value to obtain a directed graph structure, and to generate a set of candidate propagation paths by traversing the directed graph structure along the potential energy gradient direction starting from the fault source node; The fourth unit is used to construct an invariant verification graph neural network. The subgraph structure of each path in the candidate propagation path set and the node embedding representation are input into the invariant verification graph neural network, and the output is a score of the degree of satisfaction of the physical invariant constraint set. The physical invariant constraint set includes topological connectivity constraints, power conservation constraints, and temporal causality constraints. The fifth unit is used to filter paths whose satisfaction rating is higher than a preset rating threshold, and select the path with the minimum topological entropy as the fault propagation path.

Images