A power quality disturbance signal classification and identification method based on VMD and HMM

Abstract

The application provides a power quality disturbance signal classification and recognition method based on VMD and HMM, comprising the following steps: generating signals of different disturbance types through simulation, dividing the signals into training data and test data according to a certain proportion, and decomposing each signal through VMD, with the number of decomposition layers being 4; calculating the energy value E of each layer of IMF to form a feature vector, which is used as an observation sequence after discretization; training HMMs according to different disturbance types to obtain an HMM model library; inputting the observation sequence obtained after processing the test data into the trained HMM library, calculating the log-likelihood probability value corresponding to each model, and outputting the disturbance type corresponding to the HMM with the maximum log-likelihood, which is the type of the current disturbance signal. The application can accurately detect and locate power quality disturbances, extract disturbance features and identify the specific type of the disturbance.

Description

Technical Field

[0001] This invention relates to the field of power quality disturbance detection and analysis technology, specifically a power quality disturbance signal classification and identification method based on VMD and HMM. Background Technology

[0002] With the accelerating pace of smart grid construction and the increasing maturity of distributed generation technology, a large number of complex and diverse nonlinear loads, power electronic devices, and distributed power sources are being connected to the grid, leading to increasingly complex power system structures and a deteriorating power quality. Against this backdrop, power quality has become a key concern for power companies and users, a necessary prerequisite for improving and addressing power quality issues, and the foundation for the safe operation of the power system, thus possessing significant research value. Summary of the Invention

[0003] The purpose of this invention is to provide a method for identifying power quality disturbance signals based on VMD and HMM, which can accurately detect and locate power quality disturbances, extract disturbance features, and identify the specific type of disturbance.

[0004] A method for identifying power quality disturbance signals based on VMD and HMM includes the following steps:

[0005] Step 1: Generate signals with different perturbation types through simulation, divide the signals into training data and test data according to a certain ratio, and perform VMD decomposition on each type of signal, with a decomposition layer of 4 layers;

[0006] Step 2: Calculate the energy value E of each IMF layer to form an eigenvector, and discretize it to obtain the observation sequence;

[0007] Step 3: Train HMMs individually according to different perturbation types to obtain an HMM model library;

[0008] Step 4: Input the observed sequence obtained after processing the test data into the trained HMM library, solve the log-likelihood probability value of the sequence under each model, and output the perturbation type corresponding to the HMM with the maximum log-likelihood, which is the type of the current perturbation signal.

[0009] Furthermore, the waveforms of the different disturbance types include four single disturbance signals: harmonics (S1), voltage sag (S2), voltage swell (S3), and voltage interruption (S4), as well as two composite disturbance signals: sag + harmonics (S5) and swell + harmonics (S6), for a total of six disturbance signals, which are generated by simulation using MATLAB software.

[0010] Furthermore, the VMD decomposition process in step one is specifically as follows:

[0011] The expression for the modal components obtained from VMD decomposition is:

[0012] u k (t)=A k cos(φ k (t)) (1)

[0013] In the formula, A k Represented as the instantaneous amplitude of the IMF, also known as the envelope function; φ k (t) represents the phase of the modal component; u k (t) represents the modulation and demodulation signal;

[0014] The analytic signals corresponding to each IMF are obtained using the Hilbert transform, and their one-sided spectra are calculated as follows:

[0015]

[0016] In the formula, δ(t) represents the impact function, which is expressed by an exponential function. The center frequency of each IMF is adjusted by multiplying with the single-sided spectrum, and its spectrum is modulated onto the corresponding fundamental frequency band:

[0017]

[0018] Based on the Gaussian smoothing estimation method, the bandwidth of each component of the demodulated signal is estimated, which is the square norm, resulting in the constrained variational model of VMD:

[0019]

[0020] Where: K represents the number of modes, {u k}={u1,u2,...,u k} represents the set of IMFs obtained from the decomposition of signal f(t), {w k}={w1,w2,...,w k} represents the set of center frequencies corresponding to each IMF; by introducing a quadratic penalty factor α and a Lagrange operator λ, the constrained variational problem is transformed into an unconstrained variational problem to obtain the optimal solution of the model. The specific equation is as follows:

[0021]

[0022] Where <> represents the inner product operation, the variational problem above is solved using the alternating direction method of multipliers to obtain the minimum point of the model, i.e., by iteratively updating u. k n+1 w k n+1 , λ n+1 By searching for the "saddle point" in the above expression, we obtain... The expression:

[0023]

[0024] Where n represents the number of iterations, and by transforming equation (6) to the frequency domain according to Parseval / Plancherel, we can further obtain:

[0025]

[0026] Similarly, by performing a similar optimization on the center frequency, we obtain:

[0027]

[0028] In the formula, represent Wiener filtering, It is the current center frequency of the IMF. Perform the inverse Fourier transform to obtain the real part u. k (t), thus obtaining the time-domain form of each IMF.

[0029] Furthermore, step two specifically includes:

[0030] First, extract the energy characteristics of the intrinsic mode components of each layer of the signal:

[0031]

[0032] Where E(j) is the energy value of the j-th intrinsic mode component, j = 1, 2, ..., K, and K is the VMD decomposition layer number; x(t) is the signal sequence; N is the number of sampling points of the sample. The energy values ​​of each component obtained from the decomposition of different training samples are used to form a feature vector F = [E1, E2, E3, E4]. Then, the K-means clustering method is used to quantify the feature values. The specific process is as follows:

[0033] Step 1: Initialization: Randomly select K samples as the initial cluster centroids {μ1,μ2,…,μ...} K};

[0034] Step 2, Iterative Update: Calculate the distance from each point i in the space to each centroid, and assign the current point to the nearest centroid, thereby generating K clusters. The specific calculation formula is as follows:

[0035]

[0036] For each point in the cluster, recalculate its distance to the centroid, and use the point corresponding to the average distance of all points as the new centroid of the cluster:

[0037]

[0038] Where K is the number of clusters, C (i) μ represents the class closest to the sample point. j Indicates the new center of mass;

[0039] Step 3, Termination: Stop iterating when the centroid of each cluster no longer changes or the maximum number of iterations is reached, and the clustering ends.

[0040] After the above steps, the discretized features are obtained, which are then used as the observation sequence input into the HMM for training.

[0041] Furthermore, in step three, the Hidden Markov Model (HMM) is trained using the Baum-Welch algorithm. The specific process is as follows:

[0042] Define variable ζ t (i,j) is used to characterize the state S of the model at time t, given the observation sequence O and model parameters λ. i And the state corresponding to time t+1 is S j The probability is:

[0043] ζ t (i,j)=P(q t =i,q t+1 =j|O,λ) (12) Then, based on the α derived from the forward-backward algorithm t (i) and β t (i) Perform the calculation:

[0044]

[0045] Given λ and O, the state of the model at time t is S. i The probability is:

[0046] γ t (i)=P(q t =S i |O,λ) (14)

[0047] γ is calculated using a forward-backward algorithm. t (i), that is:

[0048]

[0049] Finally, the parameters of the HMM model are re-estimated using the following formula:

[0050]

[0051] Repeat the above calculations and update the model parameters until the convergence of P(O|λ) reaches the required level. The final λ obtained is the desired HMM model.

[0052] Furthermore, the convergence of P(O|λ) is defined as the change in model output being less than ε or the number of iterations reaching 30, where ε is set to 1×10. -4 .

[0053] Furthermore, in step four, the 40 sets of test signals with different disturbance signals are subjected to 4-level VMD decomposition and the energy value of the intrinsic mode components of each level are calculated. The energy values ​​are combined to obtain a 1*4 feature vector, which is also discretized and input as the observation sequence into 6 HMM models representing different disturbance types. Each HMM will output the log-likelihood probability value of the test data, which represents the degree of matching between the observation sequence and the current HMM. The disturbance type corresponding to the HMM with the largest output log-likelihood is the type of the current disturbance signal.

[0054] This invention decomposes each type of disturbance signal into VMD (Vibration Mode Decomposition) with four layers. The energy value E of each layer's IMF is then calculated to form a feature vector, which is discretized and used as an observation sequence. Subsequently, Hidden Markov Models (HMMs) are trained separately for each disturbance type. The observation sequences obtained after processing the test data are input into the trained HMM library for current disturbance signal type identification. The combined use of VMD and HMM can accurately detect and locate power quality disturbances, extract disturbance features, and identify the specific type of disturbance, thereby improving and mitigating power quality problems. Attached Figure Description

[0055] Figure 1 This is a flowchart of the VMD algorithm in this invention;

[0056] Figure 2 This is a flowchart of one embodiment of the power quality disturbance signal identification method based on VMD and HMM of the present invention;

[0057] Figure 3 The training of this invention yields HMM models corresponding to six perturbation types upon completion;

[0058] Figure 4 The log-likelihood values ​​of six test data under the harmonic HMM model of this invention. Detailed Implementation

[0059] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] like Figure 2 As shown, this embodiment of the invention provides a method for identifying power quality disturbance signals based on VMD and HMM, including the following steps;

[0061] A: By simulating and generating signals of different disturbance types, the signals are divided into training data and test data according to a certain ratio, and each signal is decomposed by VMD with 4 decomposition layers. The different disturbance types of signals mainly include four single disturbance signals: harmonics (S1), voltage sag (S2), voltage swell (S3), and voltage interruption (S4), and two composite disturbance signals: sag + harmonics (S5) and swell + harmonics (S6), for a total of 6 disturbance signal waveforms, which can be generated by simulation using MATLAB software.

[0062] The fundamental frequency of the signal is 50Hz, and the third harmonic is 150Hz. The sampling frequency is set to 10kHz, and the simulation duration is 0.5s. Considering the possibility of noise in actual situations, the signal-to-noise ratio of the simulated signals is 20dB. 100 sets of data are randomly generated for each perturbation signal, of which 60 sets are used as training data to train the HMM model, and the remaining 40 sets are used as test data to test the model's recognition accuracy.

[0063] B: Calculate the energy value E of each IMF layer to form an eigenvector, and discretize it to serve as the observation sequence;

[0064] The training samples are decomposed into 4-level VMD to obtain each intrinsic mode function and calculate its energy value. The energy values ​​of each component obtained from the decomposition of different training samples are used to form a feature vector F = [E1, E2, E3, E4]. After discretization, it is used as the observation sequence. The 6 perturbation types represent 6 hidden states, which are denoted as λ1, λ2, λ3, λ4, λ5, λ6 respectively.

[0065] C: Train HMMs individually according to different perturbation types to obtain an HMM model library;

[0066] Based on the initialization of the three main parameters of the Hidden Markov Model (HMM), the Baum-Welch algorithm is used for model training. During the training process, the log-likelihood probability value is output for each iteration. The convergence error ε is set to 1×10⁻⁶. -4 The maximum number of iterations is set to 30. The iteration stops when the change in the model output is less than ε or the number of iterations reaches 30, thus obtaining the optimal HMM parameters and obtaining the HMM models corresponding to the six types of perturbations.

[0067] D: Input the observed sequence obtained after processing the test data into the trained HMM library, solve the log-likelihood probability value of the sequence under each model, and output the perturbation type corresponding to the HMM with the maximum log-likelihood, which is the type of the current perturbation signal.

[0068] The VMD decomposition process in step A is as follows: Figure 1 As shown, the modal component expressions obtained from VMD decomposition are:

[0069] u k (t)=A k cos(φ k (t)) (1)

[0070] In the formula, A k Represented as the instantaneous amplitude of the IMF, also known as the envelope function; φ k (t) represents the phase of the modal component; u k (t) represents the modulation and demodulation signal.

[0071] The analytic signals corresponding to each IMF are obtained using the Hilbert transform, and their one-sided spectra are calculated as follows:

[0072]

[0073] In the formula, δ(t) represents the impact function. This is achieved through an exponential function. The center frequency of each IMF is adjusted by multiplying with the single-sided spectrum, and its spectrum is modulated onto the corresponding fundamental frequency band:

[0074]

[0075] Based on the Gaussian smoothing estimation method, the bandwidth of each component of the demodulated signal is estimated, which is the square norm. The constrained variational model of VMD can then be obtained as follows:

[0076]

[0077] Where: K represents the number of modes, {u k}={u1,u2,...,u k} represents the set of IMFs obtained from the decomposition of signal f(t), {w k}={w1,w2,...,w k Let} represent the set of center frequencies corresponding to each IMF. Solving the above equation directly is quite difficult, so a quadratic penalty factor α and a Lagrange operator λ are introduced. The aim is to transform the constrained variational problem into an unconstrained variational problem to obtain the optimal solution of the model. The specific equation is as follows:

[0078]

[0079] Here, <> represents the inner product operation. The purpose of introducing α is to ensure good reconstruction accuracy even with noise components; its weights are inversely proportional to the noise content in the signal. λ ensures the constraints are strictly enforced. The variational problem is solved using the alternating direction method of multipliers to obtain the model's minimum point, i.e., by iteratively updating u. k n+1 w k n+1 , λ n+1 Searching for the "saddle point" in the above expression yields... The expression:

[0080]

[0081] Where n represents the number of iterations. Transforming equation (6) to the frequency domain using Parseval / Plancherel, we can further obtain:

[0082]

[0083] Similarly, by performing a similar optimization on the center frequency, we can obtain:

[0084]

[0085] In the formula, represent Wiener filtering, It is the current center frequency of the IMF. Perform the inverse Fourier transform to obtain the real part u. k (t), thus the time-domain form of each IMF can be obtained.

[0086] Step B first extracts the energy characteristics of the intrinsic mode components of each layer of the signal:

[0087]

[0088] Where E(j) is the energy value of the j-th intrinsic mode component, j = 1, 2, ..., K, and K is the VMD decomposition layer number; x(t) is the signal sequence; N is the number of sampling points of the sample. The energy values ​​of each component obtained from the decomposition of different training samples are used to form a feature vector F = [E1, E2, E3, E4]. The K-means clustering method is then used to quantify the feature values. The specific process is as follows:

[0089] Step 1: Initialization: Randomly select K samples as the initial cluster centroids {μ1,μ2,...,μ...} K}

[0090] Step 2, Iterative Update: Calculate the distance from each point i in the space to each centroid, and assign the current point to the nearest centroid, thereby generating K clusters. The specific calculation formula is as follows:

[0091]

[0092] For each point in the cluster, recalculate its distance to the centroid, and use the point corresponding to the average distance of all points as the new centroid of the cluster:

[0093]

[0094] Where K is the number of clusters, C (i) μ represents the class closest to the sample point. j Indicates the new center of mass;

[0095] Step 3, Termination: Stop iterating when the centroid of each cluster no longer changes or the maximum number of iterations is reached, and the clustering ends.

[0096] After the above steps, the discretized features are obtained, which can then be used as the input of the observation sequence into the HMM for training.

[0097] In step C, the Hidden Markov Model (HMM) is trained using the Baum-Welch algorithm. The specific process is as follows:

[0098] Define variable ζ t (i,j) is used to characterize the state S of the model at time t, given the observation sequence O and model parameters λ. i And the state corresponding to time t+1 is S j The probability is:

[0099] ζ t (i,j)=P(q t =i,q t+1 =j|O,λ) (12)

[0100] Then, based on the α derived from the forward-backward algorithm t (i) and β t (i) can be calculated:

[0101]

[0102] Given λ and O, the state of the model at time t is S. i The probability is:

[0103] γ t (i)=P(q t =S i |O,λ) (14)

[0104] γ can be calculated using the forward-backward algorithm. t (i), that is:

[0105]

[0106] Finally, the parameters of the HMM model are re-estimated using the following formula:

[0107]

[0108] Repeat the above calculations and update the model parameters until the convergence of P(O|λ) reaches the required level. The final λ obtained is the desired HMM model.

[0109] In step D, 4-level VMD decomposition and energy value calculation of the intrinsic mode components at each level are performed on 40 sets of test signals with different perturbation signals. The energy values ​​are combined to obtain a 1*4 feature vector, which is also discretized and input as the observation sequence into 6 HMM models representing different perturbation types. Each HMM outputs the log-likelihood probability value of the test data, representing the degree of matching between the observation sequence and the current HMM. The perturbation type corresponding to the HMM with the highest output log-likelihood is the type of the current perturbation signal. Figure 3 The output results for 40 sets of test data with various perturbations under different Hidden Markov Models (HMMs) are shown. The horizontal axis represents the test data for six perturbation types, and the vertical axis represents the log-likelihood probability of the test data under that HMM. Figure 4 The values ​​represent the log-likelihood of six test data sets under the harmonic HMM model.

[0110] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for identifying power quality disturbance signals based on VMD and HMM, The features include the following steps: Step one, generate signals of different disturbance types by simulation, divide the signals into training data and test data according to a certain proportion, and decompose each signal by VMD with 4 layers; Step two, calculate the energy values E of each layer IMF to form a feature vector, which is discretized as an observation sequence; Step three, train HMM separately according to different disturbance types to obtain an HMM model library; Step four, input the observation sequence obtained by processing the test data into the trained HMM library, solve the log-likelihood probability value of the sequence under each model, and output the disturbance type corresponding to the HMM with the maximum log-likelihood, which is the type of the current disturbance signal; Step two specifically includes: First, extract the energy features of each layer of the intrinsic modal component of the signal: (9)； wherein, is the energy value of the jth intrinsic modal component, K is the number of VMD decomposition layers; is the signal sequence; N is the number of sampling points of the sample, and the component energy values obtained by decomposing different training samples form a feature vector The K-means clustering method is then used to quantize the feature values, and the specific process is as follows: Step 1, initialization: randomly select K samples as initial cluster centroids ; Step 2, Iterative update: Calculate the distance from each point in the space to each centroid and assign the current point to the nearest centroid, thereby generating K clusters. The formula is as follows: to each centroid and assign the current point to the nearest centroid, thereby generating K clusters. The formula is as follows: (10)； For each point in the cluster, recalculate the distance from the centroid, and take the point corresponding to the average value of all points as the new centroid of the cluster: (11)； where K is the number of clusters, denotes the class closest to the sample point, denotes the new centroid; Step 3, termination: when the centroid in each cluster no longer changes or reaches the maximum number of iterations, stop iteration and clustering; After the above steps, the discretized feature quantity is obtained, which is input into the HMM for training; in step three, the HMM is trained by the Baum-welch algorithm, and the specific process is as follows: Define variables to characterize the conditional probability of a given observation sequence and model parameters under the condition that the model is in state at time and in state at time , i.e. (12)； Further according to the forward-backward algorithm derived and The calculation is performed: (13)； When and Given the situation at time t, the model's probability of being in state s at time t+1 is: P(s | s, t)​ (14)； The forward-backward algorithm is used to calculate i.e. (15)； Finally, re-estimate the parameters of the HMM model according to the following formula: (16)； The above calculation is repeated, the model parameters are updated, and the calculation is stopped until the convergence degree reaches the requirement, and the finally obtained is the HMM model to be solved. In step four, 40 groups of test signals of different disturbance signals are decomposed by 4-layer VMD and the energy values of each layer of intrinsic modal component are calculated, the energy values are combined to obtain a 1*4 feature vector, which is also discretized and input into the 6 HMM models representing different disturbance types, each HMM outputs the log-likelihood probability value of the test data, which represents the matching degree of the observation sequence and the current HMM, and the disturbance type corresponding to the HMM with the maximum log-likelihood is the type of the current disturbance signal.

2. The method of claim 1, wherein the method is based on VMD and HMM. The different disturbance types of signals include four single disturbance signals: harmonic (S1), voltage sag (S2), voltage swell (S3), and voltage interruption (S4), and two composite disturbance signals: sag + harmonic (S5) and swell + harmonic (S6), a total of 6 disturbance signal waveforms, which are generated by MATLAB software simulation.

3. The method of claim 2, wherein the VMD and HMM-based power quality disturbance signal recognition method is characterized by: The VMD decomposition process in step one is as follows: The modal component expression obtained by VMD decomposition is: (1)； wherein denotes the instantaneous amplitude of the IMF, also referred to as the envelope function; denotes the phase of the modal component; is the modulated signal; The analytic signal corresponding to each IMF is solved by Hilbert transform, and its unilateral spectrum is calculated as follows: (2)； wherein represents the impact function, which is adjusted by an exponential function The center frequency of each IMF is adjusted by multiplying it with a monospectral in such a way that its spectrum is modulated to the corresponding base band: (3)； Based on the Gaussian smoothing estimation method, the bandwidth of each component of the demodulated signal is estimated, that is, the square norm, and the constraint variational model of VMD is obtained as: (4)； wherein: K represents the number of modes, for the signal a set of IMFs obtained by decomposition, denotes a set of center frequencies corresponding to each IMF; a quadratic penalty factor and a Lagrange operator , the above constrained variational problem is converted into an unconstrained variational problem to obtain the optimal solution of the model, and the specific equation is: (5)； wherein represents an inner product operation, the above variational problem is solved according to a multiplier alternating direction method to obtain a minimum point of the model, i.e. by iteratively updating searching for a "saddle point" of the above expression to obtain the expression of (6)； Where n represents the number of iterations, according to Parseval / Plancherel, equation (6) is transformed into the frequency domain, and further obtained: (7)； Similarly, the center frequency is solved by similar optimization: (8)； In the formula, represent Wiener filtering, It is the current center frequency of the IMF. Perform the inverse Fourier transform to obtain the real part. Thus, the time-domain form of each IMF is obtained.

4. The method of claim 1, wherein the method is characterized by: The Convergence is that the model output changes less than Or the number of iterations reaches 30 times, Set to .

Images