Method for identifying power quality disturbance characteristics of distributed photovoltaic interference power supply
By combining the improved variational mode decomposition and Hamiltonian dynamics model with meteorological data prediction methods, the accuracy and reliability problems of identifying the power quality disturbance characteristics of photovoltaic systems in the existing technology are solved, and efficient disturbance identification and early warning of distributed photovoltaic systems are achieved.
Patent Information
- Application Number
- CN202510548134.5
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-04-28
- Publication Date
- 2025-09-12
AI Technical Summary
When identifying the power quality disturbance characteristics of distributed photovoltaic systems, existing technologies are unable to take into account the complex signal characteristics and the influence of meteorological factors, resulting in insufficient monitoring and early warning accuracy and reliability, and the inability to accurately identify short-term complex disturbances.
An improved variational mode decomposition method is used in combination with information geometry optimization and Bayesian optimization for adaptive decomposition. The Hamiltonian dynamics model and Lyapunov index are used to analyze the nonlinear characteristics of the disturbance signal. The ARIMA and Transformer-LSTM hybrid model is used in combination with meteorological data for short-term and long-term forecasts.
It achieves accurate decomposition and identification of disturbance signals of photovoltaic systems, improves the interpretability and reliability of power quality monitoring, enhances the accuracy of disturbance warnings, and supports intelligent scheduling of photovoltaic systems.
Smart Images

Figure CN120632803A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of power quality monitoring, and in particular to a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply. Background Art
[0002] The widespread adoption of distributed photovoltaic systems has led to increasingly prominent power quality issues. Especially during the grid connection process, photovoltaic power sources can generate various disturbances, such as voltage fluctuations, frequency variations, and harmonic distortion. These disturbances can affect power quality and lead to system instability. Therefore, accurately identifying and providing early warning of disturbance characteristics is crucial for ensuring grid security and improving the operational reliability of power systems.
[0003] During the operation of the current system, disturbances such as low-frequency harmonics, short-term fluctuations, and instantaneous drops often occur. Existing technologies mainly rely on fixed parameter decomposition and traditional statistical models, which cannot take into account the influence of complex signal characteristics and meteorological factors, restricting the accuracy and reliability of monitoring and early warning.
[0004] Currently, most solutions use fixed-parameter variational modal decomposition. The parameters need to be set manually and cannot be adjusted dynamically. They cannot cope with the diversity of disturbance signals and the decomposition effect is unstable.
[0005] Existing technologies mostly rely on traditional Fourier transform or empirical mode decomposition when modeling disturbance signals, which cannot capture the nonlinear chaotic characteristics in the signal, lack physical connotations, and cannot accurately identify short-term complex disturbances.
[0006] Common prediction methods only use statistical regression or a single neural network model, ignoring the impact of environmental meteorological data on disturbances. The accuracy of both short-term and long-term predictions is limited, and the early warning response is delayed.
[0007] Therefore, those skilled in the art provide a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply to solve the above-mentioned problems. Summary of the Invention
[0008] In view of the deficiencies in the prior art, the present invention provides a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply to solve the problems raised in the above background technology.
[0009] To achieve the above objectives, the present invention is implemented through the following technical solutions: a method for identifying power quality disturbance characteristics of distributed photovoltaic interference power supply, comprising:
[0010] Step S1, data acquisition and preprocessing, collecting the power quality signal of the photovoltaic system and performing denoising, normalization and time alignment;
[0011] Step S2, adaptively decomposing the disturbance signal based on improved variational modal decomposition, automatically selecting the number of modes using information geometry optimization, and selecting the penalty factor through Bayesian optimization;
[0012] Step S3, performing disturbance characteristic modeling based on the Hamiltonian dynamics model, constructing a Hamiltonian system and calculating the phase space dynamic equation of the disturbance signal;
[0013] Step S4, calculating the Lyapunov exponent of the Hamiltonian dynamics equation established in step S3, performing nonlinear characteristic analysis of the disturbance signal based on the Lyapunov exponent, calculating the maximum Lyapunov exponent of the disturbance signal and determining its chaotic characteristics;
[0014] Step S5: Combined with the disturbance prediction of meteorological data, the ARIMA model is used for short-term disturbance prediction, and the Transformer-LSTM hybrid model is used for long-term disturbance prediction;
[0015] Step S6: Output and apply the results, output the disturbance identification results and issue a real-time warning.
[0016] Preferably, in step S1, data collection and preprocessing further include:
[0017] Step 1.1, collect power quality signals:
[0018] When collecting power quality signals, use data acquisition equipment to record the three-phase voltage u output by the photovoltaic system a (t),u b (t),u c (t) and three-phase current i a (t), i b (t), i c (t), and simultaneously record the power quality parameters on the grid side;
[0019] Let the sampling frequency be f s , the sampling time is T, then the discrete form of the collected signal is expressed as:
[0020] u n [k]=u n (kT s ), i n [k]=i n (kT s ), k=0,1,…,N-1,
[0021] Among them, u n [k] is the discrete value of the voltage of the nth phase at the kth sampling point, i n [k] is the discrete value of the n-th phase current at the k-th sampling point, T s =1 / f sis the sampling period, N=Tf s is the total number of sampling points;
[0022] Step 1.2, signal denoising:
[0023] Wavelet transform is used to decompose the signal at multiple scales to remove high-frequency noise components;
[0024] Assuming the number of wavelet decomposition layers is J and the mother wavelet function is ψ(t), the wavelet transform coefficient of the signal is expressed as:
[0025]
[0026] Among them, W ψ (a, b) is the wavelet transform coefficient of signal u(t) at scale a and position b, ψ(t) is the mother wavelet function, ψ * (t) is the conjugate complex number of the mother wavelet function, a is the scale factor, and b is the translation factor;
[0027] The denoising process consists of the following steps:
[0028] Calculate the wavelet transform coefficient W of the signal ψ (a,b);
[0029] Set threshold λ j , perform soft threshold filtering on high-frequency wavelet coefficients:
[0030]
[0031] Among them, λ j is the threshold of the j-th layer wavelet decomposition, is the update coefficient after soft threshold processing;
[0032] Perform wavelet reconstruction to obtain the denoised signal
[0033] Step 1.3, signal normalization:
[0034] Using the min-max normalization method, the signal is mapped to the interval [1, 2]:
[0035]
[0036] Among them, u(t) is the original voltage signal and i(t) is the original current signal.
[0037] u min is the minimum value of the voltage signal, u max is the maximum value of the voltage signal, i min is the minimum value of the current signal, i max is the maximum value of the current signal, is the normalized voltage signal, is the normalized current signal;
[0038] Step 1.4, signal time alignment:
[0039] The time delay between signals is calculated by the cross-correlation function. The calculation formula of the cross-correlation function is:
[0040]
[0041] Among them, R ui (τ) is the voltage signal u a [k] and current signal i a [k] The cross-correlation function value between u a [k] is the kth discrete value of the ath phase voltage signal, i a [k] is the kth discrete value of the ath phase current signal, and τ is the time delay;
[0042] Find the τ corresponding to the maximum mutual correlation coefficient max , and make corresponding time offset adjustments to the signals so that all signals are aligned to the same time base.
[0043] Preferably, in step S2, the adaptive decomposition of the disturbance signal based on improved variational mode decomposition further includes:
[0044] Step 2.1, mathematical modeling of variational mode decomposition:
[0045] Assume that the preprocessed disturbance signal is expressed as x(t):
[0046]
[0047] Among them, u k (t) is the kth group of modal components, x(t) is the preprocessed disturbance signal, and K is the target mode number;
[0048] Each modal component u k (t) at a specific center frequency ω k The neighborhood has limited bandwidth;
[0049] Variational mode decomposition achieves signal decomposition by solving the following variational optimization problem:
[0050]
[0051] Among them, H(u k ) is the modal component u k (t) is the signal after Hilbert transform, j is the imaginary unit, ω k is the center frequency of the modal component, is the time derivative operator, ui (t) is the modal component of group i;
[0052] The constraints are introduced using the Lagrange multiplier method to construct the Lagrange function:
[0053]
[0054] Among them, L({u k},{ω k},λ) is the Lagrangian function of variational mode decomposition, α is the penalty factor, λ is the Lagrangian multiplier, and <·,·> is the inner product operation;
[0055] Step 2.2, information geometry optimization automatically selects the number of modes:
[0056] Based on the information geometry optimization method, the spectral entropy H of each modal component is maximized. k Calculate the optimal number of modes:
[0057] in, is the modal component u k (t) at frequency f n Normalized energy distribution at :
[0058]
[0059] Let the total information entropy H total As the objective function:
[0060]
[0061] Traverse different modal numbers K and find the total The maximum optimal modal number K * :
[0062]
[0063] Among them, H k is the spectral entropy of the kth modal component, is the modal component u k (t) at frequency f n Normalized energy distribution at , U k (f n ) is the modal component u k (t) at frequency f n Fourier transform at, N is the number of frequency points, H total is information entropy, K * is the optimal mode number;
[0064] Step 2.3, Bayesian optimization to determine the penalty factor: Use the Bayesian optimization method to decompose the reconstruction error E of the signal rec As the objective function:
[0065] Gaussian process regression model is used to predict E under different α values rec , and select the optimal penalty factor based on the expected improvement criterion:
[0066] Among them, E rec is the signal reconstruction error, α is the penalty factor, is the reconstructed signal of all modal components after variational mode decomposition, To calculate the penalty factor α that minimizes the reconstruction error * ;
[0067] Step 2.4, calculate the decomposed disturbance signal:
[0068] Based on the determined K * and α * , perform the final variational mode decomposition on the preprocessed disturbance signal to obtain K * Group modal components:
[0069] Among them, each u k (t) represents the different frequency components in the photovoltaic disturbance signal;
[0070] At this point, the adaptive decomposition of the disturbance signal based on the improved variational mode decomposition is completed, and the decomposed disturbance signal u is obtained k (t),u k (t) will be used to construct the Hamiltonian dynamics model in step S3.
[0071] Preferably, in step S3, the disturbance characteristic modeling based on the Hamiltonian dynamics model further includes:
[0072] After completing the adaptive decomposition of the disturbance signal in step S2, we get K * Group modal component u k (t), the modal components represent the disturbance signals of different frequency bands in the photovoltaic system and will be used to construct the Hamiltonian dynamics model in step S3;
[0073] Step 3.1, Hamiltonian dynamics modeling:
[0074] Assume the phase space of the disturbance signal is (q k (t),p k (t)), where q k (t) and p k (t) represents the modal component uk (t) position and momentum;
[0075]
[0076] Among them, m k is the modal component u k (t) mass, V(q k ) is the potential energy function;
[0077] For the disturbance signal, the momentum p k (t) and position q k (t) satisfies the following Hamiltonian equation:
[0078]
[0079] in, and Indicates speed and force;
[0080] Step 3.2, construct the phase space model of the disturbance signal:
[0081] By calculating the modal components u k (t) is numerically integrated to obtain the motion trajectory in the phase space;
[0082] Assume that each modal component u k The potential energy function V(q k ) is a simple quadratic potential energy:
[0083]
[0084] Among them, k k is the modal component u k (t) elastic constant;
[0085] By solving the Hamiltonian equation, we obtain the dynamic equation in phase space:
[0086]
[0087] Among them, q k (t) and p k (t) represents the modal component q k (t),p k (t) position and momentum in phase space;
[0088] Step 3.3, calculate the phase space dynamic equation of the disturbance signal:
[0089] The classic fourth-order Runge-Kutta algorithm is used to numerically integrate the Hamiltonian equation, setting the time step Δt and iteratively updating q k (t) and pk The value of (t):
[0090]
[0091] p k (t+Δt)=p k (t)-Δt·k k q k (t),
[0092] Step 3.4: Combine the disturbance signal to perform Hamiltonian dynamics modeling analysis:
[0093] After step S3 is completed, the dynamic characteristics of the disturbance signal have been modeled by the Hamiltonian dynamics model, and then step S4 is entered to further analyze the nonlinear characteristics of the disturbance signal by using the Lyapunov exponent.
[0094] Preferably, in step S4: the nonlinear characteristic analysis of the disturbance signal based on the Lyapunov exponent further includes:
[0095] In step S3, a Hamiltonian dynamics model of the disturbance signal is constructed, and the phase space dynamic equation of the disturbance signal is obtained. To further analyze the nonlinear characteristics and chaotic characteristics of the disturbance signal, in step S4, the Lyapunov exponent needs to be calculated, and the chaotic characteristics of the disturbance signal are judged based on the maximum Lyapunov exponent.
[0096] Step 4.1, mathematical definition of Lyapunov exponent:
[0097] In the phase space, the phase trajectory of the disturbance signal is described by the Hamiltonian dynamics equation, and the state vector is: X(t) = [q1(t), p1(t), q2(t), p2(t), ..., q K* (t),p K* (t)] T ,
[0098] Among them, q K (t) and p K (t) is the modal component u k (t) is the phase space coordinate, K * is the optimal number of modal components;
[0099] The time evolution of the disturbance signal in phase space is described by the following Hamiltonian dynamics equation:
[0100]
[0101] Where, F(X(t)) is the evolution equation of the system;
[0102] At the initial time t0, adjacent trajectory points X(t0) and X′(t0) in the phase space are selected, and the initial small perturbation is: δX(t0)=X′(t0)-X(t0),
[0103] The evolution after time t satisfies:
[0104] Where λ is the Lyapunov exponent, which describes the exponential growth rate of the disturbance between adjacent trajectories. The rate of change of the distance between trajectories is given by the following formula:
[0105] Step 4.2, calculate the maximum Lyapunov exponent:
[0106] To calculate the maximum Lyapunov exponent, we need to perform numerical simulation on the phase space trajectory and use the small perturbation method to track the divergence of adjacent trajectories. The specific calculation steps are as follows:
[0107] Step 4.2.1, in the phase space trajectory X(t), select the initial point X(t0), and select the initial small perturbation δX(t0) in the neighborhood to satisfy: ||δX(t0)||=∈0,
[0108] Among them, ∈0 is the initial disturbance amplitude;
[0109] Step 4.2.2, solve the Hamiltonian dynamics equation by numerical integration method and calculate the distance between the perturbed trajectory X′(t) and the original trajectory X(t): d(t) = ||X′(t) - X(t)||,
[0110] Among them, d(t) reflects the changing trend of the disturbance trajectory over time.
[0111] Step 4.2.3: In the logarithmic coordinate system, plot the curve of ln(t) changing with time t, and use the linear regression method to solve the slope λ max :
[0112] Among them, N t is the number of time steps;
[0113] Step 4.3: Determine the chaotic characteristics based on the maximum Lyapunov exponent:
[0114] According to the maximum Lyapunov exponent λ max The numerical value of is used to judge the chaotic characteristics of the disturbance signal:
[0115] If λ max >0, the disturbance signal has chaotic characteristics;
[0116] If λ max =0, the system belongs to quasi-periodic motion;
[0117] If λ max <0, then the system has an attractor;
[0118] Step 4.4, calculate multiple groups of Lyapunov exponents:
[0119] Except for the maximum Lyapunov exponent λ max In addition, calculate other Lyapunov exponents λ i To further analyze the stability of the disturbance signal, the complete Lyapunov exponent spectrum Λ consists of multiple groups of exponents:
[0120] Λ=[λ1,λ2,…,λ 2K* ],
[0121] Where Λ is the optimal mode number;
[0122] Using the QR decomposition method, calculate the Lyapunov index spectrum and first construct the Jacobian matrix J(X) of the perturbation equation: Then, the exponential divergence rate of the orthogonal basis is tracked by QR decomposition to calculate all Lyapunov exponents;
[0123] Step 4.5, result output and next step prediction: After completing the Lyapunov exponent calculation, the nonlinear characteristics of the disturbance signal are obtained:
[0124] If λ max >0, it indicates that the disturbance signal has chaotic characteristics, and a prediction model that adapts to chaotic characteristics is used.
[0125] If λ max ≤0, the disturbance signal exhibits periodic and convergent behavior, and the traditional time series prediction model is used;
[0126] At this point, step S4 is completed, the Lyapunov exponent of the disturbance signal is successfully extracted, and the nonlinear characteristics are analyzed. Next, step S5 is entered to perform short-term and long-term predictions on the disturbance signal based on meteorological data.
[0127] Preferably, in step S5, the disturbance prediction based on meteorological data further includes:
[0128] In step S4, the nonlinear characteristics of the disturbance signal are analyzed by calculating the Lyapunov exponent and the chaotic behavior is determined. In step S5, the disturbance signal is predicted in the short and long term in combination with the meteorological data of the photovoltaic system.
[0129] Step 5.1, short-term disturbance forecast based on ARIMA model:
[0130] Assume that the time series of the disturbance signal is: s(t)=[s1,s2,...,s N ], where s Nis the disturbance signal value of N time steps, N is the total number of observation time points;
[0131] Step 5.1.1, establishment of ARIMA model:
[0132] The ARIMA model is determined by the parameters (p, d, q):
[0133] p is the autoregressive order, which indicates the influence of the value at the previous p moments on the current value;
[0134] d is the difference order, which indicates the number of differences required to stabilize the time series;
[0135] q is the moving average order, which represents the influence of the first q error terms on the current value;
[0136] The mathematical expression of the ARIMA model is:
[0137] Φ p (B)(1-B) d s(t)=Θ q (B)∈(t),
[0138] Among them, B is the backshift operator, Φ p (B) is the polynomial of the autoregressive part, Θ q (B) is the polynomial of the moving average part, ∈(t) is the white noise term;
[0139] Step 5.1.2, use AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) to optimize and select the optimal (p, d, q):
[0140] AIC=-2lnL+2k, BIC=-2lnL+klnN,
[0141] Where L is the log-likelihood function of the model, k = p + q + d is the number of model parameters, and N is the total number of observation time points;
[0142] The optimal parameter (p * ,d * ,q * ) is determined by the following optimization objectives:
[0143]
[0144] Step 5.1.3, predict disturbance signals: Based on the optimal ARIMA model, predict short-term disturbance signals:
[0145]
[0146] Among them, h is the prediction step size, is the predicted disturbance signal value at the future time step t+h, φi is the autoregressive coefficient, θ j is the moving average coefficient.
[0147] Preferably, in step S5, the disturbance prediction based on meteorological data further includes:
[0148] Step 5.2, long-term disturbance prediction is based on the Transformer-LSTM hybrid model:
[0149] Due to the chaotic characteristics of the disturbance signal, a deep learning model that can capture long-term dependencies is adopted, namely the Transformer-LSTM hybrid model;
[0150] Step 5.2.1, LSTM network modeling. The LSTM structure consists of an input gate, a forget gate, and an output gate. The calculation formula is as follows:
[0151] Forget gate: f t =σ(W f ·[h t-1 ,x t ]+b f )
[0152] Among them, x t is the current input, h t-1 is the hidden state at the previous moment, W f and b f are the weight and bias of the forget gate, σ(·) is the Sigmoid function;
[0153] Input gate: i t =σ(W i ·[h t-1 ,x t ]+b i ),
[0154]
[0155] Among them, i t is the output of the input gate, W i is the weight matrix of the input gate, b i is the bias vector of the input gate, is the candidate cell state, W C is the weight matrix of the candidate cell state, b C is the bias vector of the candidate cell state;
[0156] Cell status update:
[0157] Among them, C t is the current cell state, C t-1 is the cell state at the previous moment;
[0158] Output gate: o t =σ(W o ·[h t-1 ,x t ]+b o ),
[0159] h t =o t tanh(C t ),
[0160] Among them, t is the output of the output gate, W o is the weight matrix of the output gate, b o is the bias vector of the output gate, h t is the final output of LSTM;
[0161] Step 5.2.2, Transformer processes global features: Transformer uses the self-attention mechanism to extract global features. Given the input perturbation signal sequence S, define the query Q, key K and value V:
[0162] Q=W Q S, K = W K S, V = W V S,
[0163] Among them, W Q 、W K 、W K is the linear transformation matrix;
[0164] Calculate attention weights:
[0165] Among them, d k is the dimension of the key vector;
[0166] Step 5.2.3, predict the disturbance signal: use Transformer to extract long-term dependency features and use LSTM to predict:
[0167]
[0168] Among them, f Transformer is the feature extracted by Transformer, f LSTM is the LSTM prediction function.
[0169] Preferably, in step S5, the disturbance prediction based on meteorological data further includes:
[0170] Step 5.3, perform disturbance correction based on meteorological data:
[0171] The disturbance signal of the photovoltaic system is affected by meteorological factors. Meteorological data include:
[0172] Solar radiation intensity G(t), temperature T(t), wind speed T(t) and cloud cover T(t);
[0173] Construct a disturbance correction model:
[0174]
[0175] Among them, β1, β2, β3, and β4 are the weights of meteorological factors. is the forecast value after the meteorological data is corrected, is the original predicted value;
[0176] Step 5.4, result output and application:
[0177] Short-term forecasting is used for real-time monitoring of photovoltaic systems;
[0178] Long-term forecasts are used for power grid dispatching;
[0179] The prediction results will be used for real-time warning in step S6 to prevent the impact of photovoltaic disturbances on the power grid. At this point, step S5 is completed, and the short-term and long-term evolution trends of the disturbance signal are successfully predicted.
[0180] A terminal device includes a memory, a processor, and a computer program stored in the memory and running on the processor, wherein the computer program is configured to execute a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply.
[0181] A storage medium stores a computer program, which, when executed by a processor, implements a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply.
[0182] The present invention provides a method for identifying the power quality disturbance characteristics of a distributed photovoltaic interference power supply. It has the following beneficial effects:
[0183] 1. The present invention adopts an improved variational modal decomposition method, combined with information geometry optimization and Bayesian optimization, to achieve adaptive modal decomposition of disturbance signals. The improved variational modal decomposition can automatically determine the optimal number of modes and penalty factors, avoiding the decomposition error caused by manual parameter setting in traditional variational modal decomposition. Compared with the fixed-parameter variational modal decomposition method in the prior art, the present invention improves the analytical ability of complex disturbance signals and accurately decomposes the disturbance modes in photovoltaic systems.
[0184] 2. The present invention adopts Hamiltonian dynamics modeling to model and analyze photovoltaic disturbance signals from the perspective of energy conservation. It further combines Lyapunov exponent calculation to extract the chaotic characteristics of the signal, which can accurately identify short-term composite disturbances and effectively distinguish disturbances from systematic disturbances. Compared with Fourier analysis or empirical mode decomposition in the prior art, the dynamic modeling method of the present invention can intuitively reveal the physical characteristics of the disturbance signal and improve the interpretability and reliability of power quality monitoring.
[0185] 3. The present invention combines the meteorological data of the photovoltaic system and adopts a hybrid prediction model to achieve short-term and long-term prediction of disturbance signals. The ARIMA model is suitable for short-term trend prediction, while the Transformer-LSTM can capture long-term evolution patterns. The combination of the two can effectively improve the accuracy of disturbance warning. Compared with the existing prediction methods based on statistical regression or single LSTM, the present invention can accurately predict the power quality disturbances caused by light fluctuations and meteorological mutations, making the intelligent scheduling of the photovoltaic system efficient. BRIEF DESCRIPTION OF THE DRAWINGS
[0186] Figure 1 Flowchart of the present invention. DETAILED DESCRIPTION
[0187] To help those skilled in the art understand the present invention, the following will provide a clear and complete description of the technical solutions in the embodiments of the present invention, in conjunction with the accompanying drawings. Obviously, the described embodiments are only partial embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.
[0188] The present invention is described in detail below with reference to the accompanying drawings:
[0189] Example:
[0190] Please see the attached Figure 1 The embodiment of the present invention provides a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply, comprising:
[0191] Step S1, data acquisition and preprocessing, collecting the power quality signal of the photovoltaic system and performing denoising, normalization and time alignment;
[0192] Step 1.1, collect power quality signals:
[0193] When collecting power quality signals, use data acquisition equipment to record the three-phase voltage u output by the photovoltaic system a (t),u b (t),u c (t) and three-phase current i a (t), ib (t), i c (t), and simultaneously record the power quality parameters on the grid side;
[0194] Let the sampling frequency be f s , the sampling time is T, then the discrete form of the collected signal is expressed as:
[0195] u n [k]=u n (kT s ), i n [k]=i n (kT s ), k=0,1,…,N-1,
[0196] Among them, u n [k] is the discrete value of the voltage of the nth phase at the kth sampling point, i n [k] is the discrete value of the n-th phase current at the k-th sampling point, T s =1 / f s is the sampling period, N=Tf s is the total number of sampling points;
[0197] Step 1.2, signal denoising:
[0198] Wavelet transform is used to decompose the signal at multiple scales to remove high-frequency noise components;
[0199] Assuming the number of wavelet decomposition layers is J and the mother wavelet function is ψ(t), the wavelet transform coefficient of the signal is expressed as:
[0200]
[0201] Among them, W ψ (a, b) is the wavelet transform coefficient of signal u(t) at scale a and position b, ψ(t) is the mother wavelet function, ψ * (t) is the conjugate complex number of the mother wavelet function, a is the scale factor, and b is the translation factor;
[0202] The denoising process consists of the following steps:
[0203] Calculate the wavelet transform coefficient W of the signal ψ (a,b);
[0204] Set threshold λ j , perform soft threshold filtering on high-frequency wavelet coefficients:
[0205]
[0206] Among them, λ j is the threshold of the j-th layer wavelet decomposition, is the update coefficient after soft threshold processing;
[0207] Perform wavelet reconstruction to obtain the denoised signal
[0208] Step 1.3, signal normalization:
[0209] Using the min-max normalization method, the signal is mapped to the interval [1, 2]:
[0210]
[0211] Among them, u(t) is the original voltage signal, i(t) is the original current signal, u min is the minimum value of the voltage signal, u max is the maximum value of the voltage signal, i min is the minimum value of the current signal, i max is the maximum value of the current signal, is the normalized voltage signal, is the normalized current signal;
[0212] Step 1.4, signal time alignment:
[0213] The time delay between signals is calculated by the cross-correlation function. The calculation formula of the cross-correlation function is:
[0214]
[0215] Among them, R ui (τ) is the voltage signal u a [k] and current signal i a [k] The cross-correlation function value between u a [k] is the kth discrete value of the ath phase voltage signal, i a [k] is the kth discrete value of the ath phase current signal, and τ is the time delay;
[0216] Find the τ corresponding to the maximum mutual correlation coefficient max , and make corresponding time offset adjustments to the signals so that all signals are aligned to the same time base;
[0217] Step S2, adaptively decomposing the disturbance signal based on improved variational modal decomposition, automatically selecting the number of modes using information geometry optimization, and selecting the penalty factor through Bayesian optimization;
[0218] Step 2.1, mathematical modeling of variational mode decomposition:
[0219] Assume that the preprocessed disturbance signal is expressed as x(t):
[0220]
[0221] Among them, u k (t) is the kth group of modal components, x(t) is the preprocessed disturbance signal, and K is the target mode number;
[0222] Each modal component u k (t) at a specific center frequency ω k The neighborhood has limited bandwidth;
[0223] Variational mode decomposition achieves signal decomposition by solving the following variational optimization problem:
[0224]
[0225] Among them, H(u k ) is the modal component u k (t) is the signal after Hilbert transform, j is the imaginary unit, ω k is the center frequency of the modal component, is the time derivative operator, u i (t) is the modal component of group i;
[0226] The constraints are introduced using the Lagrange multiplier method to construct the Lagrange function:
[0227]
[0228]
[0229] Among them, L({u k},{ω k},λ) is the Lagrangian function of variational mode decomposition, α is the penalty factor, λ is the Lagrangian multiplier, and <·,·> is the inner product operation;
[0230] Step 2.2, information geometry optimization automatically selects the number of modes:
[0231] Based on the information geometry optimization method, the spectral entropy H of each modal component is maximized. k Calculate the optimal number of modes:
[0232]
[0233] in, is the modal component u k (t) at frequency f n Normalized energy distribution at :
[0234]
[0235] Let the total information entropy H total As the objective function:
[0236]
[0237] Traverse different modal numbers K and find the total The maximum optimal modal number K * :
[0238]
[0239] Among them, H k is the spectral entropy of the kth modal component, is the modal component u k (t) at frequency f n Normalized energy distribution at , U k (f n ) is the modal component u k (t) at frequency f n Fourier transform at, N is the number of frequency points, H total is information entropy, K * is the optimal mode number;
[0240] Step 2.3, Bayesian optimization to determine the penalty factor: Use the Bayesian optimization method to decompose the reconstruction error E of the signal rec As the objective function:
[0241]
[0242] Gaussian process regression model is used to predict E under different α values rec , and select the optimal penalty factor based on the expected improvement criterion:
[0243]
[0244] Among them, E rec is the signal reconstruction error, α is the penalty factor, is the reconstructed signal of all modal components after variational mode decomposition, To calculate the penalty factor α that minimizes the reconstruction error * ;
[0245] Step 2.4, calculate the decomposed disturbance signal:
[0246] Based on the determined K * and α * , perform the final variational mode decomposition on the preprocessed disturbance signal to obtain K * Group modal components:
[0247]
[0248] Among them, each u k(t) represents the different frequency components in the photovoltaic disturbance signal;
[0249] At this point, the adaptive decomposition of the disturbance signal based on the improved variational mode decomposition is completed, and the decomposed disturbance signal u is obtained k (t),u k (t) will be used to construct the Hamiltonian dynamics model in step S3;
[0250] Step S3, performing disturbance characteristic modeling based on the Hamiltonian dynamics model, constructing a Hamiltonian system and calculating the phase space dynamic equation of the disturbance signal;
[0251] Step 3.1, Hamiltonian dynamics modeling:
[0252] Assume the phase space of the disturbance signal is (q k (t),p k (t)), where q k (t) and p k (t) represents the modal component u k (t) position and momentum;
[0253]
[0254] Among them, m k is the modal component u k (t) mass, V(q k ) is the potential energy function;
[0255] For the disturbance signal, the momentum p k (t) and position q k (t) satisfies the following Hamiltonian equation:
[0256]
[0257] in, and Indicates speed and force;
[0258] Step 3.2, construct the phase space model of the disturbance signal:
[0259] By calculating the modal components u k (t) is numerically integrated to obtain the motion trajectory in the phase space;
[0260] Assume that each modal component u k The potential energy function V(q k ) is a simple quadratic potential energy:
[0261]
[0262] Among them, k k is the modal component uk (t) elastic constant;
[0263] By solving the Hamiltonian equation, we obtain the dynamic equation in phase space:
[0264]
[0265] Among them, q k (t) and p k (t) represents the modal component q k (t),p k (t) position and momentum in phase space;
[0266] Step 3.3, calculate the phase space dynamic equation of the disturbance signal:
[0267] The classic fourth-order Runge-Kutta algorithm is used to numerically integrate the Hamiltonian equation, setting the time step Δt and iteratively updating q k (t) and p k The value of (t):
[0268]
[0269] p k (t+Δt)=p k (t)-Δt·k k q k (t),
[0270] Step 3.4: Combine the disturbance signal to perform Hamiltonian dynamics modeling analysis:
[0271] After step S3 is completed, the dynamic characteristics of the disturbance signal have been modeled by the Hamiltonian dynamics model, and then step S4 is entered to further analyze the nonlinear characteristics of the disturbance signal by Lyapunov exponents;
[0272] Step S4, calculating the Lyapunov exponent of the Hamiltonian dynamics equation established in step S3, performing nonlinear characteristic analysis of the disturbance signal based on the Lyapunov exponent, calculating the maximum Lyapunov exponent of the disturbance signal and determining its chaotic characteristics;
[0273] Step 4.1, mathematical definition of Lyapunov exponent:
[0274] In the phase space, the phase trajectory of the disturbance signal is described by the Hamiltonian dynamics equation, and the state vector is: X(t) = [q1(t), p1(t), q2(t), p2(t), ..., q K* (t),p K* (t)] T ,
[0275] Among them, qK (t) and p K (t) is the modal component u k (t) is the phase space coordinate, K * is the optimal number of modal components;
[0276] The time evolution of the disturbance signal in phase space is described by the following Hamiltonian dynamics equation:
[0277]
[0278] Where, F(X(t)) is the evolution equation of the system;
[0279] At the initial time t0, adjacent trajectory points X(t0) and X′(t0) in the phase space are selected, and the initial small perturbation is: δX(t0)=X′(t0)-X(t0),
[0280] The evolution after time t satisfies:
[0281] Where λ is the Lyapunov exponent, which describes the exponential growth rate of the disturbance between adjacent trajectories. The rate of change of the distance between trajectories is given by the following formula:
[0282] Step 4.2, calculate the maximum Lyapunov exponent:
[0283] To calculate the maximum Lyapunov exponent, we need to perform numerical simulation on the phase space trajectory and use the small perturbation method to track the divergence of adjacent trajectories. The specific calculation steps are as follows:
[0284] Step 4.2.1, in the phase space trajectory X(t), select the initial point X(t0), and select the initial small perturbation δX(t0) in the neighborhood to satisfy: ||δX(t0)||=∈0,
[0285] Among them, ∈0 is the initial disturbance amplitude;
[0286] Step 4.2.2, solve the Hamiltonian dynamics equation by numerical integration method and calculate the distance between the perturbed trajectory X′(t) and the original trajectory X(t): d(t) = ||X′(t) - X(t)||,
[0287] Among them, d(t) reflects the changing trend of the disturbance trajectory over time.
[0288] Step 4.2.3: In the logarithmic coordinate system, plot the curve of ln(t) changing with time t, and use the linear regression method to solve the slope λ max :
[0289] Among them, N t is the number of time steps;
[0290] Step 4.3: Determine the chaotic characteristics based on the maximum Lyapunov exponent:
[0291] According to the maximum Lyapunov exponent λ max The numerical value of is used to judge the chaotic characteristics of the disturbance signal:
[0292] If λ max >0, the disturbance signal has chaotic characteristics;
[0293] If λ max =0, the system belongs to quasi-periodic motion;
[0294] If λ max <0, then the system has an attractor;
[0295] Step 4.4, calculate multiple groups of Lyapunov exponents:
[0296] Except for the maximum Lyapunov exponent λ max In addition, calculate other Lyapunov exponents λ i To further analyze the stability of the disturbance signal, the complete Lyapunov exponent spectrum Λ consists of multiple groups of exponents:
[0297] Λ=[λ1,λ2,…,λ 2K* ],
[0298] Where Λ is the optimal mode number;
[0299] Using the QR decomposition method, calculate the Lyapunov index spectrum and first construct the Jacobian matrix J(X) of the perturbation equation: Then, the exponential divergence rate of the orthogonal basis is tracked by QR decomposition to calculate all Lyapunov exponents;
[0300] Step 4.5, result output and next step prediction: After completing the Lyapunov exponent calculation, the nonlinear characteristics of the disturbance signal are obtained:
[0301] If λ max >0, it indicates that the disturbance signal has chaotic characteristics, and a prediction model that adapts to chaotic characteristics is used.
[0302] If λ max ≤0, the disturbance signal exhibits periodic and convergent behavior, and the traditional time series prediction model is used;
[0303] At this point, step S4 is completed, the Lyapunov exponent of the disturbance signal is successfully extracted, and the nonlinear characteristics are analyzed. Next, step S5 is entered to perform short-term and long-term predictions of the disturbance signal based on meteorological data;
[0304] Step S5: Combined with the disturbance prediction of meteorological data, the ARIMA model is used for short-term disturbance prediction, and the Transformer-LSTM hybrid model is used for long-term disturbance prediction;
[0305] Step 5.1, short-term disturbance forecast based on ARIMA model:
[0306] Assume that the time series of the disturbance signal is: s(t)=[s1,s2,...,s N ], where s N is the disturbance signal value of N time steps, N is the total number of observation time points;
[0307] Step 5.1.1, establishment of ARIMA model:
[0308] The ARIMA model is determined by the parameters (p, d, q):
[0309] p is the autoregressive order, which indicates the influence of the value at the previous p moments on the current value;
[0310] d is the difference order, which indicates the number of differences required to stabilize the time series;
[0311] q is the moving average order, which represents the influence of the first q error terms on the current value;
[0312] The mathematical expression of the ARIMA model is:
[0313] Φ p (B)(1-B) d s(t)=Θ q (B)∈(t),
[0314] Among them, B is the backshift operator, Φ p (B) is the polynomial of the autoregressive part, Θ q (B) is the polynomial of the moving average part, ∈(t) is the white noise term;
[0315] Step 5.1.2, use AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) to optimize and select the optimal (p, d, q):
[0316] AIC=-2lnL+2k, BIC=-2lnL+klnN,
[0317] Where L is the log-likelihood function of the model, k = p + q + d is the number of model parameters, and N is the total number of observation time points;
[0318] The optimal parameter (p * ,d * ,q * ) is determined by the following optimization objectives:
[0319]
[0320] Step 5.1.3, predict disturbance signals: Based on the optimal ARIMA model, predict short-term disturbance signals:
[0321]
[0322] Among them, h is the prediction step size, is the predicted disturbance signal value at the future time step t+h, φ i is the autoregressive coefficient, θ j is the moving average coefficient;
[0323] Step 5.2, long-term disturbance prediction is based on the Transformer-LSTM hybrid model:
[0324] Due to the chaotic characteristics of the disturbance signal, a deep learning model that can capture long-term dependencies is adopted, namely the Transformer-LSTM hybrid model;
[0325] Step 5.2.1, LSTM network modeling. The LSTM structure consists of an input gate, a forget gate, and an output gate. The calculation formula is as follows:
[0326] Forget gate: f t =σ(W f ·[h t-1 ,x t ]+b f )
[0327] Among them, x t is the current input, h t-1 is the hidden state at the previous moment, W f and b f are the weight and bias of the forget gate, σ(·) is the Sigmoid function;
[0328] Input gate: i t =σ(W i ·[h t-1 ,x t ]+b i ),
[0329]
[0330] Among them, i t is the output of the input gate, W i is the weight matrix of the input gate, b i is the bias vector of the input gate, is the candidate cell state, W C is the weight matrix of the candidate cell state, b Cis the bias vector of the candidate cell state;
[0331] Cell status update:
[0332] Among them, C t is the current cell state, C t-1 is the cell state at the previous moment;
[0333] Output gate: o t =σ(W o ·[h t-1 ,x t ]+b o ),
[0334] h t =o t tanh(C t ),
[0335] Among them, t is the output of the output gate, W o is the weight matrix of the output gate, b o is the bias vector of the output gate, h t is the final output of LSTM;
[0336] Step 5.2.2, Transformer processes global features: Transformer uses the self-attention mechanism to extract global features. Given the input perturbation signal sequence S, define the query Q, key K and value V:
[0337] Q=W Q S, K = W K S, V = W V S,
[0338] Among them, W Q 、W K 、W K is the linear transformation matrix;
[0339] Calculate attention weights:
[0340] Among them, d k is the dimension of the key vector;
[0341] Step 5.2.3, predict the disturbance signal: use Transformer to extract long-term dependency features and use LSTM to predict:
[0342]
[0343] Among them, f Transformer is the feature extracted by Transformer, f LSTMis the LSTM prediction function;
[0344] Step 5.3, perform disturbance correction based on meteorological data:
[0345] The disturbance signal of the photovoltaic system is affected by meteorological factors. Meteorological data include:
[0346] Solar radiation intensity G(t), temperature T(t), wind speed T(t) and cloud cover T(t);
[0347] Construct a disturbance correction model:
[0348]
[0349] Among them, β1, β2, β3, and β4 are the weights of meteorological factors. is the forecast value after the meteorological data is corrected, is the original predicted value;
[0350] Step 5.4, result output and application:
[0351] Short-term forecasting is used for real-time monitoring of photovoltaic systems;
[0352] Long-term forecasts are used for power grid dispatching;
[0353] The prediction results will be used for real-time warning in step S6 to prevent the impact of photovoltaic disturbances on the power grid. At this point, step S5 is completed, and the short-term and long-term evolution trends of the disturbance signal are successfully predicted.
[0354] Step S6: Output and apply the results, output the disturbance identification results and issue a real-time warning.
[0355] This embodiment adopts a multi-stage processing process to deal with the low-frequency harmonics and short-term fluctuations commonly seen in photovoltaic systems. The technical solution of the present invention starts with data acquisition and preprocessing, performs adaptive signal decomposition through improved variational mode decomposition, then uses Hamiltonian dynamics and Lyapunov exponents to extract nonlinear features, and finally combines meteorological data to build a hybrid prediction model.
[0356] Step S1, data collection and preprocessing:
[0357] Data acquisition uses a dedicated data acquisition device to simultaneously record the three-phase voltage, current and grid-side parameters output by the photovoltaic system; the sampling frequency and sampling duration are set to form a discrete signal sequence; signal denoising uses wavelet transform to perform multi-scale decomposition to eliminate high-frequency noise;
[0358] In the signal denoising process, the number of decomposition layers J and the mother wavelet function are set, the coefficients of each layer are soft-threshold filtered, and then the signal is restored through wavelet reconstruction;
[0359] Signal normalization and time alignment use the minimum-maximum normalization method to map the original data to a predetermined interval, then calculate the time delay by calculating the cross-correlation function, align the signals of each channel to ensure synchronization.
[0360] Step S2, adaptive decomposition based on improved variational mode decomposition:
[0361] An improved variational modal decomposition model is established to establish a mathematical model for the preprocessed disturbance signal, and each modal component in the model is set to meet the finite bandwidth characteristics; the automatic parameter selection uses the information geometry optimization method to automatically select the optimal mode number K based on the spectral entropy maximization; through Bayesian optimization, the reconstruction error is used as the objective function to dynamically determine the penalty factor α; the decomposition calculation decomposes the signal into several modes, and each mode reflects disturbances in different frequency bands; in the numerical calculation process, the iterative method is used to update the Lagrange multiplier to ensure that the constraints are met; after using the improved variational modal decomposition, there is no need to manually fix the parameters, so that the decomposition results fit the actual changes of the signal, thereby improving the accuracy of disturbance pattern recognition.
[0362] Step S3, Hamiltonian dynamics model construction:
[0363] Phase space construction constructs position and momentum variables for each modal component, establishes the corresponding Hamiltonian function, and simply uses quadratic potential energy to describe the vibration characteristics; numerical integration and dynamic equations use the fourth-order Runge-Kutta algorithm to numerically integrate the Hamiltonian equation to obtain the motion trajectory of each mode in the phase space, truly reflecting the disturbance energy transfer, and through the Hamiltonian model, analyzes the disturbance signal from the perspective of energy conservation, intuitively revealing the inherent physical mechanism of the signal.
[0364] Step S4, Lyapunov exponent calculation and chaos feature extraction:
[0365] The initial perturbation setting selects adjacent points in the phase space as the initial perturbation to ensure that the perturbation is small; after the exponential calculation is numerically integrated, the divergence rate of adjacent trajectories is tracked, the maximum Lyapunov exponent is calculated, and the Lyapunov exponent spectrum is simultaneously obtained using the QR decomposition method; characteristic judgment: if the maximum Lyapunov exponent is greater than zero, it proves that the perturbation signal has chaotic characteristics; a deep interpretation of the nonlinear behavior of the perturbation signal is achieved, improving the physical interpretability of the signal and the monitoring reliability.
[0366] Step S5, disturbance prediction based on meteorological data:
[0367] For short-term forecasting, the ARIMA model is used to model the time series characteristics of the disturbance signal, and the AIC and BIC criteria are used to determine the optimal model parameters to achieve short-term trend forecasting;
[0368] For long-term prediction, a Transformer-LSTM hybrid model is constructed. Transformer is used to extract global features, while LSTM captures long-term dependencies, effectively grasping the long-term evolution of disturbances.
[0369] Meteorological data fusion synchronously collects meteorological data, builds a correction model, corrects the prediction results, and takes into account the influence of the external environment. At the same time, by combining meteorological data, the prediction model is close to reality and the accuracy of the warning is significantly improved.
[0370] Step S6, result output and application:
[0371] The disturbance identification result integrates the decomposition results, dynamic model data and predicted values output from each step to form a real-time disturbance warning signal;
[0372] The application scenario can be used for real-time monitoring of photovoltaic systems; the power grid dispatching center can respond quickly based on early warning signals to ensure power grid stability.
[0373] A terminal device includes a memory, a processor, and a computer program stored in the memory and running on the processor, wherein the computer program is configured to execute a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply.
[0374] A storage medium stores a computer program. When the computer program is executed by a processor, a method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply is implemented.
[0375] The terminal equipment automatically collects multiple power quality signals and performs signal preprocessing, denoising, normalization, and time alignment. It also uses improved variational mode decomposition to adaptively decompose disturbance signals, constructs a Hamiltonian dynamics model, calculates Lyapunov exponents, and extracts signal nonlinear characteristics. It also combines meteorological data for short-term and long-term disturbance forecasting.
[0376] The terminal equipment has a simple structure and is easy to integrate. It can process and feedback disturbance information in real time on site, reduce the misjudgment rate, and ensure monitoring accuracy.
[0377] The storage medium can be non-volatile solid-state memory and flash memory; the computer program includes a complete signal acquisition, preprocessing, decomposition, modeling and prediction process; it can be easily ported to different hardware platforms to ensure technology reproduction and system upgrades.
[0378] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.
Claims
1. A method for identifying power quality disturbance characteristics of distributed photovoltaic interference power supply, characterized in that: include: Step S1, data acquisition and preprocessing, collecting the power quality signal of the photovoltaic system and performing denoising, normalization and time alignment; Step S2, adaptively decomposing the disturbance signal based on improved variational modal decomposition, automatically selecting the number of modes using information geometry optimization, and selecting the penalty factor through Bayesian optimization; Step S3, performing disturbance characteristic modeling based on the Hamiltonian dynamics model, constructing a Hamiltonian system and calculating the phase space dynamic equation of the disturbance signal; Step S4, calculating the Lyapunov exponent of the Hamiltonian dynamics equation established in step S3, performing nonlinear characteristic analysis of the disturbance signal based on the Lyapunov exponent, calculating the maximum Lyapunov exponent of the disturbance signal and determining its chaotic characteristics; Step S5: Combined with the disturbance prediction of meteorological data, the ARIMA model is used for short-term disturbance prediction, and the Transformer-LSTM hybrid model is used for long-term disturbance prediction; Step S6: Output and apply the results, output the disturbance identification results and issue a real-time warning.
2. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S1, data collection and preprocessing further include: Step 1.1, collect power quality signals: When collecting power quality signals, use data acquisition equipment to record the three-phase voltage u output by the photovoltaic system a (t),u b (t),u c (t) and three-phase current i a (t), i b (t), i c (t), and simultaneously record the power quality parameters on the grid side; Let the sampling frequency be f s , the sampling time is T, then the discrete form of the collected signal is expressed as: you n [k]=u n (kT s ),i n [k]=i n (kT s ), k=0,1,…,N-1, Among them, u n [k] is the discrete value of the voltage of the nth phase at the kth sampling point, i n [k] is the discrete value of the n-th phase current at the k-th sampling point, T s =1 / f s is the sampling period, N=Tf s is the total number of sampling points; Step 1.2, signal denoising: Wavelet transform is used to decompose the signal at multiple scales to remove high-frequency noise components; Assuming the number of wavelet decomposition layers is J and the mother wavelet function is ψ(t), the wavelet transform coefficient of the signal is expressed as: Among them, W ψ (a, b) is the wavelet transform coefficient of signal u(t) at scale a and position b, ψ(t) is the mother wavelet function, ψ * (t) is the conjugate complex number of the mother wavelet function, a is the scale factor, and b is the translation factor; The denoising process consists of the following steps: Calculate the wavelet transform coefficient W of the signal ψ (a,b); Set threshold λ j , perform soft threshold filtering on high-frequency wavelet coefficients: Among them, λ j is the threshold of the j-th layer wavelet decomposition, is the update coefficient after soft threshold processing; Perform wavelet reconstruction to obtain the denoised signal Step 1.3, signal normalization: Using the min-max normalization method, the signal is mapped to the interval [1, 2]: Among them, u(t) is the original voltage signal, i(t) is the original current signal, u min is the minimum value of the voltage signal, u max is the maximum value of the voltage signal, i min is the minimum value of the current signal, i max is the maximum value of the current signal, is the normalized voltage signal, is the normalized current signal; Step 1.4, signal time alignment: The time delay between signals is calculated by the cross-correlation function. The calculation formula of the cross-correlation function is: Among them, R ui (τ) is the voltage signal u a [k] and current signal i a [k] The cross-correlation function value between u a [k] is the kth discrete value of the ath phase voltage signal, i a [k] is the kth discrete value of the ath phase current signal, and τ is the time delay; Find the τ corresponding to the maximum mutual correlation coefficient max , and make corresponding time offset adjustments to the signals so that all signals are aligned to the same time base.
3. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S2, the adaptive decomposition of the disturbance signal based on the improved variational modal decomposition further includes: Step 2.1, mathematical modeling of variational mode decomposition: Assume that the preprocessed disturbance signal is expressed as x(t): Among them, u k (t) is the kth group of modal components, x(t) is the preprocessed disturbance signal, and K is the target mode number; Each modal component u k (t) at a specific center frequency ω k The neighborhood has limited bandwidth; Variational mode decomposition achieves signal decomposition by solving the following variational optimization problem: Among them, H(u k ) is the modal component u k (t) is the signal after Hilbert transform, j is the imaginary unit, ω k is the center frequency of the modal component, is the time derivative operator, u i (t) is the modal component of group i; The constraints are introduced using the Lagrange multiplier method to construct the Lagrange function: Among them, L({u k },{ω k },λ) is the Lagrangian function of variational mode decomposition, α is the penalty factor, λ is the Lagrangian multiplier, and <·,·> is the inner product operation; Step 2.2, information geometry optimization automatically selects the number of modes: Based on the information geometry optimization method, the spectral entropy H of each modal component is maximized. k Calculate the optimal number of modes: in, is the modal component u k (t) at frequency f n Normalized energy distribution at : Let the total information entropy H total As the objective function: Traverse different modal numbers K and find the H total The maximum optimal modal number K * : Among them, H k is the spectral entropy of the kth modal component, is the modal component u k (t) at frequency f n Normalized energy distribution at , U k (f n ) is the modal component u k (t) at frequency f n Fourier transform at, N is the number of frequency points, H total is information entropy, K * is the optimal modal number; Step 2.3, Bayesian optimization to determine the penalty factor: Use the Bayesian optimization method to decompose the reconstruction error E of the signal rec As the objective function: Gaussian process regression model is used to predict E under different α values rec , and select the optimal penalty factor based on the expected improvement criterion: Among them, E rec is the signal reconstruction error, α is the penalty factor, is the reconstructed signal of all modal components after variational mode decomposition, To calculate the penalty factor α that minimizes the reconstruction error * ; Step 2.4, calculate the decomposed disturbance signal: Based on the determined K * and α * , perform the final variational mode decomposition on the preprocessed disturbance signal to obtain K * Group modal components: Among them, each u k (t) represents the different frequency components in the photovoltaic disturbance signal; At this point, the adaptive decomposition of the disturbance signal based on the improved variational mode decomposition is completed, and the decomposed disturbance signal u is obtained k (t),u k (t) will be used to construct the Hamiltonian dynamics model in step S3.
4. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S3, the disturbance characteristic modeling based on the Hamiltonian dynamics model further includes: After completing the adaptive decomposition of the disturbance signal in step S2, we get K * Group modal component u k (t), the modal components represent the disturbance signals of different frequency bands in the photovoltaic system and will be used to construct the Hamiltonian dynamics model in step S3; Step 3.1, Hamiltonian dynamics modeling: Assume the phase space of the disturbance signal is (q k (t),p k (t)), where q k (t) and p k (t) represents the modal component u k (t) position and momentum; Among them, m k is the modal component u k (t) mass, V(q k ) is the potential energy function; For the disturbance signal, the momentum p k (t) and position q k (t) satisfies the following Hamiltonian equation: in, and Indicates speed and force; Step 3.2, construct the phase space model of the disturbance signal: By calculating the modal components u k (t) is numerically integrated to obtain the motion trajectory in the phase space; Assume that each modal component u k The potential energy function V(q k ) is a simple quadratic potential energy: Among them, k k is the modal component u k (t) elastic constant; By solving the Hamiltonian equation, we obtain the dynamic equation in phase space: Among them, q k (t) and p k (t) represents the modal component q k (t),p k (t) position and momentum in phase space; Step 3.3, calculate the phase space dynamic equation of the disturbance signal: The classic fourth-order Runge-Kutta algorithm is used to numerically integrate the Hamiltonian equation, setting the time step Δt and iteratively updating q k (t) and p k The value of (t): p k (t+Δt)=p k (t)-Δt·k k q k (t), Step 3.4: Combine the disturbance signal to perform Hamiltonian dynamics modeling analysis: After step S3 is completed, the dynamic characteristics of the disturbance signal have been modeled by the Hamiltonian dynamics model, and then step S4 is entered to further analyze the nonlinear characteristics of the disturbance signal by using the Lyapunov exponent.
5. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S4, the nonlinear characteristic analysis of the disturbance signal based on the Lyapunov exponent further includes: In step S3, a Hamiltonian dynamics model of the disturbance signal is constructed, and the phase space dynamic equation of the disturbance signal is obtained. To further analyze the nonlinear characteristics and chaotic characteristics of the disturbance signal, in step S4, the Lyapunov exponent needs to be calculated, and the chaotic characteristics of the disturbance signal are judged based on the maximum Lyapunov exponent. Step 4.1, mathematical definition of Lyapunov exponent: In the phase space, the phase trajectory of the disturbance signal is described by the Hamiltonian dynamics equation, and the state vector is: X(t) = [q1(t), p1(t), q2(t), p2(t), ..., q K* (t),p K* (t)] T , Among them, q K (t) and p K (t) is the modal component u k (t) is the phase space coordinate, K * is the optimal number of modal components; The time evolution of the disturbance signal in phase space is described by the following Hamiltonian dynamics equation: Where, F(X(t)) is the evolution equation of the system; At the initial time t0, select adjacent trajectory points X(t0) and X ′ (t0), the initial small perturbation is: δX(t0)=X ′ (t0)-X(t0), The evolution after time t satisfies: Where λ is the Lyapunov exponent, which describes the exponential growth rate of the disturbance between adjacent trajectories. The rate of change of the distance between trajectories is given by the following formula: Step 4.2, calculate the maximum Lyapunov exponent: To calculate the maximum Lyapunov exponent, we need to perform numerical simulation on the phase space trajectory and use the small perturbation method to track the divergence of adjacent trajectories. The specific calculation steps are as follows: Step 4.2.1, in the phase space trajectory X(t), select the initial point X(t0), and select the initial small perturbation δX(t0) in the neighborhood to satisfy: ||δX(t0)||=∈0, Among them, ∈0 is the initial disturbance amplitude; Step 4.2.2, solve the Hamiltonian dynamics equation by numerical integration method and calculate the distance between the perturbed trajectory X′(t) and the original trajectory X(t): d(t) = ||X′(t) - X(t)||, Among them, d(t) reflects the changing trend of the disturbance trajectory over time; Step 4.2.3: In the logarithmic coordinate system, plot the curve of ln(t) changing with time t, and use the linear regression method to solve the slope λ max : Among them, N t is the number of time steps; Step 4.3: Determine the chaotic characteristics based on the maximum Lyapunov exponent: According to the maximum Lyapunov exponent λ max The numerical value of is used to judge the chaotic characteristics of the disturbance signal: If λ max >0, the disturbance signal has chaotic characteristics; If λ max =0, the system belongs to quasi-periodic motion; If λ max <0, then the system has an attractor; Step 4.4, calculate multiple groups of Lyapunov exponents: Except for the maximum Lyapunov exponent λ max In addition, calculate other Lyapunov exponents λ i To further analyze the stability of the disturbance signal, the complete Lyapunov exponent spectrum Λ consists of multiple groups of exponents: Λ=[λ1,λ2,…,λ 2K* ], Where Λ is the optimal mode number; Using the QR decomposition method, calculate the Lyapunov index spectrum and first construct the Jacobian matrix J(X) of the perturbation equation: Then, the exponential divergence rate of the orthogonal basis is tracked by QR decomposition to calculate all Lyapunov exponents; Step 4.5, result output and next step prediction: After completing the Lyapunov exponent calculation, the nonlinear characteristics of the disturbance signal are obtained: If λ max >0, it indicates that the disturbance signal has chaotic characteristics, and a prediction model that adapts to chaotic characteristics is used; If λ max ≤0, the disturbance signal exhibits periodic and convergent behavior, and the traditional time series prediction model is used; At this point, step S4 is completed, the Lyapunov exponent of the disturbance signal is successfully extracted, and the nonlinear characteristics are analyzed. Next, step S5 is entered to perform short-term and long-term predictions on the disturbance signal based on meteorological data.
6. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S5, the disturbance prediction based on meteorological data further includes: In step S4, the nonlinear characteristics of the disturbance signal are analyzed by calculating the Lyapunov exponent and the chaotic behavior is determined. In step S5, the disturbance signal is predicted in the short and long term in combination with the meteorological data of the photovoltaic system. Step 5.1, short-term disturbance forecast based on ARIMA model: Assume that the time series of the disturbance signal is: s(t)=[s1,s2,...,s N ], where s N is the disturbance signal value of N time steps, where N is the total number of observation time points; Step 5.1.1, establishment of ARIMA model: The ARIMA model is determined by the parameters (p, d, q): p is the autoregressive order, which indicates the influence of the value at the previous p moments on the current value; d is the difference order, which indicates the number of differences required to stabilize the time series; q is the moving average order, which represents the impact of the first q error terms on the current value; The mathematical expression of the ARIMA model is: F p (B)(1-B) d s(t)=Θ q (B)∈(t), Among them, B is the backshift operator, Φ p (B) is the polynomial of the autoregressive part, Θ q (B) is the polynomial of the moving average part, ∈(t) is the white noise term; Step 5.1.2, use AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) to optimize and select the optimal (p, d, q): AIC=-2lnL+2k, BIC=-2lnL+klnN, Where L is the log-likelihood function of the model, k = p + q + d is the number of model parameters, and N is the total number of observation time points; The optimal parameter (p * ,d * ,q * ) is determined by the following optimization objectives: Step 5.1.3, predict disturbance signals: Based on the optimal ARIMA model, predict short-term disturbance signals: Among them, h is the prediction step size, is the predicted disturbance signal value at the future time step t+h, φ i is the autoregressive coefficient, θ j is the moving average coefficient.
7. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S5, the disturbance prediction based on meteorological data further includes: Step 5.2, long-term disturbance prediction is based on the Transformer-LSTM hybrid model: Due to the chaotic characteristics of the disturbance signal, a deep learning model that can capture long-term dependencies is adopted, namely the Transformer-LSTM hybrid model; Step 5.2.1, LSTM network modeling. The LSTM structure consists of an input gate, a forget gate, and an output gate. The calculation formula is as follows: Forget gate: f t =σ(W f ·[h t-1 ,x t ]+b f ) Among them, x t is the current input, h t-1 is the hidden state of the previous moment, W f and b f are the weight and bias of the forget gate, σ(·) is the Sigmoid function; Input gate: i t =σ(W i ·[h t-1 ,x t ]+b i ), Among them, i t is the output of the input gate, W i is the weight matrix of the input gate, b i is the bias vector of the input gate, is the candidate cell state, W C is the weight matrix of the candidate cell state, b C is the bias vector of the candidate cell state; Cell status update: Among them, C t is the current cell state, C t-1 is the cell state at the previous moment; Output gate: o t =σ(W o ·[h t-1 ,x t ]+b o ), h t =o t ·tanh(C t ), Among them, t is the output of the output gate, W o is the weight matrix of the output gate, b o is the bias vector of the output gate, h t is the final output of LSTM; Step 5.2.2, Transformer processes global features: Transformer uses the self-attention mechanism to extract global features. Given the input perturbation signal sequence S, define the query Q, key K and value V: Q=W Q S,K=W K S,V=W V S, Among them, W Q 、W K 、W K is the linear transformation matrix; Calculate attention weights: Among them, d k is the dimension of the key vector; Step 5.2.3, predict the disturbance signal: use Transformer to extract long-term dependency features and use LSTM to predict: Among them, f Transformer is the feature extracted by Transformer, f LSTM is the LSTM prediction function.
8. The method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to claim 1, characterized in that: In step S5, the disturbance prediction based on meteorological data further includes: Step 5.3, perform disturbance correction based on meteorological data: The disturbance signal of the photovoltaic system is affected by meteorological factors. Meteorological data include: Solar radiation intensity G(t), temperature T(t), wind speed T(t) and cloud cover T(t); Construct a disturbance correction model: Among them, β1, β2, β3, and β4 are the weights of meteorological factors. is the forecast value after the meteorological data is corrected, is the original predicted value; Step 5.4, result output and application: Short-term forecasting is used for real-time monitoring of photovoltaic systems; Long-term forecasts are used for power grid dispatch; The prediction results will be used for real-time warning in step S6 to prevent the impact of photovoltaic disturbances on the power grid. At this point, step S5 is completed, and the short-term and long-term evolution trends of the disturbance signal are successfully predicted.
9. A terminal device, characterized in that: The invention comprises a memory, a processor and a computer program stored in the memory and running on the processor, wherein the computer program is configured to execute the method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to any one of claims 1 to 8.
10. A storage medium, characterized in that: The storage medium stores a computer program, and when the computer program is executed by a processor, the method for identifying power quality disturbance characteristics of a distributed photovoltaic interference power supply according to any one of claims 1 to 8 is implemented.
Citation Information
Cited By
Voltage monitoring method for small power generation equipment
CN121036358A
Rapid fault positioning method and system for photovoltaic power station
CN121308675A
System control method and device based on power system, storage medium and computer equipment
CN121529994A
Power quality composite disturbance collaborative modeling method based on multi-source data
CN121859602A