Intelligent power distribution network robust dynamic state estimation method based on two-cycle data association

Abstract

The invention relates to a robust dynamic state estimation method for a two-cycle data associated intelligent power distribution network, which relates to the field of power distribution networks and comprises the following steps of: 1, carrying out measurement transformation on the power distribution network, establishing a unified measurement model and turning to the step 2; 2, ultra-short-term power prediction is carried out, parallel acceleration is carried out on matrix calculation, and the step 3 is carried out; 3, judging whether the sampling time is an SCADA sampling time or not, if so, forming a full-measurement model by using mu M-PMU and SCADA measurement, and turning to the step 4; otherwise, forming a linear measurement model according to the [mu] M-PMU amount measurement; if the step is the first mu M-PMU sampling moment in one period, turning to the step 6; 4, performing static robust state estimation; turning to the step 5; 5, updating the internal quantity of dynamic filtering, and associating data of two periods; turning to step 6; 6, performing dynamic filtering, and turning to the step 7; 7, performing state prediction according to the dynamic model; according to the method, the data in the previous and later periods are associated, and the estimation result of the static robust state estimation is effectively utilized.

Description

Technical Field

[0001] The present invention relates to a method for estimating the dynamic state of a smart distribution network, and in particular to a method for estimating the dynamic state of a smart distribution network with two-cycle data association. Background Art

[0002] As the link directly facing power users in the "generation, transmission, distribution and use" of the power system, the safety and economy of the distribution network are very important. The safe and economical operation of the distribution network depends on real-time, reliable and complete distribution network operation status data. Due to errors or failures in measurement devices, communication system errors, potential data attacks and other reasons, the data obtained by the distribution network control center often contains errors or bad data. Therefore, distribution network state estimation is required to pre-process the data to improve the accuracy and quantity of the data. Among them, dynamic state estimation relies on iterative estimation based on the measurement data of the previous and next two time sections, which can provide both the estimated value and the predicted value of the system state.

[0003] Micro Phasor Measurement Unit (μM-PMU) is an important data source for monitoring, measurement and control of distribution management systems in distribution networks. Due to construction planning and cost control, the SCADA system measurement and μM-PMU measurement will coexist for a long time in the current smart distribution network. In the hybrid measurement of distribution networks, μM-PMU is mainly used to measure the voltage phasor and current phasor of the deployed nodes. For nodes without μM-PMU deployment, SCADA can collect branch power, node injection power and node voltage amplitude. In order to increase the redundancy of state estimation, a unified measurement model is established for the results of SCADA and μM-PMU measurements. For this purpose, measurement transformation is required to uniformly convert the measurements of the two into a rectangular coordinate system.

[0004] As μM-PMU is gradually deployed in the distribution network, μM-PMU and SCADA form a long and short two-cycle sampling mode. The comprehensive use of μM-PMU and SCADA data to estimate the dynamic state of the power system is a demand and hot spot in the industry. The existing two-cycle or multi-cycle state estimation method cannot use the redundant information provided by SCADA to improve the robustness, and there is no data association between the two long cycles before and after, resulting in insufficient data utilization efficiency.

[0005] Smart distribution networks have various types of user loads, distributed photovoltaics, and wind turbines. Their power changes will bring many challenges to the operation of the distribution network: after the access of new energy output, the distribution network becomes a multi-power system, making it difficult to converge the power flow calculation based on forward and backward generation; the volatility and intermittent nature of wind power and photovoltaic power generation make the operating point of the distribution system constantly changing, which brings challenges to state estimation. For smart distribution networks, static state estimation can no longer meet the demand for mastering the dynamic trend of the system. Dynamic state estimation provides the predicted value of injected power by introducing ultra-short-term power prediction, thereby further mastering the changing trend of node voltage amplitude and voltage phase angle in the distribution network.

[0006] The methods for solving linear algebraic equations can be roughly divided into two categories: direct method and iterative method. The advantage of the direct method is that it can obtain an accurate solution without taking into account rounding errors. However, when the condition number of the coefficient matrix is ​​large, that is, when the coefficient matrix is ​​ill-conditioned, the influence of rounding errors often causes the obtained solution to be far from the accurate solution. In addition, when the scale of the equation system is large, the direct method generally requires a large amount of memory and a long time because there may be many fill-in elements after the matrix decomposition. The iterative method is one of the hot topics in the mathematical research on solving equation systems. Its advantage is that the storage requirement is smaller than that of the direct method, and it has a better speed in solving large linear equation systems. Therefore, the CPU-based solution of the equation system is divided into the following three steps, which is also the common practice of large solvers:

[0007] 1) LU symbolic decomposition of the matrix

[0008] Symbolic decomposition is a simulated LU decomposition of the matrix without performing decomposition calculations on the specific element values. Its purpose is to record the positions of the elements involved in the decomposition and the positions of the filled elements, as well as the positions of the elements eliminated during the decomposition process, so as to establish a storage structure for accessing data for numerical LU decomposition. In this way, in the process of numerical LU decomposition, based on these records, a simple and direct addressing method is used, which greatly reduces the amount of data query and significantly improves the calculation speed. When the matrix structure remains unchanged, symbolic decomposition only needs to be performed once.

[0009] 2) LU numerical decomposition of the matrix

[0010] For the distribution network, the network structure is likely to change, so LU decomposition is generally performed when the structure of the power grid changes, and the number of decompositions is relatively small. In the continuous dynamic state estimation process, the number of LU decompositions is not large, so the LU decomposition speed is not the main time bottleneck for solving the system of equations.

[0011] 3) Solve the equation by back-substitution

[0012] This is a conventional factor table solution method. The solution is calculated based on the LU matrix and the known vector on the right side, and the solution is calculated through forward and back substitution processes.

[0013] Therefore, according to the inventor's research and actual computing experience, static robust state estimation is introduced into the two-cycle dynamic state estimation, and a two-cycle data fusion method is established, which can effectively use SCADA data for robust estimation and use the estimation results to associate the data of the two cycles before and after. Summary of the invention

[0014] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and provide a dynamic state estimation method of a smart distribution network with two-cycle data association.

[0015] The purpose of the present invention can be achieved by the following technical solutions:

[0016] A two-cycle data-associated robust dynamic state estimation method for a smart distribution network, the method comprising the following steps:

[0017] (1) Perform distribution network measurement transformation, establish a unified measurement model, and go to step (2);

[0018] (2) Perform ultra-short-term power prediction and accelerate the matrix calculation in parallel, and then go to step (3);

[0019] (3) Determine whether it is a SCADA sampling moment. If so, use μM-PMU and SCADA measurements to form a full measurement model and go to step (4); otherwise, form a linear measurement model based on μM-PMU measurements; if this step is the first μM-PMU sampling moment in a cycle, go to step (6);

[0020] (4) Perform static robustness state estimation; go to step (5)

[0021] (5) Update the internal quantity of dynamic filtering and associate the two-cycle data; go to step (6)

[0022] (6) Perform dynamic filtering and go to step (7);

[0023] (7) State prediction based on dynamic models;

[0024] Furthermore, the distribution network measurement transformation in step (1) mainly includes the following steps:

[0025] In the rectangular coordinate system, after the μM-PMU measurement is transformed, we can get:

[0026]

[0027]

[0028]

[0029] Where i and j are node numbers, is the phase number, is the node i Phase voltage amplitude, is the node i Phase voltage phase angle, For branch ij Phase current amplitude, and In the rectangular coordinate system The real and imaginary parts of and for Real and imaginary measurements, is the phase angle difference between the first and the last end of branch ij, is the node i Phase injection current amplitude, is the corresponding injection current phase angle, and is the node i Real and imaginary parts of phase injected current.

[0030] In the rectangular coordinate system, after the measurement transformation of the SCADA branch power measurement and the injected power measurement, we can get:

[0031]

[0032]

[0033] In the formula, and SCADA measures the branch ij Phase active and reactive power. and is the SCADA measurement node i Phase active injection power and reactive injection power.

[0034] After the measurement transformation, the mixed measurement vector is obtained:

[0035]

[0036] Where z is composed of the μM-PMU measurement and SCADA measurement. Specifically refers to the voltage amplitude measured by SCADA, and It is the result of voltage phase measured by μM-PMU after measurement transformation.

[0037] Consider the following power system measurement model:

[0038] z=h(x)+v

[0039] Where z is the mixed quantity measurement vector, x is the state quantity vector, v is the measurement error, and h(x) is the system measurement function.

[0040] In order to utilize the voltage amplitude measurement of SCADA, the real part, imaginary part and voltage amplitude of the phase node voltage are selected as state variables:

[0041]

[0042] In the formula is an additional state variable, together with the other two state variables and There are quadratic form constraints:

[0043]

[0044] According to the above measurement model, detailed measurement equations are established for phase branch current, node injection current and node voltage measurement respectively.

[0045] 2.1SCADA Measurement

[0046] The voltage amplitude measurement equation is:

[0047]

[0048] In the formula and are the true value and measured value of the node phase voltage amplitude, and Δ represents the measurement error.

[0049] The measurement equation of the phase branch ij is:

[0050]

[0051] In the formula and is the phase conductance and susceptance of branch ij to ground, and are the phase conductance and susceptance of the corresponding branch.

[0052] The node injection current measurement equation is:

[0053]

[0054] 2.2μM-PMU measurement

[0055] The voltage measurement equation is:

[0056]

[0057] The measurement equation of the phase branch ij is:

[0058]

[0059] The node injection current measurement equation is:

[0060]

[0061] After measurement transformation, there is a linear relationship between the measured quantity and the state quantity, and the hybrid measurement model of the distribution network is:

[0062]

[0063] In short:

[0064]

[0065] Where g(x) is a quadratic function constraint on x.

[0066] After the above measurement transformation, a quadratic constrained linear measurement model about the state variables is obtained, wherein the measurement model is linear, but a quadratic constraint needs to be added to the state estimation model.

[0067] Furthermore, the parallelization acceleration of the ultra-short-term power prediction matrix calculation in step (2) mainly includes the following steps:

[0068] The ultra-short-term power prediction matrix can open up multiple threads for centralized computing. If distributed computing is used, the computing tasks of this node in the model are assigned from the upper computing node to the lower sub-nodes for computing.

[0069] If there is only one computing node in the distribution network, the node needs to obtain all the information and perform ultra-short-term power forecasting through the network data and status data in the shared memory. When there are computing nodes in multiple substations in the distribution network, distributed ultra-short-term power forecasting can also be performed in multiple substations to reduce the computing pressure of the core computing nodes. In this process, two-way communication between the upper and lower sites is essential.

[0070] In the distribution network power flow calculation, the key iterative steps are:

[0071] v (k+1) =(Y+Y L ) -1 (i PQ (v)+i I (v) -Y NS v S )

[0072] Rewrite the above formula as:

[0073] v(k+1) =A -1 f(v)

[0074] in

[0075]

[0076] The iterative method for solving nonlinear equations transforms the ultra-short-term power forecasting problem into an iterative solution of a series of linear equations. Its characteristics are that the matrix A keeps the matrix structure and elements unchanged, and only the known vector f(v) on the right side changes.

[0077] In order to achieve efficient matrix solution, the CUDA numerical computing library cuSparse, cuBLAS and cuSolver are used in the programming implementation. These three libraries can realize matrix-vector calculation based on GPU+CPU. cuSolver has a three-layer structure, in which the first layer is LU decomposition, QR decomposition, SVD decomposition and LDL decomposition of dense matrices. T Decomposition, as well as multiplication operations of dense matrices and vectors. This layer takes advantage of the massive parallel computing capabilities of the GPU to accelerate the calculation of dense matrices.

[0078] The second layer is a set of sparse matrix decomposition library functions, among which the representative one is sparse QR decomposition. For distribution networks, due to the weak ring network or radial nature, the linear equations associated with their admittance matrix are very suitable for sparse parallel solving. In this case, cuSolver can provide efficient acceleration.

[0079] The third layer is the library functions related to matrix refactoring, which provide efficient solutions when only non-zero elements change while the sparsity characteristics remain unchanged.

[0080] For linear equations

[0081] Ax=b

[0082] Where A is an n-dimensional sparse matrix, and A is decomposed into LU with sorting:

[0083] PA=LU

[0084] The LU decomposition step can be parallelized or executed sequentially. cuSparse provides the csrlsvlu function for LU decomposition of the matrix and the csrlsvqr function for QR decomposition of the matrix. The QR decomposition is:

[0085] A=QR

[0086] Decompose the matrix A into an orthogonal matrix Q and an upper triangular matrix R.

[0087] In practice, the cusolverSpDcsrlsvqr function is used to perform matrix calculations on the GPU. If the matrix parallelization effect is not good, the cusolverSpDcsrlsvluHost function is used to perform multi-threaded matrix calculations on the CPU.

[0088] Furthermore, the construction of the static model and dynamic model of the system in step (3) mainly includes the following steps:

[0089] 4.1 Construct static and dynamic equations, the compact form is:

[0090]

[0091] x is the state variable of the system dynamic process, including the active and reactive state variables P in the dynamic load model. 0,i , Q 0,i Or the virtual angular velocity ω in the Virtual Synchronous Machine (VSM) model i , that is, x=[P 0,i ,Q 0,i ,ω i ]; y is the algebraic variable of the system, including the voltage amplitude or current amplitude of the system node, that is, y = [V m,i ,V a,i ,I m,i ,I a,i ,P i ,Q i ]Where: V m,i ,V a,i is the magnitude and phase angle of the voltage phasor at node i, I m,i ,I a,i is the amplitude and phase angle of the current phasor at node i, P i ,Q i is the active and reactive power injected into node i; u is the disturbance in the system, including random fluctuations in the output of renewable energy generation or load power, i.e., ΔP i ,ΔQ i , is the static equation of the system, is the dynamic equation of the system.

[0092] 4.2 Constructing the state transfer equation

[0093] Discretize the formula in 2.1 to obtain the state transfer equation:

[0094] x k+1 =f(x k )+w k

[0095] where wk is the perturbation vector.

[0096] Further, step (4) includes:

[0097] 5.1 Forming a full measurement model:

[0098] The full measurement model refers to the measurement model formed by μM-PMU and SCADA after step (1). Its compact form is:

[0099] z=h(x)+v

[0100] st

[0101] In the formula is the set of μM-PMU and SCADA measurement variables, h(x) is the formed measurement model, which represents the function mapping from the system state variables to the measurement variables, and v represents the measurement error. Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between .

[0102] 5.2 Forming a linear measurement model based on μM-PMU

[0103] Since the measurement model formed by μM-PMU is a linear equation, it is called a linear measurement model, and its compact form is:

[0104] z=Hx+v

[0105] st

[0106] Where z can be expressed as is the set of μM-PMU measurement variables. H is a real number matrix, which is jointly determined by the μM-PMU measurement variables and the system state quantity according to the physical quantity relationship. Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between them.

[0107] Furthermore, the static robust state estimation described in step (4) includes the following steps:

[0108] 6.1 Establish a static robust state estimation model. The compact form of the model is:

[0109]

[0110] Where J(x) is the optimization objective function, x is the system state, r i is the i-th element in the residual vector, w iis the corresponding weight. The corresponding constraint is the rectangular coordinate system variable and polar coordinate variables The quadratic constraint relationship between and the definition constraint of the residual vector.

[0111] L δ For a given robust loss function, the Pseudo-Huber loss function is used:

[0112]

[0113] Where: r is the residual vector, δ is an adjustable parameter, usually 1.5.

[0114] 6.2 Solving the Given Robust Estimation Problem

[0115] For a given problem, use the interior point method to solve and obtain the optimal solution The optimal estimated value of each state quantity of the corresponding system.

[0116] Furthermore, the dynamic filtering internal quantity update described in step (5) includes the following steps:

[0117] 7.1 Dynamic Filter State Update

[0118]

[0119] Where x k|k is the correction value of the state variable in the dynamic filtering at that moment, is the optimal solution obtained in step (3).

[0120] 7.2 Fading Memory Residual Covariance Matrix Update

[0121]

[0122] Where V k is the fading memory residual covariance matrix in the dynamic filter at that moment, z k is the measurement variable at the current moment, and h(·) is the measurement function.

[0123] 7.3 Prediction Covariance Matrix Update

[0124]

[0125] In the formula is the prediction covariance matrix in the dynamic filtering at the moment, x k+1|k is the predicted value of the system state quantity, χ i,k+1|k is the Sigma point set after unscented transformation (UT), is the corresponding covariance weight obtained by UT transformation. k is the covariance of the noise matrix.

[0126] Furthermore, the dynamic filtering described in step (6) mainly includes the following steps:

[0127] 8.1 Computing the Autovariance and Cross-Variance Matrices

[0128] According to x k+1|k and P k+1|k Get the Sigma point set {ξ i,k}, and the Sigma point {y i,k+1}, calculate the autovariance matrix and the cross-variance matrix.

[0129]

[0130] In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix, R k is the variance of the measurement error.

[0131] 8.2 Calculating the residual covariance V of the fading memory k .

[0132]

[0133] Where ρ is the forgetting factor. The larger ρ is, the better k The stronger the memory of historical residual covariance, r is the residual.

[0134] 8.3 Calculation of the fading factor λ k

[0135]

[0136] M k =S k+1 -V k+1 +N k

[0137]

[0138] λ k =max(1,τ k )

[0139] Where β is the weakening factor and tr is the trace of the matrix.

[0140] 8.4 Modify the prediction covariance matrix.

[0141]

[0142] According to x k+1|k and Calculate the Sigma point set again

[0143]

[0144] In the formula { i} is the required Sigma point set, generated using a symmetric sampling strategy, L is the total number of sampling points, is the corresponding mean weight of each Sigma sampling point, is the variance weight of the corresponding Sigma sampling points. i}, the transformed mean and variance are:

[0145]

[0146] 8.5 Update:

[0147]

[0148] According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k} can be obtained:

[0149]

[0150] Furthermore, the state prediction according to the dynamic model in step (7) mainly includes the following steps:

[0151] 9.1 Prediction based on (UnscentedKalmanFilter, UKF) prediction step

[0152] Construct a given Sigma point set {χ i,k}, based on the UKF prediction step, predict the state variables at the next moment:

[0153]

[0154] Where x k+1|k is the state variable at the next predicted μM-PMU sampling time. k+1|k is the predicted error covariance matrix, Q k is the model error covariance.

[0155] 9.2 Update the predicted values ​​of the self-variance matrix and the cross-variance matrix

[0156] According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k} can be obtained:

[0157]

[0158]

[0159] In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix.

[0160] The beneficial effects of the present invention are as follows: the method of the present invention adopts robust estimation, and the influence of bad data is significantly reduced; at the same time, the present invention associates the data in the two previous and next cycles, and effectively utilizes the estimation results of static robust state estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0161] Figure 1 A diagram of the architecture of a two-cycle data-associated intelligent distribution network robust dynamic state estimation method provided by the present invention;

[0162] Figure 2 A flow chart of a two-cycle data-associated intelligent distribution network robust dynamic state estimation method provided by the present invention;

[0163] Figure 3 A parallel computing architecture diagram of an ultra-short-term prediction module of a two-cycle data-associated intelligent distribution network robust dynamic state estimation method provided by the present invention.

[0164] Figure 4 The present invention provides a wiring diagram of an IEEE33 node system in a two-cycle data-associated intelligent distribution network robust dynamic state estimation method.

[0165] Figure 5 This is a diagram of estimation results when there is no bad data in a two-cycle data associated intelligent distribution network robust dynamic state estimation method provided by the present invention.

[0166] Figure 6 A diagram of estimation results when bad data exists in a two-cycle data-associated intelligent distribution network robust dynamic state estimation method provided by the present invention. DETAILED DESCRIPTION

[0167] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0168] A two-cycle data-associated robust dynamic state estimation method for a smart distribution network, the method comprising the following steps:

[0169] (1) Perform distribution network measurement transformation, establish a unified measurement model, and go to step (2);

[0170] (2) Perform ultra-short-term power prediction and accelerate the matrix calculation in parallel, and then go to step (3);

[0171] (3) Determine whether it is a SCADA sampling moment. If so, use μM-PMU and SCADA measurements to form a full measurement model and go to step (4); otherwise, form a linear measurement model based on μM-PMU measurements; if this step is the first μM-PMU sampling moment in a cycle, go to step (6);

[0172] (4) Perform static robustness state estimation; go to step (5)

[0173] (5) Update the internal quantity of dynamic filtering and associate the two-cycle data; go to step (6)

[0174] (6) Perform dynamic filtering and go to step (7);

[0175] (7) State prediction based on dynamic models;

[0176] Furthermore, the distribution network measurement transformation in step (1) mainly includes the following steps:

[0177] In the rectangular coordinate system, after the μM-PMU measurement is transformed, we can get:

[0178]

[0179]

[0180]

[0181] Where i and j are node numbers, φ∈{a,b,c} is the phase number, is the φ phase voltage amplitude at node i, is the φ phase voltage phase angle of node i, is the φ-phase current amplitude of branch ij, and In the rectangular coordinate system The real and imaginary parts of and for Real and imaginary measurements, is the phase angle difference between the first and the last end of branch ij, is the φ-phase injected current amplitude of node i, is the corresponding injection current phase angle, and The real and imaginary parts of the current are injected into the φ phase of node i.

[0182] In the rectangular coordinate system, after the measurement transformation of the SCADA branch power measurement and the injected power measurement, we can get:

[0183]

[0184]

[0185] In the formula, and It is used to measure the active and reactive power of the φ phase of the lower branch ij for SCADA. and is the φ-phase active injection power and reactive injection power of node i under SCADA measurement.

[0186] After the measurement transformation, the mixed measurement vector is obtained:

[0187]

[0188] Where z is composed of the μM-PMU measurement and SCADA measurement. Specifically refers to the voltage amplitude measured by SCADA, and It is the result of voltage phase measured by μM-PMU after measurement transformation.

[0189] Consider the following power system measurement model:

[0190] z=h(x)+v

[0191] Where z is the mixed quantity measurement vector, x is the state quantity vector, v is the measurement error, and h(x) is the system measurement function.

[0192] In order to utilize the voltage amplitude measurement of SCADA, the real part, imaginary part and voltage amplitude of the phase node voltage are selected as state variables:

[0193]

[0194] In the formula is an additional state variable, together with the other two state variables and There are quadratic form constraints:

[0195]

[0196] According to the above measurement model, detailed measurement equations are established for phase branch current, node injection current and node voltage measurement respectively.

[0197] 2.1SCADA Measurement

[0198] The voltage amplitude measurement equation is:

[0199]

[0200] In the formula and are the true value and measured value of the node phase voltage amplitude, and Δ represents the measurement error.

[0201] The measurement equation of the phase branch ij is:

[0202]

[0203] In the formula and is the phase conductance and susceptance of branch ij to ground, and are the phase conductance and susceptance of the corresponding branch.

[0204] The node injection current measurement equation is:

[0205]

[0206] 2.2μM-PMU measurement

[0207] The voltage measurement equation is:

[0208]

[0209] The measurement equation of the phase branch ij is:

[0210]

[0211] The node injection current measurement equation is:

[0212]

[0213] After measurement transformation, there is a linear relationship between the measured quantity and the state quantity, and the hybrid measurement model of the distribution network is:

[0214]

[0215] In short:

[0216] z=Hx+v

[0217] stg(x)=0(0-3)

[0218] Where g(x) is a quadratic function constraint on x.

[0219] After the above measurement transformation, a quadratic constrained linear measurement model about the state variables is obtained, wherein the measurement model is linear, but a quadratic constraint needs to be added to the state estimation model.

[0220] Furthermore, the parallelization acceleration of the ultra-short-term power prediction matrix calculation in step (2) mainly includes the following steps:

[0221] The ultra-short-term power prediction matrix can open up multiple threads for centralized computing. If distributed computing is used, the computing tasks of this node in the model are assigned from the upper computing node to the lower sub-nodes for computing.

[0222] If there is only one computing node in the distribution network, the node needs to obtain all the information and perform ultra-short-term power forecasting through the network data and status data in the shared memory. When there are computing nodes in multiple substations in the distribution network, distributed ultra-short-term power forecasting can also be performed in multiple substations to reduce the computing pressure of the core computing nodes. In this process, two-way communication between the upper and lower sites is essential.

[0223] In the distribution network power flow calculation, the key iterative steps are:

[0224] v (k+1) =(Y+Y L ) -1 (i PQ (v)+i I (v) -Y NS v S )

[0225] Rewrite the above formula as:

[0226] v (k+1) =A -1 f(v)

[0227] in

[0228]

[0229] The iterative method for solving nonlinear equations transforms the ultra-short-term power forecasting problem into an iterative solution of a series of linear equations. Its characteristics are that the matrix A keeps the matrix structure and elements unchanged, and only the known vector f(v) on the right side changes.

[0230] In order to achieve efficient matrix solution, the CUDA numerical computing library cuSparse, cuBLAS and cuSolver are used in the programming implementation. These three libraries can realize matrix-vector calculation based on GPU+CPU. cuSolver has a three-layer structure, in which the first layer is LU decomposition, QR decomposition, SVD decomposition and LDL decomposition of dense matrices. T Decomposition, as well as multiplication operations of dense matrices and vectors. This layer takes advantage of the massive parallel computing capabilities of the GPU to accelerate the calculation of dense matrices.

[0231] The second layer is a set of sparse matrix decomposition library functions, among which the representative one is sparse QR decomposition. For distribution networks, due to the weak ring network or radial nature, the linear equations associated with their admittance matrix are very suitable for sparse parallel solving. In this case, cuSolver can provide efficient acceleration.

[0232] The third layer is the library functions related to matrix refactoring, which provide efficient solutions when only non-zero elements change while the sparsity characteristics remain unchanged.

[0233] For linear equations

[0234] Ax=b

[0235] Where A is an n-dimensional sparse matrix, and A is decomposed into LU with sorting:

[0236] PA=LU

[0237] The LU decomposition step can be parallelized or executed sequentially. cuSparse provides the csrlsvlu function for LU decomposition of the matrix and the csrlsvqr function for QR decomposition of the matrix. The QR decomposition is:

[0238] A=QR

[0239] Decompose the matrix A into an orthogonal matrix Q and an upper triangular matrix R.

[0240] In practice, the cusolverSpDcsrlsvqr function is used to perform matrix calculations on the GPU. If the matrix parallelization effect is not good, the cusolverSpDcsrlsvluHost function is used to perform multi-threaded matrix calculations on the CPU.

[0241] Furthermore, the construction of the static model and dynamic model of the system in step (3) mainly includes the following steps:

[0242] 4.1 Construct static and dynamic equations, the compact form is:

[0243]

[0244] x is the state variable of the system dynamic process, including the active and reactive state variables P in the dynamic load model. 0,i , Q 0,i Or the virtual angular velocity ω in the Virtual Synchronous Machine (VSM) model i , that is, x=[P 0,i ,Q 0,i ,ω i]; y is the algebraic variable of the system, including the voltage amplitude or current amplitude of the system node, that is, y = [V m,i ,V a,i ,I m,i ,I a,i ,P i ,Q i ]Where: V m,i ,V a,i is the magnitude and phase angle of the voltage phasor at node i, I m,i ,I a,i is the amplitude and phase angle of the current phasor at node i, P i ,Q i is the active and reactive power injected into node i; u is the disturbance in the system, including random fluctuations in the output of renewable energy generation or load power, i.e., ΔP i ,ΔQ i , is the static equation of the system, is the dynamic equation of the system.

[0245] 4.2 Constructing the state transfer equation

[0246] Discretize the formula in 2.1 to obtain the state transfer equation:

[0247] x k+1 =f(x k )+w k

[0248] where w k is the perturbation vector.

[0249] Further, step (4) includes:

[0250] 5.1 Forming a full measurement model:

[0251] The full measurement model refers to the measurement model formed by μM-PMU and SCADA after step (1). Its compact form is:

[0252] z=h(x)+v

[0253] st

[0254] In the formula is the set of μM-PMU and SCADA measurement variables, h(x) is the formed measurement model, which represents the function mapping from the system state variables to the measurement variables, and v represents the measurement error. Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between them.

[0255] 5.2 Forming a linear measurement model based on μM-PMU

[0256] Since the measurement model formed by μM-PMU is a linear equation, it is called a linear measurement model, and its compact form is:

[0257] z=Hx+v

[0258] st

[0259] Where z can be expressed as is the set of μM-PMU measurement variables. H is a real number matrix, which is jointly determined by the μM-PMU measurement variables and the system state quantity according to the physical quantity relationship. Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between them.

[0260] Furthermore, the static robust state estimation described in step (4) includes the following steps:

[0261] 6.1 Establish a static robust state estimation model. The compact form of the model is:

[0262]

[0263] Where J(x) is the optimization objective function, x is the system state, r i is the i-th element in the residual vector, w i is the corresponding weight. The corresponding constraint is the rectangular coordinate system variable and polar coordinate variables The quadratic constraint relationship between and the definition constraint of the residual vector.

[0264] L δ For a given robust loss function, the Pseudo-Huber loss function is used:

[0265]

[0266] Where: r is the residual vector, δ is an adjustable parameter, usually 1.5.

[0267] 6.2 Solving the Given Robust Estimation Problem

[0268] For a given problem, use the interior point method to solve and obtain the optimal solution The optimal estimated value of each state quantity of the corresponding system.

[0269] Furthermore, the dynamic filtering internal quantity update described in step (5) includes the following steps:

[0270] 7.1 Dynamic Filter State Update

[0271]

[0272] Where x k|k is the correction value of the state variable in the dynamic filtering at that moment, is the optimal solution obtained in step (3).

[0273] 7.2 Fading Memory Residual Covariance Matrix Update

[0274]

[0275] Where V k is the fading memory residual covariance matrix in the dynamic filter at that moment, z k is the measurement variable at the current moment, and h(·) is the measurement function.

[0276] 7.3 Prediction Covariance Matrix Update

[0277]

[0278] In the formula is the prediction covariance matrix in the dynamic filtering at the moment, x k+1|k is the predicted value of the system state quantity, χ i,k+1|k is the Sigma point set after unscented transformation (UT), is the corresponding covariance weight obtained by UT transformation. k is the covariance of the noise matrix.

[0279] Furthermore, the dynamic filtering described in step (6) mainly includes the following steps:

[0280] 8.1 Computing the Autovariance and Cross-Variance Matrices

[0281] According to x k+1|k and P k+1|k Get the Sigma point set {ξ i,k}, and the Sigma point {y i,k+1}, calculate the autovariance matrix and the cross-variance matrix.

[0282]

[0283] In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix, and Rk is the variance of the measurement error.

[0284] 8.2 Calculate the residual covariance Vk of the fading memory.

[0285]

[0286] Where ρ is the forgetting factor. The larger ρ is, the stronger Vk’s memory of the historical residual covariance is, and r is the residual.

[0287] 8.3 Calculation of the fading factor λk

[0288]

[0289] M k =S k+1 -V k+1 +N k

[0290]

[0291] λ k =max(1,τ k )

[0292] Where β is the weakening factor and tr is the trace of the matrix.

[0293] 8.4 Modify the prediction covariance matrix.

[0294]

[0295] According to x k+1|k and Calculate the Sigma point set again

[0296]

[0297] In the formula { i} is the required Sigma point set, generated using a symmetric sampling strategy, L is the total number of sampling points, is the corresponding mean weight of each Sigma sampling point, is the variance weight of the corresponding Sigma sampling points. i}, the transformed mean and variance are:

[0298]

[0299] 8.5 Update:

[0300]

[0301] According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k} can be obtained:

[0302]

[0303] Furthermore, the state prediction according to the dynamic model in step (7) mainly includes the following steps:

[0304] 9.1 Prediction based on the (Unscented Kalman Filter, UKF) prediction step

[0305] Construct a given Sigma point set {χ i,k}, based on the UKF prediction step, predict the state variables at the next moment:

[0306]

[0307] Where x k+1|k is the state variable at the next predicted μM-PMU sampling time. k+1|k is the predicted error covariance matrix, and Qk is the model error covariance.

[0308] 9.2 Update the predicted values ​​of the self-variance matrix and the cross-variance matrix

[0309] According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k} can be obtained:

[0310]

[0311]

[0312] In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix.

[0313] Example

[0314] This embodiment is described by taking the IEEE 33-node system as an example. The system wiring diagram is as follows: Figure 4 As shown in , when the measurement is normal and there is no bad data, the state estimation is performed using this method and compared with the state estimation achieved by the conventional weighted least squares method, as shown in Figure 5 As shown, WLS represents the conventional weighted least squares state estimation, and Robust represents the estimation result when there is no bad data.

[0315] When there is no bad data, the calculation error of this method is similar to that of the least squares solution.

[0316] 2% of SCADA measurements are randomly selected and set to 5 times the true value to simulate the bad data scenario. The estimated results are as follows: Figure 6 As shown, WLS represents the conventional weighted least squares state estimation, and Robust represents the estimation result when bad data exists in this method.

[0317] When bad data exists in the system, the least square method is significantly affected by gross errors, and large errors appear in both the voltage amplitude and the voltage phase angle. The gross errors have little effect on this method, which effectively suppresses the adverse effects of bad data on state estimation.

[0318] This embodiment is only one case where the method is applied. Under other system conditions involving dynamic state estimation, the method can be used to improve the state estimation accuracy.

Claims

1. A two-cycle data association method for robust dynamic state estimation of smart distribution network, characterized in that: The method comprises the following steps: (1) Perform distribution network measurement transformation, establish a unified measurement model, and go to step (2); (2) Perform ultra-short-term power prediction and accelerate the matrix calculation in parallel, and then go to step (3); (3) Determine whether it is a sampling moment of the Supervisory Control And Data Acquisition (SCADA) system. If so, use the micro-phasor measurement unit (μM-PMU) and SCADA measurement to form a full measurement model, and go to step (4); otherwise, form a linear measurement model based on the μM-PMU measurement; if this step is the first μM-PMU sampling moment in a cycle, go to step (6); (4) Perform static robustness state estimation; go to step (5); (5) Update the internal quantity of dynamic filtering and associate the data of two cycles; Go to step (6); (6) Perform dynamic filtering and go to step (7); (7) Perform state prediction based on the dynamic model.

2. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The distribution network measurement transformation in step (1) comprises the following steps: In the rectangular coordinate system, after the μM-PMU measurement is transformed, we can get: Where i and j are node numbers, is the phase number, is the node i Phase voltage amplitude, is the node i Phase voltage phase angle, For branch ij Phase current amplitude, and In the rectangular coordinate system The real and imaginary parts of and for Real and imaginary measurements, is the phase angle difference between the first and the last end of branch ij, is the node i Phase injection current amplitude, is the corresponding injection current phase angle, and is the node i Real and imaginary parts of phase injection current; In the rectangular coordinate system, after the measurement transformation of the SCADA branch power measurement and the injected power measurement, we can get: In the formula, and SCADA measures the branch ij Phase active and reactive power, and is the SCADA measurement of node i Phase active injection power and reactive injection power; After the measurement transformation, the mixed measurement vector is obtained: Where z is composed of the μM-PMU measurement and SCADA measurement. Specifically refers to the voltage amplitude measured by SCADA, and It is the result of voltage phase measured by μM-PMU after measurement transformation; Consider the following power system measurement model: z=h(x)+v Where z is the mixed quantity measurement vector, x is the state quantity vector, v is the measurement error, and h(x) is the system measurement function; In order to utilize the voltage amplitude measurement of SCADA, the real part, imaginary part and voltage amplitude of the phase node voltage are selected as state variables: In the formula is an additional state variable, and the other two state variables and There are quadratic form constraints: According to the above measurement model, detailed measurement equations are established for phase branch current, node injection current and node voltage measurement respectively; SCADA Measurements: The voltage amplitude measurement equation is: In the formula and are the true value and measured value of the node phase voltage amplitude, Δ represents the measurement error; The measurement equation of the phase branch ij is: In the formula and is the phase conductance and susceptance of branch ij to ground, and is the phase conductance and susceptance of the corresponding branch; The node injection current measurement equation is: μM-PMU measurement: The voltage measurement equation is: The measurement equation of the phase branch ij is: The node injection current measurement equation is: After measurement transformation, there is a linear relationship between the measured quantity and the state quantity, and the hybrid measurement model of the distribution network is: i∈{1,…,N}, In short: Where g(x) is a quadratic function constraint about x; After the above measurement transformation, a quadratic constrained linear measurement model about the state variables is obtained, wherein the measurement model is linear, but a quadratic constraint needs to be added to the state estimation model.

3. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The parallel acceleration of the ultra-short-term power prediction matrix calculation in step (2) comprises the following steps: The ultra-short-term power prediction matrix can open up multiple threads for centralized computing. If distributed computing is used, the computing task of this node in the ultra-short-term power prediction matrix is ​​assigned by the upper computing node to the lower sub-node for computing; If there is only one computing node in the distribution network, the node needs to obtain all the information and perform ultra-short-term power forecasting through the network data and status data in the shared memory. When there are computing nodes in multiple substations in the distribution network, distributed ultra-short-term power forecasting is performed in multiple substations to reduce the computing pressure of the core computing nodes. In this process, two-way communication is carried out between the upper and lower sites. In the distribution network power flow calculation, the key iterative steps are: v (k+1) =(Y+Y L ) -1 (i PQ (v)+i I (v)-Y NS v S ) Rewrite the above formula as: v (k+1) =A -1 f(v) in The iterative method for solving nonlinear equations transforms the ultra-short-term power forecasting problem into an iterative solution of a series of linear equations. Its characteristics are that the matrix A keeps the matrix structure and elements unchanged, and only the known vector f(v) on the right side changes; In order to achieve efficient matrix solution, the CUDA numerical computing library cuSparse, cuBLAS and cuSolver are used in the programming implementation. These three libraries can realize matrix-vector calculation based on GPU+CPU; cuSolver has a three-layer structure, in which the first layer is LU decomposition, QR decomposition, SVD decomposition and LDL decomposition of dense matrices. T Decomposition and multiplication of dense matrices and vectors. This layer uses the large-scale parallel computing capabilities of the GPU to accelerate the calculation of dense matrices. The second layer is a set of sparse matrix decomposition library functions, among which the most representative one is sparse QR decomposition. For distribution networks, due to the weak ring network or radial nature, the linear equations associated with its admittance matrix are very suitable for sparse parallel solving. In this case, cuSolver provides efficient acceleration; The third layer is the library functions related to matrix refactoring, which provide efficient solutions when only non-zero elements change while the sparsity characteristics remain unchanged; For linear equations Ax=b Where A is an n-dimensional sparse matrix, and A is decomposed into LU with sorting: PA=LU The LU decomposition step is performed in parallel or sequentially. cuSparse provides the csrlsvlu function for LU decomposition of the matrix and the csrlsvqr function for QR decomposition of the matrix. The QR decomposition is: A=QR Decompose the matrix A into an orthogonal matrix Q and an upper triangular matrix R; The cusolverSpDcsrlsvqr function is used to perform matrix calculations on the GPU. If the matrix parallelization effect is not good, the cusolverSpDcsrlsvluHost function is used to perform multi-threaded matrix calculations on the CPU.

4. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: In step (3), building the static model and dynamic model of the system includes the following steps: 4.1 Construct static and dynamic equations, the compact form is: x is the state variable of the system dynamic process, including the active and reactive state variables P in the dynamic load model. 0,i , Q 0,i Or the virtual angular velocity ω in the Virtual Synchronous Machine (VSM) model i , that is, x=[P 0,i ,Q 0,i ,ω i ]; y is the algebraic variable of the system, including the voltage amplitude or current amplitude of the system node, that is, y = [V m,i ,V a,i ,I m,i ,I a,i ,P i ,Q i ]Where: V m,i ,V a,i is the amplitude and phase angle of the voltage phasor at node i, I m,i ,I a,i is the amplitude and phase angle of the current phasor at node i, P i ,Q i is the active and reactive power injected into node i; u is the disturbance in the system, including random fluctuations in the output of renewable energy generation or load power, i.e., ΔP i ,ΔQ i , is the static equation of the system, is the dynamic equation of the system; 4.2 Constructing the state transfer equation Discretize the formula in 4.1 to obtain the state transfer equation: x k+1 =f(x k )+w k where w k is the perturbation vector.

5. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: Step (4) comprises: 5.1 Forming a full measurement model: The full measurement model refers to the measurement model formed by μM-PMU and SCADA after step (1). Its compact form is: z=h(x)+v In the formula is the set of μM-PMU and SCADA measurement variables, h(x) is the measurement model formed, representing the function mapping from the system state variables to the measurement variables, v represents the measurement error, Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between them; 5.2 Forming a linear measurement model based on μM-PMU Since the measurement model formed by μM-PMU is a linear equation, it is called a linear measurement model. Its compact form is: z=Hx+v Where z is represented by is the set of μM-PMU measurement variables, H is a real number matrix, which is jointly determined by the μM-PMU measurement variables and the system state quantity according to the physical quantity relationship. Represents the rectangular coordinate system variable in the measurement and polar coordinate variables The quadratic constraint relationship between .

6. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The static robust state estimation described in step (4) includes the following steps: 6.1 Establish a static robust state estimation model. The compact form of the model is: Where J(x) is the optimization objective function, x is the system state, r i is the i-th element in the residual vector, w i is the corresponding weight, and the corresponding constraint is the rectangular coordinate system variable and polar coordinate variables The quadratic constraint relationship between and the definition constraint of the residual vector; L δ For a given robust loss function, the Pseudo-Huber loss function is used: Where: r is the residual vector, δ is an adjustable parameter, which is 1.5; 6.2 Solving the Given Robust Estimation Problem For a given robust estimation problem, the interior point method is used to solve it and obtain the optimal solution The optimal estimated value of each state quantity of the corresponding system.

7. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The dynamic filtering internal quantity update described in step (5) includes the following steps: 7.1 Dynamic Filter State Update Where x k|k is the correction value of the state variable in the dynamic filtering at that moment, is the optimal solution obtained in step (3); 7.2 Fading Memory Residual Covariance Matrix Update Where V k is the fading memory residual covariance matrix in the dynamic filter at that moment, z k is the measurement variable at the current moment, h(·) is the measurement function; 7.3 Prediction Covariance Matrix Update In the formula is the prediction covariance matrix in the dynamic filtering at the moment, x k+1|k is the predicted value of the system state quantity, χ i,k+1|k is the Sigma point set after unscented transformation (UT), is the corresponding covariance weight obtained by UT transformation, Q k is the covariance of the noise matrix.

8. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The dynamic filtering described in step (6) comprises the following steps: 8.1 Computing the Autovariance and Cross-Variance Matrices According to x k+1|k and P k+1|k Get the Sigma point set {ξ i,k }, and the Sigma point {y i,k+1 }, calculate the autovariance matrix and the cross-variance matrix; In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix, R k is the variance of the measurement error; 8.2 Calculating the residual covariance V of the fading memory k Where ρ is the forgetting factor. The larger ρ is, the better k The stronger the memory of the historical residual covariance, r is the residual; 8.3 Calculation of the fading factor λ k M k =S k+1 -V k+1 +N k l k =max(1,τ k ) Where β is the weakening factor, tr is the trace of the matrix; 8.4 Modified prediction covariance matrix According to x k+1|k and Calculate the Sigma point set again In the formula { i } is the required Sigma point set, generated using a symmetric sampling strategy, L is the total number of sampling points, is the corresponding mean weight of each Sigma sampling point, For the corresponding variance weights of each Sigma sampling point, obtain the Sigma point set {χ i }, the transformed mean and variance are: 8.5 Update: According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k } can be obtained:

9. The method for estimating the robust dynamic state of a smart distribution network based on two-cycle data association according to claim 1, characterized in that: The state prediction according to the dynamic model in step (7) includes the following steps: 9.1 Prediction based on the prediction step of the Unscented Kalman Filter (UKF) Construct a given Sigma point set {χ i,k }, based on the UKF prediction step, predict the state variables at the next moment: Where x k+1|k is the state variable at the next μM-PMU sampling time, P k+1|k is the predicted error covariance matrix, Q k is the model error covariance; 9.2 Update the predicted values ​​of the self-variance matrix and the cross-variance matrix According to x k+1|k With P k+1|k Construct a new Sigma point set {ξ i,k } can be obtained: In the formula is the predicted mean value of the measurement, S k+1 is the variance matrix, C k+1 is the mutual variance matrix.

Images