Adaptive integrated wind power prediction method based on probability density function shape control

Abstract

The invention provides an adaptive integrated wind power prediction method based on probability density function shape control, and the method is realized based on the following steps: collecting multi-dimensional fan time sequence data, and decomposing a multi-dimensional signal into K intrinsic mode components through multivariate variational mode decomposition; evaluating the complexity of each IMF by using a sample entropy (SE), merging sub-modules, eliminating noise interference and reducing dimensionality; training an adaptive fuzzy neural network based on the reconstructed features, wherein the parameter updating rule of the adaptive fuzzy neural network ensures global convergence through a Lyapunov function; a PDF shape controller is constructed, error distribution and a target form are matched in real time, and error distribution steady-state convergence is achieved through reference predictor weight resetting and PDF controller bandwidth self-adaptive updating adjustment. The method can be widely applied to new energy management core links such as power system dynamic safety evaluation, wind power plant cluster power dispatching optimization and wind-light-storage combined frequency modulation control, and is particularly suitable for prediction reliability enhancement under extreme meteorological conditions such as typhoon and thunderstorm.

Description

Technical Field

[0001] The present invention belongs to the technical field of wind power prediction, and in particular relates to an adaptive integrated wind power prediction method based on probability density function shape control. Background Art

[0002] With the rapid development of renewable energy, wind power, as a key component of clean energy, has seen its large-scale grid integration place higher demands on the stable operation of the power grid. The accuracy of wind power prediction (WPP) is directly related to the security and economic viability of the power grid. However, wind power data exhibits significant nonlinearity and nonstationarity, resulting in an asymmetric distribution of prediction errors. Traditional integrated models based on symmetric loss functions (such as mean squared error (MSE)) struggle to effectively capture the spatial distribution of these errors, limiting further improvements in prediction accuracy.

[0003] At present, wind power prediction methods are mainly divided into two categories: physical model driven methods and data driven methods. Data driven methods have become a research hotspot due to their strong adaptability, especially those based on deep learning models. However, existing deep learning models still have problems such as difficulty in parameter initialization and easy to fall into local optimality. Specifically, they are as follows: 1. The existing method directly weights the integrated sub-model prediction value, but does not consider the asymmetry of the error distribution, resulting in insufficient prediction ability of the model for high error intervals. 2. Although the improved variational mode decomposition (VMD) combined with the bidirectional GRU in the existing deep learning model can extract time series features to a certain extent, its parameter initialization relies on trial and error, and the optimization process is prone to fall into local optimality. It lacks theoretical support for global convergence, resulting in poor prediction stability. Summary of the Invention

[0004] The purpose of the present invention is to overcome the obvious deficiencies of existing wind power prediction technologies in processing asymmetric errors, global optimization of models and dynamic integration, and to provide an adaptive integrated wind power prediction method based on probability density function shape control.

[0005] To achieve the above objectives, the technical solutions provided by the present invention are:

[0006] The adaptive integrated wind power prediction method based on probability density function shape control includes the following steps:

[0007] Step 1: Obtain the original wind power data and the corresponding feature data to form a multivariate signal X(t); divide the multivariate signal X(t) into a base learning training set, a weight prediction set, and a test set, perform multivariate variational mode decomposition on each data set, obtain multiple intrinsic mode functions, calculate the sample entropy of each intrinsic mode function, merge the intrinsic mode functions whose sample entropy difference is less than the set threshold, and obtain multi-mode reconstruction features M1, M2, ..., M after merging. n , n is the number of patterns;

[0008] The characteristic data include wind speed, temperature and air pressure;

[0009] Step 2: Build a cluster of benchmark predictors:

[0010] Step 2.1: Reconstruct the multi-modal features M1, M2, ..., M n Each intrinsic mode function in the method establishes a basic fuzzy neural network as a reference predictor, wherein the reference predictor is used to output wind power according to input characteristic data;

[0011] Step 2.2: Use the base learning training set data obtained in step 1 to train each benchmark predictor. During the training process, dynamically update the weights of the benchmark predictor according to the prediction error; update the model parameters and perform gradient updates at each time step;

[0012] Step 2.3: Use the weighted prediction set data obtained in step 1 and make predictions through each benchmark predictor trained in step 2.2, and output the corresponding benchmark power prediction value;

[0013] Step 3: constructing a PDF shape controller, wherein the PDF shape controller is used to integrate the reference power prediction values ​​output by each reference predictor to obtain a final wind power prediction value;

[0014]

[0015] in, represents the benchmark power prediction value output by the kth benchmark predictor; w represents the weight matrix, which is composed of the weight of each benchmark predictor; the number of benchmark predictors is n, and the initial weight of each benchmark predictor is defined as 1 / n;

[0016] Step 4: Determine the weight matrix w k , we get the optimal benchmark predictor and the optimal PDF shape controller:

[0017] Input the benchmark power prediction value output by each benchmark predictor in step 2.3 into the constructed PDF shape controller to determine whether the change rate of the weight matrix at adjacent moments is less than the set change rate threshold. If so, the weight matrix w is obtained. kand the optimal benchmark predictor; otherwise, update the weight matrix and return to step 2.2 to update the weight of each benchmark predictor until the change rate of the weight matrix at adjacent moments meets the requirements;

[0018] The weight matrix update formula is:

[0019]

[0020] Where η is the learning rate, is the gradient of the objective function with respect to the weight matrix;

[0021] Objective function Among them, f pdf (e) is the actual integrated error probability density, f target (e) is the target integrated error probability density;

[0022] The probability density of the actual integrated error is Where: κ is the kernel bandwidth, representing the sample standard deviation; T is the total number of samples; L is the sample index, representing the Lth sample;

[0023] e is the ensemble error of the benchmark predictor; is the prediction error of the benchmark predictor at time t, P(t) is the actual value of wind power, is the baseline power prediction value;

[0024] Step 5: Input the test set data obtained in step 1 into the optimal benchmark predictor for prediction, and use the optimal PDF shape controller to output the final prediction result to verify the prediction error of the optimal benchmark predictor. If the prediction error meets the design requirements, execute step 6; otherwise, return to step 2 to adjust the model parameters of the benchmark predictor.

[0025] Step 6: Obtain actual wind farm characteristic data and use the optimal benchmark predictor and the optimal PDF shape controller to predict the output wind power.

[0026] Furthermore, the specific process of step 1 is as follows:

[0027] Step 1.1: Obtain raw wind power data and corresponding feature data to form a multivariate signal X(t), where the time resolution of the multivariate signal X(t) is less than 15 minutes; divide the multivariate signal X(t) into a base learning training set, a weight prediction set, and a test set;

[0028] Step 1.2: Establish the solution objective and constraints of minimizing the total bandwidth of the modal components in the multivariate mode decomposition, perform multivariate mode decomposition on each data set, and obtain K intrinsic mode functions;

[0029] The solution objective and constraint formula for minimizing the total bandwidth of the modal components are:

[0030]

[0031] Where X(t) represents the original input multivariate signal, u k (t) is the kth eigenmode function, ω k is the center frequency of the kth eigenmode function; Represents the derivative operation with respect to time t; the goal of multivariate variational mode decomposition is to decompose the original input multivariate signal into multiple modal signals;

[0032] Step 1.3: Use the alternating direction multiplier method to iteratively update the frequency domain modal components and frequency domain centers of each eigenmode function, and output the optimized eigenmode function;

[0033] Step 1.4: Calculate the sample entropy of each optimized intrinsic mode function, compare the size of each sample entropy, merge the optimized intrinsic mode functions with sample entropy difference less than 0.1 into the same mode, and take the optimized intrinsic mode functions with sample entropy difference not less than 0.1 as separate modes, and obtain multi-mode reconstruction features M1, M2, ..., M after merging. n , n is the number of patterns.

[0034] Furthermore, the process of step 1.3 is as follows:

[0035] The frequency domain modal component and center frequency of each intrinsic mode function are iteratively updated using the alternating direction multiplier method. When the frequency domain modal component change rate and center frequency change rate of each intrinsic mode function are less than 1%, the iterative update ends and the optimized intrinsic mode function is output.

[0036] The frequency domain modal component update formula is:

[0037] in, is the frequency domain of the kth eigenmode function in the l+1th iteration, is the synthetic signal of the original frequency domain signal, which serves as the reference input for modal update; is the frequency domain of the i-th intrinsic mode function in the l-th iteration; λ(ω) is the frequency domain form of the Lagrange multiplier, which corresponds to the weight term of the constraint condition and is the adjustment parameter; β * is the bandwidth penalty parameter, which is a positive real number and is used to control the distribution width of the mode in the frequency domain; is the center frequency of the kth eigenmode function after the lth iteration; Penalize the frequency domain components that deviate from the center frequency and force the modal energy to aggregate; ω is the independent variable of frequency domain analysis, which represents the frequency component of the signal after Fourier transform; represents the frequency domain superposition of all eigenmode functions except the kth eigenmode function in the lth iteration;

[0038] The center frequency update formula is:

[0039]

[0040] Where, represents the center frequency of the kth eigenmode function after the l+1th iteration update; represents the frequency domain obtained by the last iterative update, that is, the frequency domain of the kth intrinsic mode function in the l+1th iteration.

[0041] Furthermore, in step 1.4, the conditional probability estimate of the data subsequence matching in the time series of each optimized intrinsic mode function is calculated respectively, and the obtained conditional probability estimate is used as the sample entropy corresponding to each optimized intrinsic mode function;

[0042] The conditional probability estimate is calculated as follows:

[0043] Where N is the signal sequence length of the optimized intrinsic mode function, m is the embedding dimension, r is the tolerance threshold, which is 0.15, A is the logarithmic count of matching m+1-dimensional templates, and B is the logarithmic count of matching m-dimensional templates.

[0044] Furthermore, in step 2.1, the multi-mode reconstruction features M1, M2, ..., M n The process of establishing a baseline predictor for each eigenmode function in is:

[0045] Step 2.1.1: Construct a fuzzy neural network model and set the model parameters of the fuzzy neural network model; the fuzzy neural network includes an input layer, a membership function layer, a rule layer and an output layer;

[0046] The input layer is used to receive the historical window data x corresponding to the multi-modal reconstruction feature k (t)=[P k (t-τ),v k (t-τ),...,P k (t),v k (t)], where τ is the window length; P k (t) is the power component in the reconstructed feature of the kth mode, V k (t) is the feature component in the reconstructed feature set of the kth pattern;

[0047] The membership function layer adopts Gaussian membership function to convert clear input data into fuzzy sets and calculate the membership of input variables of different modes in the input data under different fuzzy rules;

[0048] The rule layer is used to calculate the trigger strength of each fuzzy rule and normalize all fuzzy rules to obtain the final trigger strength;

[0049] The output layer weights and fuses the output results of multiple fuzzy rules to obtain a wind power prediction result;

[0050] Step 2.1.2: Introduce auxiliary error and construct Lyapunov function; the Lyapunov function is used to judge the convergence degree of the prediction error during the training process of the fuzzy neural network model, so that the model parameters can be continuously updated until the prediction error converges to a minimum.

[0051] Furthermore, in step 2.2, the weight of the auxiliary error in each benchmark predictor is updated using the following rule:

[0052]

[0053] Among them, ΔΘ k (t) is the update amount of the auxiliary error weight of the k-th reference predictor; k (t) is the trigger strength vector of the k-th predictor, is the auxiliary error; η k (t) is the learning rate, e k (t) is the prediction error, is the current weight estimate of the k-th base predictor; is the final trigger strength vector.

[0054] Furthermore, in step 2.2, the model parameters are updated during the training process using the following formula, and gradient updates are performed at each time step:

[0055]

[0056] Where: represents the model parameters of the benchmark predictor at time step t, and all model parameters are updated after each iteration; L(θ t ) is the loss function; λ is the dynamic weight coefficient of the loss function; η is the learning rate; is the gradient of the loss function with respect to the error distribution parameter θ at time step t; m t is the bias-corrected first-order moment estimate, representing the exponentially weighted average of the gradient; To prevent the formula from having an extremely small constant with a denominator of zero, the value is taken as 0.01.

[0057] v tis the second-order moment estimate, which represents the exponentially weighted average of the square of the gradient; β1 is the decay rate of the first-order moment, ranging from (0,1), which is used to control the influence of momentum; is the bias correction term; β2 is the decay rate of the second-order moment, ranging from (0,1), which is used to control the influence of the square of the gradient.

[0058] The advantages of the present invention are:

[0059] 1. The method of the present invention synchronously decomposes wind speed and power data based on multivariate variational mode decomposition (MVMD) and sample entropy (SE), extracts multi-mode features from the original wind power sample data, merges sub-modes with similar complexity, eliminates noise interference and reduces data dimensionality, reduces the modeling difficulty of the benchmark predictor, and solves the frequency mismatch problem of single variable decomposition in traditional methods.

[0060] 2. In the present invention, an adaptive parameter update rule is designed based on Lyapunov theory to avoid falling into local optimum when updating model parameters during the benchmark predictor training process, thereby ensuring the global convergence effect.

[0061] 3. The present invention constructs a probability density function (PDF) shape controller, which matches the prediction error distribution and the target error shape in real time through the PDF shape controller, and realizes the steady-state convergence of the prediction error distribution by adaptively updating the weight of the baseline predictor and the weight matrix of the PDF shape controller, breaking through the limitation of the traditional use of mean square error (MSE) that only captures the second-order moment information, and improving the robustness of the baseline predictor.

[0062] Additional aspects and advantages of the present invention will be set forth in part in the description which follows and, in part, will be obvious from the description which follows, or may be learned by practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments with reference to the accompanying drawings, in which:

[0064] Figure 1 This is a flow chart of the adaptive integrated wind power prediction method based on probability density function shape control of the present invention;

[0065] Figure 2 It is a data preprocessing flow chart in the method of the present invention;

[0066] Figure 3 1 is a four-layer network structure diagram of the adaptive global convergence benchmark predictor constructed in the method of the present invention;

[0067] Figure 4 This is a flowchart of the integrated prediction output process of the wind power prediction method of the present invention;

[0068] Figure 5 Schematic diagram of the wind power sequence and the corresponding intrinsic mode function obtained by decomposition in the present invention;

[0069] Figure 6 is the membership function of the hub height wind speed characteristic in the embodiment of the present invention;

[0070] Figure 7 2. Schematic diagram of the distribution of fuzzy rule triggering strength corresponding to temperature characteristics in an embodiment of the present invention;

[0071] Figure 8 It is a global convergence discriminant diagram based on the Lyapunov function in an embodiment of the present invention;

[0072] Figure 9 is a comparison diagram of the original prediction error distribution and the prediction error distribution using the method of the present invention;

[0073] Figure 10 This is a comparison chart of the prediction results of the prediction method of the present invention, the traditional random forest method, and the LSTM prediction method with the actual wind power. DETAILED DESCRIPTION

[0074] The following describes in detail embodiments of the present invention. The embodiments are exemplary and intended to explain the present invention, but are not to be construed as limiting the present invention.

[0075] An embodiment of the present invention provides an adaptive integrated wind power prediction method based on shape control of a probability density function, comprising the following steps:

[0076] Step S1: Obtain the original wind power data and the corresponding feature data, and segment the multivariate signal X(t) composed of the original wind power and wind speed data into a base learner training set, a weight prediction set, and a test set according to a ratio of 0.6:0.2:0.2; perform multivariate variational mode decomposition on each data set to obtain multiple intrinsic mode functions, calculate the sample entropy of each intrinsic mode function, merge the intrinsic mode functions whose sample entropy difference is less than the set threshold, and obtain multi-mode reconstruction features M1, M2, ..., M after merging. n , n is the number of patterns. The specific process includes the following:

[0077] The multivariate variational mode decomposition (MVMD) process is:

[0078] The original data of wind power and wind speed, temperature, air pressure and other M-dimensional time series data are used to form a multivariate signal X(t), where M is the number of features. The multivariate signal time resolution is required to be less than fifteen minutes to ensure the capture of high-frequency dynamic features. Decompose the multivariate signal X(t) into k intrinsic mode functions (IMFs), and establish the solution objective and constraint formula for minimizing the total bandwidth of the modal components as follows:

[0079] Where: u k (t) is the kth eigenmode function, ω k is the center frequency of the kth eigenmode function, represents the derivative operation with respect to time t, X(t) represents the input multivariate signal, and the solution goal is to decompose the multivariate signal X(t) into multiple modal signals.

[0080] Iterative optimization of intrinsic mode functions: The frequency domain and center frequency of each intrinsic mode function are iteratively calculated using the alternating direction multiplier method during the multivariate variational mode decomposition (MVMD) process, and Lagrange multipliers are introduced in the calculation process to ensure convergence.

[0081] The frequency domain modal component update formula is:

[0082]

[0083] in, is the frequency domain representation of the kth mode function in the l+1th iteration, The frequency domain synthesis form of the original signal (which may be the frequency domain superposition of multi-channel signals in a multivariate scenario) is used as the reference input for modal update. is the frequency domain representation of the i-th mode function in the l-th iteration, λ(ω) is the frequency domain form of the Lagrange multiplier, corresponding to the weight term of the constraint condition, used for the constraint condition, is the adjustment parameter, β * is the bandwidth penalty parameter (positive real number), which is used to control the distribution width of the mode in the frequency domain, that is, the degree of separation of the components. ω is the independent variable of frequency domain analysis, which represents the frequency component of the signal after Fourier transform. is the center frequency of the kth eigenmode function after the lth iteration, through Penalize frequency domain components that deviate from the center frequency and force the modal energy to concentrate. represents the frequency domain superposition of all modes except the kth mode in the lth iteration.

[0084] The center frequency update formula is:

[0085]

[0086] in, It represents the center frequency of the kth intrinsic mode function after the l+1th iterative update, reflects the dominant frequency of the current intrinsic mode function component, and is an important parameter for synchronous optimization during the MVMD decomposition process. Represents the frequency domain obtained in the previous step, and its value when the number of iterations is k. It updates the center frequency through frequency domain energy calculation to ensure the frequency alignment of multivariate signals. The generation of the final intrinsic mode function depends on ADMM (alternating direction multiplier method). ADMM iteration is used to minimize the total bandwidth of the modal components. By updating the modal components and center frequencies multiple times, as the iteration proceeds, the compensation effect of the Lagrange multiplier λ(ω) multiplier on the constraints gradually saturates. When the sensitivity of the solution update to the parameter adjustment is reduced to a negligible range (that is, the rate of change of the frequency domain modal components and the rate of change of the center frequency are both less than 1% when the parameters are adjusted), the convergence boundary is reached. The output IMF is the result after multiple optimization iterations, that is, the optimized intrinsic mode function.

[0087] Sample Entropy (SE) calculation and merging: By calculating the sample entropy of each optimized intrinsic mode function to measure the sub-mode complexity, the number of mode components is reduced and the effect of signal decomposition is optimized. The specific process is as follows:

[0088] Calculate sample entropy: Calculate sample entropy SE for each optimized IMF k , which measures its complexity. Sample entropy is calculated according to the formula Calculate the conditional probability estimate of the matching of data subsequences in the time series of each optimized intrinsic mode function, where if the difference between the corresponding data points of the two subsequences is less than the tolerance threshold, they are considered to be matched.

[0089] Where N is the length of the signal sequence of the optimized intrinsic mode function, m is the embedding dimension, r is the tolerance threshold, which is 0.15; A is the logarithmic count of matching m+1-dimensional templates, and B is the logarithmic count of matching m-dimensional templates.

[0090] Merging strategy: If the difference in sample entropy between two IMFs is less than a threshold of 0.1, they are merged into a single mode. If the difference in sample entropy between two IMFs is not less than a threshold of 0.1, they are treated as separate modes. The purpose of merging sub-modes is to reduce the number of sub-modes, avoid excessive component division, improve the overall stability and interpretability of the signal decomposition, and reduce the complexity of modeling.

[0091] After merging, all modes are summarized to obtain multi-mode reconstruction features M1, M2, ..., M n , n is the number of patterns.

[0092] Step S2: Construct an adaptive global convergence benchmark predictor cluster.

[0093] Step S2.1: Reconstruct the multi-modal features M1, M2, ..., M nEach intrinsic mode function in the model establishes a basic fuzzy neural network as a benchmark predictor, which is used to output wind power based on the input characteristic data. The specific construction process includes:

[0094] Step S2.1.1: Construct the fuzzy neural network structure and set the model parameters. The fuzzy neural network is used to map the nonlinear relationship of input features (such as wind speed and power history) to the predicted output. The fuzzy neural network includes an input layer, a membership function layer, a rule layer, and an output layer. The input layer is used to receive the multi-mode reconstruction feature corresponding to the historical window data, which is a feature data sequence over a period of time in the past. Specifically, it is the power component and feature component sequence of each reconstruction mode for τ consecutive time steps before the prediction time t. The historical data window x k (t)=[P k (t-τ),v k (t-τ),...,P k (t),v k (t)], the window length τ is 24, P k (t) is mode M k The power component, V k (t) is the characteristic component corresponding to the kth mode (such as wind speed, temperature, and the original multi-source data collected). Before input, it is necessary to set the window length τ for the base learner training set and independently intercept the historical window of each reconstructed mode; the window contains the power component and characteristic component at the current time t and the previous τ-1 time, and the historical data window x is obtained. k (t)=[P k (t-τ),v k (t-τ),...,P k (t),v k (t)].

[0095] The membership function layer uses Gaussian membership function to transform the clear input data into fuzzy sets and calculate the membership of different modes of input variables in the input data under different fuzzy rules. Specifically, the input variable x is calculated according to the following Gaussian membership function ki The membership degree μ under the jth fuzzy rule kj (x ki ):

[0096]

[0097] Among them, kj is the center (mean parameter) of the jth fuzzy rule corresponding to the kth eigenmode function and is updated with time t. kj is the width (standard deviation) of the jth fuzzy rule corresponding to the kth intrinsic mode function, x ki is the value of the kth intrinsic mode function of the input data at the i-th moment.

[0098] The rule layer is used to calculate the trigger strength of each fuzzy rule and normalize all fuzzy rules to obtain the final trigger strength. Specifically, the trigger strength μ of each rule is calculated as follows: kj (x k (t)), and normalize all fuzzy rules to obtain the final trigger intensity vector The dimension of the final trigger strength vector is equal to the number of fuzzy rules, and the sum of its elements is always 1, which can quantify the contribution of each rule to the current input. The calculation formula is:

[0099]

[0100] Among them, ji (t) is the membership function center of the jth rule to the i-th feature, δ ji (t) is the membership function width of the j-th rule to the i-th feature, x i (t) is the i-th characteristic component of the historical data window. R is the total number of fuzzy rules.

[0101] And use the following formula to weight the output of different fuzzy rules to obtain the prediction results corresponding to different modes

[0102]

[0103] Where: θ kj is the output weight of the jth rule of the kth intrinsic mode function.

[0104] Step S2.1.2: Introduce auxiliary error and construct Lyapunov function. Lyapunov function is used to determine the convergence degree of prediction error during the training process of fuzzy neural network model, so that the model parameters can be continuously updated until the prediction error converges to the minimum. The specific instructions are as follows:

[0105] First, the auxiliary error is introduced:

[0106] definition represents the raw error of the kth base predictor (the base learner prediction value and the true value P k (t)), the auxiliary error for:

[0107]

[0108] Among them, η k (t) is the time-varying learning rate, which represents the historical error weight, It is the fuzzy rule trigger strength vector, which is calculated by the input data through the membership function. Its function is to adaptively adjust the model parameter update direction according to the fuzzy activation state of the current input feature to ensure the global convergence of the model weight. The weight parameter vector of the fuzzy rules updated during online learning. The introduction of the auxiliary error function realizes the function of the loss function in the network, which can minimize the prediction error in each training, thereby adjusting the network parameters of the benchmark predictor and making each base learner reach the optimal model. k (t) is the trigger strength vector associated with the kth reference predictor.

[0109] Then, according to formula e k (t+1)=αe k (t)+Δ k (t)-d(t) constructs the original error dynamic system of the kth benchmark predictor, where d(t) is the approximation error; α is the error dynamic propagation coefficient, which is used to describe the time decay characteristics of the error. Its value is usually set according to the system stability analysis, and the error convergence is ensured by Lyapunov theory. is the auxiliary error correction term, Ξ k (t) is a positive definite matrix that controls the step size of model parameter adjustment (equivalent to the adaptive learning rate).

[0110] Finally, the Lyapunov function is constructed. The predictor adjusts the weight of the auxiliary error by minimizing the prediction error. The role of the Lyapunov function is to ensure the stable convergence of the prediction system by guiding the parameter update direction and to achieve stability by controlling the attenuation of the energy. To verify the global convergence of the error. The update equation of the auxiliary error weight of the baseline predictor is:

[0111]

[0112] Θ k (t+1)=Θ k (t)+ΔΘ k (t)

[0113] Among them, ΔΘ k (t) is the update amount of the auxiliary error weight of the k-th reference predictor, k (t) is the trigger strength vector associated with the kth reference predictor, is the auxiliary error, used to adjust the weight parameter Θ of the baseline predictor k (t), η k (t) is the learning rate, e k (t) is the prediction error, Θ k (t) is the weight vector of the kth benchmark predictor at time t, Θk (t+1) is the weight vector of the kth base predictor at time t+1.

[0114] Step S2.2: Use the base learning training set data obtained in step 1 to train each benchmark predictor. During the training process, dynamically update the weight of the benchmark predictor according to the prediction error; update the model parameters and perform gradient update at each time step;

[0115] During the training process, the model parameters are updated using the following formula, and the gradient is updated at each time step:

[0116]

[0117] Where: represents the model parameters of the benchmark predictor at time step t, and all model parameters are updated after each iteration; L(θ t ) is the loss function; λ is the dynamic weight coefficient of the loss function; η is the learning rate; is the gradient of the loss function with respect to the error distribution parameter θ at time step t; m t is the bias-corrected first-order moment estimate, representing the exponentially weighted average of the gradient; To prevent the formula from having an extremely small constant with a denominator of zero, the value is taken as 0.01;

[0118] v t is the second-order moment estimate, which represents the exponentially weighted average of the square of the gradient; β1 is the decay rate of the first-order moment, ranging from (0,1), which is used to control the influence of momentum; is the bias correction term; β2 is the decay rate of the second-order moment, ranging from (0,1), which is used to control the influence of the square of the gradient.

[0119] Step S2.3: Use the weighted prediction set data obtained in step 1 and make predictions through each benchmark predictor trained in step 2.2, and output the corresponding benchmark power prediction value

[0120] Step S3: constructing a PDF shape controller, which is used to integrate the reference power prediction values ​​output by each reference predictor to obtain a final wind power prediction value.

[0121] The integrated calculation formula is:

[0122] in, represents the benchmark power prediction value output by the kth benchmark predictor; w represents the weight matrix, which is composed of the weight of each benchmark predictor; the number of benchmark predictors is n, and the initial weight of each benchmark predictor is defined as 1 / n.

[0123] Step S4: Based on the error PDF shape control, the probability density function of the integrated error is calculated, and the weight matrix is ​​updated by minimizing the objective function until the weight matrix change rate is less than the set change rate threshold of 0.01 and convergence is reached to obtain the final weight matrix w k Obtain the optimal benchmark predictor and the optimal PDF shape controller. The specific process is:

[0124] Input the reference power prediction value output by each reference predictor in step S2.3 into the constructed PDF shape controller to determine whether the change rate of the weight matrix at adjacent moments is less than the set change rate threshold. If so, obtain the weight matrix w k and the optimal benchmark predictor; otherwise, update the weight matrix and return to step 2.2 to update the weight of each benchmark predictor until the change rate of the weight matrix at adjacent moments meets the requirements. Follow the steps below to calculate and determine the change of the weight matrix at adjacent moments:

[0125] (1) The probability density of the actual integrated error of the PDF shape controller is calculated according to the kernel density formula. The probability density calculation formula of the actual integrated error is:

[0126]

[0127] Among them, the actual integrated error of the PDF shape controller is κ is the kernel bandwidth, representing the sample standard deviation; T is the total number of samples; L is the sample index, representing the Lth sample; P(t) is the actual value of wind power, is the baseline wind power forecast value.

[0128] (2) Design objective function: The objective function is Use the objective function to minimize the probability density f of the actual integration error pdf (e) Probability density of the target integration error f target (e) Differences (e.g., Gaussian distribution).

[0129] (3) Weight optimization of the benchmark predictor: The weight matrix update formula is used to perform sequential quadratic programming (SQP) to update the weight of the benchmark predictor. The weight matrix update formula is:

[0130]

[0131] Where η is the learning rate, is the gradient of the objective function with respect to the weight matrix, calculated using the chain rule. Convergence occurs when the weight change rate is less than 0.01, resulting in the final weight matrix. The final weight matrix is ​​used to determine the masses of each benchmark predictor, the optimal benchmark predictor, and the optimal PDF shape controller.

[0132] Step 5: Input the test set data obtained in step 1 into the optimal benchmark predictor for prediction, and use the optimal PDF shape controller to output the final prediction result to verify the prediction error of the optimal benchmark predictor. If the prediction error meets the design requirements, execute step 6; otherwise, return to step 2 to adjust the model parameters of the benchmark predictor.

[0133] Step 6: Obtain the actual wind farm characteristic data and use the trained and predicted optimal benchmark predictor and optimal PDF shape controller to predict the output wind power.

[0134] The following is a specific example to illustrate the training and prediction process of the optimal benchmark predictor and the optimal PDF shape controller constructed by the method of the present invention, so as to reflect the prediction effect of the method of the present invention.

[0135] (1) Data partitioning: The meteorological and power data of a wind farm from January to December 2019 were selected. The characteristic data included: wind speed (m / s) at the height of 10m, 30m, 50m, and 70m of the wind tower; wind speed (m / s) at the hub height; wind direction at the height of 10m, 30m, 50m, and 70m of the wind tower; wind direction at the hub height; temperature, air pressure, and humidity. The target variable was the wind turbine output power (kW). The sampling interval was 15 minutes, with a total of 35,041 data items. The input data included a 14-dimensional original data set of meteorological information such as power, wind speed, and air pressure, which was divided into a base learning training set, a weight prediction set, and a test set according to the ratio of 7:2:1.

[0136] (2) Hyperparameter tuning: Determine the number of multivariate mode decomposition (MVMD) decompositions K = 7 through grid search and set it according to the signal complexity (recommended value range: 5-12). Set the Lagrange multiplier λ to the default range of [0.001, 100], Figure 5 The eigenmode functions (IMFs) after power signal decomposition are shown. The alternating direction multiplier method (ADMM) is used to iteratively update the IMFs to ensure their frequency-division characteristics. The sample entropy of the components is calculated, with the embedding dimension set to 3 and the tolerance threshold set to 0.02. If the difference in sample entropy between two IMFs is less than the threshold of 0.1, they are merged into the same mode. Otherwise, they are saved as separate modes.

[0137] (3) Model verification: Match the historical 24-hour data window and calculate the Gaussian membership function parameters. Figure 6The core mechanism of fuzzification of the benchmark predictor input is demonstrated. The membership function uses a dynamic Gaussian distribution, where the parameter μ represents the center position of the fuzzy rule (corresponding to typical wind speed values ​​of 5, 15, and 20 m / s), and the parameter σ reflects the rule coverage (set to 2, 3, and 1.5, respectively). The triggering probability is calculated for wind speed inputs within different intervals, achieving flexible partitioning of the feature space. For example, when the input wind speed is 10 m / s, two fuzzy rules, μ = 5 (membership 0.22) and μ = 15 (membership 0.21), are simultaneously triggered and updated online via the Lyapunov auxiliary error. Figure 7 The figure shows the trigger intensity distribution of fuzzy rules in the temperature-wind speed joint feature space. The red area in the figure indicates that when the wind speed is close to 10m / s and the temperature is in the range of 10-20℃, the activation intensity of the specific fuzzy rule reaches above 0.7. This multivariable collaborative activation mechanism effectively expresses the nonlinear coupling effect of wind speed and temperature on power output. For example, in the low wind speed and low temperature zone (wind speed <5m / s, temperature <5℃), it exhibits weak activation characteristics (blue area), which meets the actual wind turbine cut-in wind speed limit conditions. The performance of different weight update strategies was tested on the validation set to select the optimal parameters.

[0138] (4) Online adaptive update of the baseline predictor: Set the kernel bandwidth to 0.06 and the learning rate to 0.01 to obtain the baseline predictor output. Figure 8 The verification of parameter convergence characteristics based on Lyapunov theory is demonstrated. The blue curve represents the decay process of the Lyapunov function value with the number of iterations, and the red dashed line is the convergence threshold line of 0.1. The results show that the function value drops below the threshold when the training reaches the 27th iteration, proving that the adaptive update rule of the present invention has strict stability guarantees. Compared with the traditional gradient descent method (in which stability verification is not performed), the present invention solves the local optimum problem caused by parameter oscillation.

[0139] (5) If the prediction error is greater than the set threshold of 0.05, the weight matrix w of the baseline predictor is updated. Figure 9 This study reveals the optimization effect of the PDF shape controller-based control strategy on the error distribution. The original error (blue) before optimization exhibits asymmetric and heavy-tailed characteristics, with a left-side skewness of 0.85 and a kurtosis of 5.7. After optimization using the PDF shape controller (red), the error distribution effectively converges to near zero mean, with the skewness optimized to 0.12 and the kurtosis reduced to 3.1. In particular, the probability density ratio in the error interval [0.8, 1.5] decreases from 22.3% to 6.7%, verifying that the model constructed in this paper can achieve asymmetry suppression capabilities.

[0140] (6) Fine-tune the parameters of the baseline predictor and the PDF shape controller integration weights. Figure 10As shown in the figure, after applying the error dynamic optimization method of the present invention, the prediction effect is greatly improved. Figure 10 There are four curves, corresponding to the actual observation value, namely the true value (black curve), the prediction result of the method of the present invention (blue curve), the prediction result of the traditional random forest algorithm (orange curve) and the prediction result of the traditional LSTM neural network model (green curve). Figure 10 It can be seen that the wind power curve predicted by the method of the present invention has a higher degree of fit with the true value, the fluctuation trend is basically consistent, and it can better track the actual changes in wind power; while the random forest algorithm and LSTM neural network model prediction methods can reflect the changing law of wind power in the overall trend, but the prediction deviation is large at some moments, especially in the interval where wind power changes sharply. In this interval, the prediction result of the method of the present invention is closer to the true value.

[0141] The above description is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with this technical field can easily think of various equivalent modifications or replacements within the technical scope disclosed in the present invention, and these modifications or replacements should all be included in the scope of protection of the present invention.

Claims

1. An adaptive integrated wind power prediction method based on probability density function shape control, characterized in that: The following steps are involved: Step 1: Obtain the original wind power data and the corresponding feature data to form a multivariate signal X(t); divide the multivariate signal X(t) into a base learning training set, a weight prediction set, and a test set, perform multivariate variational mode decomposition on each data set, obtain multiple intrinsic mode functions, calculate the sample entropy of each intrinsic mode function, merge the intrinsic mode functions whose sample entropy difference is less than the set threshold, and obtain multi-mode reconstruction features M1, M2, ..., M after merging. n , n is the number of patterns; The characteristic data include wind speed, temperature and air pressure; Step 2: Build a cluster of benchmark predictors: Step 2.1: Reconstruct the multi-modal features M1, M2, ..., M n Each intrinsic mode function in the method establishes a basic fuzzy neural network as a reference predictor, wherein the reference predictor is used to output wind power according to input characteristic data; Step 2.2: Use the base learning training set data obtained in step 1 to train each benchmark predictor. During the training process, the weight of the benchmark predictor is dynamically updated according to the prediction error. Update model parameters and perform gradient updates at each time step; Step 2.3: Use the weighted prediction set data obtained in step 1 and make predictions through each benchmark predictor trained in step 2.2, and output the corresponding benchmark power prediction value; Step 3: constructing a PDF shape controller, wherein the PDF shape controller is used to integrate the reference power prediction values ​​output by each reference predictor to obtain a final wind power prediction value; in, represents the benchmark power prediction value output by the kth benchmark predictor; w represents the weight matrix, which is composed of the weight of each benchmark predictor; the number of benchmark predictors is n, and the initial weight of each benchmark predictor is defined as 1 / n; Step 4: Determine the weight matrix w k , we get the optimal benchmark predictor and the optimal PDF shape controller: Input the benchmark power prediction value output by each benchmark predictor in step 2.3 into the constructed PDF shape controller to determine whether the change rate of the weight matrix at adjacent moments is less than the set change rate threshold. If so, the weight matrix w is obtained. k and the optimal benchmark predictor; otherwise, update the weight matrix and return to step 2.2 to update the weight of each benchmark predictor until the change rate of the weight matrix at adjacent moments meets the requirements; The weight matrix update formula is: Where η is the learning rate, is the gradient of the objective function with respect to the weight matrix; Objective function Among them, f pdf (e) is the actual integrated error probability density, f target (e) is the target integrated error probability density; The probability density of the actual integrated error is Where: κ is the kernel bandwidth, representing the sample standard deviation; T is the total number of samples; L is the sample index, representing the Lth sample; e is the ensemble error of the benchmark predictor; is the prediction error of the benchmark predictor at time t, P(t) is the actual value of wind power, is the baseline power prediction value; Step 5: Input the test set data obtained in step 1 into the optimal benchmark predictor for prediction, and use the optimal PDF shape controller to output the final prediction result to verify the prediction error of the optimal benchmark predictor. If the prediction error meets the design requirements, execute step 6; otherwise, return to step 2 to adjust the model parameters of the benchmark predictor. Step 6: Obtain actual wind farm characteristic data and use the optimal benchmark predictor and the optimal PDF shape controller to predict the output wind power.

2. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 1 is characterized in that: The specific process of step 1 is: Step 1.1: Obtain raw wind power data and corresponding feature data to form a multivariate signal X(t), where the time resolution of the multivariate signal X(t) is less than 15 minutes; divide the multivariate signal X(t) into a base learning training set, a weight prediction set, and a test set; Step 1.2: Establish the solution objective and constraints of minimizing the total bandwidth of the modal components in the multivariate mode decomposition, perform multivariate mode decomposition on each data set, and obtain K intrinsic mode functions; The solution objective and constraint formula for minimizing the total bandwidth of the modal components are: Where X(t) represents the original input multivariate signal, u k (t) is the kth eigenmode function, ω k is the center frequency of the kth eigenmode function; Represents the derivative operation with respect to time t; the goal of multivariate variational mode decomposition is to decompose the original input multivariate signal into multiple modal signals; Step 1.3: Use the alternating direction multiplier method to iteratively update the frequency domain modal components and frequency domain centers of each eigenmode function, and output the optimized eigenmode function; Step 1.4: Calculate the sample entropy of each optimized intrinsic mode function, compare the size of each sample entropy, merge the optimized intrinsic mode functions with sample entropy difference less than 0.1 into the same mode, and take the optimized intrinsic mode functions with sample entropy difference not less than 0.1 as separate modes, and obtain multi-mode reconstruction features M1, M2, ..., M after merging. n , n is the number of patterns.

3. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 2 is characterized in that: The process of step 1.3 is as follows: The frequency domain modal component and center frequency of each intrinsic mode function are iteratively updated using the alternating direction multiplier method. When the frequency domain modal component change rate and center frequency change rate of each intrinsic mode function are less than 1%, the iterative update ends and the optimized intrinsic mode function is output. The frequency domain modal component update formula is: in, is the frequency domain of the kth eigenmode function in the l+1th iteration, is the synthetic signal of the original frequency domain signal, which serves as the reference input for modal update; is the frequency domain of the i-th intrinsic mode function in the l-th iteration; λ(ω) is the frequency domain form of the Lagrange multiplier, which corresponds to the weight term of the constraint condition and is the adjustment parameter; β * is the bandwidth penalty parameter, which is a positive real number and is used to control the distribution width of the mode in the frequency domain; is the center frequency of the kth eigenmode function after the lth iteration; Penalize the frequency domain components that deviate from the center frequency and force the modal energy to aggregate; ω is the independent variable of frequency domain analysis, which represents the frequency component of the signal after Fourier transform; represents the frequency domain superposition of all eigenmode functions except the kth eigenmode function in the lth iteration; The center frequency update formula is: Where, represents the center frequency of the kth intrinsic mode function after the l++1th iteration update; represents the frequency domain obtained by the last iterative update, that is, the frequency domain of the kth intrinsic mode function in the l+1th iteration.

4. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 3 is characterized in that: In step 1.4, the conditional probability estimates of the data subsequence matching in the time series of each optimized intrinsic mode function are calculated respectively, and the obtained conditional probability estimates are used as the sample entropy corresponding to each optimized intrinsic mode function; The conditional probability estimate is calculated as follows: Where N is the signal sequence length of the optimized intrinsic mode function, m is the embedding dimension, r is the tolerance threshold, which is 0.15, A is the logarithmic count of matching m+1-dimensional templates, and B is the logarithmic count of matching m-dimensional templates.

5. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 1 is characterized in that: In step 2.1, the multi-modal reconstruction features M1, M2, ..., M n The process of establishing a baseline predictor for each eigenmode function in is: Step 2.1.1: Construct a fuzzy neural network model and set the model parameters of the fuzzy neural network model; the fuzzy neural network includes an input layer, a membership function layer, a rule layer and an output layer; The input layer is used to receive the historical window data x corresponding to the multi-modal reconstruction feature k (t)=[P k (t-τ),v k (t-τ),...,P k (t),v k (t)], where τ is the window length; P k (t) is the power component in the reconstructed feature of the kth mode, V k (t) is the feature component in the reconstructed feature set of the kth pattern; The membership function layer adopts Gaussian membership function to convert clear input data into fuzzy sets and calculate the membership of input variables of different modes in the input data under different fuzzy rules; The rule layer is used to calculate the trigger strength of each fuzzy rule and normalize all fuzzy rules to obtain the final trigger strength; The output layer weights and fuses the output results of multiple fuzzy rules to obtain a wind power prediction result; Step 2.1.2: Introduce auxiliary error and construct Lyapunov function; the Lyapunov function is used to judge the convergence degree of the prediction error during the training process of the fuzzy neural network model, so that the model parameters can be continuously updated until the prediction error converges to a minimum.

6. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 5 is characterized in that: In step 2.2, the weight of the auxiliary error in each benchmark predictor is updated using the following rule: Among them, ΔΘ k (t) is the update amount of the auxiliary error weight of the k-th reference predictor; k (t) is the trigger strength vector with the kth reference predictor, is the auxiliary error; η k (t) is the learning rate, e k (t) is the prediction error, is the current weight estimate of the k-th base predictor; is the final trigger strength vector.

7. The adaptive integrated wind power prediction method based on probability density function shape control according to claim 6, characterized in that: In step 2.2, the following formula is used to update the model parameters during training, and the gradient is updated at each time step: Where: represents the model parameters of the benchmark predictor at time step t, and all model parameters are updated after each iteration; L(θ t ) is the loss function; λ is the dynamic weight coefficient of the loss function; η is the learning rate; is the gradient of the loss function with respect to the error distribution parameter θ at time step t; m t is the bias-corrected first-order moment estimate, representing the exponentially weighted average of the gradient; To prevent the formula from having an extremely small constant with a denominator of zero, the value is taken as 0.01; v t is the second-order moment estimate, which represents the exponentially weighted average of the square of the gradient; β1 is the decay rate of the first-order moment, ranging from (0,1), which is used to control the influence of momentum; is the bias correction term; β2 is the decay rate of the second-order moment, ranging from (0,1), which is used to control the influence of the square of the gradient.

Images