Harmonic impedance calculation method and system based on data fusion

Abstract

The invention discloses a harmonic impedance calculation method and system based on data fusion. The method comprises the following steps: respectively using a FastICA algorithm and a fluctuation quantity method to solve the harmonic impedance of a measurement signal; a prior state prediction value # imgabs0 # and a measurement prediction value # imgabs1 # of a measurement signal are solved, and an input noise covariance matrix Q and a measurement noise covariance matrix R are obtained through a hyper-parameter estimation method; a priori error covariance matrix # imgabs2 is calculated; a Kalman gain Kk is calculated; updating the error covariance matrix Pk; returning to the prior error covariance matrix calculation step for iterative calculation, stopping iteration when a preset iteration value is reached, and outputting a final prediction value; according to the method, the hyper-parameter estimation method of the covariance matrix of the input noise and the measurement noise is provided so as to achieve accurate evaluation of the covariance matrix of the input noise and the measurement noise, the advantages of the FastICA and the fluctuation quantity method are effectively combined, the application range is widened, and the calculation precision is improved.

Description

Technical Field

[0001] The present invention relates to a harmonic impedance calculation method and system, and in particular to a harmonic impedance calculation method and system based on data fusion. Background Art

[0002] As harmonic sources in the power grid become more diverse and complex, simply locating the dominant harmonic source can no longer meet industrial needs. Accurately quantifying the responsibility of each harmonic source is conducive to formulating a harmonic reward and punishment mechanism and guiding the development of harmonic governance work.

[0003] Most of the existing methods are based on the fact that the background harmonics are relatively stable and the user-side harmonic impedance |Z c |》System measurement of harmonic impedance | Z u |Assumption. These two assumptions are usually valid for traditional nonlinear users, but for the harmonic sources of new energy sites, on the one hand, the system side contains high-power harmonic sources, which greatly increases the risk of background harmonic instability; on the other hand, the new energy site is equipped with an AC filter, which reduces the amplitude of the harmonic impedance on the new energy site side under certain subharmonics and is no longer much larger than the harmonic impedance on the system side. As the harmonic impedance on the new energy site side decreases, it becomes more difficult to solve the harmonic impedance on the system side, and at this time, Z can no longer be ignored in the quantification of harmonic responsibility. c In order to solve the problem of harmonic source tracing in new energy stations, the harmonic impedance on both sides needs to be calculated simultaneously. In summary, the assumptions required by traditional algorithms are no longer valid in the problem of harmonic source tracing in new energy stations, and it is difficult for existing algorithms to quantify the harmonic responsibility of new energy stations. Therefore, how to overcome the limitations of single-point harmonic responsibility quantification methods in practical applications, and consider the instability of background harmonics of harmonic sources in new energy stations, |Z c |>>|Z u | is not true, and so on. It is an urgent problem to propose a harmonic responsibility quantification algorithm that adapts to new scenarios. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to provide a harmonic impedance calculation method and system based on data fusion that can improve calculation accuracy.

[0005] Technical solution: A harmonic impedance calculation method based on data fusion described in the present invention includes:

[0006] (1) The FastICA algorithm and the wave quantity method are used to obtain the harmonic impedance of the measurement signal;

[0007] (2) Obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method;

[0008] (3) Calculate the prior error covariance matrix

[0009] (4) Calculate the Kalman gain K k ;

[0010] (5) Update the error covariance matrix P k ;

[0011] (6) Return to step (3) to iterate the calculation of the prior error covariance matrix. When the preset iteration value is reached, stop the iteration and output the final prediction value.

[0012] Furthermore, the step (1) of obtaining the harmonic impedance of the measurement signal using the FastICA algorithm includes:

[0013] Using entropy to represent the non-Gaussian measure, negative entropy is defined as:

[0014] N g (Y) = H (Y G )-H(Y)(1)

[0015] in

[0016] H(Y)=-∫P Y (ξ)lgP Y (ξ)dξ(2)

[0017] Where Y G and Y are Gaussian random variables with equal variances, P Y (ξ) is the probability density function of the random variable Y, ξ represents the value of the random variable Y; N g (Y) is a measure of non-Gaussianity; when independent component separation is completed, N g (Y) has the largest value; the approximate calculation formula for negative entropy is:

[0018] N g (Y) = {E[g(Y)] - E[g(Y G )]} 2 (3)

[0019] in

[0020] g=tan(a1y),1≤a1≤2(4)

[0021] In the formula, E is the average value, y is the sample value; g is the nonlinear function;

[0022] The measurement signal is preprocessed, and the preprocessing includes decentralization and whitening processing; the decentralization is expressed as:

[0023]

[0024] Among them, x(t i) represents the original signal at time point t i The value of represents the signal value after decentralized processing, and N represents the total number of signal samples;

[0025] The mathematical expression of the whitening process is:

[0026]

[0027] in

[0028]

[0029] In the formula, Z(t i ) is the signal after whitening processing, matrix R is the eigenvector matrix of the covariance matrix, matrix V is the whitening matrix, matrix U is the eigenvalue matrix of the covariance matrix, and the measurement signal is whitened to obtain matrix V, so that the measurement signal satisfies E[ZZ T ]=I, I is the identity matrix, T represents transpose;

[0030] Perform iterative calculations to maximize the negative entropy; the maximum approximate value of the negative entropy must meet the following conditions:

[0031] R[g(V T Z)]+βV=0(8)

[0032] in

[0033] β=R[V1 T Zg(V1 T Z)](9)

[0034] In the formula, V1 is the optimized V value; the Newton iteration method is used to solve the above formula to obtain the optimal solution, and the mathematical expression is:

[0035]

[0036] In the formula, i is the number of iterations, V i represents the value of V at the i-th iteration, ‖‖ represents the norm, and g′ is the derivative of g; when ‖V i+1 ‖When it tends to be constant, the separation of independent components is completed;

[0037] The harmonic analysis Norton equivalent equation is converted into complex form, as shown in equation (11);

[0038]

[0039] In the formula, ΔU pcc-x and ΔI pcc-x are the real parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common busbar; ΔU pcc-y and ΔI pcc-yare the imaginary parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common bus; ΔI s-x and ΔI c-x are the real parts of the harmonic source current vector fluctuations on the system side and the user side respectively; ΔI s-y and ΔI c-y are the imaginary parts of the harmonic source current vector fluctuations on the system side and the user side respectively; a mn (m,n=1,2,3,4) is the harmonic source current fluctuation coefficient;

[0040] Convert equation (11) into matrix form, expressed as:

[0041] X=AS (12)

[0042] In the formula, X = [ΔU pcc-x , ΔU pcc-y , ΔI pcc-x , ΔI pcc-y ] T ; S = [ΔI c-x , ΔI c-y , ΔI s-x , ΔI s-y ] T ; A is the coefficient matrix;

[0043] According to formula (12), after decoupling, four independent components I are obtained: n shape information, n = 1, 2, 3, 4;

[0044] set up in is a non-zero real number, n = 1, 2, 3, 4; use the least squares method to find I n The mixing coefficient K m =[k 1m ,k 2m ,k 3m ,k 4m ], mathematically expressed as:

[0045] K m =(I T I) -1 I T X m (13)

[0046] Among them, X m is the mth vector in the X matrix,

[0047]

[0048] k mn The relationship with the harmonic impedance on the system side and the user side is as follows:

[0049]

[0050] The solution of formula (15) is as follows:

[0051]

[0052] In formula (16), Z s-x is the real part of the harmonic impedance on the system side; Z s-y is the imaginary part of the harmonic impedance on the system side; Z c-x is the real part of the user-side harmonic impedance; Z c-y is the imaginary part of the harmonic impedance at the user side;

[0053] By judging the positive or negative of the real part of the harmonic impedance, the solution on the system side is inferred. When j = 1, solve the coefficient equation and obtain the solution corresponding to I1, which is set as b 11 and b 21 , when j = 2, solve the coefficient equation and obtain the corresponding solution of I2, set it as b 12 and b 22 ; Then the harmonic impedance on the system side is expressed as:

[0054]

[0055] Furthermore, the step (1) of obtaining the harmonic impedance of the measurement signal using the wave quantity method includes:

[0056] Considering the changes of harmonic current on the system side and harmonic current on the user side, the mathematical expression is:

[0057]

[0058]

[0059] is the voltage fluctuation at the PCC point, is the current fluctuation at the PCC point, The fluctuation amount of current injected into the user side, The fluctuation of the injected current is measured by the system, Z s represents the system impedance, Z c Indicates user impedance;

[0060] Definition of PCC point harmonic fluctuation volt-ampere characteristic K hVA It is the ratio of the harmonic voltage fluctuation at the PCC point to the harmonic current fluctuation at the common bus:

[0061]

[0062] Analyze the PCC point fluctuations caused by fluctuations on the system side or user side:

[0063]

[0064] K hVA The slope of the harmonic volt-ampere characteristic curve is corresponding to the harmonic impedance on both sides of the PPC point estimated by the ratio of the harmonic current and harmonic voltage fluctuations at the PCC point; when the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is positive, the result is the harmonic impedance on the system side, and when the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is negative, the result is the negative number of the harmonic impedance on the user side;

[0065] Since the load capacity on the user side is much smaller than the short-circuit capacity on the system side, there is Z s Z c , transform formula (19) into:

[0066]

[0067] In the formula, e sys is the error caused by background harmonic fluctuation.

[0068] Furthermore, the step (2) of obtaining the prior state prediction value of the measurement signal is and measured predicted values include:

[0069] Consider the state space model of the linear discrete dynamic system as:

[0070] X(k)=AX(k-1)+ΓW(k-1) (23)

[0071] Z(k)=HX(k)+V(k)(24)

[0072] Wherein, k is the discrete time; X(k) and Z(k) are the state variable and the measurement variable respectively; W(k-1) and V(k) are the input noise and the measurement noise respectively; Γ is the noise driving matrix; A is the state transfer matrix; H is the observation matrix; Equation (23) is the state equation, and Equation (24) is the observation equation;

[0073] Assume that the input noise and measurement noise are independent of each other, their mean is zero, their covariance matrices are Q and R respectively, and their statistical characteristics are:

[0074]

[0075] Where k is the discrete time, j is the discrete time that is not synchronized with k, W(k) and V(k) are the input noise and measurement noise, Q and R are the covariance matrix of the input noise and the covariance matrix of the measurement noise, P(W) and P(V) represent the probability distribution of the input noise and the output noise, δ kj is the Kronecker symbol, which is 1 when k = j, indicating that the noise at the same time has non-zero correlation and the noise at different times is independent of each other;

[0076] Considering X(k) and Z(k) as the calculation results of the two methods respectively, the Kalman filter recursive formula is as follows:

[0077] When the noise error is not considered, the prior prediction value is obtained, which is mathematically expressed as:

[0078]

[0079]

[0080] In the formula, is the k-th prior state prediction value; To measure the predicted value.

[0081] Furthermore, the step (3) comprises:

[0082] The noise is eliminated by Kalman filtering, and the final prediction value of Kalman filtering at time k is expressed as:

[0083]

[0084] In the formula, K k is the Kalman gain, and K k ∈[0,H -1 ]; According to formula (28), when K k = 0, The final prediction value is equal to the prior state prediction value. k =H -1 hour, The final predicted value is equal to the measured predicted value;

[0085] Take the prediction error e k and prior error for:

[0086]

[0087]

[0088] When the variance of the error is minimized, the final predicted value is closest to the true value x k , the error covariance matrix P k for:

[0089]

[0090] Combining equation (29) and equation (31), we get:

[0091]

[0092] make:

[0093]

[0094] in

[0095]

[0096] Combining formula (25) to obtain the prior error covariance matrix:

[0097]

[0098] Furthermore, the step (4) comprises:

[0099] Calculate the error covariance matrix P k :

[0100]

[0101] The covariance matrix P of the error k The trace is:

[0102]

[0103] When tr(P k ) is minimum, the prediction error e k The variance of tr(P k ) is derived, we have:

[0104]

[0105] Pick

[0106]

[0107] Then there is

[0108]

[0109] Furthermore, the step (2) of obtaining the input noise covariance matrix Q and the measurement noise covariance matrix R by the hyperparameter estimation method includes:

[0110] Assume that the calculation error e is:

[0111] e=Z s-j -Z s-z (41)

[0112] In the formula, Z s-j is the calculated value of harmonic impedance on the system side; Z s-z is the true value of the harmonic impedance on the system side; the variance of the error e is:

[0113] VAR(e)=VAR(Z s-j -Z s-z )(42)

[0114] Where VAR() represents the calculation of variance; in the process of harmonic impedance calculation, Z s-z is a constant; according to the variance operation theorem:

[0115] VAR(e)=VAR(Z s-j -Z s-z )=VAR(Z s-j )(43)

[0116] According to formula (43), the values ​​of the input noise covariance matrix Q and the measurement noise covariance matrix R are determined by the variance of the results calculated by FastICA and the fluctuation method.

[0117] Furthermore, the step (5) comprises:

[0118] Update the error covariance matrix P k , the formula is:

[0119]

[0120] Based on the same inventive concept, the present invention also provides a harmonic impedance calculation system based on data fusion, comprising:

[0121] A preliminary calculation module for obtaining the harmonic impedance of the measurement signal using the FastICA algorithm and the wave quantity method respectively;

[0122] Parameter calculation module, used to obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method;

[0123] Correction module for calculating the prior error covariance matrix

[0124] Gain calculation module, used to calculate the Kalman gain K k ;

[0125] Update module, used to update the error covariance matrix P k ;

[0126] The prediction output module is used to return the iterative calculation of the parameter calculation module, stop the iteration when the preset iteration value is reached, and output the final prediction value.

[0127] Based on the same inventive concept, the present invention also provides a computing device, comprising: one or more processors, one or more memories, and one or more programs, wherein the programs are stored in the memories and configured to be executed by the processors, and when the programs are loaded into the processors, the steps of the harmonic impedance calculation method based on data fusion according to any one of the above items are implemented.

[0128] Beneficial effects: Compared with the prior art, the present invention combines the characteristics of FastICA with the wave quantity method, and still has higher calculation accuracy when the background harmonic fluctuates violently and the system side and the user side are correlated; the present invention proposes a hyperparameter estimation method for the covariance matrix of input noise and measurement noise, which can accurately estimate the covariance matrix of input noise and measurement noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0129] Figure 1 is a flow chart of a method according to an embodiment of the present invention;

[0130] Figure 2 A schematic diagram of the solution steps of the FastICA method according to an embodiment of the present invention;

[0131] Figure 3 A schematic diagram of the calculation steps of the FastICA method according to an embodiment of the present invention;

[0132] Figure 4 Schematic diagram of iterative calculation according to an embodiment of the present invention. DETAILED DESCRIPTION

[0133] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0134] As attached Figure 1 As shown, the harmonic impedance calculation method based on data fusion of this embodiment includes:

[0135] (1) The FastICA algorithm and the wave quantity method are used to obtain the harmonic impedance of the measurement signal;

[0136] (2) Obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method;

[0137] (3) Calculate the prior error covariance matrix

[0138] (4) Calculate the Kalman gain K k ;

[0139] (5) Update the error covariance matrix P k ;

[0140] (6) Return to step (3) to iterate the calculation of the prior error covariance matrix. When the preset iteration value is reached, stop the iteration and output the final prediction value.

[0141] Specifically, in step (1), the FastICA algorithm uses the Newton iteration method to analyze and process big data. It is a fixed-point iterative optimization algorithm that minimizes the mutual information of estimated components. When the independent components are separated, the non-Gaussianity measure of the signal is maximized. It has fast convergence speed, high robustness, no step size parameter, and can separate non-Gaussian independent components.

[0142] Entropy is used to represent non-Gaussian measurement. The modified form of entropy, namely negative entropy, is often used, which is defined as:

[0143] N g (Y) = H (Y G )-H(Y)(1)

[0144] in

[0145] H(Y)=-∫P Y (ξ)lgP Y (ξ)dξ(2)

[0146] Where Y G and Y are Gaussian random variables with equal variances, P Y (ξ) is the probability density function of the random variable Y; N g (Y) is the non-Gaussian measure. When the independent component separation is completed, N g (Y) has the largest value. The approximate calculation formula for negative entropy is:

[0147] N g (Y) = {E[g(Y)] - E[g(Y G )]} 2 (3)

[0148] in

[0149] g=tan(a1y),1≤a1≤2(4)

[0150] In the formula, E is the average value; g is a nonlinear function. Usually, a1=1; g=yexp(-y 2 / 2); g = y 3 .

[0151] In order to reduce the amount of calculation and the speed of signal separation, the measured signal needs to be preprocessed before using FastICA, mainly including decentralization and whitening. Decentralization reduces the complexity of calculation, speeds up the separation and convergence of FastICA, and avoids errors caused by large data differences and different dimensions. After the signal is decentralized, a new set of data with a mean of zero is obtained, which can be expressed as:

[0152]

[0153] Whitening can eliminate the correlation between the measurement signals, so that the measurement signals have unit variance, which is mathematically expressed as:

[0154]

[0155] in

[0156]

[0157] Where, matrix R is the eigenvector matrix of the covariance matrix; matrix V is the whitening matrix; and matrix U is the eigenvalue matrix of the covariance matrix. The purpose of signal whitening is to obtain matrix V so that the measured signal satisfies E[ZZ T ]=I, I is the unit matrix.

[0158] After being decentralized and whitened, the measured signal becomes a variable with unit variance and mean value of zero, and each variable is orthogonal to each other. Then iterative calculation is performed to maximize the negative entropy. The maximum approximate value of the negative entropy must meet the following conditions:

[0159] R[g(V T Z)]+βV=0 (8)

[0160] in

[0161] β=R[V1 T Zg(V1 T Z)] (9)

[0162] In the formula, V1 is the optimized V value. Then the Newton iteration method is used to solve the above formula to obtain the optimal solution, and its mathematical expression is:

[0163]

[0164] Where i is the number of iterations. i+1 When ‖ tends to be constant, the separation of independent components is completed. In summary, the solution steps of FastICA are as follows: Figure 2 shown.

[0165] When FastICA is used for harmonic impedance calculation, the harmonic analysis Norton equivalent equation is converted into a complex form for the convenience of calculation, as shown in formula (11).

[0166] In formula (11), ΔU pcc-x and ΔI pcc-x are the real parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common busbar; ΔU pcc-y and ΔI pcc-y are the imaginary parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common bus; ΔI s-x and ΔI c-xare the real parts of the harmonic source current vector fluctuations on the system side and the user side respectively; ΔI s-y and ΔI c-y are the imaginary parts of the harmonic source current vector fluctuations on the system side and the user side respectively; a ij (i,j=1,2,3,4) is the harmonic source current fluctuation coefficient.

[0167]

[0168] Converting the above formula into matrix form, it can be expressed as:

[0169] X=AS (12)

[0170] In the formula, X = [ΔU pcc-x , ΔU pcc-y , ΔI pcc-x , ΔI pcc-y ] T ; S = [ΔI c-x , ΔI c-y , ΔI s-x , ΔI s-y ] T ; A is the coefficient matrix.

[0171] According to formula (12), after decoupling, four independent components I are obtained: n shape information, n = 1, 2, 3, 4; but its amplitude, phase and corresponding relationship cannot be determined, so the decoupled independent components need to be processed further.

[0172] set up in is a non-zero real number, n = 1, 2, 3, 4; use the least squares method to find I n The mixing coefficient K m =[k 1m ,k 2m ,k 3m ,k 4m ], mathematically expressed as:

[0173] K m =(I T I) -1 I T X m (13)

[0174] Among them, X m is the mth vector in the X matrix,

[0175]

[0176] Comparison found that k mn The relationship with the harmonic impedance on the system side and the user side is as follows:

[0177]

[0178] The solution of equation (15) can be expressed by equation (16).

[0179] In formula (16), Z s-x is the real part of the harmonic impedance on the system side; Z s-y is the imaginary part of the harmonic impedance on the system side; Z c-x is the real part of the user-side harmonic impedance; Z c-y is the imaginary part of the harmonic impedance on the user side.

[0180]

[0181] Generally speaking, the real part of harmonic impedance is positive, so the solution on the system side can be inferred by judging the sign of the real part of harmonic impedance. When j = 1, solve the coefficient equation and get the solution corresponding to I1, let it be b 11 and b 21 , the corresponding solution of I2 can be obtained by the same logic, let it be b 12 and b 22 . Then the harmonic impedance on the system side can be expressed as:

[0182]

[0183] Combined with the above, the calculation steps are as follows Figure 3 shown.

[0184] The Fluctuation Volume Method:

[0185] The principle of the fluctuation quantity method is to obtain the harmonic impedance on the system side and the user side by analyzing the fluctuation quantity ratio characteristics of the harmonic voltage at the PCC point and the harmonic current at the common bus. In the actual operation of the power grid, the operation mode of the system will change, and the user load will also fluctuate, which will cause the harmonic voltage and harmonic current at the PCC point to fluctuate. The harmonic impedance on the system side depends on the short-circuit capacity of the system, and its size can be regarded as constant in a short time. Changes in the harmonic current on the user side can reflect changes in the harmonic impedance on the user side. Therefore, when analyzing the fluctuations on the system side and the user side, only the changes in the harmonic current on the system side and the harmonic current on the user side can be considered, and the mathematical expression is:

[0186]

[0187]

[0188] is the voltage fluctuation at the PCC point, is the current fluctuation at the PCC point, The fluctuation amount of current injected into the user side, The fluctuation of the injected current is measured by the system, Zs represents the system impedance, Z c Indicates user impedance;

[0189] Definition of PCC point harmonic fluctuation volt-ampere characteristic K hVA It is the ratio of the harmonic voltage fluctuation at the PCC point to the harmonic current fluctuation at the common bus, and its value is:

[0190]

[0191] When analyzing PCC point fluctuations caused by fluctuations on the system side or user side:

[0192]

[0193] K hVA The geometric meaning of is the slope of the corresponding harmonic volt-ampere characteristic curve. From the above formula, it can be seen that the harmonic impedance on both sides of the PPC point can be estimated by the ratio of the harmonic current and harmonic voltage fluctuations at the PCC point. According to the sign of the harmonic fluctuation volt-ampere characteristic at the PCC point, it can be determined whether the result is the harmonic impedance on the system side or the harmonic impedance on the user side. When the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is positive, the result is the harmonic impedance on the system side, and when the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is negative, the result is the negative of the harmonic impedance on the user side.

[0194] Since the load capacity on the user side is generally much smaller than the short-circuit capacity on the system side, there is Z s <<Z c , then formula (19) can be transformed into:

[0195]

[0196] In the formula, e sys is the error caused by background harmonic fluctuation.

[0197] In step (2), Kalman filtering technology is used. Kalman filtering technology is a kind of optimal estimation method. The optimal estimation value is obtained by minimizing the variance of the estimation error. The principle is as follows.

[0198] Consider the state space model of the linear discrete dynamic system as:

[0199] X(k)=AX(k-1)+ΓW(k-1) (23)

[0200] Z(k)=HX(k)+V(k) (24)

[0201] Wherein, k is the discrete time; X(k) and Z(k) are the state variable and the measurement variable respectively; W(k-1) and V(k) are the input noise and the measurement noise respectively; Γ is the noise driving matrix; A is the state transfer matrix; H is the observation matrix; Equation (23) is the state equation, and Equation (24) is the observation equation;

[0202] Assume that the input noise and measurement noise are independent of each other, their mean is zero, and their covariance matrices are Q and R respectively, and their statistical characteristics are:

[0203]

[0204] Where k is the discrete time, j is the discrete time that is not synchronized with k, W(k) and V(k) are the input noise and measurement noise, Q and R are the covariance matrix of the input noise and the covariance matrix of the measurement noise, P(W) and P(V) represent the probability distribution of the input noise and the output noise, δ kj is the Kronecker symbol, which is 1 when k = j, indicating that the noise at the same time has non-zero correlation and the noise at different times is independent of each other;

[0205] When considering the calculation of harmonic impedance, X(k) and Z(k) can be regarded as the calculation results of two methods respectively. The Kalman filter recursive formula is as follows.

[0206] When the noise error is not considered, the prior prediction value can be obtained, which is mathematically expressed as:

[0207]

[0208]

[0209] In the formula, is the k-th prior state prediction value; is the measured predicted value. However, due to the existence of input noise and measurement noise in the system, the estimated values ​​have errors, so Kalman filtering is needed to eliminate the noise.

[0210] Step (3) The final prediction value of the Kalman filter at time k can be expressed as:

[0211]

[0212] In the formula, K k is the Kalman gain, and K k ∈[0,H -1 ]. It can be seen from formula (28) that when K k = 0, That is, the final prediction value is equal to the prior state prediction value. k =H -1 hour, That is, the final predicted value is equal to the measured predicted value. Therefore, it is necessary to find the optimal Kalman gain so that the final predicted value is infinitely close to the true value x k .

[0213] Take the prediction error e k and prior error for:

[0214]

[0215]

[0216] When the variance of the error is minimized, the final predicted value is closest to the true value. The covariance matrix of the error is P k for:

[0217]

[0218] Combining equation (29) and equation (31), we can get:

[0219]

[0220] For the convenience of calculation, let

[0221]

[0222] in

[0223]

[0224] Combining formula (25), we can get:

[0225]

[0226] Step (4) Calculate the error covariance matrix P k :

[0227]

[0228] The covariance matrix P of the error k The trace is:

[0229]

[0230] Find a suitable K k So that tr(P k ) is minimum, the prediction error e k The variance of tr(P k ) is derived, we have:

[0231]

[0232] Pick

[0233]

[0234] Then there is

[0235]

[0236] It is worth noting that in the practical application of Kalman filtering, whether the input noise covariance matrix Q and the measurement noise covariance matrix R can be correctly estimated has a great influence on the calculation accuracy and convergence speed of Kalman filtering. However, in practical engineering applications, since the true values ​​are unknown, the values ​​of Q and R are usually empirical values, which will increase the error and reduce the convergence speed.

[0237] However, when calculating harmonic impedance, in order to avoid errors caused by changes in the harmonic impedance on the system side due to long samples, the PCC point harmonic voltage and current data in a short period of time are usually selected. At this time, it can be considered that the harmonic impedance on the system side is constant in a short period of time. Based on this, in step (2), the present invention proposes a hyperparameter estimation method for the covariance matrix of input noise and measurement noise to determine the input noise covariance matrix Q and the covariance matrix R of measurement noise, as follows:

[0238] Assume that the calculation error e is:

[0239] e=Z s-j -Z s-z (41)

[0240] In the formula, Z S-j is the calculated value of harmonic impedance on the system side; Z S-z is the true value of the harmonic impedance on the system side. The variance of the error e is:

[0241] VAR(e)=VAR(Z s-j -Z s-z )(42)

[0242] Where VAR() represents the calculation of variance. In the process of harmonic impedance calculation, Z S-z is constant and can be regarded as a constant. According to the variance operation theorem:

[0243] VAR(e)=VAR(Z s-j -Z s-z )=VAR(Z s-j )(43)

[0244] It can be seen from formula (43) that the values ​​of the input noise covariance matrix Q and the measurement noise covariance matrix R can be determined by the variance of the results calculated by FastICA and the fluctuation method.

[0245] In step (5), the error covariance matrix P is updatedk , the formula is:

[0246]

[0247] In step (6), let k = k + 1, repeat step (3), reach the maximum number of iterations 100 times, and output the final prediction value, such as Figure 4 shown.

[0248] The updated error covariance matrix is ​​used to filter the signal at time k+1. In subsequent iterative calculations, the predicted value will gradually approach the true value. The smaller the error between the calculated result and the true value, the smaller the error generated by the Kalman filter.

[0249] Based on the same inventive concept, this embodiment also provides a harmonic impedance calculation system based on data fusion, including:

[0250] A preliminary calculation module for obtaining the harmonic impedance of the measurement signal using the FastICA algorithm and the wave quantity method respectively;

[0251] Parameter calculation module, used to obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method;

[0252] Correction module for calculating the prior error covariance matrix

[0253] Gain calculation module, used to calculate the Kalman gain K k ;

[0254] Update module, used to update the error covariance matrix P k ;

[0255] The prediction output module is used to return the iterative calculation of the parameter calculation module, stop the iteration when the preset iteration value is reached, and output the final prediction value.

[0256] Based on the same inventive concept, this embodiment also provides a computing device, including: one or more processors, one or more memories, and one or more programs, wherein the programs are stored in the memories and configured to be executed by the processors, and when the programs are loaded into the processors, the steps of the harmonic impedance calculation method based on data fusion according to any one of the above items are implemented.

Claims

1. A harmonic impedance calculation method based on data fusion, characterized in that: include: (1) The FastICA algorithm and the wave quantity method are used to obtain the harmonic impedance of the measurement signal; (2) Obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method; (3) Calculate the prior error covariance matrix (4) Calculate the Kalman gain K k ; (5) Update the error covariance matrix P k ; (6) Return to step (3) to iterate the calculation of the prior error covariance matrix. When the preset iteration value is reached, stop the iteration and output the final prediction value.

2. The harmonic impedance calculation method based on data fusion according to claim 1 is characterized in that: The method of using the FastICA algorithm to obtain the harmonic impedance of the measurement signal in step (1) includes: Using entropy to represent the non-Gaussian measure, negative entropy is defined as: N g (Y)=H(Y G )-H(Y) (1) in H(Y)=-∫P Y (ξ) lg P Y (ξ)dξ (2) Where Y G and Y are Gaussian random variables with equal variances, P Y (ξ) is the probability density function of the random variable Y, ξ represents the value of the random variable Y; N g (Y) is a measure of non-Gaussianity; when independent component separation is completed, N g (Y) has the largest value; the approximate calculation formula for negative entropy is: N g (Y)={E[g(Y)]-E[g(Y G )]} 2 (3) in g=tan(a1y),1≤a1≤2 (4) In the formula, E is the average value, y is the sample value; g is the nonlinear function; The measurement signal is preprocessed, and the preprocessing includes decentralization and whitening processing; the decentralization is expressed as: Among them, x(t i ) represents the original signal at time point t i The value of represents the signal value after decentralized processing, and N represents the total number of signal samples; The mathematical expression of the whitening process is: in In the formula, Z(t i ) is the signal after whitening processing, matrix R is the eigenvector matrix of the covariance matrix, matrix V is the whitening matrix, matrix U is the eigenvalue matrix of the covariance matrix, and the measurement signal is whitened to obtain matrix V, so that the measurement signal satisfies E[ZZ T ]=I, I is the identity matrix, T represents transpose; Perform iterative calculations to maximize the negative entropy; the maximum approximate value of the negative entropy must meet the following conditions: R[g(V T Z)]+βV=0 (8) in β=R[V1 T Zg(V1 T Z)] (9) In the formula, V1 is the optimized V value; the Newton iteration method is used to solve the above formula to obtain the optimal solution, and the mathematical expression is: In the formula, i is the number of iterations, V i represents the value of V at the i-th iteration, ‖‖ represents the norm, and g′ is the derivative of g; when ‖V i+1 ‖When it tends to be constant, the separation of independent components is completed; The harmonic analysis Norton equivalent equation is converted into complex form, as shown in equation (11); In the formula, ΔU pcc-x and ΔI pcc-x are the real parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common busbar; ΔU pcc-y and ΔI pcc-y are the imaginary parts of the harmonic voltage at the PCC point and the harmonic current vector fluctuation at the common bus; ΔI s-x and ΔI c-x are the real parts of the harmonic source current vector fluctuations on the system side and the user side respectively; ΔI s-y and ΔI c-y are the imaginary parts of the harmonic source current vector fluctuations on the system side and the user side respectively; a mn (m,n=1,2,3,4) is the harmonic source current fluctuation coefficient; Convert equation (11) into matrix form, expressed as: X=AS (12) In the formula, X = [ΔU pcc-x , ΔU pcc-y , ΔI pcc-x , ΔI pcc-y ] T ; S = [ΔI c-x , ΔI c-y , ΔI s-x , ΔI s-y ] T ; A is the coefficient matrix; According to formula (12), after decoupling, four independent components I are obtained: n shape information, n = 1, 2, 3, 4; set up in is a non-zero real number, n = 1, 2, 3, 4; use the least squares method to find I n The mixing coefficient K m =[k 1m ,k 2m ,k 3m ,k 4m ], mathematically expressed as: K m =(I T I) -1 I T X m (13) Among them, X m is the mth vector in the X matrix, k mn The relationship with the harmonic impedance on the system side and the user side is as follows: The solution of formula (15) is as follows: In formula (16), Z s-x is the real part of the harmonic impedance on the system side; Z s-y is the imaginary part of the harmonic impedance on the system side; Z c-x is the real part of the user-side harmonic impedance; Z c-y is the imaginary part of the harmonic impedance at the user side; By judging the positive or negative of the real part of the harmonic impedance, the solution on the system side is inferred. When j = 1, solve the coefficient equation and obtain the solution corresponding to I1, which is set as b 11 and b 21 , when j = 2, solve the coefficient equation and obtain the corresponding solution of I2, set it as b 12 and b 22 ; Then the harmonic impedance on the system side is expressed as:

3. The harmonic impedance calculation method based on data fusion according to claim 1 is characterized in that: The step (1) of obtaining the harmonic impedance of the measurement signal by using the wave quantity method comprises: Considering the changes of harmonic current on the system side and harmonic current on the user side, the mathematical expression is: is the voltage fluctuation at the PCC point, is the current fluctuation at the PCC point, The fluctuation amount of current injected into the user side, The fluctuation of the injected current is measured by the system, Z s represents the system impedance, Z c Indicates user impedance; Definition of PCC point harmonic fluctuation volt-ampere characteristic K hVA It is the ratio of the harmonic voltage fluctuation at the PCC point to the harmonic current fluctuation at the common bus: Analyze the PCC point fluctuations caused by fluctuations on the system side or user side: K hVA The slope of the harmonic volt-ampere characteristic curve is corresponding to the harmonic impedance on both sides of the PPC point estimated by the ratio of the harmonic current and harmonic voltage fluctuations at the PCC point; when the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is positive, the result is the harmonic impedance on the system side, and when the real part of the harmonic fluctuation volt-ampere characteristic at the PCC point is negative, the result is the negative number of the harmonic impedance on the user side; Since the load capacity on the user side is much smaller than the short-circuit capacity on the system side, there is Z s < <Z c , transform formula (19) into: In the formula, e sys is the error caused by background harmonic fluctuation.

4. The harmonic impedance calculation method based on data fusion according to claim 1, characterized in that: Step (2) of obtaining the prior state prediction value of the measurement signal and measured predicted values include: Consider the state space model of the linear discrete dynamic system as: X(k)=AX(k-1)+ΓW(k-1) (23) Z(k)=HX(k)+V(k) (24) Wherein, k is the discrete time; X(k) and Z(k) are the state variable and the measurement variable respectively; W(k-1) and V(k) are the input noise and the measurement noise respectively; Γ is the noise driving matrix; A is the state transfer matrix; H is the observation matrix; Equation (23) is the state equation, and Equation (24) is the observation equation; Assume that the input noise and measurement noise are independent of each other, their mean is zero, their covariance matrices are Q and R respectively, and their statistical characteristics are: Where k is the discrete time, j is the discrete time that is not synchronized with k, W(k) and V(k) are the input noise and measurement noise, Q and R are the covariance matrix of the input noise and the covariance matrix of the measurement noise, P(W) and P(V) represent the probability distribution of the input noise and the output noise, δ kj is the Kronecker symbol, which is 1 when k = j, indicating that the noise at the same time has non-zero correlation and the noise at different times is independent of each other; Considering X(k) and Z(k) as the calculation results of the two methods respectively, the Kalman filter recursive formula is as follows: When the noise error is not considered, the prior prediction value is obtained, which is mathematically expressed as: In the formula, is the k-th prior state prediction value; To measure the predicted value.

5. The harmonic impedance calculation method based on data fusion according to claim 4 is characterized in that: The step (3) comprises: The noise is eliminated by Kalman filtering, and the final prediction value of Kalman filtering at time k is expressed as: In the formula, K k is the Kalman gain, and K k ∈[0,H -1 ]; According to formula (28), when K k = 0, The final prediction value is equal to the prior state prediction value. k =H -1 hour, The final predicted value is equal to the measured predicted value; Take the prediction error e k and prior error for: When the variance of the error is minimized, the final predicted value is closest to the true value x k , the error covariance matrix P k for: Combining equation (29) and equation (31), we get: make: in Combining formula (25) to obtain the prior error covariance matrix:

6. The harmonic impedance calculation method based on data fusion according to claim 5, characterized in that: The step (4) comprises: Calculate the error covariance matrix P k : The covariance matrix P of the error k The trace is: When tr(P k ) is minimum, the prediction error e k The variance of tr(P k ) is derived, we have: Pick Then there is 7. The harmonic impedance calculation method based on data fusion according to claim 6, characterized in that: The step (2) of obtaining the input noise covariance matrix Q and the measurement noise covariance matrix R by the hyperparameter estimation method includes: Assume that the calculation error e is: e=Z s-j -Z s-z (41) In the formula, Z s-j is the calculated value of harmonic impedance on the system side; Z s-z is the true value of the harmonic impedance on the system side; the variance of the error e is: VAR(e)=VAR(Z s-j -Z s-z ) (42) Where VAR() represents the calculation of variance; in the process of harmonic impedance calculation, Z s-z is a constant; according to the variance operation theorem: VAR(e)=VAR(Z s-j -Z s-z )=VAR(Z s-j ) (43) According to formula (43), the values ​​of the input noise covariance matrix Q and the measurement noise covariance matrix R are determined by the variance of the results calculated by FastICA and the fluctuation method.

8. The harmonic impedance calculation method based on data fusion according to claim 7 is characterized in that: The step (5) comprises: Update the error covariance matrix P k , the formula is:

9. A harmonic impedance calculation system based on data fusion, characterized in that: include: A preliminary calculation module for obtaining the harmonic impedance of the measurement signal using the FastICA algorithm and the wave quantity method respectively; Parameter calculation module, used to obtain the prior state prediction value of the measurement signal and measured predicted values The input noise covariance matrix Q and the measurement noise covariance matrix R are obtained by the hyperparameter estimation method; Correction module for calculating the prior error covariance matrix Gain calculation module, used to calculate the Kalman gain K k ; Update module, used to update the error covariance matrix P k ; The prediction output module is used to return the iterative calculation of the parameter calculation module, stop the iteration when the preset iteration value is reached, and output the final prediction value.

10. A computing device, characterized in that: include: One or more processors, one or more memories, and one or more programs, wherein the programs are stored in the memories and configured to be executed by the processors, and when the programs are loaded into the processors, the steps of the harmonic impedance calculation method based on data fusion according to any one of claims 1 to 8 are implemented.

Images