Harmonic power flow probability evaluation method and system based on skew gaussian mixture model

Abstract

The application discloses a harmonic power flow probability evaluation method and system based on a skew Gaussian mixture model, belongs to the field of power system power quality analysis, and comprises the following steps: obtaining active power time series of various harmonic sources according to historical data, respectively performing day cycle clustering and minute level clustering, and completing working condition division; establishing skew Gaussian mixture models corresponding to various working conditions, which are used for describing probability distribution of variable parameters in a three-branch parallel harmonic source model corresponding to the harmonic sources; when power flow evaluation is performed, at each evaluation time in an evaluation day, all possible joint working conditions and corresponding probabilities are determined according to the harmonic sources accessed; under each joint working condition, harmonic data of each node at the time is calculated by means of the corresponding skew Gaussian mixture model, corresponding harmonic evaluation indexes are calculated, the probability distribution characteristics of the harmonic evaluation indexes of the nodes are calculated, and harmonic power flow probability evaluation is completed. The application can improve the accuracy of harmonic power flow probability evaluation.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of power system power quality analysis, and more particularly, to a harmonic power flow probability evaluation method and system based on a skew Gaussian mixture model. BACKGROUND

[0002] With large-scale grid connection of renewable energy and wide application of power electronic equipment, the harmonic problem of power system is increasingly prominent and presents complex time-varying characteristics. Devices such as wind power generation, photovoltaic power generation, and electric vehicle charging piles inject rich harmonics into the power grid, intensifying the uncertainty of harmonic power flow in time and spatial distribution, and causing serious impact on the safe and stable operation of the power grid and power quality. Accurate evaluation of the probability distribution characteristics of time-varying harmonic power flow is an important prerequisite for realizing fine control of the safe and stable operation of the power grid and power quality under high proportion of renewable energy access.

[0003] To accurately evaluate the probability of harmonic power flow, the working condition needs to be divided first. The existing probability evaluation methods mainly consider the variation law of harmonic sources in the daily cycle scale, and adopt a single time scale to divide the working condition, which leads to rough working condition division, making it difficult to build an accurate segmented probability model, and further restricting the spatial resolution of the evaluation results. Secondly, due to the influence of factors such as load mutation and power fluctuation, the actual distribution of harmonic source model parameters often presents asymmetric and multi-peak characteristics. However, the existing methods mostly use normal distribution or simple mixed distribution to describe the harmonic parameters, which is difficult to accurately fit the complex statistical characteristics such as asymmetry and multi-peak state existing in the actual harmonic data, leading to distortion of the probability density function fitting, affecting the accuracy of the probability harmonic model. In addition, time-varying probability evaluation based on Monte Carlo simulation needs a large number of repeated calculations of deterministic harmonic power flow, and the traditional deterministic harmonic power flow algorithm has certain degree of convergence difficulty or precision deficiency problem, which is difficult to meet the online evaluation demand of large-scale new-type power system.

[0004] Overall, there is an urgent need for an efficient time-varying harmonic power flow probability evaluation method that can integrate multi-scale working condition division and complex distribution modeling to improve the quantitative accuracy of harmonic uncertainty in new-type power system. SUMMARY

[0005] In view of the defects and improvement needs of the prior art, the present application provides a harmonic power flow probability evaluation method and system based on a skew Gaussian mixture model, which aims to fine working condition division, improve harmonic data fitting accuracy, and thus improve the accuracy of harmonic power flow probability evaluation.

[0006] To achieve the above-mentioned purpose, according to one aspect of the present application, a harmonic power flow probability evaluation method based on a skew Gaussian mixture model is provided, comprising:

[0007] For each type of harmonic source, after sampling the historical voltage and current data of the harmonic source at a minute level sampling interval, the active power of each sampling point is calculated to obtain an active power time series; the active power time series is divided into multiple day power time series in a day cycle, and the multiple day power time series are clustered to obtain multiple day cycle clusters; the sampling points in each day cycle cluster are clustered respectively to obtain multiple minute-level clusters; the sampling points in the same minute-level cluster under the same day cycle cluster are divided into a working condition;

[0008] For each working condition, a corresponding skew Gaussian mixture model is established to describe the probability distribution of the variable parameters in the corresponding three-branch parallel harmonic source model of the harmonic source; the skew Gaussian mixture model is linearly combined by multiple skew Gaussian components; the variable parameters include the voltage source amplitude , impedance modulus and impedance phase ;

[0009] When performing power flow assessment, at each assessment time in the assessment day, all possible joint working conditions and corresponding probabilities are determined according to the accessed harmonic sources; the joint working condition is composed of a combination of working conditions of each harmonic source;

[0010] Under each joint working condition, the probability distribution of the variable parameters of the three-branch parallel harmonic source model is obtained based on the skew Gaussian mixture model corresponding to the working condition of each harmonic source, and is sampled to obtain the variable parameters of each three-branch parallel harmonic source model, the harmonic data of each node at this time is calculated, and the corresponding harmonic assessment index is calculated; the harmonic data includes: the voltage and current of each frequency of the node;

[0011] According to the probability of each joint working condition and the harmonic assessment index of each node under the corresponding joint working condition, the harmonic assessment index mean and standard deviation of each node at the assessment time is calculated to obtain the probability distribution characteristics of the harmonic assessment index of each node at the assessment time , and the harmonic power flow probability assessment is completed.

[0012] Further, for each working condition, a corresponding skew Gaussian mixture model is established, including:

[0013] A1: the voltage source amplitude , impedance modulus and impedance phase of each sampling point under the current working condition are calculated as the corresponding variable parameters of the corresponding sampling point;

[0014] A2: according to the preset number of skew Gaussian components K , the variable parameters of each sampling point are clustered to obtain K variable parameter clusters;K is a positive integer;

[0015] A3: initialize the weight of the first skew Gaussian component as the proportion of sampling points corresponding to the first variable parameter cluster, and initialize the mean, standard deviation and skew parameter of the first skew Gaussian component as the mean, standard deviation and skew parameter of the first variable parameter cluster respectively; wherein, k k k k k = 1, 2, …… K ;

[0016] A4: calculate the posterior probability of each sample of the variable parameter belonging to the first skew Gaussian component, and update the weight of the first skew Gaussian component as: x i k k ; N x denotes the total number of samples of the variable parameter;

[0017] A5: for each skew Gaussian component, update the values of the mean, standard deviation and skew parameter in order to minimize the negative log-likelihood function thereof;

[0018] A6: iteratively execute A4-A5 until a preset convergence condition is met;

[0019] A7: linearly combine each skew Gaussian component according to the corresponding weight to obtain a skew Gaussian mixture model corresponding to the current working condition.

[0020] Further, the determination manner of the number of skew Gaussian components includes: K

[0021] B1: initialize ;

[0022] B2: execute steps A1-A7 to obtain a skew Gaussian mixture model corresponding to the value of ;

[0023] B3: calculate a candidate score BIC( ) corresponding to the value of according to K' ;

[0024] B4: if​​​​​​​​​​​​​​​​​ Then If the value is increased by 1, proceed to B2; otherwise, proceed to B5.

[0025] B5: Select the candidate with the lowest score. K' Determined as the number of skewed Gaussian components K ;

[0026] in, K max This indicates the preset maximum number of skewed Gaussian components; For the present The value of corresponds to the maximum likelihood value of the skewed Gaussian mixture model.

[0027] Furthermore, the voltage source amplitude corresponding to each sampling point under the current operating condition. The calculation expression is:

[0028]

[0029] impedance magnitude and impedance phase The calculation methods include:

[0030] according to Calculate the impedance in a three-branch parallel harmonic source model ;

[0031] impedance Decompose it into magnitude and phase to obtain the impedance magnitude. and impedance phase ;

[0032] in, h Indicates the current harmonic order; H Indicates the highest harmonic order. n h This indicates the total frequency of the voltage being considered; V N Indicates the rated voltage of the point of common coupling. V 1 represents the fundamental voltage amplitude. V q express q Second harmonic voltage amplitude; and Representing common connection points h Subharmonic voltage vector and current phasor; This represents the phasor of the voltage source parameters, and its magnitude is... The phase is 0°; Indicates the parameters of the current source branch. This represents the admittance branch parameters.

[0033] Furthermore, the parameters of the current source branch and admittance branch parameters The calculation methods include:

[0034] Establish the following equation:

[0035]

[0036] The parameter matrix is ​​fitted using the complex partial least squares method. ;

[0037] Extracting the parameter matrix The first column of elements yields the parameters of the current source branch. Extracting the parameter matrix The admittance branch parameters are obtained from the main diagonal elements. ;

[0038] in, Indicates the first l sampling points h Second harmonic current Indicates the first l sampling points h Subharmonic voltage; h= 1,2…… H , l =1,2…… L , L This indicates the number of sampling points selected.

[0039] Furthermore, the calculation methods for harmonic data at each node include:

[0040] according to The node where the harmonic source is connected in the computing system i The h Subharmonic current ;

[0041] Calculate the first value of each node in the system according to the following equation. h Subharmonic voltage:

[0042]

[0043] The system contains a total of n There are nodes, of which the first m Each node is a node where a harmonic source is connected; Represents a node i Accessed harmonic sources h Secondary equivalent admittance; and Representing nodes respectively i The parameters of the current source branch and the admittance branch in the three-branch parallel harmonic source model corresponding to the harmonic source are connected. and Representing nodes respectively i The impedance and voltage sources in the three-branch parallel harmonic source model corresponding to the harmonic source are connected; Represents a node i Accessed harmonic sources h Secondary equivalent current source Represents a node i and nodes j Between h Secondary network admittance i =1,2…… n , j =1,2…… n ; and Representing nodes respectively i of h Secondary harmonic voltage and current.

[0044] Furthermore, the sampling points in each daily periodic cluster were clustered separately to obtain multiple minute-level clusters, including:

[0045] E1: Calculate the mean of the daily power time series of each daily period cluster as a typical daily power curve;

[0046] E2: For each typical daily power curve, detect the points of change and use the detected points of change to divide the typical daily power curve into multiple intervals;

[0047] E3: Extract the time-domain statistical features, waveform morphology features, trend features, and fluctuation features of each interval based on the sampling points contained within the interval, and combine them into the interval features of the corresponding interval.

[0048] E4: Cluster each interval based on its characteristics to obtain multiple minute-level clusters.

[0049] Furthermore, for each typical daily power curve, before detecting the points of change, the following steps are performed to determine the optimal number of points of change. and the corresponding change point positions :

[0050] F1: Typical daily power curve Mapping to the reproducing kernel Hilbert space via Gaussian radial basis functions , to obtain the feature sequence ; T This indicates the number of sampling points included in the daily power time series. The first in a typical daily power curve t One sampling point, Indicates the mapping Corresponding feature points;

[0051] F2: Set the number of change points ;

[0052] F3: Initialization B The position of each change point and set , ;

[0053] F4: For each interval obtained by dividing the points of change, extract the time-domain statistical features, waveform morphology features, trend features and fluctuation features of each interval based on the sampling points contained in the interval, and combine them into the interval features of the corresponding interval.

[0054] F5: Construct the following cost function:

[0055]

[0056] in, Representing an interval The mean of the characteristic sequence, b =1,2,…… B ; Represents the regenerated nucleus Hilbert space norm, ;

[0057] F6: If B If the number of change points is less than the preset maximum number, then... B If the value is increased by 1, proceed to F3; otherwise, proceed to F7.

[0058] F7: Minimizes the cost function B Value as the optimal number of change points and the B The location of the change point corresponding to the value is determined as the optimal number of change points. The corresponding change point location.

[0059] Furthermore, when clustering the multiple daily power time series obtained from the division, and when clustering each interval in step E3, the improved elbow rule is used to determine the optimal number of clusters; the improved elbow rule includes:

[0060] Within the candidate clustering range Traverse within the range, for each r After clustering the samples to be clustered using this as the clustering quantity, according to... Calculate the inertia value and construct an inertia value curve that shows the change of inertia value with the number of clusters. Determine the number of clusters corresponding to the point with the largest rate of change of curvature as the optimal number of clusters.

[0061] Cluster the samples to be clustered according to the determined optimal number of clusters;

[0062] When clustering multiple daily power time series obtained from the division, the samples to be clustered are daily power time series; when clustering each interval, the samples to be clustered are interval features. For the first c A sample set of clusters, Indicates the first c The centroid sequence of each cluster; Representing the i One sample; DTW() represents the dynamic time warp distance; r max This indicates the preset maximum number of clusters.

[0063] According to another aspect of the present invention, a harmonic power flow probability assessment system based on a skewed Gaussian mixture model is provided, comprising:

[0064] A computer-readable storage medium for storing computer programs;

[0065] And a processor for reading a computer program stored in a computer-readable storage medium and executing the above-described harmonic power flow probability assessment method based on a skewed Gaussian mixture model provided by the present invention.

[0066] In summary, the above-described technical solutions conceived in this invention can achieve the following beneficial effects:

[0067] (1) This invention divides the operating conditions of harmonic sources into multiple scales by combining daily cycle time scale and minute-level dynamic time scale. It takes into account the overall operating law of harmonic sources in a long period and captures their dynamic fluctuation characteristics in a short period of time, which effectively improves the accuracy of operating condition division and provides reliable basic data support for constructing a segmented skewed Gaussian mixture model of harmonic sources. It can effectively improve the accuracy of harmonic power flow probability assessment.

[0068] (2) This invention utilizes a skewed Gaussian mixture model to describe the probability distribution of variable parameters in the three-branch parallel harmonic source model corresponding to each harmonic source. This skewed Gaussian mixture model is composed of multiple skewed Gaussian components in a linear combination, which can effectively fit the probability density function with complex distribution characteristics such as asymmetry and multimodality, and is more in line with the statistical law of actual harmonic parameters, thereby improving the accuracy of harmonic source probability modeling. In a further preferred embodiment of this invention, the optimal number of skewed Gaussian components in the model is determined based on the maximum likelihood value of the model corresponding to different numbers of skewed Gaussian components, which can further improve the accuracy of model fitting.

[0069] (3) Based on the three-branch parallel harmonic source model, the present invention derives the power flow equation between harmonic voltage and harmonic current at each node of the system. The relevant parameters can be accurately calculated without iteration, thereby realizing non-iterative deterministic harmonic power flow calculation, avoiding convergence problems and improving calculation speed. Attached Figure Description

[0070] Figure 1 For the existing three-branch parallel harmonic source model h Schematic diagram of the single-phase equivalent circuit for subharmonics.

[0071] Figure 2 The overall implementation flowchart of a harmonic power flow probability assessment method based on a piecewise skewed Gaussian mixture model provided by the present invention is shown below.

[0072] Figure 3 This is a schematic diagram of a typical daily power curve of a wind farm in an embodiment of the present invention.

[0073] Figure 4 This is a schematic diagram of a typical daily power curve of a photovoltaic power station in an embodiment of the present invention.

[0074] Figure 5 This is the result of the minute-level dynamic timescale division of the typical daily power curve of the wind farm in the embodiments of the present invention.

[0075] Figure 6 This is the result of the minute-level dynamic timescale division of the typical daily power curve of a photovoltaic power station in this embodiment of the invention.

[0076] Figure 7 This is a schematic diagram of the construction of the probabilistic harmonic model in an embodiment of the present invention.

[0077] Figure 8 Wind farm in an embodiment of the present invention 5th harmonic parameters under operating conditions , and The histogram and the probability density function of the corresponding skewed Gaussian mixture model.

[0078] Figure 9 This is the modified IEEE 33-node system topology in this embodiment of the invention.

[0079] Figure 10 The mean and standard deviation of the fifth harmonic voltage amplitude of node 18 in this embodiment of the invention are time-varying results within the evaluation day.

[0080] Figure 11 The mean and standard deviation of the total harmonic distortion rate of the voltage at node 18 in this embodiment of the invention are time-varying results within the evaluation day. Detailed Implementation

[0081] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0082] In this invention, the terms "first," "second," etc. (if present) in the invention and the accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence.

[0083] To effectively improve the accuracy of harmonic power flow probability assessment, this invention provides a method and system for harmonic power flow probability assessment based on a skewed Gaussian mixture model. First, it utilizes harmonic source voltage and current measurement data to classify harmonic source operating conditions using a combination of daily periodic timescales and minute-level dynamic timescales. Second, based on the operating condition classification results, it constructs a probabilistic harmonic model corresponding to each operating condition using a skewed Gaussian mixture distribution. Finally, it completes the harmonic power flow probability assessment based on the established skewed Gaussian mixture model.

[0084] Harmonic power flow probability assessment involves calculations related to harmonic sources. In this invention, a three-branch parallel harmonic source model is used to complete the relevant calculations. For a description of the three-branch parallel harmonic source model, please refer to the paper "A three parallelbranches harmonic model with adaptive capability of operating state variation" (ISSN: 0142-0615). Figure 1 The image shows a model of a three-branch parallel harmonic source. h The physical meaning of each branch in the subharmonic single-phase equivalent circuit is as follows: Current source branch Characterizing the fundamental voltage at the steady-state operating point caused by h Secondary current component; admittance branch Characterizing voltage caused by the same frequency h Secondary current component; voltage source and impedance series branch - Characterizing the changes in non-co-frequency harmonic voltages, harmonic source characteristic parameters relative to the steady-state operating point, and first-order Peano remainder pairs. h The combined effects of secondary current.

[0085] Example 1:

[0086] A method for assessing the probability of harmonic power flow based on a skewed Gaussian mixture model, such as Figure 2 As shown, it includes the following steps:

[0087] S1, acquire historical voltage and current data of the harmonic source, and classify the operating conditions of the harmonic source by comprehensively considering the daily periodic time scale and the minute-level dynamic time scale, specifically:

[0088] S1.1 Harmonic Measurement Data Acquisition. This embodiment uses wind farms and photovoltaic power plants as typical harmonic sources, deploying power quality monitoring devices at their point of common connection (PCC) with the grid. Voltage and current waveform data for phases A, B, and C are simultaneously acquired at 3-minute sampling intervals. Fast Fourier Transform (FFT) is used to extract voltage and current measurement data for each frequency. This embodiment collects 29 days of harmonic measurement data from wind farms and photovoltaic power plants, selecting phase A as the target analysis phase. By adjusting the target analysis phase, this method can be implemented equivalently. The measurement data includes:

[0089] (1) Fundamental component: voltage amplitude Current amplitude and its phase angle , ;

[0090] (2) Harmonic components: H Voltage amplitude of harmonics within the first order Current amplitude and its phase angle , ;

[0091] (3) Timestamp information: Time stamp information corresponding to each sampling point , , N This represents the total number of sampling points.

[0092] S1.2, Calculation of active power of harmonic sources and construction of time series. The active power of harmonic sources is calculated based on the fundamental frequency and harmonic components, and its expression is as follows:

[0093]

[0094] In the formula, H In this embodiment, the highest harmonic order is considered during the calculation process. H =25, in other embodiments H Other values ​​can be set; , They represent Phase A at time h RMS values ​​of second harmonic voltage and current , The corresponding phase angle is given. The time-power two-dimensional matrix is ​​constructed as follows:

[0095]

[0096] Using sampling time as an index, a time series of active power of harmonic sources is constructed within the statistical period. .

[0097] S1.3, for active power time series The three-level data preprocessing is performed sequentially, with the specific steps as follows:

[0098] S1.3.1, Abnormal Data Detection and Repair. The sliding window method is used to calculate the active power time series. rolling mean with standard deviation Construct dynamic threshold range .Will The data points were identified as outliers and corrected using the forward padding method, with the following correction rules:

[0099]

[0100] S1.3.2, Smoothing Filtering. A polynomial fitting filter is used to smooth the repaired active power time series. Smoothing is performed to obtain a smoothed sequence. To address the potential negative power values ​​after smoothing filtering, a non-negativity constraint is implemented, with the following rules:

[0101]

[0102] in, Indicates to Active power after applying nonnegativity constraints;

[0103] S1.3.3, Normalization processing. Active power is normalized through extreme value normalization. Mapping to the [0, 1] interval yields the normalized active power time series. .

[0104] S1.4, Normalized active power time series from a daily periodic time scale Clustering is performed to obtain typical daily power curves of harmonic sources. The specific steps are as follows:

[0105] S1.4.1, Daily Dimension Data Reconstruction. Normalized active power time series data is reconstructed with a daily period. Partitioning and constructing a three-dimensional tensor .in, (Number of days) (Number of sampling points per day) (Power characteristics only)

[0106] S1.4.2, the improved elbow rule is used to determine the number of typical daily power curves. (This refers to the number of candidate clusters.) Traversing within the range, a K-means clustering model is constructed using Dynamic Time Warping (DTW) as the distance metric, and an inertia value function is defined:

[0107]

[0108] In the formula, For the first c A set of active power time series data for each cluster; Indicates the first c The centroid sequence of each cluster; Representing the i Given a time series sample, such as the active power curve for a specific day, calculate the number of different candidate clusters. r Corresponding inertia value The inertia value curve is obtained, and the optimal number of clusters, i.e. the number of typical daily power curves, is determined by the point with the maximum rate of curvature change in the inertia value curve.

[0109] S1.4.3, perform time series clustering based on the number of typical daily power curves, and calculate the probability corresponding to each typical daily power curve. The K-means model based on the DTW metric is used for the three-dimensional tensor. Perform clustering prediction to obtain the cluster label to which each daily power time series belongs. Statistics on the labels of typical solar clusters. The corresponding number of sample days is used as the proportion of the total number of daily power time series samples to determine the corresponding probability. In the original active power time series... In the process, a "daily periodic cluster label" attribute is added to each sampling point, the value of which is the cluster label corresponding to the date of that sampling point. .

[0110] S1.4.4, Extraction of typical daily power curves for harmonic sources. The mean of the daily power time series within each cluster is calculated as the typical daily power curve. In this embodiment, the typical daily power curves for wind farms and photovoltaic power plants are as follows: Figure 3 and Figure 4 As shown in the figure, wind farms include three typical daily power curves with probabilities of 0.5862, 0.2069, and 0.2069, respectively; photovoltaic power plants include three typical daily power curves with probabilities of 0.1379, 0.5172, and 0.3449, respectively.

[0111] S1.5, based on typical daily power curves, classifies the operating conditions of harmonic sources from a minute-level dynamic time scale. The specific steps are as follows:

[0112] S1.5.1, typical daily power curve Mapped to the reproducing kernel Hilbert space via Gaussian radial basis functions (RBF). , to obtain the feature sequence The RBF kernel function is... Bandwidth parameters Determined using the median heuristic: The cost function is constructed by minimizing the sum of squared distances between the piecewise signals and the segment means within the kernel space:

[0113]

[0114] In the formula, The location of the change point ( , ); Representing the b The number is an interval The mean of the feature sequence. The maximum number of variation points in the feature sequence is set to... , in turn B equal to 1 to The optimal point of change is found using dynamic programming:

[0115]

[0116] Plot a line graph of the cost function and the number of points of change, and record the number of points of change when the cost function tends to stabilize. and corresponding positions .

[0117] S1.5.2, Typical daily power curve interval feature extraction. Based on the detected... The typical daily power curve is divided into several points of change. This embodiment extracts the time-domain statistical features, waveform morphology features, trend features, and fluctuation features of each interval, and concatenates them into the corresponding interval feature vector. Simultaneously extracting these features from each interval eliminates redundant change points and achieves more accurate operating condition division. The time-domain statistical features include mean, standard deviation, median, skewness, and kurtosis; the waveform morphology features include peak-to-peak value, number of peaks, and waveform area; the trend feature is the slope of the first-order linear fit of the typical daily power curve within the corresponding interval; and the fluctuation features include the mean of the absolute values ​​of the first-order differences between adjacent points of the typical daily power curve within the corresponding interval and the proportion of positive differences.

[0118] S1.5.3, Construction of the Feature Matrix of the Typical Daily Power Curve and Determination of the Optimal Operating Conditions. The extracted feature vectors from each interval are arranged row-wise to construct the feature matrix of the typical daily power curve, and Z-score standardization is performed. The standardized feature matrix is ​​clustered using a K-means model, and the inertia values ​​for different numbers of clusters are recorded. Based on the inertia values ​​corresponding to different numbers of clusters, inertia value variation curves are plotted, and the cluster number corresponding to the point with the maximum rate of curvature change is selected as the optimal operating condition number for the typical daily power curve. .

[0119] S1.5.4, based on the K-means model, divides the operating conditions of harmonic sources into minute-level dynamic time scales. Under the determined optimal operating condition number... Below, the K-means model is used to cluster the feature matrix to obtain the cluster labels corresponding to each interval of the typical daily power curve. In the original active power time series In the process, a "minute-level cluster label" attribute is added to each sampling point, the value of which is the cluster label corresponding to the sampling time of that point. In this embodiment, the minute-level operating condition division results of typical daily power curves for wind farms and photovoltaic power plants are as follows: Figure 5 and Figure 6 As shown in the figure, the dashed lines represent the locations of the points of change. Figure 5 The display shows a typical daily power curve of a wind farm. Divided into 5 operating conditions ( ); Divided into 6 operating conditions ( ); Divided into 5 operating conditions ( ). Figure 6 The display shows a typical daily power curve of a photovoltaic power plant. Divided into 6 operating conditions ( ); Divided into 5 operating conditions ( ); Divided into 5 operating conditions ( ).

[0120] S2, based on the daily cycle and minute-level multi-scale operating condition division results, establishes a piecewise skewed Gaussian mixture model of the harmonic source for each typical day, specifically:

[0121] S2.1, Multi-scale operating condition data grouping. Clustering is performed based on the daily periodic clustering labels corresponding to the sampling points. and minute-level clustering tags The voltage and current data at different frequencies are grouped to obtain harmonic data for different operating conditions. In this embodiment, the voltage and current data at different frequencies of the wind farm are divided into 16 groups, corresponding to 16 operating conditions; the voltage and current data at different frequencies of the photovoltaic power station are also divided into 16 groups, corresponding to 16 operating conditions.

[0122] S2.2 Calculate the parameters of the three-branch parallel harmonic source model under various operating conditions. The harmonic source model consists of a voltage source and an impedance series branch under each harmonic. - Current source branch Admittance branch Parallel connection constitutes its h The equation for the second harmonic current is:

[0123]

[0124] In the formula, and For PCC points h Second harmonic voltage and current phasors. Based on a specific operating condition. L ( L (Indicates the number of selected sampling points) Group PCC point harmonic voltage and current measurement data ( , , Establish the following equation:

[0125]

[0126] The parameter matrix is ​​fitted using the complex partial least squares method. The parameters of the current source branch can be obtained by extracting its first column and main diagonal elements. and admittance branch parameters Voltage source branch parameters The phase is set to 0°, and its amplitude is set to the average of the fundamental voltage fluctuation and the non-harmonic voltages of the same frequency:

[0127]

[0128] In the formula, For the total frequency of voltage to be considered, This is the rated voltage at point PCC. Impedance branch parameters. Derived from Kirchhoff's laws:

[0129]

[0130] S2.3, Establish a piecewise skewed Gaussian mixture model of the harmonic source on each typical day. The specific steps are as follows:

[0131] S2.3.1, Determine the variable parameters of the probabilistic harmonic model. Under operating conditions... Below is a schematic diagram of the probabilistic harmonic model construction. Figure 7 As shown. (Regarding operating conditions) Current source parameters and admittance parameters The parameters are constant, while the voltage source parameters are constant. and impedance parameters It changes dynamically with the operating state of the harmonic source. (Through decomposition) and In terms of magnitude and phase, we can obtain the following respectively: and In the parameter calculation of S2.2, It is set to 0° and therefore not included in the range of variable parameters. Thus, under specific operating conditions, the variable parameter of the probabilistic harmonic model is the voltage source amplitude. Impedance modulus and impedance phase .

[0132] S2.3.2, Construction of the skewed Gaussian mixture model. Let the variable parameters be... The probability distribution is given by K The linear combination of skewed Gaussian components represents the probability density function as follows:

[0133]

[0134] in, Represents the skewed Gaussian components; For the first k The weights of the skewed Gaussian components satisfy the following conditions: ,and ; and These are the probability density function and cumulative distribution function of the standard normal distribution, respectively; , , The first k The mean, standard deviation, and skewness parameters of each component.

[0135] S2.3.3, Initialize the parameters of the skewed Gaussian mixture model. The K-means algorithm is used to initialize the variable parameters. x The samples are clustered according to the preset number of skewed Gaussian components. K Generate cluster labels. Initialize model parameters based on the statistical characteristics of samples within each cluster: from the first... k Initial weights based on class sample proportion The mean is initialized using the mean, standard deviation, and skewness of the within-class samples. Standard deviation and skewness parameters .

[0136] S2.3.4, Calculate the variable parameters Each sample Belongs to the k The posterior probability of each component:

[0137]

[0138] S2.3.5, Update the weights of each skewed Gaussian component:

[0139]

[0140] In the formula, For variable parameters x The total number of samples. For each skewed Gaussian component, the L-BFGS-B numerical optimization algorithm is used to minimize the negative log-likelihood function, and the calculation is performed. , and Parameter estimates:

[0141]

[0142] S2.3.6, iterate through S2.3.4 and S2.3.5 until the convergence condition is met. Optional convergence conditions are: the change in log-likelihood is less than the threshold 1e-4 or the maximum number of iterations (100) is reached.

[0143] S2.3.7, Determining the optimal number of components. The optimal number of components is automatically selected based on the Bayesian information criterion. :

[0144]

[0145] In the formula, L () represents the maximum likelihood value of the model. Optionally, variable parameters can be set. x The probability density function consists of at most 5 skewed Gaussian components. The number of candidate components is then traversed. Then execute steps S2.3.3 to S2.3.7 sequentially, selecting the option that minimizes BIC. Value as the number of optimal components The value of .

[0146] S2.3.8, Construction of the piecewise skewed Gaussian mixture model of the harmonic source. Operating conditions... Lower variable parameter The sample input skewed Gaussian mixture model is used, and parameter estimation and component number optimization are performed through S2.3.3~S2.3.7 to calculate the harmonic source in a typical day. Piecewise probabilistic harmonic model:

[0147]

[0148] In the formula, Indicates the first y Variable parameters under minute-level operating conditions x The skewed Gaussian mixture probability distribution.

[0149] Based on wind farm operating conditions Taking the 5th harmonic as an example, Figure 8 The parameters were displayed. , and The histogram and its corresponding probability density function of the skewed Gaussian mixture model (obtained from S2.3.3~S2.3.7). Among them, and The optimal number of components is 2. The optimal number of components is 3. The specific parameters (weights) of the skewed Gaussian mixture model are as follows: mean Standard deviation and skewness parameters As shown in Table 1.

[0150] Table 1 Parameters , and Parameters of skewed Gaussian mixture model

[0151]

[0152] S3 uses a combined enumeration method to determine the joint operating conditions and corresponding probabilities at each time point, and performs time-varying harmonic power flow probability assessment based on a piecewise skewed Gaussian mixture model, specifically:

[0153] S3.1 System Parameter Initialization and Timing Settings. Identify the harmonic source access nodes in the system (e.g., wind power, photovoltaic access nodes). Set the evaluation time point sequence within the evaluation day. ,in T To assess the total number of time points, this embodiment modifies the IEEE 33-node system: photovoltaic power is integrated into nodes 17 and 18, and wind power is integrated into nodes 32 and 33. The modified IEEE 33-node system topology is as follows: Figure 9 As shown. The set evaluation time sequence is the hourly times within the evaluation day, i.e. , T =24.

[0154] S3.2, determine the combined operating conditions and their probabilities at the evaluation time. For the evaluation time... Extract the tags of all harmonic sources connected at that moment. Combined into system joint operating conditions The combined enumeration method is used to count... All combined operating conditions And calculate the probability of its occurrence. (satisfy ).

[0155] S3.3, Time-varying operating condition matching and harmonic source model parameter extraction. For combined operating conditions... Based on the working condition labels it contains The constant parameters of the harmonic source model are extracted from the results of S2.2: current source parameters. and admittance parameters In addition, according to the operating condition label The piecewise skewed Gaussian mixture model established from S2.3.8 The skewed Gaussian mixture distribution parameters of the variable parameters are extracted. The variable parameters are the voltage source amplitude. Impedance modulus and impedance phase The parameters of a skewed Gaussian mixture distribution include the weights of each component. mean Standard deviation and skewness parameters ( , (The optimal number of components).

[0156] S3.4, Based on the probability distribution and correlation of the variable parameters of the harmonic source, a random sample is generated using improved Latin hypercube sampling, and a random sample matrix with a specified correlation is obtained through the principle of equal probability transformation. The specific steps are as follows:

[0157] S3.4.1, Read target system parameters. Read system network parameters (at least including bus, line, transformer, and generator parameters) and the harmonic sources extracted in S3.3 connected to the combined operating condition. The harmonic source model parameters (probability distribution parameters of constant and variable parameters) are defined. Based on historical data or engineering experience, target correlation matrices are set between variable parameters of different harmonic source connections and between different parameters of the same node. .

[0158] S3.4.2 generates approximately orthogonal independent standard normal samples. An improved Latin hypercube sampling method is employed, iteratively optimizing the method to minimize the correlation between variables, resulting in samples with the following dimensions. ( z For the sample size, m Independent uniform sample matrix (number of nodes connected to harmonic sources) .matrix Each column is uniformly distributed in the interval [0,1], and the columns are approximately orthogonal to each other (correlation coefficient close to 0). For the matrix... Applying the inverse cumulative distribution function of the standard normal distribution to each column yields an independent standard normal distribution sample matrix. .

[0159] S3.4.3 generates normally distributed samples with specified correlations. This involves generating the target correlation matrix from the original random variable space (where each parameter follows a skewed Gaussian mixture distribution). Transform to the standard normal space, and denote the correlation matrix after transformation as follows: .right Perform Cholesky decomposition to obtain the lower triangular matrix. Using matrices Independent standard normal samples Transform into a normal sample with a specified correlation. .in, .

[0160] S3.4.4, equal probability transformation mapping back to the original variable space. Based on the cumulative distribution function of the standard normal distribution. The inverse cumulative distribution function of the skewed Gaussian mixture distribution of each variable parameter extracted in S3.3 , related normal samples Map back to the original random variable space to generate a random parameter sample matrix with specified correlation. :

[0161]

[0162] In the formula, the matrix Each row represents a complete random scenario sample, containing all harmonic sources connected in the current operating condition. Sampled values ​​of all variable parameters.

[0163] S3.5, for the random sample matrix For each sample, fundamental wave power flow calculation and deterministic non-iterative harmonic power flow calculation are performed sequentially, with the specific steps as follows:

[0164] S3.5.1, perform fundamental power flow calculation to obtain the fundamental voltage amplitude and phase of all nodes.

[0165] S3.5.2, for the harmonic order under consideration Deterministic non-iterative harmonic power flow calculations are performed sequentially.

[0166] For each node connected to the harmonic source i ,calculate h Subharmonic current :

[0167]

[0168] The system contains n There are nodes, of which the first m The nodes are those where harmonic sources are connected. For ease of description, the corresponding nodes are referred to by their node numbers. The nth node of each node in the system is calculated according to the following equation. h Subharmonic voltage:

[0169]

[0170] in, and Representing nodes respectively i The parameters of the current source branch and the admittance branch in the three-branch parallel harmonic source model corresponding to the harmonic source are connected. and Representing nodes respectively i The impedance and voltage sources in the three-branch parallel harmonic source model corresponding to the harmonic source are connected; Represents a node i Accessed harmonic sources h Secondary equivalent admittance Represents a node i Accessed harmonic sources h Secondary equivalent current source; Represents a node i and nodes j Between h The secondary network admittance can be calculated based on the harmonic characteristics of the generator, load, and line. i =1,2…… n , j =1,2…… n ; and Representing nodes respectively i of h Secondary harmonic voltage and current.

[0171] S3.5.3, Based on the solved harmonic current and harmonic voltage, solve for the corresponding harmonic evaluation index. Optionally, in this embodiment, the harmonic evaluation index specifically includes the node harmonic voltage amplitude and the total harmonic distortion rate of the node voltage. Specifically, the total harmonic distortion rate of the node voltage can be calculated based on the fundamental voltage and harmonic voltage of each node. It should be noted that in practical applications, the selected harmonic evaluation index can be flexibly set. In some other embodiments of the present invention, it can also be set to include the total harmonic loss of the system. Specifically, the total harmonic loss of the system can be calculated by summing the harmonic losses of each branch.

[0172] S3.6, for the evaluation time All combined operating conditions Then execute steps S3.3 to S3.5 in sequence.

[0173] S3.7, for all evaluation moments Then execute steps S3.2 through S3.6 in sequence.

[0174] S3.8, using the law of total probability to evaluate the time step. The overall probability distribution characteristics of the system's harmonic indices are calculated by performing a probabilistic weighted synthesis of all combined operating conditions. The specific steps are as follows:

[0175] S3.8.1, Probabilistic weighted calculation of harmonic indices at a single time step. Determine the system at time step [missing information]. Existing C Combined working conditions and their corresponding probabilities are For each harmonic evaluation index Calculate at time Weighted mean and weighted second-order raw moments:

[0176]

[0177]

[0178] In the formula, and They represent the operating conditions respectively. Lower harmonic evaluation index The sample mean and sample second-order raw moments. Optionally, harmonic indices include nodal harmonic voltage amplitudes. and voltage total harmonic distortion Harmonic evaluation indices are calculated based on the sample mean and the sample second-order raw moment. At any moment The overall standard deviation is as follows:

[0179]

[0180] S3.8.2, Comprehensive Time-Varying Characteristics Across All Time Periods. This applies to all evaluation times. Repeat step S3.8.1. Arrange the weighted mean and overall standard deviation of each harmonic assessment index calculated at each time point, along with other probability statistics, in a time series as the probability assessment result of the time-varying harmonic power flow. Use S3 to calculate the time-varying harmonic power flow. Figure 10 and Figure 11The time-varying results of the fifth harmonic voltage amplitude, the mean and standard deviation of the total harmonic distortion rate of voltage at representative node 18 in the system during the evaluation day are shown and compared with the Monte Carlo simulation. The comparison results show that the harmonic power flow probability assessment method provided in this embodiment can accurately perform time-varying probability harmonic power flow calculation.

[0181] Example 2:

[0182] A harmonic power flow probability assessment system based on a skewed Gaussian mixture model includes:

[0183] A computer-readable storage medium for storing computer programs;

[0184] And a processor for reading a computer program stored in a computer-readable storage medium and executing the harmonic power flow probability assessment method based on the skewed Gaussian mixture model provided in Embodiment 1 above.

[0185] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for assessing the probability of harmonic power flow based on a skewed Gaussian mixture model, characterized in that, include: For each type of harmonic source, the historical voltage and current data of the harmonic source are sampled at minute-level sampling intervals. The active power at each sampling point is calculated to obtain the active power time series. The active power time series is divided into daily periods, and the multiple daily power time series are clustered to obtain multiple daily period clusters. The sampling points in each daily period cluster are clustered to obtain multiple minute-level clusters. The sampling points in the same minute-level cluster under the same daily period cluster are classified into one operating condition. For each operating condition, a corresponding skewed Gaussian mixture model is established to describe the probability distribution of the variable parameters in the three-branch parallel harmonic source model corresponding to the harmonic source; the skewed Gaussian mixture model is composed of a linear combination of multiple skewed Gaussian components; the variable parameters include the voltage source amplitude. Impedance modulus and impedance phase ; When conducting a current flow assessment, at each assessment moment within the assessment day Based on the accessed harmonic sources, all possible combined operating conditions and their corresponding probabilities are determined; the combined operating conditions are composed of one operating condition of each harmonic source. Under each combined operating condition, the probability distribution of the variable parameters of the corresponding three-branch parallel harmonic source model is obtained based on the skewed Gaussian mixture model under the corresponding operating condition of each harmonic source, and the variable parameters of each three-branch parallel harmonic source model are obtained by sampling. The harmonic data of each node are calculated at this time, and the corresponding harmonic evaluation index is calculated. The harmonic data includes: the voltage and current of each frequency of the node. Based on the probability of each joint operating condition and the harmonic evaluation index of each node under the corresponding joint operating condition, the values ​​of each node at the evaluation time are calculated. The mean and standard deviation of the harmonic evaluation index are used to obtain the values ​​of each node at the evaluation time. Based on the probability distribution characteristics of the lower harmonic evaluation index, complete the probability assessment of harmonic power flow. The three-branch parallel harmonic source model includes: current source branch. Admittance branch And voltage source and impedance series branch - Current source branch Characterizing the fundamental voltage at the steady-state operating point caused by h Secondary current component; admittance branch Characterizing voltage caused by the same frequency h Secondary current component; voltage source and impedance series branch - Characterizing the changes in non-co-frequency harmonic voltages, harmonic source characteristic parameters relative to the steady-state operating point, and first-order Peano remainder pairs. h The combined effects of secondary current.

2. The harmonic power flow probability assessment method based on a skewed Gaussian mixture model as described in claim 1, characterized in that, For each working condition, a corresponding skewed Gaussian mixture model is established, including: A1: Calculate the voltage source amplitude corresponding to each sampling point under the current operating condition. Impedance modulus and impedance phase , which serves as the variable parameter corresponding to the respective sampling point; A2: Based on the preset number of skewed Gaussian components K Clustering is performed on the variable parameters corresponding to each sampling point to obtain K Clustering with variable parameters; K It is a positive integer; A3: The first k Weights of each skewed Gaussian component Initialized to the first k The proportion of sampling points corresponding to the i-th variable parameter cluster, and the i-th variable parameter clustering corresponding to the i-th sampling point proportion .... k The mean of the skewed Gaussian components Standard deviation and skewness parameters Initialize to the first k The mean, standard deviation, and skewness parameters of each variable parameter cluster; among them... k =1,2,…… K ; A4: Calculate the variable parameters Each sample x i Belongs to the k The posterior probability of a skewed Gaussian component and will the k The weights of the skewed Gaussian components are updated as follows: ; N x This represents the total number of samples for the variable parameter; A5: For each skewed Gaussian component, update the mean with the objective of minimizing its negative log-likelihood function. Standard deviation and skewness parameters The possible values ​​of ; A6: Iterate through A4-A5 until the preset convergence condition is met; A7: The skewed Gaussian components are linearly combined according to their respective weights to obtain the skewed Gaussian mixture model corresponding to the current working condition.

3. The harmonic power flow probability assessment method based on a skewed Gaussian mixture model as described in claim 2, characterized in that, skewed Gaussian component number K The methods for determining include: B1: Initialization ; B2: Execute steps A1-A7 to obtain the current... The value corresponds to the skewed Gaussian mixture model; B3: According to Calculate the current The candidate scores corresponding to the values ​​of . ; B4: If Then If the value is increased by 1, proceed to B2; otherwise, proceed to B5. B5: Select the candidate with the lowest score. Determined as the number of skewed Gaussian components K ; in, K max This indicates the preset maximum number of skewed Gaussian components; For the present The value of corresponds to the maximum likelihood value of the skewed Gaussian mixture model.

4. The harmonic power flow probability assessment method based on the skewed Gaussian mixture model as described in claim 2, characterized in that, Voltage source amplitude corresponding to each sampling point under the current operating condition The calculation expression is: impedance magnitude and impedance phase The calculation methods include: according to Calculate the impedance in a three-branch parallel harmonic source model ; impedance Decompose it into magnitude and phase to obtain the impedance magnitude. and impedance phase ; in, h Indicates the current harmonic order; H Indicates the highest harmonic order. n h This indicates the total frequency of the voltage being considered; V N Indicates the rated voltage of the point of common coupling. V 1 represents the fundamental voltage amplitude. V q express q Second harmonic voltage amplitude; and Representing common connection points h Second harmonic voltage phasor and current phasor; This represents the phasor of the voltage source parameters, and its magnitude is... The phase is 0°; Indicates the parameters of the current source branch. This represents the admittance branch parameters.

5. The harmonic power flow probability assessment method based on the skewed Gaussian mixture model as described in claim 4, characterized in that, Current source branch parameters and admittance branch parameters The calculation methods include: Establish the following equation: The parameter matrix is ​​fitted using the complex partial least squares method. ; Extracting the parameter matrix The first column of elements yields the parameters of the current source branch. Extracting the parameter matrix The admittance branch parameters are obtained from the main diagonal elements. ; in, Indicates the first l sampling points h Secondary harmonic current Indicates the first l sampling points h Subharmonic voltage; , l =1,2…… L , L This indicates the number of sampling points selected.

6. The harmonic power flow probability assessment method based on a skewed Gaussian mixture model as described in any one of claims 1-5, characterized in that, The calculation methods for harmonic data at each node include: according to Harmonic source access node in computing system i The h Subharmonic current ; Calculate the first value of each node in the system according to the following equation. h Subharmonic voltage: The system contains a total of n There are nodes, of which the first m Each node is a node where a harmonic source is connected; and Representing nodes respectively i The parameters of the current source branch and the admittance branch in the three-branch parallel harmonic source model corresponding to the harmonic source are connected. and Representing nodes respectively i The impedance and voltage sources in the three-branch parallel harmonic source model corresponding to the harmonic source are connected; Represents a node i Accessed harmonic sources h Secondary equivalent admittance Represents a node i Accessed harmonic sources h Secondary equivalent current source Represents a node i and nodes j Between h Secondary network admittance i =1,2…… n , j =1,2…… n ; and Representing nodes respectively i of h Secondary harmonic voltage and current.

7. The harmonic power flow probability assessment method based on a skewed Gaussian mixture model as described in any one of claims 1-5, characterized in that, Clustering was performed on the sampling points in each daily period cluster to obtain multiple minute-level clusters, including: E1: Calculate the mean of the daily power time series of each daily period cluster as a typical daily power curve; E2: For each typical daily power curve, detect the points of change and use the detected points of change to divide the typical daily power curve into multiple intervals; E3: Extract the time-domain statistical features, waveform morphology features, trend features, and fluctuation features of each interval based on the sampling points contained within the interval, and combine them into the interval features of the corresponding interval. E4: Cluster each interval based on its characteristics to obtain multiple minute-level clusters.

8. The harmonic power flow probability assessment method based on the skewed Gaussian mixture model as described in claim 7, characterized in that, For each typical daily power curve, before detecting the points of change, the following steps are performed to determine the optimal number of points of change. and the corresponding change point positions : F1: Typical daily power curve Mapping to the reproducing kernel Hilbert space via Gaussian radial basis functions , to obtain the feature sequence ; T This indicates the number of sampling points included in the daily power time series. The first in a typical daily power curve t One sampling point, Indicates the mapping Corresponding feature points; F2: Set the number of change points ; F3: Initialization B The position of each change point and set , ; F4: For each interval obtained by dividing the points of change, extract the time-domain statistical features, waveform morphology features, trend features and fluctuation features of each interval based on the sampling points contained in the interval, and combine them into the interval features of the corresponding interval. F5: Construct the following cost function: in, Representing an interval The mean of the characteristic sequence, b =1,2,…… B ; Represents the regenerated nucleus Hilbert space norm, ; F6: If B If the number of change points is less than the preset maximum number, then... B If the value is increased by 1, proceed to F3; otherwise, proceed to F7. F7: Minimizes the cost function B Value as the optimal number of change points and the B The location of the change point corresponding to the value is determined as the optimal number of change points. The corresponding change point location.

9. The harmonic power flow probability assessment method based on the skewed Gaussian mixture model as described in claim 7, characterized in that, The improved elbow rule is used to determine the optimal number of clusters when clustering the multiple daily power time series obtained from the division, as well as when clustering each interval in step E3. The improved elbow rule includes: Within the candidate clustering range Traverse within the range, for each r After clustering the samples to be clustered using this as the clustering quantity, according to... Calculate the inertia value and construct an inertia value curve that shows the change of inertia value with the number of clusters. Determine the number of clusters corresponding to the point with the largest rate of change of curvature as the optimal number of clusters. Cluster the samples to be clustered according to the determined optimal number of clusters; When clustering multiple daily power time series obtained from the division, the samples to be clustered are daily power time series; when clustering each interval, the samples to be clustered are interval features. For the first c A sample set of clusters, Indicates the first c The centroid sequence of each cluster; Representing the i One sample; DTW() represents the dynamic time warp distance; r max This indicates the preset maximum number of clusters.

10. A harmonic power flow probability assessment system based on a skewed Gaussian mixture model, characterized in that, include: A computer-readable storage medium for storing computer programs; And a processor for reading a computer program stored in the computer-readable storage medium and executing the harmonic power flow probability assessment method based on the skewed Gaussian mixture model as described in any one of claims 1-9.

Images