An improved dynamic harmonic estimation method and system

Abstract

The present invention provides an improved dynamic harmonic estimation method and system, comprising: collecting harmonic voltages of synchronous phasor measurement devices of each branch or bus of a power grid; substituting the harmonic voltages into a nonlinear state equation and a measurement equation model; using an improved Sage‑Husa square root unscented Kalman filter algorithm, and estimating the harmonic current of the branch to be estimated based on the nonlinear state equation and the measurement equation model, to obtain an estimated value of the harmonic current of the branch to be estimated; wherein the improved Sage‑Husa square root unscented Kalman filter algorithm includes: updating the state of the nonlinear state equation and the measurement equation model based on measurement update and forgetting factor. The present invention addresses the problem that accurate estimation results cannot be obtained from unknown noise when using square root unscented Kalman filter to estimate the harmonic state of a power system. The present invention introduces the idea of ​​maximum a posteriori estimation based on the Sage‑Husa filter algorithm and utilizes the forgetting factor of the Sage‑Husa filter to achieve real-time dynamic estimation of unknown noise and obtain high-precision harmonic estimation results.

Description

Technical Field

[0001] The present invention belongs to the technical field of power system state estimation, and in particular relates to an improved dynamic harmonic estimation method and system. Background Art

[0002] Harmonic state estimation (HSE) is an effective harmonic monitoring and analysis technology that has developed in recent years. Harmonic state estimation utilizes data provided by synchronized phasor measurement units (PMUs) installed on some of the system's buses and lines to infer the branch harmonic current and node harmonic voltage states of the entire power grid based on appropriate estimation criteria. Dynamic harmonic state estimation algorithms estimate the state quantities at the next moment based on the motion equations of the power system harmonics and the measured data at a certain moment as initial values. Compared to static harmonic state estimation, dynamic harmonic state estimation can perform both filtering and prediction.

[0003] The Kalman filter (KF) is a highly efficient autoregressive filter that can estimate the state of a dynamic system from a combination of information in the presence of numerous uncertainties. The traditional KF algorithm is only suitable for linear system state estimation. To address the process noise and observation noise estimation issues in nonlinear system state estimation, derivative algorithms based on the KF algorithm, such as the Extended Kalman Filter (EKF) and the Unscented Kalman Filter (UKF), have been proposed. However, the EKF algorithm approximates the nonlinear model to a linear model by extracting low-order terms through a Taylor series expansion, which introduces approximation errors. The UKF offers improved estimation accuracy compared to the EKF, but suffers from the drawback of difficulty in determining the positive definiteness of the state noise covariance matrix.

[0004] Considering that the power system is a time-varying nonlinear system and its harmonics are also dynamically changing, the process noise variance and observation noise variance of the system are unknown when estimating harmonics, resulting in strong uncertainty. Estimation can usually only be done through experience, and incorrect parameter estimation often leads to filter divergence, which in turn makes it impossible to obtain accurate estimation results. The Square-Root Unscented Kalman Filter (SRUKF) uses the square root of the covariance instead of the covariance in the recursive operation, which can to some extent solve the problem of filter divergence caused by the negative definiteness of the covariance matrix. However, its essence is to use the normal distribution to approximate the posterior probability density of the system state, which does not effectively solve the problem that the strong uncertainty of power system noise will affect the accuracy of the estimation results. To improve the robustness and accuracy of power system harmonic state estimation, the shortcomings of existing state estimation algorithms should be scientifically improved so that they can filter out unknown noise and obtain accurate harmonic estimation. Summary of the Invention

[0005] To overcome the above-mentioned deficiencies of the prior art, the present invention proposes an improved dynamic harmonic estimation method, the improvement of which includes:

[0006] Collect harmonic voltages of synchronized phasor measurement devices on each branch or bus of the power grid;

[0007] Substituting the harmonic voltage into a pre-established nonlinear state equation and measurement equation model;

[0008] An improved Sage-Husa square root unscented Kalman filter algorithm is used to perform harmonic current estimation calculation on the branch to be estimated based on the nonlinear state equation and the measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated;

[0009] The improved Sage-Husa square root unscented Kalman filter algorithm includes: performing state updates on the nonlinear state equation and the measurement equation model based on measurement updates and forgetting factors.

[0010] The first preferred technical solution provided by the present invention is improved in that the nonlinear state equation and measurement equation model are as follows:

[0011]

[0012] Wherein, subscript k represents the next sampling time, subscript k-1 represents the current sampling time; X k Represents the estimated value of harmonic current at the next sampling moment, X k-1 Indicates the estimated value of harmonic current at the current sampling moment, U krepresents the harmonic voltage at the next sampling moment, F is the system state function, H is the system observation function, A k Represents the system matrix at the next sampling moment, B k Represents the control matrix at the next sampling moment, C k Represents the observation matrix at the next sampling moment, D k represents the direct transfer matrix at the next sampling moment, W k represents the system process noise at the next sampling moment, V k represents the system observation noise at the next sampling moment, Z k Indicates the harmonic voltage observation value at the next sampling moment.

[0013] The second preferred technical solution provided by the present invention is improved in that the improved Sage-Husa square root unscented Kalman filter algorithm is used to perform harmonic current estimation calculation on the preset branch to be estimated according to the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated, including:

[0014] Based on the nonlinear state equation and measurement equation model, an unscented Kalman filter algorithm is used for time update;

[0015] The nonlinear state equation and measurement equation model are updated based on the time update result and the Sage-Husa method; the initial value of the measurement sequence in the nonlinear state equation and measurement equation model is initialized by the harmonic voltage collected at the first sampling moment;

[0016] Based on the measurement update result and the preset forgetting factor, the state of the nonlinear state equation and the measurement equation model is updated to obtain an estimated value of the harmonic current.

[0017] The third preferred technical solution provided by the present invention is improved in that the initial values ​​of the measurement sequence in the nonlinear state equation and measurement equation model are initialized as shown in the following formula:

[0018]

[0019] Wherein, subscript 0 represents the initial sampling time, X0 represents the initial value of harmonic current, Q0 represents the covariance vector of system process noise at the initial sampling time, R0 represents the covariance vector of system observation noise at the initial sampling time, is the square root of Q0, is the square root of R0, is the mathematical expectation of X0; E represents the mathematical expectation, Chol function is the matrix Cholesky factor decomposition function, S0 is The covariance with X0;

[0020] The relationship between the initial value X0 of the harmonic current and the harmonic voltage is as follows:

[0021] Z0=H[X0,U0]=C0X0+D0U0+V0

[0022] Where Z0 represents the harmonic voltage after measurement transformation at the initial sampling time, U0 represents the harmonic voltage at the initial sampling time, V0 represents the system observation noise at the initial sampling time, H is the observation function, C is the observation matrix, and D is the direct transfer matrix.

[0023] The fourth preferred technical solution provided by the present invention is improved in that the time update is performed using an unscented Kalman filter algorithm based on the nonlinear state equation and measurement equation model, including:

[0024] Performing an untraceable transformation on the harmonic current in the initialized nonlinear state equation and measurement equation model to obtain a Sigma matrix containing a Sigma point set;

[0025] The Sigma matrix is ​​nonlinearly transformed according to the nonlinear state equation, and the square root of the covariance of the system process noise is used to participate in the recursive operation of the unscented Kalman filter algorithm to predict the harmonic current and the harmonic current covariance.

[0026] The fifth preferred technical solution provided by the present invention is improved in that the calculation formula of the recursive operation is as follows:

[0027]

[0028] Where, Indicates the harmonic current intermediate value at the next sampling moment, χ i,k-1 Indicates the value of the i-th Sigma point corresponding to the harmonic current at the current sampling moment. The value range of i is 0 to 2n, and n is the number of dimensions of the harmonic current. U k-1 Indicates the harmonic voltage at the current sampling moment, A k-1 Represents the system matrix at the current sampling moment, B k-1 represents the control matrix at the current sampling moment, F is the system state function, Indicates the weight of the i-th Sigma point when calculating the mean, Indicates the weight of the i-th Sigma point when calculating the covariance, χ * i,k-1 Indicates the corresponding χ i,k-1 The intermediate process variable, S xx,k Indicates the covariance of the harmonic current at the next sampling moment, Q k represents the system process noise covariance vector at the next sampling moment, Indicates the weight of the 0th Sigma point in calculating the covariance, χ * 0,k-1 Indicates the intermediate process variable corresponding to the 0th Sigma point corresponding to the harmonic current at the current sampling moment, Indicates the corresponding S xx,k The intermediate process variables, qr is the orthogonal triangular matrix decomposition operation, Cholupdata is the Cholesky decomposition operation;

[0029] The value of the i-th Sigma point corresponding to the harmonic current at the current sampling moment is χ i,k-1 As shown in the following formula:

[0030]

[0031] Where σ k-1 Indicates the traceless change parameter at the current sampling moment, Indicates the mean value of harmonic current at the current sampling moment;

[0032] The calculation formula of the traceless change parameter at the current sampling moment is as follows:

[0033]

[0034] Where S xx,k-1 represents the covariance of the harmonic current at the current sampling moment, and λ is the scaling factor.

[0035] The sixth preferred technical solution provided by the present invention is improved in that the nonlinear state equation and measurement equation model are measured and updated based on the time-updated result and the preset forgetting factor based on the Sage-Husa method, including:

[0036] Based on the result of time update, the Sigma point is resampled to obtain the Sigma point;

[0037] Perform nonlinear transformation on Sigma points through measurement equation and calculate measurement residuals;

[0038] Based on the measurement residual and taking into account the preset forgetting factor, the square root of the system observation noise covariance is updated using the square root unscented Kalman filtering method based on Sage-Husa;

[0039] Based on the square root of the system observation noise covariance, the covariance between the harmonic voltage and harmonic current and the covariance of the harmonic voltage are updated;

[0040] The Kalman filter gain is calculated based on the covariance between the harmonic voltages and harmonic currents and the covariance of the harmonic voltages.

[0041] The seventh preferred technical solution provided by the present invention is improved in that the resampled Sigma point set is nonlinearly transformed and the measurement residual is calculated by the measurement equation, as shown in the following formula:

[0042]

[0043] Where H is the observation function, χ i,k Indicates the value of the i-th Sigma point corresponding to the harmonic current at the next sampling moment of resampling, U k-1 Indicates the harmonic voltage at the current sampling moment, The weight of the i-th Sigma point when calculating the mean, It represents the harmonic voltage after adding weights, correction and measurement transformation at the next sampling moment, Z k Indicates the harmonic voltage observation value collected at the next sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, Z * i,k Indicates correspondence The intermediate process variable of the i-th Sigma point.

[0044] The eighth preferred technical solution provided by the present invention is improved in that, based on the measurement residual and taking into account the preset forgetting factor, the square root of the system observation noise covariance is updated using the square root unscented Kalman filtering method based on Sage-Husa, as shown in the following formula:

[0045]

[0046] Where, Represents the square root of the system observation noise covariance at the next sampling moment, R ** k Indicates correspondence The first intermediate process variable, R * k Indicates correspondence The second intermediate process variable, represents the square root of the system observation noise covariance at the current sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, d k is the forgetting correction factor at the next sampling moment calculated based on the forgetting factor, Represents the weight of the i-th Sigma point when calculating the covariance, It represents the harmonic voltage after the weighted correction and measurement transformation of the i-th Sigma point at the next sampling moment. It represents the harmonic voltage of the i-th Sigma point after adding weights, correction and measurement transformation at the current sampling moment. Represents R * k The transpose of , Cholupdata is the Cholesky decomposition operation, and diag is the operation of constructing a diagonal matrix;

[0047] The forgetting correction factor d at the next sampling moment k The calculation formula is as follows:

[0048] d k =(1-b) / (1-b k+1 )

[0049] Where b is the forgetting factor.

[0050] The ninth preferred technical solution provided by the present invention is improved in that the covariance between the harmonic voltage and the harmonic current and the covariance of the harmonic voltage are updated based on the square root of the system observation noise covariance, as shown in the following formula:

[0051]

[0052] Where, P xz,k It represents the covariance between harmonic voltage and harmonic current at the next sampling moment, S zz,k Represents the covariance of the harmonic voltage at the next sampling moment, S * zz,k Indicates the corresponding S zz,k The intermediate process variables, It represents the harmonic voltage after adding weights, correction and measurement transformation at the next sampling moment, Z * i,k Indicates correspondence The intermediate process variable of the i-th Sigma point, Indicates the intermediate value of harmonic current obtained at the next sampling moment, χ * i,k Indicates the corresponding χ i,k The intermediate process variable, χ i,k Indicates the value of the i-th Sigma point corresponding to the harmonic current at the next sampling moment, Indicates the weight of the i-th Sigma point when calculating the covariance, R * k Indicates correspondence The second intermediate process variable, Indicates the weight of the 0th Sigma point when calculating the covariance, Z * 0,k-1 Indicates correspondence The intermediate process variable of the 0th Sigma point, It represents the harmonic voltage after weighted correction and measurement transformation at the current sampling moment, qr is the orthogonal triangular matrix decomposition operation, and Cholupdata is the Cholesky decomposition operation.

[0053] The tenth preferred technical solution provided by the present invention is improved in that the calculation formula of the Kalman filter gain is as follows:

[0054]

[0055] Where K k Represents the Kalman filter gain at the next sampling moment, S T zz,k Indicates S zz,k The transpose of .

[0056] The eleventh preferred technical solution provided by the present invention is improved in that, based on the measurement update result, the state of the nonlinear state equation and the measurement equation model is updated to obtain the estimated value of the harmonic current, including:

[0057] Calculation is performed based on the harmonic current intermediate value, measurement residual and Kalman filter gain in the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current.

[0058] The twelfth preferred technical solution provided by the present invention is improved in that the calculation formula of the estimated value of the harmonic current is as follows:

[0059]

[0060] Among them, X k Represents the estimated value of the harmonic current at the next sampling moment, Indicates the harmonic current intermediate value obtained at the next sampling moment, K k represents the Kalman filter gain at the next sampling moment, e k Indicates the measurement residual of the harmonic voltage at the next sampling moment.

[0061] The thirteenth preferred technical solution provided by the present invention is improved in that, after obtaining the estimated value of the harmonic current, the method further includes:

[0062] updating the covariance estimate of the harmonic current and updating the square root of the system process noise covariance vector;

[0063] The parameters of the nonlinear state equation and measurement equation model at the current sampling moment are transferred to the nonlinear state equation and measurement equation model at the next sampling moment.

[0064] The fourteenth preferred technical solution provided by the present invention is improved in that the calculation formula of the covariance estimation value of the harmonic current is as follows:

[0065] S k =cholupdate{S k-1 ,K k S zz,k ,-1}

[0066] Among them, S k It represents the covariance estimate of the harmonic current at the next sampling moment, S k-1 Represents the covariance estimate of the harmonic current at the current sampling moment, K k Represents the Kalman filter gain at the next sampling moment, S zz,k Represents the covariance of the harmonic voltage at the next sampling moment, and Cholupdata is the Cholesky decomposition operation;

[0067] The square root of the system process noise covariance vector is calculated as follows:

[0068]

[0069] in, represents the square root of the system process noise covariance vector at the next sampling moment, Q k-1 represents the system process noise covariance vector after the k-1th recursion, d k is the forgetting correction factor at the next sampling moment, K k represents the Kalman filter gain at the next sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, S zz,k represents the covariance of the harmonic voltage after the kth recursion, Q ** k Indicates correspondence The first intermediate process variable, Q * k Indicates correspondence The second intermediate process variable, Cholupdata is the Cholesky decomposition operation, and diag is the diagonal matrix construction operation.

[0070] The fifteenth preferred technical solution provided by the present invention is improved in that the method of transferring the parameters of the nonlinear state equation and measurement equation model at the current sampling moment to the nonlinear state equation and measurement equation model at the next sampling moment includes:

[0071] The intermediate value of the harmonic current at the current sampling moment is used as the intermediate value of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, the covariance estimate of the harmonic current at the current sampling moment is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, the square root of the system process noise covariance vector at the current sampling moment is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, and the square root of the system observation noise covariance at the current sampling moment is used as the square root of the system observation noise covariance of the nonlinear state equation and measurement equation model at the next sampling moment.

[0072] Based on the same inventive concept, the present invention also provides an improved dynamic harmonic estimation system, which is improved in that it includes: a data acquisition module, a data input module and a harmonic estimation module;

[0073] The data acquisition module is used to collect the harmonic voltages of the synchronized phasor measurement devices of each branch or bus of the power grid;

[0074] The data input module is used to bring the harmonic voltage into a pre-established nonlinear state equation and measurement equation model;

[0075] The harmonic estimation module is used to use an improved Sage-Husa square root unscented Kalman filter algorithm to perform harmonic current estimation calculation on the branch to be estimated based on the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated;

[0076] The improved Sage-Husa square root unscented Kalman filter algorithm includes: performing state updates on the nonlinear state equation and the measurement equation model based on measurement updates and forgetting factors.

[0077] The sixteenth preferred technical solution provided by the present invention is improved in that the harmonic estimation module includes: a time update module, a measurement update module and a state update module;

[0078] The time update module is used to perform time update based on the nonlinear state equation and measurement equation model using an unscented Kalman filter algorithm; the initial value of the measurement sequence in the nonlinear state equation and measurement equation model is initialized by the harmonic voltage collected at the first sampling moment;

[0079] The measurement update module is used to perform measurement update on the nonlinear state equation and measurement equation model based on the time update result and the Sage-Husa method;

[0080] The state updating module is used to update the state of the nonlinear state equation and the measurement equation model based on the measurement update result and the preset forgetting factor to obtain the estimated value of the harmonic current.

[0081] Compared with the closest prior art, the present invention has the following beneficial effects:

[0082] The present invention provides an improved dynamic harmonic estimation method and system, comprising: collecting harmonic voltages from synchronized phasor measurement devices on each branch or bus of a power grid; subjecting the harmonic voltages to a pre-established nonlinear state equation and measurement equation model; and employing an improved Sage-Husa square root unscented Kalman filter algorithm to perform harmonic current estimation calculations on the branch to be estimated based on the nonlinear state equation and measurement equation model, thereby obtaining an estimated value of the harmonic current of the branch to be estimated. The improved Sage-Husa square root unscented Kalman filter algorithm includes performing a state update on the nonlinear state equation and measurement equation model based on measurement update and a forgetting factor. To address the problem of the inability to obtain accurate estimation results from unknown noise when using a square root unscented Kalman filter for power system harmonic state estimation, the present invention introduces the concept of maximum a posteriori estimation in the Sage-Husa filter algorithm to propose an improved Sage-Husa square root unscented Kalman filter harmonic dynamic estimation algorithm. By utilizing the forgetting factor in the Sage-Husa filter, the algorithm can achieve real-time dynamic estimation of unknown noise and obtain high-precision harmonic estimation results.

[0083] The improved dynamic harmonic estimation method and system provided by the present invention can further continuously correct the statistical characteristics of process noise and observation noise through recursive filtering, thereby obtaining more accurate harmonic estimation results and realizing dynamic tracking of harmonic states. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 A schematic flow chart of an improved dynamic harmonic estimation method provided by the present invention;

[0085] Figure 2 A simplified flowchart of an improved dynamic harmonic estimation method provided by the present invention;

[0086] Figure 3 A general flow chart of state estimation using an improved Sage-Husa square root unscented Kalman filter provided by the present invention;

[0087] Figure 4 A detailed flow chart of an improved dynamic harmonic estimation method provided by the present invention;

[0088] Figure 5 A schematic diagram of the basic structure of an improved dynamic harmonic estimation system provided by the present invention;

[0089] Figure 6 A detailed structural diagram of an improved dynamic harmonic estimation system provided by the present invention. DETAILED DESCRIPTION

[0090] The specific embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0091] Example 1:

[0092] The flow chart of an improved dynamic harmonic estimation method provided by the present invention is as follows: Figure 1 As shown, including:

[0093] Step 1: Collect the harmonic voltages of the synchronized phasor measurement devices of each branch or bus of the power grid;

[0094] Step 2: Substitute the harmonic voltage into the pre-established nonlinear state equation and measurement equation model;

[0095] Step 3: Using the improved Sage-Husa square root unscented Kalman filter algorithm, based on the nonlinear state equation and measurement equation model, the harmonic current of the estimated branch is estimated and the estimated value of the harmonic current of the estimated branch is obtained;

[0096] Among them, the improved Sage-Husa square root unscented Kalman filter algorithm includes: updating the state of the nonlinear state equation and measurement equation model based on measurement update and forgetting factor.

[0097] Specifically, the improved dynamic harmonic estimation method provided in this embodiment is a dynamic harmonic estimation method of an improved Sage-Husa square root unscented Kalman filter, and its implementation steps are as follows: Figure 2 As shown, including:

[0098] Step 101: Establish a nonlinear state equation and measurement equation model based on the power grid topology;

[0099] Step 102: Input the current measurement data;

[0100] Step 103: Using the improved Sage-Husa square root unscented Kalman filter algorithm to perform dynamic harmonic estimation calculation;

[0101] Step 104: Parameter transfer, repeat steps 102 to 104.

[0102] In the present invention, a general flow chart of an improved Sage-Husa square root unscented Kalman filter to achieve state estimation is as follows: Figure 3 As shown, the method is applied to the specific dynamic harmonic estimation to provide an improved dynamic harmonic estimation method. The detailed process is as follows Figure 4 shown.

[0103] In the improved Sage-Husa square root unscented Kalman filter dynamic harmonic estimation method, the measurement equation in step 101 is to establish a nonlinear harmonic measurement equation. According to the power system network architecture, the measurement quantity is set to the bus harmonic voltage, and the state estimation quantity is the harmonic current injected into the node to be estimated. The general form of the state equation and the measurement equation is:

[0104]

[0105] Where, I is the harmonic current state estimation sequence;

[0106] Y is the admittance matrix function;

[0107] I a Inject vectors for additional harmonic currents, typically transformer magnetizing currents;

[0108] Z is the harmonic voltage measurement vector after measurement transformation;

[0109] U is the harmonic voltage measurement vector;

[0110] M is the harmonic voltage measurement matrix function;

[0111] ε and η are the process noise vector and observation noise vector respectively;

[0112] The subscripts k and k-1 of the distribution represent the k-th and k-1-th sampling moments, or the k-th and k-1-th recursions.

[0113] The improved Sage-Husa square root unscented Kalman filter dynamic harmonic estimation method mentioned above requires that the measurement equation (1) be changed to the discrete system standard model shown in equation (15) before applying the improved Sage-Husa square root unscented Kalman filter for dynamic harmonic estimation:

[0114]

[0115] Where X and Z are the n-dimensional state variables and m-dimensional measurement variables, respectively. F and H are the system state function and observation function, respectively. Since the power system is nonlinear, F and H are nonlinear functions. W and V are the n-dimensional system process noise and m-dimensional system observation noise sequences, respectively. A is the system matrix, B is the control matrix, C is the observation matrix, and D is the direct transfer matrix.

[0116] This embodiment provides a dynamic harmonic estimation method based on an improved Sage-Husa square root unscented Kalman filter. After the power system network architecture, branch transmission line parameters, and load are determined, the state variable sequence X in the harmonic current to be estimated is selected to form the formula (15). k; Select the appropriate bus voltage to form the measurement sequence Z in formula (15) k ; Function Y is Z k Map to X k The configured admittance matrix must satisfy the requirement that its rank is full rank and rank(Y)=n. If not, the measurement sequence needs to be adjusted until the full rank requirement is met. Where n is the number of dimensions of the state vector. In this embodiment, the state variable sequence X k That is the wave current state estimation sequence I k .

[0117] A dynamic harmonic estimation method based on an improved Sage-Husa square root unscented Kalman filter is proposed. The measurement data comes from PMUs deployed on different branches or buses in the power system network. Typically, PMUs are not deployed on all branches or buses in a power system network. Therefore, it is necessary to establish system measurement equations to estimate the harmonic state of branches or buses without PMUs.

[0118] In this embodiment, a dynamic harmonic estimation method of an improved Sage-Husa square root unscented Kalman filter is implemented as follows in step 102:

[0119] a. Initialization.

[0120]

[0121] Where, subscript 0 represents the initial sampling time, X0 represents the initial value of the state vector, i.e., the harmonic current, Q0 represents the covariance vector of the system process noise W0 at the initial sampling time, and R0 represents the covariance vector of the system observation noise V0 at the initial sampling time. is the square root of Q0, is the square root of R0, is the mathematical expectation of X0; E represents the mathematical expectation, Chol function is the matrix Cholesky factor decomposition function, S0 is and the covariance of X0.

[0122] Initialization is performed only at the initial sampling time when k=0. When k>0, no initialization is performed.

[0123] b. Time update: Iterate for k=1, 2, ...

[0124] First, perform an untraceable transformation, which includes two parts: Sigma calculation and its weight calculation. Calculate the value of the Sigma sampling point: Let the sampling point dimension be the dimension n of the state quantity. 2n+1 Sigma sampling points need to be calculated. The calculation method is shown in Equation (3).

[0125]

[0126] Where, χ i represents the value of the i-th Sigma point, and S are the mean and covariance of the n-dimensional state vector X, respectively. λ is the scaling factor, and λ = α 2 (n+ρ)-n. α controls the distribution of Sigma points, and its value range is (10 -4 ,1). The premise of the value of ρ scale parameter is to ensure that (n+λ)S is a semi-positive definite matrix.

[0127] Calculate the weight ω of the Sigma point:

[0128]

[0129] Where the superscripts m and c denote the weights used in the mean and covariance calculations, respectively, and the subscript i denotes the dimension. The parameter β is ≥ 0, and its value depends on whether the influence of higher-order terms needs to be considered.

[0130] After the unscented transformation, the matrix containing the Sigma point set is obtained:

[0131]

[0132] Among them, σ k-1 represents the traceless change parameter after the k-1th recursion, Where S xx,k-1 represents the covariance of the harmonic current after the k-1th recursion.

[0133] According to the system formula (15), the Sigma matrix is ​​transformed nonlinearly, and the state quantity and covariance are predicted in one step. The square root of the system process noise covariance is taken to replace the covariance in the ordinary UKF algorithm to participate in the recursive operation.

[0134]

[0135] Where, represents the intermediate value of harmonic current obtained after the kth recursion, χ i,k-1 It represents the value of the i-th Sigma point corresponding to the harmonic current after the k-1th recursion. The value range of i is 0 to 2n, n is the dimension of the harmonic current, U k-1 Represents the harmonic voltage after the k-1th recursion, A k-1 represents the system matrix after the k-1th recursion, B k-1 represents the control matrix after the k-1th recursion, F is the system state function, Indicates the weight of the i-th Sigma point when calculating the mean, Indicates the weight of the i-th Sigma point when calculating the covariance, χ * i,k-1 Indicates the corresponding χ i,k-1 The intermediate process variable, S xx,k represents the covariance of the harmonic current after the kth recursion, Q k represents the system process noise covariance vector after the k-th recursion, Indicates the weight of the 0th Sigma point in calculating the covariance, χ * 0，k-1 It represents the intermediate process variable corresponding to the 0th Sigma point of the harmonic current after the k-1th recursion, Indicates the corresponding S xx,k The intermediate process variables, qr is the orthogonal triangular matrix decomposition operation, and Cholupdata is the Cholesky decomposition operation.

[0136] c. Measurement update.

[0137] First, the k-th sampling period, that is, the k-th recursive Sigma point is resampled.

[0138]

[0139] in,

[0140] Perform nonlinear transformation on the Sigma point set through the measurement equation and calculate the measurement residual e k :

[0141]

[0142] Where H is the observation function, χ i,k Indicates the value of the i-th Sigma point corresponding to the harmonic current after the k-th recursion, U k-1 represents the harmonic voltage after the k-1th recursion, The weight of the i-th Sigma point when calculating the mean, It represents the harmonic voltage after adding weights and correction and measurement transformation after the kth recursion, Z k represents the harmonic voltage observation value collected after the kth recursion, e k Represents the measurement residual of harmonic voltage after the kth recursion, Z * i,k Indicates correspondence The intermediate process variable of the i-th Sigma point.

[0143] The Sage-Husa filtering idea is introduced, and the forgetting factor is added when updating the estimated measurement noise statistical characteristics to improve the SRUKF. The square root of the system observation noise covariance is Updated to

[0144]

[0145] Where, represents the square root of the system observation noise covariance after the kth recursion, R ** k Indicates correspondence The first intermediate process variable, R * k Indicates correspondence The second intermediate process variable, represents the square root of the system observation noise covariance after the k-1th recursion, e k represents the measurement residual of harmonic voltage after the kth recursion, d k is the forgetting correction factor after the k-th recursion, Represents the weight of the i-th Sigma point when calculating the covariance, It represents the harmonic voltage after the k-th recursion, the correction of the i-th Sigma point after adding weights and after measurement and transformation, It represents the harmonic voltage of the i-th Sigma point after the k-1th recursion, after adding weights and correction, and after measurement transformation. Represents R * k The transpose of , Cholupdata is the Cholesky decomposition operation, diag is the construction of diagonal matrix operation; d k =(1-b) / (1-b k+1 ), b is the forgetting factor, which ranges from 0.9 to 1.

[0146] The covariance P between the update quantity measurement and the state quantity xz,k and the measurement variance matrix S zz,k :

[0147]

[0148] Where, P xz,k It represents the covariance between harmonic voltage and harmonic current after the kth recursion, S zz,k represents the covariance of harmonic voltage after the kth recursion, S * zz,k Indicates the corresponding S zz,k The intermediate process variables, It represents the harmonic voltage after adding weights and correction and measurement transformation after the kth recursion, Z * i,kIndicates correspondence The intermediate process variable of the i-th Sigma point, represents the intermediate value of harmonic current obtained after the kth recursion, χ * i,k Indicates the corresponding χ i,k The intermediate process variable, χ i,k It represents the value of the i-th Sigma point corresponding to the harmonic current after the k-th recursion, Indicates the weight of the i-th Sigma point when calculating the covariance, R * k Indicates correspondence The second intermediate process variable, Indicates the weight of the 0th Sigma point when calculating the covariance, Z * 0,k-1 Indicates correspondence The intermediate process variable of the 0th Sigma point, represents the harmonic voltage after the k-1th recursion, after weight addition and correction, and after measurement transformation. qr is the orthogonal triangular matrix decomposition operation, and Cholupdata is the Cholesky decomposition operation.

[0149] Calculate the Kalman filter gain K k :

[0150]

[0151] Where K k represents the Kalman filter gain after the kth recursion, S T zz,k Indicates S zz,k The transpose of .

[0152] d. Status update.

[0153] Estimate the corrected state value X k and the covariance estimate S of the state value k :

[0154]

[0155] Among them, X k represents the estimated value of harmonic current after the kth recursion, It represents the intermediate value of harmonic current obtained after the k-th recursion, S k represents the covariance estimate of the harmonic current after the kth recursion, S k-1 It represents the estimated value of the covariance of the harmonic current after the k-1th recursion.

[0156] Noise Statistical Characteristics of the New Estimation Process

[0157]

[0158] in, represents the square root of the system process noise covariance vector after the kth recursion, Q k-1 represents the system process noise covariance vector after the k-1th recursion, d k is the forgetting correction factor after the kth recursion, K k represents the Kalman filter gain after the kth recursion, e k Represents the measurement residual of harmonic voltage after the kth recursion, S zz,k represents the covariance of the harmonic voltage after the kth recursion, Q ** k Indicates correspondence The first intermediate process variable, Q * k Indicates correspondence The second intermediate process variable.

[0159] In step 104, the parameter transfer is shown in formula (14):

[0160]

[0161] That is, after the arrival of the new moment k+1, the intermediate value of the harmonic current at the previous moment k is used as the intermediate value of the harmonic current of the nonlinear state equation and measurement equation model at the current moment k+1, the covariance estimate of the harmonic current at the previous moment k is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the current moment k+1, the square root of the system process noise covariance vector at the previous moment k is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the current moment k+1, and the square root of the system observation noise covariance at the previous moment k is used as the square root of the system observation noise covariance of the nonlinear state equation and measurement equation model at the current moment k+1.

[0162] Example 2:

[0163] Based on the same inventive concept, the present invention also provides an improved dynamic harmonic estimation system. Since the principles of these devices in solving technical problems are similar to those of the improved dynamic harmonic estimation method, the repeated parts will not be repeated.

[0164] The basic structure of the system is as follows Figure 5 As shown, it includes: a data acquisition module, a data input module and a harmonic estimation module;

[0165] A data acquisition module for collecting harmonic voltages from synchronized phasor measurement devices on each branch or bus of the power grid;

[0166] Data input module, used to bring harmonic voltage into pre-established nonlinear state equation and measurement equation models;

[0167] The harmonic estimation module is used to use the improved Sage-Husa square root unscented Kalman filter algorithm to perform harmonic current estimation calculation on the estimated branch according to the nonlinear state equation and measurement equation model to obtain the estimated value of the harmonic current of the estimated branch;

[0168] Among them, the improved Sage-Husa square root unscented Kalman filter algorithm includes: updating the state of the nonlinear state equation and measurement equation model based on measurement update and forgetting factor.

[0169] The detailed structure of an improved dynamic harmonic estimation system is as follows: Figure 6 shown.

[0170] Among them, the harmonic estimation module includes: a time update module, a measurement update module and a state update module;

[0171] A time update module is used to perform time update based on a nonlinear state equation and measurement equation model using an unscented Kalman filter algorithm; the initial value of the measurement sequence in the nonlinear state equation and measurement equation model is initialized by the harmonic voltage collected at the first sampling moment;

[0172] The measurement update module is used to update the nonlinear state equation and measurement equation model based on the time update results and the Sage-Husa method;

[0173] The state update module is used to update the state of the nonlinear state equation and the measurement equation model based on the measurement update result and the preset forgetting factor to obtain the estimated value of the harmonic current.

[0174] The time update module includes: a traceless transformation unit and a time update unit;

[0175] An untraceable transformation unit is used to perform an untraceable transformation on the harmonic current in the initialized nonlinear state equation and measurement equation model to obtain a Sigma matrix containing a Sigma point set;

[0176] The time update unit is used to perform nonlinear transformation on the Sigma matrix according to the nonlinear state equation, and use the square root of the covariance of the system process noise to participate in the recursive operation of the unscented Kalman filter algorithm to predict the harmonic current and the harmonic current covariance.

[0177] The measurement update module includes: a resampling unit, a measurement residual unit, an observation noise update unit, a measurement update unit and a filter gain unit;

[0178] A resampling unit is used to resample the Sigma point based on the result of time update to obtain the Sigma point;

[0179] The measurement residual unit is used to perform nonlinear transformation on the Sigma point through the measurement equation and calculate the measurement residual;

[0180] An observation noise update unit is used to update the square root of the system observation noise covariance based on the measurement residual and taking into account the preset forgetting factor, using the square root unscented Kalman filter method based on Sage-Husa;

[0181] a measurement update unit for updating the covariance between harmonic voltages and harmonic currents and the covariance of harmonic voltages based on the square root of the system observation noise covariance;

[0182] The filter gain unit is used to calculate the Kalman filter gain according to the covariance between the harmonic voltage and the harmonic current and the covariance of the harmonic voltage.

[0183] Among them, the system also includes a current and process noise update module and a parameter transfer module;

[0184] A current and process noise update module is used to update the covariance estimate of the harmonic current and the square root of the system process noise covariance vector;

[0185] The parameter transfer module is used to transfer the parameters of the nonlinear state equation and measurement equation model at the current sampling moment to the nonlinear state equation and measurement equation model at the next sampling moment.

[0186] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0187] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0188] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0189] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0190] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application and are not intended to limit its scope of protection. Although the present application has been described in detail with reference to the above embodiments, ordinary technicians in the relevant field should understand that after reading this application, those skilled in the art may still make various changes, modifications or equivalent substitutions to the specific implementation methods of the application, but these changes, modifications or equivalent substitutions are all within the scope of protection of the pending claims of the application.

Claims

1. An improved dynamic harmonic estimation method, characterized in that: include: Collect harmonic voltages of synchronized phasor measurement devices on each branch or bus of the power grid; Substituting the harmonic voltage into a pre-established nonlinear state equation and measurement equation model; An improved Sage-Husa square root unscented Kalman filter algorithm is used to perform harmonic current estimation calculation on the branch to be estimated based on the nonlinear state equation and the measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated; The improved Sage-Husa square root unscented Kalman filter algorithm includes: updating the state of the nonlinear state equation and the measurement equation model based on measurement update and forgetting factor; The improved Sage-Husa square root unscented Kalman filter algorithm is used to perform harmonic current estimation calculation on the preset branch to be estimated based on the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated, including: Based on the nonlinear state equation and measurement equation model, an unscented Kalman filter algorithm is used for time update; the initial value of the measurement sequence in the nonlinear state equation and measurement equation model is initialized by the harmonic voltage collected at the first sampling moment; Based on the time update result and the preset forgetting factor, the nonlinear state equation and measurement equation model are measured and updated based on the Sage-Husa method; Based on the measurement update result, the nonlinear state equation and the measurement equation model are updated to obtain an estimated value of the harmonic current; The method of performing measurement update on the nonlinear state equation and measurement equation model based on the time-based update result and the preset forgetting factor and the Sage-Husa method includes: Based on the result of time update, the Sigma point is resampled to obtain the Sigma point; Perform nonlinear transformation on Sigma points through measurement equation and calculate measurement residuals; Based on the measurement residual and taking into account the preset forgetting factor, the square root of the system observation noise covariance is updated using the square root unscented Kalman filtering method based on Sage-Husa; Based on the square root of the system observation noise covariance, the covariance between the harmonic voltage and harmonic current and the covariance of the harmonic voltage are updated; Calculate the Kalman filter gain based on the covariance between the harmonic voltage and the harmonic current and the covariance of the harmonic voltage; According to the measurement residual and considering the preset forgetting factor, the square root of the system observation noise covariance is updated using the square root unscented Kalman filtering method based on Sage-Husa, as shown in the following formula: Where, Represents the square root of the system observation noise covariance at the next sampling moment, R ** k Indicates correspondence The first intermediate process variable, R * k Indicates correspondence The second intermediate process variable, represents the square root of the system observation noise covariance at the current sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, d k is the forgetting correction factor at the next sampling moment calculated based on the forgetting factor, Represents the weight of the i-th Sigma point when calculating the covariance, It represents the harmonic voltage after the weighted correction and measurement transformation of the i-th Sigma point at the next sampling moment. It represents the harmonic voltage of the i-th Sigma point after adding weights and correction and measurement transformation at the current sampling moment, R * k T Represents R * k The transpose of , Cholupdata is the Cholesky decomposition operation, and diag is the operation of constructing a diagonal matrix; The forgetting correction factor d at the next sampling moment k The calculation formula is as follows: d k =(1-b) / (1-b k+1 ) Where b is the forgetting factor.

2. The method according to claim 1, wherein The nonlinear state equation and measurement equation model are shown below: Wherein, subscript k represents the next sampling time, subscript k-1 represents the current sampling time; X k Represents the estimated value of harmonic current at the next sampling moment, X k-1 Indicates the estimated value of harmonic current at the current sampling moment, U k represents the harmonic voltage at the next sampling moment, F is the system state function, H is the system observation function, and A k Represents the system matrix at the next sampling moment, B k Represents the control matrix at the next sampling moment, C k Represents the observation matrix at the next sampling moment, D k represents the direct transfer matrix at the next sampling moment, W k represents the system process noise at the next sampling moment, V k represents the system observation noise at the next sampling moment, Z k Indicates the harmonic voltage observation value at the next sampling moment.

3. The method according to claim 1, wherein Initialize the initial values ​​of the measurement sequence in the nonlinear state equation and measurement equation model as shown below: Wherein, subscript 0 represents the initial sampling time, X0 represents the initial value of harmonic current, Q0 represents the covariance vector of system process noise at the initial sampling time, R0 represents the covariance vector of system observation noise at the initial sampling time, is the square root of Q0, is the square root of R0, is the mathematical expectation of X0; E represents the mathematical expectation, Chol function is the matrix Cholesky factor decomposition function, S0 is The covariance with X0; The relationship between the initial value X0 of the harmonic current and the harmonic voltage is as follows: Z0=H[X0,U0]=C0X0+D0U0+V0 Where Z0 represents the harmonic voltage after measurement transformation at the initial sampling time, U0 represents the harmonic voltage at the initial sampling time, V0 represents the system observation noise at the initial sampling time, H is the observation function, C is the observation matrix, and D is the direct transfer matrix.

4. The method according to claim 1, wherein The method of using an unscented Kalman filter algorithm to perform time updating based on the nonlinear state equation and measurement equation model includes: Performing an untraceable transformation on the harmonic current in the initialized nonlinear state equation and measurement equation model to obtain a Sigma matrix containing a Sigma point set; The Sigma matrix is ​​nonlinearly transformed according to the nonlinear state equation, and the square root of the covariance of the system process noise is used to participate in the recursive operation of the unscented Kalman filter algorithm to predict the harmonic current and the harmonic current covariance.

5. The method according to claim 4, wherein The calculation formula of the recursive operation is as follows: Where, Indicates the harmonic current intermediate value at the next sampling moment, χ i,k-1 Indicates the value of the i-th Sigma point corresponding to the harmonic current at the current sampling moment. The value range of i is 0 to 2n, and n is the number of dimensions of the harmonic current. U k-1 Indicates the harmonic voltage at the current sampling moment, A k-1 Represents the system matrix at the current sampling moment, B k-1 represents the control matrix at the current sampling moment, F is the system state function, Indicates the weight of the i-th Sigma point when calculating the mean, Represents the weight of the i-th Sigma point when calculating the covariance, Indicates the corresponding χ i,k-1 The intermediate process variable, S xx,k Indicates the covariance of the harmonic current at the next sampling moment, Q k represents the system process noise covariance vector at the next sampling moment, Indicates the weight of the 0th Sigma point in calculating the covariance, χ * 0,k-1 Indicates the intermediate process variable corresponding to the 0th Sigma point corresponding to the harmonic current at the current sampling moment, Indicates the corresponding S xx,k The intermediate process variables, qr is the orthogonal triangular matrix decomposition operation, Cholupdata is the Cholesky decomposition operation; The value of the i-th Sigma point corresponding to the harmonic current at the current sampling moment is χ i,k-1 As shown in the following formula: Where σ k-1 Indicates the traceless change parameter at the current sampling moment, Indicates the mean value of harmonic current at the current sampling moment; The calculation formula of the traceless change parameter at the current sampling moment is as follows: Where S xx,k-1 represents the covariance of the harmonic current at the current sampling moment, and λ is the scaling factor.

6. The method according to claim 1, wherein The resampled Sigma point set is nonlinearly transformed and the measurement residual is calculated by the measurement equation, as shown in the following formula: Where H is the observation function, χ i,k Indicates the value of the i-th Sigma point corresponding to the harmonic current at the next sampling moment of resampling, U k-1 Indicates the harmonic voltage at the current sampling moment, The weight of the i-th Sigma point when calculating the mean, It represents the harmonic voltage after adding weights, correction and measurement transformation at the next sampling moment, Z k Indicates the harmonic voltage observation value collected at the next sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, Z * i,k Indicates correspondence The intermediate process variable of the i-th Sigma point.

7. The method according to claim 1, wherein The covariance between the harmonic voltage and the harmonic current and the covariance of the harmonic voltage are updated based on the square root of the system observation noise covariance, as shown in the following formula: Where, P xz,k It represents the covariance between harmonic voltage and harmonic current at the next sampling moment, S zz,k Represents the covariance of the harmonic voltage at the next sampling moment, S * zz,k Indicates the corresponding S zz,k The intermediate process variables, It represents the harmonic voltage after adding weights, correction and measurement transformation at the next sampling moment, Z * i,k Indicates correspondence The intermediate process variable of the i-th Sigma point, Indicates the intermediate value of harmonic current obtained at the next sampling moment, χ * i,k Indicates the corresponding χ i,k The intermediate process variable, χ i,k Indicates the value of the i-th Sigma point corresponding to the harmonic current at the next sampling moment, Indicates the weight of the i-th Sigma point when calculating the covariance, R * k Indicates correspondence The second intermediate process variable, Indicates the weight of the 0th Sigma point when calculating the covariance, Z * 0,k-1 Indicates correspondence The intermediate process variable of the 0th Sigma point, It represents the harmonic voltage after weighted correction and measurement transformation at the current sampling moment, qr is the orthogonal triangular matrix decomposition operation, and Cholupdata is the Cholesky decomposition operation.

8. The method according to claim 7, wherein The calculation formula of the Kalman filter gain is as follows: Where K k Represents the Kalman filter gain at the next sampling moment, S T zz,k Indicates S zz,k The transpose of .

9. The method according to claim 1, wherein The step of updating the state of the nonlinear state equation and the measurement equation model based on the measurement update result to obtain an estimated value of the harmonic current includes: Calculation is performed based on the harmonic current intermediate value, measurement residual and Kalman filter gain in the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current.

10. The method according to claim 9, wherein The estimated value of the harmonic current is calculated as follows: Among them, X k Represents the estimated value of the harmonic current at the next sampling moment, Indicates the harmonic current intermediate value obtained at the next sampling moment, K k Represents the Kalman filter gain at the next sampling moment, e k Indicates the measurement residual of the harmonic voltage at the next sampling moment.

11. The method according to claim 9, wherein After obtaining the estimated value of the harmonic current, the method further includes: updating the covariance estimate of the harmonic current and updating the square root of the system process noise covariance vector; The parameters of the nonlinear state equation and measurement equation model at the current sampling moment are transferred to the nonlinear state equation and measurement equation model at the next sampling moment.

12. The method according to claim 11, wherein The calculation formula of the covariance estimation value of the harmonic current is as follows: S k =cholupdate{S k-1 ,K k S zz,k ,-1} Among them, S k It represents the covariance estimate of the harmonic current at the next sampling moment, S k-1 Represents the covariance estimate of the harmonic current at the current sampling moment, K k Represents the Kalman filter gain at the next sampling moment, S zz,k Represents the covariance of the harmonic voltage at the next sampling moment, and Cholupdata is the Cholesky decomposition operation; The square root of the system process noise covariance vector is calculated as follows: in, represents the square root of the system process noise covariance vector at the next sampling moment, Q k-1 represents the system process noise covariance vector after the k-1th recursion, d k is the forgetting correction factor at the next sampling moment, K k represents the Kalman filter gain at the next sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, S zz,k represents the covariance of the harmonic voltage after the kth recursion, Q ** k Indicates correspondence The first intermediate process variable, Q * k Indicates correspondence The second intermediate process variable, Cholupdata is the Cholesky decomposition operation, and diag is the diagonal matrix construction operation.

13. The method according to claim 11, wherein The transferring of the parameters of the nonlinear state equation and measurement equation model at the current sampling moment to the nonlinear state equation and measurement equation model at the next sampling moment includes: The intermediate value of the harmonic current at the current sampling moment is used as the intermediate value of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, the covariance estimate of the harmonic current at the current sampling moment is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, the square root of the system process noise covariance vector at the current sampling moment is used as the covariance estimate of the harmonic current of the nonlinear state equation and measurement equation model at the next sampling moment, and the square root of the system observation noise covariance at the current sampling moment is used as the square root of the system observation noise covariance of the nonlinear state equation and measurement equation model at the next sampling moment.

14. An improved dynamic harmonic estimation system, characterized in that: include: Data acquisition module, data input module and harmonic estimation module; The data acquisition module is used to collect the harmonic voltages of the synchronized phasor measurement devices of each branch or bus of the power grid; The data input module is used to bring the harmonic voltage into a pre-established nonlinear state equation and measurement equation model; The harmonic estimation module is used to use an improved Sage-Husa square root unscented Kalman filter algorithm to perform harmonic current estimation calculation on the branch to be estimated based on the nonlinear state equation and measurement equation model to obtain an estimated value of the harmonic current of the branch to be estimated; The improved Sage-Husa square root unscented Kalman filter algorithm includes: updating the state of the nonlinear state equation and the measurement equation model based on measurement update and forgetting factor; The harmonic estimation module includes: a time update module, a measurement update module and a state update module; The time update module is used to perform time update based on the nonlinear state equation and measurement equation model using an unscented Kalman filter algorithm; the initial value of the measurement sequence in the nonlinear state equation and measurement equation model is initialized by the harmonic voltage collected at the first sampling moment; The measurement update module is used to perform measurement update on the nonlinear state equation and measurement equation model based on the time update result and the Sage-Husa method; The state updating module is used to update the state of the nonlinear state equation and the measurement equation model based on the measurement update result and the preset forgetting factor to obtain an estimated value of the harmonic current; The measurement update module includes: a resampling unit, a measurement residual unit, an observation noise update unit, a measurement update unit and a filter gain unit; A resampling unit is used to resample the Sigma point based on the result of time update to obtain the Sigma point; The measurement residual unit is used to perform nonlinear transformation on the Sigma point through the measurement equation and calculate the measurement residual; An observation noise update unit is used to update the square root of the system observation noise covariance based on the measurement residual and taking into account the preset forgetting factor, using the square root unscented Kalman filter method based on Sage-Husa; a measurement update unit for updating the covariance between harmonic voltages and harmonic currents and the covariance of harmonic voltages based on the square root of the system observation noise covariance; A filter gain unit, configured to calculate a Kalman filter gain based on the covariance between the harmonic voltage and the harmonic current and the covariance of the harmonic voltage; According to the measurement residual and considering the preset forgetting factor, the square root of the system observation noise covariance is updated using the square root unscented Kalman filtering method based on Sage-Husa, as shown in the following formula: Where, Represents the square root of the system observation noise covariance at the next sampling moment, R ** k Indicates correspondence The first intermediate process variable, R * k Indicates correspondence The second intermediate process variable, represents the square root of the system observation noise covariance at the current sampling moment, e k Represents the measurement residual of the harmonic voltage at the next sampling moment, d k is the forgetting correction factor at the next sampling moment calculated based on the forgetting factor, Represents the weight of the i-th Sigma point when calculating the covariance, It represents the harmonic voltage after the weighted correction and measurement transformation of the i-th Sigma point at the next sampling moment. It represents the harmonic voltage of the i-th Sigma point after adding weights and correction and measurement transformation at the current sampling moment, R * k T Represents R * k The transpose of , Cholupdata is the Cholesky decomposition operation, and diag is the operation of constructing a diagonal matrix; The forgetting correction factor d at the next sampling moment k The calculation formula is as follows: d k =(1-b) / (1-b k+1 ) Where b is the forgetting factor.

Images