Summary of the invention
Goal of the invention: aiming at the problems existing in the prior art, the present invention provides a kind of based on optimal VMD and Synchronous fluorimetry
Ultra-short term wind speed forecasting method, it can be considered that mutual shadow when model parameter and mode input variable separately optimizing between the two
It rings, while obtaining optimal model parameter and being combined with input variable, the total optimization of implementation model parameter and input variable;It can
Reduce wind speed time series it is strong it is non-stationary predicted value must be influenced, and then improve precision of prediction.
Technical solution: the present invention provides a kind of ultra-short term wind speed forecasting method based on optimal VMD and Synchronous fluorimetry, packet
It includes following steps: step 1: obtaining wind field history actual measurement air speed data, air speed data is surveyed according to history and establishes wind speed time sequence
Column, are divided into training sample and test sample for wind speed time series;Step 2: determine that VMD is decomposed using centre frequency observation
Parameter K and the parameter τ that criterion determines that optimal VMD decomposes is minimized using residual error evaluation index, by wind speed Time Series
Submodule state with different center frequency;Step 3: each submodule state resulting to step 2 is normalized to [0,1] area
Between, the partial autocorrelation function value of submodule state is calculated to determine initial candidate input variable;Step 4: it to each submodule state, uses
Mixing backtracking searching algorithm synchronizes optimization to prediction model candidate input variable and parameter, and wherein binary system backtracking search is calculated
Each candidate input variable that method is used to identify partial autocorrelation function carries out postsearch screening, and real number is recalled searching algorithm and is used for pole
The parameter of limit learning machine prediction model optimizes;Step 5: the resulting optimal model parameters of step 4, optimal input are gathered
And test sample input substitutes into extreme learning machine prediction model, the predicted value and renormalization of test phase is obtained, to all
The prediction result of submodule state carries out summation reconstruct, obtains the predicted value of wind speed time series;Step 6: root-mean-square error is used
(RMSE), mean absolute error (MAE), four evaluation indexes of average absolute percentage error (MAPE) and related coefficient (R) are judged
The accuracy of "current" model.
Further, in the step 2, the corresponding center of each submodule state under different decomposition mode number is calculated first
Frequency enables mode decomposition number K=1,2 ..., k decompose respectively to original signal, if occurring similar as K=k
Submodule state centre frequency, then optimal Decomposition mode number is selected as k-1;Then, it calculates original signal under different update parameter τ and goes
The root-mean-square error of noise cancellation signal is denoted as residual error evaluation index REI, chooses the smallest τ value of REI as VMD optimal Decomposition parameter.
Preferably, the residual error evaluation index REI can be described as follows:
In formula, N indicates the length for most indicating original signal, uk(i) k-th of submodule state, i-th of sample point, f (i) table are indicated
Show i-th of sample point of original signal.
Further, in the step 4, using the parameter of mixing backtracking searching algorithm Synchronous fluorimetry extreme learning machine
With the input matrix of prediction model, each of mixed population individual includes R real-valued parameter and B binary features mask,
Essential parameter is made of the weight and threshold value of extreme learning machine model, in the corresponding candidate input set of each binary mask
One input vector, the encoded radio of binary features mask are " 1 ", then the candidate input vector of the position it is selected finally enter to
Duration set, encoded radio are " 0 ", then delete the input vector.
Preferably, fitness function is selected as the root-mean-square error between measured value and predicted value:
In formula, fitness indicates adaptive value, NtIndicate the length of training sample, NoIndicate the output of extreme learning machine model
Vector number, FijFor model predication value, OijFor measured value.
Further, in the step 4, extreme learning machine parameter is optimized using real number backtracking searching algorithm
The specific implementation steps are as follows:
(1) initialize: enabling Population Size is N, maximum number of iterations G, initializes the number of iterations k=1.According to formula (8)
Random to generate initial population P, enabling i-th of individual in kth time iteration is Pi(k), i=1,2 ..., N, k=1,2..., G.Kind
The initialization procedure of group is as follows:
PI, j=rand* (uPj-lowj)+lowj, i=1,2 ..., N, j=1,2 ..., D (8) formula in, D indicate kind
Group's dimension, upjAnd lowjThe upper bound and the lower bound of jth dimension component are respectively indicated, rand~U (0,1), U () are to be uniformly distributed, PI, j
Indicate that the ground j of i-th of individual of initial population ties up variable.
(2) select I: history population Old P generated according to formula (9) at random, according to formula (10) and (11) to history population into
Row selection updates:
oldPI, j=rand* (upj-lowj)+lowj (9)
If a < b then oldP=P, else oldP=oldP (10)
OldP=randshuffle (oldP) (11)
In formula, oldPI, jIndicate that i-th of individual ground j of history population ties up variable, a, b, which indicate to obey on (0,1), uniformly to be divided
The random number of cloth, randshuffle () indicate random selection function;In formula (11), whether random selection is planted current parent generation
Group is substituted for some randomly selected history population, realizes the memory function of backtracking optimization algorithm.
(3) it makes a variation: it is randomly ordered to the individual progress of parent population, then make a variation according to equation (13):
OldP=Permuting (oldP) (12)
Mutant=P+F (oldP-P) (13)
In formula, Permuting () indicates that randomly ordered function, F are mutation scaling coefficient, is used to command deployment direction square
The variation amplitude of battle array, Mutant are the population after variation, and the default setting of usual F is F=3rand, rand~N (0,1), N
() is standardized normal distribution;
(4) intersect: firstly generating the binary matrix map of N*D dimension, and pass through two random numbers rand1 and rand2
Equiprobability calls two kinds of pre-set Crossover Strategies, then generates new experimental population T, works as mapijWhen=1, experimental population
In individual TijBy PijReplacement;Steps are as follows for specific intersection:
T=map.*P+ (~map) .*Mutant (15)
In formula, mixrate is crossover probability, rand, rand1 and rand2 equally distributed random to obey on (0,1)
Number, randi indicate to generate a random integers, the random integers between one 1~D of randi (D) expression generation;
(5) it selects II: calculating the fitness value f of current population at individual according to fitness functioni(k), and kind of lower generation is updated
Group, If fi(k) < fbest, then fbest=fi(k), Pbest=Pi(k);Wherein, fbestIndicate current adaptive optimal control value, Pbest
Indicate current optimal population;
(6) if k < G, enables k=k+1, go to step (2), and otherwise output limit learning machine prediction model is optimal defeated
Enter weight and threshold value and optimal input set.
Further, in the step 4, in binary system backtracking searching algorithm, the initialization of population, selection, intersection
It is identical with operating process and real number the backtracking searching algorithm such as variation, but the calculation of adaptive value is different, in population often each and every one
Body is encoded as a binary vector, and is converted the fitness function value of individual to [0,1] sky by sigmoid function
Between:
In formula, k indicates the number of iterations, fiIt (k) is i-th of ideal adaptation angle value of population, SiIt (k) is the population after conversion
Body fitness value;The encoded radio BP of binary system population at individuali(k) it updates as follows:
Preferably, in the step 2, each submodule state u is estimated by following 3 stepskBandwidth: S1. passes through
Hilbert transformation solves each submodule state ukAnalytic signal, further obtain the unilateral frequency spectrum of each submodule state;S2.
By each submodule state ukCorresponding centre frequency ωkExponential term aliasing, further the frequency spectrum of submodule state is converted to base band
Region;S3. the bandwidth of each submodule state is estimated by calculating square L2 norm of demodulated signal gradient.
Preferably, in the S3, corresponding constraint variation problem can be described as follows:
In formula,It indicates to carry out the time derivation, t is the time, and δ (t) indicates that impulse function, j indicate imaginary unit, *
Indicate convolution algorithm, ukWith wkIt respectively decomposes and obtains the time-domain signal and centre frequency of k-th of mode.Introduce secondary penalty term
Unconstrained Optimization Problem further is converted by constrained optimization problem above with Lagrange multiplier:
In formula, λ is Lagrange multiplier, and α is the punishment parameter of quadratic term;Above-mentioned nothing is solved using alternating direction multipliers method
The basic decomposition step of constrained optimization problem, VMD can be expressed as follows: (1) being initializedInitialize iteration
Frequency n=1;(2) to each submodule state uk, according to alternating direction multipliers method, ukRenewal process it is as follows:
In formula, ω is random frequency, ωkFor submodule state ukCentre frequency,Become for the Fourier of signal f (ω)
It changes,For the Fourier transformation of λ (ω),It indicatesFourier transformation;(3) centre frequency ωkRenewal process such as
Under:
In formula,Indicate ukThe Fourier transformation of (ω);(4) renewal process of Lagrange multiplier λ is as follows:
In formula, τ is that step-length updates coefficient;(5) if(ε > 0 be discrimination precision), then iteration
Stop;Otherwise, n=n+1 is enabled, gos to step 2.
Preferably, in the step 6, the calculation formula of four evaluation indexes is as follows:
In formula, fsIt (i) is the analogue value of i-th of sample;foIt (i) is the measured value of i-th of sample;WithIt respectively indicates
The averaging analog value of sample data set and average measured value;The size of N expression sample set.
The utility model has the advantages that being directed to strong randomness, fluctuation and the intermittent feature of wind speed, the invention proposes one kind to be based on
The ultra-short term wind speed forecasting method and system of optimal VMD and Synchronous fluorimetry, are firstly introduced into based on centre frequency observation and residual error
Original wind speed Time Series are had the submodule state of different center frequency by the optimal VMD that evaluation index minimizes criterion,
On the basis of this, secondary sieve is carried out to each candidate input variable that partial autocorrelation function identifies using binary system backtracking searching algorithm
Choosing, while the parameter of extreme learning machine is optimized using real number backtracking searching algorithm, realize extreme learning machine model
The Synchronous fluorimetry of candidate input variable and parameter realizes wind speed finally, carrying out summation reconstruct to the prediction result of all submodule states
The high-precision forecast of time series.
In general, through the invention it is contemplated above technical scheme is compared with the prior art, have below beneficial to effect
Fruit:
1) present invention optimizes the parameter and input variable of model using the thought of Synchronous fluorimetry simultaneously, it can be considered that
Influencing each other between the two when model parameter and mode input variable separately optimizing, at the same obtain optimal model parameter with it is defeated
Enter variable combination, the total optimization of implementation model parameter and input variable.
2) it is directed to fluctuation, randomness and the intermittence of wind speed time series, the present invention uses optimal VMD by non-stationary wind
Fast Time Series are several relatively stable submodule states, and are predicted respectively sub- mode, are finally reconstructed cumulative
The wind speed time series arrived can reduce the strong of wind speed time series and non-stationary obtain shadow to predicted value as final predicted value
It rings, and then improves precision of prediction.
Specific embodiment
The present invention is described in detail with reference to the accompanying drawing.
Present embodiment carries out example using the 10min wind speed time series of the wind field of inner mongolia two as embodiment
Emulation, to verify effect of the invention.Fig. 1 is that the ultra-short term wind speed provided by the invention based on optimal VMD and Synchronous fluorimetry is pre-
Model flow figure is surveyed, implementation steps are as follows:
Step 1: the 10min air speed data for choosing the wind field of inner mongolia two is as sample data (600 samples
Notebook data point), two wind speed time serieses are established, for the two wind speed time serieses (data set 1 and data set 2), respectively
Using preceding 470 sample number strong points as training sample, rear 130 sample number strong points are as test sample.
Step 2: the parameter K and use residual error evaluation index minimum that optimal VMD is decomposed are determined using centre frequency observation
Change criterion and determine the parameter τ that optimal VMD is decomposed, the wind speed Time Series of two wind fields are had to the son of different center frequency
Mode.Parameter τ influences optimal VMD decomposition result such as Fig. 3 as shown in Fig. 2, two air speed data collection to VMD decomposition REI value
It is shown.
Variation mode decomposition (VMD) is that one-dimensional non-stationary signal f (t) is resolved into K finite bandwidth submodule by a kind of trial
State ukThe onrecurrent signal processing technology of (k=1,2 ..., K), each variation mode have different centre frequency ωk.It is logical
It crosses following 3 steps and estimates each submodule state ukBandwidth:
1. being converted by Hilbert and solving each submodule state ukAnalytic signal, further obtain each submodule state
Unilateral frequency spectrum;
2. by each submodule state ukCorresponding centre frequency ωkExponential term aliasing, further by the unilateral of submodule state
Frequency spectrum is converted to baseband region;
3. estimating the bandwidth of each submodule state by square L2 norm for calculating demodulated signal gradient.Corresponding constraint variation
Problem can be described as follows:
In formula,It indicates to carry out the time derivation, t is the time, and δ (t) indicates that impulse function, j indicate imaginary unit, *
Indicate convolution algorithm, ukWith wkIt respectively decomposes and obtains the time-domain signal and centre frequency of k-th of mode.
It introduces secondary penalty term and Lagrange multiplier and further converts no constraint for constrained optimization problem above
Optimization problem:
In formula, λ is Lagrange multiplier, and α is the punishment parameter of quadratic term.
Above-mentioned Unconstrained Optimization Problem is solved using alternating direction multipliers method, the basic decomposition step of VMD can be stated such as
Under:
(1) it initializesInitialize the number of iterations n=1.
(2) to each submodule state uk, according to alternating direction multipliers method, ukRenewal process it is as follows:
In formula, ω is random frequency, ωkFor submodule state ukCentre frequency,Become for the Fourier of signal f (ω)
It changes,For the Fourier transformation of λ (ω),It indicatesFourier transformation.
(3) centre frequency ωkRenewal process it is as follows:
In formula,Indicate ukThe Fourier transformation of (ω),
(4) renewal process of Lagrange multiplier λ is as follows:
In formula, τ is that step-length updates coefficient.
(5) if(ε > 0 be discrimination precision), then iteration stopping;Otherwise, n=n+1 is enabled, is jumped
Go to step 2.
The conventional sub- mode decomposition number of VMD method is difficult to determine, easily falls into modal overlap and owe to decompose, for this purpose, first
The corresponding centre frequency of each submodule state under different decomposition mode number is calculated, mode decomposition number K=1,2 ..., k are enabled, respectively
Original signal is decomposed, if there is similar submodule state centre frequency, then optimal Decomposition mode number is selected as K=k
k-1;Then, the root-mean-square error for calculating original signal and denoised signal under different update parameter τ, is denoted as residual error evaluation index
(Residual Evaluation Index, REI) chooses the smallest τ value of REI as VMD optimal Decomposition parameter, residual error evaluation
Index REI can be described as follows:
In formula, N indicates the length for most indicating original signal, uk(i) k-th of submodule state, i-th of sample point, f (i) table are indicated
Show i-th of sample point of original signal.
Step 3: each submodule state resulting to step 2 is normalized to [0,1] section, calculates the inclined of submodule state
Auto-correlation function value is to determine initial candidate input variable.
Step 4: to each submodule state, using mixing backtracking searching algorithm to prediction model candidate input variable and parameter
Optimization is synchronized, each candidate input variable that wherein binary system backtracking searching algorithm is used to identify partial autocorrelation function carries out
Postsearch screening, real number backtracking searching algorithm is for optimizing the parameter of extreme learning machine prediction model.
Using the parameter of mixing backtracking searching algorithm Synchronous fluorimetry extreme learning machine and the input matrix of prediction model, mixing
Each of population individual includes R real-valued parameter and B binary features mask, and essential parameter is by extreme learning machine model
Weight and threshold value composition, the corresponding candidate input vector inputted in set of each binary mask, binary features are covered
The encoded radio of code is " 1 ", then the candidate input vector of the position is selected finally enters vector set, and encoded radio is " 0 ", then deletes
The input vector.Fitness function is selected as the root-mean-square error between measured value and predicted value:
In formula, fitness indicates adaptive value, NtIndicate the length of training sample, N.Indicate the defeated of extreme learning machine model
Outgoing vector number, FijFor model predication value, QijFor measured value.
Optimize that the specific implementation steps are as follows to extreme learning machine parameter using real number backtracking searching algorithm:
(1) it initializes.Enabling Population Size is N, maximum number of iterations G, initializes the number of iterations k=1.According to formula (8)
Random to generate initial population P, enabling i-th of individual in kth time iteration is Pi(k), i=1,2 ..., N, k=1,2..., G.Kind
The initialization procedure of group is as follows:
PI, j=rand* (upj-lowj)+lowj, i=1,2 ..., N, j=1,2 ..., D (8)
In formula, D indicates population dimension, upjAnd lowjRespectively indicate jth dimension component the upper bound and lower bound, rand~U (0,
1), U () is to be uniformly distributed, PI, jIndicate that the ground j of i-th of individual of initial population ties up variable.
(2) I is selected.History population OldP is generated at random according to formula (9), and history population is carried out according to formula (10) and (11)
Selection updates.
oldPI, j=rand* (upj-lowj)+lowj (9)
If a < b then oldP=P, else oldP=oldP (10)
OldP=randshuffle (oldP) (11)
In formula, oldPI, jIndicate that i-th of individual ground j of history population ties up variable, a, b, which indicate to obey on (0,1), uniformly to be divided
The random number of cloth, randshuffle () indicate random selection function.In formula (11), whether random selection is planted current parent generation
Group is substituted for some randomly selected history population, realizes the memory function of backtracking optimization algorithm.
(3) it makes a variation.It is randomly ordered to the individual progress of parent population, then make a variation according to equation (13):
OldP=Permuting (oldP) (12)
Mutant=P+F (oldP-P) (13)
In formula, Permuting () indicates that randomly ordered function, F are mutation scaling coefficient, is used to command deployment direction square
The variation amplitude of battle array, Mutant are the population after variation, and the default setting of usual F is F=3rand, rand~N (0,1), N
() is standardized normal distribution.
(4) intersect.The binary matrix map of N*D dimension is firstly generated, and passes through two random numbers rand1 and rand2
Equiprobability calls two kinds of pre-set Crossover Strategies, then generates new experimental population T, works as mapijWhen=1, experimental population
In individual TijBy PijReplacement.Steps are as follows for specific intersection:
T=map.*P+ (~map) .*Mutant (15)
In formula, mixrate is crossover probability, rand, rand1 and rand2 equally distributed random to obey on (0,1)
Number, randi indicate to generate a random integers, the random integers between one 1~D of randi (D) expression generation.
(5) II is selected.The fitness value f of current population at individual is calculated according to fitness functioni(k), and kind of lower generation is updated
Group, If fi(k) < fbest, then fbest=fi(k), Pbest=Pi(k).Wherein, fbestIndicate current adaptive optimal control value, Pbest
Indicate current optimal population.
(6) if k < G, enables k=k+1, go to step (2), and otherwise output limit learning machine prediction model is optimal defeated
Enter weight and threshold value and optimal input set.
In binary system backtracking searching algorithm, the operating process such as initialization, selection, intersection and variation of population and real number are returned
The searching algorithm that traces back is identical, but the calculation of adaptive value is different, and each individual is encoded as a binary vector in population,
And the fitness function value of individual is converted to [0,1] space by sigmoid function:
In formula, k indicates the number of iterations, fiIt (k) is i-th of ideal adaptation angle value of population, SiIt (k) is the population after conversion
Body fitness value.
The encoded radio BP of binary system population at individuali(k) it updates as follows:
Step 5: the resulting optimal model parameters of step 4, optimal input set and test sample input are substituted into pole
Learning machine prediction model is limited, the predicted value and renormalization of test phase is obtained, the prediction result of all submodule states is asked
And reconstruct, obtain the predicted value of the wind speed time series of two wind fields.
Step 6: the accuracy of "current" model is judged using 4 evaluation indexes.This 4 indexs include: root-mean-square error
(RMSE), mean absolute error (MAE), average absolute percentage error (MAPE) and related coefficient (R).Evaluation index calculation formula
It is as follows:
In formula, fsIt (i) is the analogue value of i-th of sample;foIt (f) is the measured value of i-th of sample;WithIt respectively indicates
The averaging analog value of sample data set and average measured value;The size of N expression sample set.
The ultra-short term forecasting wind speed based on optimal VMD Yu Synchronous fluorimetry (OVMD-HBSA-ELM) mentioned using the present invention
Model predicts the wind speed time series of two wind fields, in order to verify the validity of institute's climbing form type of the present invention, by itself and biography
System ELM model and the ELM model (HBSA-ELM) for using mixing backtracking searching algorithm to optimize compare, air speed data collection 1
The error criterion statistical value difference for testing the prediction result of phase with data set 2 is as shown in Table 1 and Table 2.
1 data set of table, 1 model prediction resultant error statistics
2 data set of table, 2 model prediction resultant error statistics
As shown in Table 1, for two datasets, RMSE, MAE and the MAPE of the mentioned HBSA-ELM model of the present invention compare
ELM model is small, and R ratio ELM model is big, illustrates that the prediction effect of HBSA-ELM model is substantially better than ELM model, using mixing back
Optimization algorithm of tracing back, which synchronizes optimization to model parameter and input variable, can be improved the precision of prediction of ELM model;Compare OVMD-
The prediction result of HBSA-ELM model and HBSA-ELM model can be seen that the evaluation of proposed OVMD-HBSA-ELM model herein
Index is substantially better than HBSA-ELM model, and illustrating that optimal VMD is decomposed can be effectively relatively flat by wind speed Time Series
Steady submodule state can significantly improve model prediction essence by the way of predict to sub- mode then summation reconstruct respectively
Degree.
Test phase prediction result residual plot difference is as shown in Figure 4 and Figure 5, by Fig. 4 and Fig. 5 it is found that using OVMD-HBSA-
The residual error that ELM model predicts wind speed time series is minimum, HBSA-ELM model secondly, the residual error of ELM model is maximum,
Further demonstrate the validity of the proposed method of the present invention.
The technical concepts and features of above embodiment only to illustrate the invention, its object is to allow be familiar with technique
People cans understand the content of the present invention and implement it accordingly, and it is not intended to limit the scope of the present invention.It is all according to the present invention
The equivalent transformation or modification that Spirit Essence is done, should be covered by the protection scope of the present invention.