CN105391057B  A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates  Google Patents
A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates Download PDFInfo
 Publication number
 CN105391057B CN105391057B CN201510809809.3A CN201510809809A CN105391057B CN 105391057 B CN105391057 B CN 105391057B CN 201510809809 A CN201510809809 A CN 201510809809A CN 105391057 B CN105391057 B CN 105391057B
 Authority
 CN
 China
 Prior art keywords
 msub
 mrow
 matrix
 nodes
 node
 Prior art date
 Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
 Expired  Fee Related
Links
 239000011159 matrix material Substances 0.000 claims abstract description 83
 238000004364 calculation method Methods 0.000 claims abstract description 32
 230000001174 ascending Effects 0.000 claims description 15
 230000000875 corresponding Effects 0.000 claims description 9
 238000003860 storage Methods 0.000 claims description 9
 238000002347 injection Methods 0.000 claims description 6
 239000007924 injection Substances 0.000 claims description 6
 210000001503 Joints Anatomy 0.000 claims description 2
 238000006073 displacement reaction Methods 0.000 claims description 2
 238000000034 method Methods 0.000 abstract description 5
 238000009826 distribution Methods 0.000 description 4
 230000001133 acceleration Effects 0.000 description 3
 230000000694 effects Effects 0.000 description 2
 238000004519 manufacturing process Methods 0.000 description 2
 238000011030 bottleneck Methods 0.000 description 1
 238000004422 calculation algorithm Methods 0.000 description 1
 238000005094 computer simulation Methods 0.000 description 1
 230000003247 decreasing Effects 0.000 description 1
 230000005611 electricity Effects 0.000 description 1
 230000005520 electrodynamics Effects 0.000 description 1
 238000005516 engineering process Methods 0.000 description 1
 238000011156 evaluation Methods 0.000 description 1
 238000005457 optimization Methods 0.000 description 1
 230000035945 sensitivity Effects 0.000 description 1
 238000004088 simulation Methods 0.000 description 1
 241000894007 species Species 0.000 description 1
 230000003068 static Effects 0.000 description 1
 230000002123 temporal effect Effects 0.000 description 1
 230000001550 time effect Effects 0.000 description 1
Classifications

 H—ELECTRICITY
 H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
 H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
 H02J3/00—Circuit arrangements for ac mains or ac distribution networks

 H—ELECTRICITY
 H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
 H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
 H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
 H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
Abstract
The invention discloses a kind of GPU threaded design methods that direction of energy Jacobi battle array calculates, it is characterised in that：Methods described includes：Input electric network data, pretreatment node admittance battle array Y；CPU distinguishes calculated submatrix nonzero element and node admittance battle array Y nonzero element position mapping tables；CPU distinguishes calculated submatrix nonzero element and Jacobian matrix nonzero element position mapping table；GPU interior joint injecting power kernel functions S calculates all node injecting powers；Jacobi submatrix calculates kernel function and distinguishes nonzero element in calculated submatrix and be stored in Jacobian matrix in GPU.The present invention is calculated in Jacobian matrix nonzero entry calculating process by submatrix to judge the branched structure of the affiliated submatrix region process of element when avoiding and directly being calculated using single kernel function, lifts execution efficiency；And by centralized calculation there is nonzero element in the submatrix of identical calculations formula to solve the unbalanced problem of threads load, lift parallel efficiency calculation.
Description
Technical field
The invention belongs to High performance computing in power system application field, more particularly to a kind of direction of energy Jacobi battle array to calculate
GPU threaded design methods.
Background technology
Load flow calculation is most widely used, most basic and most important a kind of electric computing in power system.In power train
In the research of the method for operation of uniting and programme, it is required for carrying out Load flow calculation to compare the method for operation or plan power supply plan
Feasibility, reliability and economy.Meanwhile for the running status of realtime electric power monitoring system, it is also desirable to carry out a large amount of and fast
The Load flow calculation of speed.Therefore, in the method for operation of programming and planning and schedule system, using offline Load flow calculation；In electricity
In the realtime monitoring of Force system running status, then calculated using online power flow.
And in actual production process, offline trend and online power flow calculating all have this to compare the calculating speed of trend
High requirement.It is being related to planning and designing and complexity situations such as in the offline trend for arranging the method for operation, landing scheme because of equipment, is needing
Want the species of simulation run more, Load flow calculation amount is big, and single Load flow calculation time effects integrally emulate duration；And in power system
The online power flow carried out in operation is calculated to calculating temporal sensitivity height, it is necessary to provide calculation of tidal current in real time, is such as being envisioned
In accident, the equipment Load flow calculation out of service to the influence of static security, system needs to calculate trend under a large amount of forecast accidents
Distribution, and the method for operation Adjusted Option of anticipation is made in real time.
In traditional NewtonLaphson method Load flow calculation, the calculating time of Jacobian matrix accounts for trend and amounts to evaluation time
30%, the calculating speed of Jacobian matrix influences the overall performance of program.And slowing down with the lifting of CPU calculating speeds, existing rank
The single Load flow calculation calculating time of section has reached a bottleneck.The accelerated method of Load flow calculation, which is concentrated mainly on, at present makes
More trends are carried out with coarseness acceleration with cluster and multiplecore server, it is actual into being ground to single trend internal arithmetic acceleration in production
Study carefully less.
GPU is a kind of manycore parallel processor, will be considerably beyond CPU in the quantity of processing unit.GPU traditionally is only
Responsible figure renders, and CPU has all been given in most processing.Method battle array is a kind of multinuclear to present GPU, multithreading, is had
There are powerful calculating ability and high bandwidth of memory, programmable processor.Under universal computer model, associations of the GPU as CPU
Processor works, and is decomposed by task reasonable distribution and completes highperformance calculation.
Jacobian matrix, which calculates, has concurrency, and the calculating of each of which nonzero element is separate, does not rely on pass
System, it can naturally be handled by parallel calculating, be adapted to GPU to accelerate.Therefore can be fast by rational scheduling between CPU and GPU
Speed completes the calculating of Jacobian matrix, and the method that domestic and foreign scholars have begun to carry out GPU Jacobi acceleration is studied,
But without deep optimization threaded design, computational threads design is studied from the distribution of amount of calculation merely, to thread calculating side
Formula, data directory mode are not furtherd investigate, and program can not be made to give full play to GPU advantage.
It would therefore be highly desirable to solve the above problems.
The content of the invention
Goal of the invention：In view of the shortcomings of the prior art, when Jacobian matrix calculating can be greatly decreased in present invention offer one kind
Between and can be lifted Load flow calculation speed a kind of direction of energy Jacobi battle array calculate GPU threaded design methods.
Technical scheme：The present invention proposes a kind of a kind of tide of electric power based on GPU that computing resource is distributed using mapping relations
Flow the GPU threaded design methods that Jacobian matrix calculates.
Load flow calculation：Electrodynamic noun, refer in given power system network topology, component parameters and generating, load parameter
Under the conditions of, calculate the distribution of active power, reactive power and voltage in power network.
Parallel computation：It is a kind of algorithm that once can perform multiple instruction, it is therefore an objective to improve and calculate relative to serial arithmetic
Speed, and by expanding problem solving scale, solve largescale and complicated computational problem.
GPU：Graphics processor (English：GraphicsProcessingUnit, abbreviation：GPU).
Admittance matrix：Established based on the Equivalent admittance of system element, each node voltage of description electric power networks and
The matrix of relation between Injection Current.
The invention discloses a kind of GPU threaded design methods that direction of energy Jacobi battle array calculates, methods described includes：
(1) electric network data, pretreatment node admittance battle array Y are inputted；
(2) CPU calculates the sub matrix nonzero elements of H, N, J, L tetra respectively and node admittance battle array Y nonzero element positions map
Relation table M_{H2Y}、M_{N2Y}、M_{J2Y}、M_{L2Y}, circular comprises the following steps：
(2.1) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set that Hmatrix calculates
H_{Y}, generation inquiry table M_{H2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to active reactive PQ nodes or active voltage PV node
When, remember that i=pqpv [x], wherein pqpv are PQ, the index of PV node, x is call number；Row number j is scanned again, when j belongs to PQ or PV
During node, remember j=pqpv [y], choose the element to be added to set H_{Y}；
B) element is often added to set H_{Y}When, set of records ends H_{Y}CurrentElement quantity k, and the element is matrix Y's
Position h in COO forms_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ, PV node_{k}；Set H_{Y}Middle element representation
For { k, h_{k}, x_{k}, y_{k}}；
C) k → h is generated_{k}Mapping table M_{H2Y}；
(2.2) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set that N matrix calculates
N_{Y}, generation mapping table M_{N2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ or PV node, remembers i=pqpv [x], then sweep
Row number j is retouched, when j belongs to PQ nodes, remembers that j=pq [y], wherein pq are the index of PQ nodes, y is call number, chooses the element
It is added to set N_{Y}；
B) element is often added to set N_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is matrix Y's
Position n in COO forms_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ nodes_{k}；Set N_{Y}Middle element representation is
{ k, n_{k}, x_{k}, y_{k}}；
C) k → n is generated_{k}Mapping table M_{N2Y}；
(2.3) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set corresponding to the element of J matrix computations
M_{Y}, generation mapping table M_{J2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scan row number
J, when j belongs to PQ or PV node, remember j=pqpv [y], choose the element to be added to set M_{Y}；
B) element is often added to set J_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is matrix Y's
Position m in COO forms_{k}；Record i position x in PQ nodes_{k}, j position y in PQ, PV node_{k}；Set J_{Y}Middle element representation is
{ k, j_{k}, x_{k}, y_{k}}；
C) k → j is generated_{k}Mapping table M_{J2Y}；
(2.4) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set corresponding to the element of L matrix computations
L_{Y}, generation inquiry table M_{L2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scan row number
J, when j belongs to PQ nodes, remember j=pq [y], choose the element to be added to set L_{Y}；
B) element is often added to set L_{Y}When, set of records ends L_{Y}CurrentElement quantity k, and the element is matrix Y's
Position l in COO forms_{k}；Record i position x in PQ nodes_{k}, j position y in PQ nodes_{k}；Set L_{Y}Middle element representation for k,
l_{k}, x_{k}, y_{k}}；
C) k → l is generated_{k}Mapping table M_{L2Y}；
(3) CPU calculates the sub matrix nonzero elements of H, N, J, L tetra respectively and the mapping of Jacobian matrix nonzero element position is closed
It is table M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}, circular comprises the following steps：
(3.1) according to the H obtained in step 2_{Y}、N_{Y}、J_{Y}、L_{Y}In x_{k}, y_{k}Obtain matrix H, N, J, L sparse format
(3.2) definition set Cmp { x, y, NO, zone }, wherein x are Jacobian matrix line number, y is Jacobian matrix row
Number, NO be element position in H, N, J, L Matrix C OO forms, equivalent to the k in step 2, in zone H, N, J, L subregion
The identifier of element, PQ, PV node number are npqpv；
(3.3) element in matrix H, N, J, L is inserted in Cmp, wherein element zone is set to 1 in H battle arrays；Element y values in N battle arrays
Npqpv, zone is added to be set to 2；Element x value adds npqpv, zone to be set to 3 in J battle arrays；Element x value adds npqpv in L battle arrays, and y values add
Npqpv, zone are set to 4；
(3.4) set Cmp is reordered by x ascending orders, y ascending orders, obtains the Jacobi Matrix C OO forms of row ascending order；
In set Cmp add variable num, represent Jacobian matrix in position of the element in standard COO forms, and in order assign 1,
2、3…；
(3.5) to Cmp { x, y, NO, zone, num } by zone ascending orders, NO ascending sorts, the num in each zone regions
Sequence is inquiry table M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}；
(4) GPU interior joints injecting power kernel function S calculates all node injecting powers；
(5) Jacobi submatrix calculates kernel function H, M, N, L and calculates nonzero in tetra submatrixs of H, N, J, L respectively in GPU
Element is simultaneously stored in Jacobian matrix.
Wherein, the nodes of power network are n in the step (1), the node admittance battle array Y { i, j, G, B } of power network with line number i,
Row number j, conductance G, susceptance B Coordinate sparse formats storage, by relative skew of the node admittance battle array each row of data in Y
Position is stored in line displacement vector R, and definition node voltage phase angle vector is θ, and voltage magnitude vector is V.
Wherein, the k threads that the calculation of step (4) the interior joint injecting power is kernel function S calculate k nodes
The injection of active and reactive power, calculation formula is：
Wherein,
K is parallel thread number in GPU Programmings；
P_{k}Injected for the active power of k nodes；
Q_{k}Injected for the reactive power of k nodes；
k_{b}Indexed for the starting of admittance battle array Y row k nonzero elements, k_{b}=R [k]；
k_{e}Indexed for the termination of admittance battle array Y row k nonzero elements, k_{e}=R [k+1] 1；
J, G, B are the data of lth of position storage in admittance battle array Y；
V_{k}For the voltage magnitude of k nodes；
V_{j}For the voltage magnitude of j nodes；
θ_{k}For the voltage phase angle of k nodes；
θ_{j}For the voltage phase angle of j nodes.
Wherein, in the step (5), the computational methods of nonzero element are in H submatrixs：
According to the H submatrixs nonzero element and node admittance battle array Y mapping table M obtained in step (2)_{H2Y}, H core letters
When several t threads calculate, inquiry table M_{H2Y}T positions take out Y battle arrays call number k_{H2Y}, take out Y [k_{H2Y}] corresponding in element
Data { i, j, G, B } are calculated, and the calculation formula of tth of nonzero element of Jacobian matrix neutron matrix H is：
Wherein,
T is parallel thread number in GPU Programmings；H_{t}For the nonzero element of tth of position in H submatrix COO forms；
I, j, G, B are kth in admittance battle array Y_{H2Y}The data of individual position storage；
Q_{i}For the injection reactive power of node i；
According still further to the H submatrixs nonzero element and matrix J a mapping relations M obtained in step (3)_{H2Ja}, inquire about M_{H2Ja}T
Number position take out H_{t}Storage call number k of the result of calculation in Ja matrixes_{H2Ja}, then by H_{t}Result of calculation deposit Ja
[k_{H2Ja}]。
Beneficial effect：Compared with the prior art, beneficial effects of the present invention are：First, the present invention uses GPU processors to refined
Gram accelerated than matrix computations, compared to the speedup ratio that CPU obtains 16 times；Secondly, threaded design method of the invention is to GPU
The performance of program is optimized, and is calculated in Jacobian matrix nonzero entry calculating process by submatrix to avoid using single
The branched structure of the affiliated submatrix region process of element, lifting GPU branches execution efficiency are judged when kernel function directly calculates；And lead to
Crossing centralized calculation, there is nonzero element in the submatrix of identical calculations formula to solve the unbalanced problem of threads load, and lifting is parallel
Computational efficiency.
Brief description of the drawings
Fig. 1 is the flow chart for the GPU threaded design methods that a kind of direction of energy Jacobi battle array of the present invention calculates
Embodiment
As shown in figure 1, the present invention proposes a kind of a kind of tide of electric power based on GPU that computing resource is distributed using mapping relations
The GPU threaded design methods that Jacobian matrix calculates are flowed, the method includes the steps of：
If the nodes of power network are n, the node admittance battle array Y { i, j, G, B } of power network is with line number, row number, conductance, susceptance
COO sparse formats store, and index in COO forms are converted into CompressedSparseRow (CSR) row compressed format, row is inclined
Shifting array vector is R, and node voltage phase angle vector is θ, and voltage magnitude vector is V；
Step 1：Generate inquiry table M_{H2Y}、M_{N2Y}、M_{J2Y}、M_{L2Y}
1) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set H that Hmatrix calculates_{Y}, it is raw
Into inquiry table M_{H2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to active reactive PQ nodes or active voltage PV node
When, remember that i=pqpv [x], wherein pqpv are PQ, the index of PV node, x is call number；Row number j is scanned again, when j belongs to PQ or PV
During node, remember j=pqpv [y], choose the element to be added to set H_{Y}；
B) element is often added to set H_{Y}When, set of records ends H_{Y}CurrentElement quantity k, and the element is matrix Y's
Position h in COO forms_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ, PV node_{k}；Set H_{Y}Middle element representation
For { k, h_{k}, x_{k}, y_{k}}；
C) k → h is generated_{k}Mapping table M_{H2Y}；
2) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set N that N matrix calculates_{Y}, it is raw
Into mapping table M_{N2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ or PV node, remembers i=pqpv [x], then sweep
Row number j is retouched, when j belongs to PQ nodes, remembers that j=pq [y], wherein pq are the index of PQ nodes, y is call number, chooses the element
It is added to set N_{Y}；
B) element is often added to set N_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is matrix Y's
Position n in COO forms_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ nodes_{k}；Set N_{Y}Middle element representation is
{ k, n_{k}, x_{k}, y_{k}}；
C) k → n is generated_{k}Mapping table M_{N2Y}；
3) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set M corresponding to the element of J matrix computations_{Y}, it is raw
Into mapping table M_{J2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scan row number
J, when j belongs to PQ or PV node, remember j=pqpv [y], choose the element to be added to set M_{Y}；
B) element is often added to set J_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is matrix Y's
Position m in COO forms_{k}；Record i position x in PQ nodes_{k}, j position y in PQ, PV node_{k}；Set J_{Y}Middle element representation is
{ k, j_{k}, x_{k}, y_{k}}；
C) k → j is generated_{k}Mapping table M_{J2Y}；
4) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set L corresponding to the element of L matrix computations_{Y}, it is raw
Into inquiry table M_{L2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scan row number
J, when j belongs to PQ nodes, remember j=pq [y], choose the element to be added to set L_{Y}；
B) element is often added to set L_{Y}When, set of records ends L_{Y}CurrentElement quantity k, and the element is matrix Y's
Position l in COO forms_{k}；Record i position x in PQ nodes_{k}, j position y in PQ nodes_{k}；Set L_{Y}Middle element representation for k,
l_{k}, x_{k}, y_{k}}；
C) k → l is generated_{k}Mapping table M_{L2Y}；
Step 2：Generate inquiry table M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}
1) according to the H obtained in step 1_{Y}、N_{Y}、J_{Y}、L_{Y}In x_{k}, y_{k}Obtain matrix H, N, J, L sparse format
2) definition set Cmp { x, y, NO, zone }, wherein x are Jacobian matrix line number, y is Jacobian matrix row number, NO
For element in H, N, J, L Matrix C OO forms position, equivalent to the k in step 1, zone H, N, J, L subregion interior element
Identifier, PQ, PV node number are npqpv；
3) element in matrix H, N, J, L is inserted in Cmp, wherein element zone is set to 1 in H battle arrays；Element y values add in N battle arrays
Npqpv, zone are set to 2；Element x value adds npqpv, zone to be set to 3 in J battle arrays；Element x value adds npqpv, y values plus npqpv in L battle arrays,
Zone is set to 4；
4) set Cmp is reordered by x ascending orders, y ascending orders, obtains the Jacobi Matrix C OO forms of row ascending order；Collecting
Close Cmp in add variable num, represent Jacobian matrix in position of the element in standard COO forms, and in order assign 1,2,
3…；
5) to Cmp { x, y, NO, zone, num } by zone ascending orders, NO ascending sorts, the num sequences in each zone regions
As inquiry table M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}。
Step 3：Calculate node injecting power
S kernel functions calculate the active reactive injection of all nodes, and setting kernel function startup parameter is<<<n/256+1,256>
>>, that is, start onedimensional piece of block quantity [n/256]+1, open 256 thread thread in each block, amount to n
threads.Required data are passed in GPU in advance, the k threads of S kernel functions calculate the active and reactive power note of k nodes
Enter, thread calculation formula is：
Wherein,
k_{b}Indexed for the starting of admittance battle array Y row k nonzero elements, k_{b}=R [k]；
k_{e}Indexed for the termination of admittance battle array Y row k nonzero elements, k_{e}=R [k+1] 1；
J, G, B are the data of lth of position storage in admittance battle array Y.
Step 4：Calculate nonzero element value in Jacobian matrix
Jacobian matrix is calculated and is divided into four parts, respectively calculated submatrix H, M, J, L, each subfunction in calculating process
Thread by mapping directly read calculate data calculated.
By taking the calculating of submatrix H nonzero entries value as an example,
H_{1}Kernel function is responsible for calculated submatrix H nonzero element, sets its startup parameter to be<<<H_{num}/256+1,256>>
>, that is, start onedimensional block quantity [H_{num}/ 256]+1, each block is interior to open 256 threads, amounts to H_{num}It is individual
Threads, wherein H_{num}For submatrix H nonzero entry number.Required data are passed in GPU in advance, H_{1}The t lines of kernel function
When journey calculates, inquiry table M_{H2Y}T positions take out Y battle arrays call number k_{H2Y}, take out Y [k_{H2Y}] corresponding data in element i, j,
G, B } calculated, tth of nonzero element that content is Jacobian matrix neutron matrix H is calculated, calculation formula is：
Wherein,
I, j, G, B are kth in admittance battle array Y_{H2Y}The data of individual position storage；
Q_{i}For the injection reactive power of node i.
Establish H_{2}Global the thread number t and matrix J acobi of kernel function mapping relations M_{H2Ja}, H_{2}The t threads of kernel function
After calculating result, M is inquired about_{H2Ja}T positions take out storage call number k of the result of calculation in Jacobi matrixes_{H2Ja}, then
By result of calculation deposit Jacobi [k_{H2Ja}]。
The calculating of other submatrixs is similar therewith.
Claims (4)
1. a kind of GPU threaded design methods that direction of energy Jacobi battle array calculates, it is characterised in that：Methods described includes：
(1) electric network data, pretreatment node admittance battle array Y are inputted；
(2) CPU calculates the sub matrix nonzero elements of H, N, J, L tetra and node admittance battle array Y nonzero element position mapping relations respectively
Table M_{H2Y}、M_{N2Y}、M_{J2Y}、M_{L2Y}, circular comprises the following steps：
(2.1) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set H that Hmatrix calculates_{Y}, generation
Inquiry table M_{H2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to active reactive PQ nodes or active voltage PV node,
Remember that i=pqpv [x], wherein pqpv are PQ, the index of PV node, x is call number；Row number j is scanned again, when j belongs to PQ or PV sections
During point, remember j=pqpv [y], choose the element to be added to set H_{Y}；
B) element is often added to set H_{Y}When, set of records ends H_{Y}CurrentElement quantity k, and the element is in matrix Y COO
Position h in form_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ, PV node_{k}；Set H_{Y}Middle element representation is
{ k, h_{k}, x_{k}, y_{k}}；
C) k → h is generated_{k}Mapping table M_{H2Y}；
(2.2) the node admittance battle array Y of COO forms is scanned, chooses and wherein corresponds to the element addition set N that N matrix calculates_{Y}, generation
Mapping table M_{N2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ or PV node, remembers i=pqpv [x], then scan columns
Number j, when j belongs to PQ nodes, remember that j=pq [y], wherein pq are the index of PQ nodes, y is call number, chooses the element to add
To set N_{Y}；
B) element is often added to set N_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is in matrix Y COO
Position n in form_{k}；Record i position x in PQ, PV node_{k}, j position y in PQ nodes_{k}；Set N_{Y}Middle element representation for k,
n_{k}, x_{k}, y_{k}}；
C) k → n is generated_{k}Mapping table M_{N2Y}；
(2.3) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set M corresponding to the element of J matrix computations_{Y}, generation
Mapping table M_{J2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scans row number j, when
When j belongs to PQ or PV node, remember j=pqpv [y], choose the element to be added to set M_{Y}；
B) element is often added to set J_{Y}When, set of records ends N_{Y}CurrentElement quantity k, and the element is in matrix Y COO
Position m in form_{k}；Record i position x in PQ nodes_{k}, j position y in PQ, PV node_{k}；Set J_{Y}Middle element representation for k,
j_{k}, x_{k}, y_{k}}；
C) k → j is generated_{k}Mapping table M_{J2Y}；
(2.4) the node admittance battle array Y of COO forms is scanned, chooses and wherein adds set L corresponding to the element of L matrix computations_{Y}, generation
Inquiry table M_{L2Y}；
A) the line number i of each element in COO forms is scanned, when it belongs to PQ nodes, remembers i=pq [x], then scans row number j, when
When j belongs to PQ nodes, remember j=pq [y], choose the element to be added to set L_{Y}；
B) element is often added to set L_{Y}When, set of records ends L_{Y}CurrentElement quantity k, and the element is in matrix Y COO
Position l in form_{k}；Record i position x in PQ nodes_{k}, j position y in PQ nodes_{k}；Set L_{Y}Middle element representation is { k, l_{k},
x_{k}, y_{k}}；
C) k → l is generated_{k}Mapping table M_{L2Y}；
(3) CPU calculates the sub matrix nonzero elements of H, N, J, L tetra and Jacobian matrix nonzero element position mapping table respectively
M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}, circular comprises the following steps：
(3.1) according to the H obtained in step 2_{Y}、N_{Y}、J_{Y}、L_{Y}In x_{k}, y_{k}Obtain matrix H, N, J, L sparse format
(3.2) definition set Cmp { x, y, NO, zone }, wherein x are Jacobian matrix line number, y is Jacobian matrix row number, NO
For element in H, N, J, L Matrix C OO forms position, equivalent to the k in step 2, zone H, N, J, L subregion interior element
Identifier, PQ, PV node number are npqpv；
(3.3) element in matrix H, N, J, L is inserted in Cmp, wherein element zone is set to 1 in H battle arrays；Element y values add in N battle arrays
Npqpv, zone are set to 2；Element x value adds npqpv, zone to be set to 3 in J battle arrays；Element x value adds npqpv, y values plus npqpv in L battle arrays,
Zone is set to 4；
(3.4) set Cmp is reordered by x ascending orders, y ascending orders, obtains the Jacobi Matrix C OO forms of row ascending order；Collecting
Close Cmp in add variable num, represent Jacobian matrix in position of the element in standard COO forms, and in order assign 1,2,
3…；
(3.5) to Cmp { x, y, NO, zone, num } by zone ascending orders, NO ascending sorts, the num sequences in each zone regions
As inquiry table M_{H2Ja}、M_{N2Ja}、M_{J2Ja}、M_{L2Ja}；
(4) GPU interior joints injecting power kernel function S calculates all node injecting powers；
(5) Jacobi submatrix calculates kernel function H, M, N, L and calculates nonzero element in tetra submatrixs of H, N, J, L respectively in GPU
And it is stored in Jacobian matrix.
2. the GPU threaded design methods that direction of energy Jacobi battle array according to claim 1 calculates, it is characterised in that：Institute
The nodes for stating power network in step (1) are n, and the node admittance battle array Y { i, j, G, B } of power network is with line number i, row number j, conductance G, susceptance
B Coordinate sparse formats storage, line displacement is stored in by relative deviation post of the node admittance battle array each row of data in Y
In vectorial R, definition node voltage phase angle vector is θ, and voltage magnitude vector is V.
3. the GPU threaded design methods that direction of energy Jacobi battle array according to claim 1 calculates, it is characterised in that：
The k threads that the calculation of step (4) the interior joint injecting power is kernel function S calculate k nodes active and
Reactive power is injected, and calculation formula is：
<mrow>
<msub>
<mi>P</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msub>
<mi>V</mi>
<mi>k</mi>
</msub>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>l</mi>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>b</mi>
</msub>
</mrow>
<msub>
<mi>k</mi>
<mi>e</mi>
</msub>
</munderover>
<msub>
<mi>V</mi>
<mi>j</mi>
</msub>
<mo>&lsqb;</mo>
<mi>G</mi>
<mi> </mi>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mi>k</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>B</mi>
<mi> </mi>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mi>k</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<msub>
<mi>Q</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msub>
<mi>V</mi>
<mi>k</mi>
</msub>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>l</mi>
<mo>=</mo>
<msub>
<mi>k</mi>
<mi>b</mi>
</msub>
</mrow>
<msub>
<mi>k</mi>
<mi>e</mi>
</msub>
</munderover>
<msub>
<mi>V</mi>
<mi>j</mi>
</msub>
<mo>&lsqb;</mo>
<mi>G</mi>
<mi> </mi>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mi>k</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo></mo>
<mi>B</mi>
<mi> </mi>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&theta;</mi>
<mi>k</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
Wherein,
K is S kernel functions parallel thread number in GPU Programmings；
P_{k}Injected for the active power of k nodes；
Q_{k}Injected for the reactive power of k nodes；
k_{b}Indexed for the starting of admittance battle array Y row k nonzero elements, k_{b}=R [k]；
k_{e}Indexed for the termination of admittance battle array Y row k nonzero elements, k_{e}=R [k+1] 1；
J, G, B are the data of lth of position storage in admittance battle array Y；
V_{k}For the voltage magnitude of k nodes；
V_{j}For the voltage magnitude of j nodes；
θ_{k}For the voltage phase angle of k nodes；
θ_{j}For the voltage phase angle of j nodes.
4. the GPU threaded design methods that direction of energy Jacobi battle array according to claim 1 calculates, it is characterised in that：
In the step (5), the computational methods of nonzero element are in H submatrixs：
According to the H submatrixs nonzero element and node admittance battle array Y mapping table M obtained in step (2)_{H2Y}, the t of H kernel functions
When number thread calculates, inquiry table M_{H2Y}T positions take out Y battle arrays call number k_{H2Y}, take out Y [k_{H2Y}] corresponding data in element
{ i, j, G, B } is calculated, and the calculation formula of tth of nonzero element of Jacobian matrix neutron matrix H is：
<mrow>
<msub>
<mi>H</mi>
<mi>t</mi>
</msub>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mo></mo>
<msub>
<mi>V</mi>
<mi>i</mi>
</msub>
<msub>
<mi>V</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<mi>G</mi>
<mi> </mi>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>&theta;</mi>
<mi>i</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mo></mo>
<mi>B</mi>
<mi> </mi>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>&theta;</mi>
<mi>i</mi>
</msub>
<mo></mo>
<msub>
<mi>&theta;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>i</mi>
<mo>&NotEqual;</mo>
<mi>j</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>V</mi>
<mi>i</mi>
<mn>2</mn>
</msubsup>
<mi>B</mi>
<mo>+</mo>
<msub>
<mi>Q</mi>
<mi>i</mi>
</msub>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>j</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
Wherein,
T is H kernel functions parallel thread number in GPU Programmings；
H_{t}For the nonzero element of tth of position in H submatrix COO forms；
I, j, G, B are kth in admittance battle array Y_{H2Y}The data of individual position storage；
Q_{i}For the injection reactive power of node i；
V_{i}For the voltage magnitude of i nodes；
V_{j}For the voltage magnitude of j nodes；
θ_{i}For the voltage phase angle of i nodes；
θ_{j}For the voltage phase angle of j nodes；
According still further to the H submatrixs nonzero element and matrix J a mapping relations M obtained in step (3)_{H2Ja}, inquire about M_{H2Ja}T positions
Put and take out H_{t}Storage call number k of the result of calculation in Ja matrixes_{H2Ja}, then by H_{t}Result of calculation deposit Ja [k_{H2Ja}]。
Priority Applications (1)
Application Number  Priority Date  Filing Date  Title 

CN201510809809.3A CN105391057B (en)  20151120  20151120  A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

CN201510809809.3A CN105391057B (en)  20151120  20151120  A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates 
Publications (2)
Publication Number  Publication Date 

CN105391057A CN105391057A (en)  20160309 
CN105391057B true CN105391057B (en)  20171114 
Family
ID=55423022
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

CN201510809809.3A Expired  Fee Related CN105391057B (en)  20151120  20151120  A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates 
Country Status (1)
Country  Link 

CN (1)  CN105391057B (en) 
Families Citing this family (6)
Publication number  Priority date  Publication date  Assignee  Title 

CN106157176B (en) *  20160726  20190712  东南大学  A kind of LU decomposition method for the direction of energy Jacobian matrix that GPU accelerates 
CN106684863B (en) *  20161229  20180227  华中科技大学  A kind of discrimination method of power distribution network bus admittance matrix 
CN107368455A (en) *  20170622  20171121  东南大学  Trigonometric equation group back substitution method on the direction of energy that a kind of GPU accelerates 
CN108879691B (en) *  20180621  20200904  清华大学  Largescale continuous power flow calculation method and device 
CN109344361B (en) *  20180827  20220520  南昌大学  Method for quickly forming Jacobian matrix in power system load flow calculation 
CN111478333B (en) *  20200414  20211130  广东电网有限责任公司广州供电局  Parallel static security analysis method for improving power distribution network recovery after disaster 
Citations (5)
Publication number  Priority date  Publication date  Assignee  Title 

CN101976835A (en) *  20101011  20110216  重庆大学  Parallel computation method for Newton power flow of largescale electric power system 
WO2012061674A3 (en) *  20101104  20130704  Siemens Corporation  Stochastic state estimation for smart grids 
CN103617150A (en) *  20131119  20140305  国家电网公司  GPU (graphic processing unit) based parallel power flow calculation system and method for largescale power system 
CN103793590A (en) *  20121101  20140514  同济大学  GPUbased computation method for quickly solving power flow in distribution networks 
CN104967121A (en) *  20150713  20151007  中国电力科学研究院  Largescale electric power system node load flow computing method 

2015
 20151120 CN CN201510809809.3A patent/CN105391057B/en not_active Expired  Fee Related
Patent Citations (5)
Publication number  Priority date  Publication date  Assignee  Title 

CN101976835A (en) *  20101011  20110216  重庆大学  Parallel computation method for Newton power flow of largescale electric power system 
WO2012061674A3 (en) *  20101104  20130704  Siemens Corporation  Stochastic state estimation for smart grids 
CN103793590A (en) *  20121101  20140514  同济大学  GPUbased computation method for quickly solving power flow in distribution networks 
CN103617150A (en) *  20131119  20140305  国家电网公司  GPU (graphic processing unit) based parallel power flow calculation system and method for largescale power system 
CN104967121A (en) *  20150713  20151007  中国电力科学研究院  Largescale electric power system node load flow computing method 
NonPatent Citations (3)
Title 

GPUbased power flow analysis with Chebyshev preconditioner andconjugate gradient method;Xue Li，Fangxing Li;《Electric Power Systems Research》;20140614;全文 * 
基于GPU的电力系统并行潮流计算的实现;夏俊峰 等;《电力系统保护与控制》;20100916;第38卷(第18期);全文 * 
基于道路树分层的大电网潮流并行算法及其GPU优化实现;陈德扬;《电力系统自动化》;20141125;第38卷(第22期);全文 * 
Also Published As
Publication number  Publication date 

CN105391057A (en)  20160309 
Similar Documents
Publication  Publication Date  Title 

CN105391057B (en)  A kind of GPU threaded design methods that direction of energy Jacobi battle array calculates  
CN106874113A (en)  A kind of many GPU heterogeneous schemas static security analysis computational methods of CPU+  
CN106157176B (en)  A kind of LU decomposition method for the direction of energy Jacobian matrix that GPU accelerates  
CN105576648A (en)  Static security analysis doublelayer parallel method based on GPUCUP heterogeneous computing platform  
CN102842917B (en)  A kind of general parallel net type photovoltaic generating system machineelectricity transient model  
CN104967121B (en)  A kind of tidal current computing method of largescale electrical power system node  
CN103455948B (en)  A kind of distribution system multidimensional multiresolution Modeling and the method for analysis  
CN103337861A (en)  Power distribution network reactive power optimization method based on gold chaotic ecological niche particle swarm algorithm  
CN106026107B (en)  A kind of QR decomposition method for the direction of energy Jacobian matrix that GPU accelerates  
CN106407158A (en)  GPU accelerated method for performing batch processing of isomorphic sparse matrixes multiplied by full vectors  
CN106058863A (en)  Random optimal trend calculation method based on random response surface method  
CN101976838B (en)  Newtonprocess power flow calculation method for study purpose  
CN105680473A (en)  Physical fusion modeling method for general electromechanical transient information of photovoltaic power generation system  
CN103279824A (en)  Modeling method for relay protection setting calculation system  
CN108258725B (en)  Doublyfed wind turbine dynamic equivalence method based on equivalent power angle coherence  
CN106354479B (en)  A kind of GPU acceleration QR decomposition method of a large amount of isomorphism sparse matrixes  
CN104734148A (en)  Threephrase powerdistributing network continuation power flow analysis of distributed power supply  
CN105896558B (en)  VSCbased UPFC electromechanical transient modular modeling method  
CN106294022B (en)  A kind of Jacobian matrix redundancy storage method for static security analysis  
CN107368454A (en)  A kind of GPU of the sparse lower trigonometric equation group of a large amount of isomorphisms pushes away method before accelerating  
CN106410811B (en)  Iteration small impedance branches endpoint changes the tidal current computing method of Jacobian matrix for the first time  
CN109698516A (en)  The maximum capacity computing system and method for renewable energy access power distribution network  
Mao et al.  Improved fast shortterm wind power prediction model based on superposition of predicted error  
Bai et al.  Heuristic optimization for wind energy integrated optimal power flow  
CN104182909A (en)  Multicore parallel successive approximation method of hydropower system optimal scheduling 
Legal Events
Date  Code  Title  Description 

C06  Publication  
PB01  Publication  
C10  Entry into substantive examination  
SE01  Entry into force of request for substantive examination  
GR01  Patent grant  
GR01  Patent grant  
CF01  Termination of patent right due to nonpayment of annual fee 
Granted publication date: 20171114 Termination date: 20181120 

CF01  Termination of patent right due to nonpayment of annual fee 