Voltage quality composite disturbance identification method considering the fusion of time-frequency features of complex signals

Abstract

This invention relates to a method for identifying complex voltage quality disturbances by considering the fusion of time-frequency features of complex signals, comprising the following steps: Step S1, acquiring the time-frequency domain distribution information of the voltage quality disturbance signal; Step S2, dividing the disturbance identification input feature index into time-domain statistical features, frequency-domain statistical features, and energy features; Step S3, extracting deep features from the voltage quality disturbance information and then fusing them with the disturbance identification input feature index; Step S4, assigning weights to the fused feature set; Step S5, using an ensemble learner XGBoost based on a Boosting strategy as the disturbance classifier, inputting the weighted fused feature set into XGBoost to achieve the identification of complex disturbance signals. This invention can solve the technical problems of tracing the source of power quality disturbances and improving power system stability.

Description

Technical Field

[0001] This invention belongs to the field of voltage quality disturbance identification technology, and relates to a method for identifying composite voltage quality disturbances, particularly a method for identifying composite voltage quality disturbances that considers the fusion of time-frequency features of complex signals. Background Technology

[0002] Against the backdrop of the energy revolution driven by the "dual carbon" goal, the massive influx of heterogeneous distributed energy sources and nonlinear loads into the power grid, while enhancing the grid's flexibility, has also triggered more complex and variable voltage quality disturbances, seriously affecting the stability and economy of power system operation.

[0003] Voltage quality disturbance identification is a key component of voltage quality monitoring systems in smart grid paradigms and a crucial prerequisite for managing and improving voltage quality. However, in the complex process of building new power systems, voltage quality disturbances have evolved from traditional single disturbances to hybrid disturbances, which are composed of various basic voltage quality disturbances of different types, intensities, and start and end times.

[0004] Moreover, the composite voltage quality disturbance characteristics are not simply a linear superposition of the basic voltage quality disturbance characteristics, but may be overlapping and complexly intersecting between time-frequency domain characteristics, posing a huge challenge to traditional voltage quality disturbance identification techniques.

[0005] Therefore, extracting appropriate voltage quality indicators and using artificial intelligence algorithms to accurately identify voltage quality disturbances has become a major research challenge in the field of voltage quality in recent years.

[0006] A search revealed no publicly available literature of the same or similar prior art as this invention. Summary of the Invention

[0007] The purpose of this invention is to overcome the shortcomings of the prior art and propose a voltage quality composite disturbance identification method that considers the fusion of time-frequency characteristics of complex signals, which can solve the technical problems of power quality disturbance tracing and power system stability.

[0008] The present invention solves its practical problem by adopting the following technical solution:

[0009] A method for identifying complex voltage quality disturbances that considers the fusion of time-frequency features of complex signals includes the following steps:

[0010] Step S1: Obtain the time-frequency domain distribution information of the voltage quality disturbance signal using the improved Hilbert-Huang transform (IHHT) method;

[0011] Step S2: Based on the time-frequency domain distribution information of the voltage quality disturbance signal obtained in step S1, design disturbance identification input feature indexes and divide these disturbance identification input feature indexes into three categories: time-domain statistical features, frequency-domain statistical features, and energy features. Construct a more comprehensive feature set from the perspectives of amplitude, frequency, phase, and energy.

[0012] Step S3: Employ a time-frequency domain joint discrimination enhanced sparse autoencoder to extract deep features from voltage quality disturbance information through time-frequency domain joint training. Then, fuse these features with the disturbance identification input feature index manually designed in Step S2 to obtain a more comprehensive fusion feature set.

[0013] Step S4: The dynamic feature weighting method based on dual-delay deep deterministic policy gradient (TD3) and the state space, action space and reward function are used to guide TD3 to learn the optimal dynamic input feature adjustment method. In the offline application stage, the weight of the input features can be adaptively adjusted according to the voltage signal change trend to weight the fused feature set in step S3.

[0014] Step S5: The XGBoost ensemble learner based on the Boosting strategy is used as the perturbation classifier. The fusion feature set obtained in step S3 is weighted in step S4 and then input into XGBoost to identify the composite perturbation signal.

[0015] Advantages and beneficial effects of the present invention:

[0016] 1. This invention proposes a voltage quality composite disturbance identification method that considers the fusion of time-frequency features of complex signals. It involves time-frequency domain feature extraction, dynamic feature weight adjustment, and disturbance identification classification model construction, which effectively improves the accuracy of voltage quality disturbance identification.

[0017] 2. This invention proposes an improved Hilbert-Huang transform method to analyze the time-frequency domain characteristics of perturbation signals, which alleviates the mode mixing and endpoint effects problems caused by the traditional EMD method.

[0018] 3. This invention proposes a discriminative enhanced sparse autoencoder, which extracts time-frequency domain features from power quality disturbance information through joint time-frequency domain training. It achieves an effective combination of manual feature extraction and deep learning adaptive feature extraction, which can both design indicators based on expert experience and adaptively learn the complex patterns and features of disturbance signals.

[0019] 4. This invention proposes a dynamic feature weighting method based on dual-delay deep deterministic policy gradient (TD3), and uses an integrated learner XGBoost to classify perturbation signals, thereby improving the efficiency of perturbation identification while reducing redundant features and noise feature interference. Attached Figure Description

[0020] Figure 1 This is a flowchart of the voltage quality composite disturbance identification method that considers the fusion of time-frequency features of complex signals provided by the present invention;

[0021] Figure 2 This is a partial voltage quality disturbance waveform generated by the mathematical model in this invention;

[0022] Figure 3 This is a diagram illustrating the effect of signal decomposition using ICEEMDAN in this invention.

[0023] Figure 4 This is a two-dimensional spatial distribution map of different perturbation signals after feature extraction by IHHT in this invention;

[0024] Figure 5 This is a marginal spectrum curve diagram of different perturbation signals in this invention;

[0025] Figure 6 This is a comparison chart of the accuracy of perturbation identification for each model under different signal-to-noise ratios in this invention. Detailed Implementation

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

[0027] A method for identifying complex voltage quality perturbations that considers the fusion of time-frequency features of complex signals, such as... Figure 1 As shown, it includes the following steps:

[0028] Step S1: Use the improved Hilbert-Huang transform (IHHT) method to obtain the time-frequency domain distribution information of the voltage quality disturbance signal;

[0029] An improved adaptive noise fully integrated empirical mode decomposition (ICEEMDAN) is used to replace the empirical mode decomposition (EMD) in the traditional Hilbert-Huang transform. The best IMF components with clear features are selected from a series of intrinsic mode functions (IMFs) obtained from ICEEMDAN using the frequency domain cross-correlation coefficient for time-frequency domain feature extraction.

[0030] The specific steps of step S1 include:

[0031] The IHHT in step S1 includes three processes: ICEEMDAN mode decomposition, frequency domain cross-correlation coefficient feature selection, and Hilbert transform.

[0032] Step S101: Inject the modal components of the i-th group of white noise into the original time series S using ICEEMDAN to construct a new time series X. (i) :

[0033] X (i)=S+β0E1(ω (i) )

[0034] β0=ε0σ(S) / σ(E1(ω (i) ))

[0035] In the formula: β0 is the noise figure, σ(·) is the standard deviation operator, and E i (·) represents the i-th mode obtained from EMD decomposition, ω (i) Let ε0 be the white noise of the i-th group (i = 1, 2, ..., I) which follows an N(0, 1) distribution, where I is the total number of modes decomposed, and ε0 is the amplitude of the white noise.

[0036] The working principle of step S101 is as follows: Traditional HHT consists of two parts: Empirical Mode Decomposition (EMD) and Hilbert Transform. The ICEEMDAN introduced in IHHT overcomes the mode mixing and endpoint effects problems existing in traditional EMD.

[0037] Step S102: Calculate the local mean to obtain the residual R1 of the first decomposition.

[0038] R1 = <M(X (i) ).

[0039] In the formula: <·. represents the average value calculation, and M(·) is the local mean value calculation function for the electrical energy signal;

[0040] Step S103: When the modal decomposition order of the sampled signal is 1, calculate the first modal component IMF1 using the original signal S and the residual signal R1, that is:

[0041] IMF1 = S-R1

[0042] Step S104: Add white noise to the residual to construct a new signal to be decomposed, R1+β1E2(ω). (i) Furthermore, the second set of residuals R2 and its modal components IMF2 can be calculated:

[0043] IMF2 = R1 - R2 = R1 - <M(R1+β1E2(ω (i) ))>

[0044] Step S105: Calculate the residual R of the k-th decomposition, and so on. k and its modal components IMF k :

[0045] R k = <M(R k-1 +β k-1 E k ω (i) ).

[0046] IMF k =R k-1 -R k

[0047] In the formula: when k≥1, β k =ε0σ(R) k )

[0048] Step S106: Repeat step S105 to obtain all modal components and the final residual;

[0049] Step S107: Using the frequency domain correlation coefficient as the evaluation criterion, select the best IMF component with clear features from a series of IMF components obtained in step S106 for signal reconstruction and time-frequency domain feature extraction.

[0050] The working principle of the frequency domain correlation coefficient in step S107 is as follows: assuming G x and G y These are signals x i and y i spectral power, f a If the analysis frequency is the signal x, then the signal x i and y i The cross-correlation coefficient in the frequency domain can be expressed as:

[0051]

[0052] In the formula: G xi and G yi G x and G y Spectral power at frequency i; and G x and G y The average spectral power.

[0053] Step S108: Based on the dominant IMF obtained in step S107, the local amplitude and frequency information existing in the processed signal is obtained through Hilbert transform;

[0054] The IMF after Hilbert transformation can be represented as:

[0055]

[0056] In the formula: H(·) is the Hilbert transform function; P is the Cauchy principal value; IMF k (t) represents the value of the k-th eigenmode function at time t; τ represents the integration variable;

[0057] Step S109: Construct the analytic signal Z based on the Hilbert transform result obtained in step S108. k (t)=IMFk (t)+jH[IMF k [(t)], the magnitude function a of each IMF k (t) and phase function θ k (t) can be calculated as:

[0058]

[0059] In the formula: Z k (t) represents the analytic signal corresponding to the k-th IMF; arctan is the arctangent function.

[0060] Step S110: Based on the phase function θ obtained in step S109 k (t), calculate the instantaneous frequency, Hilbert spectrum, and final marginal spectrum of IHHT to obtain the time-frequency domain distribution information of the voltage quality disturbance signal, in order to design subsequent feature indicators:

[0061]

[0062] In the formula: T is the length of the signal; ω k H(ω,t) represents the instantaneous frequency of the k-th IMF at time t; H(ω,t) represents the Hilbert spectrum, which accurately reflects the trend of amplitude variation with duration and frequency; h(ω) represents the marginal spectrum, which reflects the contribution of each frequency to the amplitude.

[0063] Step S2: Based on the time-frequency domain distribution information of the voltage quality disturbance signal obtained in step S1, design disturbance identification input feature indexes and divide these disturbance identification input feature indexes into three categories: time-domain statistical features, frequency-domain statistical features, and energy features. The aim is to construct a more comprehensive feature set from the perspectives of amplitude, frequency, phase, and energy.

[0064] The disturbance identification input feature index design in step S2 includes time-domain statistical features, frequency-domain statistical features, and energy features;

[0065] The specific steps of step S2 include:

[0066] Step S201: Design the time-domain statistical features as: maximum voltage amplitude, minimum voltage amplitude, average voltage amplitude, standard deviation of voltage, peak voltage factor, and duration of disturbance.

[0067] The working principle of step S201 is that the original perturbation signal contains important information about the PQD waveform, which enhances the classifier's ability to understand and represent the data.

[0068] Step S202: After the Hilbert transform, the instantaneous characteristics of each IMF can be obtained. Therefore, the following time-domain statistical characteristics are added: maximum instantaneous amplitude of each IMF, minimum instantaneous amplitude of each IMF, average instantaneous amplitude of each IMF, standard deviation of instantaneous amplitude of each IMF, and standard deviation of instantaneous phase of each IMF.

[0069] Step S203: Design the frequency domain statistical features as follows: the top three amplitudes in the marginal spectrum and their corresponding frequencies, the minimum marginal spectrum amplitude, the average marginal spectrum amplitude, the standard deviation of the marginal spectrum amplitude, the maximum high-frequency marginal spectrum amplitude, the minimum high-frequency marginal spectrum amplitude, the average high-frequency marginal spectrum amplitude, and the standard deviation of the high-frequency marginal spectrum amplitude.

[0070] The working principle of step S203 is that traditional studies mostly focus only on the instantaneous amplitude-phase characteristics obtained after Hilbert transform, ignoring the marginal spectral characteristics.

[0071] Step S204: Design the energy characteristics as: IMF energy and total marginal spectrum energy.

[0072] IMF Energy E imf This reflects the strength of the IMF, and the specific calculation method is as follows:

[0073]

[0074] Where: IMF k (i) can be the i-th sampling point in the k-th IMF; N is the total number of sampling points.

[0075] The marginal spectrum is obtained by integrating the Hilbert spectrum over the time axis, reflecting the total energy of the signal at each frequency. Therefore, the total marginal spectral energy is defined as:

[0076]

[0077] In the formula: N represents the number of sampling points of the marginal spectrum curve; h(ω,i) represents the i-th sampling point of the marginal spectrum curve h(ω).

[0078] The working principle of step S204 is that energy characteristics can quantify the energy distribution of a signal in the entire time and frequency domain, thereby reflecting the intensity or severity of the disturbance.

[0079] Step S3: Employ a time-frequency domain joint discrimination enhanced sparse autoencoder to extract deep features from voltage quality disturbance information through time-frequency domain joint training. Then, fuse these features with the disturbance identification input feature index manually designed in Step S2 to obtain a more comprehensive fusion feature set.

[0080] The specific steps of step S3 include:

[0081] Step S301: Based on the encoding function of the autoencoder (AE), the input signal is converted into a compressed feature representation through weighting and biasing:

[0082] H = f e (W e ·X+b e )

[0083] X R =f d (W d ·H+b d )

[0084]

[0085] In the formula: W e and b e These represent the weights and biases from the input layer to the hidden layer, respectively. d and b d X, H, and X represent the weights and biases from the hidden layer to the output layer, respectively. R These represent the original input, hidden layer output, and reconstructed data, respectively; f e (·) and f d (·) is the activation function; J AE X represents the autoencoder reconstruction error; i and Let i and j represent the i-th original input and the reconstructed data, respectively.

[0086] Step S302: Based on step S301, construct a sparse autoencoder (SAE). The sparsity constraint affects the encoder's behavior by adding a sparsity penalty term to the loss function. KL divergence is used to measure the difference between the actual activation distribution of the hidden layer and the target sparsity.

[0087]

[0088] In the formula: ρ represents the sparsity parameter of the KL divergence; represents the average activation of hidden layer node j; the purpose of the penalty term is to make active hidden nodes close to 1 and inactive nodes close to 0, forcing the model to learn more generalizable feature representations.

[0089] Furthermore, the autoencoder loss function after introducing sparsity constraints can be expressed as:

[0090]

[0091] In the formula: J SAE (W,b) is the loss function after introducing sparsity constraints; β is the weighting coefficient of KL divergence.

[0092] The working principle of step S302 is that adding sparsity constraints to the autoencoder can filter out redundant information in the input data, thereby better focusing on the key features of the data.

[0093] Step S303: Based on step S302, a discriminator module is added to construct a class label-assisted discriminative enhanced sparse autoencoder (DESAE), thereby improving the quality of feature extraction. During training, the hidden layer features of the SAE are input into the discriminator, and the class label information helps to enhance the attention to perturbation signals. Assuming a given training sample {X,Y}, where Y represents the perturbation class to which the sample belongs, the discriminator calculates the perturbation class information it may contain:

[0094]

[0095] In the formula: d(·) represents the discriminator activation function; W g b represents the discriminator weight coefficients; g Represents the discriminator bias term; This represents the disturbance category information output by the discriminator.

[0096] Since composite perturbation signals may contain multiple perturbation label information, this invention uses the cross-entropy loss function to measure the similarity between the discriminator output and the actual label of the sample, and achieves training by minimizing the following loss function:

[0097]

[0098] In the formula: K D This represents the number of disturbance categories.

[0099] Therefore, the total loss function of DESAE is defined as:

[0100]

[0101] The working principle of step S303 is that since most of the voltage signal containing disturbances is a normal signal, the features extracted by SAE, which is dominated by minimizing reconstruction error, are prone to lacking discriminative power.

[0102] Step S304: Considering that both the time domain and the frequency domain contain important information reflecting the type of disturbance, the DESAE designed in step S303 simultaneously extracts features from both the time domain and the frequency domain to improve the accuracy and robustness of disturbance identification by fusing time-frequency multimodal features.

[0103] The working principle of step S304 is as follows: the time-domain DESAE and the frequency-domain DESAE are respectively responsible for reconstructing their input signals, and the shared discriminator simultaneously receives the encoder outputs from the time domain and the frequency domain and calculates the possible perturbation information:

[0104]

[0105] H TF =[H T H F ] = [f te (W te ·X T +b te ),f fe (W fe ·X F +b fe )]

[0106] Where: H TF X represents the fused time-frequency domain features; T H T f te (·), W te b te These represent the input time-domain signal of the time-domain DESA, the hidden layer output, the encoder activation function, the encoder weight matrix, and the encoder bias term, respectively; X F H F f fe (·), W fe b fe These are the input frequency domain signal of the frequency domain DESAE, the hidden layer output, the encoder activation function, the encoder weight matrix, and the encoder bias term, respectively.

[0107] Step S305: Guided by the reconstruction loss and the discriminator, the gradient of the time-frequency domain DESAE loss function with respect to the encoder weight coefficients and bias terms is calculated as follows:

[0108]

[0109] In the formula: and These represent the total loss functions of the time-domain DESAE and the frequency-domain DESAE, respectively. As can be seen from the gradient formula, the encoder extracts time-frequency signal features under the combined correction of reconstruction error, KL divergence, and discriminator error.

[0110] Step S306: The gradient of the loss function with respect to the discriminator weight coefficients and bias term is calculated as follows:

[0111]

[0112] Finally, minimizing the loss function of the joint DESAE can be achieved by combining gradient optimization algorithm with backpropagation;

[0113] Step S307: The time-frequency domain DESAE trained through steps S305 and S306 can then perform the time-frequency domain feature extraction task, thereby fusing it with the manually designed indicators obtained in step S2 to obtain a fused feature set, providing a data foundation for subsequent dynamic feature weighting methods.

[0114] Step S4: The dynamic feature weighting method based on dual-delay deep deterministic policy gradient (TD3) and the state space, action space and reward function are used to guide TD3 to learn the optimal dynamic input feature adjustment method. In the offline application stage, the weight of the input features can be adaptively adjusted according to the voltage signal change trend to weight the fused feature set in step S3.

[0115] The specific steps of step S4 include:

[0116] The dynamic feature weighting method based on dual-delay deep deterministic policy gradient (TD3) in step S4 includes the following steps:

[0117] Step S401: Convert the complex vibration signal into a simpler envelope signal through the amplitude envelope to accurately reflect the disturbance of the input signal;

[0118] Therefore, the state space of TD3 is set as follows:

[0119] S={S E ,S F}

[0120] S F =[F1...F 26 H T H F ]

[0121] In the formula: S E S represents the magnitude envelope after the Hilbert transform; F This represents the complete set of perturbation features constructed in this invention.

[0122] The working principle of step S401 is as follows: the variables affecting the control action decision are set as the system's state variables. The state space of the intelligent agent of this invention not only needs to reflect the changing trend of the disturbance signal, but also needs to include the input features to be adjusted.

[0123] Step S402: Design the set of input feature weights This represents the action space of the intelligent agent. Among them, Indicates the Nth F Changes in each hyperparameter; N F This represents the number of input features.

[0124] The working principle of step S402 is as follows: in order to obtain the maximum immediate reward, the intelligent system determines the weight coefficient of each input feature based on the observed system state S.

[0125] Step S403: Reward Function Design. This invention uses the change in classification accuracy before and after an action as the reward function to incentivize the model to perform actions that improve the accuracy of recognizing complex perturbations.

[0126]

[0127] In the formula: and These represent the classification losses in the current state and after the agent's action, respectively; BCE(·) is the cross-entropy loss calculation function; f C (·) represents a classification function; ⊙ represents the Hadamard product operation; A F This represents the actions taken by the intelligent agent.

[0128] The working principle of step S403 is that the classifier can serve as a medium for quantifying the quality of the agent's actions, and provide guidance for adjusting the agent's actions through classification loss.

[0129] Step S404: Based on the action, state, and reward functions designed in steps S401-S403, guide the TD3 reinforcement learning model to dynamically adjust the input feature weights. First, introduce the pruning double Q-learning technique, using two independent Critic networks to calculate Q-values, and select the smaller of the two estimates as the target Q-value:

[0130]

[0131] In the formula: y(t) represents the target Q value; To avoid actions that combine local optima with random noise during the agent's exploration process; The target Critic function; The parameters represent the target Critic network.

[0132] The Q-value in step S404 is interpreted as follows: The Q-value represents the expected reward that the agent can obtain after taking action a in state s. It can also be seen as an evaluation or measure of the merits of taking a specific action in a specific state.

[0133] Step S405: In step S404, a target policy smoothing technique is introduced. Clipped random noise is added to the target actions, making the Critic function smoother across different actions and reducing the likelihood of the policy exploiting Critic function errors. Therefore, It can be represented as:

[0134]

[0135] Where: π λ Represents the target policy function; The parameters represent the target policy network; ε~clip(N(0,σ),-M,M) represents the clipped normal distribution noise, and -M and M represent the upper and lower limits of the clipped noise, respectively.

[0136] Step S406: Based on step S405, a random experience replay method is used to randomly batch-sample data from the experience replay pool as training data. Then, the target Q-value is substituted into the Bellman equation to calculate the temporal difference error and loss function, and the loss function is minimized using the gradient descent algorithm to train the Critic network.

[0137]

[0138] In the formula: Q i Represents the Critic function; N SA The number of state-action pairs involved in training; θ i L(θ) represents the parameters of the Critic network. i ) represents the loss function of the Critic network; η C This represents the learning rate of the Critic network.

[0139] Step S407: During the training process described in step S406, a delayed update strategy was adopted for the Actor network, meaning the update frequency of the Actor network is lower than that of the Critic network. Unlike the Critic network update, the Actor network updates its parameters using a gradient ascent method.

[0140]

[0141] In the formula: θ π The parameters of the Actor network; J(θ) π ) is the objective function of the Actor network; The objective function is relative to θ π The gradient of η; A is the learning rate of the Actor network.

[0142] Step 408: TD3 network parameter update; to avoid error divergence, target network parameters and Approaching θ slowly through a soft update strategy π and θ i :

[0143]

[0144] In the formula: This represents the soft update factor.

[0145] Step S409: The updated TD3 from step S408 can be used to perform adaptive adjustment of feature weights online, assigning weights to the fused feature set from step S3 to provide higher quality input information for the subsequent perturbation identification model.

[0146] Step S5: The XGBoost ensemble learner based on the Boosting strategy is used as the perturbation classifier. The fusion feature set obtained in step S3 is weighted in step S4 and then input into XGBoost to realize the identification of composite perturbation signals.

[0147] The specific steps of step S5 include:

[0148] Step S501: Construct the XGBoost model. Perform a multi-round training process, in which a new CART base classifier is trained in each round to correct the errors in the previous round's predictions. These errors are the differences between the actual perturbation type and the current model's prediction.

[0149] Step S502: In each round, the new CART base classifier is fitted to the error of the previous round, minimizing the objective function, including classification error and regularization term, by adjusting the values ​​of the leaf nodes to prevent overfitting.

[0150]

[0151] In the formula: The cross-entropy loss function represents the type of predicted perturbation. and the true perturbation type y t The error between; f k The classification function represents the k-th base classifier; Ω is the regularization function, expressed as:

[0152]

[0153] In the formula: w represents the total number of leaf nodes in the k-th CART; k λ represents the weight coefficient of the k-th CART leaf node; μ and λ represent the penalty coefficient and regularization coefficient in XGBoost, respectively.

[0154] Step S503: Each CART base classifier uses a greedy algorithm, starting from the root node, to recursively select splitting features and thresholds to maximize information gain. This process continuously splits the dataset, forming multiple leaf nodes. The final XGBoost classification model is a weighted sum of multiple base classifiers.

[0155] Step S504: Power quality disturbance identification using the trained XGBoost model: The input features after dynamically adjusting the weights in step S4 are input into XGBoost. XGBoost outputs the disturbance identification results based on the features, and finally the model is evaluated.

[0156] The invention will be further illustrated below with specific examples:

[0157] To identify composite voltage quality disturbances, a suitable dataset must first be found and preprocessed to provide data support for subsequent numerical analysis. According to the IEEE STD 1159-2019 standard, voltage quality disturbances are classified into seven single disturbance types: swell, transient pulse, sag, flicker, interruption, oscillation, and harmonics. Their general mathematical model can be expressed as:

[0158] Dis=Mul(t)·[sin(ωt)+Add(t)]

[0159] In the formula: ω is the angular frequency of the voltage signal; Mul(t) represents the multiplicative disturbance; Add(t) represents the additive disturbance.

[0160] To verify the effectiveness of the method of this invention, 13 types of signals were classified, namely 7 single disturbances and 6 composite disturbances, namely, voltage sag (S1), voltage dip (S2), voltage interruption (S3), harmonics (S4), oscillation (S5), transient pulse (S6), voltage flicker (S7), voltage sag and harmonic composite disturbance (S8), voltage dip and harmonic composite disturbance (S9), voltage interruption and harmonic composite disturbance (S10), voltage dip and oscillation composite disturbance (S11), voltage flicker and oscillation composite disturbance (S12), and voltage sag, harmonics and oscillation triple composite disturbance (S13). The sample size for each disturbance type was 200, and some disturbance signals are shown below. Figure 2 As shown.

[0161] To verify the discriminative power of this invention, which extracts frequency domain features based on IHHT, in classifying power quality disturbances, the results after ICEEMDAN decomposition are as follows: Figure 3 As shown, the IMF components gradually decrease from high frequency to low frequency, exhibiting the characteristics of different frequency components. Furthermore, the residual components show a smooth trend signal, indicating that the main high-frequency and mid-frequency components have been effectively extracted. In addition, the energy distribution of each IMF is reasonable, and the correlation between them is low, ensuring the independence and orthogonality of the components, thus proving the practicality of the signal decomposition method proposed in this invention. The IHHT method is used to process each perturbation signal, and then feature extraction is performed using the feature index constructed in step S2. The features are then visualized in two-dimensional space using principal component analysis, as shown below. Figure 4As shown, the index designed in this invention effectively distinguishes the characteristics of different disturbance signals in two-dimensional space. Furthermore, Figure 5 The marginal spectral energy shows a clear fluctuation trend, indicating that using it as a feature index is reasonable in this invention. However, the distributions of S1, S9, and S13 in the feature space are still relatively similar, which also verifies the necessity of further feature extraction using DESAE in this invention.

[0162] To demonstrate the superiority of the IHHT proposed in this invention, features extracted by HHT and IHHT for 13 types of disturbance signals were input into the XGBoost model for power quality disturbance identification. The accuracy comparison results for different signal-to-noise ratios are shown in Table 1. IHHT, using ICEEMDAN, yields better signal quality than HHT, exhibiting superior feature quality and noise resistance. Therefore, its classification accuracy is higher than that of HHT under various noise environments.

[0163] Table 1. Comparison of classification accuracy (%) of features obtained from HHT and IHHT under different noise levels.

[0164]

[0165] To demonstrate the effectiveness of the DESAE developed in this invention, the deep features extracted by the autoencoder (AE), sparse autoencoder (SAE), and DESAE were respectively input into the XGBoost model for perturbation identification. The accuracy comparison results at different signal-to-noise ratios are shown in Table 2. It can be seen that the DESAE proposed in this invention exhibits significantly better classification accuracy than traditional methods under various signal-to-noise ratio conditions. At high signal-to-noise ratios (90dB and 80dB), the accuracy of DESAE reached 96.74% and 96.58%, respectively, significantly higher than the corresponding results of AE and SAE; while at lower signal-to-noise ratios (70dB and 60dB), DESAE still maintained excellent performance, with accuracies of 96.46% and 96.14%, respectively. This indicates that DESAE can effectively extract signal features not only under ideal conditions but also maintain robust performance in noisy environments, fully demonstrating the significant superiority of the DESAE developed in this invention in perturbation signal feature extraction and classification.

[0166] Table 2. Comparison of classification accuracy (%) of different deep learning methods

[0167]

[0168]

[0169] Furthermore, examples demonstrate the classification performance of the dynamically weighted XGBoost proposed in this invention, referred to herein as TD3-XGBoost, and its performance is compared with that of Backpropagation Neural Network (BPNN), Convolutional Neural Network (CNN), and Adaptive Boosting (AdaBoost) in Boosting ensemble learning. The input to each classification model is a time-frequency domain feature set obtained by fusing the features obtained in step S2 and step S3. The results of each classifier are as follows: Figure 6 As shown in the figure, BPNN performs well under noise-free conditions, but its accuracy drops significantly with increasing noise, demonstrating its sensitivity to noise. LSTM and CNN perform excellently under high signal-to-noise ratio conditions and have certain noise resistance, but LSTM's resistance to noise is slightly inferior to CNN. AdaBoost has high accuracy at high signal-to-noise ratios, but it is relatively sensitive to noise. The TD3-XGBoost proposed in this invention performs excellently under all signal-to-noise ratio conditions, maintaining a consistently high accuracy, especially under low signal-to-noise ratios, demonstrating its superior noise resistance and robustness, and proving its significant superiority in feature extraction and classification of perturbed signals.

[0170] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0171] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0172] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0173] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

Claims

1. A method for identifying composite voltage quality disturbances considering the fusion of time-frequency features of complex signals, characterized in that: Includes the following steps: Step S1: Obtain the time-frequency domain distribution information of the voltage quality disturbance signal using the improved Hilbert-Huang transform (IHHT) method; Step S2: Based on the time-frequency domain distribution information of the voltage quality disturbance signal obtained in step S1, design disturbance identification input feature indexes and divide these disturbance identification input feature indexes into three categories: time-domain statistical features, frequency-domain statistical features, and energy features. Construct a more comprehensive feature set from the perspectives of amplitude, frequency, phase, and energy. Step S3: Employ a time-frequency domain joint discrimination enhanced sparse autoencoder to extract deep features from voltage quality disturbance information through time-frequency domain joint training. Then, fuse these features with the disturbance identification input feature index manually designed in Step S2 to obtain a more comprehensive fusion feature set. Step S4: The dynamic feature weighting method based on dual-delay deep deterministic policy gradient (TD3) and the state space, action space and reward function are used to guide TD3 to learn the optimal dynamic input feature adjustment method. In the offline application stage, the weight of the input features can be adaptively adjusted according to the voltage signal change trend to weight the fused feature set in step S3. Step S5: The XGBoost ensemble learner based on the Boosting strategy is used as the perturbation classifier. The fusion feature set obtained in step S3 is weighted in step S4 and then input into XGBoost to realize the identification of composite perturbation signals. The specific steps of step S2 include: Step S201: Design the time-domain statistical characteristics as: maximum voltage amplitude, minimum voltage amplitude, average voltage amplitude, standard deviation of voltage, peak factor of voltage amplitude, and duration of disturbance; Step S202: After the Hilbert transform, the instantaneous characteristics of each IMF can be obtained, and the following time-domain statistical characteristics are added: maximum instantaneous amplitude of each IMF, minimum instantaneous amplitude of each IMF, average instantaneous amplitude of each IMF, standard deviation of instantaneous amplitude of each IMF, and standard deviation of instantaneous phase of each IMF. Step S203: Design the frequency domain statistical features as follows: the top three amplitudes in the marginal spectrum and their corresponding frequencies, the minimum marginal spectrum amplitude, the average marginal spectrum amplitude, the standard deviation of the marginal spectrum amplitude, the maximum high-frequency marginal spectrum amplitude, the minimum high-frequency marginal spectrum amplitude, the average high-frequency marginal spectrum amplitude, and the standard deviation of the high-frequency marginal spectrum amplitude. Step S204: Design the energy characteristics as: IMF energy and total marginal spectrum energy; IMF Energy E imf This reflects the strength of the IMF, and the specific calculation method is as follows: Where: IMF k (i) can be the i-th sampling point in the k-th IMF; N is the total number of sampling points; The marginal spectrum is obtained by integrating the Hilbert spectrum over the time axis, reflecting the total energy of the signal at each frequency. Therefore, the total marginal spectral energy is defined as: In the formula: N represents the number of sampling points of the marginal spectrum curve; h(ω,i) represents the i-th sampling point of the marginal spectrum curve h(ω); The specific steps of step S4 include: Step S401: Convert the complex vibration signal into a simpler envelope signal through the amplitude envelope to accurately reflect the disturbance of the input signal; Therefore, the state space of TD3 is set as follows: S={S E ,S F } S F =[F1...F 26 ,H T ,H F ] In the formula: S E S represents the magnitude envelope after the Hilbert transform; F This represents the complete set of perturbation features constructed in this invention; Step S402: Design the set of input feature weights This is the action space of the intelligent agent; where, Indicates the Nth F Changes in each hyperparameter; N F Represents the number of input features; Step S403: Reward function design: Use the change in classification accuracy before and after the action as the reward function to incentivize the model to perform actions that improve the accuracy of recognizing composite perturbations; In the formula: and These represent the classification losses in the current state and after the agent's action, respectively; BCE(·) is the cross-entropy loss calculation function; f C (·) represents a classification function; ⊙ represents the Hadamard product operation; A F Represents the actions taken by the intelligent agent; Step S404: Based on the action, state, and reward functions designed in steps S401-S403, guide the TD3 reinforcement learning model to dynamically adjust the input feature weights. First, introduce the pruning double Q learning technique, using two independent Critic networks to calculate the Q-value, and select the smaller of the two estimates as the target Q-value: In the formula: y(t) represents the target Q value; To avoid actions that combine local optima with random noise during the agent's exploration process; The target Critic function; Representative Critic network parameters; The Q value in step S404 is interpreted as follows: the Q value represents the expected reward that the agent can obtain after taking action a in state s; it is regarded as an evaluation or measure of the quality of taking an action in a certain state. Step S405: In step S404, a target policy smoothing technique is introduced. Clipped random noise is added to the target actions, making the Critic function smoother across different actions and reducing the possibility of the policy exploiting Critic function errors. Therefore... Represented as: Where: π λ Represents the target policy function; The parameters represent the target policy network; ε~clip(N(0,σ),-M,M) represents the clipped normal distribution noise, and -M and M represent the upper and lower limits of the clipped noise, respectively. Step S406: Based on step S405, the random experience replay method is used to randomly draw data in batches from the experience replay pool as training data. Then, the target Q value is substituted into the Bellman equation to calculate the temporal difference error and loss function, and the loss function is minimized through the gradient descent algorithm to train the Critic network. In the formula: Q i Represents the Critic function; N SA The number of state-action pairs involved in training; θ i L(θ) represents the parameters of the Critic network. i ) represents the loss function of the Critic network; η C This represents the learning rate of the Critic network. Step S407: During the training process in step S406, a delayed update strategy was adopted for the Actor network, meaning the update frequency of the Actor network is lower than that of the Critic network. Unlike the Critic network update, the Actor network updates its parameters using a gradient ascent method. In the formula: θ π The parameters of the Actor network; J(θ) π ) is the objective function of the Actor network; The objective function is relative to θ π The gradient of η; A is the learning rate of the Actor network; Step 408: TD3 network parameter update; to avoid error divergence, target network parameters and Approaching θ slowly through a soft update strategy π and θ i : In the formula: Represents the soft update factor; Step S409: The updated TD3 from step S408 can be used to perform adaptive adjustment of feature weights online, assigning weights to the fused feature set from step S3 to provide higher quality input information for the subsequent perturbation identification model.

2. The voltage quality composite disturbance identification method considering the fusion of time-frequency features of complex signals according to claim 1, characterized in that: The specific steps of step S1 include: Step S101: Inject the modal components of the i-th group of white noise into the original time series S using ICEEMDAN to construct a new time series X. (i) : X (i) =S+β0E1(ω (i) ) β0=ε0σ(S) / σ(E1(ω (i) )) In the formula: β0 is the noise figure, σ(·) is the standard deviation operator, and E i (·) represents the i-th mode obtained from EMD decomposition, ω (i) Let the i-th group of white noise (i = 1, 2, ..., I) follow an N(0, 1) distribution, where I is the total number of modes decomposed, and ε0 is the amplitude of the white noise; Step S102: Calculate the local mean to obtain the residual R1 of the first decomposition. R1= <M(X (i) ) In the formula: <·> is the average value calculation, and M(·) is the local mean value calculation function of the electrical energy signal; Step S103: When the modal decomposition order of the sampled signal is 1, calculate the first modal component IMF1 using the original signal S and the residual signal R1, that is: IMF1 = S-R1 Step S104: Add white noise to the residual to construct a new signal to be decomposed, R1+β1E2(ω). (i) And calculate the second set of residuals R2 and its modal components IMF2: IMF2=R1-R2=R1-<M(R1+β1E2(ω (i) ))> Step S105: Calculate the residual R of the k-th decomposition, and so on. k and its modal components IMF k : R k = <M(R k-1 +b k-1 E k oh (i) )> IMF k =R k-1 -R k In the formula: when k≥1, β k =ε0σ(R) k ) Step S106: Repeat step S105 to obtain all modal components and the final residual; Step S107: Using the frequency domain correlation coefficient as the evaluation criterion, select the best IMF component with clear features from a series of IMF components obtained in step S106 for signal reconstruction and time-frequency domain feature extraction. Step S108: Based on the dominant IMF obtained in step S107, the local amplitude and frequency information existing in the processed signal is obtained through Hilbert transform; The IMF after Hilbert transformation can be represented as: In the formula: H(·) is the Hilbert transform function; P is the Cauchy principal value; IMF k (t) represents the value of the k-th eigenmode function at time t; τ represents the integration variable; Step S109: Construct the analytic signal Z based on the Hilbert transform result obtained in step S108. k (t)=IMF k (t)+jH[IMF k [t], the magnitude function a of each IMF k (t) and phase function θ k (t) is calculated as follows: In the formula: Z k (t) represents the analytical signal corresponding to the k-th IMF; arctan is the arctangent function; Step S110: Based on the phase function θ obtained in step S109 k (t), calculate the instantaneous frequency, Hilbert spectrum, and final marginal spectrum of IHHT to obtain the time-frequency domain distribution information of the voltage quality disturbance signal, in order to design subsequent feature indicators: In the formula: T is the length of the signal; ω k H(ω,t) represents the instantaneous frequency of the k-th IMF at time t; H(ω,t) represents the Hilbert spectrum, which accurately reflects the trend of amplitude variation with duration and frequency; h(ω) represents the marginal spectrum, which reflects the contribution of each frequency to the amplitude; and K represents the number of intrinsic mode functions.

3. The voltage quality composite disturbance identification method considering the fusion of time-frequency features of complex signals according to claim 1, characterized in that: The specific steps of step S3 include: Step S301: Based on the encoding function of the autoencoder AE, the input signal is converted into a compressed feature representation through weighting and biasing: H=f e (W e ·X+b e ) X R =f d (W d ·H+b d ) In the formula: W e and b e W represents the weights and biases from the input layer to the hidden layer, respectively. d and b d X, H, and X represent the weights and biases from the hidden layer to the output layer, respectively. R These represent the original input, hidden layer output, and reconstructed data, respectively; f e (·) and f d (·) is the activation function; J AE X represents the autoencoder reconstruction error; i and These represent the i-th original input and the reconstructed data, respectively. Step S302: Based on step S301, construct a sparse autoencoder (SAE). Sparsity constraints affect the encoder's behavior by adding a sparsity penalty term to the loss function. KL divergence is used to measure the difference between the actual activation distribution of the hidden layer and the target sparsity. In the formula: ρ represents the sparsity parameter of the KL divergence; represents the average activation of hidden layer node j; the purpose of the penalty term is to make active hidden nodes close to 1 and inactive nodes close to 0, forcing the model to learn more generalizable feature representations. Furthermore, the autoencoder loss function after introducing sparsity constraints is expressed as: In the formula: J SAE (W,b) is the loss function after introducing sparsity constraints; β is the weighting coefficient of KL divergence; Step S303: Based on step S302, add a discriminator module to construct a class label-assisted discriminative enhanced sparse autoencoder (DESAE). During training, the hidden layer features of the SAE are input into the discriminator, and the class label information helps to enhance the attention to the perturbation signal. Assuming a given training sample {X,Y}, where Y represents the perturbation class to which the sample belongs, the discriminator calculates the perturbation class information it may contain: In the formula: d(·) represents the discriminator activation function; W g b represents the discriminator weight coefficients; g Represents the discriminator bias term; This represents the disturbance category information output by the discriminator; The cross-entropy loss function is used to measure the similarity between the discriminator's output and the actual labels of the samples. Training is achieved by minimizing the following loss function: Where: K D This represents the number of disturbance categories; The total loss function of DESAE is defined as: Step S304: Considering that both the time domain and the frequency domain contain important information reflecting the type of disturbance, the DESAE designed in step S303 simultaneously extracts features from both the time domain and the frequency domain to improve the accuracy and robustness of disturbance identification by fusing time-frequency multimodal features. The time-domain DESAE and frequency-domain DESAE are responsible for reconstructing their respective input signals. The shared discriminator simultaneously receives encoder outputs from both the time and frequency domains and calculates any potential perturbation information. H TF =[H T ,H F ]=[f te (W te ·X T +b te ),f fe (W fe ·X F +b fe )] In the formula: H TF X represents the fused time-frequency domain features; T H T f te (·), W te b te These represent the input time-domain signal of the time-domain DESA, the hidden layer output, the encoder activation function, the encoder weight matrix, and the encoder bias term, respectively; X F H F f fe (·), W fe b fe These are the input frequency domain signal of the frequency domain DESAE, the hidden layer output, the encoder activation function, the encoder weight matrix, and the encoder bias term, respectively. Step S305: Guided by the reconstruction loss and the discriminator, the gradient of the time-frequency domain DESAE loss function with respect to the encoder weight coefficients and bias terms is calculated as follows: In the formula: and These represent the total loss functions of the time-domain DEESAE and the frequency-domain DEESAE, respectively. As can be seen from the gradient formula, the encoder extracts time-frequency signal features under the joint correction of reconstruction error, KL divergence, and discriminator error. Step S306: The gradient of the loss function with respect to the discriminator weight coefficients and bias term is calculated as follows: Finally, the minimization of the loss function of the joint DESAE is achieved by combining gradient optimization algorithm with backpropagation; Step S307: The time-frequency domain DESAE trained through steps S305 and S306 can then perform the time-frequency domain feature extraction task, thereby fusing it with the manually designed indicators obtained in step S2 to obtain a fused feature set, providing a data foundation for subsequent dynamic feature weighting methods.

4. The voltage quality composite disturbance identification method considering the fusion of time-frequency features of complex signals according to claim 1, characterized in that: The specific steps of step S5 include: Step S501: Construct the XGBoost model; perform a multi-round training process, in which a new CART base classifier is trained in each round to correct the errors in the previous round's predictions; these errors are the differences between the actual perturbation type and the current model's prediction. Step S502: In each round, the new CART base classifier is fitted to the error of the previous round, minimizing the objective function, including classification error and regularization term, by adjusting the values ​​of the leaf nodes to prevent overfitting. In the formula: The cross-entropy loss function represents the type of predicted perturbation. and the true perturbation type y t The error between; f k The classification function represents the k-th base classifier; Ω is the regularization function, expressed as: In the formula: w represents the total number of leaf nodes in the k-th CART; k λ represents the weight coefficient of the k-th CART leaf node; μ and λ represent the penalty coefficient and regularization coefficient in XGBoost, respectively. Step S503: Each CART base classifier uses a greedy algorithm to recursively select splitting features and thresholds starting from the root node to maximize information gain; by continuously splitting the dataset to form multiple leaf nodes, the final XGBoost classification model is a weighted sum of multiple base classifiers; Step S504: Power quality disturbance identification using the trained XGBoost model: The input features after dynamically adjusting the weights in step S4 are input into XGBoost. XGBoost outputs the disturbance identification results based on the features, and finally the model is evaluated.

Images