An ultra-short-term photovoltaic output prediction method based on a bimodal weather classification

Abstract

The application discloses a kind of based on dual-mode weather classification ultra-short photovoltaic output prediction method, belong to new energy power system prediction technical field, including: obtaining meteorological data and photovoltaic output data is preprocessed;According to the data after preprocessing, construct dual-mode feature matrix;Adopt dual-mode joint distance function self-adapting distribution time-frequency domain feature weight coefficient, according to the weight coefficient after distribution, by reducing the dual-mode feature matrix of processing, then using spectral clustering realizes weather classification;According to the classification result, different weather scene uses corresponding signal decomposition method;According to signal decomposition result, construct hybrid prediction model, and pass through improved optimization algorithm optimization model hyperparameter;According to the hybrid prediction model after optimization, input fusion model is predicted, obtains the final prediction result.The application effectively solves the problem of inaccurate prediction caused by weather sensitivity, modal aliasing and parameter solidification, significantly improves the prediction accuracy under rainy weather scene.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of new energy power system prediction, and particularly relates to a super-short-term photovoltaic output prediction method based on double-mode weather classification. BACKGROUND

[0002] With the vigorous development of global economy, the demand for energy is increasing, and as an important clean energy, the accurate prediction of solar power output becomes crucial. At present, artificial intelligence methods have provided a new and potential solution for photovoltaic output prediction, but the existing hybrid model still has some limitations, resulting in insufficient prediction technology to meet the strict accuracy requirements.

[0003] The existing problems and shortcomings mainly include: (1) single weather classification: traditional methods only rely on time sequence patterns, ignoring frequency energy distribution characteristics, resulting in inaccurate rain weather identification. (2) fixed decomposition strategy: fixed methods cannot adapt to the fluctuation characteristics of different weather types, which may produce redundant modes or information loss. (3) parameter optimization limitations: traditional optimization algorithms are prone to local optimization, and the global collaborative optimization of classification-decomposition-prediction parameters has not been realized.

[0004] In view of the above problems of the prior art, it is urgent to propose a super-short-term photovoltaic output prediction method based on double-mode weather classification and dynamic signal decomposition. SUMMARY

[0005] To solve the above technical problems, the application provides a super-short-term photovoltaic output prediction method based on double-mode weather classification, comprising the following steps:

[0006] Obtain meteorological data and photovoltaic output data of the research area and perform preprocessing;

[0007] According to the preprocessed data, a double-mode feature matrix is constructed by combining dynamic time warping distance and frequency energy feature matrix;

[0008] Adaptive allocation of time-frequency domain feature weight coefficients is performed using a double-mode joint distance function, and the double-mode feature matrix is processed by dimension reduction according to the allocated weight coefficients;

[0009] According to the double-mode feature matrix after dimension reduction, spectral clustering is used to realize the classification of four types of weather, namely sunny, cloudy, rainy and special weather;

[0010] According to the classification results, corresponding signal decomposition methods are used for different weather scenarios;

[0011] According to the signal decomposition results, a hybrid prediction model is constructed, and the model hyperparameters are optimized in parallel by an improved constant star oscillation optimization algorithm;

[0012] According to the optimized hybrid prediction model, input the fusion model for prediction to obtain the final prediction result.

[0013] Optionally, the process of acquiring and preprocessing the meteorological data and photovoltaic output data of the research area comprises:

[0014] The meteorological data and photovoltaic output data of the research area are subjected to outlier cleaning, interpolation of missing data, data standardization, and key meteorological factor screening.

[0015] Optionally, the process of implementing the four types of weather classification, i.e., sunny, cloudy, rainy, and special weather, comprises:

[0016] According to the input data, time window slices are divided, and complementary features are extracted from the time domain and the frequency domain;

[0017] According to the time domain and frequency domain features, a bimodal joint distance function is defined, the bimodal feature matrix is reduced in dimension when the weight coefficient is a preset weight coefficient, the reduced bimodal feature matrix is subjected to spectral clustering, and a temporary clustering result is obtained;

[0018] According to the temporary clustering result, an optimization objective function is defined, a multi-objective optimization algorithm is used to obtain the optimal weight coefficient of the bimodal feature matrix, the optimal weight coefficient is used to repeat the process of dimension reduction and spectral clustering, a spectral clustering result is obtained, and a weather type label is dynamically mapped.

[0019] Optionally, the process of using a corresponding signal decomposition method for different weather scenarios according to the classification result comprises:

[0020] If it is a sunny scenario, CEEMDAN decomposition is used and sample entropy is calculated, and only low-frequency components with sample entropy less than 0.2 are retained;

[0021] If it is a cloudy scenario, CEEMDAN decomposition is used and sample entropy is calculated, and only medium-frequency components with sample entropy less than 0.4 are retained;

[0022] If it is a rainy scenario, CEEMDAN decomposition is first used, and then the top 5 components in terms of sample entropy value are subjected to secondary decomposition, and finally components with sample entropy less than 0.5 are retained;

[0023] If it is a special weather scenario, decomposition is used and components with sample entropy less than 0.6 are retained.

[0024] Optionally, the hybrid prediction model comprises a Transformer module and a bidirectional gated recurrent module.

[0025] The Transformer module is used to dynamically calculate the correlation weight between different input features through the self-attention mechanism of the Transformer encoder, and adaptively identify and focus on key influencing factors, including photovoltaic output, irradiance and temperature; the bidirectional gated recurrent module is used to dynamically capture the long-term dependence relationship between nonlinear time series by using the front and rear information in the sequence.

[0026] Optionally, the process of optimizing the model hyperparameters in parallel through the improved stellar oscillation optimization algorithm comprises:

[0027] According to the initial parameters, an initial oscillator population is generated, and the position of each oscillator is disturbed by using a Logistic mapping;

[0028] According to the disturbed oscillator position, each generated oscillator individual is updated by stellar seismic oscillation, and the oscillator is subjected to Levy flight disturbance for mutation;

[0029] According to the oscillator after the Levy flight disturbance, an oscillatory motion affected by a sine and cosine function is added to generate a new position;

[0030] According to the new position, the fitness is evaluated, the mirror reverse mapping of the oscillator with poor fitness is performed to generate a mirror solution and keep the better position, the particle position with the optimal fitness is kept, and the above steps are repeatedly executed until the maximum iteration number or the convergence threshold is reached, and the optimized model hyperparameters are obtained.

[0031] Optionally, the process of inputting the optimized hybrid prediction model into the fusion model for prediction to obtain the final prediction result comprises:

[0032] According to the optimized hybrid prediction model, a multi-source feature matrix is input into the finally built fusion model, and the multi-source feature matrix includes encoded and decomposed components and meteorological factors;

[0033] According to the characteristics of the model structure and the photovoltaic output data, a key hyperparameter search space is defined;

[0034] According to the defined hyperparameter search space, the hybrid prediction model is trained with an initial parameter set;

[0035] According to the training result, the above steps are repeatedly executed until a specified iteration number is reached or a performance target is reached, and the parameter configuration with the minimum error on the validation set is selected as the setting of the hybrid prediction model;

[0036] According to the final hybrid prediction model, the final photovoltaic output prediction result is obtained.

[0037] The application also provides a computer device, comprising a memory, a processor and a computer program stored in the memory, wherein the processor executes the computer program to realize the steps of the method.

[0038] The application also provides a computer readable storage medium, wherein a computer program is stored on the computer readable storage medium, and the computer program is executed by a processor to realize the steps of the method.

[0039] The application also provides a computer program product, comprising a computer program, wherein the computer program is executed by a processor to realize the steps of the method.

[0040] Compared with the prior art, the application has the following advantages and technical effects:

[0041] The prediction method provided by the application significantly improves the weather typing accuracy and physical interpretability. The bimodal feature matrix is constructed by the dynamic time warping (DTW) distance matrix and the power spectral density (PSD) frequency energy feature, so that the weather misjudgment problem caused by the fact that the traditional method only depends on the time sequence form is solved. Meanwhile, the NSGA-II multi-objective optimization algorithm is combined to adaptively distribute the time-frequency domain weight coefficients, so that the weather recognition accuracy is improved, and the physical interpretability of the typing process is enhanced.

[0042] The dynamic decomposition strategy of the application effectively suppresses the modal aliasing and noise interference. The signal decomposition method is adaptively selected based on the weather typing label, the information loss problem of the fixed parameter decomposition is overcome, and the feature purity is significantly improved.

[0043] The hybrid prediction model of the application cooperates with the improved optimization algorithm to improve the generalization ability of the model, and more accurate prediction results are obtained. BRIEF DESCRIPTION OF DRAWINGS

[0044] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The embodiments of this application and their description are used to explain the application and are not intended to limit the application. In the drawings:

[0045] Figure 1 The figure is a technical flowchart of the embodiment of the application;

[0046] Figure 2 The figure is a flowchart of the DTW-PSD bimodal weather typing of the embodiment of the application;

[0047] Figure 3 The figure is a flowchart of the signal decomposition method adaptively selected based on the dynamic weather typing label of the embodiment of the application;

[0048] Figure 4 The figure is a network structure diagram of the Transformer-BiGRU prediction model of the embodiment of the application;

[0049] Figure 5 The improved ISOO algorithm flowchart of the embodiment of the application. DETAILED DESCRIPTION

[0050] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0051] It should be noted that the steps shown in the flowchart of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and although the logical order is shown in the flowchart, in some cases, the steps shown or described herein can be executed in an order different from that shown herein.

[0052] Embodiment one

[0053] The existing technology weather typing precision is insufficient and the physical explainability is missing: the traditional method mainly relies on single time sequence shape analysis, and ignores the frequency energy distribution characteristics, resulting in high error in the recognition of rain weather scene. This defect makes the prediction model unable to adapt to the output fluctuation characteristics of different weather types, directly affecting the accuracy of the subsequent decomposition and prediction links.

[0054] The signal decomposition strategy in the prior art is fixed, resulting in modal aliasing and information loss: the prior art uses a fixed parameter decomposition method (such as single CEEMDAN decomposition or VMD decomposition), which cannot dynamically adjust the decomposition strategy according to the weather type. This fixed strategy exacerbates the modal aliasing phenomenon, reduces the purity of feature extraction, and the prediction effect is poor.

[0055] The hyperparameter optimization of the prediction model in the prior art uses a traditional optimization algorithm, which is easy to fall into local optimum and has slow convergence speed, and the model prediction effect is poor.

[0056] In order to solve the problems in the prior art described above, the present embodiment aims to provide a more accurate, efficient and economical scheme for ultra-short-term photovoltaic output prediction, so as to effectively improve the effect of photovoltaic output prediction.

[0057] As shown in Figure 1 The ultra-short-term photovoltaic output prediction method based on double-mode weather typing provided by the present application mainly includes the following steps:

[0058] Step S1, data acquisition and preprocessing;

[0059] Step S2, type weather based on time domain feature DTW and frequency domain feature PSD double mode;

[0060] Step S3, self-adaptive selection of signal decomposition method based on dynamic weather typing label;

[0061] Step S4, constructing a prediction model and optimizing hyperparameters;

[0062] Step S5, predicting photovoltaic output data.

[0063] Specifically, step S1, data acquisition and preprocessing process includes:

[0064] Step S11, multi-source data acquisition and outlier cleaning, including the following steps:

[0065] Step S111, collect multi-source data: first, collect the preset time hourly meteorological data (such as total irradiance GHI, normal direct irradiance DNI, scattered irradiance DHI, ambient temperature, wind speed, relative humidity, atmospheric pressure) and photovoltaic output data (photovoltaic output sequence every 15 minutes) and photovoltaic power station characteristic data (such as photovoltaic module parameters, photovoltaic installed capacity) from the meteorological station and photovoltaic power station in the study area.

[0066] Step S112, clean the outliers of the data, follow these two principles: (1) calculate the mean X mean and the standard deviation X std of the photovoltaic output sequence, and eliminate the outliers outside the range [X mean -3X std , X mean +3X std ]; (2) if the photovoltaic output value exceeds 120% of the installed capacity or the output is >0 at night (irradiance ≤10 W / m²), it is determined as an abnormal value and is removed.

[0067] Step S12, interpolate missing values and align multi-source data, including the following steps:

[0068] Step S121, if the data is missing within 4 time sequences, it is determined as short-time missing (≤4 time sequences), and the mean value of the adjacent time of the same photovoltaic power station is used to fill; if the data is missing within 4 time sequences, it is determined as long-time missing (>4 time sequences), and the mean value of the same time in the past 3 days is used to fill.

[0069] Step S122, linearly interpolate the meteorological data (irradiance, temperature, etc.) based on the photovoltaic output timestamp to ensure uniform time resolution and strictly align the meteorological data and photovoltaic output data timestamps.

[0070] Step S13, standardizing the data and selecting key meteorological factors, including the following steps:

[0071] Step S131, only normalize the meteorological variables, using the maximum-minimum value normalization method, the specific formula is as follows:

[0072] ;

[0073] MET new is the standardized data, ranging from [0, 1]; MET raw is the original data; MET max and MET min are the maximum and minimum values of the original data, respectively.

[0074] Step S132, photovoltaic output data standardization, the photovoltaic output keeps its physical dimension, only standardization is performed, the specific formula is as follows:

[0075] ;

[0076] MET std is the standardized photovoltaic output data; P is the original photovoltaic output; X mean is the mean value of the photovoltaic output; X std is the standard deviation of the photovoltaic output sequence.

[0077] Step S133, screening of key meteorological factors, Pearson correlation coefficient analysis is adopted to calculate the linear correlation of each meteorological variable and photovoltaic output, and variables with Pearson correlation coefficient ≥ 0.3 are reserved; the calculation formula of the Pearson correlation coefficient is as follows:

[0078] ;

[0079] MET is the Pearson correlation coefficient between the meteorological variable and the photovoltaic output; MET is the meteorological variable; P is the photovoltaic output data; i is the time point; MET mean , P mean is the average value of the variable MET, P; n is the total number of data samples.

[0080] Step S2, DTW-PSD bimodal weather classification, the contents of each part are shown in Figure 2 .

[0081] Step S21, input data division and time window slicing, including:

[0082] Step S211, input historical data set: photovoltaic output sequence P i (i=1, 2,..., N, keeping the physical dimension) and meteorological variable matrix X=[GHI, DNI, temperature, wind speed, …], here the length of the meteorological variable matrix and the length of the photovoltaic output sequence keep consistent, both are N.

[0083] Step S212, time window slicing: the photovoltaic output sequence is divided into M subsequences Pd (d = 1, 2, …, M), (data is 15 minutes, then the sequence length of a day is 96, M represents the total number of days).

[0084] Step S22, bimodal feature extraction, complementary features are extracted from time domain and frequency domain, and the physical interpretability of weather classification is enhanced, and the steps are as follows:

[0085] Step S221, calculate the time domain shape feature (DTW distance matrix); calculate the shape similarity between photovoltaic output subsequences P d The calculation formula is as follows:

[0086] ;

[0087] In the formula: route is a regular path that meets the monotonicity, continuity and boundary alignment constraints; is the kth sampling point of the dth subsequence; is the lth sampling point of the dth subsequence, ; the output D DTW is the DTW distance matrix.

[0088] Step S222, calculate the frequency energy feature matrix (PSD feature matrix), as follows:

[0089] Step S222a, first, perform fast Fourier transform on each subsequence P d :

[0090] ;

[0091] In the formula: is the frequency component, with Hz as the unit; represents the complex frequency spectrum value of the dth day subsequence at the frequency ; j is the imaginary unit; i represents the ith time point.

[0092] Step S222b, according to the calculated frequency of different subsequences, classify different sequences: the frequency range of low-frequency sequences is 0.1~1 Hz, reflecting the daily periodic change (such as sunny day); the frequency range of medium-frequency sequences is 1~3 Hz, reflecting the hourly fluctuation (such as cloudy day); the frequency range of high-frequency sequences is 3~10 Hz, reflecting the minute-level mutation (such as cloud movement).

[0093] Step S222c, calculate the energy of different frequency sequences in different frequency bands:

[0094] ;

[0095] In the formula: , , These represent the energies of the d-th day's sequence at low, medium, and high frequencies, respectively. For frequency components, the unit is Hz; Represents the frequency of the subsequence on day d. The complex spectrum value at that location;

[0096] Step S222d: Normalize the eigenvectors to obtain the PSD matrix.

[0097] ;

[0098] ;

[0099] In the formula: The frequency domain feature vector representing day d has a dimension of 1×3; each element represents the energy proportion of the corresponding frequency band, ranging from 0 to 1. The final output is the PSD feature matrix F. PSD。

[0100] Step S23: Define the bimodal joint distance function and solve for the t-SNE dimensionality reduction matrix of the fusion matrix when the weight coefficient α is 0.5, including:

[0101] Step S231: Define and calculate the bimodal joint distance function. The formula is as follows:

[0102] ;

[0103] In the formula, p and q represent the indices of the weather subsequence; α is a weighting coefficient that controls the contribution in the time and frequency domains; It is the Euclidean distance, which can be used to calculate the differences in frequency domain features; This represents the fusion feature matrix.

[0104] Step S232: Using a weight coefficient α=0.5, perform t-SNE dimensionality reduction on the fusion matrix. By minimizing the KL divergence, the high-dimensional features are mapped to a 2-dimensional space, including the following steps:

[0105] Step S232a: Input the fusion feature matrix Initialize the dimensionality reduction matrix Z0, and perform iterative optimization with the objective function of minimizing the KL divergence. Calculate the weather feature similarity (HDS) in the high-dimensional space. p,q The formula is as follows:

[0106] ;

[0107] In the formula, p and q represent the indices of the weather subsequence; It is the index for traversing all data points except p; The perplexity is calculated by the weather-specific perplexity = 30, ensuring that each data point retains about 30 neighboring points, fitting the local fluctuation characteristics of the photovoltaic sequence.

[0108] Step S232b, calculating the low-dimensional space similarity LDS p,q , the formula is as follows:

[0109] ;

[0110] In the formula, is the second index traversing all point pairs, and are not equal; Z is the dimension reduction matrix.

[0111] Step S232c, calculating the divergence using the gradient formula, the formula is as follows:

[0112] ;

[0113] Step S232d, updating the dimension reduction matrix Z using the momentum acceleration principle, the formula is as follows:

[0114] ;

[0115] ;

[0116] In the formula, represents the momentum vector of the dimension reduction matrix Z t at time t, is the momentum vector of the dimension reduction matrix Z t-1 at the previous time; is the learning rate, taking 200.

[0117] Step S232e, when the KL divergence changes less than 0.0001 for 1000 iterations or 50 consecutive generations, stop the iteration, and output the final dimension reduction matrix Z.

[0118] Step S24, performing spectral clustering on the dimension reduction matrix obtained when the weight coefficient a = 0.5 to obtain a temporary clustering result, including the following steps:

[0119] Step S241, constructing a similarity matrix, the formula is as follows:

[0120] ;

[0121] In the formula, p and q represent the indexes of the weather sub-sequence; represents the width of the Gaussian kernel, controlling the decay rate of the similarity, taking 0.5, fitting the discrete characteristics of photovoltaic weather.

[0122] Step S242, symmetrically normalizing the Laplacian matrix, the formula is as follows:

[0123] ;

[0124] In the formula, L sym represents a symmetric normalized Laplacian matrix, D is a degree matrix, W is a similarity matrix, and I is a unit matrix.

[0125] Step S243, mapping the original nonlinear data to a spectral space by eigen-decomposing the Laplacian matrix, and then selecting a cluster number K, comprising:

[0126] Step S243a, solving an eigenvalue problem, and the formula is as follows:

[0127] ;

[0128] In the formula, V is an eigenvector corresponding to the first K smallest eigenvalues; represents an eigenvalue diagonal matrix.

[0129] Step S243b, selecting a cluster number K, traversing K from 2 to 10, calculating a silhouette coefficient, the formula is as follows, selecting the K with the largest silhouette coefficient as the number of weather types of clustering.

[0130] ;

[0131] In the formula: a x represents the average value of the distance between the xth data X x and other data in the same cluster; b x represents the nearest average distance obtained from any cluster other than X x . For data X i , Silhouette(K) represents the overall classification effect, and the value range is [-1, 1], and the closer to 1, the better the classification effect. The value of Silhouette(K) is calculated when K takes different values, and the value of K is obtained when Silhouette(K) is the largest.

[0132] Step S244, performing Kmeans clustering on the data on the spectral space, comprising:

[0133] Step S244a, initializing a cluster center, and randomly selecting K centers from the eigenvectors;

[0134] Step S244b, calculating the minimum distance between each eigenvector and each center eigenvector, and assigning each eigenvector to the class represented by the nearest cluster center;

[0135] Step S244c, for each cluster, recalculating the eigenvector of the cluster center;

[0136] Step S244d, when the center change threshold is less than 0.0001, output the current clustering center and the new K-class feature vector.

[0137] Step S25, define an optimization objective function, and obtain the optimal weight coefficient a when fusing the matrix by using the NSGA-II algorithm, dynamically balance the weights of the time domain and frequency domain features, and improve the classification robustness, including:

[0138] Step S251, define the optimization objective function:

[0139] Objective function 1: maximize the contour coefficient, and the objective function is as follows:

[0140] ;

[0141] In the formula, a x represents the average distance of the xth data X x to other data in the same cluster; b x represents the nearest average distance obtained by selecting any cluster other than X x ; Silhouette(K) represents the overall classification effect, and M represents the total number of days.

[0142] Objective function 2: minimize the variance of the DTW distance in the same class, and the objective function is as follows:

[0143] ;

[0144] In the formula, p and q represent the indices of the weather subsequence; represents a variance calculation function.

[0145] Step S252, obtain the optimal weight coefficient a by using the NSGA-II algorithm, and the flow of the algorithm is as shown in Figure 1 , and the specific steps are as follows:

[0146] Step S252a, randomly generate 50 groups of initial populations a, a [0, 1], which are the parent populations. Among them, the evolution generation counter is set to Gen = 1.

[0147] Step S252b, calculate the values of the two objective functions for each a, then sort the solution set according to the Pareto dominance relationship, front 1 is the Pareto optimal solution, which is not dominated by any other individual, front 2 is the solution only dominated by the front 1 individual, and so on. Sort the individuals in the same front layer according to the objective function value. Randomly select two different individuals to simulate binary crossover to generate offspring, and perform polynomial mutation on the generated offspring. Output the offspring population with a size of 50; update the evolution generation: Gen = Gen + 1 = 2.

[0148] Step S252c, population merging and loop control, merging the parent population and the child population, if a new parent population is generated, executing step S252b; otherwise, executing step S253d.

[0149] Step S252d, calculating the value of two objective functions for each individual in the merged population, marking the front level and calculating the crowding degree, screening the top 50 solutions according to the rule of ascending front level + descending crowding degree, and obtaining a new parent population.

[0150] Step S252e, selecting the optimal weight coefficient α after 100 iterations to obtain the best contribution degree of the control time domain and frequency domain.

[0151] Step S26: using the optimal weight coefficient α, repeating steps S232-S24 to obtain the results of spectral clustering and dynamically mapping the weather type label, including the following:

[0152] Step S261, marking each sub-sequence after spectral clustering as , inputting the different clusters after clustering, and calculating the intra-class DTW distance variance and the energy proportion of the frequency band , , ;

[0153] Step S262, dynamic label mapping, if and , marking the cluster as sunny; if and , marking the cluster as rainy weather; if and , marking the cluster as cloudy weather. If , marking the cluster as special weather.

[0154] Step S263, mapping the cluster label to all sub-sequences, outputting the weather type label set .

[0155] Step S3, adaptively selecting a signal decomposition method based on the dynamic weather classification label (sunny / cloudy / rainy / special weather) to achieve physically meaningful and accurate feature extraction, and the network structure of the model is shown as follows: Figure 3

[0156] Step S31, if it is a sunny scene (L d = sunny), CEEMDAN decomposition is used, and the specific steps are as follows:

[0157] Step S311, adding times of Gaussian white noise with mean value 0 to the original photovoltaic power sequence Y(t) to generate a disturbance signal set: ​

[0158] ;

[0159] where t is time; is the signal-to-noise ratio, and controls the amplitude of the added noise; is the first white noise sequence, The value range of M is [1, M].

[0160] Step S312, each is decomposed by applying the Empirical Mode Decomposition (EMD) algorithm to obtain the set average I1(t) of the first modal component (IMF1) and the residual component s1(t).

[0161] ;

[0162] ;

[0163] where t is time; is the first modal component of the sequence to be decomposed.

[0164] Step S313, noise is added to the residual component, and EMD is continued to be decomposed, and the residual is updated:

[0165] ;

[0166] ;

[0167] ;

[0168] where t is time; num represents the modal order, is the signal-to-noise ratio of the first num-1 order white noise; is the num-1 order modal component obtained by EMD decomposition; I num is the average value of the num order modal component after decomposition; s num is the num-1 order residual component.

[0169] Step S314, when the extreme points of the residual s num (t) do not exceed 2, stop decomposing, and the final signal is represented as:

[0170] ;

[0171] where t is time; num' represents that the original photovoltaic power sequence is decomposed into num' modal components, represents the residual term.

[0172] Step S315, sample entropy SE is calculated by SE, and the specific steps are as follows:

[0173] Step S315a, the modal component I num (t) is constructed as a set of k-dimensional vectors:

[0174] ;

[0175] ;

[0176] In the formula, k em is the embedding dimension, N em is the sequence length, and z represents the number of vectors in the vector set.

[0177] Step S315b, set the similarity tolerance ST, and calculate the maximum absolute distance between any two vectors Y z (k em ) and Y z’ (k em ):

[0178] ;

[0179] In the formula, is an index variable for traversing each element in the k em dimensional vector, ; if , it is considered that the two vectors are similar.

[0180] Step S315c, the number of vector pairs satisfying the similarity condition B z (k em ) and B z (k em +1) are counted, and the sample entropy is defined as:

[0181] ;

[0182] In the formula, B z (k em ) and B z (k em +1) represent the number of vector pairs when the dimension is k em and k em +1, respectively.

[0183] Step S316, only the ultra-low frequency component with SE<0.2 is retained, reflecting the daily cycle trend, and the noise is removed.

[0184] Step S32, if it is a cloudy scene (L d = cloudy), SEEMDAN decomposition and sample entropy SE calculation are adopted, and only the medium frequency component with SE<0.4 is retained.

[0185] Step S33, if it is a rainy day scene (L d = rainy day), first decompose using CEEMDAN, and then perform VMD secondary decomposition on the SE entropy values of the first 5 IMF components. The specific steps are as follows:

[0186] Step S331, assuming that the original photovoltaic power sequence Y(t) can be decomposed into J modal components, the goal is to minimize the sum of the estimated bandwidths of all modes and ensure that the sum of the components is equal to the original signal. The constrained variational problem can be expressed as:

[0187] ;

[0188] In the formula, t is the time; u β (t) represents the βth modal component, w β is the center frequency; J represents the total number of modal components, represents the impulse function; * represents convolution operation, j represents the imaginary unit, and is used to construct a complex frequency modulation kernel function.

[0189] Step S332, convert the original problem into an unconstrained optimization problem:

[0190] ;

[0191] In the formula, t is the time; γ represents a quadratic penalty factor, and the quadratic penalty term can improve the convergence of the limited weight of the penalty term; is the Lagrange operator, used to perform constraints.

[0192] Step S333, update , and alternately until the convergence condition is met:

[0193] ;

[0194] ;

[0195] In the formula, is the Fourier transform of ; is the Fourier transform of the original signal, l num is the number of iterations.

[0196] Step S334, after decomposition, L components u1(t), u2(t), , u L (t) are obtained, and the center frequencies are arranged from high to low.

[0197] Step S335, calculate the sample entropy value of each sequence, and retain the IMF components with SE <0.5.

[0198] Step S34, if it is a special weather scenario (L d = special weather), CEEMDAN decomposition is used and IMF components with SE < 0.6 are retained.

[0199] Step S35, the weather type label L d is One-Hot encoded, and the steps are as follows:

[0200] If L d = sunny, the One-Hot encoding is [1, 0, 0, 0], if L d = cloudy, the One-Hot encoding is [0, 1, 0, 0], if L d = rainy, the One-Hot encoding is [0, 0, 1, 0], and if L d = special weather, the One-Hot encoding is [0, 0, 0, 1].

[0201] Step S36, the One-Hot encoding is spliced with the IMF components obtained by decomposition and the weather factor matrix to form a three-dimensional feature tensor and output.

[0202] Step S4, a Transformer-BiGRU hybrid prediction model fused with ISOO hyperparameter optimization is constructed;

[0203] Step S41, a Transformer-BiGRU hybrid prediction model is constructed, and the model structure is as shown in Figure 4 . The specific steps are as follows:

[0204] Step S411, a Transformer module is built. The core mechanism of the Transformer model is the self-attention mechanism, which associates and pays attention to each part of the input sequence with other parts, so that the model can better capture long-distance dependencies. The self-attention mechanism is implemented by the following formula:

[0205] ;

[0206] In the formula, Q matrix is the query matrix, K matrix is the key matrix, V matrix is the value matrix, all of which are obtained by linear transformation; softmax is a normalization function, d k is the dimension of the key matrix K matrix , is a scaling factor, and T is the length of the input sequence.

[0207] Step S412: Construct a bidirectional gated recurrent (BiGRU) module. Based on the Transformer output, BiGRU models the time series data from both forward and backward directions. Given the same time series as input, it processes the context information of each position in the input sequence from both forward and backward directions, and then concatenates the hidden vectors as the final output. At each time step, the network hidden state is iterated jointly by two GRU networks in different directions, as shown in the following expression:

[0208] ;

[0209] In the formula, t represents time. , These represent the output states of the forward and reverse GRU networks at time t, respectively. , x represents the output state of the forward and reverse GRU networks at the previous time step, respectively. t w represents the input feature at the current time step. f,t w b,t These are the weights corresponding to the forward and backward hidden layers, b. t Let be the bias matrix at time t.

[0210] Step S413: Construct a Transformer-BiGRU hybrid prediction model. First, the Transformer encoder's self-attention mechanism dynamically calculates the correlation weights between different input features, adaptively identifying and focusing on key influencing factors such as photovoltaic output, irradiance, and temperature, deeply exploring the complex nonlinear relationship between meteorological characteristics and photovoltaic output. Then, the extracted global features are input into the BiGRU network, utilizing the preceding and following information in the sequence to dynamically capture long-term dependencies between nonlinear time-series sequences, enhancing the modeling capability of the data. By combining the advantages of the Transformer's multi-head attention mechanism and the BiGRU network, the constructed combined model can effectively extract deep-level information from multi-feature data, significantly improving the prediction accuracy and generalization ability of complex sequence data, thus making it more suitable for photovoltaic output prediction tasks with multiple feature inputs.

[0211] Step S42: For the Transformer-BiGRU hybrid prediction model, the improved Stellar Oscillation Optimization (ISOO) algorithm is used to search for the global optimum of the hyperparameters. The ISOO algorithm process is as follows: Figure 5 As shown. The specific steps are as follows:

[0212] Step S421: Generate the initial oscillator population, using the following formula:

[0213] ;

[0214] In the formula, represents the position of each particle in the multi-dimensional solution space, wherein represents the particle currently being updated, represents the parameter dimension currently calculated; C is the size of the candidate solution contained in the population, and Dim represents the expandable dimension of the solution space; and respectively represent the upper and lower bounds of the dimensional solution space, represents a random matrix uniformly distributed on [0, 1] C×Dim .

[0215] Step S422, the position of each oscillator is disturbed by using the Logistic mapping, breaking the uniformity of the traditional random generation, effectively improving the diversity of the population and enhancing the global search ability of the algorithm. The Logistic mapping formula is as follows:

[0216] ;

[0217] In the formula, is the position of the th oscillator at the f th iteration; is the new position of the th oscillator after disturbance at the f +1 th iteration; h is the parameter of the Logistic mapping, and the value range is (3.57, 4].

[0218] Step S423, each generated oscillator individual is subjected to star seismic oscillation update, two candidate solutions are generated and the average position is calculated, and the final updated position is the average value, and the formula is as follows:

[0219] ;

[0220] ;

[0221] ;

[0222] In the formula, t is the time; r1, r2, r3 are random numbers generated in the range of [0, 1] at each iteration of the oscillator; and respectively are two candidate solutions generated in the dimensional solution space; w(t) represents a weight decay function; S(t) represents an amplitude scaling function; represents the best solution generated when the oscillator is iterated; represents the latest oscillator at this oscillation update.

[0223] Step S424, Lévy flight disturbance is applied to the oscillator for mutation, and the formula of Lévy flight disturbance is as follows:

[0224] ;

[0225] In the formula, r6 is a random number in the range [0, 1] generated by the oscillator in each iteration; This is the probability threshold that can trigger Levi's flight, set to 0.3; if r6 is higher than... If the value is small, then execute the Levi flight disturbance; This represents the globally optimal position up to the f-th iteration; It is a dynamic step size factor that decays with the number of iterations. r4 and r5 are independent random variables that follow a standard normal distribution. This is the Levi index, set to 1.5.

[0226] Step S425: For the oscillator after the Lévy flight disturbance, add oscillatory motion influenced by sine and cosine functions to generate a new position, as shown in the following formula:

[0227] ;

[0228] In the formula, r7 is a random number in the range [0, 1] generated by the oscillator in each iteration; Representing the The first oscillator The position of the solution space in the (f+1)th iteration; , , These are randomly selected reference oscillator positions.

[0229] Step S426: Evaluate fitness. For oscillators with poor fitness, perform mirror inversion mapping to generate a mirror solution and retain the better position, as shown in the following formula:

[0230] ;

[0231] In the formula, Indicates the first The oscillator in the first... The mirror position generated by the dimension Indicates the first The oscillator is at the f-th iteration. The position of the dimension, Fitness represents the fitness value of the oscillator.

[0232] Step S427: Retain the position of the particle with the best fitness, and repeat steps S424-S426 until the maximum number of iterations is reached to obtain the optimized model hyperparameters.

[0233] Step S5: Utilize the Transformer-BiGRU hybrid prediction model with ISOO hyperparameter optimization to perform predictions, obtain the final prediction results, and employ commonly used prediction evaluation metrics to measure model accuracy, including:

[0234] Step S51, input the multi-source feature matrix to the final built fusion ISOO hyperparameter optimized Transformer-BiGRU model, the multi-source feature matrix including One-Hot encoding, decomposed IMF component and meteorological factors.

[0235] Step S52, based on the structure of the Transformer-BiGRU hybrid prediction model and the characteristics of the photovoltaic output data, define the key hyperparameter search space;

[0236] Step S53, train the Transformer-BiGRU hybrid prediction model with the initial parameter set, train the model on the training set (80% data), and record the root mean square error (RMSE) of the validation set (20% data) as the fitness value; the calculation formula of RMSE is as follows:

[0237] ;

[0238] In the formula, y e is the actual value of the photovoltaic power generation power; is the predicted value of the photovoltaic power generation power, and n is the sample number.

[0239] Step S54, repeat steps S52-S53 until the specified number of iterations is reached or the performance target is reached. In each iteration, the ISOO algorithm will recommend a new parameter set for evaluation, gradually optimizing the hyperparameters of the model. Select the parameter configuration with the smallest error on the validation set as the setting of the Transformer-BiGRU hybrid prediction model.

[0240] Step S55, obtain the final photovoltaic output prediction result, and use common prediction evaluation indexes to measure the model precision, i.e. mean absolute error (MAE), root mean square error (RMSE), symmetric average absolute percentage error (SMAPE) and determination coefficient (R 2 ). The calculation formula of RMSE is shown in step S53; the calculation formulas of MAE, SMAPE and R 2 are as follows:

[0241] ;

[0242] ;

[0243] ;

[0244] In the formula, y e is the actual value of the photovoltaic power generation power; is the predicted value of the photovoltaic power generation power, and n is the sample number.

[0245] Example analysis:

[0246] This embodiment uses the photovoltaic output data of a 20MW photovoltaic power station in Gonghe County, Qinghai Province in 2016, and the data time resolution is 15 minutes. The abnormal values are detected and corrected, the missing values are filled by the interpolation method, and the meteorological features highly correlated with the photovoltaic power output are selected through correlation analysis. The Pearson correlation coefficients of the meteorological features and the actual output data are shown in Table 1. The temperature, scattered irradiance DHI, normal direct irradiance DNI, and total irradiance GHI are retained, and the maximum-minimum extreme value normalization is performed.

[0247] Table 1

[0248]

[0249] Input the 15-minute photovoltaic output data, temperature, scattered irradiance DHI, normal direct irradiance DNI, and total irradiance GHI data of the whole year of 2016; divide the photovoltaic output into sub-sequences by day, i.e. 365 photovoltaic output sub-sequences, 15 minutes per data point, and the sequence length of each day is 96. Calculate the DTW distance between the photovoltaic output sequences of the whole year, and calculate the 365x365 DTW distance matrix. At the same time, extract the frequency domain features of the daily sequence, calculate the photovoltaic output frequency band of each day, including the low frequency band, the medium frequency band and the high frequency band, i.e. the FSD matrix of each day is 1x3, and the FSD matrix of a year is 365x3.

[0250] Select 100 values of weight coefficients between (0, 1) through the NSGA-II algorithm, calculate the silhouette coefficient and intra-class distance under different weight coefficients, and obtain the best weight coefficient as 0.78. At this time, the maximum silhouette coefficient is 0.78, and the corresponding clustering number is 4; the minimum intra-class distance is 1.9. Obtain the fusion matrix, reduce it to 2 dimensions by t-SNE, and select the perplexity as 0 and the learning rate as 200 during dimension reduction, and iterate 1000 times. The dimension of the reduced matrix is 365x2. Perform spectral clustering on the reduced matrix, and select the Gaussian kernel width as =0.5 and K=4 (the maximum point of the silhouette coefficient) at this time. The final clustering result is 4 classes, including sunny days, rainy days, cloudy days, and special weather. The results after statistical analysis of the whole year data are shown in Table 2.

[0251] Table 2

[0252]

[0253] For different scenes, different decomposition methods are adopted. For sunny days, CEEMDAN decomposition is adopted and low-frequency components with SE < 0.2 are retained. For cloudy days, CEEMDAN decomposition is adopted and medium-frequency components with SE < 0.4 are retained. For rainy days, CEEMDAN decomposition is first adopted, and then VMD secondary decomposition is performed on the components with SE in the top 5, and finally the components with SE less than 0.5 are retained.

[0254] Taking the rainy day scene as an example, CEEMDAN decomposition is first performed to obtain IMF1-IMF11, and the sample entropy SE value is calculated, as shown in Table 3. VMD secondary decomposition is performed on the IMF (IMF1-IMF5) with SE in the top 5, and finally the components with SE < 0.5 are retained, obtaining 4 effective components, as shown in Table 4.

[0255] Table 3

[0256]

[0257] Table 4

[0258]

[0259] The hyperparameters of the improved ISOO algorithm optimized Transformer-BiGRU hybrid prediction model are defined as follows: the number of multi-head attention heads of the Transformer module is a discrete set of {2, 4, 6, 8}; the L2 regularization coefficient is a logarithmic space of [1e-4, 1e-1]; the Dropout is {0.1, 0.2, 0.3}; during model training, the training batch is set to 16 and the iteration is 100. The optimizer uses the Adam algorithm, and the initial learning rate is set to 0.1, then it is decayed by 0.1 times every 100 training periods until it reaches the minimum learning rate of 0.0001. The initial parameter set is used to train the Transformer-BiGRU hybrid prediction model, the model is trained on the training set (80% data), and the RMSE of the validation set (20% data) is recorded as the fitness value until the specified number of iterations is reached. In each iteration, the ISOO algorithm will recommend a new parameter set for evaluation, and gradually optimize the hyperparameters of the model. The parameter configuration with the smallest error on the validation set is selected as the setting of the Transformer-BiGRU hybrid prediction model, as shown in Table 5.

[0260] Table 5

[0261]

[0262] The four IMF components obtained by twice VMD decomposition of the rainy day scene and the corresponding meteorological factors are input into the trained mixed prediction model of Transformer-BiGRU to obtain the final prediction result under the rainy day scene. At the same time, the same operation is performed for sunny, cloudy and special meteorological scenes to obtain the prediction result under each weather type; at the same time, without weather typing and signal decomposition of the annual data, the original prediction model without ISOO algorithm optimization of super parameters is used for prediction. The prediction accuracy is measured by using RMSE, MAE, SMAPE and R 2 The prediction evaluation indexes are used to measure the prediction accuracy, and the results are shown in Table 6.

[0263] Table 6

[0264]

[0265] Example Two

[0266] The embodiment also discloses a computer device, including a memory, a processor and a computer program stored in the memory, the processor executes the computer program to realize the steps of the method in the embodiment one.

[0267] Example Three

[0268] The embodiment also discloses a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to realize the steps of the method in the embodiment one.

[0269] Example Four

[0270] The embodiment also discloses a computer program product, which includes a computer program, and the computer program is executed by a processor to realize the steps of the method in the embodiment one.

[0271] The above is only the preferred specific embodiment of the present application, but the protection scope of the present application is not limited to this, any person skilled in the art can easily think of changes or replacements within the technical range disclosed in the present application, which should be covered in the protection scope of the present application. Therefore, the protection scope of the present application should be limited by the protection scope of the claims.

Claims

1. A method for ultra-short-term photovoltaic power output prediction based on bimodal weather classification, characterized in that, The method comprises the following steps: acquiring meteorological data and photovoltaic output data of a research area and preprocessing the data; constructing a bimodal feature matrix according to the preprocessed data and combining dynamic time warping distance and frequency energy feature matrix; adopting a bimodal joint distance function to adaptively assign time-frequency domain feature weight coefficients, and processing the bimodal feature matrix through dimension reduction according to the assigned weight coefficients; realizing the classification of four types of weather, i.e., sunny, cloudy, rainy and special weather, according to the dimension-reduced bimodal feature matrix and using spectral clustering; adopting a corresponding signal decomposition method for different weather scenes according to the classification result; constructing a hybrid prediction model according to the signal decomposition result, and optimizing model hyperparameters in parallel through an improved star oscillation optimization algorithm; inputting the optimized hybrid prediction model into a fusion model for prediction to obtain a final prediction result; the process of optimizing model hyperparameters in parallel through the improved star oscillation optimization algorithm comprises: generating an initial oscillator population according to initial parameters, and perturbing the position of each oscillator through Logistic mapping; updating each generated oscillator individual through star earthquake oscillation according to the perturbed oscillator position, and performing mutation through Lévy flight perturbation; adding oscillatory motion influenced by sine and cosine functions to the oscillator after Lévy flight perturbation to generate a new position; evaluating the fitness according to the new position, performing mirror reverse mapping on the oscillator with poor fitness to generate a mirror solution and keeping the better position, keeping the particle position with the best fitness, and repeating the above steps until a maximum iteration number or a convergence threshold is reached to obtain the optimized model hyperparameters.

2. The method of claim 1, wherein the process of acquiring meteorological data and photovoltaic output data of a research area and preprocessing the data comprises: performing outlier cleaning, interpolating missing data, data standardization and key meteorological factor screening on the meteorological data and photovoltaic output data of the research area.

3. The method of claim 1, wherein the process of realizing the classification of four types of weather, i.e., sunny, cloudy, rainy and special weather, comprises: extracting complementary features from time domain and frequency domain according to time window slicing of input data; defining a bimodal joint distance function, performing dimension reduction on the bimodal feature matrix when solving the weight coefficient as a preset weight coefficient, performing spectral clustering on the dimension-reduced bimodal feature matrix to obtain a temporary clustering result, according to the time domain and frequency domain features; defining an optimization objective function according to the temporary clustering result, obtaining the optimal weight coefficient of the bimodal feature matrix through a multi-objective optimization algorithm, and repeatedly performing dimension reduction and spectral clustering using the optimal weight coefficient to obtain the spectral clustering result and dynamically mapping the weather type label.

4. The method of claim 1, wherein the process of adopting a corresponding signal decomposition method for different weather scenes according to the classification result comprises: if it is a sunny scene, performing CEEMDAN decomposition and calculating sample entropy, and only keeping the low-frequency component with sample entropy less than 0.2; if it is a cloudy scene, performing CEEMDAN decomposition and calculating sample entropy, and only keeping the medium-frequency component with sample entropy less than 0.

4. If it is a rainy day scene, CEEMDAN is used for decomposition first, and then the first five components of sample entropy are decomposed by VMD again, and finally the components with sample entropy less than 0.5 are reserved; If it is a special weather scene, decomposition is used and the components with sample entropy less than 0.6 are reserved.

5. The method of claim 1, wherein, The hybrid prediction model comprises a Transformer module and a bidirectional gated recurrent module.

6. The method of claim 1, wherein, The process of inputting the optimized hybrid prediction model into the fusion model for prediction to obtain the final prediction result comprises: According to the optimized hybrid prediction model, input the multi-source feature matrix including the encoded and decomposed components and the weather factors into the finally built fusion model; According to the characteristics of the model structure and the photovoltaic output data, define the hyperparameter search space; According to the defined hyperparameter search space, train the hybrid prediction model with the initial parameter set; According to the training result, repeat the above steps until a specified number of iterations is reached or a performance target is reached, and select the parameter configuration with the smallest error on the validation set as the setting of the hybrid prediction model; According to the final hybrid prediction model, the final photovoltaic output prediction result is obtained.

7. A computer apparatus comprising a memory, a processor, and a computer program stored on the memory, wherein the computer program, when executed by the processor, causes the processor to perform the method of any one of claims 1 to 6. The processor executes the computer program to implement the steps of the method of any one of claims 1-6.

8. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the steps of the method of any one of claims 1-6.

9. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the steps of the method of any one of claims 1-6.

Images