CN111711186B - Power system PQ decomposition state estimation method based on GPU parallel computing - Google Patents
Power system PQ decomposition state estimation method based on GPU parallel computing Download PDFInfo
- Publication number
- CN111711186B CN111711186B CN202010470768.0A CN202010470768A CN111711186B CN 111711186 B CN111711186 B CN 111711186B CN 202010470768 A CN202010470768 A CN 202010470768A CN 111711186 B CN111711186 B CN 111711186B
- Authority
- CN
- China
- Prior art keywords
- node
- branch
- active
- matrix
- est
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 28
- 239000011159 matrix material Substances 0.000 claims abstract description 161
- 238000005259 measurement Methods 0.000 claims abstract description 151
- 239000013598 vector Substances 0.000 claims abstract description 59
- 238000004364 calculation method Methods 0.000 claims description 36
- 238000002347 injection Methods 0.000 claims description 26
- 239000007924 injection Substances 0.000 claims description 26
- 230000008569 process Effects 0.000 claims description 25
- 230000008859 change Effects 0.000 claims description 23
- 239000000243 solution Substances 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 6
- 230000005611 electricity Effects 0.000 claims description 4
- 230000006835 compression Effects 0.000 claims description 3
- 238000007906 compression Methods 0.000 claims description 3
- DGAQECJNVWCQMB-PUAWFVPOSA-M Ilexoside XXIX Chemical compound C[C@@H]1CC[C@@]2(CC[C@@]3(C(=CC[C@H]4[C@]3(CC[C@@H]5[C@@]4(CC[C@@H](C5(C)C)OS(=O)(=O)[O-])C)C)[C@@H]2[C@]1(C)O)C)C(=O)O[C@H]6[C@@H]([C@H]([C@@H]([C@H](O6)CO)O)O)O.[Na+] DGAQECJNVWCQMB-PUAWFVPOSA-M 0.000 claims 1
- 229910052708 sodium Inorganic materials 0.000 claims 1
- 239000011734 sodium Substances 0.000 claims 1
- 238000003491 array Methods 0.000 description 4
- 238000013480 data collection Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005538 encapsulation Methods 0.000 description 2
- 238000007667 floating Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
Images
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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
- Y02E40/30—Reactive power compensation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Power Engineering (AREA)
- Supply And Distribution Of Alternating Current (AREA)
Abstract
The invention provides a power system PQ decomposition state estimation method based on GPU parallel computing, which comprises the following steps: collecting measurement parameters and network parameters of a power system; generating a node admittance matrix and a branch admittance list; initializing state variable vectors, calculating an active jacobian matrix, a reactive jacobian matrix, an active residual equation matrix and a reactive residual equation matrix, and respectively performing LLT decomposition on the active residual equation matrix and the reactive residual equation matrix; copying an active jacobian matrix, a reactive jacobian matrix, an active lower triangular matrix, a reactive lower triangular matrix and a branch admittance list into a main memory of the GPU; and (3) performing iterative computation in the GPU, returning to obtain a converged state estimation value, and returning to an error if iteration is not converged. The method provided by the invention can improve the speed of power system state estimation, ensure that the power dispatching system timely reflects the real running state of the power grid, and ensure the safe and efficient running of the power system.
Description
Technical Field
The invention relates to the technical field of power system state estimation, in particular to a power system PQ decomposition state estimation method based on GPU parallel computing.
Background
The power system state estimation function is to estimate the current running state of the power system according to various measurement information of the power system. The power system state estimation is the basis of each large-scale application in the power grid dispatching system, and is the basis of maintaining the stability of the power system and ensuring the safe and efficient operation of the power system, so that the calculation response speed of the power system state estimation determines whether the power grid dispatching system can timely reflect the real state of the power grid, and also influences the control effect of other power grid automatic control systems, and the stable operation of the whole power system is important.
At present, a power dispatching system often adopts a least square method and a PQ decomposition method to perform state estimation on the power system. In the calculation process of the two methods, a large number of mutually independent calculation processes and matrix operation exist, the capability of the power dispatching system for timely reflecting the real state of the power grid cannot be ensured, the parallelization processing is necessary in the calculation process of the state estimation for improving the real-time property of the power dispatching system for reflecting the real state of the power grid, the parallel calculation method for the state estimation of the power system is disclosed in the China patent No. CN1070696A in 2017 month 8, the two processes of state estimation parallel model establishment and parallel decomposition are integrated, the complexity of the state estimation by adopting an optimization algorithm is overcome, the parallelization processing is conveniently realized by combining the related parallel platform, and the calculation response speed of the state estimation can be improved.
Disclosure of Invention
In order to overcome the defect that the traditional power system state estimation method cannot guarantee that a power dispatching system timely reflects the real state of a power grid, the currently disclosed power system state estimation parallel calculation method is complex in related process and has the defect of non-uniform encapsulation integration improvement; the invention provides a power system PQ decomposition state estimation method based on GPU parallel computing, which improves the speed of power system state estimation, ensures that a power dispatching system timely reflects the real running state of a power grid, and ensures the safe and efficient running of the power system.
In order to achieve the technical effects, the technical scheme of the invention is as follows:
the invention provides a power system PQ decomposition state estimation method based on GPU parallel computing, which at least comprises the following steps:
s1, collecting measurement parameters of an electric power system and network parameters of the electric power system;
s2, generating a node admittance matrix and a branch admittance list according to network parameters of the power system;
s3, using node voltage amplitude and voltage phase angle of the power system as state variable vectors, initializing the node voltage amplitude to 1.0pu, and initializing the node voltage phase angle to 0rad;
s4, respectively calculating an active jacobian matrix H according to the measurement parameters of the power system and the network parameters of the power system P Reactive jacobian matrix H Q ；
S5, calculating an active residual equation matrix A P And reactive residual equation matrix A Q For the active residual equation matrix A P And reactive residual equation matrix A Q Respectively performing LLT decomposition to obtain an active lower triangular matrix A PL Reactive lower triangular matrix A QL ；
S6, an active jacobian matrix H P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL Copying the branch admittance list into a main memory of the GPU;
s7, based on GPU parallel calculation, carrying out active iteration in the GPU, judging whether an active iteration convergence condition is met, if yes, executing a step S8, otherwise, updating a node voltage phase angle, and executing the step S8;
s8, performing reactive iteration in the GPU, judging whether reactive iteration convergence conditions are met, and if yes, executing a step S9; otherwise, updating the node voltage amplitude, and executing step S9;
s9, judging whether the active iteration and the reactive iteration are converged or not, and if yes, outputting a state estimation result; otherwise, executing step S10;
s10, judging whether the iteration number kk is greater than the maximum iteration number k th If yes, the state estimation iteration is not converged; otherwise, the process returns to step S7.
The measurement parameters of the power system are collected through a data collection system, the network parameters of the power system are collected through an Energy Management System (EMS), and the node admittance matrix is generated according to the network parameters of the power system in step S2 according to a traditional calculation method; in addition, for the active residual equation matrix A P And reactive residual equation matrix A Q Respectively performing LLT decomposition is a relatively mature traditional technology; in the whole state estimation, if the current iteration times are compared with the set maximum iteration times when the active iteration convergence condition and the reactive iteration convergence condition are not met, if the current iteration times exceed the maximum iteration times, but the state estimation process still does not meet the active iteration convergence condition and the reactive iteration convergence condition, the state estimation process is not converged; when the active iteration convergence condition and the reactive iteration convergence condition are both satisfied, the output state estimation result is node voltage amplitude and voltage phase angle data of the power system, namely the current value of the state variable vector.
Preferably, the measurement parameters of the power system in step S1 include: node voltage U and node injection active power P inj_M Node injection reactive power Q inj_M Branch active power P ft_M 、P tf_M Branch reactive power Q ft_M 、Q tf_M A current magnitude; active value vector Z formed by measuring parameters P_M Vector Z of reactive power measurement values Q_M And a weight matrix R for measuring the power P And reactive power measurement weight matrix R Q The method comprises the steps of carrying out a first treatment on the surface of the The collected invalid or actually undelivered measurement parameters are set to 0, and the weight coefficient is also 0.
The network parameters of the power system include: node number nodeNum, branch number lnNum, branch head node number f, branch end node number t, branch resistance R, branch reactance X, branch nonstandard transformation ratio K, earth admittance B and node parallel conductance G sh Node parallel susceptance B sh Node parallel conductance G of non-parallel element node sh Node parallel susceptance B sh Is 0.
The node voltage U information in the acquired network parameters comprises the amplitude value and the voltage phase angle of the node voltage; in addition, in order to adapt to the subsequent parallel computation performed in the GPU, the CPU memory is considered to read the parameters in units of 128 floating point numbers or integers, so the parameter data collected in step S1 are all aligned to be processed: the length of the measurement vector is extended to an integer multiple of 128, and the node number nodeNum and the branch number lnNum are also extended to an integer multiple of 128.
The branch admittance list in step S2 sequentially comprises the following admittance storage structures from one side to the other side: branch admittance column: y is Y s [1]～Y s [lnNum]Branch f side-to-ground admittance column: y is Y f [1]～Y f [lnNum]Branch t side to ground admittance column: y is Y t [1]～Y t [lnNum]；
Any branch k admittance Y in branch admittance column s (k) The solution formula of (2) is:
wherein ,Ys (k) Representing the branch admittance, R (k) representing the branch resistance, X (k) representing the branch reactance; k (K) represents the non-standard transformation ratio of the branch; g s (k) Conductance for the branch; b s (k) Is thatThe branch susceptance; j is the imaginary part mark in complex calculation;
the kth branch f side-to-ground admittance Y in the branch f side-to-ground admittance column f (k) The solution formula of (2) is:
wherein B (k) represents the capacitance to ground of the line, g f (k) Conducting ground for the branch f side; b f (k) A susceptance to ground for the branch f side;
the kth branch t side-to-ground admittance Y in the branch t side-to-ground admittance column t (k) The solution formula of (2) is:
wherein ,gt (k) Conducting electricity to the side t of the branch; b t (k) And susceptance is applied to the t side of the branch.
In order to improve the speed of acquiring parameter data by the GPU, prevent access conflict among threads of the GPU, ensure the speed of parallel computation of the GPU, store the branch admittance list in sequence according to the admittance storage structure from one side to the other side, and the admittance parameter of the branch which does not exist actually is 0.
Preferably, the length of the state variable vector in step S3 is: nodeNum x 2.
Preferably, the active jacobian matrix H described in step S4 P For a sparse matrix of nodenum+lnnum row by 2 nodeNum column, an active jacobian matrix H P The element calculation formula of (a) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pinj (i) Injecting active power for point i; b (B) ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, θ, of node i i A voltage phase angle theta of the node i j The voltage phase angle of the node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) Is the voltage phase angle of the last node of branch k. f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j.
Reactive jacobian matrix H Q A sparse matrix of nodeNum x 2+lnnum x 2 rows and nodeNum columns, and a reactive jacobian matrix H Q The element calculation formula of (a) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qinj (i) Injecting reactive power for node i, B ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, U, of node i i Voltage amplitude, U, for node i j The voltage amplitude at node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) Is the voltage phase angle of the last node of branch k. f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j.
Preferably, the active residual equation matrix A described in step S5 P The calculation formula of (2) is as follows:
wherein ,
as an active jacobian matrix H P Transpose of->
For measuring weight matrix R of active quantity P An inverse matrix of (a); reactive residual equation matrix A Q The calculation formula of (2) is as follows:
wherein ,
is a reactive jacobian matrix H Q Transpose of->
For measuring weight matrix R of reactive power Q Is a matrix of inverse of (a).
Here, due to the active weight matrix R P And reactive power measurement weight matrix R Q Are all diagonal arrays, and when the measurement parameters of the power system are collected, the collected invalid or actually undelivered measurement parameters are set to 0, so that the active measurement weight matrix R P And reactive power measurement weight matrix R Q In which a part of diagonal elements are set to 0, and there is a case of no inverse matrix, for the active weight matrix R P And reactive power measurement weight matrix R Q The non-zero elements in (1) are inverted to perform equivalent inversion.
Preferably, the active jacobian matrix H described in step S6 P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL When copying to the main memory of the GPU, adopting a column priority sparse matrix compression storage format; copying the branch admittance list into the main memory of the GPU, and adopting a column-priority storage format as an active jacobian matrix H P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL The number of lines of any one parameter column in the branch admittance list is unequalAnd when the number of the GPU threads is integer multiple, inserting all 0 parameter empty lines to complete, thereby ensuring the combined access of the triggering GPU to the parameters copied to the main memory.
Preferably, the active iteration process described in step S7 is:
s701, calculating an active power measurement estimated value;
firstly, calculating the branch active power measurement estimated value, wherein the calculation formula of the branch active power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pft_est (k) Measuring an estimated value for active power of the branch from the f side to the t side; p (P) tf_est (k) Measuring an estimated value for active power of the branch from the t side to the f side; u (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j.
S702, traversing all nodes, and setting an active injection measurement estimated value variable P for any point i inj_est (i) And initialized to 0.
S703, traversing all branches, and for any branch k, superposing the active power measurement estimated values on two sides of the branch to the injection active power measurement estimated value of the measurement node. For example, for the kth branch
P inj_est (f(k))＝P inj_est (f(k))+P ft_est (k)
P inj_est (t(k))＝P inj_est (t(k))+P tf_est (k)
wherein ,Pft_est (f) Injecting a power measurement estimated value, P, into the node of the branch f-side node tf_est And (t) injecting a useful measuring estimated value for the node of the node at the t side of the branch.
S704, traversing all nodes, superposing the node parallel admittance active power measurement estimation value to the node injection active power measurement estimation value, and for any node i:
wherein ,Gsh (i) Parallel conductance for node i
S705. P inj_est 、P ft_est 、P tf_est Sequentially arranged to form a vector z of the estimated value of the active power measurement P_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z P_est And active measured value vector Z P_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s706, according to an active residual equation matrix A P Error vector z of active measurement estimation P_err The change amount delta theta of the state variable theta is obtained, and the change amount delta theta of the state variable theta satisfies the following conditions:
preferably, the reactive iteration process described in step S7 is:
s711, calculating a reactive power measurement estimated value;
firstly, calculating a branch reactive power measurement estimated value, wherein the calculation formula of the branch reactive power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est And (t) injecting reactive power measurement estimated values into the node of the node at the t side of the branch. U (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j.
S712, traversing all nodes, and setting an active injection measurement estimated value variable Q for any point i inj_est (i) And initialized to 0.
S713, traversing all branches, and for any branch k, superposing reactive power measurement estimated values at two sides of the branch to the injection reactive power measurement estimated value of the measurement node. For example, for the kth branch
Q inj_est (f(k))＝Qinj_est(f(k))+Q ft_est (k)
Q inj_est (t(k))＝Qinj_est(t(k))+Q tf_est (k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est And (t) injecting reactive power measurement estimated values into the node of the node at the t side of the branch.
S714, traversing all nodes, superposing the node parallel admittance reactive power measurement estimation value to the node injection reactive power measurement estimation value, and for any node i:
wherein ,Bsh (i) And susceptance is connected in parallel for the node i.
S715. Will Q inj_est 、Q ft_est 、Q tf_est U is arranged in turnThe columns form a reactive power measurement estimate vector z Q_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z Q_est And active measured value vector z Q_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s716, according to an active residual equation matrix A Q Error vector z of active measurement estimation Q_err The change quantity delta U of the state variable U is obtained, and the change quantity delta U of the state variable U meets the following conditions:
in the process of active iteration or reactive iteration, as the parallel calculation is realized by each thread of the GPU in the form of single-instruction multi-data stream, the number of branches connected by each node is different, and the traditional method for calculating the active or reactive power measurement estimated value according to the nodes can cause waiting among the threads, so that the parallelism is seriously influenced, the active and reactive power measurement estimated values of the branches and the parallel admittance active and reactive power measurement estimated values of all the nodes are calculated at first, the power of the branches is summed to obtain the node injection power, and the calculated parallelism is improved; in addition, in the process of solving the change amount Δu of the state variable U or the change amount Δθ of the state variable θ, the state estimation is based on the PQ decomposition method, and the active residual equation matrix a is then P Reactive residual equation matrix A Q All are constant matrixes, and do not change in the iterative process, so that the lower triangular arrays APL and AQL respectively obtained by LLT decomposition are pushed forward for generation, and delta theta and delta U are further obtained, so that the solving process is simplified.
Preferably, the active iteration convergence condition is:
max(Δθ)<θ th
wherein Δθ represents the amount of change in the state variable θ; θ th Representing an active iteration convergence thresholdThe method comprises the steps of carrying out a first treatment on the surface of the When the active iteration convergence condition is not satisfied, updating a node voltage phase angle, wherein an updating formula of the node voltage phase angle is as follows:
θ (kk+1) ＝θ (kk) +Δθ
wherein ,θ(kk) Representing a node voltage phase angle state variable at a kth iteration; θ (kk+1) Representing the node voltage phase angle state variable at the kk+1th iteration; the reactive iteration convergence condition is:
max(ΔU)<U th
wherein DeltaU represents the variation of the voltage amplitude U of the state variable node; u (U) th Representing a reactive iteration convergence threshold; when the reactive iteration convergence condition is not satisfied, updating the node voltage amplitude, wherein the updating formula of the node voltage amplitude is as follows:
U( kk+1 )＝U( kk )+ΔU
wherein ,U(kk) Representing a node voltage magnitude state variable at a kth iteration; u (U) (kk+1) Representing the node voltage magnitude state variable at the kk+1th iteration.
Compared with the prior art, the technical scheme of the invention has the beneficial effects that:
the invention provides a power system PQ decomposition state estimation method based on GPU parallel computing, which is based on state estimation of a traditional power system PQ decomposition method, uniformly copies and packages parameters related in the state estimation process into a main memory of a GPU, triggers the GPU to merge and access the parameters copied into the main memory, and then carries out parallel computing based on the GPU.
Drawings
FIG. 1 is a flow chart of a method for estimating the PQ decomposition state of a power system based on GPU parallel computing according to the present invention;
fig. 2 shows a branch admittance list encapsulation structure diagram proposed in an embodiment of the present invention.
Detailed Description
The drawings are for illustrative purposes only and are not to be construed as limiting the present patent;
for better illustration of the present embodiment, some parts of the drawings may be omitted, enlarged or reduced, and do not represent actual dimensions;
it will be appreciated by those skilled in the art that some well known descriptions in the figures may be omitted.
The technical scheme of the invention is further described below with reference to the accompanying drawings and examples.
Example 1
Fig. 1 is a flow chart of a power system PQ decomposition state estimation method based on GPU parallel computing according to the present invention, see fig. 1, wherein H in the illustration P Representing an active jacobian matrix; h Q Representing a reactive jacobian matrix; a is that P Representing an active residual equation matrix; a is that Q Representing a reactive residual equation matrix; a is that PL Representing an active lower triangular matrix; a is that QL Representing a reactive lower triangular matrix, the method comprising:
s1, collecting measurement parameters of an electric power system and network parameters of the electric power system; in specific implementation, measurement parameters of the power system are collected through a data collection system, and network parameters of the power system are collected through an Energy Management System (EMS);
the measurement parameters of the power system include: node voltage U and node injection active power P inj_M Node injection reactive power Q inj_M Active power P on both sides of branch ft_M 、P tf_M Reactive power Q at two sides of branch ft_M 、Q tf_M The current amplitude value, the collected invalid or actually undelivered measurement parameter is set to 0, and the weight coefficient is also 0; active value vector Z formed by measuring parameters P_M Vector Z of reactive power measurement values Q_M And a weight matrix R for measuring the power P And reactive power measurement weight matrix R Q The method comprises the steps of carrying out a first treatment on the surface of the The node voltage U information in the acquired network parameters comprises the amplitude value and the voltage phase angle of the node voltage;
network parameter packet of the power systemThe method comprises the following steps: node number nodeNum, branch number lnNum, branch head node number f, branch end node number t, branch resistance R, branch reactance X, branch nonstandard transformation ratio K, earth admittance B and node parallel conductance G sh Node parallel susceptance B sh Node parallel conductance G of non-parallel element node sh Node parallel susceptance B sh Is 0. In addition, in order to adapt to the parallel computation performed in the GPU in the subsequent implementation, the CPU memory is considered to read the parameters in units of 128 floating point numbers or integers, so the parameter data collected in step S1 are all aligned to be processed: the length of the measurement vector is extended to an integer multiple of 128, and the node number nodeNum and the branch number lnNum are also extended to an integer multiple of 128.
S2, generating a node admittance matrix and a branch admittance list according to network parameters of the power system; in this embodiment, the node admittance matrix is generated according to a conventional calculation method;
s3, using node voltage amplitude and voltage phase angle of the power system as state variable vectors, initializing the node voltage amplitude to 1.0pu, and initializing the node voltage phase angle to 0rad;
s4, respectively calculating an active jacobian matrix H according to the measurement parameters of the power system and the network parameters of the power system P Reactive jacobian matrix H Q ；
S5, calculating an active residual equation matrix A P And reactive residual equation matrix A Q For the active residual equation matrix A P And reactive residual equation matrix A Q Respectively performing LLT decomposition to obtain an active lower triangular matrix A PL Reactive lower triangular matrix A QL ；
S6, an active jacobian matrix H P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL Copying the branch admittance list into a main memory of the GPU;
s7, based on GPU parallel calculation, carrying out active iteration in the GPU, judging whether an active iteration convergence condition is met, if yes, executing a step S8, otherwise, updating a node voltage phase angle, and executing the step S8;
s8, performing reactive iteration in the GPU, judging whether reactive iteration convergence conditions are met, and if yes, executing a step S9; otherwise, updating the node voltage amplitude, and executing step S9;
s9, judging whether the active iteration and the reactive iteration are converged or not, and if yes, outputting a state estimation result; otherwise, executing step S10;
s10, judging whether the iteration number kk is greater than the maximum iteration number k th If yes, the state estimation iteration is not converged; otherwise, the process returns to step S7.
In the whole state estimation, if the current iteration times are compared with the set maximum iteration times when the active iteration convergence condition and the reactive iteration convergence condition are not met, if the current iteration times exceed the maximum iteration times, but the state estimation process still does not meet the active iteration convergence condition and the reactive iteration convergence condition, the state estimation process is not converged; when the active iteration convergence condition and the reactive iteration convergence condition are both satisfied, the output state estimation result, namely the current value of the state variable vector, is node voltage amplitude and voltage phase angle data of the power system.
In this embodiment, as shown in fig. 2, in order to increase the speed of the GPU obtaining parameter data, access conflict between threads of the GPU is prevented, the speed of GPU parallel computing is ensured, and the admittance storage structure of the branch admittance list from one side to the other side is as follows: branch admittance column: y is Y s [1]～Y s [lnNum]Branch f side-to-ground admittance column: y is Y f [1]～Y f [lnNum]Branch t side to ground admittance column: y is Y t [1]～Y t [lnNum]The method comprises the steps of carrying out a first treatment on the surface of the The admittance parameter of the branch that is not actually present is 0.
Any branch k admittance Y in branch admittance column s (k) The solution formula of (2) is:
wherein ,Ys (k) Representing the branch admittance, R (k) represents the branchA path resistance, X (k) representing the path reactance; k (K) represents the non-standard transformation ratio of the branch; g s (k) Conductance for the branch; b s (k) Susceptance for the branch; j is the imaginary part mark in complex calculation;
the kth branch f side-to-ground admittance Y in the branch f side-to-ground admittance column f (k) The solution formula of (2) is:
wherein B (k) represents the capacitance to ground of the line, g f (k) Conducting ground for the branch f side; b f (k) A susceptance to ground for the branch f side;
the kth branch t side-to-ground admittance Y in the branch t side-to-ground admittance column t (k) The solution formula of (2) is:
wherein ,gt (k) Conducting electricity to the side t of the branch; b t (k) And susceptance is applied to the t side of the branch.
The length of the state variable vector in step S3 is as follows: nodeNum x 2.
In this embodiment, the active jacobian matrix H described in step S4 P As a sparse matrix of (nodenum+lnnum) row nodeNum columns, an active jacobian matrix H P The element calculation formula of (a) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pinj (i) Injecting active power for point i; b (B) ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, θ, of node i i A voltage phase angle theta of the node i j The voltage phase angle of the node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) Is the voltage phase angle of the last node of branch k. f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j.
Reactive jacobian matrix H Q Is a sparse matrix of (nodeNum x 2+lnnum x 2) row nodeNum columns, reactive jacobian matrix H Q The element calculation formula of (a) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qinj (i) Injecting reactive power for node i, B ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, U, of node i i Voltage amplitude, U, for node i j The voltage amplitude at node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) Is the voltage phase angle of the last node of branch k. f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j.
In this embodiment, the active residual equation matrix A described in step S5 P The calculation formula of (2) is as follows:
wherein ,
as an active jacobian matrix H P Transpose of->
For measuring weight matrix R of active quantity P An inverse matrix of (a); reactive residual equation matrix A Q The calculation formula of (2) is as follows:
wherein ,
is a reactive jacobian matrix H Q Transpose of->
Weight matrix R for reactive power measurement Q Is a matrix of inverse of (a). In the concrete implementation, the weight matrix R is measured due to the active quantity P And reactive power measurement weight matrix R Q Are all diagonal arrays, and when the measurement parameters of the power system are collected, the collected invalid or actually undelivered measurement parameters are set to 0, so that the active measurement weight matrix R P And reactive power measurement weight matrix R Q In which a part of diagonal elements are set to 0, and there is a case of no inverse matrix, for the active weight matrix R P And reactive power measurement weight matrix R Q The non-zero elements in (1) are inverted to perform equivalent inversion.
In this embodiment, an active jacobian matrix H as described in step S6 P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL When copying to the main memory of the GPU, adopting a column priority sparse matrix compression storage format; copying the branch admittance list into the main memory of the GPU, and adopting a column-priority storage format as an active jacobian matrix H P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL And when the number of lines of any parameter column in the branch admittance list is not equal to the integer multiple of the number of GPU threads, inserting all 0 parameter blank lines to complete, thereby ensuring the combined access of triggering the GPU to the parameters copied to the main memory.
In this embodiment, the active iteration process described in step S7 is:
s701, calculating an active power measurement estimated value;
firstly, calculating the branch active power measurement estimated value, wherein the calculation formula of the branch k active power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pft_est (k) Measuring an estimated value for active power of the branch from the f side to the t side; p (P) tf_est (k) Measuring an estimated value for active power of the branch from the t side to the f side; u (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j.
S702, traversing all nodes, and setting an active injection measurement estimated value variable P for any point i inj_est (i) And initialized to 0.
S703, traversing all branches, and for any branch k, superposing the active power measurement estimated values on two sides of the branch to the injection active power measurement estimated value of the measurement node. For example, for the kth branch
P inj_est (f(k))＝P inj_est (f(k))+P ft_est (k)
P inj_est (t(k))＝P inj_est (t(k))+P tf_est (k)
wherein ,Pft_est (f) Injecting a power measurement estimated value, P, into the node of the branch f-side node tf_est And (t) injecting a useful measuring estimated value for the node of the node at the t side of the branch.
S704, traversing all nodes, superposing the node parallel admittance active power measurement estimation value to the node injection active power measurement estimation value, and for any node i:
wherein ,Gsh (i) Parallel conductance for node i
S705. P inj_est 、P ft_est 、P tf_est Sequentially arranged to form a vector z of the estimated value of the active power measurement P_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z P_est And active measured value vector Z P_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s706, according to an active residual equation matrix A P Error vector z of active measurement estimation P_err The change amount delta theta of the state variable theta is obtained, and the change amount delta theta of the state variable theta satisfies the following conditions:
will P i_est 、P ij_est 、P ji_est Sequentially arranged to form a vector z of the estimated value of the active power measurement P_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z P_est And active measured value vector Z P_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s706, according to an active residual equation matrix A P Error vector z of active measurement estimation P_err The change amount delta theta of the state variable theta is obtained, and the change amount delta theta of the state variable theta satisfies the following conditions:
preferably, the reactive iteration process described in step S7 is:
s711, calculating a reactive power measurement estimated value;
firstly, calculating a branch reactive power measurement estimated value, wherein the calculation formula of the branch reactive power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est And (t) injecting reactive power measurement estimated values into the node of the node at the t side of the branch. U (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j.
S712, traversing all nodes, and setting an active injection measurement estimated value variable Q for any point i inj_est (i) And initialized to 0.
S713, traversing all branches, and for any branch k, superposing reactive power measurement estimated values at two sides of the branch to the injection reactive power measurement estimated value of the measurement node. For example, for the kth branch
Q inj_est (f(k))＝Q inj_est (f(k))+Q ft_est (k)
Q inj_est (t(k))＝Q inj_est (t(k))+Q tf_est (k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est And (t) injecting reactive power measurement estimated values into the node of the node at the t side of the branch.
S714, traversing all nodes, superposing the node parallel admittance reactive power measurement estimation value to the node injection reactive power measurement estimation value, and for any node i:
wherein ,Bsh (i) And susceptance is connected in parallel for the node i.
S715. Will Q i_est 、Q ij_est 、Q ji_est Sequentially arranged to form a reactive power measurement estimated value vector z Q_est The method comprises the steps of carrying out a first treatment on the surface of the Using the reactive power measurement estimate vector z Q_est And reactive power measurement vector z Q_M Calculating a reactive power measurement estimation error vector:
z Q_err ＝z Q_M -z Q_est
wherein ,zQ_err Representing a reactive power measurement estimation error vector;
s713, according to reactive residual equation matrix A Q Reactive power measurement estimation error vector z Q_err The change quantity delta U of the state variable U is obtained, and the change quantity delta U of the state variable U meets the following conditions:
in the process of active iteration or reactive iteration, when calculating the active or reactive power measurement estimated value, as each thread of the GPU realizes parallel calculation in the form of single-instruction multi-data stream, the number of branches connected by each node is different, and the traditional method for calculating the active or reactive power measurement estimated value by node can cause waiting among each thread and seriously affect the parallelism, firstly calculating the active and reactive power measurement estimated values of the branches and the parallel admittance active and reactive power measurement estimated values of all nodes, and solving the branch powerAnd obtaining node injection power, so as to improve the parallelism of calculation; in addition, in the process of solving the change amount Δu of the state variable U or the change amount Δθ of the state variable θ, the state estimation is based on the PQ decomposition method, and the active residual equation matrix a is then P Reactive residual equation matrix A Q All are constant matrixes, and do not change in the iterative process, so that the lower triangular arrays APL and AQL respectively obtained by LLT decomposition are pushed forward for generation, and delta theta and delta U are further obtained, so that the solving process is simplified.
The active iteration convergence condition is:
max(Δθ)<θ th
wherein Δθ represents the amount of change in the state variable θ; θ th Representing an active iteration convergence threshold; when the active iteration convergence condition is not satisfied, updating a node voltage phase angle, wherein an updating formula of the node voltage phase angle is as follows:
θ (kk+1) ＝θ (kk) +Δθ
wherein ,θ(kk) Representing a node voltage phase angle state variable at a kth iteration; θ (kk+1) Representing the node voltage phase angle state variable at the kk+1th iteration; the reactive iteration convergence condition is:
max(ΔU)<U th
wherein DeltaU represents the variation of the voltage amplitude U of the state variable node; u (U) th Representing a reactive iteration convergence threshold; when the reactive iteration convergence condition is not satisfied, updating the node voltage amplitude, wherein the updating formula of the node voltage amplitude is as follows:
U( kk+1 )＝U( kk )+ΔU
wherein ,U(kk) Representing a node voltage magnitude state variable at a kth iteration; u (U) (kk+1) Representing the node voltage magnitude state variable at the kk+1th iteration.
The relationships depicted in the drawings are for illustrative purposes only and are not to be construed as limiting the present patent;
it is to be understood that the above examples of the present invention are provided by way of illustration only and are not intended to limit the scope of the invention. Other variations or modifications of the above teachings will be apparent to those of ordinary skill in the art. It is not necessary here nor is it exhaustive of all embodiments. Any modification, equivalent replacement, improvement, etc. which come within the spirit and principles of the invention are desired to be protected by the following claims.
Claims (7)
1. The utility model provides a power system PQ decomposition state estimation method based on GPU parallel computing, which is characterized by at least comprising the following steps:
s1, collecting measurement parameters of an electric power system and network parameters of the electric power system;
the measurement parameters of the power system in step S1 include: node voltage U and node injection active power P inj_M Node injection reactive power Q inj_M Active power P on both sides of branch ft_M 、P tf_M Reactive power Q at two sides of branch ft_M 、Q tf_M A current magnitude; active value vector Z formed by measuring parameters P_M Vector Z of reactive power measurement values Q_M And a weight matrix R for measuring the power P And reactive power measurement weight matrix R Q The method comprises the steps of carrying out a first treatment on the surface of the The network parameters of the power system include: node number nodeNum, branch number lnNum, branch head node number f, branch end node number t, branch resistance R, branch reactance X, branch nonstandard transformation ratio K, earth admittance B and node parallel conductance G sh Node parallel susceptance B sh ；
S2, generating a node admittance matrix and a branch admittance list according to network parameters of the power system;
the branch admittance list in step S2 sequentially comprises the following admittance storage structures from one side to the other side: branch admittance column: y is Y s [1]～Y s [lnNum]Branch f side-to-ground admittance column: y is Y f [1]～Y f [lnNum]Branch t side to ground admittance column: y is Y t [1]～Y t [lnNum]；
Any branch k admittance Y in branch admittance column s (k) The solution formula of (2) is:
wherein ,Ys (k) Representing the branch admittance, R (k) representing the branch resistance, X (k) representing the branch reactance; k (K) represents the non-standard transformation ratio of the branch; g s (k) Conductance for the branch; b s (k) Susceptance for the branch; j is the imaginary part mark in complex calculation;
the kth branch f side-to-ground admittance Y in the branch f side-to-ground admittance column f (k) The solution formula of (2) is:
wherein B (k) represents the capacitance of the branch to ground, g f (k) Conducting ground for the branch f side; b f (k) A susceptance to ground for the branch f side;
the kth branch t side-to-ground admittance Y in the branch t side-to-ground admittance column t (k) The solution formula of (2) is:
wherein ,gt (k) Conducting electricity to the side t of the branch; b t (k) A sodium is oppositely arranged for the t side of the branch;
s3, using node voltage amplitude and voltage phase angle of the power system as state variable vectors, initializing the node voltage amplitude to 1.0pu, and initializing the node voltage phase angle to 0.0rad;
s4, respectively calculating an active jacobian matrix H according to the measurement parameters of the power system and the network parameters of the power system P Reactive jacobian matrix H Q ；
An active jacobian matrix H as described in step S4 P For a sparse matrix of nodenum+lnnum row by 2 nodeNum column, an active jacobian matrix H P The element calculation formula of (a) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pinj (i) Injecting active power for point i; b (B) ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, θ, of node i i A voltage phase angle theta of the node i j The voltage phase angle of the node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) The voltage phase angle of the k end node of the branch is; f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j;
reactive jacobian matrix H Q A sparse matrix of nodeNum, 3+lnnum, 2 rows and nodeNum columns, and a reactive jacobian matrix H Q In (a) elementsThe element calculation formula is:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qinj (i) Injecting reactive power for node i, B ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, U, of node i i Voltage amplitude, U, for node i j The voltage amplitude at node j; u (U) f(k) The voltage amplitude of the starting node of the branch k; u (U) t(k) The voltage amplitude of the k-terminal node of the branch is; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase of the k-terminal node of the branchA corner; f (k) is the number of the starting node of the branch k, t (k) is the number of the last node of the branch k, and θ ft Representing the voltage phase angle difference between node i and node j;
s5, calculating an active residual equation matrix A P And reactive residual equation matrix A Q For the active residual equation matrix A P And reactive residual equation matrix A Q Respectively performing LLT decomposition to obtain an active lower triangular matrix A PL Reactive lower triangular matrix A QL ；
An active residual equation matrix A described in the step S5 P The calculation formula of (2) is as follows:
wherein ,
as an active jacobian matrix H P Transpose of->
For measuring weight matrix R of active quantity P An inverse matrix of (a);
reactive residual equation matrix A Q The calculation formula of (2) is as follows:
wherein ,
is a reactive jacobian matrix H Q Transpose of->
For reactive power measurement weight matrix R Q An inverse matrix of (a);
s6, an active jacobian matrix H P Reactive jacobian matrix H Q Active powerLower triangular matrix A PL Triangle matrix A under reactive power QL Copying the branch admittance list into a main memory of the GPU;
s7, based on GPU parallel calculation, carrying out active iteration in the GPU, judging whether an active iteration convergence condition is met, if yes, executing a step S8, otherwise, updating a node voltage phase angle, and executing the step S8;
s8, performing reactive iteration in the GPU, judging whether reactive iteration convergence conditions are met, and if yes, executing a step S9; otherwise, updating the node voltage amplitude, and executing step S9;
s9, judging whether the active iteration and the reactive iteration are converged or not, and if yes, outputting a state estimation result; otherwise, executing step S10;
s10, judging whether the iteration number kk is greater than the maximum iteration number k th If yes, the state estimation iteration is not converged; otherwise, the process returns to step S7.
2. The GPU-parallel computing-based power system PQ decomposition state estimation method according to claim 1, wherein the collected invalid or actually untransmitted measurement parameters are set to 0, and the weight coefficient is also 0; node parallel conductance G of non-parallel element node sh Node parallel susceptance B sh Is 0.
3. The GPU-parallel computing-based power system PQ decomposition state estimation method according to claim 1, wherein the length of the state variable vector in step S3 is: nodeNum x 2.
4. The method for estimating a state of power system PQ decomposition based on GPU parallel computing according to claim 1, wherein the active jacobian matrix H of step S6 P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL When copying to the main memory of the GPU, adopting a column priority sparse matrix compression storage format; copying the branch admittance list into the main memory of the GPU, and adopting a column-priority storage format as an active jacobian matrix H P Reactive jacobian matrix H Q Triangle matrix A under active power PL Triangle matrix A under reactive power QL And inserting all 0 parameter blank columns to complete when the number of lines of any parameter column in the branch admittance list is not equal to the integer multiple of the number of GPU threads.
5. The GPU-parallel computing-based power system PQ decomposition state estimation method according to claim 1, wherein the active iteration process of step S7 is:
s701, calculating an active power measurement estimated value;
firstly, calculating the branch active power measurement estimated value, wherein the calculation formula of the branch k active power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Pft_est (k) Measuring an estimated value for active power of the branch from the f side to the t side; p (P) tf_est (k) Measuring an estimated value for active power of the branch from the t side to the f side; u (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j;
s702, traversing all nodes, and setting an active injection measurement estimated value variable P for any point i inj_est (i) And initialized to 0;
s703, traversing all branches, and for any branch k, superposing active power measurement estimated values on two sides of the branch to the injection active power measurement estimated value of the measurement node; for the kth branch
P inj_est (f(k))＝P inj_est (f(k))+P ft_est (k)
P inj_est (t(k))＝P inj_est (t(k))+P tf_est (k)
wherein ,Pft_est (f) Injecting a power measurement estimated value, P, into the node of the branch f-side node tf_est (t) injecting a useful-load measurement estimated value for a node of the node at the t side of the branch;
s704, traversing all nodes, superposing the node parallel admittance active power measurement estimation value to the node injection active power measurement estimation value, and for any node i:
wherein ,Gsh (i) Parallel conductance for node i
S705. P inj_est 、P ft_est 、P tf_est Sequentially arranged to form a vector z of the estimated value of the active power measurement P_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z P_est And active measured value vector Z P_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s706, according to an active residual equation matrix A P Error vector z of active measurement estimation P_err The change amount delta theta of the state variable theta is obtained, and the change amount delta theta of the state variable theta satisfies the following conditions:
6. the GPU-parallel computing-based power system PQ decomposition state estimation method according to claim 5, wherein the reactive iteration process of step S7 is:
s711, calculating a reactive power measurement estimated value;
firstly, calculating a branch reactive power measurement estimated value, wherein the calculation formula of the branch reactive power measurement estimated value is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est (t) injecting a reactive power measurement estimated value for a node of the branch t-side node; u (U) f(k) Representing the node voltage magnitude on the branch f side; u (U) t(k) Representing the node voltage amplitude on the t side of the branch; θ f(k) The voltage phase angle of the starting node of branch k, theta t(k) For the voltage phase angle, θ, of the branch k-terminal node ft Representing the voltage phase angle difference between node i and node j;
s712, traversing all nodes, and setting an active injection measurement estimated value variable Q for any point i inj_est (i) And initialized to 0;
s713, traversing all branches, and for any branch k, superposing reactive power measurement estimated values at two sides of the branch to the injection reactive power measurement estimated value of the measurement node; for the kth branch
Q inj_est (f(k))＝Q inj_est (f(k))+Q ft_est (k)
Q inj_est (t(k))＝Q inj_est (t(k))+Q tf_est (k)
wherein ,Qft_est (f) Injecting reactive power measurement estimated value, Q for the node of the branch f-side node tf_est (t) injecting a reactive power measurement estimated value for a node of the branch t-side node;
s714, traversing all nodes, superposing the node parallel admittance reactive power measurement estimation value to the node injection reactive power measurement estimation value, and for any node i:
wherein ,Bsh (i) A susceptance is connected in parallel for the node i;
s715. Will Q inj_est 、Q ft_est 、Q tf_est U are sequentially arranged to form a reactive power measurement estimated value vector z Q_est The method comprises the steps of carrying out a first treatment on the surface of the Using active measurement estimate vector z Q_est And active measured value vector z Q_M Calculating an active measurement estimation error vector:
z P_err ＝z P_M -z P_est
wherein ,zP_err Representing a useful-quantity-measurement estimation error vector;
s716, according to an active residual equation matrix A Q Error vector z of active measurement estimation Q_err The change quantity delta U of the state variable U is obtained, and the change quantity delta U of the state variable U meets the following conditions:
7. the GPU-parallel computing-based power system PQ decomposition state estimation method of claim 6, wherein the active iteration convergence condition is:
max(Δθ)<θ th
wherein Δθ represents the amount of change in the state variable θ; θ th Representing an active iteration convergence threshold; when the active iteration convergence condition is not satisfied, updating a node voltage phase angle, wherein an updating formula of the node voltage phase angle is as follows:
θ (kk+1) ＝θ (kk) +Δθ
wherein ,θ(kk) Representing a node voltage phase angle state variable at a kth iteration; θ (kk+1) Representing node electricity at kk+1th iterationA phase angle state variable; the reactive iteration convergence condition is:
max(ΔU)<U th
wherein DeltaU represents the variation of the voltage amplitude U of the state variable node; u (U) th Representing a reactive iteration convergence threshold; when the reactive iteration convergence condition is not satisfied, updating the node voltage amplitude, wherein the updating formula of the node voltage amplitude is as follows:
U( kk+1 )＝U( kk )+ΔU
wherein ,U(kk) Representing a node voltage magnitude state variable at a kth iteration; u (U) (kk+1) Representing the node voltage magnitude state variable at the kk+1th iteration.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010470768.0A CN111711186B (en) | 2020-05-28 | 2020-05-28 | Power system PQ decomposition state estimation method based on GPU parallel computing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010470768.0A CN111711186B (en) | 2020-05-28 | 2020-05-28 | Power system PQ decomposition state estimation method based on GPU parallel computing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111711186A CN111711186A (en) | 2020-09-25 |
| CN111711186B true CN111711186B (en) | 2023-05-02 |
Family
ID=72538687
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010470768.0A Active CN111711186B (en) | 2020-05-28 | 2020-05-28 | Power system PQ decomposition state estimation method based on GPU parallel computing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111711186B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112749369A (en) * | 2021-01-19 | 2021-05-04 | 东方电子股份有限公司 | Power system state estimation method based on Givens orthogonal transformation |
| CN113315118B (en) * | 2021-04-26 | 2022-12-30 | 中国南方电网有限责任公司 | Power system state estimation method based on parallel computing and particle swarm optimization |
| CN114372235A (en) * | 2021-12-30 | 2022-04-19 | 东方电子股份有限公司 | Power system state estimation method for preventing singular transformation matrix |
| CN117060373B (en) * | 2023-06-27 | 2024-03-26 | 国网信息通信产业集团有限公司 | Active power distribution network state estimation method and device based on measurement and alignment |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2000059994A (en) * | 1998-08-05 | 2000-02-25 | Mitsubishi Electric Corp | Optimum power flow calculation device |
| CN103902814A (en) * | 2014-03-10 | 2014-07-02 | 中国南方电网有限责任公司 | Electric power system operation state detecting method based on dynamic partitioning |
| CN107016489A (en) * | 2017-03-09 | 2017-08-04 | 中国电力科学研究院 | A kind of electric power system robust state estimation method and device |
| CN107196306A (en) * | 2017-07-10 | 2017-09-22 | 大连海事大学 | Algorithm quicksort tidal current computing method based on Matlab sparse matrixes |
-
2020
- 2020-05-28 CN CN202010470768.0A patent/CN111711186B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2000059994A (en) * | 1998-08-05 | 2000-02-25 | Mitsubishi Electric Corp | Optimum power flow calculation device |
| CN103902814A (en) * | 2014-03-10 | 2014-07-02 | 中国南方电网有限责任公司 | Electric power system operation state detecting method based on dynamic partitioning |
| CN107016489A (en) * | 2017-03-09 | 2017-08-04 | 中国电力科学研究院 | A kind of electric power system robust state estimation method and device |
| CN107196306A (en) * | 2017-07-10 | 2017-09-22 | 大连海事大学 | Algorithm quicksort tidal current computing method based on Matlab sparse matrixes |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111711186A (en) | 2020-09-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111711186B (en) | Power system PQ decomposition state estimation method based on GPU parallel computing | |
| CN110110413B (en) | Structural topology optimization method based on material field reduction progression expansion | |
| CN105449675B (en) | The electric power networks reconstructing method of Optimum distribution formula energy access point and access ratio | |
| CN104933528B (en) | A kind of method that Jacobian matrix during electric power system tide calculates quickly is formed based on sparse matrix technology | |
| CN107591807B (en) | Optimization method for power transmission network planning under new energy access | |
| CN113420401A (en) | Optimal arrangement method for bias current blocking devices of power system | |
| CN107171365A (en) | Multiple target stochastic and dynamic economic load dispatching method with asynchronous iteration is decoupled based on scene | |
| CN113191105A (en) | Electrical simulation method based on distributed parallel operation method | |
| CN113659604A (en) | Electromechanical transient simulation method and device for LCC-VSC hybrid direct-current power grid and storage medium | |
| CN114336635B (en) | Full-pure embedded power flow calculation method and device based on constant term value and priori node | |
| CN111327048A (en) | Robust operation optimization method for power distribution network containing three-terminal SNOP | |
| CN113051796B (en) | Structural topology optimization design method applied to additive manufacturing | |
| CN107959287B (en) | Method for constructing two-voltage-level power grid growth evolution model | |
| CN109830987A (en) | The active distribution network Probabilistic Stability method of meter and distributed photovoltaic randomness | |
| CN116578830A (en) | Photovoltaic array reconstruction method based on improved marine predator algorithm | |
| CN111371125B (en) | Splitting and grouping optimization method for improving system coherence under condition of considering fan access | |
| CN113315118B (en) | Power system state estimation method based on parallel computing and particle swarm optimization | |
| CN108649585A (en) | Direct method for quickly searching static voltage stability domain boundary of power system | |
| CN111697607A (en) | Multi-terminal flexible direct-current transmission receiving-end power grid access method and system | |
| CN110930263B (en) | Medium-voltage distribution network short-circuit current calculation method containing photovoltaic power supply and induction motor based on black hole particle swarm algorithm | |
| CN114389261A (en) | Distributed photovoltaic power cluster division method considering grid-connected management mode | |
| CN114204613A (en) | Reactive compensation method and system for offshore wind farm access power system | |
| CN116341394B (en) | Hybrid driving model training method, device, computer equipment and storage medium | |
| CN115642597B (en) | Distributed photovoltaic bearing capacity calculation method and device | |
| CN113033024B (en) | Fine-grained parallel electromagnetic transient simulation method, system, terminal and medium for power transmission network |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |