CN113315118B - Power system state estimation method based on parallel computing and particle swarm optimization - Google Patents
Power system state estimation method based on parallel computing and particle swarm optimization Download PDFInfo
- Publication number
- CN113315118B CN113315118B CN202110455081.4A CN202110455081A CN113315118B CN 113315118 B CN113315118 B CN 113315118B CN 202110455081 A CN202110455081 A CN 202110455081A CN 113315118 B CN113315118 B CN 113315118B
- Authority
- CN
- China
- Prior art keywords
- node
- state estimation
- div
- reactive
- active
- 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
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 OR CALCULATING; COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Molecular Biology (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Power Engineering (AREA)
- Supply And Distribution Of Alternating Current (AREA)
Abstract
The invention provides a power system state estimation method based on parallel computing and a particle swarm algorithm, which overcomes the defect that the traditional power system state estimation method is not easy to converge, and simultaneously solves the problem that the operation time of the traditional power system state estimation optimizing algorithm is too long.
Description
Technical Field
The invention relates to the technical field of power system state estimation, in particular to a power system state estimation method based on parallel computing and particle swarm optimization.
Background
The power system state estimation refers to estimating the current operation state of the power system according to various measurement information of the power system. The state estimation of the power system is the basis of various high-grade applications in the power grid dispatching system and is the basis for 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 state estimation of the power system determines whether the power grid dispatching system can reflect the real state of the power grid in time, the control effect of other automatic power grid control systems is influenced, and the stable operation of the whole power system is very important.
At present, most of state estimation programs of actual power system operation are realized by taking measurement residual least squares as a target function based on a Newton Raphson method or a PQ decomposition method, the two methods are mainly applied to a convex optimization problem, but the actual power system is not necessarily a convex optimization problem completely, and under the condition, the problem of easy non-convergence exists, and the safety and stability analysis of the power system is greatly influenced. In 2010, limei, zhang Shiyuan, ren Weijian and other people published an article of "power system state estimation research based on particle swarm evolutionary algorithm" in the journal of power system protection and control, and the improved particle swarm evolutionary algorithm is applied to state estimation, so that the problem of non-convergence of traditional state estimation is solved, but the method is not combined with the traditional method, has long calculation time, cannot meet the requirement of a power system on real-time performance, and is only suitable for offline analysis.
Disclosure of Invention
In order to overcome the defect that the traditional power system state estimation method is not easy to converge and solve the problem that the operation time of the existing power system state estimation optimization algorithm is too long, the invention provides a power system state estimation method based on parallel computing and a particle swarm algorithm, so that the state estimation convergence is improved and the operation time is greatly reduced.
In order to achieve the technical effects, the technical scheme of the invention is as follows:
a power system state estimation method based on parallel computing and particle swarm optimization at least comprises the following steps:
s1, collecting measurement parameters of an electric power system and network parameters of the electric power system, and generating a branch admittance list according to the network parameters of the electric power system;
s2, respectively calculating an active Jacobian matrix H P And reactive Jacobian matrix H Q Active residual equation matrix A P And a reactive residual equation matrix A Q ；
S3, active residual equation matrix A P And reactive residual equation matrix A Q Respectively carrying out LLT decomposition to obtain an active lower triangular matrix A PL And a reactive lower triangular matrix A QL ；
S4, initializing node voltage amplitude and node voltage phase angle, carrying out state estimation iterative computation by a PQ decomposition method in a CPU, and simultaneously utilizing an active lower triangular matrix A in a GPU PL And a reactive lower triangular matrix A QL Calculating an active basis solution combination and a reactive basis solution combination;
s5, judging whether the state estimation iterative computation of the PQ decomposition method is converged, if so, finishing the iteration, and outputting a state estimation result; otherwise, confirming the unconverged node number after the last iterative computation of the PQ decomposition method, and executing the step S6;
s6, performing state estimation iterative optimization calculation on unconverged nodes in the GPU based on the active basis solution combination and the reactive basis solution combination and in combination with a particle swarm algorithm;
s7, judging whether the state estimation iterative optimization calculation is converged, if so, finishing the iterative optimization calculation, and outputting a state estimation result; otherwise, an error is prompted and the operation is exited.
According to the technical scheme, parameters required by state estimation of a traditional PQ decomposition method are prepared firstly, then power system state estimation is carried out in a CPU based on the traditional PQ decomposition method, iterative optimization calculation is carried out in a GPU based on a particle swarm algorithm when the traditional PQ decomposition method is not converged, state estimation convergence is improved, in addition, parallel calculation is carried out in the GPU simultaneously during iterative calculation of the traditional PQ decomposition method, the operation time of the particle swarm algorithm is greatly reduced, and the real-time performance of the power system state estimation is kept.
Preferably, the measurement parameters of the power system in step S1 include:
node voltage U, node injection active power P inj_M Node injection reactive power Q inj_M Active power P on two sides of branch ft_M 、P tf_M And reactive power Q at two sides of branch ft_M 、Q tf_M And current amplitude, wherein t and f respectively represent the serial numbers of the head end and the tail end of the branch; forming an active measurement value vector Z by measuring parameters P_M Vector Z of the measured value of the reactive power Q_M And an active power measurement weight matrix R P And a reactive power measurement weight matrix R Q The collected invalid or actual untransmitted measurement parameters are set to be 0, and the weight coefficient is also 0;
the network parameters of the power system include: node number nodeNum, branch number lnNum, branch first node number f, branch last node number t, branch resistance R, branch reactance X, branch non-standard transformation ratio K, ground admittance B and node parallel conductance G sh Node parallel susceptance B sh Node parallel conductance G without parallel element nodes sh And node parallel connectionNano B sh Is 0.
Preferably, the admittance Y of any branch k in said branch admittance list in step S1 s (k) The solving formula of (2) is as follows:
wherein, Y s (k) Representing the branch admittance, R (k) representing the branch resistance, X (k) representing the branch reactance; k (K) represents the branch nonstandard transformation ratio; g s (k) Conducting electricity to the branch; b s (k) Susceptance for the branch; j is the imaginary part identification in the complex number calculation.
Preferably, the active residual equation matrix a in step S2 P The calculation formula of (2) is as follows:
wherein,
is an active Jacobian matrix H P The transpose of (a) is performed,
is an active measurement weight matrix R P The inverse matrix of (d);
reactive residual equation matrix A Q The calculation formula of (c) is:
wherein,
is a reactive Jacobian matrix H Q The transpose of (a) is performed,
is a reactive power measurement weight matrix R Q Contrary to (2)A matrix;
active residual equation matrix A P Reactive residual equation matrix A Q All are positive definite matrixes, LLT decomposition is carried out to obtain the active lower triangular matrix A in the step S3 PL And a reactive lower triangular matrix A QL 。
Preferably, the active lower triangular matrix a is used in step S4 PL The method for calculating the active base solution combination comprises the following steps:
s401, according to a forward-backward substitution method, utilizing an active lower triangular matrix A PL Solving equation set A P ·x Pi ＝b Pi Wherein b is Pi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Pi Is the solution of the equation set;
s402, combining the solutions of the equation set to obtain an active radical solution combination G xP ＝{x P1 ，x P2 ，……，x PnodeNUm }；
Using a reactive lower triangular matrix A QL The method for calculating the reactive base solution combination comprises the following steps:
s411, according to the forward-backward substitution method, utilizing a reactive lower triangular matrix A QL Solving equation A Q ·x Qi ＝b Qi Wherein b is Qi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Qi Is the solution of the equation set;
s412, combining the solutions of the equation set to obtain a reactive base solution combination G xQ ＝{x Q1 ，x Q2 ，……，x QnodeNUm }；
The calculation processes of each base solution in the active base solution combination and the reactive base solution combination are independent from each other, and the solving process of each base solution is processed in a single GPU thread.
Preferably, in step S4, the node voltage amplitude is initialized to 1.0pu, the node voltage phase angle is initialized to 0.0rad, and the iterative calculation of PQ decomposition state estimation performed in the CPU and the GPU using the active lower triangular matrix A PL And a reactive lower triangular matrix A QL The process of calculating the active base solution combination and the reactive base solution combination is carried out simultaneously, and the active lower triangular moment is utilizedArray A PL And a reactive lower triangular matrix A QL And calculating the power base solution combination and the reactive base solution combination, preparing for carrying out state estimation iterative optimization calculation on the non-convergence nodes in the GPU, reducing the operation time and ensuring the real-time performance of the state estimation of the power system.
Preferably, the process of confirming the unconverged node number after the last iterative computation of the PQ decomposition method in step S5 is as follows:
confirming a convergence threshold E of the PQ decomposition method, and extracting a node voltage amplitude correction quantity and a node voltage phase angle correction quantity after the last iterative calculation of the PQ decomposition method;
and taking the node number when the node voltage amplitude correction quantity and the node voltage phase angle correction quantity are both larger than the convergence threshold E as the unconverged node number after the last iterative calculation of the PQ decomposition method.
Preferably, the step S6, in combination with the particle swarm algorithm, of performing state estimation iterative optimization calculation on the non-converged node in the GPU includes:
s61, taking the node voltage amplitude and the node voltage phase angle of the unconverged node after the last iterative computation of the PQ decomposition method as initial values, setting the iterative optimization-seeking computation times of state estimation as m1 and the upper limit value as m1max; setting the iteration number of the particle swarm algorithm as m2, setting the upper limit value as m2max, and initializing both m1 and m2 as 0;
s62, performing state estimation iterative optimization calculation by taking the node voltage amplitude and the node voltage phase angle of the unconverged node as state variables and taking the total network measurement error least square in the particle swarm optimization as a target function;
s63, judging whether the particles are converged in the state estimation iterative optimization calculation, if so, executing a step S64, otherwise, returning to the step S62;
s64, judging whether the iteration number m2 of the particle swarm algorithm reaches an upper limit value m2max, if so, taking an optimal state, and confirming the node voltage amplitude state correction dU of the unconverged node corresponding to the optimal state div Sum node voltage phase angle state correction d theta div Step S65 is executed; otherwise, returning to step S62;
s65, judging the node voltage of the unconverged nodeAmplitude state correction dU div Sum node voltage phase angle state correction d theta div Whether a state estimation convergence condition is met or not is judged, and if yes, a result is output; otherwise, go to step S66;
s66, judging whether the state estimation iteration optimizing calculation times m1 reach an upper limit value m1max, and if so, prompting an error exit; otherwise, correcting the node voltage amplitude state of the non-convergence node by a correction amount dU div Sum node voltage phase angle state correction d theta div And (4) superimposing the state variable on the state variable, and returning to the step S62.
Preferably, in step S62, the process of performing state estimation iterative optimization calculation by using the total network measurement error least square in the particle swarm algorithm as the objective function is as follows:
s621, setting m3 particles, and generating a non-convergence node voltage velocity vector v for each particle o u (o) elements of which obey evenly distributed random quantities, the length being the number of unconverged node voltages; generating an unconverged node voltage phase angle velocity vector v θ (o) its elements are random quantities subject to an even distribution, the length being the number of unconverged phase angles;
s622, after the last iterative computation of the PQ decomposition method, the node voltage amplitude of the unconverged node and the voltage velocity vector v of the unconverged node u (o) adding to obtain a voltage vector u of the particle div (o); after the last iterative computation of the PQ decomposition method, the node voltage phase angle of the unconverged node and the speed vector v of the unconverged node voltage phase angle θ (o) adding to obtain a voltage vector u of the particle div (o)；
S623, parallel computing the measurement estimation value, the measurement error and the objective function value corresponding to each particle in the GPU, selecting the particle with the minimum objective function, and recording the serial number p of the particle best ；
Calculating a direction vector of each particle
d u ＝u div (o)-u div (p best )
d θ ＝θ div (o)-θ div (p best )
Updating the velocity vector of each particle
v u (o)＝ω·v u (o)+l 1 ·rand(0,1)d u
v θ (o)＝ω·v θ (o)+l 2 ·rand(0,1)d θ
Wherein w is a non-negative inertia factor, l 1 、l 2 For the acceleration constant, rand (0,1) is a random number between 0 and 1;
updating the vector of voltage values for each particle
u div (o)＝u div (o)+v u (o)
θ div (o)＝θ div (o)+θ u (o)
S623, if the direction vectors d of all the particles u Is less than 3 times of node voltage amplitude value d θ If the maximum value of the phase angle is less than 3 times of the node voltage phase angle, the particle swarm iteration is optimized and converged, and the node voltage amplitude u is output div (pbest), node voltage phase angle θ div (pbest); otherwise, the particle swarm iterative optimization is not converged, and step S624 is executed;
s624, judging whether the iteration number m2 of the particle swarm algorithm reaches an upper limit value m2max, if so, optimizing the particle swarm iteration without convergence, and prompting an error exit; otherwise, increasing the iteration number of the particle swarm algorithm by one, and returning to the step C622.
The number m3 of the particles is set according to an operation unit of the GPU, the range of the elements of the voltage velocity vector vu (o) of the unconverged node subject to the average distribution is-0.2, and the range of the elements of the voltage phase angle velocity vector v theta (o) of the unconverged node subject to the average distribution is-0.3.
Preferably, step S623 outputs a node voltage magnitude u div (pbest), node voltage phase angle θ div (pbest) after which the output u is used div (pbest)、θ div (pbest) calculation measurement estimation value and active and reactive measurement residual error r P 、r Q (ii) a Based on this, the variation of the optimal particle is solved:
wherein, b P 、b Q Is a residual vector, b Pi 、b Qi Are respectively b P 、b Q The ith element of (1), x Pi ,x Qi Is the ith base solution vector of the active and reactive base solution combination.
Judging max (abs (d θ) div ) Max (abs (dU)) div ) Whether the convergence criterion of the node voltage amplitude and the node voltage phase angle in the state estimation of the PQ decomposition method is met or not respectively, if yes, the convergence is considered, and a state estimation result is output; otherwise, judging the number of iterative optimization searching calculation times by entering state estimation:
if the state estimation iterative optimization searching calculation time m1 does not reach the upper limit value m1max, increasing one by the state estimation iterative optimization searching calculation time m1, and returning to the state estimation iterative optimization searching calculation which is carried out by taking the whole network measurement error least square in the particle swarm optimization as a target function; otherwise, the convergence is not performed, and an error exit is prompted.
Compared with the prior art, the technical scheme of the invention has the beneficial effects that:
the invention provides a power system state estimation method based on parallel computing and a particle swarm algorithm.
Drawings
Fig. 1 is a schematic flow chart of a power system state estimation method based on parallel computing and particle swarm optimization according to an embodiment of the present invention.
Detailed Description
The drawings are for illustrative purposes only and are not to be construed as limiting the patent;
for better illustration of the present embodiment, certain parts of the drawings may be omitted, enlarged or reduced, and do not represent actual dimensions;
it will be understood by those skilled in the art that certain well-known descriptions of the figures may be omitted.
The positional relationships depicted in the drawings are for illustrative purposes only and are not to be construed as limiting the present patent;
the technical solution of the present invention is further described with reference to the drawings and the embodiments.
Example 1
The flow chart of the power system state estimation method based on parallel computing and particle swarm optimization as shown in fig. 1 is schematic, and referring to fig. 1, the method comprises the following steps:
s1, collecting measurement parameters of an electric power system and network parameters of the electric power system, and generating a branch admittance list according to the network parameters of the electric power system;
in specific implementation, measurement parameters of the power system are collected by a data collection system, network parameters of the power system are collected by an Energy Management System (EMS), and the measurement parameters of the power system include:
node voltage U, node injection active power P inj_M Node injection reactive power Q inj_M Active power P on two sides of branch ft_M 、P tf_M Reactive power Q at two sides of branch ft_M 、Q tf_M And current amplitude, wherein t and f respectively represent the serial numbers of the head end and the tail end of the branch; forming an active measurement value vector Z by measuring parameters P_M Vector Z of the measured value of the reactive power Q_M And an active measurement weight matrix R P And a reactive power measurement weight matrix R Q The collected invalid or actual untransmitted measurement parameters are set to be 0, and the weight coefficient is also 0;
the network parameters of the power system include: node number nodeNum, branch number lnNum, branch first node number f, branch last node number t, branch resistance R, branch reactance X, branch non-standard transformation ratio K, ground admittance B and node parallel conductance G sh Node parallel susceptance B sh Node parallel conductance G without parallel element nodes sh And node parallel susceptance B sh Is 0, the admittance Y of any branch k in the branch admittance list s (k) The solving formula of (2) is as follows:
wherein, Y s (k) Representing the branch admittance, R (k) representing the branch resistance, X (k) representing the branch reactance; k (K) represents the branch nonstandard transformation ratio; g is a radical of formula s (k) Conducting electricity to the branch; b s (k) Susceptance for the branch; j is the imaginary part identification in the complex number calculation, and the acquisition and processing of the above related parameters are all necessary preparations for the subsequent state estimation.
S2, respectively calculating active power Jacobian matrix H P Reactive Jacobian matrix H Q Active residual equation matrix A P And a reactive residual equation matrix A Q ；
In the present embodiment, the active Jacobian matrix H P The calculation formula of the elements in (1) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein, P inj (i) Active power is injected for point i; b ij Representing the conductance of the branch between node i and node j, B ii Representing the self-conductance, theta, of node i i Is the phase angle of the voltage at node i, θ j Is the voltage phase angle of node j; u shape f(k) Is the voltage amplitude of the starting node of branch k; u shape t(k) The voltage amplitude of the k end node of the branch circuit; theta f(k) Phase angle of voltage at starting node of branch k, theta t(k) Is the voltage phase angle of the k end node of the branch. f (k) is the number of the initial node of the branch k, t (k) is the number of the last node of the branch k, and theta ft Representing the voltage phase angle difference between node i and node j.
Reactive Jacobian matrix H Q Sparse matrix of nodeNum + 2+ lnNum 2 rows of nodeNum columns, reactive Jacobian matrix H Q The calculation formula of the elements in (1) is as follows:
θ ft ＝θ f(k) -θ t(k)
wherein Q is inj (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 Is the voltage amplitude, U, of node i j Is the voltage amplitude of node j; u shape f(k) Is the voltage amplitude of the starting node of branch k; u shape t(k) The voltage amplitude of the k end node of the branch circuit; theta f(k) Phase angle of voltage at starting node of branch k, theta t(k) Is the voltage phase angle of the k end node of the branch. f (k) is the number of the initial node of the branch k, t (k) is the number of the last node of the branch k, and theta ft Representing the voltage phase angle difference between node i and node j.
Active residual equation matrix A P The calculation formula of (2) is as follows:
wherein,
is an active Jacobian matrix H P The transpose of (a) is performed,
is an active measurement weight matrix R P The inverse matrix of (d); reactive residual equation matrix A Q The calculation formula of (2) is as follows:
wherein,
is a reactive Jacobian matrix H Q The transpose of (a) is performed,
measure the weight matrix R for the reactive power Q The inverse matrix of (c). Due to the active power measurement weight matrix R P And a reactive power measurement weight matrix R Q Are diagonal arrays, and when collecting the measurement parameters of the power system, the collected invalid or actually untransmitted measurement parameters are set to 0, so that the active measurement weight matrix R P And a reactive power measurement weight matrix R Q In some cases, the active measurement weight matrix R is set to 0, and there may be a case without an inverse matrix P And a reactive power measurement weight matrix R Q The medium non-zero elements are inverted to obtain equivalent inversion.
S3, to active residual equation matrix A P And reactive residual equation matrix A Q Respectively carrying out LLT decomposition to obtain an active lower triangular matrix A PL And a reactive lower triangular matrix A QL (ii) a In specific implementation, the active residual equation matrix A P Reactive residual equation matrix A Q All are positive definite matrixes, LLT decomposition is carried out to obtain the active lower triangular matrix A in the step S3 PL And a reactive lower triangular matrix A QL For active residual equation matrix A P And reactive residual equation matrix A Q Performing LLT decomposition separately is a relatively mature conventional technique.
S4, initializing node voltage amplitude and node voltage phase angle, and performing PQ decomposition in a CPUAnd (3) state estimation iterative calculation, specifically, initializing the node voltage amplitude to 1.0pu, initializing the node voltage phase angle to 0.0rad, and using an active lower triangular matrix A in the step S4 PL The method for calculating the active base solution combination comprises the following steps:
s401, according to a forward-backward substitution method, utilizing an active lower triangular matrix A PL Solving equation set A P ·x Pi ＝b Pi Wherein b is Pi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Pi Is the solution of the equation set;
s402, combining the solutions of the equation set to obtain an active radical solution combination G xP ＝{x P1 ，x P2 ，……，x PnodeNUm }；
Using a reactive lower triangular matrix A QL The method for calculating the reactive base solution combination comprises the following steps:
s411, according to the forward-backward substitution method, utilizing a reactive lower triangular matrix A QL Solving equation A Q ·x Qi ＝b Qi Wherein b is Qi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Qi Is the solution of the equation set;
s412, combining the solutions of the equation set to obtain a reactive base solution combination G xQ ＝{x Q1 ，x Q2 ，……，x QnodeNUm }；
The calculation processes of each base solution in the active base solution combination and the reactive base solution combination are independent from each other, and the solving process of each base solution is processed in a single GPU thread.
S5, judging whether the state estimation iterative computation of the PQ decomposition method is converged, if so, finishing the iteration, and outputting a state estimation result; otherwise, confirming the unconverged node number after the last iterative computation of the PQ decomposition method, and executing the step S6; the process of state estimation by the PQ decomposition method is not repeated herein, and in specific implementation, if the PQ decomposition method does not converge, the node voltage amplitude correction and the node voltage phase angle correction after the last iteration calculation by the PQ decomposition method are extracted according to the convergence threshold E of the PQ decomposition method; and taking the node number when the node voltage amplitude correction quantity and the node voltage phase angle correction quantity are both larger than the convergence threshold E as the unconverged node number after the last iterative calculation of the PQ decomposition method.
S6, performing state estimation iterative optimization calculation on unconverged nodes in the GPU based on the active basis solution combination and the reactive basis solution combination and in combination with a particle swarm algorithm;
s7, judging whether the state estimation iterative optimization calculation is converged, if so, finishing the iterative optimization calculation, and outputting a state estimation result; otherwise, an error is prompted and the operation is exited.
In this embodiment, in combination with the particle swarm algorithm, the process of performing state estimation iterative optimization calculation on the non-converged node in the GPU in step S6 includes:
s61, taking the node voltage amplitude and the node voltage phase angle of the unconverged node after the last iterative computation of the PQ decomposition method as initial values, setting the iterative optimization-seeking computation times of state estimation as m1 and the upper limit value as m1max; setting the iteration number of the particle swarm algorithm as m2, setting the upper limit value as m2max, and initializing both m1 and m2 as 0;
s62, performing state estimation iterative optimization calculation by taking the node voltage amplitude and the node voltage phase angle of the unconverged node as state variables and taking the total network measurement error least square in the particle swarm optimization as a target function;
s63, judging whether the particles are converged in the state estimation iterative optimization calculation, if so, executing a step S64, otherwise, returning to the step S62;
s64, judging whether the iteration number m2 of the particle swarm algorithm reaches an upper limit value m2max, if so, taking an optimal state, and confirming the node voltage amplitude state correction dU of the unconverged node corresponding to the optimal state div Sum node voltage phase angle state correction d theta div Step S65 is executed; otherwise, returning to step S62;
s65, judging node voltage amplitude state correction dU of unconverged node div Sum node voltage phase angle state correction d theta div Whether a state estimation convergence condition is met or not is judged, and if yes, a result is output; otherwise, go to step S66;
s66, judging whether the iterative optimization calculation times m1 of state estimation reaches an upper limit value m1max, if yes, prompting an error exit; otherwise, correcting the node voltage amplitude state of the unconverged node by a correction amount dU div Sum node voltage phase angle state correction d theta div And (4) superimposing the state variable on the state variable, and returning to the step S62.
In step S62, the process of performing state estimation iterative optimization calculation with the least square of the measurement error of the whole network in the particle swarm optimization as the objective function is as follows:
s621, setting m3 particles, and generating a non-convergence node voltage velocity vector v for each particle o u (o) elements of which obey evenly distributed random quantities, the length being the number of unconverged node voltages; generating an unconverged node voltage phase angle velocity vector v θ (o) its elements are random quantities subject to an even distribution, the length being the number of unconverged phase angles; the number m3 of the particles is set according to an operation unit of the GPU, the range of the elements of the voltage velocity vector vu (o) of the unconverged node subject to the average distribution is-0.2, and the range of the elements of the voltage phase angle velocity vector v theta (o) of the unconverged node subject to the average distribution is-0.3.
S622, after the last iterative computation of the PQ decomposition method, the node voltage amplitude of the unconverged node and the voltage velocity vector v of the unconverged node u (o) adding to obtain a voltage vector u of the particle div (o); after the last iterative computation of the PQ decomposition method, the node voltage phase angle of the unconverged node and the speed vector v of the unconverged node voltage phase angle θ (o) adding to obtain a voltage vector u of the particle div (o)；
S623, parallel computing the measurement estimation value, the measurement error and the objective function value corresponding to each particle in the GPU, selecting the particle with the minimum objective function, and recording the serial number p of the particle best ；
Calculating a direction vector of each particle
d u ＝u div (o)-u div (p best )
d θ ＝θ div (o)-θ div (p best )
Updating the velocity vector of each particle
v u (o)＝ω·v u (o)+l 1 ·rand(0,1)d u
v θ (o)＝ω·v θ (o)+l 2 ·rand(0,1)d θ
Wherein w is a non-negative inertia factor, l 1 、l 2 For the acceleration constant, rand (0,1) is a random number between 0 and 1;
updating the vector of voltage values for each particle
u div (o)＝u div (o)+v u (o)
θ div (o)＝θ div (o)+θ u (o)
S623, if the direction vectors d of all the particles u Is less than 3 times of node voltage amplitude value d θ The maximum value of the output node is less than 3 times of the node voltage phase angle, the particle swarm iteration is optimized and converged, and the node voltage amplitude u is output div (pbest), node voltage phase angle θ div (pbest); otherwise, the particle swarm iterative optimization is not converged, and step S624 is executed;
s624, judging whether the iteration number m2 of the particle swarm algorithm reaches an upper limit value m2max, if so, optimizing the particle swarm iteration without convergence, and prompting an error exit; otherwise, increasing the iteration number of the particle swarm algorithm by one, and returning to the step C622.
In the present embodiment, the step S623 outputs the node voltage amplitude u div (pbest), node voltage phase angle θ div (pbest) followed by the use of the output u div (pbest)、θ div (pbest) calculation measurement estimation value and active and reactive measurement residual error r P 、r Q (ii) a Based on this, the variation of the optimal particle is solved:
wherein, b P 、b Q Is a residual vector, b Pi 、b Qi Are respectively b P 、b Q The ith element of (1), x Pi ,x Qi Is the ith base solution vector of the active and reactive base solution combination.
Judging max (abs (d θ) div ) Max (abs (dU)) div ) Whether the convergence criteria of the node voltage amplitude and the node voltage phase angle in the state estimation of the PQ decomposition method are met respectively, if yes, convergence is considered, and a state estimation result is output; otherwise, judging the number of iterative optimization searching calculation times by entering state estimation:
if the state estimation iterative optimization searching calculation time m1 does not reach the upper limit value m1max, increasing one by the state estimation iterative optimization searching calculation time m1, and returning to the state estimation iterative optimization searching calculation which is carried out by taking the whole network measurement error least square in the particle swarm optimization as a target function; otherwise, the convergence is not achieved, and an error exit is prompted.
The positional relationships depicted in the drawings are for illustrative purposes only and are not to be construed as limiting the present patent;
it should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.
Claims (9)
1. A power system state estimation method based on parallel computing and particle swarm optimization 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, and generating a branch admittance list according to the network parameters of the electric power system;
s2, respectively calculating an active Jacobian matrix H P Reactive Jacobian matrix H Q Active residual equation matrix A P And a reactive residual equation matrix A Q ；
S3, to active residual equation matrix A P And reactive residual equation matrix A Q Respectively carrying out LLT decomposition to obtain an active lower triangular matrix A PL And a reactive lower triangular matrix A QL ；
S4, initializing a node voltage amplitude value and a node voltage phase angle, and performing state estimation iterative calculation of a PQ decomposition method in a CPU; meanwhile, an active lower triangular matrix A is utilized in the GPU PL And a reactive lower triangular matrix A QL Calculating an active basis solution combination and a reactive basis solution combination;
in step S4, an active lower triangular matrix A is utilized PL The method for calculating the active base solution combination comprises the following steps:
s401, according to a forward-backward substitution method, utilizing an active lower triangular matrix A PL Solving equation set A P ·x Pi ＝b Pi Wherein b is Pi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Pi Is the solution of the equation set;
s402, combining solutions of the equation set to obtain an active radical solution combination G xP ＝{x P1 ，x P2 ，……，x PnodeNUm }；
Using a reactive lower triangular matrix A QL The method for calculating the reactive base solution combination comprises the following steps:
s411, according to the forward-backward substitution method, utilizing a reactive lower triangular matrix A QL Solving equation A Q ·x Qi ＝b Qi Wherein b is Qi Is length nodeNum, is the column vector with the ith element being 1 and the remaining elements being 0; x is the number of Qi Is the solution of the equation set;
s412, combining the solutions of the equation set to obtain a reactive base solution combination G xQ ＝{x Q1 ，x Q2 ，……，x QnodeNUm }；
The calculation process of each base solution in the active base solution combination and the reactive base solution combination is independent, and the solving process of each base solution is processed in a single GPU thread;
s5, judging whether the state estimation iterative computation of the PQ decomposition method is converged, if so, finishing the iteration, and outputting a state estimation result; otherwise, confirming the unconverged node number after the last iterative computation of the PQ decomposition method, and executing the step S6;
s6, performing state estimation iterative optimization calculation on unconverged nodes in the GPU based on the active basis solution combination and the reactive basis solution combination and in combination with a particle swarm algorithm;
s7, judging whether the state estimation iterative optimization calculation is converged, if so, finishing the iterative optimization calculation, and outputting a state estimation result; otherwise, an error is prompted and the operation is exited.
2. The method for estimating the state of a power system based on parallel computing and particle swarm optimization according to claim 1, wherein the measurement parameters of the power system in step S1 comprise:
node voltage U, node injection active power P inj_M Node injection reactive power Q inj_M Active power P on two sides of branch ft_M 、P tf_M Reactive power Q at two sides of branch ft_M 、Q tf_M And current amplitude, wherein t and f respectively represent the serial numbers of the head end and the tail end of the branch; forming an active measurement value vector Z by measuring parameters P_M Vector Z of the measured value of the dead time Q_M And an active power measurement weight matrix R P And a reactive power measurement weight matrix R Q The collected invalid or actual untransmitted measurement parameters are set to be 0, and the weight coefficient is also 0;
the network parameters of the power system include: node number nodeNum, branch number lnNum, branch first node number f, branch last node number t, branch resistance R, branch reactance X, branch non-standard transformation ratio K, ground admittance B and node parallel conductance G sh Node parallel susceptance B sh Node parallel conductance G without parallel element nodes sh And node parallel susceptance B sh Is 0.
3. The parallel computing and particle swarm algorithm-based power system state estimation method according to claim 2, wherein in step S1, the admittance Y of any branch k in the branch admittance list is s (k) The solving formula of (2) is as follows:
wherein Y is s (k) Representing the branch admittance, R (k) representing the branch resistance, X (k) representing the branch reactance; k (K) represents the branch nonstandard transformation ratio; g s (k) Conducting electricity to the branch; b is a mixture of s (k) Susceptance for the branch; j is the imaginary part identification in the complex number calculation.
4. The parallel computing and particle swarm algorithm-based power system state estimation method according to claim 3, wherein the active residual equation matrix A in the step S2 P The calculation formula of (2) is as follows:
wherein,
is an active Jacobian matrix H P The transpose of (a) is performed,
is an active measurement weight matrix R P The inverse matrix of (d);
reactive residual equation matrix A Q The calculation formula of (2) is as follows:
wherein，
Is a reactive Jacobian matrix H Q The transpose of (a) is performed,
is a reactive power measurement weight matrix R Q The inverse matrix of (d);
active residual equation matrix A P Reactive residual equation matrix A Q All are positive definite matrixes, LLT decomposition is carried out to obtain the active lower triangular matrix A in the step S3 PL And a reactive lower triangular matrix A QL 。
5. The parallel computing and particle swarm algorithm-based power system state estimation method according to claim 4, wherein in step S4, the node voltage amplitude is initialized to 1.0pu, the node voltage phase angle is initialized to 0.0rad, and the PQ decomposition state estimation iterative computation performed in CPU and GPU utilize the active lower triangular matrix A PL And a reactive lower triangular matrix A QL The process of calculating the active and reactive base solution combinations is carried out simultaneously, and an active lower triangular matrix A is utilized PL And a reactive lower triangular matrix A QL And calculating a power base solution combination and a reactive base solution combination, and preparing for performing state estimation iterative optimization calculation on the non-convergence nodes in the GPU.
6. The method for estimating state of power system based on parallel computing and particle swarm optimization according to claim 5, wherein the step S5 of confirming the unconverged node number after the last iteration of the PQ decomposition method comprises:
confirming a convergence threshold E of the PQ decomposition method, and extracting a node voltage amplitude correction quantity and a node voltage phase angle correction quantity after the last iteration calculation of the PQ decomposition method;
and taking the node number when the node voltage amplitude correction quantity and the node voltage phase angle correction quantity are both larger than a convergence threshold E as the node number which is not converged after the last iterative calculation of the PQ decomposition method.
7. The method for estimating the state of the power system based on the parallel computing and the particle swarm optimization algorithm according to claim 6, wherein the step S6 of performing the state estimation iterative optimization calculation on the non-converged nodes in the GPU in combination with the particle swarm optimization algorithm comprises:
s61, taking the node voltage amplitude and the node voltage phase angle of the unconverged node after the last iterative computation of the PQ decomposition method as initial values, setting the iterative optimization-seeking computation times of state estimation as m1 and the upper limit value as m1max; setting the iteration number of the particle swarm algorithm as m2, setting the upper limit value as m2max, and initializing both m1 and m2 as 0;
s62, performing state estimation iterative optimization calculation by taking the node voltage amplitude and the node voltage phase angle of the unconverged node as state variables and taking the total network measurement error least square in the particle swarm optimization as a target function;
s63, judging whether the particles are converged in the state estimation iterative optimization calculation, if so, executing a step S64, otherwise, returning to the step S62;
s64, judging whether the iteration number m2 of the particle swarm algorithm reaches an upper limit value m2max, if so, taking an optimal state, and confirming the node voltage amplitude state correction dU of the unconverged node corresponding to the optimal state div Sum node voltage phase angle state correction d theta div Step S65 is executed; otherwise, returning to step S62;
s65, judging node voltage amplitude state correction dU of unconverged node div Sum node voltage phase angle state correction d theta div Whether a state estimation convergence condition is met or not is judged, and if yes, a result is output; otherwise, go to step S66;
s66, judging whether the state estimation iteration optimizing calculation times m1 reach an upper limit value m1max, and if so, prompting an error exit; otherwise, correcting the node voltage amplitude state of the unconverged node by a correction amount dU div Sum node voltage phase angle state correction d theta div And (4) superimposing the state variable on the state variable, and returning to the step S62.
8. The method for estimating the state of the power system based on the parallel computing and the particle swarm optimization according to claim 7, wherein the iterative state estimation optimization computing process with the least square of the measurement error of the whole network in the particle swarm optimization as the objective function in step S62 comprises:
s621, setting m3 particles, and generating a non-convergence node voltage velocity vector v for each particle o u (o) elements of which obey evenly distributed random quantities, the length being the number of unconverged node voltages; generating an unconverged node voltage phase angle velocity vector v θ (o) its elements are random quantities subject to an even distribution, the length being the number of unconverged phase angles;
s622, after the last iterative computation of the PQ decomposition method, the node voltage amplitude of the unconverged node and the voltage velocity vector v of the unconverged node u (o) adding to obtain a voltage vector u of the particle div (o); after the last iterative computation of the PQ decomposition method, the node voltage phase angle of the unconverged node and the speed vector v of the unconverged node voltage phase angle θ (o) adding to obtain a voltage phase angle vector θ of the particle div (o)；
S623, parallel computing the measurement estimation value, the measurement error and the objective function value corresponding to each particle in the GPU, selecting the particle with the minimum objective function, and recording the serial number p of the particle best ；
Calculating a direction vector of each particle
d u ＝u div (o)-u div (p best )
d θ ＝θ div (o)-θ div (p best )
Updating the velocity vector of each particle
v u (o)＝ω·v u (o)+l 1 ·rand(0，1)d u
v θ (o)＝ω·v θ (o)+l 2 ，rand(0，1)d θ
Where ω is a non-negative inertia factor, l 1 、l 2 For the acceleration constant, rand (0,1) is a random number between 0 and 1;
updating the vector of voltage values for each particle
u div (o)＝u div (o)+v u (o)
θ div (o)＝θ div (o)+θ u (o)
S623, if the direction vectors d of all the particles u Is less than 3 times of node voltage amplitude value d θ The maximum value of the output node is less than 3 times of the node voltage phase angle, the particle swarm iteration is optimized and converged, and the node voltage amplitude u is output div (pbest), node voltage phase angle θ div (pbest); otherwise, the particle swarm iterative optimization is not converged, and step S624 is executed;
s624, judging whether the iteration times m2 of the particle swarm algorithm reaches an upper limit value m2max, if yes, optimizing the particle swarm iteration and failing to converge, and prompting an error exit; otherwise, increasing the iteration times of the particle swarm algorithm by one, and returning to the step C622;
the number m3 of the particles is set according to an operation unit of the GPU, the range of the elements of the voltage velocity vector vu (o) of the unconverged node subject to the average distribution is-0.2, and the range of the elements of the voltage phase angle velocity vector v theta (o) of the unconverged node subject to the average distribution is-0.3.
9. The parallel computing and particle swarm algorithm-based power system state estimation method according to claim 8, wherein step S623 outputs a node voltage amplitude u div (pbest), node voltage phase angle θ div (pbest) followed by the use of the output u div (pbest)、θ div (pbest) calculation measurement estimation value and active and reactive measurement residual error r P 、r Q ；
Based on this, the variation of the optimal particle is solved:
wherein, b P 、b Q Is a residual vector, b Pi 、b Qi Are respectively b P 、b Q The ith element of (1), x Pi ，x Qi The ith base solution vector is the combination of active and reactive base solutions;
judging max (abs (d θ) div ) Max (abs (dU)) div ) Whether the convergence criterion of the node voltage amplitude and the node voltage phase angle in the state estimation of the PQ decomposition method is met or not respectively, if yes, the convergence is considered, and a state estimation result is output; otherwise, judging the number of iterative optimization searching calculation times by entering state estimation:
if the state estimation iterative optimization searching calculation time m1 does not reach the upper limit value m1max, increasing one by the state estimation iterative optimization searching calculation time m1, and returning to the state estimation iterative optimization searching calculation which is carried out by taking the whole network measurement error least square in the particle swarm optimization as a target function; otherwise, the convergence is not achieved, and an error exit is prompted.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110455081.4A CN113315118B (en) | 2021-04-26 | 2021-04-26 | Power system state estimation method based on parallel computing and particle swarm optimization |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110455081.4A CN113315118B (en) | 2021-04-26 | 2021-04-26 | Power system state estimation method based on parallel computing and particle swarm optimization |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113315118A CN113315118A (en) | 2021-08-27 |
| CN113315118B true CN113315118B (en) | 2022-12-30 |
Family
ID=77371172
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110455081.4A Active CN113315118B (en) | 2021-04-26 | 2021-04-26 | Power system state estimation method based on parallel computing and particle swarm optimization |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113315118B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115202992B (en) * | 2022-09-15 | 2022-11-22 | 中国空气动力研究与发展中心计算空气动力研究所 | CFD operation convergence monitoring method for slurm scheduling system |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107425520A (en) * | 2017-06-12 | 2017-12-01 | 东南大学 | A kind of probabilistic active distribution network three-phase section method for estimating state of injecting power containing node |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9645558B2 (en) * | 2012-09-29 | 2017-05-09 | Operation Technology, Inc. | Dynamic parameter tuning using particle swarm optimization |
| CN105699849A (en) * | 2016-03-18 | 2016-06-22 | 国网上海市电力公司 | A voltage sag estimation method based on quantum behavior particle swarm optimization algorithm |
| CN107528312B (en) * | 2017-07-14 | 2020-04-28 | 同济大学 | Power system state estimation method |
| CN108462184B (en) * | 2018-03-26 | 2021-06-18 | 昆明理工大学 | An optimal configuration method for series compensation of power system lines |
| CN110504676B (en) * | 2019-07-05 | 2021-01-26 | 国网浙江平湖市供电有限公司 | Power distribution network state estimation method based on APSO-BP |
| CN111711186B (en) * | 2020-05-28 | 2023-05-02 | 中国南方电网有限责任公司 | Power system PQ decomposition state estimation method based on GPU parallel computing |
| CN112149273B (en) * | 2020-08-18 | 2024-06-18 | 浙江工业大学 | Fast-convergence alternating-current power grid distributed state estimation method |
-
2021
- 2021-04-26 CN CN202110455081.4A patent/CN113315118B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107425520A (en) * | 2017-06-12 | 2017-12-01 | 东南大学 | A kind of probabilistic active distribution network three-phase section method for estimating state of injecting power containing node |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113315118A (en) | 2021-08-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111711186B (en) | Power system PQ decomposition state estimation method based on GPU parallel computing | |
| Saidi et al. | Precise parameter identification of a PEMFC model using a robust enhanced salp swarm algorithm | |
| CN105512502B (en) | One kind is based on the normalized weight function the least square estimation method of residual error | |
| CN107994567B (en) | A Generalized Fast Decomposition State Estimation Method | |
| CN113315118B (en) | Power system state estimation method based on parallel computing and particle swarm optimization | |
| CN105116346A (en) | Series-connected battery system and method for estimating state of charge thereof | |
| CN103094900A (en) | Distributed generation power distribution network three-phase load flow calculation method taking phase sequence mixing method into consideration | |
| CN108462184B (en) | An optimal configuration method for series compensation of power system lines | |
| CN109038545B (en) | A Distribution Network Reconfiguration Method Based on Differential Evolution Invasion Weed Algorithm | |
| CN112564084B (en) | Method for rapidly determining voltage stability of large-scale distributed photovoltaic access power distribution network | |
| WO2025162502A1 (en) | Fuel cell parameter identification method and apparatus, and device and storage medium | |
| CN110601172A (en) | Multi-direct-current coordination controller design method based on convex polyhedron uncertainty | |
| CN105514977B (en) | A kind of hyperbolic cosine type robust state estimation method of POWER SYSTEM STATE | |
| CN105914752B (en) | Pilot bus selection method based on fast search and density peaks cluster | |
| CN108923431A (en) | The AC-DC hybrid power grid of the power grid containing flexible direct current decouples tidal current computing method | |
| CN108988393B (en) | Micro-source time sequence optimization method for micro-grid black start | |
| CN113222095B (en) | Fractional order prediction auxiliary state estimation method for power system based on evolutionary computation | |
| CN109698505A (en) | The regulation of the online prevention and control of bulk power grid static voltage stability quantifies mapping calculation method | |
| CN113506896A (en) | Predictive control method of SOFC temperature distribution based on space-time separation modeling | |
| CN118137592A (en) | Reactive power optimization method for distributed photovoltaic access distribution network based on cluster division | |
| CN105303471B (en) | A Hyperbolic Cosine Maximum Exponential Square Robust Estimation Method for Power System State | |
| Djari et al. | Voltage stability assessment using different methods | |
| Wang et al. | Calculation of maximum permitted capacity of photovoltaic based on flower pollination algorithm combined with genetic algorithm (GA-FPA) in distribution network | |
| CN113162063A (en) | Design method of multi-direct-current coordination controller for inhibiting ultralow frequency oscillation | |
| Wang et al. | Identification of thermal process using fractional-order transfer function based on intelligent optimization |
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 |