A wind speed vector prediction method based on a Gaussian process and a combination model

Abstract

The application discloses a wind speed vector prediction method based on a Gaussian process and a combination model, which predicts future average wind speed and wind direction according to a wind speed sequence and related features. The method designs a wind speed vector decomposition method, decomposes the wind speed sequence through singular spectrum analysis, fits a periodic term by using a difference autoregressive moving average model, predicts a trend term by using a combination model, and finds model weights by using a weight-corrected particle swarm algorithm, predicts a residual term by using a Gaussian process based on a time window sliding, and finally obtains predicted wind speed and wind direction through a vector synthesis method. When evaluating the performance of the model, the periodicity of the direction and the inapplicability of the traditional indicators are considered, and a new indicator is used to evaluate the performance of the model. The application adopts the idea of "vector" to process wind speed and wind direction, adopts the idea of "combination" to make prediction, and has good prediction ability and generalization ability. The new indicator can better reflect the prediction ability of the model and has stronger physical meaning.

Description

TECHNICAL FIELD

[0001] The application discloses a wind speed vector prediction method based on a Gaussian process and a combination model, which predicts the average wind speed and wind direction size in the future 1 minute according to a wind speed sequence and related features. The method is based on a data set collected from a supervisory control and data acquisition (SCADA) system during the operation of a wind turbine, a wind speed vector decomposition method is designed and related features are constructed, real-time short-term wind speed prediction is realized through time series decomposition, combination prediction, residual prediction and vector synthesis, and when evaluating the performance of the model, the periodicity of the direction and the inadaptability of the traditional indicators are considered, and new indicators RMSE θ , MAPE θ and MAPE' are used to evaluate the performance of the model. BACKGROUND

[0002] With the increasing maturity of new energy generation technology, renewable and clean energy such as wind energy and solar energy have attracted more and more researchers worldwide and gradually replaced part of the traditional fossil energy. Due to the uncertainty, volatility and suddenness of wind energy, the application of wind power generation technology also faces various challenges. In the short term, the rapid and frequent change characteristics of wind speed will affect the dynamic load of mechanical parts of wind turbines, especially when the control system cannot track the change of wind speed in real time, resulting in reduced service life and frequent failures of the unit. On the other hand, the delay of the control system in tracking the change of wind speed leads to the decline of the power generation efficiency of the wind farm. Finally, accurate and rapid wake modeling, efficient and accurate wind farm level power optimization scheduling, and even wind farm site selection and layout all rely on the support of short-term wind speed prediction technology.

[0003] In the background of wind farm, a common method for real-time short-term wind speed prediction is to use AR, MA, ARIMA and other statistical learning related knowledge to make regression prediction only for wind speed and wind direction, but it cannot use more related features and the feature dimension is limited. Or use some common machine learning models for fast wind speed and wind direction regression fitting, but the sudden change characteristics of wind direction at the two edges of 0 and 360 degrees lead to poor fitting effect. Or use some common deep time series prediction networks such as LSTM, but the fixed training set may lead to larger and larger prediction errors in engineering use. Therefore, it is of great significance to design a model with high prediction accuracy and good generalization ability for real-time short-term wind speed prediction. SUMMARY

[0004] The method is based on a data set collected from a supervisory control and data acquisition (SCADA) system during operation of a wind turbine, a wind speed vector decomposition method is designed and relevant features are constructed, a singular spectrum analysis (SSA) is used to decompose the wind speed sequence into a trend item and a periodic item, a difference autoregressive moving average model (ARIMA) is used to fit the periodic item, a combination model composed of an extreme gradient boosting (XGB), a gated recurrent unit (GRU) and a support vector regression (SVR) is used to predict the trend item, a weight-corrected particle swarm algorithm is used to find the weight of the combination model, a time window is constructed, a Gaussian process based on the sliding of the time window is used to predict the residual item, the periodic item, the trend item and the residual item are added to obtain the prediction result of the model, and finally the prediction result is converted into the predicted wind speed and wind direction by a vector synthesis method. When evaluating the performance of the model, considering the periodicity of the direction and the inadaptability of the traditional indicators, new indicators RMSE θ , MAPE θ , MAPE' are used to evaluate the performance of the model. The application uses the idea of "vector" to process the wind speed and wind direction, uses the idea of "combination" to make predictions, and uses the "self-learning" method to predict the residual item, has good prediction ability and generalization ability, and the new indicators can better reflect the prediction ability of the model and have stronger physical meaning.

[0005] The object of the application is achieved by the following technical solutions: a wind speed vector prediction method based on a Gaussian process and a combination model, the method comprising the following steps:

[0006] Step 1, selecting a wind turbine, obtaining N pieces of operation data recorded in the SCADA system of the wind turbine, and dividing the operation data into a training set, a validation set and a test set;

[0007] Step 2, vector decomposition, converting the wind speed and wind direction obtained at each time point in the SCADA system into X-direction wind speed and Y-direction wind speed;

[0008] Step 3, constructing features, obtaining the wind speed sequence processed in step 2, calculating the average wind speed before each time node and taking it as a new feature, calculating the average wind speed 1 min after each time node and taking it as a prediction label. The new features, the prediction label and other features obtained by the SCADA are standardized on the training set, wherein the other features obtained by the SCADA include 30s average wind alignment error, 3s average wind alignment error and current nacelle position. The mean and variance of each feature standardization are saved, and then input to the model for training;

[0009] Step 4, training phase, using singular spectrum analysis (SSA) to decompose the prediction label of the training set into different component sequences, since the frequency of the decomposed component sequence is basically in an increasing relationship, the first term is selected as the trend term, and the sum of the second term and the subsequent terms is selected as the periodicity, and the remaining terms are ignored as noise;

[0010] Step 5, model training for periodic terms, a difference autoregressive moving average model (ARIMA) is selected for fitting, when selecting the number of autoregressive terms p and the number of moving average terms q, first use ACF and PACF plots to select the range of p and q, and then according to the BIC criterion, grid search is used to determine the specific values of p and q to prevent the model from being too complex. After determining the parameters, input the training set data to obtain the fitted periodic term data ytrain arima-X and ytrain arima-Y ;

[0011] Step 6, model training for trend terms, using the standardized input features in step 3, three models of extreme gradient boosting (XGB), gated recurrent unit (GRU), and support vector regression (SVR) are selected for fitting, and some hyperparameters are determined using the validation set. After determining the parameters, input the training set data to obtain the fitted trend term data ytrain XGB-X , ytrain GRU-X , ytrain SVR-X , and the fitted trend term data ytrain XGB-Y , ytrain GRU-Y , ytrain SVR-Y ;

[0012] Step 7, determine the weight of the combined model, use the fitted data obtained in step 5 and step 6, and adopt the weight correction particle swarm algorithm to find the model weight on the validation set, while ensuring that the weight ω ARIMA of the difference autoregressive moving average model is 1, the weight ω XGB of the extreme gradient boosting model, the weight ω GRU of the gated recurrent unit model, and the weight ω SVR of the support vector regression model satisfy ω XGB + ω XRU + ω SVR = 1. In the process of finding the model weight, the modulus of the difference between the prediction vector and the fitted vector is used as the loss function in the finding process, and the weight corresponding to the minimum is selected. Through the combined model calculation formula, the fitted data ytrain X , ytrain Y of X and Y independent models are obtained.

[0013] Step 8, using the mean and variance in the standardization process in step 2, the fitting data obtained in step 7 is de-standardized, and the residual of the de-standardized data and the predicted label data in step 3 is calculated, and the input features standardized in step 3 are used to construct a time window of 60s before each time node, and the data in the time window is used to fit the Gaussian process, and the model parameters are directly determined in the training set, the mean absolute error (MAE) of the residual fitting is selected as the loss function, and the hyperparameters are searched by grid method;

[0014] Step 9, test phase, using the mean and variance in the standardization process in step 3, the input feature data of the test set is standardized and input into the model for prediction;

[0015] Step 10, using the periodic term model trained in step 5, predict the periodic term data of the test set X, Y independent model Using the trend term model trained in step 6, predict the three trend term data of the test set X independent model And the three trend term data of the Y independent model Using the combined model trained in step 7, the combined model data is obtained by the combined model calculation formula, and the combined model prediction output is obtained by de-standardizing the combined model data using the mean and variance in the standardization process in step 3

[0016] Step 11, using the time window sliding Gaussian model trained in step 8, predict the residual data e of the test set X, Y independent model X Y Add the combined model prediction output obtained in step 10 to the residual data to obtain the wind speed prediction value in X direction And the wind speed prediction value in Y direction That is

[0017] Step 12, using the X, Y direction prediction values calculated in step 11, the calculation method of vector composition is used to obtain the predicted wind speed v and wind direction θ;

[0018] Step 13, using the predicted wind speed and wind direction obtained in step 12 and the real wind speed and wind direction of the test set, calculate the new index RMSE θ And MAPE θ , the modified index MAPE' to avoid zero division phenomenon, and the traditional regression evaluation index MAPE, RMSE, all indexes are better and better, the model prediction performance is measured by each index, and the advantages of the new index are shown by comparing the new index with the traditional index. ​

[0019] Furthermore, in step 2, the vector decomposition calculation method is as follows:

[0020] v X =vcosθ

[0021] v Y =vsinθ

[0022] Where v is the wind speed obtained by the SCADA system at each moment, and θ is the wind direction obtained by the SCADA system at each moment. X v Y These represent the wind speed in the X direction and the wind speed in the Y direction, respectively.

[0023] Furthermore, in step 3, the feature standardization calculation method is as follows:

[0024]

[0025]

[0026]

[0027] Where g represents the data for a single feature, u represents the mean of that feature, σ represents the standard deviation of that feature, g' represents the standardized data for that feature, and n represents the total number of samples. By standardizing each feature individually, we obtain the standardized input data for the training set, which is then fed into the model for training.

[0028] Furthermore, in step 5, the ARIMA calculation process is as follows:

[0029]

[0030] Where m t Let p and q be the values ​​at time t, C be a constant, p and q be the orders, β and α be the corresponding coefficients, and ε be the value at time t. t Let be the error at time t.

[0031] When searching for parameters p and q, the BIC criterion is calculated as follows:

[0032] BIG = Kln(n) - 2ln(L)

[0033] Where n is the number of samples, L is the likelihood function, and K is the number of model parameters.

[0034] Furthermore, in step 6, the calculation process of the objective function of XGB and the model of GRU is as follows:

[0035] a1,XGB is an additive model, and the output of our t-th iteration is as follows:

[0036]

[0037] in, Let i be the predicted value of sample i after the t-th iteration. For the prediction results of the first t-1 trees, f t (.) represents the model of the t-th tree;

[0038] The prediction accuracy of a model is determined by both its bias and variance. The loss function represents the model's bias, and to minimize variance, a regularization term needs to be added to the objective function to prevent overfitting. Therefore, the objective function consists of the model's loss function and a regularization term to suppress model complexity. The objective function is defined as follows:

[0039]

[0040] Where l is the loss function, Ω is the model complexity function, n represents the sample size, and y i For the true value, This is a predicted value;

[0041] Since the structure of the first t-1 trees is already determined, the sum of the complexities of the first t-1 trees can be represented by a constant during the t-th iteration. This constant term can be ignored when calculating the minimization. The objective function of the model for the t-th iteration is as follows:

[0042]

[0043] Expanding the loss function using the second-order Taylor formula yields:

[0044]

[0045] Where g i h is the first derivative of the loss function. i This is the second derivative of the loss function. Note that the derivative here is... Differentiate;

[0046] The regularization term of the objective function is determined by the normal form of the vector composed of the number of leaf nodes in all generated decision trees and the weights of all nodes. The regularization term is expanded as follows:

[0047]

[0048] Where T is the number of leaf nodes, ωXGB 2 Let L2 be the L2 norm of the leaf node weight vector;

[0049] Since the first term of the second-order Taylor expansion of the loss function is a constant, it will not affect the optimization of the objective function. Therefore, the approximate value of the objective function is:

[0050]

[0051] The second line of the above equation calculates the loss function for each sample after traversing all samples. However, the samples will eventually fall on the leaf nodes, so we can also traverse all leaf nodes, obtain the sample set on the leaf nodes, and finally calculate the loss function for each leaf node. To minimize the objective function, we can minimize the objective function of each leaf node. Using the maximum and minimum value formula of a quadratic function, we can obtain the weight of each leaf node by taking the first derivative of the objective function and setting it equal to 0. Then the third line of the above equation can be simplified to the maximum and minimum value form of a quadratic function, which simplifies the calculation.

[0052] a2, GRU is a simplified model of LSTM, memory h t By integrating long-term and short-term memory, it can effectively solve nonlinear time series problems. The model calculation process is as follows:

[0053]

[0054]

[0055]

[0056]

[0057] in For the sigmoid function, Let x be the candidate state at time t. t This is a sample at time t. t The update gate at time t controls the amount of data that can be retained up to the current time from previously memorized information. When the update gate z... t The closer the output is to 1, the closer it is to... When updating the gate z t When the value is close to 0, only the old state is retained, which comes from x. t The information was largely ignored. t The gate is reset at time t to control how much past information is forgotten. When the gate is reset, r... t When the value is close to 1, the above model is a basic RNN recurrent neural network; when the gate r is reset... t When approaching 0, It is x t The result of the multilayer perceptron as input.

[0058] Furthermore, the specific calculation method for step 7 is as follows:

[0059] b1, The update steps of the weight-adjusted particle swarm optimization algorithm are as follows:

[0060]

[0061]

[0062]

[0063] in, Let be the velocity variable of the i-th particle in the d-th iteration. Let c1 and c2 be the position variables of the i-th particle in the d-th iteration, c1 and c2 be the individual learning factors and social learning factors, and γ1 and γ2 be random numbers between zero and one. gbest is the optimal solution for the i-th particle after d-1 iterations. d-1 For the optimal solution after d-1 iterations for all particles, w start w represents the initial weight of the velocity term during the optimization process. end Here, d represents the final weight value of the velocity term during the optimization process, d represents the current iteration number, and D represents the total number of iterations.

[0064] b2, the weighted particle swarm optimization loss function is calculated as follows:

[0065]

[0066] in, v represents the prediction after running two combined models on the i-th sample. zX-i v zY-i Let represent the true component of the i-th sample after standardization, and n be the total number of samples;

[0067] b3, the calculation formula for the combined model is as follows:

[0068]

[0069]

[0070] in, The output of the combined model is ytrain, which is the fitted data obtained during the training process, and ω is the weight of the corresponding model. Therefore, the subscripts in the formula all correspond to the model.

[0071] Furthermore, in step 8, the denormalization and Gaussian process are calculated as follows:

[0072] c1, destandardized calculation method:

[0073] g=σ*g′+u

[0074] Where g' is the output of the combined model, g is the output after inverse standardization, and u and σ are the mean and variance of the predicted labels saved in step 3, respectively;

[0075] c2, the Gaussian process is calculated as follows:

[0076] The prior of the Gaussian process (the number (x,y) within the time window) is expressed as f(x) ~ N(u f ,K ff If we now observe some data (x) * ,y * ), that is, the input features of the current time model, and assuming y * f(x) follows a joint Gaussian distribution, i.e.

[0077]

[0078] Among them, K ff =k(x,x), K fy =k(x,x * ), k yy =k(x * ,x * ), Then there is

[0079]

[0080] mean The output predicted by the Gaussian process is a linear function of y, taking into account the prior mean, where y is the residual of the time window data. Based on the variance... The confidence interval for the predicted residuals can be calculated, let the mean and standard deviation of the Gaussian process prediction be u, respectively. g and σ g Then the 95% confidence interval for predicting the residuals is [u g -1.96*σ g u g +1.96*σ g ].

[0081] Furthermore, in step 10, the calculation method for the combined model prediction output is as follows:

[0082]

[0083]

[0084] Where o is the inverse standardization function.

[0085] Furthermore, in step 12, the vector composition calculation method is as follows:

[0086]

[0087]

[0088] in, The wind speeds in the X and Y directions are predicted by two independent models, respectively. For the final predicted wind speed, This is the final predicted wind direction.

[0089] Furthermore, it also includes step 13: using the predicted wind speed and direction obtained in step 12 and the actual wind speed and direction in the test set, calculate a new index RMSE to eliminate the influence of directional periodicity. θ and MAPE θ MAPE' is a correction indicator to avoid division by zero, and traditional indicators MAPE and RMSE are also mentioned. θ MAPE θ The calculation methods for MAPE, MAPE, and RMSE indices are as follows, where RMSE... θ and MAPE θ Evaluation based solely on the angle:

[0090]

[0091]

[0092] Where, Δθ i Let θ be the angle deviation during the prediction process of the i-th sample. i The true angle of the i-th sample;

[0093]

[0094]

[0095]

[0096] Among them, y i Let i be the true value of the i-th sample. Let be the predicted value of the i-th sample. When evaluating the model performance, both MAPE and MAPE' can be considered simultaneously. When the difference between the two is small, more attention should be paid to MAPE; when the difference between the two is very large, more attention should be paid to MAPE'.

[0097] Compared with the prior art, the present invention has the following innovative advantages and significant effects:

[0098] 1) By performing feature preprocessing in a standardized manner, the scaling scale of features in the training set and the test set is consistent, and the standardized features are used for prediction to improve the model's wind speed prediction adaptability.

[0099] 2) By processing the wind speed and direction obtained at each moment in a "vector" manner, it is convenient to perform operations such as averaging, calculating differences, and synthesizing the wind speed. There is no need to consider the impact of angle jump on the model. The two models are basically independent, and the errors do not affect each other. Moreover, the processing method is more in line with the physical meaning.

[0100] 3) The time series decomposition and model selection fully consider the problem of model overfitting: noise terms with minimal impact are ignored when decomposing the time series; in order to avoid the problem of model overfitting, the differential autoregressive moving average model adopts the BIC criterion for parameter determination, the extreme gradient boosting fully considers the model complexity, support vector regression itself has the lowest structural risk, and the gated recurrent unit can prevent overfitting by reducing the number of training rounds, adding L1 regularization terms, and adopting the "early stopping" strategy; the independent model has a strong ability to prevent overfitting.

[0101] 4) The differential autoregressive moving average model is used to predict periodic data. Since periodic data is generally very stable, the fitting accuracy is high, and it is less likely to experience the problem of accuracy decline when the differential autoregressive moving average model predicts more steps. The parameters of the differential autoregressive moving average model are determined by using ACF and PACF plots to determine the approximate range of parameters and using the BIC criterion for grid search. This reduces the time for parameter acquisition and prevents the model from becoming too complex.

[0102] 5) The combined model is used to predict trend data. The error of the predicted data is smaller than that of the data predicted by the individual model. The error of the base model is reduced. Furthermore, the particle swarm algorithm with linear weight correction is used to find the weights of the combined model, which speeds up the convergence speed in the weight finding process.

[0103] 6) The residual term is predicted using a Gaussian process based on a sliding time window. As the time window moves over time, the parameters of the model also change, continuously fitting the changes in the wind speed sequence. Especially when predicting wind speed sequences that the model has not learned, it can greatly compensate for the model's predictive ability and has strong self-learning and generalization capabilities.

[0104] 7) Use RMSE θ and MAPE θ The evaluation of angle takes into account the physical meaning of direction, which more realistically reflects the predictive ability of the model. The evaluation using MAPE' avoids the phenomenon of dividing the denominator by zero in the MAPE index. At the same time, MAPE' is small in difference from MAPE, and MAPE' can more accurately reflect the predictive ability of the model. Attached Figure Description

[0105] Figure 1This is a flowchart of the wind speed vector prediction model based on Gaussian process and combined model of the present invention, which illustrates the data processing process and implementation method of the X and Y independent models and the entire model;

[0106] Figure 2 This is a flowchart of the combined model of the present invention, which shows the base model and data processing method included in the combined model;

[0107] Figure 3 This is a schematic diagram of the Limiting Gradient Boosting (XGB) method of the present invention, illustrating the training process of XGB;

[0108] Figure 4 This is a schematic diagram of the gated loop unit (GRU) of the present invention, showing the specific implementation details of the GRU unit;

[0109] Figure 5 This is a schematic diagram of the Gaussian process based on time window sliding of the present invention, illustrating how the model performs residual prediction;

[0110] Figure 6 This is a graph showing the prediction results of the X and Y independent models of the present invention, illustrating the visualization results of the predictions under the independent models;

[0111] Figure 7 This is a diagram showing the wind speed and direction prediction results of the present invention, illustrating the visualization results of the entire model prediction after vector synthesis;

[0112] Figure 8 These are the ACF and PACF diagrams referenced when determining the parameters of this invention, illustrating the correlation between previous time nodes and the current time node, and used to select the order range. Detailed Implementation

[0113] The specific implementation method and working principle of the present invention will be described in detail below with reference to the accompanying drawings:

[0114] Example

[0115] This embodiment performs real-time short-term forecasting for a specific wind turbine in a wind farm. Data from the SCADA system of this wind turbine from 00:00:00 to 23:59:59 on April 11, 2022, is selected, with a data sampling interval of 1 second. The specific implementation steps are as follows:

[0116] Step 1: Obtain the operating data recorded in the SCADA system of the wind turbine for one day, which is a total of 86,400 samples. Select the first 45,060 samples as the training set, the last 5,000 samples of the training set as the validation set, and the last 30,000 samples of the validation set as the test set.

[0117] Step 2: Perform vector decomposition and use the wind speed and direction data in the system to construct new features: X-direction velocity and Y-direction velocity.

[0118] Step 3: Calculate the average wind speed in the 60 seconds before and the average wind speed in the 10 seconds before each time point, and calculate the average wind speed in the 60 seconds after each time point as the prediction label. Standardize the features (30-second average wind error, 3-second average wind error, cabin position), the average wind speed in the 60 seconds before, the average wind speed in the 10 seconds before, and the average wind speed in the 60 seconds after the training set obtained by the SCADA system, and save the mean and variance of the standardization process. Then, send the training set data into the model for training.

[0119] Step 4: Perform singular spectrum analysis (SSA) on the predicted labels. Select 20 decomposition components, select the first component as the trend term, select the sum of the second to fifteenth components as the period term, and ignore the remaining components as noise.

[0120] Step 5: Fit the periodic data using the ARIMA model. The ACF and PACF plots are shown below. Figure 8 As shown, when selecting parameters, the p-parameter range for the X model is [5,9], and the q-parameter range is [15,20]. The parameter selection range for the Y model is the same as that for the X model. The final determined p-value for the X model is 7, and the q-value is 17; the p-value for the Y model is 7, and the q-value is 17, and the fitted periodic term data are obtained.

[0121] Step 6: Use the training set data to train three models: Extreme Gradient Boosting (XGB), Gated Recurrent Unit (GRU), and Support Vector Regression (SVR). Each independent model will generate three trend prediction data.

[0122] Step 7: Use a weighted particle swarm optimization algorithm to find the weights for the combined model, while ensuring the weights ω of the differential autoregressive moving average model are maintained. ARIMA The weight ω of the limiting gradient boosting model is 1. XGB The weights ω of the gated recurrent unit model GRU The weights ω of the support vector regression model SVR Satisfy ω XGB +ω GRU +ω SVR =1. During the model weight search process, the magnitude of the difference between the predicted vector and the fitted vector is used as the loss function. Ultimately, the XGB, GRU, and SVR model weights for the X-independent model are 0.564, 0.616, and -0.180, respectively; the XGB, GRU, and SVR model weights for the Y-independent model are 0.876, -0.004, and 0.128, respectively. The combined model fitting data for the two independent models is also output.

[0123] Step 8: Destandardize the data fitted by the combined model and calculate the residuals between the model and the predicted labels on the training set. Construct a time window 60 seconds before each time point, fit the data within the time window using a Gaussian process, and directly determine the model parameters from the training set. Select the mean absolute error (MAE) of the residual fitting as the loss function, and search for hyperparameters using a grid method. The final time window length for the X-independent model is 5, and the Gaussian kernel parameter is 0.3; the time window length for the Y-independent model is 3, and the Gaussian kernel parameter is 0.1.

[0124] Step 9: Standardize the test set data using the mean and variance from the training set standardization process, and then feed it into the model for prediction.

[0125] Step 10: Using the features of the test set, calculate the predicted data of one periodic term and three trend terms for each independent model on the test set, and use the weights of the combined model to perform a weighted average. The result of the weighted average is then de-standardized and used as the prediction output of the combined model.

[0126] Step 11: Calculate the residual for each sample in the test set using a sliding Gaussian process within the time window, and add the residual to the combined model's prediction output as the final output for each independent model. The final prediction results for the X and Y independent models are as follows: Figure 6 As shown.

[0127] Step 12: Convert the prediction results of the X and Y independent models into actual wind speed and direction. The predicted wind speed and direction results are as follows: Figure 7 As shown.

[0128] Step 13, calculate the RMSE of wind speed and direction. θ MAPE θ The model performance was evaluated using the MAPE, MAPE', and RMSE indices. The calculation results are shown in Table 1. It can be seen that RMSE... θ and MAPE θ The result is smaller than the result of directly using MAPE and RMSE indicators for directional evaluation, because RMSE θ and MAPE θ The angle jump problem has been solved, which is more in line with the physical meaning. At the same time, it can be seen that MAPE' is much smaller than MAPE, thus avoiding the phenomenon of dividing the denominator by zero.

[0129] Indicator Wind speed magnitude Wind direction RMSE θ ]]> / 22.3094 MAPE θ ]]> / 17.1972 MAPE' 18.6830 / MAPE 9411066 154.3923 RMSE 1.0856 54.7810

[0130] Table 1

[0131] This invention presents a wind speed vector prediction method based on Gaussian processes and combined models. It mainly includes feature construction and preprocessing, vector decomposition, synthesis, averaging and difference calculation, time series decomposition, particle swarm optimization with weight correction for combined model prediction, Gaussian process residual prediction with time window sliding, and index modification. Experimental results show that this invention has high prediction accuracy and strong generalization ability.

[0132] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A wind speed vector prediction method based on Gaussian processes and a combined model, characterized in that, Includes the following steps: Step 1: Select a wind turbine and obtain N operational data records from the SCADA system of the wind turbine, and divide them into training set, validation set and test set; Step 2: Perform vector decomposition to convert the wind speed and direction obtained at each moment in the SCADA system into wind speed in the X direction and wind speed in the Y direction, and then construct independent X and Y models to predict the wind speed in these two directions respectively. Step 3: Construct features. Obtain the wind speed sequence after Step 2. Calculate the average wind speed before each time point and use it as a new feature. Calculate the average wind speed after each time point and use it as a prediction label. Standardize the new features, prediction labels and other features obtained by SCADA on the training set. At the same time, save the mean and variance of the standardized features. Then input them into the X and Y independent models for training. Step 4, training phase: Use singular spectrum analysis to decompose the predicted labels of the training set into different component sequences, select the first item as the trend item, select the second item and the remaining items to add together as the period item, and ignore the rest as noise. Step 5: Train the periodic term model. Select the differential autoregressive moving average model for fitting. When selecting the number of autoregressive terms p and the number of moving average terms q, first use the ACF and PACF plots to select the range of p and q, with the selection criterion being the neighborhood that produces a truncation or tailing effect. Then, according to the BIC criterion, perform grid optimization to determine the specific values ​​of p and q. After the parameters are determined, input the training set data to obtain the fitted periodic term data ytrain of the X and Y independent model. arima-X and ytrain arima-Y ; Step 6: Train the trend term model. Using the standardized input features from Step 3, select three models—limited gradient boosting, gated recurrent unit, and support vector regression—for fitting. Use the validation set to determine some hyperparameters. After determining the hyperparameters, input the training set data to obtain the fitted trend term ytrain for the X-independent model. XGB-X ytrain GRU-X ytrain SVR-X The trend term ytrain of the Y-independent model XGB-Y ytrain GRU-Y ytrain SVR-Y ; Step 7: Determine the weights of the combined model. Using the fitted data obtained in Steps 5 and 6, employ a weighted particle swarm optimization algorithm on the validation set to find the model weights, while ensuring the weights ω of the differential autoregressive moving average model are maintained. ARIMA The weight ω of the limiting gradient boosting model is 1. XGB The weights ω of the gated recurrent unit model GRU The weights ω of the support vector regression model SVR Satisfy ω XGB +ω GRU +ω SVR =1. In the model weight search process, the magnitude of the difference between the predicted vector and the fitted vector is used as the loss function, and the weight corresponding to the minimum difference is selected. By combining the model calculation formulas, the fitted data ytrain of the independent X and Y models is obtained. X ytrain Y ; Step 8: Train the Gaussian process model with a sliding time window. Using the mean and variance from the standardization process in Step 2, destandardize the fitted data obtained in Step 7 and calculate the residual between the destandardized data and the predicted label data in Step 3. For the residual data, use the input features standardized in Step 3 to construct a time window before each time node, fit the Gaussian process with the data within the time window, and directly determine the model parameters in the training set. Select the mean absolute error (MAE) of the residual fitting as the loss function and search for hyperparameters using the grid method. Step 9, Testing Phase: Using the mean and variance from the standardization process in Step 3, standardize the input feature data of the test set and input it into the X and Y independent models for prediction. Step 10: Using the periodic term model trained in Step 5, predict the periodic term data of the independent models for the test set X and Y. Using the trend term model trained in step 6, predict the three trend term data of the independent model on the test set X. and the three trend terms data of the Y-independent model Using the combined model weights determined in step 7, the combined model data is obtained through the combined model calculation formula. The combined model data is then de-standardized using the mean and variance from the standardization process in step 3 to obtain the combined model prediction output. Step 11: Using the Gaussian process model with a sliding time window trained in Step 8, predict the residual data e of the independent models of X and Y in the test set. X e Y The predicted wind speed in the X direction is obtained by adding the combined model prediction output obtained in step 10 to the residual data. Wind speed forecast in the Y direction Right now Step 12: Using the predicted X and Y values ​​calculated in Step 11, the predicted wind speed is obtained using a vector synthesis method. and wind direction 2. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 2, the vector decomposition calculation method is as follows: v X =vcosθ v Y =vsinθ Where v is the wind speed obtained by the SCADA system at each moment, and θ is the wind direction obtained by the SCADA system at each moment. X v Y These represent the wind speed in the X direction and the wind speed in the Y direction, respectively.

3. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 3, the feature standardization calculation method is as follows: Where g is the data for a feature, u is the mean of the feature, σ is the standard deviation of the feature, g′ is the standardized data for the feature, and n is the total number of samples. By standardizing each feature, the standardized input data of the training set can be obtained, and the standardized data can be fed into the model for training.

4. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 5, the calculation process of the differential autoregressive moving average model is as follows: Where, m t Let p and q be the values ​​at time t, C be a constant, p and q be the orders, β and α be the corresponding coefficients, and ε be the value at time t. t Let be the error at time t. When searching for parameters p and q, the BIC criterion is calculated as follows: BIC = Kln(n) - 2ln(L) Where n is the number of samples, L is the likelihood function, and K is the number of model parameters.

5. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 6, the calculation process of the objective function of the limiting gradient boosting and the model of the gated recurrent unit is as follows: a1, the limiting gradient boosting model is additive, and the output of the t-th iteration is as follows: in, Let i be the predicted value of sample i after the t-th iteration. f represents the prediction results of the first t-1 trees for sample i. t (.) represents the model of the t-th tree, x i Let be the i-th sample; the objective function consists of the model's loss function and a regularization term to suppress model complexity. The objective function is defined as follows: Where l is the loss function, Ω is the model complexity function, n represents the sample size, and y i For the true value, This is a predicted value; The objective function of the model in the t-th iteration is as follows: Expanding the loss function using the second-order Taylor formula yields: Among them, g i h is the first derivative of the loss function. i The second derivative of the loss function is taken here as a function of... Differentiate; The regularization term of the objective function is determined by the normal form of the vector composed of the number of leaf nodes in all generated decision trees and the weights of all nodes. The regularization term is expanded as follows: Where T is the number of leaf nodes, ωXGB 2 Let L2 be the L2 norm of the leaf node weight vector; Further, the approximate value of the objective function is obtained as follows: The second line of the above formula is the loss function for each sample after iterating through all the samples; a2, The calculation process of the gated cyclic unit model is as follows: in For the sigmoid function, Let x be the candidate state at time t. t For the sample at time t, z t The update gate at time t controls the amount of data that can be retained up to the current time from previously memorized information. When the update gate z... t The closer the output is to 1, the closer it is to... When updating the gate z t When the value is close to 0, only the old state is retained, which comes from x. t The information was ignored; t The gate is reset at time t to control how much past information is forgotten. When the gate is reset, r... t When the value is close to 1, the above model is a basic RNN recurrent neural network; when the gate r is reset... t When approaching 0, It is x t The result of the multilayer perceptron as input.

6. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, The specific calculation method for step 7 is as follows: b1, The update steps of the weight-adjusted particle swarm optimization algorithm are as follows: in, Let be the velocity variable of the i-th particle in the d-th iteration. Let c1 and c2 be the position variables of the i-th particle in the d-th iteration, c1 and c2 be the individual learning factors and social learning factors, and γ1 and γ2 be random numbers between zero and one. gbest is the optimal solution for the i-th particle after d-1 iterations. d-1 For the optimal solution after d-1 iterations for all particles, w start w represents the initial weight of the velocity term during the optimization process. end Here, d represents the final weight value of the velocity term during the optimization process, d represents the current iteration number, and D represents the total number of iterations. b2, the weighted particle swarm optimization loss function is calculated as follows: in, v represents the prediction after running two combined models on the i-th sample. zX-i v zY-i Let represent the true component of the i-th sample after standardization, and n be the total number of samples; b3, the calculation formula for the combined model is as follows: in, The output of the combined model is ytrain, which is the fitted data obtained during training, and ω is the weight of the corresponding model. The subscripts in the formula all correspond to the model.

7. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 8, the inverse standardization and Gaussian process are calculated as follows: c1, destandardized calculation method: g=σ*g′+u Where g′ is the output of the combined model, g is the output after inverse standardization, and u and σ are the mean and variance of the predicted labels saved in step 3, respectively; c2, the Gaussian process is calculated as follows: The prior of the Gaussian process is expressed as f(x) ~ N(u f K ff If some data (x) are observed now * y * ), that is, the input features of the current time model, and assuming y * f(x) follows a joint Gaussian distribution, i.e. Among them, K ff =k(x, x), K fy =k(x, x) * ), K yy =k(x * x * ), Then there is mean The output predicted by the Gaussian process is a linear function of y, where y is the residual of the time window data, calculated based on the variance. The confidence interval for the predicted residuals can be calculated, let the mean and standard deviation of the Gaussian process prediction be u, respectively. g and σ g Then the 95% confidence interval for predicting the residuals is [u g -1.96*σ g u g +1.96*σ g ].

8. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 10, the calculation method for the prediction output of the combined model is as follows: Where o is the inverse standardization function.

9. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, In step 12, the vector composition calculation method is as follows:

10. The wind speed vector prediction method based on Gaussian process and combined model according to claim 1, characterized in that, The process also includes step 13: using the predicted wind speed and direction obtained in step 12 and the actual wind speed and direction in the test set, a new index, RMSE, is calculated to eliminate the influence of directional periodicity. θ and MAPE θ MAPE′ is a correction indicator to avoid division by zero, and traditional indicators MAPE and RMSE are also mentioned. θ MAPE θ The calculation methods for MAPE, MAPE, and RMSE indices are as follows, where RMSE... θ and MAPE θ This is an evaluation based solely on the angle: Where, Δθ i Let θ be the angle deviation during the prediction process of the i-th sample. i The true angle of the i-th sample; Where y i Let i be the true value of the i-th sample. Let be the predicted value for the i-th sample.

Images