Power grid equation graph processor parallel back substitution method based on reverse order path tree

Abstract

The application belongs to the technical field of power grid equation graphic processor parallel back substitution, and discloses a power grid equation graphic processor parallel back substitution method based on an inverse order path tree, which comprises a data compression storage step, wherein the index of an elder node set is differentially encoded and compressed and stored by utilizing the local cluster characteristics of a power grid topology; when back substitution calculation is performed on a graphic processor, the sparse structure pointer of a Jacobian matrix is locked through a unified memory addressing mechanism. In the back substitution calculation stage, a mixed precision adaptive back substitution algorithm is adopted, double precision floating point numbers are used for operation in a high level stage, and automatic switching to single precision floating point numbers or half precision floating point numbers is performed for operation in a low level stage, the number of thread of a parallel back substitution kernel function is equal to the number of the elder node set of the current level, and each thread independently takes charge of the calculation task of all nodes in the elder node set. The application can effectively reduce the memory occupation and improve the calculation efficiency.

Description

Technical Field

[0001] This invention belongs to the field of parallel back-substitution technology of power grid equation graphics processors, and relates to a parallel back-substitution method for power grid equation graphics processors based on reverse path trees. Background Technology

[0002] In large-scale power system simulation, power flow calculation and state estimation are core components, both essentially involving the solution of large-scale sparse linear equation systems. As the power grid scales continuously expand, the dimensionality of the equation systems increases dramatically, placing extremely high demands on computational efficiency. Graphical processors (GPUs), with their high parallelism, are often used to accelerate such computational tasks. However, current technologies still face many challenges when using GPUs to handle the back-substitution process of power grid equations. Traditional back-substitution algorithms typically traverse based on node physical numbers or simple hierarchical order. This approach fails to fully consider the complex dependencies between nodes in the power grid topology, leading to insufficient thread utilization during parallel computation. Especially in the construction of the reverse path tree, the lack of an effective hierarchical partitioning strategy can easily result in uneven distribution of computational load among the GPU's computing units, thus limiting the overall speedup.

[0003] Furthermore, existing technologies often overlook specific structural features present in power grid topologies. For example, the situation where multiple parent nodes point to the same child node is quite common in power grid structures; these nodes can be considered a set of elder nodes. Failure to implement dedicated thread allocation and memory management for such sets of elder nodes leads to a large number of redundant memory accesses and atomic operations, severely slowing down computation. Meanwhile, while high-precision double-precision floating-point arithmetic ensures accurate results, it consumes significant computational resources. Existing solutions lack flexible precision adjustment mechanisms and cannot dynamically adjust computational precision based on computational level or residual feedback, making it difficult to achieve a balance between computational speed and numerical accuracy.

[0004] Therefore, there is an urgent need for a parallel back-substitution method that can combine the characteristics of reverse path trees, optimize the processing of elder node sets, and support mixed-precision computation. Summary of the Invention

[0005] To address the problems existing in the background technology, this invention proposes a parallel back-substitution method for power grid equations using a graphical processor based on a reverse path tree.

[0006] The first aspect of the present invention provides a parallel back-substitution method for power grid equations using a graphical processor based on a reverse path tree, comprising:

[0007] Obtain the network topology and Jacobian matrix of the power system, and establish a dependency graph between nodes based on the sparse structure of the Jacobian matrix;

[0008] A reverse hierarchical strategy is used to divide all nodes in the dependency graph into levels. The root node with an in-degree of zero is divided into level zero, and the level of the remaining nodes is determined by adding one to the maximum value of the levels of all its parent nodes.

[0009] Based on the topological characteristic that multiple nodes in the same level depend on the same child node, nodes with the same child node are divided into a set of elder nodes.

[0010] The hierarchical information and the set of elder nodes are mapped to the graphics processor, and the parallel back-generation kernel function is started.

[0011] The graphics processor performs parallel back-substitution calculations on the linear equation system in ascending order of hierarchy, using the set of elder nodes as the thread allocation unit, until all levels are traversed to obtain the solution to the equation system.

[0012] Optionally, the step of using a reverse hierarchical strategy to divide all nodes in the dependency graph into hierarchical levels specifically includes:

[0013] Calculate the in-degree of each node and the number of non-zero elements in the corresponding matrix, and define the number of non-zero elements in the matrix as the computational weight of the node;

[0014] A computational load threshold is set. During the reverse layering process, when the cumulative computational weight of a certain layer exceeds the computational load threshold, some nodes with high computational weights are forcibly moved to the next layer to balance the computational load of each layer.

[0015] Optionally, the step of mapping the partitioned hierarchical information and the set of elder nodes to the graphics processor and starting the parallel back-generation kernel function further includes:

[0016] The size of the thread blocks in the graphics processor is dynamically adjusted based on the number of nodes and computational weights at each level.

[0017] For multiple consecutive levels with a computational load less than a preset threshold, a multi-level fusion strategy is adopted to control the same thread block to continuously process the back-substitution computation of the multiple consecutive levels.

[0018] Optionally, the parallel back-substitution calculation of the linear equation system specifically includes:

[0019] Using the shared memory or registers of the graphics processor as accumulators, perform in-situ calculations on nodes within the set of elder nodes;

[0020] After all nodes in the set of elder nodes have been calculated, the calculation results are written to the global memory of the graphics processor all at once.

[0021] Optionally, the parallel back-substitution calculation of the linear equation system adopts a mixed-precision adaptive back-substitution algorithm;

[0022] Double-precision floating-point numbers are used for calculations in the higher-level stages of back-substitution computation, while single-precision floating-point numbers or half-precision floating-point numbers are automatically switched to in the lower-level stages of back-substitution computation.

[0023] Optionally, the method further includes a residual feedback compensation step:

[0024] After each level of calculation is completed, the root node of the elder node set is extracted for residual verification.

[0025] If the residual exceeds a preset threshold, a precision rollback mechanism is triggered, which either re-executes the subsequent calculations of the branch using double-precision floating-point numbers or initiates an iterative correction process based on double-precision floating-point numbers.

[0026] Optionally, the method for constructing the set of elder nodes is as follows:

[0027] Traverse all nodes at the current level and identify multiple parent nodes that depend on the same child node;

[0028] Multiple parent nodes that depend on the same child node are grouped into the same set of elder nodes to ensure that update operations on the same child node are executed within the same thread processing scope.

[0029] Optionally, the graphics processor performs parallel back-substitution calculations in ascending order of hierarchy, as shown in the following formula:

[0030] ;

[0031] Where x(j) represents the state variable of the j-th node to be solved, y(j) represents the j-th result vector element after the forward operation of the linear equation system, u(j,i) represents the element in the j-th row and i-th column of the upper triangular transformation matrix, x(i) represents the state variable of the i-th node that has been solved, and u(j,j) represents the diagonal element in the j-th row and j-th column of the upper triangular transformation matrix.

[0032] Optionally, the method further includes a data compression and storage step:

[0033] By leveraging the local clustering characteristics of the power grid topology, the index of the elder node set is differentially encoded and compressed for storage.

[0034] When the graphics processor performs back-substitution calculations, the sparse structure pointer of the Jacobian matrix is ​​locked through a unified memory addressing mechanism to avoid redundant copying of data between host memory and device memory.

[0035] Optionally, the thread allocation logic for the parallel back-substitution kernel function is as follows:

[0036] The number of threads in a graphics processor is equal to the number of parent nodes in the current level set;

[0037] Each thread is independently responsible for the computation tasks of all nodes within a set of elder nodes, and processes the nodes sequentially within the set to eliminate memory write conflicts.

[0038] Compared with the prior art, the present invention has the following beneficial effects:

[0039] This invention discloses a parallel back-substitution method for power grid equations using a graphics processor based on a reverse path tree. Through a data compression and storage step, it leverages the local clustering characteristics of the power grid topology to differentially encode and compress the indices of the elder node set, significantly reducing GPU memory usage. Furthermore, a unified memory addressing mechanism locks the sparse structure pointers of the Jacobian matrix, avoiding redundant data copying between host memory and device memory, thus greatly improving data access efficiency. This invention employs a hybrid precision adaptive back-substitution algorithm, using double-precision floating-point operations at higher levels and automatically switching to single-precision or half-precision floating-point operations at lower levels. This ensures numerical stability while fully utilizing the high throughput advantage of low-precision computation. In addition, this invention introduces a residual feedback compensation step. After each level of computation, the root node of the elder node set is checked for residuals. When the residual exceeds a preset threshold, a precision fallback mechanism is triggered, re-executing the calculation using double-precision floating-point operations or initiating iterative correction, effectively preventing solution divergence caused by precision loss. Finally, this invention, through optimized thread allocation logic, enables each thread to independently handle the calculation of a set of elder nodes and process them in order internally, completely eliminating memory write conflicts and achieving efficient parallel computing without atomic operations. Attached Figure Description

[0040] Figure 1 This is a flowchart of the parallel back-substitution method for power grid equations based on a reverse path tree in an embodiment of the present invention.

[0041] Figure 2 This is a schematic diagram of the elimination operation with p as the axis in an embodiment of the present invention;

[0042] Figure 3 This is a schematic diagram of the original matrix in an embodiment of the present invention;

[0043] Figure 4 This is a schematic diagram of the multi-threaded concurrent computation process in an embodiment of the present invention;

[0044] Figure 5 This is a schematic diagram illustrating the construction of the reverse path tree and the set of elder nodes in an embodiment of the present invention. Detailed Implementation

[0045] 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.

[0046] In one embodiment, such as Figure 1 As shown, a parallel back-substitution method for power grid equations based on a reverse path tree graphical processor is provided, which is then applied to... Figure 1 Taking China as an example, the following specific steps will be used:

[0047] S10: Obtain the network topology and Jacobian matrix of the power system, and establish a dependency graph between nodes based on the sparse structure of the Jacobian matrix.

[0048] S20: The reverse hierarchical strategy is used to divide all nodes in the dependency graph into levels. The root node with an in-degree of zero is divided into level zero, and the level of the remaining nodes is determined by adding one to the maximum value of the levels of all its parent nodes.

[0049] S30: Based on the topological characteristic that multiple nodes in the same level depend on the same child node, nodes with the same child node are divided into a set of elder nodes.

[0050] S40: Map the partitioned hierarchical information and the set of elder nodes to the graphics processor, and start the parallel back-generation kernel function.

[0051] S50: The graphics processor performs parallel back-substitution calculations on the linear equation system in descending order of hierarchy, using the set of elder nodes as the thread allocation unit, until all levels are traversed to obtain the solution to the equation system.

[0052] Specifically, such as Figure 2 As shown, the traditional back-substitution calculation process is as follows: The points on the weighted directed factor graph are updated with the values ​​obtained from previous iterations and normalization. In this graph, node numbers j start from n and decrease in descending order. The points originating from node i of all edges dependent on j are corrected using the following formula: In the sparse decomposition path graph, there are data dependencies between parent and child nodes. In back-generation computation, the child nodes can only begin computation after the parent node, i.e., the receiving node of node i, has been completed.

[0053] Traditional back-substitution computation based on sequential path tree processes computation from the perspective of the originating point. For example, when processing node i, the computation process uses the position values ​​of the already processed receiving points j (j>i) connected to node i, which are the upper triangular elements of node i. On the one hand, all receiving points of node i may be at different levels of the path tree, which brings processing difficulties to parallel computing based on the graphics processor. On the other hand, there is also a serious problem of overlapping computation. That is, in the multi-threaded computation of parallel computing, if all receiving points of node i are assigned to different threads, the parallel computation will modify the position of node i, resulting in multiple threads rewriting the same memory, which will cause serious deviations in the calculation results.

[0054] like Figure 3 As shown, back-generation computation based on the ascending path thread tree suffers from parallel overlapping computation issues. Taking the back-generation computation of node 2 as an example, when processing node 2, node 2 is connected to node 3. For the right-hand element of node 2 There has been an update calculation; node 2 is now connected to node 6. For the right-hand element of node 2 There has been an update calculation; node 2 is now connected to node 8. For the right-hand element of node 2 There are updated calculations, and the processing results are cumulatively updated, that is... .

[0055] In GPU-based matrix parallel computation, according to the traditional computation order of the ascending sparse path thread tree, nodes 3 and 6 both belong to the first layer and have no interdependence, theoretically allowing for parallel computation. In GPU-based matrix parallel computation, threads are allocated based on the number of nodes within a layer, allowing parallel computation of nodes within the same layer. Logically, the results of thread 1 (to which node 3 belongs) and thread 2 (to which node 6 belongs) should be added to the right-hand element of node 2. However, node 8 belongs to a different layer, and nodes in different layers cannot be computed in parallel. Therefore, using the ascending sparse path thread tree to perform back-substitution computation on node 2 is logically flawed. Furthermore, the back-substitution computation of node 2 by nodes 3 and 6, which can be parallelized, results in overlapping computations.

[0056] like Figure 4 As shown, the graphics processing unit (GPU) employs a single-instruction, multi-threaded execution model. The same trigonometric decomposition statement or previous generation calculation statement is executed concurrently by multiple threads. Thread 1 performs the update. After obtaining the value, no notification was sent to thread 2. Therefore, the operations performed by both thread 1 and thread 2 were to retrieve the value. The initial value is used, and the final calculation result is also the result of a single thread, rather than the expected cumulative result. In other words, there is an issue of overlapping calculations between threads. That is, when processing unrelated parent nodes at the same level in parallel, different threads operate on the same memory block, i.e., the information of the child nodes, at the same time, and the final calculation result deviates from the expectation.

[0057] In CUDA, a general-purpose parallel computing architecture for graphics processing units, a mutex can be implemented using a device function. This device function contains atomic compare-and-swap functions, and the mutex is then implemented by calling this device function through a parallel back-substitute kernel function, ensuring the independence of memory access during multi-threaded parallelism. However, the introduction of CUDA mutexes increases the execution time of the parallel back-substitute kernel function, as each thread must execute these instructions within the function. Furthermore, once a memory block is locked, all threads executing the current statement must wait for the memory block to be unlocked before continuing execution, significantly reducing the efficiency of parallel computing.

[0058] In this invention, firstly, the network topology and Jacobian matrix of the power system are obtained. Specifically, the connection relationships of each node and branch in the power grid can be obtained by reading the power system model file, thus forming the network topology. Based on this network topology and system operating parameters, a Jacobian matrix describing the nonlinear equations of the power system is constructed. The Jacobian matrix typically has a highly sparse characteristic, meaning that most elements in the matrix are zero. Next, a dependency graph between nodes is established based on the sparse structure of the Jacobian matrix. In the dependency graph, each node represents a row or column of the Jacobian matrix. If the element in the i-th row and j-th column of the Jacobian matrix is ​​not zero, and i is not equal to j, then a directed edge from node i to node j is established in the dependency graph, or a corresponding dependency relationship is established according to the specific elimination order. This dependency graph intuitively reflects the order in which the calculated values ​​of the node variables depend on each other during the back-substitution solution process.

[0059] Subsequently, a reverse-order hierarchical strategy is used to partition all nodes in the dependency graph into different levels. The specific execution process of the reverse-order hierarchical strategy is as follows: First, the dependency graph is traversed to identify all nodes with an in-degree of zero; these nodes are the root nodes, and they are assigned to level zero. The root node represents the variable at the starting point of the dependency relationship during the back-generation process. For the remaining nodes, their level is determined by their parent nodes. Specifically, for any non-root node, all its directly connected parent nodes are retrieved, and the level numbers of these parent nodes are obtained. The maximum value is selected, and this maximum value is incremented by one to become the level of the current node. Through this top-down traversal method, all nodes in the dependency graph can be divided into different levels.

[0060] After hierarchical partitioning, based on the topological characteristic that multiple nodes in the same level depend on the same child node, nodes with the same child node are grouped into a set of elder nodes. In the dependency graph, there may be multiple nodes in the same level that jointly depend on the same node or the same group of nodes in the next level. These parent nodes that depend on the same child node constitute the set of elder nodes. The purpose of partitioning into elder node sets is to optimize and align thread bundles in subsequent parallel computing. Grouping nodes with the same data access patterns or computational logic into one category is beneficial for the graphics processor to efficiently schedule threads.

[0061] Next, the partitioned hierarchical information and the set of elder nodes are mapped to the graphics processing unit (GPU), and the parallel back-substitution kernel function is initiated. Specifically, the CPU transfers an array containing node hierarchy numbers, elder node set indices, and non-zero elements of the Jacobian matrix to the GPU's video memory. Subsequently, the CPU calls the pre-compiled parallel back-substitution kernel function. This kernel function configures the corresponding thread grid and thread blocks.

[0062] In response to the kernel function's startup instruction, the graphics processing unit (GPU) performs parallel back-substitution calculations on the linear equations, using the set of elder nodes as the thread unit, in descending order of hierarchy. Multiple stream processors within the GPU work simultaneously. When calculating a specific level, the GPU allocates thread bundles to process the set of elder nodes at that level. Different threads within the same thread bundle correspond to different nodes in the elder node set, executing back-substitution calculation instructions in parallel. The calculation process typically involves reading the value of the corresponding element in the constant term vector on the right, subtracting the product of the calculated child node variable value and the corresponding Jacobian matrix element, and finally dividing by the diagonal element to obtain the solution for the current node. The GPU continuously repeats this process, with the level number decreasing progressively, until all levels in the dependency graph have been traversed. When the level number decreases to level zero, the complete solution vector of the linear equations is obtained.

[0063] It is important to note that when searching the sparse decomposition path tree, the parent nodes of nodes at the same level are called elder nodes, and a set of elder nodes is established. The back-generation calculations of nodes in the elder node set are executed sequentially in the same thread.

[0064] This embodiment accurately resolves data dependencies in the Jacobian matrix back-substitution process by constructing a dependency graph and performing reverse hierarchical layering. Utilizing the reverse hierarchical layering strategy, the originally serial back-substitution process is transformed into a quasi-parallel process performed hierarchically. Furthermore, by identifying the set of elder nodes, nodes with the same topological characteristics are aggregated, enabling the graphics processor to fully utilize the parallel execution characteristics of thread bundles when scheduling threads, reducing thread divergence and improving instruction execution efficiency. This method significantly improves the solution speed of linear equation systems for large-scale power systems, meeting the stringent computational efficiency requirements of real-time simulation and online analysis of power systems.

[0065] In one embodiment, step S20, namely, the hierarchical division of all nodes in the dependency graph using the reverse hierarchical strategy, further includes the following steps:

[0066] S21: Calculate the in-degree of each node and the number of non-zero elements in the corresponding matrix, and define the number of non-zero elements in the matrix as the computational weight of the node.

[0067] S22: Set a computational load threshold. During the reverse layering process, when the cumulative computational weight of a certain layer exceeds the computational load threshold, some nodes with high computational weights will be forcibly moved to the next layer to balance the computational load of each layer.

[0068] Specifically, to further optimize parallel computing efficiency, this invention introduces a dynamic adjustment mechanism based on computational load. First, the system traverses each node in the dependency graph, counting the in-degree of each node and the number of non-zero elements in the Jacobian matrix corresponding to that node. The number of non-zero elements directly reflects the number of floating-point operations required by the node during back-substitution computation. Therefore, this embodiment defines the number of non-zero elements as the computational weight of the node, which quantifies the computational complexity of a single node during the solution process.

[0069] After determining the computational weight of each node, the system pre-sets a computational load threshold. This threshold can be configured based on the ideal parallel processing capability of a single thread block or thread bundle of the graphics processor. During the process of dividing the system layer by layer from the root node according to the aforementioned reverse hierarchical principle, the system continuously accumulates the computational weights of all nodes within the current layer to obtain the cumulative computational weight of that layer.

[0070] When the cumulative computational weight of a certain layer does not exceed the computational load threshold, all nodes within that layer maintain their original hierarchical division. However, when it is detected that the cumulative computational weight of a certain layer is about to exceed or has already exceeded the set computational load threshold, it indicates that the computational task at that layer is too heavy. Forcing parallel processing at the same layer may overload some of the graphics processor's computational resources while leaving others idle, thus causing a load imbalance problem. At this time, the system will trigger a load balancing mechanism to select some nodes with high computational weights from the current layer and forcibly relocate these nodes to the next layer. The nodes forcibly relocated to the next layer still satisfy the logic of reverse hierarchical layering, that is, they will still perform computations after their child nodes have completed their computations, but their execution is postponed by one step at the logical level.

[0071] In this way, a large number of computational tasks that were originally piled up at a single level are rationally split into adjacent levels. When the graphics processing unit (GPU) executes the parallel back-substitution kernel function, the computational load of each level is more balanced, avoiding thread bundle waiting or resource contention caused by excessive tasks at a single level. This hierarchical dynamic adjustment strategy based on computational weights effectively improves the hardware utilization of the GPU when processing large-scale sparse matrix back-substitution, ensuring the smoothness and stability of the parallel computing process, thereby further shortening the overall solution time of the linear equations of the power system.

[0072] In one embodiment, step S40, which involves mapping the partitioned hierarchical information and the set of elder nodes to the graphics processor and starting the parallel back-generation kernel function, further includes the following steps:

[0073] S41: Dynamically adjust the size of the thread block in the graphics processor based on the number of nodes and the computational weight at each level.

[0074] S42: For multiple consecutive levels with a computational load less than a preset threshold, a multi-level fusion strategy is adopted to control the same thread block to continuously process the back-substitution computation of the multiple consecutive levels.

[0075] Specifically, in the process of mapping the partitioned hierarchical information and the set of elder nodes to the graphics processor and initiating the parallel back-substitution kernel function, this invention also includes dynamic optimization configuration of the graphics processor execution parameters to further improve the utilization efficiency of hardware resources. Specifically, the system dynamically adjusts the size of the thread blocks in the graphics processor based on the number of nodes and computational weights at each level. During actual execution, before calling the parallel back-substitution kernel function, the central processing unit (CPU) traverses the metadata of all levels, analyzing the total number of nodes and the total computational weight of each level. For levels with a large number of nodes and high computational weights, the system configures a larger thread block size to fully utilize the streaming multiprocessor resources of the graphics processor and achieve high-concurrency parallel computing. Conversely, for levels with a small number of nodes or low computational weights, the system automatically reduces the thread block size to avoid allocating too many idle threads and wasting resources. This mechanism of dynamically adjusting the thread block size based on hierarchical characteristics ensures that the graphics processor can maintain a high thread utilization rate and execution efficiency when handling computational tasks of different scales.

[0076] Furthermore, to address the potential issue of uneven hierarchical granularity in the dependency graph, this embodiment employs a multi-level fusion strategy. Specifically, the system presets a computational threshold to identify levels with lower computational requirements. During the hierarchical mapping phase, the system scans the hierarchical sequence and identifies multiple consecutive levels with computational requirements less than the preset threshold. These levels typically contain fewer nodes; starting a kernel function or performing global synchronization for each level separately would incur significant startup overhead and synchronization latency, leading to reduced computational efficiency. Therefore, the system uses a multi-level fusion strategy to logically merge these multiple consecutive levels with lower computational requirements into a single execution unit.

[0077] During parallel back-substitution computation, the graphics processing unit (GPU) controls the same thread block to continuously process multiple consecutive levels of back-substitution computation. This means that after completing the computation of one small level, the thread block does not need to return control to the central processing unit (CPU) or wait for a global synchronization signal. Instead, it directly reads the data for the next level from the GPU's registers or shared memory and continues execution. This approach effectively masks the synchronization overhead between levels and reduces the number of kernel function launches. Through this multi-level fusion strategy, not only is the scheduling burden on the GPU reduced, but the data hit rate in the cache is also improved, thereby significantly enhancing the overall throughput and response speed of parallel back-substitution solution of linear equations.

[0078] In one embodiment, step S40, i.e., the method for constructing the set of elder nodes, further includes the following steps:

[0079] S43: Traverse all nodes at the current level and identify multiple parent nodes that depend on the same child node.

[0080] S44: Group multiple parent nodes that depend on the same child node into the same set of elder nodes to ensure that update operations on the same child node are executed within the same thread processing scope.

[0081] Specifically, such as Figure 2 As shown, constructing the set of parent nodes is a key step in optimizing thread scheduling in parallel back-substitution computation. First, the system traverses all nodes at the current level, deeply analyzing the topological connection characteristics in the dependency graph. During this traversal, the algorithm identifies multiple parent nodes that depend on the same child node. In the context of back-substitution of the Jacobian matrix in power systems, this means that multiple equation variables at the current level collectively depend on the same variable in the next level.

[0082] After identifying nodes that share a common dependency, the system groups multiple parent nodes that depend on the same child node into the same set of elder nodes. This grouping process is based on the principle of locality of reference in data dependencies. By aggregating parent nodes that jointly affect the same child node, the system ensures that update operations on the same child node are executed within the same thread.

[0083] Specifically, when the graphics processing unit (GPU) performs parallel computation, nodes belonging to the same set of parent nodes are assigned to the same thread bundle or adjacent computing units for processing. Since these parent nodes need to read and update the same child node's data during back-generation computation, grouping them into the same set and assigning them to the same thread maximizes the utilization of the GPU's cache or shared memory. When a thread processes the first node in the set, it loads the child node's data into the cache. Subsequently, when processing other nodes in the set, it can directly reuse the cached data without repeatedly reading from low-bandwidth global memory.

[0084] This construction method effectively reduces global memory access conflicts and redundant reads, avoids pipeline stalls caused by different thread bundles competing for the same memory address at the same time, and significantly improves the data locality and parallel execution efficiency of the graphics processor when processing sparse matrix back substitution by spatially aggregating computational tasks with the same data access patterns, thereby accelerating the solution speed of the overall system of equations.

[0085] In one embodiment, step S50, namely the parallel back-substitution calculation of the linear equation system, further includes the following steps:

[0086] S51: Utilizes the shared memory or registers of the graphics processor as an accumulator to perform in-situ calculations on nodes within the elder node set.

[0087] S52: After all nodes in the set of elder nodes have been calculated, the calculation results are written to the global memory of the graphics processor all at once.

[0088] Specifically, in the process of mapping the partitioned hierarchical information and the set of elder nodes to the graphics processor and initiating the parallel back-substitution kernel function, this embodiment also includes dynamic optimization configuration of the graphics processor execution parameters to further improve the utilization efficiency of hardware resources. Specifically, the system dynamically adjusts the size of the thread blocks in the graphics processor based on the number of nodes and computational weights at each level. During actual execution, before calling the parallel back-substitution kernel function, the central processing unit traverses the metadata of all levels, analyzing the total number of nodes in each level and the total computational weight of that level. For levels with a large number of nodes and high computational weights, the system configures a larger thread block size to fully utilize the streaming multiprocessor resources of the graphics processor and achieve high-concurrency parallel computing. Conversely, for levels with a small number of nodes or low computational weights, the system automatically reduces the thread block size to avoid allocating too many idle threads and wasting resources. This mechanism of dynamically adjusting the thread block size based on hierarchical characteristics ensures that the graphics processor can always maintain a high thread utilization rate and execution efficiency when handling computational tasks of different scales.

[0089] Furthermore, to address the potential issue of uneven hierarchical granularity in the dependency graph, this embodiment employs a multi-level fusion strategy. Specifically, the system presets a computational threshold to identify levels with lower computational requirements. During the hierarchical mapping phase, the system scans the hierarchical sequence and identifies multiple consecutive levels with computational requirements less than the preset threshold. These levels typically contain fewer nodes; starting a kernel function or performing global synchronization for each level separately would incur significant startup overhead and synchronization latency, leading to reduced computational efficiency. Therefore, the system uses a multi-level fusion strategy to logically merge these multiple consecutive levels with lower computational requirements into a single execution unit.

[0090] During parallel back-substitution computation, the graphics processing unit (GPU) controls the same thread block to continuously process multiple consecutive levels of back-substitution computation. This means that after completing the computation of one small level, the thread block does not need to return control to the central processing unit (CPU) or wait for a global synchronization signal. Instead, it directly reads the data for the next level from the GPU's registers or shared memory and continues execution. This approach effectively masks the synchronization overhead between levels and reduces the number of kernel function launches. Through this multi-level fusion strategy, not only is the scheduling burden on the GPU reduced, but the data hit rate in the cache is also improved, thereby significantly enhancing the overall throughput and response speed of parallel back-substitution solution of linear equations.

[0091] In one embodiment, in step S50, the parallel back-substitution calculation of the linear equation system is performed using a mixed-precision adaptive back-substitution algorithm; double-precision floating-point numbers are used for calculation in the high-level stage of the back-substitution calculation, and the calculation is automatically switched to single-precision floating-point numbers or half-precision floating-point numbers in the low-level stage of the back-substitution calculation.

[0092] Specifically, the parallel back-substitution calculation of the linear equation system employs a mixed-precision adaptive back-substitution algorithm. The core of this algorithm lies in dynamically adjusting the precision format of floating-point operations based on the numerical sensitivity and error propagation characteristics of different levels during the back-substitution calculation. Specifically, in the higher-level stages of the back-substitution calculation, the system enforces the use of double-precision floating-point numbers for operations. Higher levels typically correspond to levels near the root node in the dependency graph. The calculation results at these levels often directly determine the final magnitude and overall precision of the equation system's solution and are prone to accumulating rounding errors from lower levels. Therefore, utilizing the high significant digits of double-precision floating-point numbers ensures extremely high numerical stability along the critical path of the back-substitution calculation, preventing solution divergence or result distortion due to insufficient precision.

[0093] As the back-substitution process progresses towards the leaf nodes, when the computation enters the lower-level stage of the back-substitution calculation, the algorithm automatically switches to single-precision floating-point or half-precision floating-point operations. Lower-level nodes typically correspond to variables with low coupling or small numerical contributions in the equation system, and their relative impact on the overall solution error is weak. At this stage, the system detects the level depth or a preset precision switching threshold, triggering a data type conversion instruction to migrate subsequent computational tasks to the tensor cores or stream processors in the graphics processing unit (GPU) specifically designed for low-precision operations.

[0094] This hybrid precision adaptive back-substitution algorithm achieves two key benefits. First, by retaining double-precision operations at higher levels, it ensures the convergence of the linear equations and the accuracy of the final results, meeting the stringent high-precision requirements of power system computation. Second, by introducing single-precision or half-precision operations at lower levels, it fully leverages the significantly increased throughput of modern graphics processors in low-precision mode. This not only substantially reduces the bandwidth consumption of the graphics processor's video memory and lowers the power consumption for data transfer, but also significantly shortens the total back-substitution computation time, achieving an optimal balance between computational accuracy and efficiency.

[0095] In one embodiment, step S50, namely the parallel back-substitution calculation of the linear equation system, further includes a residual feedback compensation step:

[0096] S53: After each level of calculation is completed, extract the root node of the elder node set for residual verification.

[0097] S54: If the residual exceeds the preset threshold, the precision rollback mechanism is triggered, and the subsequent calculations of the branch are re-executed in double-precision floating-point numbers or an iterative correction process based on double-precision floating-point numbers is started.

[0098] Specifically, the method further includes a residual feedback compensation step, which aims to address potential accuracy loss during mixed-precision computation and ensure the reliability of the final solution. In practice, the system performs real-time quality monitoring after each level of computation is completed. During this process, the computation unit extracts the root node from the set of elder nodes for residual verification. The root node typically represents the key state variable of that computational branch, and the accuracy of its value is crucial to the convergence of the overall equation system. The specific operation of residual verification involves substituting the currently calculated nodal solution into the original linear equation system and calculating the resulting imbalance, i.e., the residual norm.

[0099] The system compares the calculated residuals with a preset threshold. The preset threshold is a small positive number pre-set based on the accuracy requirements of power system state estimation or power flow calculation. If the residual does not exceed the preset threshold, it indicates that the calculation accuracy at the current level meets the requirements, and the system will continue to execute subsequent calculations according to the original strategy. If the residual exceeds the preset threshold, it indicates that excessive rounding errors have accumulated in previous lower-level or current-level low-precision calculations, leading to decreased numerical stability. In this case, the system will immediately trigger the accuracy rollback mechanism.

[0100] Once the precision fallback mechanism is activated, the graphics processor will take corrective measures for that computational branch. Specifically, the system will re-execute subsequent calculations of that branch using double-precision floating-point numbers. This means that multiplication and addition operations originally planned to use single-precision or half-precision will be forcibly upgraded to double-precision operations to utilize the wider mantissa width of double-precision floating-point numbers to suppress further spread of errors. Alternatively, in another implementation, the system initiates an iterative correction process based on double-precision floating-point numbers. This process uses double-precision arithmetic to perform local Newton-Raphson iterations or Gauss-Seidel iterations on the current calculation results to quickly correct the deviation value until the residual converges to an acceptable range.

[0101] By introducing a residual feedback compensation step, this embodiment achieves a dynamic balance between computational accuracy and efficiency. In regions with good numerical performance, the system maintains low-precision computation to ensure high throughput; while in numerically sensitive or ill-conditioned regions, the system automatically reverts to a high-precision double-precision mode as a fallback. This mechanism effectively avoids the risk of overall solution divergence caused by a single low-precision computation, maximizing the parallel computing performance of the graphics processor while ensuring the robustness of solving linear equation systems.

[0102] Furthermore, the parallel solution steps for the power grid power flow linear equations based on the sparse decomposition path reverse thread tree of elder nodes are as follows: Initially, the level (i) of all nodes is set to 0; then, all nodes are traversed, and nodes with an in-degree of 0 are assigned the level 0 value (level(i) = 0); for nodes with an in-degree greater than 0, the level of node i is defined as level(i) = max(level(j) + 1), where j is the column number of the non-zero, non-diagonal element in the i-th row of the upper triangular matrix U; the sparse decomposition path thread tree is established according to the above method; nodes at the same level are traversed, and nodes with the same child nodes are placed into the elder node set; back-substitution calculation is performed from level i = 0 to level i = n-1, and the number of graphics processor threads is allocated according to the number of elder node sets at the current level at each level; parallel calculation is performed at the current level.

[0103] The basic steps of back-substitution calculation are as follows: For n-1 layers, back-substitution calculation is performed from i=0 to i=-(n-1). At each layer, a corresponding GPU thread is allocated according to the number of nodes in that layer for parallel computation. The back-substitution calculation approach for GPU-based reverse road trees is to start from the perspective of the receiving node. When processing node j, the influence of receiving node j on the sending node i (j>i) connected to node j is calculated using the lower triangular elements of row j.

[0104] Following the reverse path thread tree, the parent nodes 3 and 6 of node 2 are at the same level, while node 8 is at a different level. According to the calculation order, the influence of node 8 on all its child nodes is processed first; then the influence of the parent node set (3, 6) on its child nodes is processed.

[0105] In one embodiment, in step S50, the graphics processor performs parallel back-substitution computation in descending order of hierarchy.

[0106] Specifically, the graphics processor performs parallel back-substitution calculations in descending order of hierarchy, the core of which lies in executing mathematical operations based on the upper triangular transformation matrix. In particular, this calculation process follows a specific back-substitution formula, namely… In this formula, x(j) represents the state variable of the j-th node to be solved, which corresponds to key physical quantities such as voltage amplitude or phase angle of a specific node in the power system. y(j) represents the j-th result vector element after the linear equations are forward-propagated. This result vector element is intermediate data transmitted from the central processing unit to the graphics processing unit after completing the forward operations of the lower triangular matrix. u(j, i) represents the element in the j-th row and i-th column of the upper triangular transformation matrix, which reflects the electrical coupling relationship between node j and node i. x(i) represents the state variable of the i-th node that has been solved, i.e., the variable value of the child node that is at a lower level in the dependency graph and has been calculated first. u(j, j) represents the diagonal element in the j-th row and j-th column of the upper triangular transformation matrix, serving as the normalization factor for the division operation.

[0107] In actual parallel computing, multiple threads within the graphics processor simultaneously read their respective upper triangular transformation matrix elements and forward calculation results. Each thread first performs a summation operation, multiplying the state variables x(i) of all solved child nodes i of the current node j with the corresponding matrix element u(j, i) and accumulating the results. Then, the thread subtracts the accumulated value from the forward calculation result y(j) to obtain the net injection amount or residual of the current node. Finally, the net injection amount is divided by the diagonal element u(j, j) to accurately calculate the state variable x(j) of the j-th node to be solved.

[0108] By strictly adhering to the aforementioned back-substitution formula, this embodiment can fully utilize the parallel computing architecture of the graphics processor to quickly solve the sparse linear equations of the power system. This method transforms complex matrix operations into a large number of independent thread-level arithmetic and logical operations, significantly improving the computational throughput of the back-substitution stage while ensuring the accuracy of node state calculations, thus providing a reliable mathematical foundation for real-time power system simulation.

[0109] In one embodiment, the parallel back-substitution method for power grid equations based on reverse path trees further includes a data compression and storage step: utilizing the local clustering characteristics of the power grid topology, the index of the elder node set is differentially encoded and compressed for storage; when the graphics processor performs back-substitution calculations, the sparse structure pointer of the Jacobian matrix is ​​locked through a unified memory addressing mechanism to avoid redundant copying of data between host memory and device memory.

[0110] Specifically, to address the issues of excessive data storage and transmission delays in large-scale power grid topology data, this invention first utilizes the local clustering characteristic of power grid topologies. In power systems, tightly connected electrical nodes often exhibit a continuous or near-continuous distribution in their numbering; this characteristic is known as local clustering. Based on this, the system uses differential encoding and compressed storage for the index of the elder node set.

[0111] The specific operation of differential encoding in this invention is that, when storing the index numbers of each node in the set of elder nodes, the absolute number of the node is not stored directly. Instead, the difference between the current node index and the previous node index is stored. Due to the local clustering characteristics, these differences are usually very small integers. Therefore, the system can use data types with shorter bit widths to record these differences, thereby significantly reducing the storage space required to describe the power grid topology.

[0112] To further eliminate data transfer bottlenecks during back-substitution computation by the graphics processing unit (GPU), the system uses a unified memory addressing mechanism to lock the sparse structure pointer of the Jacobian matrix. This unified memory addressing mechanism allows the CPU and GPU to share the same virtual address space. Through this mechanism, the system fixes the sparse structure pointer of the Jacobian matrix at a specific location in physical memory and establishes a unified mapping relationship.

[0113] When the graphics processor executes the back-substitute kernel function, it can directly access the sparse data of the Jacobian matrix located in host memory through this pointer, without having to explicitly copy the entire matrix data from host memory to the graphics processor's device memory beforehand. This mechanism effectively avoids redundant data copying between host memory and device memory.

[0114] This invention significantly reduces the GPU memory resources required to store the power grid topology and Jacobian matrix, enabling the graphics processor to handle larger-scale power grid models. Simultaneously, by utilizing a unified memory addressing mechanism, it avoids cumbersome and time-consuming memory copy operations, greatly reducing latency in the data preparation phase and ensuring that the graphics processor's computing core can quickly acquire the required sparse matrix data, thereby improving the overall real-time performance of parallel back-substitution computation.

[0115] In one embodiment, the thread allocation logic of the parallel back-substitution kernel function is as follows: the number of threads in the graphics processor is equal to the number of elder nodes in the current level; each thread is independently responsible for the calculation tasks of all nodes in an elder node set, and processes the nodes sequentially within the set to eliminate memory write conflicts.

[0116] Specifically, to address the memory write conflict problem caused by data dependencies in traditional parallel computing, this invention first determines the size of the thread grid based on the topological characteristics of the current layer when the graphics processor initiates parallel back-generation computation. This thread allocation logic stipulates that the number of threads in the graphics processor is equal to the number of elder node sets in the current layer. This means that each elder node set is mapped to an independent execution thread in the graphics processor, achieving a one-to-one correspondence between sets and threads.

[0117] During the actual computation execution phase, each thread independently handles the computation tasks for all nodes within a set of elder nodes. Since the elder node set is constructed based on a reverse path tree, nodes within the set have sequential dependencies on the back-generation path. Therefore, when performing computation, the thread processes nodes sequentially within the set. The thread first reads the data of the lowest-level node in the set, completes the computation, and then uses the result to compute the next node with a slightly higher level within the same set. This serial processing method strictly adheres to the order of data dependencies.

[0118] Through the above allocation logic, this embodiment effectively eliminates memory write conflicts. In traditional fine-grained parallel schemes, multiple threads may simultaneously attempt to update the same parent node, leading to data contention. However, in this invention, since multiple parent nodes that depend on the same child node are grouped into the same set of elder nodes and updated sequentially and exclusively by the same thread, there is no situation where multiple threads simultaneously write to the same memory address. This mechanism eliminates the need for complex atomic operation instructions or locking mechanisms, ensuring the correctness of the calculation results while avoiding the performance overhead of atomic operations, significantly improving the parallel efficiency of the graphics processor during the back-generation phase.

[0119] For example, taking a large power grid in a certain region as an example, the power grid has 3956 computing nodes. Through the sparse decomposition path thread layering described above, the -29th layer has 1005 nodes. The calculations between these nodes are independent of each other. Therefore, the matrix decomposition calculation kernel function can theoretically be allocated to 1005 threads. However, through analysis, it was found that among these 1005 nodes, 132 nodes have the same child nodes. If a thread mutex lock implemented by the atomic compare and swap function in CUDA is used, although it can guarantee the independence of memory access during multi-threaded parallelism, after the CUDA mutex locks a certain memory, all threads executing the current statement must wait for the memory to be unlocked before they can continue execution, which greatly reduces the efficiency of parallel computing. In this invention, the 132 mutex lock threads are serial; after their execution is completed, other threads can proceed in parallel. The test was conducted on a Linux server with a 2.1GHz CPU, 8GB of memory, and an NVIDIA A40 GPU. The computation time for layer 0 was compared using the average of multiple calculations: the mutex function took 2059ms, while the method of this invention took 32ms. This invention significantly improves computational efficiency while maintaining accuracy.

[0120] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided by this invention can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and RAMbus dynamic RAM (RDRAM), etc. Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is used as an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the system can be divided into different functional units or modules to complete all or part of the functions described above.

[0121] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A parallel back-substitution method for power grid equations based on a reverse path tree graphical processor, characterized in that, include: Obtain the network topology and Jacobian matrix of the power system, and establish a dependency graph between nodes based on the sparse structure of the Jacobian matrix; A reverse hierarchical strategy is used to divide all nodes in the dependency graph into levels. The root node with an in-degree of zero is divided into level zero, and the level of the remaining nodes is determined by adding one to the maximum value of the levels of all its parent nodes. Based on the topological characteristic that multiple nodes in the same level depend on the same child node, nodes with the same child node are divided into a set of elder nodes. The hierarchical information and the set of elder nodes are mapped to the graphics processor, and the parallel back-generation kernel function is started. The graphics processor performs parallel back-substitution calculations on the linear equation system in descending order of hierarchy, using the set of elder nodes as the thread allocation unit, until all levels are traversed to obtain the solution to the equation system.

2. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor, as described in claim 1, is characterized in that... The step of using a reverse hierarchical strategy to divide all nodes in the dependency graph into hierarchical levels specifically includes: Calculate the in-degree of each node and the number of non-zero elements in the corresponding matrix, and define the number of non-zero elements in the matrix as the computational weight of the node; A computational load threshold is set. During the reverse layering process, when the cumulative computational weight of a certain layer exceeds the computational load threshold, some nodes with high computational weights are forcibly moved to the next layer to balance the computational load of each layer.

3. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The step of mapping the partitioned hierarchical information and the set of elder nodes to the graphics processor and starting the parallel back-generation kernel function also includes: The size of the thread blocks in the graphics processor is dynamically adjusted based on the number of nodes and computational weights at each level. For multiple consecutive levels with a computational load less than a preset threshold, a multi-level fusion strategy is adopted to control the same thread block to continuously process the back-substitution computation of the multiple consecutive levels.

4. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor, as described in claim 1, is characterized in that... The parallel back-substitution calculation of the linear equation system specifically includes: Using the shared memory or registers of the graphics processor as accumulators, perform in-situ calculations on nodes within the set of elder nodes; After all nodes in the set of elder nodes have been calculated, the calculation results are written to the global memory of the graphics processor all at once.

5. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The parallel back-substitution calculation of the linear equation system is performed using a mixed-precision adaptive back-substitution algorithm. Double-precision floating-point numbers are used for calculations in the higher-level stages of back-substitution computation, while single-precision floating-point numbers or half-precision floating-point numbers are automatically switched to in the lower-level stages of back-substitution computation.

6. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 5, characterized in that, The method also includes a residual feedback compensation step: After each level of calculation is completed, the root node of the elder node set is extracted for residual verification. If the residual exceeds a preset threshold, a precision rollback mechanism is triggered, which either re-executes the subsequent calculations of the branch using double-precision floating-point numbers or initiates an iterative correction process based on double-precision floating-point numbers.

7. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The specific method for constructing the set of elder nodes is as follows: Traverse all nodes at the current level and identify multiple parent nodes that depend on the same child node; Multiple parent nodes that depend on the same child node are grouped into the same set of elder nodes to ensure that update operations on the same child node are executed within the same thread processing scope.

8. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The graphics processor performs parallel back-substitution calculations in descending order of hierarchy, using the following formula: ； Where x(j) represents the state variable of the j-th node to be solved, y(j) represents the j-th result vector element after the forward operation of the linear equation system, u(j,i) represents the element in the j-th row and i-th column of the upper triangular transformation matrix, x(i) represents the state variable of the i-th node that has been solved, and u(j,j) represents the diagonal element in the j-th row and j-th column of the upper triangular transformation matrix.

9. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The method also includes a data compression and storage step: By leveraging the local clustering characteristics of the power grid topology, the index of the elder node set is differentially encoded and compressed for storage. When the graphics processor performs back-substitution calculations, the sparse structure pointer of the Jacobian matrix is ​​locked through a unified memory addressing mechanism to avoid redundant copying of data between host memory and device memory.

10. The parallel back-substitution method for power grid equations based on a reverse path tree using a graphical processor according to claim 1, characterized in that, The thread allocation logic of the parallel back-substitution kernel function is as follows: The number of threads in a graphics processor is equal to the number of parent nodes in the current level set; Each thread is independently responsible for the computation tasks of all nodes within a set of elder nodes, and processes the nodes sequentially within the set to eliminate memory write conflicts.

Images