Natural gas load prediction method and device

Abstract

The invention provides a natural gas load prediction method and device. The natural gas load prediction method comprises the steps of obtaining a natural gas load sequence, a date sequence and exogenous variable data; decomposing a pre-obtained natural gas load sequence, and extracting a first-order period feature and a first-order trend feature of the natural gas load sequence; extracting second-order periodic features of the natural gas load sequence in a time period corresponding to the first-order periodic features; extracting second-order trend characteristics of the natural gas load sequence in a time period corresponding to the first-order trend characteristics; calculating the mutation probability of each time point in the natural gas load sequence to generate a mutation probability sequence; and aggregating and splicing results after the time domain mixed residual block processing, and inputting the results as input variables into an Informer generation model to carry out load prediction so as to generate a prediction result. According to the method, the processing capability of exogenous variables and the generalization capability of the model can be enhanced, and the prediction precision and robustness are improved.

Description

Technical Field

[0001] The present application relates to the technical field of natural gas load forecasting, in particular to the technical field related to machine learning in this technical field, and specifically to a natural gas load forecasting method and device. Background Art

[0002] In order for the natural gas industry to cope with seasonal peak regulation, emergency supply guarantee and supply shortages due to emergencies, it is necessary to accurately predict the natural gas load, aiming to accurately capture future fluctuations in natural gas load, thereby optimizing the supply and scheduling of natural gas. It also serves as an important boundary condition for the natural gas pipeline network simulation mechanism model, supporting the accuracy and real-time performance of the mechanism model solution results.

[0003] In recent years, with the rapid development of deep learning and time series forecasting technology, natural gas load forecasting based on neural networks has significantly improved the forecasting accuracy to a certain extent. However, despite the technological progress, the following key problems still need to be solved in the forecasting process:

[0004] (1) Insufficient consideration of exogenous variables: Natural gas load is affected by many exogenous factors, such as temperature and holidays. However, existing models often fail to fully consider the complexity and diversity of these exogenous variables, resulting in the prediction results of natural gas load may not reflect the actual demand fluctuations. How to effectively integrate and deal with these exogenous factors remains a challenging problem.

[0005] (2) Response to extreme weather and emergencies: Extreme weather conditions (such as cold waves, high temperatures, etc.) and emergencies (such as natural disasters, gas supply interruptions) have a huge impact on natural gas demand. Such abnormal events are often uncertain and random, and traditional forecasting models are difficult to capture or predict these extreme situations in advance, resulting in poor performance of forecasting models in actual emergencies.

[0006] (3) Model parameters are large and complex: Deep learning models, especially prediction models based on neural networks, usually contain a large number of parameters and are highly complex. This complexity not only increases the computational overhead of the model, but also easily leads to overfitting problems.

[0007] (4) Low prediction robustness: Natural gas load data often contains noise, missing values ​​or outliers, and load demand is affected by multiple highly correlated factors, which makes model predictions susceptible to fluctuations in data quality. Summary of the invention

[0008] In order to solve the above problems, the present invention provides a natural gas load forecasting method and device considering the fusion of mutation points and collaborative variables, so as to better adapt to the dynamic and mutation changes of natural gas load data containing more noise or abnormal values, and have better generalization.

[0009] In order to achieve the above objectives, in a first aspect, this party provides a natural gas load forecasting method, comprising:

[0010] Obtain natural gas load series, date series and exogenous variable data;

[0011] Decomposing the pre-acquired natural gas load sequence, and extracting the first-order periodic characteristics and the first-order trend characteristics of the natural gas load sequence;

[0012] Extracting the second-order periodic characteristics of the natural gas load sequence in the time period corresponding to the first-order periodic characteristics;

[0013] Extracting the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature;

[0014] Calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence;

[0015] splicing the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and inputting them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result;

[0016] Performing hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result;

[0017] Record the location information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate time series enhanced attention calculation results;

[0018] Inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the temporal pattern characteristics in the natural gas load data;

[0019] The results of the time domain mixed residual block processing are aggregated and spliced, and are input as input variables into the Informer generation model for load forecasting to generate a forecast result.

[0020] In some embodiments of the present invention, the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data are spliced ​​and input into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result, including:

[0021] Performing layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data;

[0022] Inputting the first normalized data into a first fully connected layer of an adaptive dynamic activation function, wherein the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence;

[0023] Perform Dropout processing on the output data of the first fully connected layer.

[0024] In some embodiments of the present invention, the step of inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the temporal pattern characteristics in the natural gas load data includes:

[0025] Performing layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data;

[0026] The second normalized data is processed by a second fully connected layer using a Tanh activation function and regularized by a Dropout layer.

[0027] In some embodiments of the present invention, the hybrid attention calculation includes:

[0028] Perform one-dimensional convolution operations on each incoming variable sequence to generate query matrix, value matrix, and keyword matrix, and input these matrices into the self-attention mechanism for calculation.

[0029] In a second aspect, the present application provides a natural gas load prediction device, the device comprising:

[0030] Data acquisition module, used to obtain natural gas load series, date series and exogenous variable data;

[0031] A first-order data extraction module, used to decompose the pre-acquired natural gas load sequence and extract the first-order periodic characteristics and first-order trend characteristics of the natural gas load sequence;

[0032] A second-order periodic feature extraction module, used to extract the second-order periodic feature of the natural gas load sequence in the time period corresponding to the first-order periodic feature;

[0033] A second-order trend feature extraction module, used to extract the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature;

[0034] A mutation probability sequence generation module is used to calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence;

[0035] A feature fusion result generation module is used to splice the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and input them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result;

[0036] A hybrid attention calculation result generation module, used for performing hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result;

[0037] A time series enhanced attention calculation result generation module is used to record the position information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate a time series enhanced attention calculation result;

[0038] A time pattern feature integration module, used for inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the time pattern features in the natural gas load data;

[0039] The prediction result generation module is used to aggregate and splice the results of the time domain mixed residual block processing, and input them as input variables into the Informer generation model for load prediction to generate prediction results.

[0040] In some embodiments of the present invention, the feature fusion result generating module includes:

[0041] A first normalized data generating unit, used for performing layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data;

[0042] A first normalized data input unit, used to input the first normalized data into a first fully connected layer of an adaptive dynamic activation function, where the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence;

[0043] An output data Dropout processing unit is used to perform Dropout processing on the output data of the first fully connected layer.

[0044] In some embodiments of the present invention, the temporal pattern feature integration module includes:

[0045] A second normalized data generating unit, used for performing layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data;

[0046] The second normalized data regularization unit is used to process the second normalized data using a second fully connected layer of a Tanh activation function and perform regularization through a Dropout layer.

[0047] In some embodiments of the present invention, the hybrid attention calculation result generation module includes:

[0048] The hybrid attention calculation result generation unit is used to perform a one-dimensional convolution operation on each incoming variable sequence to generate a query matrix, a value matrix, and a keyword matrix, and input these matrices into the self-attention mechanism for calculation.

[0049] In a third aspect, the present application provides a computer program product, including a computer program / instruction, which implements the steps of a natural gas load forecasting device when executed by a processor.

[0050] In a fourth aspect, the present application provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of a natural gas load forecasting method when executing the program.

[0051] In a fifth aspect, the present application provides a computer-readable storage medium having a computer program stored thereon, which implements the steps of a natural gas load forecasting method when executed by a processor.

[0052] Compared with the prior art, the present invention has the following beneficial effects:

[0053] The present invention can more fully consider the impact of exogenous factors (such as temperature, holidays, etc.) on natural gas load through the collaborative variable fusion module, thereby greatly enhancing the processing capability of exogenous variables. Through the synergistic effect of multiple factors, the scenario adaptability of the prediction model is improved, the problem of insufficient processing of exogenous variables in the existing model is solved, and the prediction results are more in line with actual needs.

[0054] The present invention adopts the Bayesian change point detection algorithm, which can effectively respond to emergencies and extreme weather. The algorithm can capture the mutation points in the natural gas load data in advance, thereby improving the rapid response capability of the prediction model under extreme weather or emergencies, and enhancing the robustness and emergency management capability of the natural gas supply system.

[0055] The present invention significantly improves the generalization ability of the prediction model by introducing the time domain mixed residual block and the feature domain mixed residual block. This technology can better handle the noise, missing values ​​and outliers in the natural gas load data, making the prediction model more robust and able to adapt to complex data environments, ensuring the reliability and stability of the prediction results.

[0056] The present invention effectively reduces the computational complexity of the model by combining feature fusion and adaptive dynamic activation function. The combination of GeLU function and the probability of load data mutation at a time point reduces the risk of overfitting of the model, while ensuring that the model can run efficiently and maintain high prediction accuracy under limited computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0058] Figure 1 A flow chart of a natural gas load forecasting method provided in an embodiment of the present application;

[0059] Figure 2 A flowchart of step 600 provided in an embodiment of the present application;

[0060] Figure 3 A flowchart of step 900 provided in an embodiment of the present application;

[0061] Figure 4 A schematic diagram of the principle of a natural gas load forecasting system provided in a specific application example of this application;

[0062] Figure 5 A schematic diagram of the inner loop of the STL data decomposition algorithm provided in the specific application example of this application;

[0063] Figure 6 This is a diagram of the overall prediction effect of the data set provided in the specific application example of this application;

[0064] Figure 7 A schematic diagram of the structure of a natural gas load prediction device provided in a specific application example of this application;

[0065] Figure 8 This is a structural diagram of the feature fusion result generation module 06 provided in the specific application example of this application;

[0066] Fig. 9 This is a schematic diagram of the structure of the time pattern feature integration module 09 provided in the specific application example of this application;

[0067] Fig.10 This is a schematic diagram of the structure of an electronic device provided in a specific application example of this application. DETAILED DESCRIPTION

[0068] In order to make the purpose, technical solution and advantages of the embodiments of the present application clearer, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0069] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application may adopt the form of a computer program product implemented in one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that include computer-usable program code.

[0070] It should be noted that the terms "including" and "having" and any variations thereof in the specification and claims of the present application and the above-mentioned drawings are intended to cover non-exclusive inclusions. For example, a process, method, system, product or apparatus comprising a series of steps or units is not necessarily limited to those steps or units clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or apparatuses.

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

[0072] The embodiment of the present application provides a specific implementation method of a natural gas load prediction method, see Figure 1 , the method specifically includes the following contents:

[0073] Step 100: Obtain natural gas load series, date series and exogenous variable data;

[0074] Step 200: Decomposing the pre-acquired natural gas load sequence, and extracting the first-order periodic characteristics and the first-order trend characteristics of the natural gas load sequence;

[0075] Specifically, the STL (Seasonal Trend Decomposition) algorithm can be used to decompose the natural gas load sequence to extract the first-order periodic characteristics and the first-order trend characteristics.

[0076] Step 300: extracting the second-order periodic characteristics of the natural gas load sequence in the time period corresponding to the first-order periodic characteristics;

[0077] Specifically, the first-order periodic features can be extracted using Fourier transform and converted into second-order periodic features.

[0078] Step 400: extracting the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature;

[0079] Specifically, the Markov chain can be used to extract the first-order trend features and convert them into second-order trend features.

[0080] Step 500: Calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence;

[0081] Based on the Bayesian change point detection algorithm, the mutation probability of each time point in the natural gas load sequence is calculated to generate a mutation probability sequence;

[0082] Step 600: splicing the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and inputting them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result;

[0083] Step 700: performing hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result;

[0084] It can be understood that step 700 can effectively improve the correlation between complex data.

[0085] Step 800: Record the location information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate a time series enhanced attention calculation result;

[0086] Preferably, the position information of each data point in the natural gas load sequence is recorded by absolute position coding.

[0087] Step 900: inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the temporal pattern characteristics in the natural gas load data;

[0088] Step 1000: Aggregate and splice the results of the time domain mixed residual block processing, and input them as input variables into the Informer generation model for load forecasting to generate a forecast result.

[0089] From the above description, it can be seen that the embodiment of the present application provides a natural gas load forecasting method, which starts from the two perspectives of mutation point and synergistic variable fusion, and through synergistic variable fusion, can more fully consider the impact of exogenous factors (such as temperature, holidays, etc.) on natural gas load, thereby greatly enhancing the processing ability of exogenous variables. Through the synergistic effect of multiple factors, the scenario adaptability of the prediction model is improved, and the problem of insufficient processing of exogenous variables by the existing model is solved, so that the prediction results are more in line with actual needs.

[0090] In some embodiments of the present invention, see Figure 2 , step 600 includes:

[0091] Step 601: performing layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data;

[0092] Specifically, the input features are normalized layer by layer. The features involved include second-order periodic features, second-order trend features, location information of sequence data points, and exogenous variables. This step aims to standardize different features to ensure the stability and effectiveness of subsequent calculations.

[0093] Step 602: input the first normalized data into a first fully connected layer of an adaptive dynamic activation function, where the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence;

[0094] The normalized data is input into the fully connected layer with the adaptive dynamic activation function. In this layer, the adaptive dynamic activation function uses the product of the GeLU function and the average value of the mutation probability sequence, which can dynamically scale the output value of the activation function, thereby improving the model's ability to handle complex data features.

[0095] Step 603: Perform Dropout processing on the output data of the first fully connected layer.

[0096] In order to prevent the model from overfitting, the processed data passes through the Dropout layer, which is regularized by randomly blocking some neurons to enhance the generalization ability of the model. Subsequently, the data passes through the fully connected layer without activation function, which does not contain an activation function and is only used for linear transformation. Finally, the data passes through the Dropout layer again for further regularization, and then the final result is output.

[0097] In some embodiments of the present invention, see Figure 3 , step 900 includes:

[0098] Step 901: performing layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data;

[0099] The input hybrid attention calculation results and temporal enhanced attention calculation results are layer normalized. The purpose of this step is to standardize the input different attention calculation results to ensure the stability and consistency of subsequent calculations.

[0100] Step 902: Process the second normalized data using a second fully connected layer with a Tanh activation function, and regularize the data using a Dropout layer.

[0101] Specifically, the data after the second normalization is input to the fully connected layer containing the activation function. In this layer, the Tanh activation function is used. Unlike the activation function in the feature domain mixed residual block, the Tanh activation function can handle nonlinear features and compress the output value to the range of -1 to 1.

[0102] In order to prevent the model from overfitting, the data is regularized by the Dropout layer, which randomly blocks some neurons and enhances the generalization ability of the model. The data then passes through a fully connected layer without an activation function, which only performs linear transformations and does not contain any nonlinear activations. Finally, the data is further regularized by the Dropout layer and the processing results are output.

[0103] In some embodiments of the present invention, the hybrid attention calculation includes:

[0104] Perform one-dimensional convolution operations on each incoming variable sequence to generate query matrix, value matrix, and keyword matrix, and input these matrices into the self-attention mechanism for calculation.

[0105] Specifically, a one-dimensional convolution operation is performed on each variable sequence passed in. The purpose of this step is to extract local features in the sequence through the convolution operation and generate the first query matrix and the first value matrix. These matrices are used for subsequent attention calculations. The last data point is extracted from the sequence as the first keyword matrix. The keyword matrix represents the important information in the sequence and is used to calculate the correlation between sequences in the attention mechanism. The generated first query matrix, first value matrix, and first keyword matrix are input into the self-attention mechanism together. In the self-attention mechanism, the system calculates the relationship and weight between each sequence, and outputs the result after attention calculation. This process can capture the global dependency between sequences, thereby improving the prediction performance of the model.

[0106] From the above description, it can be seen that the embodiment of the present application provides a natural gas load forecasting method, which starts from the two perspectives of mutation point and synergistic variable fusion, and through synergistic variable fusion, can more fully consider the impact of exogenous factors (such as temperature, holidays, etc.) on natural gas load, thereby greatly enhancing the processing ability of exogenous variables. Through the synergistic effect of multiple factors, the scenario adaptability of the prediction model is improved, and the problem of insufficient processing of exogenous variables by the existing model is solved, so that the prediction results are more in line with actual needs.

[0107] To further illustrate the present solution, the present application also provides a specific application example of a natural gas load forecasting system, which specifically includes the following contents.

[0108] A natural gas load forecasting system provided in the present application includes: a data acquisition module, a collaborative variable fusion module and a prediction model prediction module, wherein the collaborative variable fusion network is used for first-order period feature decomposition, first-order trend feature decomposition, time series change point probability detection, one-dimensional space domain mixing, second-order period feature extraction, second-order trend feature extraction, feature domain mixing, mixed attention calculation, time series enhanced attention calculation, time domain mixing and aggregation splicing; the prediction model is used to predict the future load changes as an Informer generative model.

[0109] See also Figure 4 The working principle of the above natural gas load forecasting system is as follows: First, the natural gas load sequence, date sequence and exogenous variables are obtained. The natural gas load sequence is decomposed into two sequences, the first-order periodic feature and the first-order trend feature, by the time series STL data decomposition algorithm. Then, the first-order periodic feature is extracted by Fourier transform and transformed into the second-order periodic feature, and the first-order trend feature is extracted by Markov chain and transformed into the second-order trend feature; the probability of mutation at any time point in the natural gas load sequence is calculated by the Bayesian change point detection algorithm, and the mutation probability sequence of any time point in the sequence is obtained, which is used as the prior information for load forecasting; the second-order periodic feature and the second-order trend feature after the above calculation are used as the prior information for load forecasting. The features, mutation probability sequences, position information of sequence data points and exogenous variables are spliced ​​and input into the feature domain mixing layer together to obtain the result after feature fusion, and then the mixed attention calculation is performed; the position information of each data point in the natural gas load sequence is recorded by absolute position encoding to form a one-dimensional spatial domain, which is used as input together with the date sequence corresponding to the natural gas load data to calculate the time series enhanced attention, and the calculated results and the mixed attention calculation results are input into the time domain mixing layer together to fuse the previous features with the time pattern of the complex natural gas load data; the results of the time domain mixing layer are aggregated and spliced ​​and used as input variables to enter the Informer generative model for prediction, thereby obtaining the prediction results.

[0110] The feature domain hybrid residual block mentioned above includes: layer normalization, fully connected layer, adaptive dynamic activation function, Dropout, etc. First, the feature domain hybrid residual block performs layer normalization on the second-order periodic features, second-order trend features, location information of sequence data points, and exogenous variables, and then passes through a fully connected layer containing an adaptive dynamic activation function, where the adaptive dynamic activation function is a GeLU function multiplied by the average value of the mutation probability sequence, which can scale the activation function value, and then passes through Dropout, a fully connected layer without an activation function, and Dropout calculation before output.

[0111] The time-domain mixed residual block mentioned above includes: layer normalization, fully connected layer, activation function, Dropout, etc. First, the time-domain mixed residual block performs layer normalization on the mixed attention calculation results of the input variables and the temporal enhanced attention calculation results, and then passes through the fully connected layer with the activation function. The activation function here is different from the feature domain mixed residual block, which is the Tanh activation function. Then it passes through Dropout, the fully connected layer without activation function, and the Dropout calculation before output.

[0112] The hybrid attention calculation mentioned above includes: one-dimensional convolution and self-attention calculation. Perform one-dimensional convolution operation on each variable sequence passed in, generate the first query matrix and the first value matrix, extract the last data point in the sequence as the first keyword matrix, and input the three matrices into the self-attention mechanism for calculation and output.

[0113] Furthermore, each module in a natural gas load forecasting system has the following definitions:

[0114] Data acquisition module: used to obtain natural gas load series, date series and exogenous variable data, and pre-process the acquired data. The natural gas load series is defined as X, X = [x1, x2, L, x T ]∈R T represents the natural gas load data obtained by a gas company after sampling over a period of time T, X i represents the natural gas load data sampled by a gas company at the i-th time point, and T represents the number of sampling times, that is, the total number of sampling time points. The date series is concatenated with other exogenous variables. The exogenous variable is defined as Z, It means there are N exogenous variables, and for any exogenous variable Its data length is also T, corresponding to the information of natural gas load data at T sampling times. The exogenous variables include but are not limited to: date sequence, temperature (maximum temperature, minimum temperature).

[0115] After obtaining the data, data preprocessing is required, which includes normalization, data set division, and sliding window division.

[0116] The natural gas load data and exogenous variables are normalized, and the calculation method is shown in formula (1) and formula (2).

[0117]

[0118] Formula (1) represents the normalization of natural gas load data, where μ is the mean of the input natural gas load data and is the variance of the input natural gas load data. Formula (2) represents the normalization of exogenous variables. Since there are N exogenous variables, the normalization operation needs to be performed on each exogenous variable separately. The normalization operation can eliminate the dimension of the data and avoid the impact of data size differences.

[0119] To ensure that the model does not leak data, the preprocessed data set is first divided into training set, validation set and test set in a ratio of 8:1:1. The training set is used to train the network parameters of the collaborative variable fusion network and the Informer prediction model, while the validation set is used to evaluate and tune the trained network. The final test set is to verify the generalization ability of the prediction method for unknown data sets.

[0120] In order to improve the training efficiency, the results after aggregation and splicing are divided into sliding windows and entered into the model for training in batches. The length of the window is recorded as w, and the number of steps the window moves forward each time is recorded as s. i is the starting position of each sliding window, that is, the i-th time point. Taking the natural gas load sequence as an example, the sliding window processing can be expressed as, X i =[x i ,x i+1 ,...,x i+w-1 ], where window i contains w data points starting from the i-th time point. The divided data are input into the collaborative variable fusion module in the order of the length of the sliding window and the time step. For the convenience of explanation, the entry of each sliding window is represented in the following, and the natural load shedding data is recorded as OT and the exogenous variable is recorded as M.

[0121] The collaborative variable fusion module is used to achieve the following functions: First-order period feature decomposition: used to extract the most basic periodic patterns in the data, helping to capture short-term periodic changes in time series such as natural gas load. First-order trend feature decomposition: used to capture long-term trends in the data, and by decomposing the trend part of the data, it helps to identify the macro trend of load changes. Time series change point probability detection: detect possible change points in the time series (such as abnormal fluctuations or trend changes) to identify sudden changes or turning points in the load data as prior information for load data prediction. One-dimensional spatial domain mixing: identify and encode the front and back positions of different data points in the variable. Second-order period feature extraction: based on the extraction of first-order period features, further capture more complex periodic patterns and identify long-term or multi-level periodic changes. Second-order trend feature extraction: based on the first-order trend, further extract deeper trend information to help accurately describe long-term change trends. Feature domain mixing: fuse the features extracted from different decompositions to improve the comprehensive understanding and processing capabilities of the prediction model for various features. Hybrid attention calculation: weight different features through the attention mechanism to highlight important features, help the model better focus on key information, and reduce the interference of irrelevant information. Time series enhanced attention calculation: weight the sequence information in the time dimension to capture the importance of specific time periods in the time series to more accurately identify the changes in patterns over time. Time domain mixing: fuse feature information at different time scales to help better capture the interaction and impact of long-term and short-term changes. Aggregation splicing: splice the processed multiple features together as the final input feature to ensure that all relevant feature information is included in the model prediction and enhance the accuracy of the prediction.

[0122] First-order period feature decomposition and first-order trend feature decomposition: Perform time series STL data decomposition on the input natural gas load data to extract first-order period and trend features. The original natural gas load data sequence can be considered as the accumulation of trend terms, seasonal terms and residual terms, as shown in the following formula:

[0123] OT t =T t +S t +R t (3)

[0124] In the formula, OT t is the original sequence; T t is the trend item; S t is the seasonal term; R t is the residual term.

[0125] The STL decomposition process consists of two parts: an inner loop and an outer loop. The inner loop uses Loess to smooth the trend term and seasonal term, and the outer loop calculates the robust weight based on the results of the inner loop to reduce the impact of noise on the next inner loop.

[0126] Assume T t (k) is the trend item sequence after the kth inner cycle, S t (k) is the seasonal term sequence after the kth inner loop, with the initial value T t (0) is 0, and the inner loop steps are as follows Figure 5 As shown, there are six process steps, namely, detrending, seasonal subsequence smoothing, seasonal subsequence low-pass filtering, desmoothing seasonal subsequence trend, deseasonalization and trend smoothing, and residual component calculation, which are continuously iterated until the residual component converges.

[0127] At the end of the inner loop, there may be outliers in the residual term, which will affect the decomposition result. Therefore, it is necessary to use the outer loop to assign weights for noise reduction. The calculation method of robust weights is as follows:

[0128] h=6median(R t |) (4)

[0129]

[0130] ρ t =B(R t | / h) (6)

[0131] In the next cycle, the robust weights are used to adjust the weights in Loess in the following way:

[0132] v i (x) new =ρ t v i (x) old (7)

[0133] In the inner and outer loops, Loess is used to smooth the data and assign weights. The weight calculation principle is as follows:

[0134] First, select the q points closest to the fitting point x as the range of the field, and then assign weights to each point according to the distance between these q points and the fitting point x. Define the function W:

[0135]

[0136] Then calculate the distances between q points and the fitting point x respectively, and record the maximum distance as λ q (x), then for a point x near the fitting point x i , whose weight v i The calculation method of (x) is as follows:

[0137]

[0138] Therefore, points closer to the fit have larger weights, and the weights decrease as they move away from the fit, until the farthest point, the qth point, has a weight decaying to 0.

[0139] Time series change point probability detection: Perform Bayesian change point detection on the input natural gas load data to obtain the change point probability sequence in the sequence. Define the change point position as t, and the change point probability sequence P can be obtained by formula (10).

[0140]

[0141] Among them, P(OT|t) can be expressed as: μ1, and μ2, are the mean and variance of the data before and after the change point, respectively, and P(t) means that the natural gas load data follows a normal distribution.

[0142] One-dimensional spatial domain mixing: Perform corresponding encoding combinations on each group of data point positions, where the position embedding formulas are shown in formulas (11) and (12). Sum the encoded sequences to obtain the result L after one-dimensional spatial domain mixing.

[0143]

[0144] Among them, pos is the absolute position index of the feature sequence, d model is the dimension of the input feature vector, j is the position index of the feature vector, 2j represents an even position, and 2j+1 represents an odd position.

[0145] Second-order periodic feature extraction: The first-order periodic feature T after decomposition of STL data t Perform Fourier transform to extract the second-order periodic features. t Apply the fast Fourier transform, as shown in formula (13), where N is S t The data length is k, and k is the frequency index. Set the threshold τ and identify the main frequency components in the Fourier transform result through formula (14). Use the inverse Fourier transform to convert the frequency domain information back to the time domain as shown in formula (15), denoted as T t (2) .

[0146]

[0147] Second-order trend feature extraction: The first-order trend feature S after decomposition of STL data t Perform Markov extraction of second-order trend features. By formula (16), the first-order trend feature S t Discretization forms the state space S, S = {s1, s2, ..., s N}, where t represents each moment, and each s i Represents a state of trend characteristics, and constructs the state transfer matrix as shown in formula (17), where the matrix element p ij Indicates that from state s i At the next moment, it transfers to state s j The probability of the state transition is 1, and the sum of the state probabilities of all possible transitions is 1. The steady-state distribution π of the system is calculated by formula (18), where π=(π1,π2,...,π N ) is a steady-state probability vector and satisfies the normalization condition The new steady-state distribution π is the global feature of the first-order trend feature, that is, the second-order trend feature, denoted by S t (2) .

[0148]

[0149] πP=π (18)

[0150] Feature domain mixed residual block: The feature domain mixed residual block combines multiple levels of feature processing to improve the model's modeling ability for different features. t (2) , second-order trend characteristics S t (2) and the exogenous variable M, and then input them into the feature domain to form F = concat(T t (2) ,S t (2) ,M), first, a similar form of normalization is performed on each row of the input multidimensional array through formula (1), and then it enters the fully connected layer with an adaptive dynamic activation function. The adaptive dynamic activation function is: in, The average value of the mutation probability sequence P(t|y) introduces the mutation probability information within the sliding window length, so that the output of the activation function can dynamically respond to the mutation of the sequence. It helps to enhance the model's sensitivity to mutation points in the time series, and is especially suitable for capturing non-stationary changes in time series. The processed results are then passed through Dropout, a fully connected layer without an activation function, and Dropout to output the results, which are recorded as a multidimensional array C.

[0151] Hybrid attention calculation: Perform a one-dimensional convolution attention operation on each variable sequence passed in. The formula is (19), where O Conis the output matrix after convolution. W is the convolution kernel weight matrix, with dimension (k, d), k is the size of the convolution kernel, d is the depth of the convolution kernel, and b is the bias vector. Generate the first query matrix Q and the first value matrix V, extract the last data point in the sequence as the first keyword matrix K, and input the three matrices into the self-attention mechanism for calculation and output. The self-attention calculation formula is formula (20), where d k is the dimension of the keyword matrix, used for normalization to stabilize the gradient.

[0152] O Con =Conv1D(C,W,b) (19)

[0153]

[0154] Temporal enhanced attention calculation: The result L after mixing in the one-dimensional spatial domain is used to calculate the query and key vectors with temporal order information through formula (), Q′=(X+P)W Q ,K′=(X+P)W K , and then the time-series enhanced features are calculated using formula (21), where B(t) is the time dependency term used to capture the dependency between time steps.

[0155]

[0156] Temporal hybrid residual block: The hybrid attention calculation result Attention(Q, K, V) and the temporal enhanced attention calculation result TemporalAttention(Q′, K′, V) are concatenated and input into the temporal domain for fusion. First, layer normalization is performed on it using formula (22), where LayerNorm represents the layer normalization operation. The normalized result passes through the fully connected layer and then is processed using formula (23), where W1 and W2 are the weight matrices of the fully connected layer, b1 and b2 are biases, and Tanh is the activation function. After the activation function, the Dropout operation is applied to reduce the risk of overfitting, and then the output is passed through the fully connected layer and the Dropout operation.

[0157]

[0158]

[0159] Aggregation splicing: mainly multiple variable sequences in time domain mixing results Perform full-connected linear layer calculations to obtain aggregated results Where W o ∈R 1×K , is the weight of the Kth variable sequence.

[0160] Prediction model prediction module: used to predict the future changes in natural gas load based on the output data of the collaborative variable fusion module. The aggregated and spliced ​​data are input into the Informer model in the order of the length of the sliding window and the time step. After multiple rounds of network parameter training, the optimal network parameters are obtained to predict the results. i The Informer model is fed in sequence for batch training. Each subsequence is trained multiple times, and the model updates its parameters until it converges to the optimal parameters. The objective loss function is: The network parameters of the model are continuously optimized through back propagation and gradient descent. After multiple rounds of network parameter training, the optimal parameter set is obtained. The Informer model can then accurately predict subsequent input data and output prediction results.

[0161] In addition, in order to verify the prediction effect of the prediction method of the present invention compared with other existing prediction methods, this embodiment uses real natural gas load data to predict the above method and other prediction methods.

[0162] The historical data of natural gas load of a city's gas company from May 16, 2015 to August 1, 2023 were collected. After investigation, it was found that weather and other related factors have the greatest impact and the highest correlation in the short-term natural gas demand forecast. Therefore, the historical data of local weather and other related influencing factors were crawled from the weather network, totaling 3,000 sets of data.

[0163] In order to evaluate the prediction effect of the model, the following four evaluation indicators are introduced:

[0164] (1) Mean absolute percentage error (MAPE)

[0165]

[0166] (2) Root mean square error (RMSE)

[0167]

[0168] (3) Mean absolute error (MAE)

[0169]

[0170] (4) Mean Square Error (MSE)

[0171]

[0172] (5) Coefficient of determination (R2)

[0173]

[0174] Where n is the number of samples; y i is the true value of the sample; Y i is the predicted value of the sample.

[0175] The smaller the values ​​of MAPE, RMSE, MAE, and MSE, the better the prediction effect of the model and the more accurate the prediction. The larger the value of R2, the better the data fits the model.

[0176] The acquired data is trained, tested, and predicted to obtain the prediction results as shown in Table 1 and Figure 6 shown.

[0177] Table 1 Comparison of evaluation indicators

[0178]

[0179] Based on the same inventive concept, the embodiments of the present application also provide a natural gas load forecasting device, which can be used to implement the method described in the above embodiments, such as the following embodiments. It should be noted here that the various modules in this embodiment correspond one-to-one to the various steps in Example 1, and the specific implementation process is the same, which will not be repeated here. As used below, the term "module" can implement a combination of software and / or hardware for a predetermined function. Although the system described in the following embodiments is preferably implemented in software, the implementation of hardware, or a combination of software and hardware, is also possible and conceivable.

[0180] The embodiment of the present application provides a specific implementation of a natural gas load prediction device capable of implementing a natural gas load prediction method, see Figure 7 , a natural gas load forecast specifically includes the following contents:

[0181] Data acquisition module 01, used to obtain natural gas load series, date series and exogenous variable data;

[0182] The first-order data extraction module 02 is used to decompose the pre-acquired natural gas load sequence and extract the first-order periodic characteristics and the first-order trend characteristics of the natural gas load sequence;

[0183] A second-order periodic feature extraction module 03, used to extract the second-order periodic feature of the natural gas load sequence in the time period corresponding to the first-order periodic feature;

[0184] A second-order trend feature extraction module 04, used to extract the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature;

[0185] The mutation probability sequence generation module 05 is used to calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence;

[0186] The feature fusion result generating module 06 is used to splice the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and input them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result;

[0187] A hybrid attention calculation result generating module 07 is used to perform hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result;

[0188] The time series enhanced attention calculation result generation module 08 is used to record the position information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate a time series enhanced attention calculation result;

[0189] A time pattern feature integration module 09, used for inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the time pattern features in the natural gas load data;

[0190] The prediction result generation module 10 is used to aggregate and splice the results of the time domain mixed residual block processing, and input them as input variables into the Informer generation model for load prediction to generate prediction results.

[0191] In some embodiments of the present invention, see Figure 8 , the feature fusion result generating module 06 comprises:

[0192] A first normalized data generating unit 61 is used to perform layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data;

[0193] A first normalized data input unit 62, used to input the first normalized data into a first fully connected layer of an adaptive dynamic activation function, wherein the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence;

[0194] The output data Dropout processing unit 63 is used to perform Dropout processing on the output data of the first fully connected layer.

[0195] In some embodiments of the present invention, see Fig. 9 , the time pattern feature integration module 09 includes:

[0196] A second normalized data generating unit 91 is used to perform layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data;

[0197] The second normalized data regularization unit 92 is used to process the second normalized data using a second fully connected layer of a Tanh activation function and perform regularization through a Dropout layer.

[0198] In some embodiments of the present invention, the mixed attention calculation result generation module 07 includes:

[0199] The hybrid attention calculation result generation unit is used to perform a one-dimensional convolution operation on each incoming variable sequence to generate a query matrix, a value matrix, and a keyword matrix, and input these matrices into the self-attention mechanism for calculation.

[0200] From the above description, it can be seen that the embodiment of the present application provides a natural gas load forecasting device, which starts from the two perspectives of mutation point and synergistic variable fusion, and through synergistic variable fusion, can more fully consider the impact of exogenous factors (such as temperature, holidays, etc.) on natural gas load, thereby greatly enhancing the processing ability of exogenous variables. Through the synergistic effect of multiple factors, the scenario adaptability of the prediction model is improved, and the problem of insufficient processing of exogenous variables by the existing model is solved, so that the prediction results are more in line with actual needs.

[0201] The embodiment of the present application also provides a specific implementation of an electronic device capable of implementing all steps in the natural gas load forecasting method considering the fusion of mutation points and cooperative variables in the above embodiment, see Fig.10 , electronic equipment specifically includes the following:

[0202] Processor (processor) 1201, memory (memory) 1202, communication interface (CommunicationsInterface) 1203 and bus 1204;

[0203] The processor 1201 , the memory 1202 , and the communication interface 1203 communicate with each other via the bus 1204 ; the communication interface 1203 is used to realize information transmission between relevant devices such as the server-side device, the power measurement device, and the user-side device.

[0204] The processor 1201 is used to call the computer program in the memory 1202, and when the processor executes the computer program, all the steps in the natural gas load forecasting method in the above embodiment are implemented.

[0205] The embodiments of the present application also provide a computer-readable storage medium capable of implementing all the steps of the natural gas load forecasting method in the above embodiments. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements all the steps of the natural gas load forecasting method in the above embodiments.

[0206] Each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. In particular, for the hardware + program embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiment.

[0207] The above is a description of a specific embodiment of the specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in an order different from that in the embodiments and still achieve the desired results. In addition, the processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0208] Although the present application provides method operation steps such as embodiments or flow charts, more or fewer operation steps may be included based on conventional or non-creative labor. The order of steps listed in the embodiments is only one way of executing the order of many steps and does not represent the only execution order. When the actual device or client product is executed, it can be executed in the order of the method shown in the embodiments or the drawings or in parallel (for example, in a parallel processor or multi-threaded processing environment).

[0209] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0210] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0211] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0212] Specific embodiments are used in this application to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core idea of ​​this application. At the same time, for those skilled in the art, according to the idea of ​​this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as a limitation on this application.

Claims

1. A natural gas load forecasting method, characterized in that: Includes: Obtain natural gas load series, date series and exogenous variable data; Decomposing the pre-acquired natural gas load sequence, and extracting the first-order periodic characteristics and the first-order trend characteristics of the natural gas load sequence; Extracting the second-order periodic characteristics of the natural gas load sequence in the time period corresponding to the first-order periodic characteristics; Extracting the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature; Calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence; splicing the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and inputting them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result; Performing hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result; Record the location information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate time series enhanced attention calculation results; Inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the temporal pattern characteristics in the natural gas load data; The results of the time domain mixed residual block processing are aggregated and spliced, and are input as input variables into the Informer generation model for load forecasting to generate a forecast result.

2. The natural gas load forecasting method according to claim 1, characterized in that: The second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data are spliced ​​and input into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result, including: Performing layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data; Inputting the first normalized data into a first fully connected layer of an adaptive dynamic activation function, wherein the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence; Perform Dropout processing on the output data of the first fully connected layer.

3. The natural gas load forecasting method according to claim 1 or 2, characterized in that: The step of inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the temporal pattern characteristics in the natural gas load data includes: Performing layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data; The second normalized data is processed by a second fully connected layer using a Tanh activation function and regularized by a Dropout layer.

4. The natural gas load forecasting method according to claim 1, characterized in that: The hybrid attention calculation includes: Perform one-dimensional convolution operations on each incoming variable sequence to generate query matrix, value matrix, and keyword matrix, and input these matrices into the self-attention mechanism for calculation.

5. A natural gas load forecasting device, characterized in that: Includes: Data acquisition module, used to obtain natural gas load series, date series and exogenous variable data; A first-order data extraction module, used to decompose the pre-acquired natural gas load sequence and extract the first-order periodic characteristics and first-order trend characteristics of the natural gas load sequence; A second-order periodic feature extraction module, used to extract the second-order periodic feature of the natural gas load sequence in the time period corresponding to the first-order periodic feature; A second-order trend feature extraction module, used to extract the second-order trend feature of the natural gas load sequence in the time period corresponding to the first-order trend feature; A mutation probability sequence generation module is used to calculate the mutation probability at each time point in the natural gas load sequence to generate a mutation probability sequence; A feature fusion result generation module is used to splice the second-order periodic feature, the second-order trend feature, the mutation probability sequence, the position information of the sequence data points and the exogenous variable data and input them into the feature domain mixed residual block to perform feature fusion to generate a feature fusion result; A hybrid attention calculation result generation module, used for performing hybrid attention mechanism calculation on the feature fusion result to generate a hybrid attention calculation result; A time series enhanced attention calculation result generation module is used to record the position information of each data point in the natural gas load sequence, and perform time series enhanced attention calculation in combination with the date sequence to generate a time series enhanced attention calculation result; A time pattern feature integration module, used for inputting the temporal enhanced attention calculation result and the hybrid attention calculation result into a time domain hybrid residual block to integrate the time pattern features in the natural gas load data; The prediction result generation module is used to aggregate and splice the results of the time domain mixed residual block processing, and input them as input variables into the Informer generation model for load prediction to generate prediction results.

6. The natural gas load prediction device according to claim 5, characterized in that: The feature fusion result generation module includes: A first normalized data generating unit, used for performing layer normalization processing on the input second-order periodic features, the second-order trend features, the position information of the sequence data points and the exogenous variables to generate first normalized data; A first normalized data input unit, used to input the first normalized data into a first fully connected layer of an adaptive dynamic activation function, where the adaptive dynamic activation function is a GeLU function multiplied by an average value of a mutation probability sequence; An output data Dropout processing unit is used to perform Dropout processing on the output data of the first fully connected layer.

7. The natural gas load prediction device according to claim 5 or 6, characterized in that: The time pattern feature integration module includes: A second normalized data generating unit, used for performing layer normalization processing on the input hybrid attention calculation result and the temporal enhanced attention calculation result to generate second normalized data; The second normalized data regularization unit is used to process the second normalized data using a second fully connected layer of a Tanh activation function and perform regularization through a Dropout layer.

8. The natural gas load prediction device according to claim 5, characterized in that: The hybrid attention calculation result generation module includes: The hybrid attention calculation result generation unit is used to perform a one-dimensional convolution operation on each incoming variable sequence to generate a query matrix, a value matrix, and a keyword matrix, and input these matrices into the self-attention mechanism for calculation.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the natural gas load forecasting method according to any one of claims 1 to 4 are implemented.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the program implements the steps of the natural gas load forecasting method according to any one of claims 1 to 4.

Images