Dynamic state estimation method of photovoltaic distribution network based on improved volumetric Kalman

Abstract

The present invention discloses a method for estimating the dynamic state of a photovoltaic distribution network based on improved volumetric Kalman, which belongs to the technical field of power systems. First, a general model for estimating the dynamic state of a photovoltaic distribution network is established, including a state transfer equation and a measurement equation; secondly, for the state transfer equation in the model, a photovoltaic side state quantity is added on the basis of the original state quantity of the distribution network, and a one-step prediction of the dynamic state trajectory of the photovoltaic distribution network is performed based on the Holt's two-parameter exponential smoothing equation; next, for the measurement equation in the model, a photovoltaic side quantity measurement is added on the basis of the original PMU quantity measurement of the distribution network, and the measurement equation of the photovoltaic distribution network is derived based on the constraint relationship between the quantity measurement and the state quantity; then an improved volumetric Kalman filtering method is designed; finally, the improved volumetric Kalman filtering method is used to estimate the dynamic state of the photovoltaic distribution network. The method provided by the present invention has good adaptability to the system operation conditions of photovoltaic grid-connected systems.

Description

Technical Field

[0001] The invention belongs to the technical field of power systems, and in particular relates to a method for estimating the dynamic state of a photovoltaic power distribution network based on an improved volumetric Kalman algorithm. Background Art

[0002] State estimation can solve the system operating status online and provide reliable data guarantee for real-time status monitoring, dispatching, control and fault analysis of the distribution system. Dynamic state estimation provides estimated values ​​and forecast values ​​of the state quantities of the two time sections before and after, which can accurately capture the dynamic changes of the system status in real time and is an important part of the safe and stable operation of the active distribution network. However, with the continuous access of distributed photovoltaics, the operation mode and state perception of the active distribution network have added more uncertainties, which puts forward higher requirements for the dynamic state estimation of the active distribution network.

[0003] The cubature Kalman filter (CKF) is based on the spherical-radial rule. It can solve the problem of high-order nonlinear filtering without selecting any parameters and is widely used in dynamic state estimation of distribution networks. However, the covariance matrix of CKF is not positive definite during iteration, which leads to filter interruption or filter divergence. In addition, the fluctuation of distributed photovoltaic output causes the injection power of system nodes to change continuously, which makes the voltage phasor of distribution network nodes change frequently, aggravating the uncertainty of distribution network state quantity and reducing the accuracy of dynamic state estimation results of photovoltaic distribution network. Summary of the invention

[0004] The purpose of the present invention is to provide a method for estimating the dynamic state of a photovoltaic distribution network based on an improved volumetric Kalman filter, so as to solve the problems of filtering interruption or divergence in the prior art and reduced accuracy of dynamic state estimation results of a photovoltaic distribution network due to frequent changes in voltage phasors.

[0005] To achieve the above object, the present invention provides a method for estimating the dynamic state of a photovoltaic power distribution network based on an improved volumetric Kalman filter, comprising the following steps:

[0006] Step 1: First, establish a general model for dynamic state estimation of photovoltaic distribution network, including state transfer equation and measurement equation;

[0007] Step 2: For the state transfer equation in the model, add the photovoltaic side state quantity on the basis of the original state quantity of the distribution network, and perform a one-step forecast of the dynamic state trajectory of the photovoltaic distribution network based on Holt's two-parameter exponential smoothing equation;

[0008] Step 3: For the measurement equation in the model, add the PV side measurement on the basis of the original PMU measurement of the distribution network, and derive the measurement equation of the PV distribution network based on the constraint relationship between the measurement and the state quantity;

[0009] Step 4: design an improved cubature Kalman filter method;

[0010] Step 5: Use the improved cubature Kalman filter method to estimate the dynamic state of the photovoltaic distribution network.

[0011] Preferably, the specific expression of the general model for dynamic state estimation of the photovoltaic distribution network established in step 1 is as follows:

[0012]

[0013] Where: x is the n-dimensional state vector; z is the m-dimensional measurement vector; the subscript k represents the time point; f is the state transfer function, which is the state transfer equation, and expresses the relationship between the state quantity at time k-1 and the state quantity at time k; h is the measurement function, which is the measurement equation, and expresses the relationship between the state quantity at time k and the quantity measurement at time k; q k-1 is the n-dimensional process noise at time k-1, and its covariance matrix is ​​Q k-1 ; r k is the m-dimensional measurement noise at time k, and its covariance matrix is ​​R k .

[0014] Preferably, in step 2, for the state transfer equation in the model, the photovoltaic side state quantity is added on the basis of the original state quantity of the distribution network, and the specific process of one-step prediction of the dynamic state trajectory of the photovoltaic distribution network is performed based on Holt's two-parameter prediction equation as follows:

[0015] S21. Establish a distributed photovoltaic model;

[0016] S22. Select the DC / AC inverter output voltage amplitude U from the established distributed photovoltaic model F and phase angle ∠θ F As the photovoltaic side state quantity x PV , take the distribution network node voltage amplitude U and phase angle ∠θ as the distribution network side state quantity x net , the photovoltaic side state quantity x PV Combined with the distribution network side state quantity x net Together they constitute the state quantity x of the photovoltaic distribution network, where the photovoltaic side state quantity x PV , grid side state quantity x net The specific expression of the state quantity x of the photovoltaic distribution network is as follows:

[0017] x PV =[U F ,∠θ F ] T

[0018] x net =[U,∠θ] T

[0019] x=[x net ,x PV ] T =[U1,…U j ,∠θ1,…∠θ j ,U F ,∠θ F ] T

[0020] In the formula, j represents the total number of nodes in the distribution network;

[0021] S23, the k-1 time state quantity x of the photovoltaic distribution network k-1 Substituting into Holt's two-parameter exponential smoothing equation, the one-step forecast x of the dynamic state trajectory of the photovoltaic distribution network is obtained. k|k-1 , which includes the state quantity x of the photovoltaic distribution network at time k-1 k-1 The specific expression is as follows:

[0022] x k-1 =[x net ,x PV ] T =[U 1,k-1 ,…U j,k-1 ,∠θ 1,k-1 ,…∠θ j,k-1 ,U F,k-1 ,∠θ F,k-1 ] T .

[0023] Preferably, the distributed photovoltaic model in S21 is a distributed photovoltaic power generation connected to a 10kV medium-voltage distribution network, adopts a two-pole three-phase structure, consists of two power links, a DC / DC direct current boost converter and a DC / AC inverter, and uses LC filtering in the filtering link.

[0024] Preferably, the specific expression of Holt's two-parameter exponential smoothing equation in S23 is as follows:

[0025]

[0026] In the formula, x k|k-1 represents the predicted value of the state quantity x at time k obtained by predicting at time k-1; L k-1 and b k-1 are the horizontal component and tilt component at time k-1 respectively; x k-1 k-2 represents the predicted value of the state quantity x at time k-1 obtained by predicting at time k-2, b k-2 represents the tilt component at time k-2, L k-2 represents the horizontal component at time k-2, α H and β HThey are smoothing parameters for calculating the horizontal component L and the tilt component b, both of which take values ​​in [0,1].

[0027] Preferably, step 3 adds photovoltaic side measurement to the original PMU measurement of the distribution network for the measurement equation in the model, and based on the constraint relationship between the measurement and the state quantity, the specific process of deriving the measurement equation containing the photovoltaic distribution network is as follows:

[0028] S31, obtain the photovoltaic side measurement and distribution network side measurement; the distributed photovoltaic directly uses its own measurement information to select the DC / AC inverter output active power P from the distributed photovoltaic model established in step 2 F , reactive power Q F And the active power P injected into the grid in and reactive power Q in As a photovoltaic side measurement, the specific expression is as follows:

[0029] z PV =[P F ,Q F ,P in ,Q in ] T

[0030] The voltage amplitude U measured by the distribution network node i where the PMU is installed i and phase angle θ U,i and the current amplitude of the connected branch I ij and phase angle θ I,ij The data is measured as the distribution network side quantity, and the specific expression is as follows:

[0031] z net =[U i ,θ U,i ,I ij ,θ I,ij ] T ;

[0032] S32, combining the photovoltaic side quantity measurement obtained in S31 with the distribution network side quantity measurement to form the quantity measurement of the photovoltaic distribution network, and the specific expression is as follows:

[0033] z=[z net ,z PV ] T =[U i ,θ U,i ,I ij ,θ I,ij ,…,P F ,Q F ,P in ,Q in ] T ;

[0034] S33. According to the quantity measurement of the photovoltaic distribution network formed in S32, a measurement equation of the photovoltaic distribution network is generated, including a distribution network side measurement equation and a photovoltaic side measurement equation. The specific expression is as follows:

[0035] The measurement equation on the distribution network side is:

[0036]

[0037] In the formula, g ij and b ij are the conductance and susceptance on the ij branch respectively; G ij and B ij are the real and imaginary parts of the corresponding elements of the node admittance matrix respectively; the superscript z indicates the quantity measurement;

[0038] The photovoltaic side measurement equation is:

[0039]

[0040] Where g and b are the conductance and susceptance of the impedance branch of the LC filter, respectively; g y and b y are the conductance and susceptance on the admittance branch; k is the transformation ratio of the transformer.

[0041] Preferably, in step 4, an improved cubature Kalman filter method is designed, including using singular value SVD decomposition to replace Choleskey decomposition in traditional CKF and introducing adaptive filtering to correct process noise statistical parameters in real time; the specific expression is as follows:

[0042] The specific expression of using singular value SVD decomposition to replace the Choleskey decomposition in traditional CKF is:

[0043]

[0044] The specific expression of the real-time correction of process noise statistical parameters by introducing adaptive filtering is:

[0045]

[0046] In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix, U k-1 Indicates the k-1 moment for P k-1 The orthogonal matrix obtained by SVD decomposition, U k|k-1 represents the k-time response to P k|k-1The orthogonal matrix obtained by SVD decomposition is represents the estimated value of the n-dimensional process noise at time k, d k represents the weighting coefficient at time k, q k-1 represents the n-dimensional process noise at time k-1, represents the estimated value of the state at the previous moment, Represents the estimated value of the state at the current moment.

[0047] Preferably, the specific process of using the improved volumetric Kalman filter method to estimate the dynamic state of the photovoltaic distribution network in step 5 is as follows:

[0048] S51, using the improved volumetric Kalman filter method to solve the general model of dynamic state estimation of the photovoltaic distribution network constructed in step 1, at time k, based on the state estimation value of the previous time and the estimated covariance matrix P k-1 Generate volume sampling points Substituting into the state transfer equation f we get Get the state prediction vector x k|k-1 and its corresponding covariance matrix P k|k-1 , the specific expression is as follows:

[0049]

[0050] In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, q k-1 represents the n-dimensional process noise at time k-1, Q k-1 Represents the covariance matrix of the n-dimensional process noise at time k-1;

[0051] S52: Generate a forecast volume point based on the value obtained in S51 Substitute into the measurement equation h to calculate the measurement prediction volume point And weighted calculation to get the measurement prediction vector z k|k-1 Finally, the measurement prediction error covariance matrix is ​​obtained and the cross covariance matrix The specific expression is as follows:

[0052]

[0053] In the formula, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix of Indicates that the volume point will be predicted Substitute the measured predicted volume point calculated by the measurement equation h, R krepresents the covariance matrix of the measurement noise at time k;

[0054] S53, using the filter gain matrix and the innovation vector e k =z k -z k|k-1 Corrected state prediction vector x k|k-1 , where z k is the measurement vector, and the state estimation vector at time k is obtained and the estimated covariance matrix P k , while updating q in real time k , the specific expression is as follows:

[0055]

[0056] Where P k|k-1 represents the forecast covariance matrix at time k, b represents the forgetting factor, which takes values ​​between [0.95, 0.995], and d k represents the weighting coefficient at time k, b k represents b to the power of k, represents the estimated value of the n-dimensional process noise at time k, q k-1 represents the n-dimensional process noise at time k-1.

[0057] Therefore, the present invention adopts the above-mentioned dynamic state estimation method of photovoltaic distribution network based on improved cubature Kalman, and uses singular value decomposition to solve the problem of filtering interruption or filtering divergence caused by improper error covariance matrix in traditional CKF; introduces adaptive filtering, and corrects the statistical parameters of process noise in real time, so that the estimation results of the amplitude and phase angle of the distribution network node are more accurate, and the accuracy of the dynamic state estimation results is improved; it can also better track the dynamic changes of the connected photovoltaic nodes, and has better adaptability to the system operating conditions after photovoltaic grid connection.

[0058] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 A schematic flow chart of a method for estimating the dynamic state of a photovoltaic power distribution network based on improved volumetric Kalman provided in an embodiment of the present invention;

[0060] Figure 2 A schematic diagram of a distributed photovoltaic topology structure according to an embodiment of the present invention;

[0061] Figure 3 A schematic diagram of the system topology of an improved IEEE33 node in an embodiment of the present invention;

[0062] Figure 4This is a schematic diagram of the change of irradiance G in an embodiment of the present invention;

[0063] Figure 5 Schematic diagram of comparison between the estimated value and the true value of the voltage of node 3 obtained by the method of this invention in an embodiment of the present invention, wherein (a) is a schematic diagram of comparison between the estimated value and the true value of the voltage amplitude of node 3; (b) is a schematic diagram of comparison between the estimated value and the true value of the voltage phase angle of node 3;

[0064] Figure 6 Schematic diagram for comparing the estimated values ​​and true values ​​of the voltages of nodes 3 and 18 obtained by using the method in this article and different methods in an embodiment of the present invention, wherein (a) is a schematic diagram for comparing the estimated value and true value of the voltage amplitude of the third node; (b) is a schematic diagram for comparing the estimated value and true value of the voltage phase angle of the third node; (c) is a schematic diagram for comparing the estimated value and true value of the voltage amplitude of the 18th node; (d) is a schematic diagram for comparing the estimated value and true value of the voltage phase angle of the 18th node. DETAILED DESCRIPTION

[0065] The following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0066] See also Figure 1-6 , a dynamic state estimation method for photovoltaic distribution network based on improved volumetric Kalman includes the following steps:

[0067] Step 1: First, establish a general model for dynamic state estimation of photovoltaic distribution network, including state transfer equation and measurement equation; the specific expression is as follows:

[0068]

[0069] Where: x is the n-dimensional state vector; z is the m-dimensional measurement vector; the subscript k represents the time point; f is the state transfer function, which is the state transfer equation, and expresses the relationship between the state quantity at time k-1 and the state quantity at time k; h is the measurement function, which is the measurement equation, and expresses the relationship between the state quantity at time k and the quantity measurement at time k; q k-1 is the n-dimensional process noise at time k-1, and its covariance matrix is ​​Q k-1 ; r k is the m-dimensional measurement noise at time k, and its covariance matrix is ​​R k In traditional CKF, q k-1 and r k All are assumed to be Gaussian noise with mean 0.

[0070] Step 2: For the state transfer equation in the model, add the photovoltaic side state quantity on the basis of the original state quantity of the distribution network, and perform a one-step forecast of the dynamic state trajectory of the photovoltaic distribution network based on Holt's two-parameter exponential smoothing equation; the specific process is as follows:

[0071] S21. Establish a distributed photovoltaic model; the distributed photovoltaic model is a distributed photovoltaic power generation connected to a 10kV medium-voltage distribution network, which adopts a two-pole three-phase structure and consists of two power links: a DC / DC boost converter and a DC / AC inverter. The filtering link adopts LC filtering. The distributed photovoltaic power generation process is as follows: when the external irradiance is G and the temperature is T, the photovoltaic array composed of solar panels receives sunlight and performs MPPT maximum power point tracking to maximize the output power of the photovoltaic array. At this time, the DC voltage and current output by the photovoltaic array are u pv ,I pv ,u pv ,I pv The DC / DC boost output capacitor voltage u dc , and then output AC voltage u through DC / AC inverter link f , active power is P F , reactive power is Q F , and finally injected into the distribution network node i through LC filtering and transformer. The injected active power and reactive power are P in , Q in Its topological structure is as follows Figure 2 shown.

[0072] S22. Select the DC / AC inverter output voltage amplitude U from the established distributed photovoltaic model F and phase angle ∠θ F As the photovoltaic side state quantity x PV , take the distribution network node voltage amplitude U and phase angle ∠θ as the distribution network side state quantity x net , the photovoltaic side state quantity x PV Combined with the distribution network side state quantity x net Together they constitute the state quantity x of the photovoltaic distribution network, where the photovoltaic side state quantity x PV , grid side state quantity x net The specific expression of the state quantity x of the photovoltaic distribution network is as follows:

[0073] x PV =[U F ,∠θ F ] T

[0074] x net =[U,∠θ] T

[0075] x=[xnet ,x PV ] T =[U1,…U j ,∠θ1,…∠θ j ,U F ,∠θ F ] T

[0076] In the formula, j represents the total number of nodes in the distribution network;

[0077] S23, the k-1 time state quantity x of the photovoltaic distribution network k-1 Substituting into Holt's two-parameter exponential smoothing equation, the one-step forecast x of the dynamic state trajectory of the photovoltaic distribution network is obtained. k|k-1 , which includes the state quantity x of the photovoltaic distribution network at time k-1 k-1 The specific expression is as follows:

[0078] x k-1 =[x net ,x PV ] T =[U 1,k-1 ,…U j,k-1 ,∠θ 1,k-1 ,…∠θ j,k-1 ,U F,k-1 ,∠θ F,k-1 ] T .

[0079] The specific expression of Holt's two-parameter exponential smoothing equation is as follows:

[0080]

[0081] In the formula, x k|k-1 represents the predicted value of the state quantity x at time k obtained by predicting at time k-1; L k-1 and b k-1 are the horizontal component and tilt component at time k-1 respectively; x k-1|k-2 represents the predicted value of the state quantity x at time k-1 obtained by predicting at time k-2, b k-2 represents the tilt component at time k-2, L k-2 represents the horizontal component at time k-2, α H and β H are smoothing parameters for calculating the horizontal component L and the tilt component b, both of which take values ​​in [0,1]. In the present invention, α H =0.85,β H =0.05.

[0082] Step 3: For the measurement equation in the model, add the PV side measurement on the basis of the original PMU measurement of the distribution network, and derive the measurement equation of the PV distribution network based on the constraint relationship between the measurement and the state quantity; due to the large number of branches and economic cost constraints, the actual distribution network is only equipped with synchronized phasor measurement units (PMUs) at some nodes. The data that the PMU can measure are the voltage amplitude and phase angle of the installation node, and the current amplitude and phase angle of the connected branches.

[0083] The specific process is as follows:

[0084] S31, obtain the photovoltaic side measurement and distribution network side measurement; the distributed photovoltaic directly uses its own measurement information to select the DC / AC inverter output active power P from the distributed photovoltaic model established in step 2 F , reactive power Q F And the active power P injected into the grid in and reactive power Q in As a photovoltaic side measurement, the specific expression is as follows:

[0085] z PV =[P F ,Q F ,P in ,Q in ] T

[0086] The voltage amplitude U measured by the distribution network node i where the PMU is installed i and phase angle θ U,i and the current amplitude of the connected branch I ij and phase angle θ I,ij The data is measured as the distribution network side quantity, and the specific expression is as follows:

[0087] z net =[U i ,θ U,i ,I ij ,θ I,ij ] T ;

[0088] S32, combining the photovoltaic side quantity measurement obtained in S31 with the distribution network side quantity measurement to form the quantity measurement of the photovoltaic distribution network, and the specific expression is as follows:

[0089] z=[z net ,z PV ] T =[U i ,θ U,i ,I ij ,θ I,ij ,…,P F,Q F ,P in ,Q in ] T ;

[0090] S33. According to the quantity measurement of the photovoltaic distribution network formed in S32, a measurement equation of the photovoltaic distribution network is generated, including a distribution network side measurement equation and a photovoltaic side measurement equation. The specific expression is as follows:

[0091] The measurement equation on the distribution network side is:

[0092]

[0093] In the formula, g ij and b ij are the conductance and susceptance on the ij branch respectively; G ij and B ij are the real and imaginary parts of the corresponding elements of the node admittance matrix respectively; the superscript z indicates the quantity measurement;

[0094] The photovoltaic side measurement equation is:

[0095]

[0096] Where g and b are the conductance and susceptance of the impedance branch of the LC filter, respectively; g y and b y are the conductance and susceptance on the admittance branch; k is the transformation ratio of the transformer.

[0097] Step 4: Design an improved cubature Kalman filter method, including using singular value SVD decomposition to replace Choleskey decomposition in traditional CKF and introducing adaptive filtering to correct process noise statistical parameters in real time. The specific expression is as follows:

[0098] The specific expression of using singular value SVD decomposition to replace the Choleskey decomposition in traditional CKF is:

[0099]

[0100] The specific basis is as follows:

[0101] The covariance matrix P loses positive definiteness k In form, it is still an n-order symmetric matrix, that is, P T =P, then P T P=PP T =PP is a real symmetric matrix of order n, whose rank is r, so there exists an orthogonal matrix V of order n such that P T P is similar to the diagonal matrix ∑ n :

[0102]

[0103] where λ1,λ2,…λ n YesP T The eigenvalues ​​of P, and Assume λ1≥λ2≥…≥λ r ,λ r+1 =…=λ n =0.

[0104] Let V = [V1, V2] = [v1, v2, …, v r ,v r+1 ,…v n ], then:

[0105]

[0106] Get V1 T P T PV1=Σ r , So PV2=0.

[0107] make but Therefore, the r column vectors of U1 are a unit orthogonal vector group. By expanding it, we obtain the n-order orthogonal matrix U=[U1,U2].

[0108]

[0109] ∵

[0110] ∴

[0111] ∵

[0112] ∴

[0113] ∵PV2=0

[0114] ∴

[0115] ∴

[0116] ∴

[0117] Again∴P T =P

[0118] ∴

[0119]

[0120] ∴

[0121] ∴

[0122] That is, after SVD decomposition,

[0123] The subtraction in the unbiased recursive noise estimator will result in The filter loses semi-positive definiteness and tends to diverge when the system order is high. After the biased estimate is introduced, some terms are omitted, resulting in Large errors reduce the estimation accuracy of the algorithm. Therefore, in the present invention, only the mean value of the process noise of the reaction error is updated in real time, and adaptive filtering is introduced to correct the specific expression of the process noise statistical parameters in real time:

[0124]

[0125] In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix, U k-1 Indicates the k-1 moment for P k-1 The orthogonal matrix obtained by SVD decomposition, U k|k-1 represents the k-time P k|k-1 The orthogonal matrix obtained by SVD decomposition is represents the estimated value of the n-dimensional process noise at time k, d k represents the weighting coefficient at time k, q k-1 represents the n-dimensional process noise at time k-1, represents the estimated value of the state at the previous moment, Represents the estimated value of the state at the current moment.

[0126] Step 5: Use the improved cubature Kalman filter method to estimate the dynamic state of the photovoltaic distribution network. The specific process is as follows:

[0127] S51, using the improved volumetric Kalman filter method to solve the general model of dynamic state estimation of the photovoltaic distribution network constructed in step 1, at time k, based on the state estimation value of the previous time and the estimated covariance matrix P k-1 Generate volume sampling points Substituting into the state transfer equation f we get Get the state prediction vector x k|k-1 and its corresponding covariance matrix P k|k-1 , the specific expression is as follows:

[0128]

[0129] In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, q k-1 represents the n-dimensional process noise at time k-1, Q k-1 Represents the covariance matrix of the n-dimensional process noise at time k-1;

[0130] S52: Generate a forecast volume point based on the value obtained in S51 Substitute into the measurement equation h to calculate the measurement prediction volume point And weighted calculation to get the measurement prediction vector z k|k-1 Finally, the measurement prediction error covariance matrix is ​​obtained and the cross covariance matrix The specific expression is as follows:

[0131]

[0132] In the formula, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix of Indicates that the volume point will be predicted Substitute the measured predicted volume point calculated by the measurement equation h, R k represents the covariance matrix of the measurement noise at time k;

[0133] S53, using the filter gain matrix and the innovation vector e k =z k -z k|k-1 Corrected state prediction vector x k|k-1 , where z k is the measurement vector, and the state estimation vector at time k is obtained and the estimated covariance matrix P k , while updating q in real time k , the specific expression is as follows:

[0134]

[0135] Where P k|k-1 represents the forecast covariance matrix at time k, b represents the forgetting factor, which takes values ​​between [0.95, 0.995], and d k represents the weighting coefficient at time k, b k represents b to the power of k, represents the estimated value of the n-dimensional process noise at time k, q k-1 represents the n-dimensional process noise at time k-1.

[0136] Example

[0137] In order to more clearly demonstrate the technical solution and technical effects provided by the present invention, the performance of the proposed method is verified and analyzed by using a modified IEEE 33-node system. A distributed photovoltaic system with a capacity of 5×550kW is connected to nodes 9, 18, 22, 25, and 33 through transformers. Figure 3 The figure shows the schematic diagram of the system topology of the improved IEEE33 node according to the present invention. In addition, PMUs are configured at 11 nodes, namely nodes 2, 5, 8, 11, 14, 17, 21, 24, 27, 30 and 32, to ensure the system is viable.

[0138] The irradiance G mutation was set to simulate the photovoltaic fluctuations that caused the state quantity mutation. The 18-node photovoltaic system was connected, and G was increased from 1000W / m in 1.2s-1.4s. 2 Reduced to 800W / m 2 ; Connected to the 25-node photovoltaic, G from 1000W / m in 1.6s-1.8s 2 Reduced to 800W / m 2 ; The rest of the PVs are operated under standard conditions (G = 1000W / m 2 , T=298k). G change curve is as follows Figure 4 As shown in the figure. The measurement sampling interval is 0.01s and the simulation time is 10s. The state value obtained by the system time domain simulation is taken as the true value, and the covariance matrix of the process and measurement noise is set to 1×10 -6 .

[0139] The root mean square error (RMSE) of the estimated results of the system state variables at each moment is defined as a quantitative evaluation index of the estimation effect. The RMSE at the kth discrete moment is defined as:

[0140]

[0141] In the formula, and are the estimated value and true value of the i-th state quantity at time k, and n is the number of state quantities.

[0142] like Figure 5 The figure shows a schematic diagram comparing the estimated value and the true value of the state quantity of the third node in the example of the present invention, wherein (a) is a schematic diagram comparing the estimated value and the true value of the voltage amplitude of the third node; (b) is a schematic diagram comparing the estimated value and the true value of the voltage phase angle of the third node. Due to the change in irradiance, the voltage of the node of the distribution network with photovoltaic access fluctuates to varying degrees, but the state estimation value obtained by using the improved cubature Kalman filter method proposed in the present invention for dynamic state estimation can accurately track the change of the true value and ensure a high state estimation accuracy.

[0143] On this basis, the method proposed in the present invention is compared with other methods. Figure 6 The diagram shown is a comparison diagram of the voltage state estimation values ​​of the 3rd and 18th nodes in the example of the present invention and other methods, wherein (a) is a comparison diagram of the estimated value and the true value of the voltage amplitude of the 3rd node; (b) is a comparison diagram of the estimated value and the true value of the voltage phase angle of the 3rd node; (c) is a comparison diagram of the estimated value and the true value of the voltage amplitude of the 18th node; (d) is a comparison diagram of the estimated value and the true value of the voltage phase angle of the 18th node. It can be seen that the accuracy of the estimation value of the method proposed in this article is significantly higher than that of other methods. The following Table 1 further compares the mean RMSE of the voltage estimation values ​​of all nodes by two methods, namely the traditional method SVD-CKF and the method SVD-ACKF adopted in this application. It can be seen that after taking into account the dynamic component state and adopting the improved volumetric Kalman filtering method proposed by the present invention, the node estimation value is closer to the true value and the error is smaller.

[0144] Table 1 Root mean square error comparison

[0145] algorithm Average RMSE of node voltage amplitude estimation Average RMSE of node voltage phase angle estimation SVD-CKF 0.0136 0.0864 SVD-ACKF 0.0061 0.0378

[0146] Therefore, the present invention adopts the above-mentioned dynamic state estimation method of photovoltaic distribution network based on improved cubature Kalman, integrates the constraint relationship of dynamic and algebraic state of photovoltaic distribution network into the node voltage state transition equation, reflects the main driving factors of node voltage change, and makes the node voltage transfer model have clear physical meaning; at the same time, SVD decomposition is used to replace Choleskey decomposition in CKF and adaptive filtering is introduced to update the process noise mean in real time to improve the traditional CKF, thereby improving the accuracy of dynamic state estimation results, and having good adaptability to the system operation conditions of photovoltaic grid-connected.

[0147] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A dynamic state estimation method for photovoltaic distribution network based on improved volumetric Kalman, characterized in that: The following steps are involved: Step 1: First, establish a general model for dynamic state estimation of photovoltaic distribution network, including state transfer equation and measurement equation; Step 2: For the state transfer equation in the model, add the photovoltaic side state quantity on the basis of the original state quantity of the distribution network, and perform a one-step forecast of the dynamic state trajectory of the photovoltaic distribution network based on Holt's two-parameter exponential smoothing equation; Step 3: For the measurement equation in the model, add the PV side measurement on the basis of the original PMU measurement of the distribution network, and derive the measurement equation of the PV distribution network based on the constraint relationship between the measurement and the state quantity; Step 4: design an improved cubature Kalman filter method; Step 5: Use the improved volumetric Kalman filter method to estimate the dynamic state of the photovoltaic distribution network; In step 2, for the state transfer equation in the model, the PV side state quantity is added to the original state quantity of the distribution network. Based on Holt's two-parameter prediction equation, the specific process of one-step prediction of the dynamic state trajectory of the PV distribution network is as follows: S21. Establish a distributed photovoltaic model; S22. Select the DC / AC inverter output voltage amplitude U from the established distributed photovoltaic model F and phase angle ∠θ F As the photovoltaic side state quantity x PV , take the distribution network node voltage amplitude U and phase angle ∠θ as the distribution network side state quantity x net , the photovoltaic side state quantity x PV Combined with the distribution network side state quantity x net Together they constitute the state quantity x of the photovoltaic distribution network, where the photovoltaic side state quantity x PV , grid side state quantity x net The specific expression of the state quantity x of the photovoltaic distribution network is as follows: x PV =[U F ,∠θ F ] T x net =[U,∠θ] T x=[x net ,x PV ] T =[U1,…U j ,∠θ1,…∠θ j ,U F ,∠θ F ] T In the formula, j represents the total number of nodes in the distribution network; S23, the k-1 time state quantity x of the photovoltaic distribution network k-1 Substituting into Holt's two-parameter exponential smoothing equation, the one-step forecast x of the dynamic state trajectory of the photovoltaic distribution network is obtained. k|k-1 , which includes the state quantity x of the photovoltaic distribution network at time k-1 k-1 The specific expression is as follows: x k-1 =[x net ,x PV ] T =[U 1,k-1 ,…U j,k-1 ,∠θ 1,k-1 ,…∠θ j,k-1 ,U F,k-1 ,∠θ F,k-1 ] T ; Step 3: For the measurement equation in the model, add the PV side measurement on the basis of the original PMU measurement of the distribution network. Based on the constraint relationship between the measurement and the state quantity, the specific process of deriving the measurement equation containing the PV distribution network is as follows: S31, obtain the photovoltaic side measurement and distribution network side measurement; the distributed photovoltaic directly uses its own measurement information to select the DC / AC inverter output active power P from the distributed photovoltaic model established in step 2 F , reactive power Q F And the active power P injected into the grid in and reactive power Q in As a photovoltaic side measurement, the specific expression is as follows: z PV =[P F ,Q F ,P in ,Q in ] T The voltage amplitude U measured by the distribution network node i where the PMU is installed i and phase angle θ U,i and the current amplitude of the connected branch I ij and phase angle θ I,ij The data is measured as the distribution network side quantity, and the specific expression is as follows: z net =[U i ,θ U,i ,I ij ,θ I,ij ] T ; S32, combining the photovoltaic side quantity measurement obtained in S31 with the distribution network side quantity measurement to form the quantity measurement of the photovoltaic distribution network, and the specific expression is as follows: z=[z net ,z PV ] T =[U i ,θ U,i ,I ij ,θ I,ij ,…,P F ,Q F ,P in ,Q in ] T ; S33. According to the quantity measurement of the photovoltaic distribution network formed in S32, a measurement equation of the photovoltaic distribution network is generated, including a distribution network side measurement equation and a photovoltaic side measurement equation. The specific expression is as follows: The measurement equation on the distribution network side is: In the formula, g ij and b ij are the conductance and susceptance on the ij branch respectively; G ij and B ij are the real and imaginary parts of the corresponding elements of the node admittance matrix respectively; the superscript z indicates the quantity measurement; The photovoltaic side measurement equation is: Where g and b are the conductance and susceptance of the impedance branch of the LC filter, respectively; g y and b y are the conductance and susceptance on the admittance branch; k is the transformation ratio of the transformer; In step 4, an improved cubature Kalman filter method is designed, which includes using singular value SVD decomposition to replace Choleskey decomposition in CKF and introducing adaptive filtering to correct process noise statistical parameters in real time; the specific expression is as follows: The specific expression of using singular value SVD decomposition to replace Choleskey decomposition in CKF is: The specific expression of the real-time correction of process noise statistical parameters by introducing adaptive filtering is: In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix, U k-1 Indicates the k-1 moment for P k-1 The orthogonal matrix obtained by SVD decomposition, U kk-1 represents the k-time P k|k-1 The orthogonal matrix obtained by SVD decomposition is represents the estimated value of the n-dimensional process noise at time k, d k represents the weighting coefficient at time k, q k-1 represents the n-dimensional process noise at time k-1, represents the estimated value of the state at the previous moment, Represents the estimated value of the state at the current moment.

2. The method for dynamic state estimation of photovoltaic power distribution network based on improved volumetric Kalman according to claim 1, characterized in that: The specific expression of the general model for dynamic state estimation of photovoltaic distribution network established in step 1 is as follows: Where: x is the n-dimensional state vector; z is the m-dimensional measurement vector; the subscript k represents the time point; f is the state transfer function, which is the state transfer equation, and expresses the relationship between the state quantity at time k-1 and the state quantity at time k; h is the measurement function, which is the measurement equation, and expresses the relationship between the state quantity at time k and the quantity measurement at time k; q k-1 is the n-dimensional process noise at time k-1, and its covariance matrix is ​​Q k-1 ; r k is the m-dimensional measurement noise at time k, and its covariance matrix is ​​R k .

3. The method for dynamic state estimation of photovoltaic power distribution network based on improved volumetric Kalman according to claim 2, characterized in that: The distributed photovoltaic model in S21 is a distributed photovoltaic power generation connected to a 10kV medium-voltage distribution network. It adopts a two-pole three-phase structure and consists of two power links: a DC / DC boost converter and a DC / AC inverter. The filtering link adopts LC filtering.

4. The method for dynamic state estimation of photovoltaic power distribution network based on improved volumetric Kalman according to claim 3, characterized in that: The specific expression of Holt's two-parameter exponential smoothing equation in S23 is as follows: In the formula, x kk-1 represents the predicted value of the state quantity x at time k obtained by predicting at time k-1; L k-1 and b k-1 are the horizontal component and tilt component at time k-1 respectively; x k-1k-2 represents the predicted value of the state quantity x at time k-1 obtained by predicting at time k-2, b k-2 represents the tilt component at time k-2, L k-2 represents the horizontal component at time k-2, α H and β H They are smoothing parameters for calculating the horizontal component L and the tilt component b, both of which take values ​​in [0,1].

5. The method for dynamic state estimation of photovoltaic power distribution network based on improved volumetric Kalman according to claim 4, characterized in that: The specific process of using the improved cubature Kalman filter method to estimate the dynamic state of the photovoltaic distribution network in step 5 is as follows: S51, using the improved volumetric Kalman filter method to solve the general model of dynamic state estimation of photovoltaic distribution network constructed in step 1, at time k, based on the state estimation value of the previous time and the estimated covariance matrix P k-1 Generate volume sampling points Substituting into the state transfer equation f we get Get the state prediction vector x kk-1 and its corresponding covariance matrix P k|k-1 , the specific expression is as follows: In the formula, S k-1 Represents the k-1 time error covariance matrix P k-1 The square root matrix, n represents the dimension of the state quantity x, q k-1 represents the n-dimensional process noise at time k-1, Q k-1 Represents the covariance matrix of the n-dimensional process noise at time k-1; S52: Generate a forecast volume point based on the value obtained in S51 Substitute into the measurement equation h to calculate the measurement prediction volume point And weighted calculation to get the measurement prediction vector z kk-1 Finally, we get the measurement prediction error covariance matrix and the cross covariance matrix The specific expression is as follows: In the formula, S k|k-1 represents the k-time forecast covariance matrix P k|k-1 The square root matrix of Indicates that the volume point will be predicted Substitute the measured predicted volume point calculated by the measurement equation h, R k represents the covariance matrix of the measurement noise at time k; S53, using the filter gain matrix and the innovation vector e k =z k -z k|k-1 Corrected state prediction vector x kk-1 , where z k is the measurement vector, and the state estimation vector at time k is obtained and the estimated covariance matrix P k , while updating q in real time k , the specific expression is as follows: Where P k|k-1 represents the forecast covariance matrix at time k, b represents the forgetting factor, which takes values ​​between [0.95, 0.995], and d k represents the weighting coefficient at time k, b k represents b to the power of k, represents the estimated value of the n-dimensional process noise at time k, q k-1 represents the n-dimensional process noise at time k-1.

Images