A harmonic state estimation method for non-gaussian distributed measurement errors

Abstract

The application discloses a harmonic state estimation method suitable for non-Gaussian distribution measurement error, and comprises the following steps: obtaining real-time measurement data of a power distribution network, and constructing a linear harmonic state estimation model, wherein the measurement error of the real-time measurement data is modeled as a mixed Gaussian distribution; using historical measurement data, extracting cross-harmonic correlation between different harmonic orders through a vine Copula self-encoder, generating a prior probability distribution, and calculating a prior expectation and a prior covariance matrix; based on Bayes theorem, combining the mixed Gaussian distribution and the prior probability distribution, constructing a maximum a posteriori target function; and using an expectation maximization algorithm to iteratively solve the maximum a posteriori target function, and outputting a harmonic state estimation value. Through the introduction of the cross-harmonic correlation prior and the Bayes inference mechanism, the estimation accuracy and numerical stability in the measurement sparse and strong noise environment are fundamentally improved.

Description

Technical Field

[0001] This invention relates to the field of power system automation and state estimation technology, and in particular to a harmonic state estimation method that adapts to non-Gaussian distribution measurement errors. Background Technology

[0002] In an ideal AC power system, the waveforms of voltage and current should be pure sine waves of a single frequency. However, with the large-scale integration of nonlinear loads such as high proportions of renewable energy and power electronic devices into modern power grids, the original sine waves become distorted and warped. According to the Fourier transform principle, this distorted non-sine wave can be decomposed into a fundamental frequency component and a series of sine wave components with frequencies that are integer multiples of the fundamental frequency. These components with frequencies that are integer multiples of the fundamental frequency are called harmonics.

[0003] Harmonic state estimation refers to the process of using limited measurement data from the power grid to calculate and estimate the actual operating state of voltage and current at various nodes in the power grid at different harmonic frequencies through mathematical algorithms. In AC power systems, phasors are used to represent the amplitude and phase angle of voltage or current at a specific frequency. Using phasors for calculation can avoid complex nonlinear power equations and achieve complete linearization of the state estimation model.

[0004] Synchronous phasor measurement units (SPMs) are advanced power measurement devices capable of providing high-precision, microsecond-level synchronous three-phase voltage and current phasor data with GPS timestamps. Miniature synchronous phasor measurement units specifically refer to miniature phasor measurement devices deployed at critical nodes in distribution networks. Due to the insufficient number of actual measurement devices deployed in distribution networks, historical data or predicted values ​​are often artificially introduced as pseudo-measurements to meet the number of equations required for state estimation.

[0005] Traditional statistical algorithms typically assume that measurement errors follow a standard normal distribution. However, under harsh real-world conditions, due to packet loss, transformer magnetic saturation, and the introduction of numerous spurious measurements, error data often becomes distorted and no longer conforms to a standard Gaussian distribution. This characteristic is known as non-Gaussian distribution measurement error. Furthermore, different types of harmonic sources possess specific spectral characteristics, and there are inherent physical correlations and statistical dependencies between harmonics of different orders. For example, there is correlation between even-order harmonics and propagation characteristics between specific harmonic multiples. These inherent physical correlations and statistical dependencies between harmonics of different orders are known as interharmonic correlation.

[0006] Existing technologies for harmonic state estimation in distribution networks have the following main drawbacks: First, they rely excessively on the Gaussian distribution error assumption, and the estimation accuracy decreases significantly when the actual error exhibits heavy-tailed non-Gaussian characteristics. Second, they treat each harmonic as a completely independent variable, ignoring the inherent correlation between different harmonics, resulting in low utilization efficiency of broadband measurement information. Third, they do not make sufficient use of historical prior information, mostly employing only a simple Gaussian prior distribution, failing to fully utilize the complex harmonic characteristics and statistical laws contained in massive historical data. Under conditions of sparse measurements or poor data quality, ill-conditioned solution problems are prone to occur, and the overall robustness of the solution is poor. Summary of the Invention

[0007] In view of this, the present invention provides a harmonic state estimation method that adapts to non-Gaussian distribution measurement errors, which can achieve accurate and robust harmonic state estimation in environments with strong noise and sparse measurement.

[0008] To achieve the above objectives, the present invention provides the following technical solution:

[0009] A harmonic state estimation method adapted to non-Gaussian distribution measurement errors includes the following steps: acquiring real-time measurement data of the distribution network and constructing a linear harmonic state estimation model, wherein the measurement error of the real-time measurement data is modeled as a mixture Gaussian distribution; using historical measurement data, extracting the cross-harmonic correlation between different harmonic orders through a Copula autoencoder to generate a prior probability distribution, and calculating the prior expectation and prior covariance matrix; based on Bayes' theorem, combining the mixture Gaussian distribution and the prior probability distribution, constructing a maximum a posteriori objective function; and using the expectation-maximization algorithm to iteratively solve the maximum a posteriori objective function to output the harmonic state estimate.

[0010] Preferably, the step of modeling the measurement error as a mixture of Gaussian distributions includes: decoupling the measurement error into Gaussian thermal noise and non-Gaussian transient interference; and quantifying the decoupled measurement error using a mixture of Gaussian impulse noise model, wherein the mixture of Gaussian impulse noise model is represented as a weighted sum of multiple Gaussian probability density functions, each Gaussian probability density function corresponding to an error component and having corresponding weight coefficients, mean, and covariance.

[0011] Preferably, the step of using historical measurement data to extract interharmonic correlations between different harmonics through a Teng Copula autoencoder, generating a prior probability distribution, and calculating the prior expectation and prior covariance matrix includes: inputting historical measurement data into the encoder network of the autoencoder to compress it into low-dimensional latent variables; mapping the low-dimensional latent variables to uniformly distributed variables through an empirical cumulative distribution function; constructing an adaptive Teng Copula decomposition to fit the uniformly distributed variables to a multivariate joint distribution to explicitly quantify interharmonic correlations; sampling new latent variables from the fitted joint distribution and mapping them back to the original high-dimensional physical space through a decoder network to reconstruct a high-dimensional sparse prior probability distribution; and calculating the prior expectation and prior covariance matrix based on the high-dimensional sparse prior probability distribution.

[0012] Preferably, the step of constructing the adaptive Copula decomposition includes: for low-order strongly coupled harmonics, Clayton Copula is selected for modeling; for high-order harmonics, Gumbel Copula is selected for modeling.

[0013] Preferably, the step of calculating the prior expectation and prior covariance matrix based on the high-dimensional sparse prior probability distribution includes: constructing an integral expression for the high-dimensional sparse prior probability distribution, wherein the integral expression is used to solve for the prior expectation and prior covariance matrix; and using a multi-dimensional iterative lattice rule to accelerate the solution of the integral expression.

[0014] Preferably, the step of accelerating the solution of the integral expression by using multi-dimensional iterative lattice rules includes: performing an affine transformation on the coordinates of the integral domain to convert the quasi-Monte Carlo lattice rules into a tensor product grid form; constructing a sign function layer by layer along the coordinate direction in the tensor product coordinate system; and reusing intermediate calculation results through iterative dimensionality reduction recursive formulas to extract the prior expectation and prior covariance matrix.

[0015] Preferably, the step of iteratively solving the maximum a posteriori objective function using the expectation-maximization algorithm includes: initializing the state vector as the prior expectation; alternately executing the E-step and M-step until the convergence condition is met: E-step: fix the current state estimate and calculate the responsibility weights corresponding to each error component in the Gaussian mixture distribution; M-step: fix the responsibility weights, solve the weighted least squares equation system containing the prior covariance matrix and the prior expectation, and update the state estimate.

[0016] Preferably, in step E, the formula for calculating the responsibility weight is: ;in, For the first The responsibility weight of each error component For the first The weighting coefficients of each error component. For the first The weighting coefficients of each error component. This represents the Gaussian probability density function. This represents the total number of Gaussian components in the mixture. For real-time measurement vectors, For the measurement matrix, For state vectors, and The first The mean and covariance of each error component, and The first The mean and covariance of each error component.

[0017] Preferably, in the M steps, the update equation for the state vector is: ;in, This represents the total number of Gaussian components in the mixture. For the measurement matrix, For real-time measurement vectors, For the first The mean of each error component, For the first The responsibility weight of each error component For the first The covariance matrix of each error component For the equivalent dynamic weight matrix, It is the inverse of the prior covariance matrix. Let be the a priori expectation.

[0018] Preferably, the real-time measurement data includes measurement data from the micro synchronous phasor measurement unit and measurement data from the smart meter.

[0019] The present invention also proposes an electronic device comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores a computer program executable by the at least one processor, the computer program being executed by the at least one processor to enable the at least one processor to perform the harmonic state estimation method adapted to non-Gaussian distribution measurement error as described in any of the above embodiments.

[0020] The present invention also proposes a computer-readable storage medium storing computer instructions for causing a processor to execute the harmonic state estimation method adapted to non-Gaussian distribution measurement errors as described in any of the above embodiments.

[0021] The present invention also proposes a computer program product, which includes a computer program that, when executed by a processor, implements a harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to any one of the above embodiments.

[0022] The beneficial effects of this invention are as follows:

[0023] First, by using a Bayesian framework and a mixture of Gaussian impulse noise model, the complex measurement error is explicitly decoupled into Gaussian thermal noise and non-Gaussian transient interference. In the iteration of the expectation-maximization algorithm, the current residual is dynamically evaluated to determine which error type it is more in order to adaptively adjust the responsibility weight. This allows for flexible and accurate handling of non-Gaussian measurement errors, completely avoiding the dependence on the Gaussian assumption and demonstrating strong robustness in scenarios with strong noise and transient interference.

[0024] Second, it innovatively introduces the Copula autoencoder to perform nonlinear compression and distribution fitting on massive historical data. It accurately matches Clayton Copula or Gumbel Copula to the physical characteristics of different harmonics, explicitly quantifies the asymmetric nonlinear statistical dependence between different harmonic orders, breaks down the information barrier of independent processing of each frequency band, makes full use of the coupling characteristics between harmonics, and greatly improves the utilization efficiency of broadband measurement information and the accuracy of overall state estimation.

[0025] Third, by using multi-dimensional iterative lattice rules, the prior expectation and cross-band correlation covariance matrix are efficiently solved, and the historical data depth is transformed into a prior regularization constraint and introduced into the state update equation. When the system has insufficient observability or faces strong noise interference, the state vector is guided to automatically converge to the pre-mined physical high-dimensional manifold. This not only significantly enhances the stability of the method in scenarios with missing data or sparse measurements, but also completely solves the underlying mathematical problem of ill-conditioned matrix inversion under harsh conditions in traditional methods. Detailed Implementation

[0026] Embodiments of the present invention are described in detail below, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to reference are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0027] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0028] The following describes a harmonic state estimation method adapted to non-Gaussian distribution measurement errors in an embodiment of the present invention.

[0029] This application discloses a harmonic state estimation method adapted to non-Gaussian distribution measurement errors, comprising the following steps: acquiring real-time measurement data of the distribution network, constructing a linear harmonic state estimation model, wherein the measurement error of the real-time measurement data is modeled as a mixture Gaussian distribution; using historical measurement data, extracting cross-harmonic correlations between different harmonic orders through a Copula autoencoder, generating a prior probability distribution, and calculating the prior expectation and prior covariance matrix; based on Bayes' theorem, combining the mixture Gaussian distribution and the prior probability distribution, constructing a maximum a posteriori objective function; using the expectation-maximization algorithm to iteratively solve the maximum a posteriori objective function, and outputting the harmonic state estimate.

[0030] Specifically, in a three-phase symmetrical distribution network, the transmission line can be simplified to a single-phase line, using... A model is constructed using an equivalent circuit. Based on Kirchhoff's laws, the linear model for harmonic state estimation can be written as follows: ,in The measurement vector contains node voltage phasors and line current phasors; The measurement matrix is ​​constructed based on network topology and line parameters; It is a state vector containing the voltage phasor of each harmonic order at each node; This is the measurement error vector. Measurement vector It integrates real-time measurement data from multi-source heterogeneous sensors at the bottom layer of the distribution network, specifically including high-precision synchronous three-phase voltage and current phasor data provided by miniature synchronous phasor measurement units deployed at key nodes, as well as asynchronous measurement data or historical pseudo-measurement data provided by widely distributed advanced metering infrastructure or smart meters.

[0031] In this step, the measurement error The model is constructed using a Gaussian mixture distribution, rather than the traditional Gaussian distribution assumption. Subsequently, using historical measurement data, the interharmonic correlations between different harmonic orders are extracted using a Copula autoencoder to generate a prior probability distribution, and the prior expectation and prior covariance matrix are calculated. Next, based on Bayes' theorem, and combining the Gaussian mixture distribution and the prior probability distribution, a maximum a posteriori objective function is constructed. Finally, the expectation-maximization algorithm is used to iteratively solve this objective function, outputting the harmonic state estimate.

[0032] This embodiment constructs a complete harmonic state estimation framework adapted to non-Gaussian errors, breaking through the limitations of traditional Gaussian error assumptions. By introducing cross-harmonic correlation priors and Bayesian inference mechanisms, it fundamentally improves estimation accuracy and numerical stability in measurement-sparse and high-noise environments. Compared to existing technologies that process each harmonic independently and rely on a single Gaussian prior, the technical solution in this embodiment can fully utilize broadband measurement information and the statistical regularities of historical data, achieving a comprehensive innovation from the underlying mathematical model to the upper-level solution algorithm.

[0033] In some embodiments, the step of modeling the measurement error as a mixture of Gaussian distributions includes: decoupling the measurement error into Gaussian thermal noise and non-Gaussian transient interference; quantizing the decoupled measurement error using a mixture of Gaussian impulse noise model, wherein the mixture of Gaussian impulse noise model is represented as a weighted sum of multiple Gaussian probability density functions, each Gaussian probability density function corresponding to an error component and having a corresponding weight coefficient, mean, and covariance.

[0034] Specifically, the measurement error is decoupled into two components: Gaussian thermal noise and non-Gaussian transient interference. Gaussian thermal noise mainly originates from the thermal disturbance of the electronic components of the measuring equipment itself, exhibiting standard Gaussian distribution characteristics; while non-Gaussian transient interference originates from adverse operating conditions such as communication packet loss, transformer magnetic saturation, and the introduction of a large amount of spurious measurement data, exhibiting heavy-tailed non-Gaussian characteristics.

[0035] A mixture of Gaussian impulse noise (JEV) model is used to quantify the decoupled measurement error. This model is expressed as a weighted sum of multiple Gaussian probability density functions, mathematically represented as: the measurement error in the nth... The probability distribution of each element is equal to the probability distribution of the elements from the 1st to the 2nd. The weighting coefficients of each Gaussian mixture component are multiplied by the sum of their respective Gaussian probability density functions, where each Gaussian probability density function has its own mean and covariance. The weighting coefficients satisfy the constraint that the sum is 1. This model can accurately capture the heavy-tailed non-Gaussian nature of measurement errors in actual distribution networks, providing an accurate likelihood function basis for subsequent Bayesian inference.

[0036] This embodiment explicitly decouples the measurement error and quantifies it using a Gaussian mixture model, enabling the algorithm to distinguish and handle different types of error sources. This avoids the model mismatch problem of traditional single Gaussian models when faced with non-Gaussian disturbances. This hybrid modeling approach endows the algorithm with adaptability to heavy-tailed noise, providing a mathematical basis for dynamic weight adjustment in the subsequent expectation-maximization algorithm, and significantly improving its robustness under harsh operating conditions.

[0037] In some embodiments, the steps of extracting interharmonic correlations between different harmonics using historical measurement data via a Teng Copula autoencoder, generating a prior probability distribution, and calculating the prior expectation and prior covariance matrix include: inputting historical measurement data into the encoder network of the autoencoder to compress it into low-dimensional latent variables; mapping the low-dimensional latent variables to uniformly distributed variables through an empirical cumulative distribution function; constructing an adaptive Teng Copula decomposition to fit a multivariate joint distribution to the uniformly distributed variables to explicitly quantify interharmonic correlations; sampling new latent variables from the fitted joint distribution and mapping them back to the original high-dimensional physical space through a decoder network to reconstruct a high-dimensional sparse prior probability distribution; and calculating the prior expectation and prior covariance matrix based on the high-dimensional sparse prior probability distribution.

[0038] Specifically, firstly, massive historical state data, such as voltage and current phasors, of each observation node under multiple harmonics are acquired and input into the encoder network of the autoencoder. The encoder network consists of a multilayer sensing mechanism, which compresses the high-dimensional historical data into a low-dimensional latent variable space through nonlinear mapping, obtaining a low-dimensional feature vector. This low-dimensional latent variable effectively removes the linear redundancy of the power grid physically, extracting a few core nonlinear driving factors that dominate the cross-frequency coupling of the system.

[0039] Then, each component of the low-dimensional latent variable is mapped to a marginal variable that follows a uniform distribution through an empirical cumulative distribution function. Based on this, an adaptive vine Copula decomposition architecture is constructed to decompose the joint probability density function into a product of a series of binary conditional Copulas, thereby explicitly quantifying the asymmetric nonlinear statistical dependence between different harmonic orders.

[0040] Next, random sampling is performed from the fitted low-dimensional Copula joint distribution to generate a new latent variable matrix with equivalent statistical properties. Finally, the new latent variables are inversely mapped back to the original high-dimensional phasor physical space using an autoencoder decoder network to generate prior states, thereby explicitly reconstructing the high-dimensional sparse joint prior probability density function containing cross-band correlations.

[0041] Finally, based on this high-dimensional sparse prior probability distribution, the prior expectation and prior covariance matrix are calculated through integration. The prior expectation is the average value of the state variable over the entire prior distribution, and the prior covariance matrix is ​​the integral of the outer product of the state variable and its expected deviation.

[0042] This embodiment achieves nonlinear compression from high-dimensional historical data to low-dimensional latent variables using the Copula autoencoder, and flexibly captures the complex dependency structure between different harmonic orders. The decoded and reconstructed prior distribution not only contains rich information on cross-harmonic correlations, but also provides accurate prior expectations and covariance matrices for Bayesian estimation through subsequent integral calculations. This mechanism overcomes the shortcomings of traditional methods that ignore the coupling relationship between harmonics, enabling prior information to truly reflect the physical characteristics of the power grid and significantly improving the effectiveness of prior constraints.

[0043] In some embodiments, the step of constructing an adaptive vine Copula decomposition includes: for low-order strongly coupled harmonics, selecting Clayton Copula for modeling; and for high-order harmonics, selecting Gumbel Copula for modeling.

[0044] Specifically, in the adaptive Copula decomposition process, the Copula function is precisely matched according to the physical characteristics of different harmonics: for low-order strongly coupled harmonics with strong lower tail correlation characteristics, Clayton Copula is adaptively selected for modeling; for high-order harmonics with burst aggregation characteristics, Gumbel Copula is selected for modeling.

[0045] Clayton Copula effectively captures the lower-tail correlations between variables, meaning that when one variable takes a small value, the other also tends to take a small value. This aligns with the strong coupling propagation characteristics of low-order harmonics under grid faults or heavy loads. Gumbel Copula, on the other hand, excels at characterizing the upper-tail correlations between variables, making it suitable for describing the aggregation characteristics of high-order harmonics under sudden disturbances. Through this targeted matching strategy, adaptive Gumbel Copula decomposition can more accurately quantify the asymmetric nonlinear statistical dependencies between different harmonic orders.

[0046] This embodiment achieves refined modeling of cross-harmonic correlation by selecting different Copula functions for different harmonic characteristics. Compared with the approach of uniformly using a single Copula type, this adaptive matching strategy can more realistically reflect the physical coupling law of different frequency harmonics in the power grid, thereby improving the accuracy of prior distribution fitting and the precision of state estimation.

[0047] In some embodiments, the steps of calculating the prior expectation and prior covariance matrix based on a high-dimensional sparse prior probability distribution include: constructing an integral expression for the high-dimensional sparse prior probability distribution, the integral expression being used to solve for the prior expectation and prior covariance matrix; and using a multi-dimensional iterative lattice rule to accelerate the solution of the integral expression.

[0048] Specifically, after obtaining the high-dimensional sparse prior probability density function, it is necessary to solve for the prior expectation and prior covariance matrix in the original high-dimensional physical space. The calculation of both involves the integration of the high-dimensional joint probability density function. To overcome the computational collapse problem caused by the extremely high dimensionality, a multi-dimensional iterative lattice rule is adopted to accelerate the solution of the integral expression.

[0049] First, we construct integral expressions for high-dimensional sparse prior probability distributions, where the integral expression for the prior expectation is the integral of the state variables and the prior probability density function over the entire space, and the integral expression for the prior covariance matrix is ​​the integral of the outer product of the state variables and the expected deviation over the prior probability density function over the entire space. Since directly solving these high-dimensional integrals is computationally infeasible, we introduce multi-dimensional iterative lattice rules to accelerate the process.

[0050] This embodiment effectively solves the computational challenge of high-dimensional prior integrals by introducing a multi-dimensional iterative lattice rule. This rule significantly reduces the originally exponentially increasing computational complexity, making it possible to efficiently extract prior expectations and covariance matrices in high-dimensional phasor spaces, thus ensuring the real-time performance of subsequent Bayesian inference.

[0051] In some embodiments, the step of accelerating the solution of the integral expression by using a multi-dimensional iterative lattice rule includes: performing an affine transformation on the coordinates of the integral domain to convert the quasi-Monte Carlo lattice rule into a tensor product grid form; constructing a sign function layer by layer along the coordinate direction in the tensor product coordinate system; and reusing intermediate calculation results through an iterative dimensionality reduction recursive formula to extract the prior expectation and prior covariance matrix.

[0052] Specifically, firstly, an affine transformation is performed on the integral domain coordinates. Traditional quasi-Monte Carlo rank-1 lattice rules use a single generator vector set of integration points. To overcome the computational redundancy caused by its non-orthogonal structure, a specific full-rank affine transformation matrix is ​​introduced to linearly remap the integral coordinate system, rigorously converting the underlying quasi-Monte Carlo lattice rules into a tensor product grid form in the transformed coordinate system. The transformed multidimensional numerical integral is approximated as multiple discrete sums along each independent coordinate direction, with each dimension having a corresponding equivalent integration weight.

[0053] Then, in the tensor product coordinate system, instead of independently calculating the function value at each grid point, a sign function is constructed layer by layer along the coordinate direction. Let the... The sign function of a layer is a function of the variables of the current layer and higher layers, and its iterative dimensionality reduction recursive formula is defined as: the first layer... The sign function of the layer is equal to that of the first layer. Layer symbol function in the first The sum of all values ​​at each grid point multiplied by their corresponding weights. Through this multi-dimensional iterative strategy, the evaluation process of the function at the integration nodes is highly clustered. The summation result of each inner-layer variable is treated as a symbolic function stored in memory and repeatedly called in the outer iterations, completely eliminating a large number of overlapping intermediate calculations.

[0054] By leveraging the efficient dimensionality reduction iteration of the sign function, the multi-dimensional iterative lattice rule successfully reduces the computational complexity of high-dimensional joint probability integrals from exponential to polynomial level, ultimately efficiently outputting analytical prior expectations and a sparse covariance matrix containing harmonic crosstalk energy.

[0055] This embodiment achieves a qualitative leap in the computational efficiency of high-dimensional integrals through affine transformation in the integral domain and a hierarchical reuse mechanism of symbolic functions. This technical solution not only ensures the accuracy of the calculation results but also enables the extraction of high-dimensional prior features with acceptable computational resource consumption, laying the foundation for the engineering practical application of the entire harmonic state estimation method.

[0056] In some embodiments, the steps of iteratively solving the maximum a posteriori objective function using the expectation-maximization algorithm include: initializing the state vector as the prior expectation; alternately executing the E-step and M-step until the convergence condition is met: E-step: fix the current state estimate and calculate the responsibility weights corresponding to each error component in the Gaussian mixture distribution; M-step: fix the responsibility weights, solve the weighted least squares equation system containing the prior covariance matrix and the prior expectation, and update the state estimate.

[0057] Specifically, after constructing the global maximum a posteriori objective function, since the objective function contains a logarithmic summation term, it is not possible to directly obtain an analytical solution by taking the derivative of the state vector. Therefore, a latent variable is introduced to construct an expectation-maximization algorithm for solving the problem.

[0058] First, the initial state iteration value is set to the calculated prior expectation, and the iteration convergence threshold and maximum number of iterations are set. Then, the process alternates between the E-step and the M-step.

[0059] In step E, the current state estimate is fixed, and the responsibility weight corresponding to each Gaussian mixture kernel is calculated and updated based on Bayes' theorem. The responsibility weight represents the residual under the current observation residual condition, which is determined by the first... The posterior probability is generated by the Gaussian error components. The physical significance of this step is that the algorithm adaptively evaluates which error type the current residual better matches, thereby dynamically adjusting the weights.

[0060] In the M-step, the responsibility weights obtained in the E-step are fixed. At this point, the originally complex nonlinear maximum a posteriori objective function can be rigorously transformed into a weighted quadratic form of the expected log-likelihood of the complete data. Taking its partial derivative with respect to the state vector and setting the derivative to zero, the adaptive weighted least squares update equation containing prior constraints can be directly derived. Solving this system of linear equations yields the updated state vector.

[0061] Calculate the state change between two consecutive iterations. If the change is less than the convergence threshold or the number of iterations reaches the maximum value, stop the iteration and output the final state estimate as the system's true posterior optimal estimate for the current measurement section.

[0062] This embodiment cleverly solves the Bayesian inference problem with models containing mixed errors using the expectation-maximization algorithm. The alternating E-step and M-step iterations allow the algorithm to adaptively adjust the confidence level for different error components in each update, while the introduction of prior constraints ensures that the state vector converges to the physical high-dimensional manifold when measurement information is insufficient. This scheme completely solves the problems of accuracy degradation and ill-conditioned matrix inversion in traditional weighted least squares methods under non-Gaussian noise.

[0063] In some embodiments, the formula for calculating the responsibility weight in step E is: ;in, For the first The responsibility weight of each error component For the first The weighting coefficients of each error component. For the first The weighting coefficients of each error component. This represents the Gaussian probability density function. This represents the total number of Gaussian components in the mixture. For real-time measurement vectors, For the measurement matrix, For state vectors, and The first The mean and covariance of each error component, and The first The mean and covariance of each error component.

[0064] Specifically, in the first In the nth iteration, the current state estimate is fixed, and the th... The responsibility weight of each error component is equal to the weight coefficient of that component multiplied by the Gaussian probability density function value, which uses the product of the real-time measurement vector and the measurement matrix and the current state estimate as variables, and the mean and covariance of that component as parameters, and then divided by all... The formula is a weighted sum of the corresponding values ​​of each component. It is derived directly from the calculation of the posterior probability in Bayes' theorem and reflects the relative contribution of each error component to the current residual given the current observations and state.

[0065] This embodiment provides a precise mathematical expression for the responsibility weights, giving the E-step a clear algorithmic basis. Through the dynamic calculation of this formula, the algorithm can sense changes in the error distribution in real time and adjust the influence of each component accordingly, thereby achieving true adaptive estimation.

[0066] In some embodiments, the update equation for the state vector in the M-step is: ;in, This represents the total number of Gaussian components in the mixture. For the measurement matrix, For real-time measurement vectors, For the first The mean of each error component, For the first The responsibility weight of each error component For the first The covariance matrix of each error component For the equivalent dynamic weight matrix, It is the inverse of the prior covariance matrix. This is the a priori expectation.

[0067] Specifically, after fixing the responsibility weights, the equivalent dynamic weight matrix is ​​defined as the product of the responsibility weights and the inverse covariance matrices of the corresponding error components. The analytical solution for this state update is to solve a system of linear equations: The coefficient matrix on the left side of the equation consists of two parts. The first part is the transpose of the measurement matrix multiplied by the sum of all equivalent dynamic weight matrices, and then multiplied by the measurement matrix. The second part is the inverse of the prior covariance matrix. The coefficient matrix as a whole is multiplied by the new state vector. The right side of the equation consists of two parts. The first part is the transpose of the measurement matrix multiplied by the sum of the products of each equivalent dynamic weight matrix and the real-time measurement vector minus the mean of the corresponding component. The second part is the inverse of the prior covariance matrix multiplied by the prior expectation.

[0068] In this update equation, the prior covariance inverse matrix term on the left side of the equation and the prior covariance inverse matrix multiplied by the prior expectation term on the right side constitute the prior regularization constraint. These two terms ensure that the state vector automatically converges to the pre-mined physical high-dimensional manifold when observability is insufficient or there is strong noise interference, fundamentally avoiding the matrix singularity or ill-conditioned inversion problem caused by insufficient information in traditional methods.

[0069] This embodiment presents a closed-form solution for M-step state updates, making the algorithm concise and efficient. The introduction of a priori regularization is the core of this embodiment; it transforms the transharmonic correlations inherent in historical data into hard constraints on state updates, greatly enhancing the stability and reliability of the algorithm under extreme conditions.

[0070] In some embodiments, real-time measurement data includes measurement data from micro-synchronous phasor measurement units and measurement data from smart meters. Specifically, real-time measurement data includes measurement data from micro-synchronous phasor measurement units and measurement data from smart meters. Micro-synchronous phasor measurement units are deployed at key nodes in the distribution network, providing high-precision, microsecond-level synchronous three-phase voltage and current phasor data with GPS timestamps. Smart meters are widely distributed in the distribution network, providing asynchronous measurement data with lower time refresh rates or historical pseudo-measurement data. The fusion of these two types of data constitutes a measurement vector, ensuring high-precision synchronous measurement at key nodes while compensating for insufficient measurement coverage through smart meter data.

[0071] This embodiment clarifies the fusion and use of multi-source heterogeneous measurement data, enabling the method of this invention to adapt to the uneven measurement configuration in actual power distribution networks. By organically combining high-precision synchronous measurements and wide-coverage asynchronous measurements, the method reduces system deployment costs while maximizing the utilization of existing measurement resources, thus improving its engineering applicability.

[0072] The present invention also proposes an electronic device comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores a computer program executable by the at least one processor, the computer program being executed by the at least one processor to enable the at least one processor to perform the harmonic state estimation method adapted to non-Gaussian distribution measurement error as described in any of the above embodiments.

[0073] The present invention also proposes a computer-readable storage medium storing computer instructions for causing a processor to execute the harmonic state estimation method adapted to non-Gaussian distribution measurement errors as described in any of the above embodiments.

[0074] The present invention also proposes a computer program product, which includes a computer program that, when executed by a processor, implements a harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to any one of the above embodiments.

[0075] Other configurations and operations of the harmonic state estimation method adapted to non-Gaussian distribution measurement errors according to embodiments of the present invention are known to those skilled in the art and will not be described in detail here.

[0076] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0077] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A harmonic state estimation method that is adaptive to non-Gaussian distribution measurement errors, characterized by, Includes the following steps: Real-time measurement data of the power distribution network is acquired, and a linear harmonic state estimation model is constructed, wherein the measurement error of the real-time measurement data is modeled as a mixture Gaussian distribution. Using historical measurement data, the cross-harmonic correlation between different harmonics is extracted by the Copula autoencoder, a prior probability distribution is generated, and the prior expectation and prior covariance matrix are calculated. Based on Bayes' theorem, and combining the Gaussian mixture distribution and the prior probability distribution, a maximum a posteriori objective function is constructed. The expectation-maximization algorithm is used to iteratively solve the maximum a posteriori objective function and output the harmonic state estimate.

2. The method of claim 1, wherein, The step of modeling the measurement error as a mixture Gaussian distribution includes: The measurement error is decoupled into Gaussian thermal noise and non-Gaussian transient interference. A Gaussian mixture impulse noise model is used to quantify the measurement error after decoupling. The Gaussian mixture impulse noise model is represented as a weighted sum of multiple Gaussian probability density functions. Each Gaussian probability density function corresponds to an error component and has corresponding weight coefficients, mean, and covariance.

3. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 1, characterized in that, The steps of using historical measurement data to extract interharmonic correlations between different harmonic orders through the Copula autoencoder, generating a prior probability distribution, and calculating the prior expectation and prior covariance matrix include: Historical measurement data is input into the encoder network of the autoencoder and compressed to obtain low-dimensional latent variables; The low-dimensional latent variables are mapped to uniformly distributed variables using an empirical cumulative distribution function; An adaptive Copula decomposition is constructed, and the uniformly distributed variables are fitted with a multivariate joint distribution to explicitly quantify the interharmonic correlation. New latent variables are generated by sampling from the fitted joint distribution and mapped back to the original high-dimensional physical space through the decoder network to reconstruct the high-dimensional sparse prior probability distribution. Based on the high-dimensional sparse prior probability distribution, the prior expectation and prior covariance matrix are calculated.

4. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 3, characterized in that, The steps for constructing the adaptive vine Copula decomposition include: For low-order strongly coupled harmonics, Clayton Copula is selected for modeling. For higher harmonics, Gumbel Copula is selected for modeling.

5. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 3, characterized in that, The steps for calculating the prior expectation and prior covariance matrix based on the high-dimensional sparse prior probability distribution include: Construct an integral expression for the high-dimensional sparse prior probability distribution, which is used to solve for the prior expectation and the prior covariance matrix; The integral expression is solved more quickly by employing a multi-dimensional iterative lattice rule.

6. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 5, characterized in that, The step of accelerating the solution of the integral expression by employing a multi-dimensional iterative lattice rule includes: An affine transformation is performed on the coordinates in the integral domain to convert the quasi-Monte Carlo lattice rule into a tensor product mesh form; In the tensor product coordinate system, a sign function is constructed layer by layer along the coordinate direction. The intermediate calculation results are reused through an iterative dimensionality reduction recursive formula to extract the prior expectation and prior covariance matrix.

7. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 1, characterized in that, The steps of iteratively solving the maximum a posteriori objective function using the expectation-maximization algorithm include: The initial state vector is the prior expectation; Alternately execute the E-step and M-step until the convergence condition is met: Step E: Fix the current state estimate and calculate the responsibility weights corresponding to each error component in the Gaussian mixture distribution; M-step: Fix the responsibility weights, solve the weighted least squares equation system containing the prior covariance matrix and the prior expectation, and update the state estimate.

8. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 7, characterized in that, In step E, the formula for calculating the responsibility weight is: ； in, For the first The responsibility weight of each error component For the first The weighting coefficients of each error component. For the first The weighting coefficients of each error component. This represents the Gaussian probability density function. This represents the total number of Gaussian components in the mixture. For real-time measurement vectors, For the measurement matrix, For state vectors, and The first The mean and covariance of each error component, and The first The mean and covariance of each error component.

9. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 7, characterized in that, In the M steps, the update equation for the state vector is: ； in, This represents the total number of Gaussian components in the mixture. For the measurement matrix, For real-time measurement vectors, For the first The mean of each error component, For the first The responsibility weight of each error component For the first The covariance matrix of each error component For the equivalent dynamic weight matrix, It is the inverse of the prior covariance matrix. Let be the a priori expectation.

10. The harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to claim 1, characterized in that, The real-time measurement data includes measurement data from the micro synchronous phasor measurement unit and measurement data from the smart meter.

11. An electronic device, characterized in that, The electronic device includes: At least one processor; and a memory communicatively connected to the at least one processor; The memory stores a computer program that can be executed by the at least one processor, which enables the at least one processor to perform the harmonic state estimation method for adapting to non-Gaussian distribution measurement errors as described in any one of claims 1-10.

12. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer instructions that, when executed by a processor, implement the harmonic state estimation method for adapting to non-Gaussian distribution measurement errors as described in any one of claims 1-10.

13. A computer program product, characterized in that, The computer program product includes a computer program that, when executed by a processor, implements the harmonic state estimation method for adapting to non-Gaussian distribution measurement errors according to any one of claims 1-10.

Images