A principal component probability modeling method for harmonic sources based on typical power intervals

Abstract

The present invention relates to the technical field of harmonic source modeling, and specifically to a harmonic source principal component probability modeling method based on a typical power interval. Based on the harmonic coupling admittance matrix model, the power interval is divided according to the rated active power of the monitoring point in the substation area, and the probability of occurrence of each power interval within the monitoring time range is calculated, and the initial matrix of each power interval is solved, and a sub-probability density optimization method is proposed. Then, based on the improved moth fire optimized particle swarm algorithm, the parameters of the mixed probability distribution function are identified; the models of different power intervals are compared for similarity, and the probability distributions of elements at the same position are compared at the element level; the obtained similar intervals are clustered at the matrix level to obtain the optimal number of models based on the typical power interval. The model parameters are thresholded through the probability distribution of the diagonal elements of the typical power interval, and the dominant elements under the typical power interval are screened out to achieve effective characterization of the node harmonic emission level.

Description

Technical Field

[0001] The present invention relates to the technical field of harmonic source modeling, and in particular to a harmonic source principal component probability modeling method based on typical power intervals. Background Art

[0002] With the rapid development of renewable energy power generation technology, power electronic components in the power grid have increased widely. They generate a large number of harmonics during operation, which reduces the power quality and causes power loss, and the harm caused is becoming more and more serious. Due to the differences in the topology, control method, and application scenarios of renewable energy grid-connected inverters, the harmonic generation mechanism is becoming more and more complex. At the same time, affected by environmental factors, access location and other factors, the harmonic characteristics of renewable energy such as photovoltaics and wind power and nonlinear loads have obvious randomness and uncertainty. In addition, the harmonic characteristics of the harmonic source are not only related to its own characteristics and operating status, but also affected by adjacent harmonic sources. There is interaction between multiple types of harmonic sources, which aggravates the mutual coupling of harmonics and makes it difficult to characterize the harmonic source model.

[0003] Existing harmonic source modeling methods are mainly divided into mechanism-based harmonic source models and data-driven harmonic source models. Harmonic source mechanism modeling is a harmonic characteristic analysis model constructed based on the circuit structure and operating parameters of harmonic sources such as power electronic converters. It comprehensively considers the influence of various types of parameters on the harmonic characteristics of harmonic sources during the operation of the power system. Its physical meaning is clear and has strong interpretability. Harmonic source mechanism modeling mainly includes: constant current source type harmonic source model, harmonic source model based on cross-frequency admittance matrix, and harmonic source model based on state space. The constant current source model is the simplest type of frequency domain harmonic source model. In this type of model, each harmonic is replaced by a constant current source based on the typical frequency spectrum of the harmonic source. However, this type of model does not consider the coupling relationship between harmonics of different frequencies, has low modeling accuracy, relies on typical harmonic spectra, is not necessarily applicable to actual scenarios, and cannot model dynamic harmonic sources; the harmonic source modeling method based on the cross-frequency admittance matrix starts from the operating mechanism and topological structure of the circuit, has a clear physical meaning, and clearly reflects the coupling relationship between the voltage and current of each harmonic. However, the cross-frequency admittance matrix is ​​based on the fact that the harmonic source can be linearized or approximately linearized. For random time-varying new energy harmonic sources, linearization will bring large errors, so high-precision modeling cannot be achieved. In addition, each element of the matrix corresponds to a specific operating condition, and each element of the matrix needs to be recalculated when the operating condition changes. Therefore, the harmonic source based on the cross-frequency admittance matrix is ​​not suitable for situations where the harmonic source condition changes frequently or the harmonic source has strong randomness; the harmonic source modeling method based on the state space is a practical power electronic equipment modeling method suitable for theoretical analysis. This type of harmonic source model has a clear physical meaning and strong interpretability, and is suitable for analyzing the coupling relationship of harmonics of each frequency. When modeling, the linearized equation is analyzed using small signals to make the model more accurate. However, this type of modeling method also has certain limitations: first, the circuit topology and the parameters of each component must be completely known, which may not be met in actual engineering; second, the constructed model is not universal and portable, and when the topology changes, the model should be re-modeled; third, the modeling steps are cumbersome and not suitable for dynamic modeling; fourth, the model is based on specific operating conditions, so it is not suitable for harmonic sources with frequently changing operating conditions.

[0004] The harmonic data-driven modeling method pays more attention to the external characteristics of the harmonic source. It does not need to consider the operating mechanism of the harmonic source itself and the parameters of the equipment topology. It uses measured data to identify the harmonic model parameters. Therefore, it is relatively universal and more suitable for the construction of harmonic source models with high monitoring levels and known equipment port data. At present, the commonly used harmonic source data-driven models include the Norton equivalent harmonic source model, the harmonic source modeling based on the least squares method, the harmonic source modeling based on the support vector machine, and the harmonic source modeling based on the neural network. The Norton model does not need to know the specific network topology structure, and only requires a few sets of measurement data to be solved, which is convenient for harmonic responsibility division and harmonic source tracing analysis. However, the Norton model is not very interpretable, and the modeling accuracy is easily affected by the measurement data. At the same time, the coupling relationship between harmonics of different frequencies is not considered, and it is impossible to characterize models with strong randomness or variable working conditions. As a classic statistical algorithm, the least squares algorithm has a good statistical theoretical basis, but this method has a more obvious disadvantage, that is, it is easily affected by outliers. In actual power systems, factors such as sudden faults or measurement noise will generate abnormal data. In order to improve modeling accuracy, the least squares algorithm used in the electrical field should improve its robustness. Partial least squares, robust least squares, complex regression least squares and other methods can reduce the impact of abnormal data on modeling results to improve the robustness of the model. Harmonic source modeling based on support vector machine only requires a small amount of data and less computing resources to achieve ideal modeling accuracy. In addition, support vector machine modeling is fast, has excellent anti-noise ability, and is robust, so it is suitable for online modeling. However, the physical meaning of such models is not clear enough, and the interpretability is poor. It is more suitable for simple harmonic source modeling with less source data. Neural network modeling essentially relies on the powerful nonlinear fitting ability of neural networks to map the coupling relationship between each harmonic voltage and current. Compared with other data-driven modeling methods, neural networks have stronger nonlinear mapping capabilities and can achieve high-precision nonlinear complex harmonic source modeling. Its nonlinear fitting ability is strong and it is suitable for complex harmonic source modeling, but such models require large-scale data sets and take a long time to train, so they are not suitable for real-time modeling. If the data set is small or the harmonic source is not highly nonlinear, it is more appropriate to use statistical methods to model it. Therefore, the neural network harmonic source modeling method is more suitable for complex harmonic sources with more historical measurement data.

[0005] In summary, the shortcomings of the prior art are as follows:

[0006] 1) Traditional modeling methods require that the topology of the circuit and the parameters of each component are completely known, but this may not be met in actual engineering;

[0007] 2) The existing model is not universal and transferable. When the topology changes, it should be remodeled;

[0008] 3) The modeling steps are cumbersome and not suitable for dynamic modeling. Moreover, the model is established based on specific operating conditions, so it is not suitable for harmonic sources with frequently changing operating conditions. Summary of the invention

[0009] In view of the above problems, the purpose of the present invention is to provide a harmonic source principal component probability modeling method based on typical power intervals, which can achieve accurate characterization of the harmonic state of each node, meet engineering reality, and have versatility and portability, and is suitable for harmonic sources with frequently changing working conditions.

[0010] The technical solution is as follows:

[0011] A harmonic source principal component probability modeling method based on a typical power interval comprises the following steps:

[0012] Step 1: Based on the harmonic coupling admittance matrix model, taking into account the randomness and stage characteristics of the collective load in the substation area, the power interval is divided according to the rated active power of the monitoring point in the substation area, and the probability of each power interval occurring within the monitoring time range is calculated. The initial matrix of each power interval is solved by the least squares method in the complex domain through the measured data;

[0013] By dividing the power range and building the model considering the load characteristics, the parameter settings are optimized, and the coupling inside the equipment is considered to enhance the system stability and efficiency.

[0014] Step 2: Establish a probability model based on the divided power intervals, and use the improved mothfire optimized particle swarm algorithm to identify the parameters of the mixed probability density distribution function of the probability model to avoid local optimal solutions and improve the convergence speed;

[0015] Step 3: Compare the similarities of the models in different power intervals. At the element level, use the Bhattacharyya distance algorithm to compare the probability distribution of elements at the same position. At the matrix level, define the matrix similarity factor to characterize the similarity of the admittance models in different power intervals, and use the mean shift algorithm to cluster the obtained similar probability intervals to obtain the optimal number of models based on the typical power interval.

[0016] Step 4: Threshold the model parameters through the probability distribution of diagonal elements in the typical power range, screen out the dominant elements in the typical power range, and achieve effective characterization of the node harmonic emission level.

[0017] Furthermore, the step 1 specifically includes:

[0018] Step 1.1: Determine the harmonic coupling admittance matrix model as a unified harmonic source model, the specific form is as follows:

[0019]

[0020] Where:

[0021] are the phasor values ​​of the lth harmonic current and the mth harmonic voltage respectively; Characterizes the contribution of the mth harmonic voltage to the lth harmonic current; l = 1, 3, ..., P; m = 1, 3, ..., P, the highest harmonic is P;

[0022] Step 1.2: Initial power interval model solution

[0023] After dividing the power interval according to the rated active power of the monitoring point in the substation, the harmonic current after the collective equipment in the substation is connected is shown as follows:

[0024]

[0025] Where n is the number of power intervals; [I] K , [Y] K and [U] K They are the contribution of harmonic current and harmonic voltage to harmonic current and the phasor value of harmonic voltage in power interval K respectively;

[0026] Step 1.3: Calculate the probability of each power interval occurring within the monitoring time range

[0027] Assume that the occurrence time corresponding to power interval K is t K , then based on the measured data, the probability of occurrence of this power interval within the monitoring time is expressed as:

[0028]

[0029] In the formula, T is the total monitoring time;

[0030] Step 1.4: Solve the initial matrix for each power interval

[0031] For the admittance matrix of each power interval, the complex domain least squares method is used to solve the initial matrix. The calculation formula is as follows:

[0032]

[0033] In the formula, the superscript “*” indicates the transpose of the matrix, and the current matrix and voltage matrix Generated through monitoring point data.

[0034] Furthermore, the step 2 specifically includes:

[0035] Step 2.1: Select a general model of mixed probability density distribution that characterizes multiple probability distributions:

[0036]

[0037] In the formula, f i (.) is the sub-probability density distribution function, λ i is the weight coefficient of the sub-probability density distribution function, indicating f i (.) is the contribution degree in the model, x is the admittance data of the Kth power interval to be solved, and M is the total number of selected sub-probability density distribution functions; that is, the M sub-probability density distribution functions with the smallest error are selected by the optimization method;

[0038] Step 2.2: Model parameter identification based on improved mothfire optimized particle swarm algorithm

[0039] The objective function is constructed by minimizing the difference between the mean and variance of the probability model and the mean and variance of the actual data, and the sub-probability density distribution function f is weighed by the weight coefficient i (.) contribution; the objective function is expressed as:

[0040]

[0041] In the formula, are the contribution of variance and mean to the objective function, respectively, and their sizes are determined by cross-validation method; Var i (x) and E i (x) represent the sub-probability density distribution function f i The variance and mathematical expectation of (.);

[0042] The weight coefficient λ of the sub-probability density distribution function i The following constraints are met:

[0043]

[0044] λ i >0 (8)

[0045] The output result of the mothfire optimization particle swarm algorithm is set as the initial position information of the particle swarm, and the initial position of the parameter to be identified is continuously updated through the mothfire optimization to avoid individuals with similar performance from being selected to improve the search efficiency; the mothfire optimization model is as follows:

[0046] S(M u ,F v )=D uv e zw cos(2πw)+F v (9)

[0047] Where S(.) is the spiral function, M u represents the u-th moth, F v represents the vth flame, D uv represents the distance D between the uth moth and the vth flameuv =|F v -M u |, z is a constant used to define the shape of the logarithmic spiral, w is a random number in the range [-1,1], and the coefficient w in the spiral function represents the distance between the moth's next position and the flame;

[0048] In order to balance the global search capability and local development capability of the algorithm in the search space, the number of flames will be adaptively reduced during the iteration process. The formula is as follows:

[0049]

[0050] Where w0 represents the current iteration number, Q represents the maximum number of flames, W0 represents the maximum iteration number, and round represents the rounding function;

[0051] Assume that the three selected f i (.) (Take M = 3 as an example) are normal distribution, t-distribution and lognormal distribution respectively; the solution f i The number of unknown parameters of (.) are 2, 3, and 2 respectively; there are also 3 unknown weight coefficients. The total number of unknown parameters j of the model is 10, so the dimension of the solution space of the model is 10. There are N particles in the j-dimensional solution space, each particle is a solution of the model, denoted by the position parameter S, S t d Represents the value of the tth solution at the dth iteration of the model.

[0052] The position of each particle is the solution of the parameters to be identified, S t d Represents the set of parameters of the tth solution at the dth iteration of the model:

[0053]

[0054] In the formula, is the jth parameter value of the tth solution at the dth iteration;

[0055] Each particle has V t d Iterate in the solution space at a speed of and calculate the value of the objective function;

[0056]

[0057]

[0058]

[0059] In the formula, v t,j d is the optimal speed in the solution space at the dth iteration; pbestt d represents the optimal position of the t-th particle in the solution space at the d-th iteration; b t,j d gbest represents the optimal position of the jth parameter of the tth solution as of the dth iteration; d represents the optimal position of all particles in the solution space up to the dth iteration; g j d represents the optimal position of the jth parameter among all solutions up to the dth iteration, that is, the value of the parameter to be identified;

[0060] The velocity and position of each particle after each iteration are expressed as:

[0061]

[0062]

[0063] In the formula, v t d represents the velocity of the t-th particle at the d-th iteration; r1 and r2 are random numbers in [0,1], Represents the solution of each parameter in the tth solution at the dth iteration of the model.

[0064] Furthermore, the neutron probability density distribution function in step 2.1 includes:

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071]

[0072] In the formula, f A (.), f B (.), f C (.), f D (.), f E (.) and f F(.) represent Weibull distribution, lognormal distribution, Beta distribution, normal distribution, t-distribution and Laplace distribution respectively; k and λ represent the shape and scale parameters of Weibull distribution; μ2 and μ4 represent the mathematical expectations of lognormal distribution function and normal distribution function respectively; σ2 and σ4 represent the variances of lognormal distribution function and normal distribution function respectively; α and β are two parameters of Beta distribution; μ5, σ5 and v1 represent the location parameter, scale parameter and shape parameter of t-distribution respectively; μ6 is the location parameter of Laplace distribution, indicating the central position of the distribution, b is the scale parameter of Laplace distribution, and Γ(.) represents the gamma function.

[0073] Furthermore, in step 2.1, the optimization method for selecting M sub-probability density distribution functions with the minimum error is as follows:

[0074] Step 2.1.1: Solve for the parameters of a single sub-probability density distribution function

[0075] The least squares method is used to calculate the probability density distribution function f in the set. i (.) Fitting is performed to obtain the parameters of each probability density distribution function according to the measured data;

[0076] Step 2.1.2: Calculate parameter errors

[0077] Through statistical calculation, the actual mathematical expectation and variance of the actual calculation data are obtained; through the calculation in step 2.1.1, f is obtained i The mathematical expectation and variance of (.); the error between the statistical value and the parameter value obtained in step 2.1.1 is calculated as follows:

[0078]

[0079] In the formula, Var i (x) and E i (x) represents the sub-probability density distribution function f i (.), Var(x) and E(x) represent the actual variance and mathematical expectation of the measured data;

[0080] Step 2.1.3: Screening the best sub-probability density distribution function

[0081] Select M sub-probability density distribution functions f with the smallest error i (.), which is used to form the general probability f(x).

[0082] Furthermore, the step 3 specifically includes:

[0083] Step 3.1: At the element level, use the Bhattacharyya distance algorithm to calculate the overlap of elements at the same position in two power intervals:

[0084]

[0085]

[0086] Where y1(x) and y2(x) are Y lm1 and Y lm2 The percentage of sampling points in the corresponding power range; BC is the Bhattacharyya coefficient, and 0≤BC≤1; W is the waveform similarity coefficient, and 0 <W≤1;Y lm1 is the first power interval in the harmonic coupling admittance matrix Y lm2 is the second power interval in the harmonic coupling admittance matrix element;

[0087] When the waveform similarity coefficient W is greater than the threshold W set When , it means that the similarity between the two waveforms at this position is high, which will be higher than the threshold W set The number and position information of the elements are stored in the set S0.

[0088] Step 3.2: At the matrix level, define the matrix similarity factor δ to quantify the similarity of the matrices:

[0089]

[0090] In the formula, S is the number of similar elements, a0 and b0 are the number of rows and columns of the matrix respectively;

[0091] The value of the matrix similarity factor δ affects the accuracy and calculation speed of the model. 0<δ≤1. The larger the δ value, the higher the accuracy of the model and the lower the calculation efficiency, and vice versa. When δ≥δ set When , the two matrices are considered similar; δ set is the similarity threshold;

[0092] Step 3.3: Use the mean shift algorithm to cluster the obtained similar power intervals

[0093] First, the calculated power interval model is divided into M n The Mean Shift algorithm is used to continuously calculate and optimize the group to select the appropriate M n The power states are taken as cluster centers, and each power interval is divided into different categories according to the distance between it and the cluster center, so as to minimize the distance to the cluster center after the power interval is aggregated, and the optimal number of models based on the typical power interval is obtained;

[0094] After the power interval model is aggregated, the parameters of the model are revalued. The value selection rules are as follows:

[0095] Y lm '=∑Y lmx ×P x (28)

[0096] Where P x is the occurrence probability of the corresponding power interval; Y lm ' is the number of similar power intervals after clustering The value of the element again, Y lmx is the matrix of similar power intervals to be determined.

[0097] Furthermore, the step 4 specifically includes:

[0098] Step 4.1: Threshold the obtained probability model, that is, reduce the order of the obtained harmonic coupling admittance matrix model; for the hth harmonic current, its expression is:

[0099]

[0100] Where Y h,n are all probability density distribution functions, is the phasor value of the nth harmonic voltage;

[0101] Step 4.2: Screening of dominant elements based on the probability density distribution of self-coupling admittance

[0102] Determine the R of the cumulative distribution function of the diagonal elements in the initial admittance matrix n % as the threshold, for each diagonal element, find the R of its cumulative distribution function n % point, for each off-diagonal element in the matrix, calculate its cumulative distribution function value in its own distribution; if this cumulative distribution function value is lower than R n %, then the value of this element is considered to be lower than the threshold, and its probability density function value is set to zero;

[0103] Apply the above process to the entire matrix, screen all off-diagonal elements and set the probability density to zero; the mathematical expression is as follows:

[0104] Let Y be the original admittance matrix, D be a diagonal matrix whose diagonal elements are the same as Y; for each off-diagonal element Y in the original admittance matrix Y lm , if its probability density f lm (Y lm ) corresponds to the cumulative distribution function value F lm (Y lm ) is lower than R n %, then do the following:

[0105]

[0106] Thus, a probabilistic coupling model of the dominant elements in a typical power range is obtained, which enables the characterization of the node harmonic emission level.

[0107] The beneficial effects of the present invention are:

[0108] 1) The present invention proposes a harmonic source probabilistic modeling method based on typical power intervals to achieve accurate characterization of the harmonic state of each node; by dividing the power intervals and building a model considering the load characteristics, using the coupling admittance matrix for modeling, considering the coupling conditions within the current new energy equipment, and enhancing the system stability and efficiency.

[0109] 2) The present invention uses probability distribution to characterize harmonic admittance, characterizes all situations in which admittance exists, and makes the modeling universal and portable. The present invention proposes a method for sub-probability density optimization, by reasonably selecting sub-probability density distribution, using improved mothfire optimized particle swarm algorithm to identify the parameters of the mixed probability density distribution function, while improving the solution speed and taking into account the accuracy, effectively avoiding local optimal solutions, and improving the convergence speed. The method of this step can be migrated and used in other fields.

[0110] 3) In order to improve the practicality of the model, the present invention compares the similarities of models in different power intervals. At the element level, the Bhattacharyya distance algorithm is used to compare the probability distribution of elements at the same position; at the matrix level, the matrix similarity factor is defined to characterize the similarity of the admittance models in different power intervals, and the mean shift algorithm is used to cluster the obtained similar intervals to obtain the optimal number of models based on the typical power interval.

[0111] 4) The present invention thresholds the model parameters through the probability distribution of diagonal elements in a typical power interval, screens out the dominant elements in the typical power interval, and achieves an effective characterization of the node harmonic emission level. BRIEF DESCRIPTION OF THE DRAWINGS

[0112] Figure 1 Flowchart for probabilistic modeling of harmonic sources based on typical power intervals.

[0113] Figure 2 This is the model parameter identification process based on the improved mothfire optimized particle swarm algorithm. DETAILED DESCRIPTION

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

[0115] The present invention proposes a harmonic source probability modeling method based on typical power intervals to achieve accurate characterization of the harmonic state of each node. The flow chart is as follows: Figure 1 As shown. First, based on the harmonic coupling admittance matrix model, taking into account the randomness and stage characteristics of the collective load in the substation, the rated active power of the substation monitoring point is evenly divided into ten power intervals, and the probability of each power interval within the monitoring time range is calculated. The initial matrix of each power interval is solved by the least square method in the complex domain through the measured data. In order to improve the fitting degree of the probability distribution function, the present invention proposes a method for optimizing the mixed probability density distribution to characterize the harmonic emission level and probability characteristics of different power intervals, and uses the improved moth fire optimization particle swarm algorithm to identify the parameters of the mixed probability distribution function, effectively avoiding the local optimal solution and improving the convergence speed; secondly, in order to improve the practicality of the model, the models of different power intervals are compared for similarity. At the element level, the Bhattacharyya distance algorithm is used to compare the probability distribution of elements at the same position. At the matrix level, the matrix similarity factor is defined to characterize the similarity of the admittance model of different power intervals, and the mean shift (MeanShift) algorithm is used to cluster the obtained similar intervals to obtain the optimal number of models based on the typical power interval. Finally, the model parameters are thresholded through the probability distribution of diagonal elements in the typical power range, and the dominant elements in the typical power range are screened out, achieving an effective characterization of the node harmonic emission level.

[0116] 1. Constructing a unified model of harmonic sources

[0117] The harmonic coupling admittance matrix model (HCAM) can intuitively show the coupling relationship between the harmonic current and harmonic voltage of the load. Compared with the traditional constant current source model or Norton model, it can more accurately reflect the harmonic generation of the load. The specific form of the HCAM model based on measured data is as follows:

[0118]

[0119] Where: are the phasor values ​​of the lth harmonic current and the mth harmonic voltage respectively; Characterizes the contribution of the mth harmonic voltage to the lth harmonic current; l=1,3,…,P; m=1,3,…,P, the highest harmonic is P.

[0120] 2. Probabilistic model establishment based on power interval division

[0121] 1. Initial power interval model solution

[0122] Considering that the collective load of the substation is composed of a variety of nonlinear devices, under the same power output, due to the different output conditions of various devices, there may be different harmonic emission characteristics, showing obvious randomness and stage characteristics. Based on this, a power interval model is constructed. The power interval is divided according to the rated active power data of the monitoring point. Assuming the rated power is P n , divided into n intervals. In this embodiment, n is 10, and the divided power interval is <10%P n , 10%~20%P n , 21%~30%P n , 31%~40%P n , 41%~50%P n , 51%~60%P n , 61%~70%P n , 71%~80%P n , 81%~90%P n , >90% P n Then the harmonic current after the collective equipment in the substation area is connected can be expressed as:

[0123]

[0124] According to the power interval division, data of multiple time periods may be obtained. Let the occurrence time corresponding to power interval K be t K , then based on the measured data, the probability of occurrence of this interval within the monitoring time can be expressed as:

[0125]

[0126] Where T is the total monitoring time.

[0127] Taking into account the efficiency and accuracy of the solution, the present invention adopts the complex domain least squares method to solve the initial matrix for the admittance matrix of each power interval. The calculation formula is as follows:

[0128]

[0129] Where: The superscript “*” indicates the transpose of the matrix. Monitor point data generation.

[0130] 2. Determination of probability model

[0131] According to the power interval division, a substation can be divided into 10 power interval states. Based on the data of the Kth interval, a harmonic coupling admittance probability density distribution model can be established to characterize the harmonic emission level of the Kth interval of the substation.

[0132] 1) Optimization of sub-probability density distribution function

[0133] In order to solve the problem that a single probability distribution function is difficult to accurately characterize the distribution of admittance parameters, the present invention proposes a general model of mixed probability density distribution that can characterize the characteristics of multiple probability distributions.

[0134]

[0135] In the formula, f i (.) is the sub-probability density distribution function, λ i is the weight coefficient of the sub-probability density distribution function, indicating f i (.) is the contribution degree in the model, x is the admittance data of the Kth power interval to be solved, and M is the total number of selected sub-probability density distribution functions.

[0136] The present invention proposes a method for optimizing sub-probability density distribution functions, selecting M f from a sub-density distribution function set i (.) to construct the general f(x).

[0137] (1) Constructing a set of sub-probability density distribution functions

[0138] First, construct a set of probability density distribution functions, including normal distribution, Poisson distribution, Laplace distribution and other common probability density distribution functions. There are M elements in the set, M = 6 in the present invention, and the six probability density distribution functions are respectively denoted as f A (.), f B (.), f C (.), f D (.), f E (.), f F (.), are Weibull distribution, lognormal distribution, Beta distribution, normal distribution, t-distribution and Laplace distribution respectively. Formulas (6)-(12) give the mathematical expressions of each distribution function.

[0139]

[0140]

[0141]

[0142]

[0143]

[0144]

[0145]

[0146] In the formula, f A (.), fB (.), f C (.), f D (.), f E (.) and f F (.) represent Weibull distribution, lognormal distribution, Beta distribution, normal distribution, t-distribution and Laplace distribution respectively; k and λ represent the shape and scale parameters of Weibull distribution; μ2 and μ4 represent the mathematical expectations of lognormal distribution function and normal distribution function respectively; σ2 and σ4 represent the variances of lognormal distribution function and normal distribution function respectively; α and β are two parameters of Beta distribution; μ5, σ5 and v1 represent the location parameter, scale parameter and shape parameter of t-distribution respectively; μ6 is the location parameter of Laplace distribution, indicating the central position of the distribution, b is the scale parameter of Laplace distribution, and Γ(.) represents the gamma function.

[0147] (2) Optimization method

[0148] M sub-probability density distribution functions were selected using an optimization method. The steps are as follows:

[0149] Step 1): Solve for the parameters of a single sub-probability density distribution function. Use the least squares method to calculate the parameters of the f in the set. i (.) is fitted and the parameters of each probability density distribution function are obtained according to the measured data.

[0150] Step 2): Calculate the parameter error. Through statistical calculation, the actual mathematical expectation and variance of the actual calculation data can be obtained. Through the calculation in step 1, f i The mathematical expectation and variance of (.). The error between the statistical value and the parameter value obtained in step 1 can be calculated according to formula (13).

[0151]

[0152] In the formula, Var i (x) and E i (x) represents f i (.), Var(x) and E(x) represent the actual variance and mathematical expectation of the measured data.

[0153] Step 3): Screen the best sub-probability density distribution function. The smaller the error value, the better the f i (.) The more accurate the description of the probability density distribution of each admittance parameter is. M is the f with the minimum error i (.), which is used to form the general probability f(x).

[0154] 2) Model parameter identification based on improved mothfire optimized particle swarm algorithm

[0155] The model parameter identification process based on the improved mothfire optimized particle swarm algorithm is as follows: Figure 2 shown.

[0156] (1) Objective function

[0157] The objective function is constructed based on the difference in statistics. The smaller the difference between the mean and variance of the probability model and the mean and variance of the actual data, the better the model fitting effect. At the same time, the contribution of the mean and variance should be considered. In addition, the weight coefficient is also considered to measure the sub-probability density distribution function f i (.) contribution, the optimization parameters minimize the objective function.

[0158]

[0159] In the formula, They are the contribution of variance and mean to the objective function respectively, and their sizes can be determined by cross-validation method.

[0160] (2) Constraints

[0161] Due to the non-negativity and normalization of the probability density distribution function, the weight coefficients must satisfy the following constraints.

[0162]

[0163] λ i >0 (16)

[0164] (3) Solution method

[0165] The particle swarm algorithm based on improved mothfire optimization is a hybrid algorithm that combines the advantages of two algorithms to improve optimization performance. The mothfire optimization algorithm is used to set the initial position of the particle swarm, avoiding the situation of local optimal solution, and also improving the convergence speed. In this embodiment, M=3 is taken as an example. It is assumed that the selected probability distributions are Weibull distribution, Laplace distribution and lognormal distribution respectively; then the number of unknown parameters to be solved for fi(.) is 2, 2, 2 respectively; and there are 3 unknown weight coefficients. There are a total of 9 parameters to be identified. The present invention sets the output result of the mothfire optimization algorithm as the initial position information of the particle swarm, and continuously updates the initial position of the parameter to be identified through mothfire optimization, avoiding the selection of individuals with similar performance, and improving the search efficiency. The mothfire optimization model is as follows:

[0166] S(M u ,F v )=D uv e zw cos(2πw)+F v (17)

[0167] Where S(.) is the spiral function, Mu represents the u-th moth, F v represents the vth flame, D uv represents the distance D between the uth moth and the vth flame uv =|F v -M u |, z is a constant that defines the shape of the logarithmic spiral, w is a random number in the range [-1,1], and the coefficient w in the spiral function represents the distance between the moth's next position and the flame (w = -1 is the moth's closest position to the flame, while w = 1 shows the farthest position).

[0168] In order to balance the global search capability and local development capability of the algorithm in the search space, the number of flames is adaptively reduced during the iteration process. The formula is as follows:

[0169]

[0170] Where w0 represents the current iteration number, Q represents the maximum number of flames, W0 represents the maximum iteration number, and round represents the rounding function;

[0171] The optimal initial position of particles can be achieved through the mothfire optimization algorithm. The optimal position of particles can be obtained through continuous iteration. The position of each particle is the solution of the parameters to be identified. t d Represents the value of the tth solution at the dth iteration of the model.

[0172]

[0173] In the formula, is the jth parameter value of the tth solution at the dth iteration, where j is the number of unknown parameters;

[0174] Each particle has V t d Iterates in the solution space at a speed of and calculates the value of the objective function. t d gbest represents the optimal position of the t-th particle in the solution space at the d-th iteration. d It represents the optimal position of all particles in the solution space at the dth iteration. It is shown in the following formula:

[0175]

[0176]

[0177]

[0178] In the formula, v t,j dis the optimal speed in the solution space at the dth iteration; b t,j d represents the optimal position of the jth parameter of the tth solution as of the dth iteration. j d It represents the optimal position of the jth parameter among all solutions up to the dth iteration, that is, the value of the parameter to be identified.

[0179] The velocity and position of each particle after each iteration are expressed as (23)-(24).

[0180]

[0181]

[0182] In the formula, v t d represents the velocity of the t-th particle at the d-th iteration; r1 and r2 are random numbers in [0,1], Represents the solution of each parameter in the tth solution at the dth iteration of the model.

[0183] 3. Determination of the dominant coupling model based on typical power range

[0184] 1. Determination of typical power range

[0185] In order to improve the practicality of the harmonic source model, the power intervals are analyzed for similarity, and the power intervals with high similarity are clustered. Since the elements of the matrix obtained in each interval obey the probability distribution, the model similarity judgment of two power intervals needs to be considered from two aspects: element level and matrix level.

[0186] 1) Elemental level

[0187] As can be seen from the previous text, the parameter Y in the coupling admittance matrix is lm Obey the probability distribution, with the interval 10%~20%P n , 21%~30%P n For example, in order to measure the similarity of elements at the same position, the present invention applies the Bhattacharyya distance algorithm to calculate the overlap of elements at the same position under two power intervals. The higher the overlap of the waveforms of the two elements, the higher the similarity of their waveforms.

[0188]

[0189]

[0190] Where y1(x) and y2(x) are Y lm1 and Y lm2The proportion of sampling points in the corresponding power interval; BC is the Bhattacharyya coefficient, and 0 ≤ BC ≤ 1. The larger BC is, the higher the overlap degree of the two waveforms. W is the waveform similarity coefficient, and 0 < W ≤ 1. When the two waveforms are more similar, the value of W is larger, and vice versa. When W is greater than the threshold W set , it indicates that the similarity of the two waveforms at this position is relatively high. W set needs to be set artificially. The number and position information of elements higher than the threshold are stored in the set S0.

[0191] 2) At the matrix level

[0192] By traversing the entire matrix for element similarity calculation, the similarity of the same positions under two power intervals can be obtained. To quantify the similarity degree of the matrix, the present invention defines the matrix similarity factor δ.

[0193]

[0194] In the formula, S is the number of similar elements, and a0 and b0 are the number of rows and columns of the matrix respectively. 0 < δ ≤ 1. The value of δ affects the accuracy and calculation speed of the model. The larger the value of δ, the higher the accuracy of the model, and the lower the calculation efficiency, and vice versa.

[0195] When δ ≥ δ set , it is considered that the two matrices are similar. The actual value of δ set can be obtained by using machine learning algorithms, and the optimal value is obtained through training with machine algorithms. The present invention will not elaborate.

[0196] The present invention uses the Mean Shift algorithm to cluster the obtained similar intervals. Mean Shift is a density-based non-parametric clustering algorithm that does not require specifying the number of clusters in advance and does not assume the distribution of data. The core idea of Mean Shift is to determine the cluster center by finding the local maximum points of the density function. First, the calculated power interval model is divided into M n groups, and through continuous calculation and optimization using the Mean Shift algorithm, appropriate M n power states are selected as the cluster centers, and they are divided into different categories according to the distance between each power interval and the cluster centers, so as to minimize the distance from the power interval aggregation to the cluster centers. At this time, the optimal number of models based on the typical power intervals can be obtained.

[0197] Since the occurrence probabilities of each power interval are different, after aggregating the power interval models, the values of each parameter of the model need to be re-determined, and the value-taking rules are as follows:

[0198] Y lm ' = ∑Y lmx ×P K(28)

[0199] Where P K is the occurrence probability of the corresponding power interval.

[0200] 2. Screening of dominant elements in typical power range

[0201] Based on the above process, the probability model of the typical power range can be obtained. Since the harmonic admittance model to be solved is a full-rank matrix state, in order to improve the calculation efficiency, the obtained probability model is thresholded, that is, the obtained matrix model is reduced in order. For the hth harmonic current, its expression can be obtained as follows:

[0202]

[0203] Where Y h,n are all probability density distribution functions.

[0204] According to existing research, the self-coupling of harmonic current and harmonic voltage is dominant, that is, Y h,h The present invention selects the dominant element according to the probability density distribution of the self-coupling admittance.

[0205] First, determine the 95% of the cumulative distribution function (CDF) of the diagonal elements in the initial admittance matrix as the threshold. For each diagonal element, find the 95% point of its CDF. For each non-diagonal element in the matrix, calculate its CDF value in its own distribution. If this CDF value is lower than 95%, it is considered that the value of this element is lower than the threshold, and its probability density function value is set to zero. Apply the above process to the entire matrix, screen all non-diagonal elements and set the probability density distribution to zero. The mathematical expression is as follows:

[0206] Let Y be the original matrix and D be a diagonal matrix whose diagonal elements are the same as Y. For each off-diagonal element Y in the matrix Y lm , if its probability density f lm (Y lm ) corresponding to the CDF value F lm (Y lm ) is less than 95%, then do the following:

[0207]

[0208] Based on the probability distribution characteristics of the diagonal elements, other elements in the matrix are screened and processed to ensure that only those elements with higher probability density distribution values ​​are retained. In this way, the probabilistic coupling model of the dominant elements in the typical power range can be obtained to characterize the node harmonic emission level.

Claims

1. A harmonic source principal component probability modeling method based on typical power intervals, characterized in that: The following steps are involved: Step 1: Based on the harmonic coupling admittance matrix model, taking into account the randomness and stage characteristics of the collective load in the substation area, the power interval is divided according to the rated active power of the monitoring point in the substation area, and the probability of each power interval occurring within the monitoring time range is calculated. The initial matrix of each power interval is solved by the least squares method in the complex domain through the measured data; Step 2: Establish a probability model based on the divided power intervals, and use the improved mothfire optimized particle swarm algorithm to identify the parameters of the mixed probability density distribution function of the probability model to avoid local optimal solutions and improve the convergence speed; Step 3: Compare the similarities of the models in different power intervals. At the element level, use the Bhattacharyya distance algorithm to compare the probability distribution of elements at the same position. At the matrix level, define the matrix similarity factor to characterize the similarity of the admittance models in different power intervals, and use the mean shift algorithm to cluster the obtained similar probability intervals to obtain the optimal number of models based on the typical power interval. Step 4: Threshold the model parameters through the probability distribution of diagonal elements in the typical power range, screen out the dominant elements in the typical power range, and achieve effective characterization of the node harmonic emission level.

2. The harmonic source principal component probability modeling method based on typical power interval according to claim 1 is characterized in that: The step 1 specifically includes: Step 1.1: Determine the harmonic coupling admittance matrix model as a unified harmonic source model, the specific form is as follows: Where: are the phasor values ​​of the lth harmonic current and the mth harmonic voltage respectively; Characterizes the contribution of the mth harmonic voltage to the lth harmonic current; l = 1, 3, ..., P; m = 1, 3, ..., P, the highest harmonic is P; Step 1.2: Initial power interval model solution After dividing the power interval according to the rated active power of the monitoring point in the substation, the harmonic current after the collective equipment in the substation is connected is shown as follows: Where n is the number of power intervals; [I] K , [Y] K and [U] K They are the contribution of harmonic current and harmonic voltage to harmonic current and the phasor value of harmonic voltage in power interval K respectively; Step 1.3: Calculate the probability of each power interval occurring within the monitoring time range Assume that the occurrence time corresponding to power interval K is t K , then based on the measured data, the probability of occurrence of this power interval within the monitoring time is expressed as: In the formula, T is the total monitoring time; Step 1.4: Solve the initial matrix for each power interval For the admittance matrix of each power interval, the complex domain least squares method is used to solve the initial matrix. The calculation formula is as follows: In the formula, the superscript "*" indicates the transpose of the matrix, and the current matrix and voltage matrix Generated through monitoring point data.

3. The harmonic source principal component probability modeling method based on typical power interval according to claim 2 is characterized in that: The step 2 specifically includes: Step 2.1: Select a general model of mixed probability density distribution that characterizes multiple probability distributions: In the formula, f i (.) is the sub-probability density distribution function, λ i is the weight coefficient of the sub-probability density distribution function, indicating f i (.) The contribution degree in the model, x is the admittance data of the Kth power interval to be solved, M is the total number of selected sub-probability density distribution functions, that is, the M sub-probability density distribution functions with the smallest error are selected by the optimization method; Step 2.2: Model parameter identification based on improved mothfire optimized particle swarm algorithm The objective function is constructed by minimizing the difference between the mean and variance of the probability model and the mean and variance of the actual data, and the sub-probability density distribution function f is weighed by the weight coefficient i (.) contribution; the objective function is expressed as: In the formula, are the contribution of variance and mean to the objective function, respectively, and their sizes are determined by cross-validation method; Var i (x) and E i (x) represent the sub-probability density distribution function f i The variance and mathematical expectation of (.); The weight coefficient λ of the sub-probability density distribution function i The following constraints are met: λ i >0(8) The output result of the mothfire optimization particle swarm algorithm is set as the initial position information of the particle swarm, and the initial position of the parameter to be identified is continuously updated through the mothfire optimization to avoid individuals with similar performance from being selected to improve the search efficiency; the mothfire optimization model is as follows: S(M u ,F v )=D uv e zw cos(2πw)+F v (9) where S(.) is the spiral function, M u represents the u-th moth, F v represents the vth flame, D uv represents the distance D between the uth moth and the vth flame uv =|F v -M u |, z is a constant used to define the shape of the logarithmic spiral, w is a random number in the range [-1,1], and the coefficient w in the spiral function represents the distance between the moth's next position and the flame; In order to balance the global search capability and local development capability of the algorithm in the search space, the number of flames will be adaptively reduced during the iteration process. The formula is as follows: Where w0 represents the current iteration number, Q represents the maximum number of flames, W0 represents the maximum iteration number, and round represents the rounding function; The position of each particle is the solution of the parameters to be identified, S t d Represents the set of parameters of the tth solution at the dth iteration of the model: In the formula, is the jth parameter value of the tth solution at the dth iteration; Each particle has V t d Iterate in the solution space at a speed of and calculate the value of the objective function; In the formula, v t,j d is the optimal speed in the solution space at the dth iteration; pbest t d represents the optimal position of the t-th particle in the solution space at the d-th iteration; b t,j d gbest represents the optimal position of the jth parameter of the tth solution as of the dth iteration; d represents the optimal position of all particles in the solution space up to the dth iteration; g j d represents the optimal position of the jth parameter among all solutions up to the dth iteration, that is, the value of the parameter to be identified; The velocity and position of each particle after each iteration are expressed as: In the formula, v t d represents the velocity of the t-th particle at the d-th iteration; r1 and r2 are random numbers in [0,1], Represents the solution of each parameter in the tth solution at the dth iteration of the model.

4. The harmonic source principal component probability modeling method based on typical power interval according to claim 3 is characterized in that: The neutron probability density distribution function in step 2.1 includes: In the formula, f A (.), f B (.), f C (.), f D (.), f E (.) and f F (.) represent Weibull distribution, lognormal distribution, Beta distribution, normal distribution, t-distribution and Laplace distribution respectively; k and λ represent the shape and scale parameters of Weibull distribution; μ2 and μ4 represent the mathematical expectations of lognormal distribution function and normal distribution function respectively; σ2 and σ4 represent the variances of lognormal distribution function and normal distribution function respectively; α and β are two parameters of Beta distribution; μ5, σ5 and v1 represent the location parameter, scale parameter and shape parameter of t-distribution respectively; μ6 is the location parameter of Laplace distribution, indicating the central position of the distribution, b is the scale parameter of Laplace distribution, and Γ(.) represents the gamma function.

5. The harmonic source principal component probability modeling method based on typical power interval according to claim 3 is characterized in that: The optimization method used in step 2.1 to select M sub-probability density distribution functions with the minimum error is as follows: Step 2.1.1: Solve for the parameters of a single sub-probability density distribution function The least squares method is used to calculate the probability density distribution function f in the set. i (.) Fitting is performed to obtain the parameters of each probability density distribution function according to the measured data; Step 2.1.2: Calculate parameter errors Through statistical calculation, the actual mathematical expectation and variance of the actual calculation data are obtained; through the calculation in step 2.1.1, f is obtained i The mathematical expectation and variance of (.); The error between the statistical value and the parameter value obtained in step 2.1.1 is calculated as follows: In the formula, Var i (x) and E i (x) represents the sub-probability density distribution function f i (.), Var(x) and E(x) represent the actual variance and mathematical expectation of the measured data; Step 2.1.3: Screening the best sub-probability density distribution function Select M sub-probability density distribution functions f with the smallest error i (.), which is used to form the general probability f(x).

6. The harmonic source principal component probability modeling method based on typical power interval according to claim 3 is characterized in that: The step 3 specifically includes: Step 3.1: At the element level, use the Bhattacharyya distance algorithm to calculate the overlap of elements at the same position in two power intervals: Where y1(x) and y2(x) are Y lm1 and Y lm2 The proportion of sampling points in the corresponding power range; BC is the Bhattacharyya coefficient, and 0≤BC≤1; W is the waveform similarity coefficient, and 0 <W≤1;Y lm1 is the first power interval in the harmonic coupling admittance matrix Y lm2 is the second power interval in the harmonic coupling admittance matrix element; When the waveform similarity coefficient W is greater than the threshold W set When , it means that the similarity between the two waveforms at this position is high, which will be higher than the threshold W set The number and position information of the elements are stored in the set S0; Step 3.2: At the matrix level, define the matrix similarity factor δ to quantify the similarity of the matrices: In the formula, S is the number of similar elements, a0 and b0 are the number of rows and columns of the matrix respectively; The value of the matrix similarity factor δ affects the accuracy and calculation speed of the model. 0<δ≤1. The larger the δ value, the higher the accuracy of the model and the lower the calculation efficiency, and vice versa. When δ≥δ set When , the two matrices are considered similar; δ set is the similarity threshold; Step 3.3: Use the mean shift algorithm to cluster the obtained similar power intervals First, the calculated power interval model is divided into M n The Mean Shift algorithm is used to continuously calculate and optimize the group to select the appropriate M n The power states are taken as cluster centers, and each power interval is divided into different categories according to the distance between it and the cluster center, so as to minimize the distance to the cluster center after the power interval is aggregated, and the optimal number of models based on the typical power interval is obtained; After the power interval model is aggregated, the parameters of the model are revalued. The value selection rules are as follows: AND lm '=∑Y lmx ×P x (28) Where P x is the occurrence probability of the corresponding power interval; Y lm ' is the number of similar power intervals after clustering The value of the element again, Y lmx is the matrix of similar power intervals to be determined.

7. The harmonic source principal component probability modeling method based on typical power interval according to claim 3 is characterized in that: The step 4 specifically includes: Step 4.1: Threshold the obtained probability model, that is, reduce the order of the obtained harmonic coupling admittance matrix model; for the hth harmonic current, its expression is: Where Y h,n are all probability density distribution functions, is the phasor value of the Pth harmonic voltage; Step 4.2: Screening of dominant elements based on the probability density distribution of self-coupling admittance Determine the R of the cumulative distribution function of the diagonal elements in the initial admittance matrix n % as the threshold, for each diagonal element, find the R of its cumulative distribution function n % point, for each off-diagonal element in the matrix, calculate its cumulative distribution function value in its own distribution; if this cumulative distribution function value is lower than R n %, then the value of this element is considered to be lower than the threshold, and its probability density function value is set to zero; Apply the above process to the entire matrix, screen all off-diagonal elements and set the probability density to zero; the mathematical expression is as follows: Let Y be the original admittance matrix, D be a diagonal matrix whose diagonal elements are the same as the diagonal elements of the original admittance matrix Y; for each off-diagonal element Y in the original admittance matrix Y lm , if its probability density f lm (Y lm ) corresponds to the cumulative distribution function value F lm (Y lm ) is lower than the threshold R n %, then do the following: Thus, a probabilistic coupling model of the dominant elements in a typical power range is obtained, which enables the characterization of the node harmonic emission level.

Images