A Krylov subspace dimension reduction power flow calculation and solution system

Abstract

The present invention relates to a Krylov subspace dimension reduction power flow calculation and solution system, comprising: a power balance module, which is used to form an n-dimensional linear power balance equation group based on a Jacobian matrix and a power imbalance matrix; a Krylov subspace dimension reduction solution module, which is used to iteratively calculate the n-dimensional linear power balance equation group using a standard orthogonal basis after orthogonalization calculation to obtain an approximate optimal solution for the power correction amount in an m-dimensional subspace; a bidirectional linked list matrix compression storage module, which is used to be responsible for the matrix operation during the iterative calculation of the Krylov subspace dimension reduction solution module, and the matrix operation is operated based on matrix block division. The present invention can speed up the operation speed of the existing power flow calculation and solution, and improve the calculation efficiency and resource utilization.

Description

Technical Field

[0001] The invention relates to the technical field of power flow calculation, and in particular to a Krylov subspace dimension reduction power flow calculation and solution system. Background Art

[0002] Power flow calculation is a very important analysis and calculation of the power system. It calculates the electrical quantities of the power system under steady-state operation according to the power system wiring mode, parameters and operating conditions, including the distribution of active power, reactive power and voltage in the power grid. Power flow calculation is also the basis for system safety, stability and reliability analysis, and is used to study various problems raised in system planning and operation. For the power system under planning, power flow calculation can be used to verify whether the proposed power system planning scheme can meet the requirements of various operating modes; for the power system under operation, power flow calculation can be used to predict whether various load changes and changes in network structure will endanger the safety of the system, whether the voltage of all buses in the system is within the allowable range, whether various components in the system (lines, transformers, etc.) will be overloaded, and what preventive measures should be taken in advance when overload may occur.

[0003] With the digital transformation and upgrading of power grids, the measured values ​​in the network are constantly increasing, and the scale of power networks is constantly expanding, which puts higher requirements on the calculation of power flow distribution in power grids. At present, some research institutions and commercial companies have developed a number of professional power flow calculation software, such as EasyPower, ETAP, Southern CASS, Matpower, PSASP, PSCAD, etc. However, the existing CPU-based software algorithms are designed based on general task-based calculations on general CPUs, and there is no performance optimization specifically for the calculation characteristics of the algorithm. Instructions are executed in clock cycle order, and only one instruction can be executed in one clock. Therefore, when performing data calculation tasks, the calculation speed is slow due to the serial calculation process, which cannot meet the effectiveness of large-scale power system power flow calculations.

[0004] The existing patent document CN111740424A proposes a method for calculating power system flow based on a computation tree GPU parallel acceleration model. This method uses a computation tree GPU parallel acceleration model to accelerate the generation of Jacobian matrices and measurement corrections in power system flow calculations, and improves calculation efficiency when the system scale is large. Although this method uses the parallel computing characteristics of GPUs to improve computing efficiency when the system data volume is large, it still has shortcomings in real-time performance.

[0005] For computationally intensive tasks such as power flow calculations, FPGAs have the characteristics of pipeline parallelism and data parallelism, which can achieve hardware-level acceleration of calculations. FPGA-based algorithms are directly implemented at the hardware circuit level. According to the characteristics of the data calculation process of a specific algorithm, the algorithm is designed using a hardware circuit description language and mapped to the underlying chip to form a customized dedicated computing circuit. The data calculation task is executed in a parallel pipeline, so it is faster. Power flow calculations involve a large number of matrix calculations and solving linear equations. FPGAs are particularly suitable for processing such multi-component computationally intensive tasks. Relying on the pipeline parallel structure system, FPGAs have technical advantages over CPUs in terms of computing speed and result return latency.

[0006] The existing open literature (Nwankpa C, Johnson J, Nagvajara P, et al.FPGA hardware results for power system computation [C] / / 2009 IEEE / PES Power Systems Conference and Exposition.0 [2023-09-26]. DOI: 10.1109 / PSCE.2009.4839953) proposed a power flow calculation FPGA system based on the Newton-Raphson method, in which the solution strategy for the linear equations is the matrix LU decomposition method, using the PCI pipeline bus architecture, and using FPGA hardware to reduce the pivot search time to accelerate the LU decomposition. However, since the LU decomposition of the Jacobian matrix formed by each iterative solution of the equation system requires a large number of matrix operations, it poses a great challenge to the computing power and storage resources of the FPGA chip. Therefore, when facing large-scale power flow calculations, this solution will have the disadvantages of high storage resource consumption and slow calculation speed, and cannot effectively exert the speed advantage of FPGA in processing intensive matrix operations. Summary of the invention

[0007] The technical problem to be solved by the present invention is to provide a Krylov subspace dimension reduction power flow calculation and solution system, which can speed up the operation speed of the existing power flow calculation and solution, and improve the calculation efficiency and resource utilization.

[0008] The technical solution adopted by the present invention to solve the technical problem is: to provide a Krylov subspace dimension reduction power flow calculation and solution system, comprising:

[0009] A power balancing module, used for forming an n-dimensional linear power balancing equation group based on a Jacobian matrix and a power imbalance matrix;

[0010] A Krylov subspace dimension reduction solution module is used to iteratively calculate the n-dimensional linear power balance equations using a standard orthogonal basis after orthogonalization calculation to obtain an approximate optimal solution for the power correction in the m-dimensional subspace;

[0011] The bidirectional linked list matrix compression storage module is used for performing matrix operation during iterative calculation of the Krylov subspace dimensionality reduction solution module, wherein the matrix operation is performed based on matrix block division.

[0012] The bidirectional linked list matrix compression storage module comprises:

[0013] A matrix block unit, used for dividing the Jacobian matrix into blocks by rows to obtain a plurality of matrix blocks;

[0014] A row-block parallel computing unit is used to obtain a product result component based on the matrix block, and the product result component is a result obtained by dividing the product result of the Jacobian matrix and the column vector into row blocks.

[0015] The matrix blocking unit completes the Jacobian matrix blocking in the following manner:

[0016] Accumulate the number of non-zero elements in the Jacobian matrix row by row to obtain the accumulated value of non-zero elements;

[0017] Determine whether the current non-zero element accumulation value is greater than s+1 times the number of average matrix division elements, where s is the row number of the Jacobi matrix row-by-row division point; if the current non-zero element accumulation value is greater than s+1 times the number of average matrix division elements, then take the current row number i as a matrix block division point, and let the Jacobi matrix row-by-row division point row number s be incremented by one;

[0018] Repeat the previous step until the current row number is the last row of the Jacobian matrix, and set the elements of the remaining Jacobian matrix as a matrix block.

[0019] The row-block parallel computing unit obtains the product result component based on the matrix block in the following manner:

[0020] Calculate the number of rows of the ith matrix block, and find out the computational components of other block parallel computing units that the jth row of the ith matrix block needs to rely on;

[0021] If the number of calculation components of other block parallel computing units that the j-th row of the i-th matrix block needs to rely on is zero, directly multiply the j-th row of the i-th matrix block with the column vector to obtain the product result component;

[0022] If the number of calculation components of other block parallel computing units that the j-th row of the i-th matrix block needs to rely on is greater than zero, then find the numbers of the calculation components of the other dependent block parallel computing units, multiply the j-th row of the i-th matrix block with the column vector to obtain a first result, then multiply the j-th row of the i-th matrix block with the calculation components of the other dependent block parallel computing units and accumulate to obtain a second result, sum the first result and the second result to obtain the product result component.

[0023] The Krylov subspace dimensionality reduction solution module includes:

[0024] The parameter initialization unit is used to set the initial iteration value x0 and the dimension m of the subspace, calculate the error vector r and the inner product beta of the error vector r, and use the Gram-Schmidt orthogonalization method to calculate the vector space and output the standard orthogonal basis V m and the Heisenberg matrix H; the error vector r is standardized by the inner product beta of the error vector r, and the standardized error vector r is used as the standard orthogonal basis V m The first element of; Set the initial value of the row list head pointer register row_ptr_head_H and the row address index register row_ptr_H of the upper Heisenberg matrix H to -1, and set the row list head pointer register cow_ptr_head_H and the column list head pointer register col_ptr_head_H, the row address index register row_ptr_H and the column address index register col_ptr_H to null pointers, and set the initial values ​​of the row element number register num_row_H, the column element number register num_col_H, the data element register data_H and the data row and column label value register data_col_row_H to 0; Wherein, M is the Jacobian matrix;

[0025] The standard orthogonal basis calculation unit is used to calculate the Jacobian matrix M and the standard orthogonal basis V m The jth element in calculates the basis vector w under dimension j, and for all the j orthogonal bases that have been generated, calculates the correction vector and performs orthogonalization calculation to complete the standard orthogonal basis V m Orthogonalization calculation of ;

[0026] The approximate optimal solution calculation unit is used to calculate the standard orthogonal basis V after orthogonalization. m Calculate the approximate optimal solution x of Mx=b m , where b is the power imbalance matrix.

[0027] The standard orthogonal basis calculation unit completes the standard orthogonal basis V in the following way m Orthogonal calculation of :

[0028] The Jacobian matrix M and the orthonormal basis V m Multiply the jth element in to get the basis vector w under dimension j;

[0029] Determine the base vector w pointer register address ptr_id_w and the generated orthogonal base pointer vector address ptr_id_V i The relationship between the size of the register and the number of non-zero elements;

[0030] When the base vector w pointer register address ptr_id_w is less than or equal to the number of non-zero elements in the base vector w register, or the orthogonal base pointer vector address ptr_id_V has been generated i When it is less than or equal to the number of non-zero elements in the orthogonal basis register, record the basis vector w and the orthogonal basis V i The element position number in , and the data value corresponding to the position number;

[0031] If the element positions of the basis vector w are numbered with the orthogonal basis V i If the element position numbers in are the same, the corresponding elements are multiplied and accumulated and stored in the inner product data_H_ij; if the element position number of the basis vector w is greater than the orthogonal basis V i The element position number in and the orthogonal base pointer vector address ptr_id_V that has been generated i When the number of non-zero elements in the orthogonal base register is less than the number of non-zero elements in the orthogonal base register, the generated orthogonal base pointer vector address ptr_id_V i Plus one; if the element position number of the basis vector w is less than the orthogonal basis V i If the element position number in , and the basic vector w pointer register address ptr_id_w is less than the number of non-zero elements in the basic vector w register, the basic vector w pointer register address ptr_id_w is increased by one;

[0032] Subtract the value in the inner product memory data_H_ij from the basis vector w and the orthogonal basis V i The product of is used as the new basis vector;

[0033] Repeat the above steps until the accumulated value of the inner product data_H_ij is calculated;

[0034] When the index of the calculated orthogonal basis vector is greater than or equal to 2, find the previous element in this row of the inner product data_H_ij, and record the address of the first element of the i-th row as the linked list head pointer register row_ptr_head_H[i]. If the address of the first element of the i-th row is empty, directly point the row linked list head pointer register row_ptr_head_H[i] to the address pointer register of the first position of the inner product data_H_ij; if the address of the first element of the i-th row is not empty, perform the following steps until the previous element position number of the inner product data_H_ij is found;

[0035] Take out the next element address next_col_id of the first element of the i-th row in the row pointer register row_ptr_H, and determine whether the next element address next_col_id is empty. If the next element address next_col_id is empty, the previous element position number of the inner product data_H_ij has been found, and the element at the next element address next_col_id points to the address pointer register of the first position of the inner product data_H_ij, and exit the process; if the next element address next_col_id is not empty, the first element address of the i-th row will be updated to the next element address next_col_id.

[0036] The accumulated value of the inner product data_H_ij is stored as follows:

[0037] The address pointer register pointing to the first position of data_H is recorded as ptr_data_H, and the initial value of the address pointer register ptr_data_H is set to zero;

[0038] Determine whether the inner product data_H_ij is zero. If the inner product data_H_ij is not zero, store the inner product data_H_ij to the address pointed to by the address pointer register ptr_data_H, and store the row and column label values ​​of the inner product data_H_ij at the same time;

[0039] Determine whether the column link list head pointer register col_ptr_head_H and the row link list head pointer register row_ptr_head_H at the position (i, j) are empty; if the column link list head pointer register col_ptr_head_H[j] of the jth column is not empty, then point the current column link list head pointer to the address pointed to by the address pointer register ptr_data_H; if the row link list head pointer register row_ptr_head_H[i] of the i-th row is not empty, then point the current row link list head pointer to the address pointed to by the address pointer register ptr_data_H;

[0040] The pointer of the non-zero element number register of the row and column where (i, j) is located is incremented by one, and the pointer of the column address index register col_ptr_H is incremented by one.

[0041] The approximate optimal solution calculation unit is Calculate the approximate optimal solution, where x m is a near optimal solution, , It represents the inverse matrix of the m-dimensional Heisenberg matrix H. unit_e is an m-dimensional column vector whose first element is 1.

[0042] The technical solution adopted by the present invention to solve its technical problem is: to provide a power flow calculation system, characterized in that it includes the above-mentioned Krylov subspace dimensionality reduction power flow calculation and solution system, a judgment module and a node voltage and power calculation module, the judgment module is used to judge whether the change of the approximate optimal solution calculated by the Krylov subspace dimensionality reduction power flow calculation and solution system within a preset time is less than a preset error; the node voltage and power calculation module is used to calculate the network node flow distribution according to the approximate optimal solution when the change of the approximate optimal solution is less than a preset error.

[0043] The power flow calculation system further comprises a Jacobian matrix calculation module and a power imbalance calculation module; the Jacobian matrix calculation module is used to generate a Jacobian matrix; the power imbalance calculation module is used to generate a power imbalance matrix.

[0044] Beneficial Effects

[0045] Due to the adoption of the above technical scheme, the present invention has the following advantages and positive effects compared with the prior art: the present invention takes into account the matrix operation characteristics in the Krylov subspace iterative algorithm, designs a matrix row block parallel calculation process and method for matrix-vector multiplication by bidirectional linked list compression storage, and the bidirectional linked list matrix compression storage method can speed up the calculation speed of the Krylov subspace iterative algorithm and save the consumption of storage resources on the chip. At the same time, according to the particularity of the dynamic change of the Jacobian matrix generated in the power flow calculation, a bidirectional linked list matrix compression storage format is designed. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 It is a block diagram of a Krylov subspace dimension reduction power flow calculation and solution system according to a first embodiment of the present invention;

[0047] Figure 2 It is a block diagram of a power flow calculation system according to a second embodiment of the present invention. DETAILED DESCRIPTION

[0048] The present invention will be further described below in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and are not intended to limit the scope of the present invention. In addition, it should be understood that after reading the content taught by the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms fall within the scope limited by the appended claims of the application equally.

[0049] The first embodiment of the present invention relates to a Krylov subspace dimension reduction power flow calculation and solution system, such as Figure 1 As shown, including:

[0050] A power balancing module, used for forming an n-dimensional linear power balance equation group based on a Jacobian matrix and a power imbalance matrix, wherein n is a positive integer greater than 2;

[0051] A Krylov subspace dimension reduction solution module is used to iteratively calculate the n-dimensional linear power balance equations using a standard orthogonal basis after orthogonalization calculation to obtain an approximate optimal solution for the power correction amount in an m-dimensional subspace, where m is a positive integer greater than 2;

[0052] The bidirectional linked list matrix compression storage module is used for the matrix operation during the iterative calculation of the Krylov subspace dimensionality reduction solution module. The matrix operation is based on matrix block operation, mainly including the data storage definition of sparse matrix elements, the addition, deletion and overwriting operations of matrix elements, and the parallel calculation process of matrix-vector multiplication, inner product and matrix multiplication.

[0053] The bidirectional linked list matrix compression storage module in this implementation manner includes a matrix blocking unit and a row blocking parallel computing unit.

[0054] In order to improve the computational efficiency of matrix-vector multiplication, this implementation has m parallel computation modules for matrix-vector multiplication inside the FPGA. , calculation module P i The amount of calculation assigned to the calculation module P i The Jacobian matrix M i In order to ensure the uniform distribution of the computational effort of each computing module, the matrix partitioning unit in this embodiment partitions the Jacobian matrix by row to obtain a number of matrix blocks, wherein the number of non-zero elements in each matrix block is as close as possible. The specific partitioning method is as follows:

[0055] Assume that the number of non-zero elements in the Jacobian matrix M is num_M, where num_M is the accumulated value of the matrix register num_row or the matrix register num_col. The average number of elements in the ideal state of each matrix block after the Jacobian matrix is ​​divided into blocks is avg_row, avg_row←num_M / m. The dividing point of the Jacobian matrix M is numbered s, s←0, and the matrix dividing point is q. s , the sum of non-zero elements of the matrix that has been currently divided into blocks is acc_num_M, acc_num_M←0.

[0056] For each row number i=1 to n of the Jacobian matrix M, the number of non-zero elements in the Jacobian matrix M is accumulated row by row to acc_num_M, acc_num_M←acc_num_M+num_row[i]. Determine the size of the accumulated value of non-zero elements in the current row acc_num_M and the number of average matrix partition elements avg_row. If the current accumulated value of non-zero elements is greater than s+1 times the number of average matrix partition elements, that is, acc_num_M>(s+1)*avg_row, then record the current row number i as a matrix block partition point, q s+1 ←i, s←s+1. Repeat the above steps until i=n, ​​and finally set the remaining matrix elements as a matrix block.

[0057] After the above online matrix block calculation process, the elements can be distributed as evenly as possible. Divide the Jacobian matrix M into m matrix blocks , the i-th matrix block M i The line in is the qth order of the original Jacobian matrix M i To q i+1 -1 row element, the total number of rows contained is q i+1 -q i , the i-th matrix block M i Stored in m independent computing modules respectively Perform parallel computing.

[0058] Calculate the product of the n*n Vyakofibrillian matrix M and the n-dimensional column vector x, y=M*x, M, x, y are divided into rows and stored in c independent calculation modules respectively In the first step, we need to ik The relative coordinate values ​​of other calculation components stored in are converted into the absolute coordinate values ​​of the original Jacobian matrix. The conversion rules are as follows:

[0059] The absolute coordinate value of the original Jacobian matrix = w ik The column number value of the component unit + (the number of this component i - Initial component number s0).

[0060] Matrix-vector multiplication y=M*x can be represented by rows:

[0061] ;

[0062] Among them, the i-th result vector y i for:

[0063] ;

[0064] in, , M i and x i Stored in the i-th computing module P i In the calculation unit, other dependent components are stored in the dependent component number register w ik And dependent component number register row_id_w ij middle.

[0065] In this embodiment, the row-block parallel computing unit is used to obtain a product result component based on the matrix block, and the product result component is the result obtained by dividing the product result of the Jacobian matrix and the column vector into row blocks. The specific calculation method is as follows:

[0066] For each calculation module i=0 to c-1, calculate the block matrix M contained in the i-th calculation module i The number of rows s i ,s i =q i+1 -q i , then the block matrix M i All rows j = 1 to s i , the jth row of the ith matrix block needs to rely on the computational components of other parallel computing units in the block. ij The number of other calculation components required for the row is num_x, num_x←row_id_w ij [j+1]- row_id_w ij [j]. Determine whether the i-th row has components that depend on other computing modules, and give y ij Assign initial value, y ij ←0, 1) If num_x = 0, no other components are needed for this row, and y is calculated directly ij ←M ij ·x i , 2) If num_x>0, it is necessary to find out the number of other components that the module depends on, and record the sequence of other components as other_x_id←row_id_w ij [j], for each dependent component k=1 to num_x, record the component position number x_id, x_id←wij [other_x_id+k-1], and M ij Multiply with other components and add to y ij , that is, y ij ←y ij +M i,x_id ·x i_id . When k=num_x, M ij The cumulative sum y of all other dependent components ij The calculation is completed. Finally, calculate M in its own block i ij and x i Multiply and add to y ij , that is, y ij ←y ij +M ij ·x i At this point, the result component y of the i-th matrix block is calculated ij Repeat the above steps until all components in all blocks are calculated and stored in y i .

[0067] It is not difficult to find that the matrix row block calculation process of the bidirectional linked list compression storage and the row block parallel calculation process of the matrix-vector multiplication of this embodiment can efficiently realize the calculation speed of matrix-vector multiplication, matrix inner product, and matrix multiplication in multiple computing modules at the same time, thereby accelerating the calculation efficiency of the Krylov subspace iterative solution.

[0068] In this implementation, the Krylov subspace dimensionality reduction solution module includes a parameter initialization unit, a standard orthogonal basis calculation unit and an approximate optimal solution calculation unit.

[0069] Among them, the parameter initialization unit is used to initialize various parameters, as follows:

[0070] Set the initial iteration value x0 and the dimension m of the subspace. To simplify the calculation, set x0 to an n-dimensional 0 vector, x0←[0,0,…,0]. According to experience, the subspace dimension m is 30%~50% of the original matrix dimension n. Calculate the error vector r←b-M·x0. At the beginning of the calculation, r=b. The error vector r is stored in data_b in column compression format, and the number of non-zero elements is num_b. At the same time, calculate the inner product beta of the error vector r, beta← , beta is calculated as follows:

[0071] First, calculate the sum of squares in the error vector r, set the address pointer ptr_id_b←0, the accumulated square sum value square_sum_b←0, and for all non-zero elements i=1 to num_b in the error vector r, perform the following steps: take out the first column vector element data in b, data←data_b[ptr_id_b], calculate the accumulated square sum square_sum_b, square_sum_b←square_sum_b+data*data, and ptr_id_b←ptr_id_b+1. When r=num_b, the accumulated square sum is completed. Finally, calculate beta, beta← .

[0072] Then, the Gram-Schmidt orthogonalization method is used to calculate the vector space and output the standard orthogonal basis V m and the Heisenberg matrix H, where the orthonormal basis V m The dimension is n*(m+1), , the dimension of the Shanghai Senburg matrix H is (m+1)*m. Normalize the error vector r and store it in V1, that is, V1←r / beta. Set the initial value of the row list head pointer register row_ptr_head_H and the row pointer register row_ptr_H of the Shanghai Senburg matrix H to -1. The doubly linked list registers of the Shanghai Senburg matrix H are set as follows: the row list head pointer register cow_ptr_head_H and the column list head pointer register col_ptr_head_H, the row address index register row_ptr_H and the column address index register col_ptr_H are set to the null pointer null, the row element number register num_row_H, the column element number register num_col_H, the data element register data_H and the data row and column label value register data_col_row_H are all set to 0, and the above registers are dynamically updated during the iteration process.

[0073] The standard orthogonal basis calculation unit in this embodiment is used to calculate the Jacobian matrix M and the standard orthogonal basis V m The jth element in calculates the basis vector w under dimension j, and for all the j orthogonal bases that have been generated, calculates the correction vector and performs orthogonalization calculation to complete the standard orthogonal basis V m The orthogonal calculation is as follows:

[0074] First, calculate the basis vector w in dimension j, w←M·V j, calculate the inner product of the jth orthogonal basis i to be data_H_ij, data_H_ij←0. Determine the basic vector w pointer register address ptr_id_w and the standard orthogonal basis pointer vector address ptr_id_V that has been generated i The relationship between the size of the number of non-zero elements in the register is that when the base vector w pointer register address ptr_id_w is less than or equal to the number of non-zero elements in the base vector w register, or the orthogonal base pointer vector address ptr_id_V has been generated i Less than or equal to the number of non-zero elements in the orthogonal base register, that is, ptr_id_w ≤ num_data_w, or ptr_id_V i ≤num_V i , perform the following steps:

[0075] Record the basis vector w and the orthogonal basis V i The element position number in, that is, this_id_w←id_w[ptr_id_w], this_id_V i ←id_V i [ptr_id_V i ], and record the data value corresponding to the position number, this_data_w←data_w[ptr_id_w], this_data_V i ←data_V i [ptr_id_V i ]. Determine the relationship between the position numbers. a) If the element position number of the basis vector w is the same as the orthogonal basis V i The element position number in is the same, that is, this_id_w = this_id_V i , then perform the product operation of the corresponding elements and accumulate and store them in data_H_ij, data_H_ij←data_H_ij+ this_data_w* this_data_V i b) If the element position number of the basis vector w is greater than the orthogonal basis V i The element position number in and the orthogonal base pointer vector address ptr_id_V that has been generated i Less than the number of non-zero elements in the orthogonal base register, that is, this_id_w>this_id_V i , and ptr_id_V i < num_V i When the orthogonal base pointer vector address ptr_id_V is generated, i Plus one, that is, ptr_id_V i ←ptr_id_V i+1. c) If the element position number of the basis vector w is less than that of the orthogonal basis V i The element position number in the base vector w pointer register address ptr_id_w is less than the number of non-zero elements in the base vector w register, that is, this_id_w <this_id_V i , and ptr_id_w<num_data_w, the address of the basic vector w pointer register ptr_id_w is incremented by 1, that is, ptr_id_w←ptr_id_w+1. At the same time, the basic vector w is subtracted from the value in the inner product memory data_H_ij and the orthogonal basis V i The product of is used as the new basis vector, that is, w←w-data_H_ij*V i .

[0076] Repeat the above steps until the accumulated value of data_H_ij is calculated. The storage process of the accumulated value of data_H_ij is as follows:

[0077] The address pointer register pointing to the first position of data_H is recorded as ptr_data_H, and the initial value of the address pointer register ptr_data_H is set to zero, that is, ptr_data_H←0. Determine whether the value of the inner product data_H_ij is zero. If the inner product data_H_ij is not zero, that is, data_H_ij≠0, perform the following steps: store the inner product data_H_ij to the address pointed to by the address pointer register ptr_data_H, that is, data_H[ptr_data_H]←data_H_ij, and store the row and column number values ​​of data_H_ij at the same time, that is, data_col_row_H[ptr_data_H][63:32]←j, data_col_row_H[ptr_data_H][31:0]←i. Determine whether the column linked list head pointer register col_ptr_head_H and the row linked list head pointer register row_ptr_head_H at the position (i, j) are empty. a) If the column linked list head pointer register col_ptr_head_H[j] of the jth column is not empty, that is, col_ptr_head_H[j]≠null, then point the current column linked list head pointer to the address pointed to by the address pointer register ptr_data_H, col_ptr_head_H[j]←ptr_data_H. b) If the row linked list head pointer register row_ptr_head_H[i] of the i-th row is not empty, that is, row_ptr_head_H[i]≠null, then point the current row linked list head pointer to the address pointed to by the address pointer register ptr_data_H, that is, row_ptr_head_H[i]←ptr_data_H. After processes a) and b) are executed, the register for the number of non-zero elements in the row and column where (i, j) is located is updated, num_col_H[j]←num_col_H[j]+1, num_row_H[i]←num_row_H[i]+1, and the column pointer index register col_ptr_H is updated, col_ptr_H[ptr_data_H]←ptr_data_H+1.

[0078] Next, store the value of the row pointer register row_ptr_H. The process is as follows:

[0079] When the calculated basis vector index j ≥ 2, perform the following steps: find the previous element in the row of the inner product data_H_ij, and record the address of the first element of the i-th row as the linked list head pointer register row_ptr_head_H[i], that is, row_id←row_ptr_head_H[i], 1) If the address of the first element of the i-th row is empty, that is, row_id=null, then directly point the row linked list head pointer register row_ptr_head_H[i] to the address pointer register of the first position of the inner product data_H_ij, that is, row_ptr_head_H[i]←ptr_data_H, 2) If the address of the first element of the i-th row is not empty, that is, row_id≠null, continue to perform the following steps until the first element of the i-th row is found. The previous element position number of data_H_ij: Take out the next element address next_col_id of the first element of the i-th row in the row pointer register row_ptr_H, that is, next_col_id←row_ptr_H[row_id], and determine whether the next element address next_col_id is empty. a) If the next element address next_col_id is empty, that is, next_col_id=null, the previous element address number has been found, and the element at the next element address next_col_id points to the address pointer register of the first position of the inner product data_H_ij, that is, row_ptr_H[next_col_id]←ptr_data_H, and exit the process. b) If the next element address next_col_id is not empty, that is, next_col_id≠null, update row_id, row_id←next_col_id.

[0080] After the above process is completed, the jth vector has completed all the generated orthogonal basis vectors V i Orthogonalization calculation is performed, and finally the data element address pointer register ptr_data_H is updated, ptr_data_H←ptr_data_H+1.

[0081] Next, calculate data_H_(j+1)j, data_H_(j+1)j← , and store it in the data matrix register. The storage process is the same as above. Determine whether data_H_(j+1)j is zero. 1) If data_H_(j+1)j=0, the m-dimensional subspace calculation is completed. The number of vectors j that have been calculated is recorded as the dimension m of the subspace, m←j, and the entire subspace calculation process ends. 2) Otherwise, the standardized correction vector w is used as the orthogonal basis V of the j-th subspace. j+1 Stored in subspace V, that is, V j+1←w / data_H_(j+1)j. Finally, repeat the above steps until the m-dimensional subspace is established.

[0082] After the above process is completed, it represents the m-dimensional subspace V m The calculation has been completed. Next, the approximate optimal solution x of Mx=b is calculated by the approximate optimal solution calculation unit. m In this embodiment, the approximate optimal solution calculation unit is Compute the approximate optimal solution, where , Represents the inverse matrix of the m-dimensional Heisenberg matrix H stored in the data_H register. unit_e is an m-dimensional column vector whose first element is 1, that is, unit_e←[1,0,…,0].

[0083] It can be seen that the linear equations solving method of Krylov subspace dimensionality reduction adopted in this embodiment can solve the approximate optimal solution of the original equation in a subspace of smaller dimension, thereby improving the efficiency of power flow calculation.

[0084] A second embodiment of the present invention relates to a power flow calculation system, such as Figure 2 As shown, it includes the Krylov subspace dimensionality reduction power flow calculation and solution system of the first embodiment, a judgment module and a node voltage and power calculation module, the judgment module is used to judge whether the change of the approximate optimal solution calculated by the Krylov subspace dimensionality reduction power flow calculation and solution system within a preset time is less than a preset error; the node voltage and power calculation module is used to calculate the network node flow distribution according to the approximate optimal solution when the change of the approximate optimal solution is less than a preset error.

[0085] The power flow calculation system also includes a Jacobian matrix calculation module and a power imbalance calculation module; the Jacobian matrix calculation module is used to generate a Jacobian matrix; the power imbalance calculation module is used to generate a power imbalance matrix.

[0086] The power flow calculation system based on the Krylov subspace dimensionality reduction power flow calculation solution system proposed in this embodiment can save on-chip storage resources required for large-scale node network power flow calculations, realize efficient compressed storage and reading of matrix elements, and improve the access efficiency of matrix operations during the Krylov subspace iteration process.

Claims

1. A Krylov subspace dimension reduction power flow calculation and solution system, characterized in that: include: A power balancing module, used for forming an n-dimensional linear power balancing equation group based on a Jacobian matrix and a power imbalance matrix; The Krylov subspace dimension reduction solution module is used to iteratively calculate the n-dimensional linear power balance equations using the standard orthogonal basis after orthogonalization calculation to obtain an approximate optimal solution for the power correction amount in the m-dimensional subspace; the Krylov subspace dimension reduction solution module includes: The parameter initialization unit is used to set the initial iteration value x0 and the dimension of the subspace m, calculate the error vector r and the inner product beta of the error vector r, and use the Gram-Schmidt orthogonalization method to calculate the vector space {r,Mr,M 2 r,...,M m-1 r}, and output the standard orthogonal basis V m and the Heisenberg matrix H; the error vector r is standardized by the inner product beta of the error vector r, and the standardized error vector r is used as the standard orthogonal basis V m The first element; Set the initial value of the row list head pointer register row_ptr_head_H and the row address index register row_ptr_H of the upper Heisenberg matrix H to -1, and set the row list head pointer register cow_ptr_head_H and the column list head pointer register col_ptr_head_H, the row address index register row_ptr_H and the column address index register col_ptr_H to null pointers, and set the initial values ​​of the row element number register num_row_H, the column element number register num_col_H, the data element register data_H and the data row and column label value register data_col_row_H to 0; Where M is the Jacobian matrix; The standard orthogonal basis calculation unit is used to calculate the Jacobian matrix M and the standard orthogonal basis V m The jth element in calculates the basis vector w under dimension j, and for all the j orthogonal bases that have been generated, calculates the correction vector and performs orthogonalization calculation to complete the standard orthogonal basis V m Orthogonalization calculation of ; The approximate optimal solution calculation unit is used to calculate the standard orthogonal basis V after orthogonalization. m Calculate the approximate optimal solution x of Mx = b m , where b is the power imbalance matrix; A bidirectional linked list matrix compression storage module is used for matrix operation during iterative calculation of the Krylov subspace dimensionality reduction solution module, wherein the matrix operation is performed based on matrix block operation; the bidirectional linked list matrix compression storage module includes: A matrix block unit, used for dividing the Jacobian matrix into blocks by rows to obtain a plurality of matrix blocks; A row-block parallel computing unit is used to obtain a product result component based on the matrix block, and the product result component is a result obtained by dividing the product result of the Jacobian matrix and the column vector into row blocks.

2. The Krylov subspace dimension reduction power flow calculation and solution system according to claim 1 is characterized in that: The matrix blocking unit completes the Jacobian matrix blocking in the following manner: Accumulate the number of non-zero elements in the Jacobian matrix row by row to obtain the accumulated value of non-zero elements; Determine whether the current non-zero element accumulation value is greater than s+1 times the number of average matrix division elements, where s is the row number of the Jacobi matrix row-by-row division point; if the current non-zero element accumulation value is greater than s+1 times the number of average matrix division elements, then take the current row number i as a matrix block division point, and let the Jacobi matrix row-by-row division point row number s be incremented by one; Repeat the previous step until the current row number is the last row of the Jacobian matrix, and set the elements of the remaining Jacobian matrix as a matrix block.

3. The Krylov subspace dimension reduction power flow calculation and solution system according to claim 1, characterized in that: The row-block parallel computing unit obtains the product result component based on the matrix block in the following manner: Calculate the number of rows of the ith matrix block, and find out the computational components of other block parallel computing units that the jth row of the ith matrix block needs to rely on; If the number of computational components of other block parallel computing units that the jth row of the i-th matrix block needs to rely on is zero, Directly multiply the j-th row of the i-th matrix block by the column vector to obtain the product result component; If the number of calculation components of other block parallel computing units that the j-th row of the i-th matrix block needs to rely on is greater than zero, then find the numbers of the calculation components of the other dependent block parallel computing units, multiply the j-th row of the i-th matrix block with the column vector to obtain a first result, then multiply the j-th row of the i-th matrix block with the calculation components of the other dependent block parallel computing units and accumulate to obtain a second result, sum the first result and the second result to obtain the product result component.

4. The Krylov subspace dimension reduction power flow calculation and solution system according to claim 1, characterized in that: The standard orthogonal basis calculation unit completes the standard orthogonal basis V in the following way m Orthogonal calculation of : The Jacobian matrix M and the orthonormal basis V m Multiply the jth element in to get the basis vector w under dimension j; determine the basis vector w pointer register address ptr_id_w and the generated orthogonal basis pointer vector address ptr_id_V i The relationship between the size of the register and the number of non-zero elements; When the base vector w pointer register address ptr_id_w is less than or equal to the number of non-zero elements in the base vector w register, or the orthogonal base pointer vector address ptr_id_V has been generated i When it is less than or equal to the number of non-zero elements in the orthogonal basis register, record the basis vector w and the orthogonal basis V i The element position number in , and the data value corresponding to the position number; If the element positions of the basis vector w are numbered with the orthogonal basis V i If the element position numbers in are the same, the product operation of the corresponding elements is performed and accumulated and stored in the inner product data_H_ij; If the element position number of the basis vector w is greater than the orthogonal basis V i The element position number in and the orthogonal base pointer vector address ptr_id_V that has been generated i When the number of non-zero elements in the orthogonal base register is less than the number of non-zero elements in the orthogonal base register, the generated orthogonal base pointer vector address ptr_id_V i plus one; If the element position number of the basis vector w is less than the orthogonal basis V i If the element position number in , and the basic vector w pointer register address ptr_id_w is less than the number of non-zero elements in the basic vector w register, the basic vector w pointer register address ptr_id_w is increased by one; Subtract the value in the inner product memory data_H_ij from the basis vector w and the orthogonal basis V i The product of is used as the new basis vector; Repeat the above steps until the accumulated value of the inner product data_H_ij is calculated; When the index of the calculated orthogonal basis vector is greater than or equal to 2, find the previous element in this row of the inner product data_H_ij, and record the address of the first element of the i-th row as the linked list head pointer register row_ptr_head_H[i]. If the address of the first element of the i-th row is empty, directly point the row linked list head pointer register row_ptr_head_H[i] to the address pointer register of the first position of the inner product data_H_ij; if the address of the first element of the i-th row is not empty, perform the following steps until the previous element position number of the inner product data_H_ij is found; Take out the next element address next_col_id of the first element of the i-th row in the row pointer register row_ptr_H, and determine whether the next element address next_col_id is empty. If the next element address next_col_id is empty, the previous element position number of the inner product data_H_ij has been found, and the element at the next element address next_col_id points to the address pointer register of the first position of the inner product data_H_ij, and exit the process; If the next element address next_col_id is not empty, the first element address of the i-th row will be updated to the next element address next_col_id.

5. The Krylov subspace dimension reduction power flow calculation and solution system according to claim 4 is characterized in that: The accumulated value of the inner product data_H_ij is stored as follows: The address pointer register pointing to the first position of data_H is recorded as ptr_data_H, and the initial value of the address pointer register ptr_data_H is set to zero; Determine whether the inner product data_H_ij is zero. If the inner product data_H_ij is not zero, store the inner product data_H_ij to the address pointed to by the address pointer register ptr_data_H, and store the row and column label values ​​of the inner product data_H_ij at the same time; Determine whether the column link list head pointer register col_ptr_head_H and the row link list head pointer register row_ptr_head_H at the position (i, j) are empty; if the column link list head pointer register col_ptr_head_H[j] of the jth column is not empty, then point the current column link list head pointer to the address pointed to by the address pointer register ptr_data_H; if the row link list head pointer register row_ptr_head_H[i] of the i-th row is not empty, then point the current row link list head pointer to the address pointed to by the address pointer register ptr_data_H; The pointer of the non-zero element number register of the row and column where (i, j) is located is incremented by one, and the pointer of the column address index register col_ptr_H is incremented by one.

6. The Krylov subspace dimension reduction power flow calculation and solution system according to claim 1, characterized in that: The approximate optimal solution calculation unit is provided by x m =x0+V m ·y m Calculate the approximate optimal solution, where x m is a near optimal solution, It represents the inverse matrix of the m-dimensional Heisenberg matrix H. unit_e is an m-dimensional column vector whose first element is 1.

7. A power flow calculation system, characterized in that: It comprises a Krylov subspace dimensionality reduction power flow calculation and solution system as described in any one of claims 1 to 6, a judgment module and a node voltage and power calculation module, wherein the judgment module is used to judge whether the change of the approximate optimal solution calculated by the Krylov subspace dimensionality reduction power flow calculation and solution system within a preset time is less than a preset error; the node voltage and power calculation module is used to calculate the network node flow distribution according to the approximate optimal solution when the change of the approximate optimal solution is less than a preset error.

8. The power flow calculation system according to claim 7, characterized in that: It also includes a Jacobian matrix calculation module and a power imbalance calculation module; the Jacobian matrix calculation module is used to generate a Jacobian matrix; the power imbalance calculation module is used to generate a power imbalance matrix.

Images