Distributed industrial energy operation optimization platform automatically constructing intelligent models and algorithms

Abstract

A distributed industrial energy operation optimization platform which is capable of automatically constructing intelligent models and algorithms, is divided into three parts: a modeling terminal, a background service and a human-computer interface. The models like data pre-processing, energy generation-consumption-storage trend forecasting and optimal scheduling decision models are encapsulated in the modeling terminal as different visualization modules facing with multiple categories production scenarios, by dragging which the complex functional models can be realized conveniently. The background service is capable of automatically constructing the training samples and the production plans/manufacturing signals series according to the device model requirements of each edge side, interacts with the trained intelligent models through corresponding interfaces, and the computing results are saved in the specified relational database. The computing results are displayed through a friendly customer human-computer interface, and the real-time state of current working condition can also be adjusted.

Claims

1. A distributed industrial energy operation optimization platform automatically constructing intelligent models and algorithms, comprising: a modeling terminal, a background service and a human-computer interface, wherein starting from business scenarios, the modeling terminal is the core part of the whole platform, and integrates mathematical models in the form of modules; in consideration of obvious difference in features of heterogeneous energy flow in the industrial production process, a general intelligent module for multiple production scenarios, multiple response scales, multiple display modes and multiple working conditions switching is provided; complex models can be built and trained by logic splicing based on services by means of drag and drop; meanwhile, an adaptive algorithm is embedded in each module to dynamically optimize parameters, thus ensuring the adaptive and generalization abilities of models; then, the built models are saved as separate “*.iail” files for repeated invocation; and the modeling terminal is divided into a data pre-processing module, an energy generation-consumption-storage trend forecasting module, an optimal scheduling decision module, a parameter optimization configuration module and a result display module, and each module comprises a plurality of models for different working condition scenarios; the modeling terminal supports the following two interactive modes: (1) on-line model computing: this interactive mode realizes the integration of model training and computing processes; in this interactive mode, a user gets access to the original data through an interface, then automatically completes model training and parameter optimization according to a built logic, and displays the results through the figure/table module; and the model is not saved, and will be retrained in the next computing; (2) off-line model computing: in this interactive mode, model training and computing are two independent processes; and in the real-time computing process, an off-line model obtained by training of historical data is invoked, the parameters of the model are fixed, and the computing process is not updated; the background service is respectively deployed on a cloud service and a edge side service; the edge side service performs real-time computing by reading and serializing information including industrial operation data and production plans/signals, then building training samples according to the standardized format and invoking the saved “*.iail” model file interface, and saves the results in a relational database of the cloud service; and the cloud service receives the foreground trigger information and the working conditions switching signals obtained by analysis of the edge side, so as to trigger the edge side model for self-updating to realize the adaptive adjustment of parameters of the model; the human-computer interface realizes an interface display of the computed result of the model; the interface is refreshed every minute to update the computed result in time; and in the case of temporary adjustment of working conditions, the temporary adjustment of working conditions is fed back to the cloud through human-machine interaction, and the cloud service will trigger the adaptive optimization of the parameters of the edge side model, so as to better meet the actual demand; a platform-based modeling terminal builds a filtering model, a filling model, a forecasting model and a scheduling model by visual drag and drop, and saves the filtering model, the filling model, the forecasting model and the scheduling model as intelligent model files; wherein the edge side service reads actual industrial data, builds training samples according to the standardized format, invokes a model file interface, performs real-time computing based on information including real-time data and production plans/signals, and saves the results in the relational database; wherein an outlier detection method is as follows: production data have aperiodic characteristics or periodic characteristics; for production data with aperiodic characteristics, an adaptive fuzzy C-means (AFCM) clustering method is for outlier detection, which can effectively identify outliers and improve the clustering quality by introducing adaptive vectors and adaptive indexes to adjust and control clustering, and steps of the outlier detection method are as follows: step 1: performing AFCM clustering on a time series to be detected to obtain c.sup.od cluster centers which are defined as v.sub.1.sup.od, . . . , v.sub.c.sub.od.sup.od; step 2: calculating the minimum distance between each data point in the time series and each cluster center; step 3: determining the threshold of an abnormal distance through a box figure, wherein distance values which are more than the threshold are determined as outliers; first, by arranging the data from small to large, and respectively calculating an upper quartile Q.sub.3.sup.od and a lower quartile Q.sub.1.sup.od; then calculating an interquartile range IQR.sup.od=Q.sub.3.sup.od−Q.sub.1.sup.od; and finally, calculating an outlier cut-off point Q.sup.od=Q.sub.3.sup.od+1.5.sub.*IQR.sup.od; the accuracy of the method is judged by a misjudgment ratio η.sub.w.sub.od.sup.od and an omission factor η.sub.ƒ.sub.od.sup.od, of which the formulae are respectively shown in formula (1) and formula (2): $\begin{array}{cc}{\eta }_{w^{\mathrm{od}}}^{\mathrm{od}}=\frac{n_{w^{\mathrm{od}}}^{\mathrm{od}}}{n_{t^{\mathrm{od}}}^{\mathrm{od}}}& \left(1\right)\end{array}$ $\begin{array}{cc}{\eta }_{f^{\mathrm{od}}}^{\mathrm{od}}=\frac{n_{f^{\mathrm{od}}}^{\mathrm{od}}}{n_{w_r^{\mathrm{od}}}^{\mathrm{od}}}& \left(2\right)\end{array}$ wherein n.sub.w.sub.od.sup.od and n.sub.t.sub.od.sup.od are respectively the number of misjudgments and the number of detection results, and n.sub.ƒ.sub.od.sup.od and n.sub.w.sub.t.sub.od.sup.od are respectively the number of omissions and the number of actual outliers; for production data with periodic characteristics, a method based on dynamic time warping is used to detect outliers; and the steps of the outlier detection method are as follows: step 1: dividing historical data samples into n.sup.od sample series [x.sub.1.sup.od, . . . , x.sub.n.sub.od.sup.od] with similar length to the data to be detected according to the periodic characteristics; step 2: calculating the sum SumD.sub.k.sub.od.sup.od of dynamic time warping (DTW) distances between the series X.sub.i.sub.od.sup.od and other n.sup.od−1 series; when SumD.sub.k.sub.od.sup.od=min{SumD.sub.k.sub.od.sup.od, . . . , SumD.sub.n.sub.od.sup.od} for the (k.sup.od).sup.th series, taking the corresponding X.sub.k.sub.od.sup.od, as a center sample, which is recorded as Y.sup.od; letting the time warping function be C.sup.od={c.sup.od(1), c.sup.od(2), . . . , c.sup.od(N)}, where N represents the path length, c.sup.od(n.sub.x.sub.i.sup.od)=(i.sup.od, j.sup.od) represents the (n.sub.x.sub.i.sup.od).sup.th matching point pair composed of the feature vector x.sub.i.sup.od of the series X.sup.od and the (j.sup.od).sup.th feature vector y.sub.j.sub.od of the series Y.sup.od, and the distance d.sup.od(i.sup.od, j.sup.od) between the two is the local matching distance; and the DTW algorithm minimizes the sum of distances through local optimization, as shown in formula (3): $\begin{array}{cc}{D}^{\mathrm{od}}\left(m_1^{\mathrm{od}},m_2^{\mathrm{od}}\right)=\min \left\{\underset{n^{\mathrm{od}}=1}{\overset{N^{\mathrm{od}}}{.Math.}}{d}^{\mathrm{od}}\left(x_{i^{\mathrm{od}}}^{\mathrm{od}},y_{j^{\mathrm{od}}}^{\mathrm{od}}\right)\right\}& \left(3\right)\end{array}$ step 3: with the series X.sup.od of a sample to be detected as the reference, scaling the center sample Y.sup.od based on DTW to obtain Y.sub.*.sup.od, wherein the scaling process meets formula (4): $\begin{array}{cc}{Y}_{*}^{\mathrm{od}}\left[j^{\mathrm{od}}\right]=\frac{1}{.Math.\left\{k^{\mathrm{od}}.Math.x_{i^{\mathrm{od}}}^{\mathrm{od}}.fwdarw.y_{j^{\mathrm{od}}}^{\mathrm{od}}\right\}.Math.}\underset{x_{i^{\mathrm{od}}}^{\mathrm{od}}.fwdarw.y_{j^{\mathrm{od}}}^{\mathrm{od}}}{.Math.}{y}_{j^{\mathrm{od}}}^{\mathrm{od}}& \left(4\right)\end{array}$ wherein k.sup.od represents the number of points corresponding to the series Y.sup.od and the series X.sup.od in the DTW process, and the symbol x.sub.j.sub.od.sup.od represents that the point x.sub.j.sub.od.sup.od in the series X.sup.od corresponds to the point y.sub.j.sub.od.sup.od in the series Y.sup.od; step 4: performing AFCM clustering on the scaled series Y.sub.*.sup.od to obtain a cluster center v.sub.*.sup.od, and calculating the minimum distance between each point in the sample to be detected and all cluster centers; and step 5: obtaining a threshold for determining an abnormal distance by calculation through a box figure; wherein the distance values which are more than the threshold are determined as outliers.

2. The distributed industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein the filtering model is as follows: for real-time data of different working condition scenarios, different intelligent methods are integrated into a filtering module; (1) for working condition scenarios with discrete features in the production process, the feature data are obvious, and a moving average filtering method is adopted; (2) for working condition scenarios with continuous production features, the data generally contain white noise, and the acquired data are processed by the empirical mode decomposition method.

3. The distributed industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein a missing value fill method is as follows: the basic features are analyzed based on the historical sample data and divided into different working condition scenarios, and different filling methods are respectively proposed for discrete single-point missing and continuous multi-point missing; case 1, for discrete single-point missing, filling is divided into two conditions according to the data standard deviation; where under the working condition where the standard deviation is less than ε, data are filled by means of moving average; and under the working condition where the standard deviation is more than ε, data are filled by means of cubic spline interpolation; case 2, for continuous missing in data, an autoregressive integrated moving average model (ARIMA) is used to fill segments; the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are taken as evaluation indexes of the autoregressive integrated moving average model, and a (p.sup.tb, q.sup.tb) combination minimizing AIC or BIC is taken as the optimum parameter of the model where p.sup.tb and q.sup.tb respectively represent the two parameters of the autoregressive moving average model, p.sup.tb is the number of autoregressive terms, and q.sup.tb is the number of moving average terms; then the complete data segments before and after the missing segment are respectively taken as training samples, the data of the missing segment are forecast based on the ARIMA model, and the average value of two forecast results is taken as the fill value to ensure full utilization of the sample information so as to improve the filling accuracy; and steps of the missing value fill method are as follows: step 1: judging the information features of time series data to be filled comprising the miss rate of data, the missing situations of missing points, the periodicity of time series and the standard deviation of time series; if all the missing points in the time series are discrete missing, executing step 2; where if all the missing points in the time series are continuous missing, executing steps 4-5; and if discrete missing points and continuous missing points both exist in the time series, executing step 2 first, and then executing steps 4-5 repeatedly for each continuous missing data segment; step 2: if the standard deviation of the series is less than ε, filling discrete points by means of moving average; otherwise, executing step 3; step 3: carrying out filling by means of cubic spline interpolation; step 4: filling data based on the autoregressive integrated moving average (ARIMA) model; and selecting complete series before and after the missing data segment as samples to forecast the data of the missing segment respectively, so as to calculate the average value of the two forecast results as the fill value; and step 5: establishing an ARIMA model to forecast the missing data segment, and taking the average value of the two forecast results as the fill value.

4. The distributed industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein an energy trend forecast method is as follows: an energy generation-consumption forecasting for class period feature working conditions comprises the following specific steps: step 1: dividing historical sample data into a training set and a validation set, wherein the training set is used for training the model and the validation set is used for validating the model accuracy; step 2: for an echo state network (ESN) method, optimizing the embedded dimensions of samples and a number of nodes in a reservoir of an ESN model by the particle swarm optimization (PSO) algorithm, the number of nodes being a parameter of the reservoir; step 3: for a least square support vector machine (LSSVM) method, optimizing the embedded dimensions of samples and the intrinsic parameters (penalty coefficient and kernel parameter) of the algorithm by the PSO algorithm; and step 4: evaluating the ESN and LSSVM methods based on the training accuracy, and taking the computed result of the method with higher accuracy as the forecast result of the model and direct modeling as follows: step a: initial filling and data set building based on phase space reconstruction; extracting time series data of a period of time from the real-time relational database of an industrial field, detecting and recording positions of missing points in the time series, and filling all missing points with zero values; wherein the function of phase space reconstruction is to build a training sample set; and a discrete time series {x.sub.1.sup.sl, x.sub.2.sup.sl, . . . , x.sub.N.sub.sl.sup.s1}, wherein x.sub.i.sub.sl.sup.sl represents the value acquired at the (i.sup.sl).sup.th sampling moment, and i.sup.sl=1, 2, . . . , N.sup.sl; wherein the (j.sup.sl).sup.th phase vector x.sub.j.sub.sl.sup.sl is shown in formula (5):
X.sub.j.sub.sl.sup.sl=(x.sub.j.sub.sl.sup.sl,x.sub.j.sub.sl.sub.+1.sup.sl, . . . ,x.sub.j.sub.sl.sub.+p.sub.sl.sub.-1.sup.sl).sup.T  (5) wherein p.sup.sl is an embedded dimension of the phase vector, and the matrix form of N.sup.sl input samples is X.sup.sl=[x.sub.1.sup.sl, x.sub.2.sup.s1, . . . , x.sub.N.sub.sl.sup.sl].sup.T; step b: building of relevance vector machine (RVM) model based on incomplete data set; building an initial RVM model, by building a functional relationship between input and output variables; and wherein the mapping relationship between the (i.sup.sl).sup.th function value y.sub.i.sub.sl.sup.sl∈R and an input vector x.sub.j.sub.sl.sup.sl as shown in formula (6): $\begin{array}{cc}{x}_{i^{\mathrm{sl}}}^{\mathrm{sl}}=\underset{j^{\mathrm{sl}}=1}{\overset{N^{\mathrm{sl}}}{.Math.}}{w}_{j^{\mathrm{sl}}}^{\mathrm{sl}}K\left(x_{i^{\mathrm{sl}}}^{\mathrm{sl}},x_{j^{\mathrm{sl}}}^{\mathrm{sl}}\right)+w_0^{\mathrm{sl}}+{\epsilon }_{i^{\mathrm{sl}}}^{\mathrm{sl}}={\left({\varphi }_{i^{\mathrm{sl}}}^{\mathrm{sl}}\right)}^T{w}^{\mathrm{sl}}+{\epsilon }_{i^{\mathrm{sl}}}^{\mathrm{sl}}& \left(6\right)\end{array}$ wherein ϕ.sub.i.sub.sl.sup.sl=[1, K(x.sub.i.sub.sl.sup.sl, x.sub.1.sup.sl), K(x.sub.i.sub.sl.sup.sl, x.sub.2.sup.sl), . . . , K(x.sub.i.sub.sl.sup.sl, x.sub.N.sub.sl.sup.sl)].sup.T and i.sup.sl=1, 2, . . . , N.sup.sl; w.sup.sl=[w.sub.0.sup.sl, w.sub.1.sup.sl, . . . , w.sub.N.sub.sl.sup.sl].sup.T represents a weight vector, and K(x.sub.i.sub.sl.sup.sl, x.sub.j.sub.sl.sup.sl) is a kernel function which is a Gaussian kernel function; ε.sub.i.sub.sl.sup.sl represents a Gaussian distribution noise term with the mean value of 0 and the variance of σ.sup.2; and the weight vector w.sup.sl is considered as a continuous random variable; and each element obeys normal distribution with the mean value of 0 and different variances, as shown in formula (7): $\begin{array}{cc}p\left(w^{\mathrm{sl}}.Math.\alpha \right)=\underset{i^{\mathrm{sl}}=0}{\overset{N^{\mathrm{sl}}}{.Math.}}N\left(w_{i^{\mathrm{sl}}}^{\mathrm{sl}}.Math.0,{\alpha }_i^{-1}\right)& \left(7\right)\end{array}$ wherein α=[α.sub.0, α.sub.1, . . . , α.sub.N].sup.T, and A=diag(α); parameters α and β are considered as random variables which obey Gamma distribution; a noise variance vector is written as β=[β.sub.[m], β.sub.[m], . . . , β.sub.[m], β.sub.[obs], β.sub.[obs]].sup.T, and β.sub.[obs]≡σ.sup.−2; in addition, β also can be written as β=[β.sub.[m].sup.T, β.sub.[obs]].sup.T in the form of diagonal matrices B.sub.[m]=diag(β.sub.[m]), B.sub.[obs]=diag(β.sub.[obs]) and B=diag(β); a training output set is composed of missing points and observable points, wherein x.sup.sl=[(x.sub.[m].sup.sl).sup.T, (x.sub.[obs].sup.sl).sup.T].sup.T, x.sup.sl.sub.[m] and x.sup.sl.sub.[obs] respectively represent the vectors of missing points and observable points in the output set, the number of output missing points and that of output observable points are respectively N.sup.sl.sub.ymiss and N.sup.sl.sub.yobs, and N.sup.sl=N.sup.sl.sub.ymiss+N.sup.sl.sub.yobs; the element β.sub.[m] in β.sub.[m] is defined as 10.sup.8 to ensure that the noise variance of the corresponding missing output point is not affected; and the likelihood function is shown in formula (8): $\begin{array}{cc}p\left(x_{\left[m\right]}^{\mathrm{sl}},x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}.Math.w,\beta \right)=N\left(x^{\mathrm{sl}}.Math.{\Phi }^{\mathrm{sl}}{w}^{\mathrm{sl}},B^{-1}\right)=1/\left[{\left(2\pi \right)}^{N^{\mathrm{sl}}/2}{.Math.{\left(B^{\mathrm{sl}}\right)}^{-1}.Math.}^{1/2}\right]\times \exp \left\{-{\frac{1}{2}\left[x^{\mathrm{sl}}-\Phi w^{\mathrm{sl}}\right]}^{T}B\left[x^{\mathrm{sl}}-\Phi w^{\mathrm{sl}}\right]\right\}& \left(8\right)\end{array}$ wherein Φ.sup.sl=[ϕ.sup.sl.sub.1, ϕ.sup.sl.sub.2, . . . , ϕ.sup.sl.sub.N.sub.sl].sup.T is a kernel function matrix; aiming at the optimization problem of a missing input point x.sup.sl.sub.[m], and considering the characteristics of the sample set built by phase space reconstruction, each missing input point has a corresponding equivalent missing point in the missing output set; and therefore, the posterior probability of each element in x.sup.sl.sub.[m] is replaced with the posterior probability of the corresponding missing output point, the posterior average of x.sup.sl.sub.[m] is filled in the corresponding position of the missing input point, and the kernel matrix Φ.sup.sl is further updated, as shown in formula (9): $\begin{array}{cc}{\mu }_{x_{\left[m\right]},i_{x_{\left[m\right]}}^{\mathrm{sl}}}=\int x_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}p\left(x_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}.Math.x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}},\alpha ,{\beta }_{\left[\mathrm{obs}\right]}\right){\mathrm{dx}}_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}& \left(9\right)\end{array}$ wherein i.sub.x.sub.[m].sup.sl=1, 2, . . . , N.sup.sl.sub.x.sub.[m], the set x.sup.sl.sub.[m] contains N.sup.sl.sub.x.sub.[m] missing input points, x.sub.[m],i.sub.sl.sup.sl represents the (i.sup.sl).sup.th element in x.sup.sl.sub.[m], and μ.sub.x.sub.[m].sub.,i.sub.sl represents the mean posterior probability of x.sub.[m],i.sub.sl.sup.sl; x.sup.sl.sub.[m],corr_i represents the corresponding variable of the (i.sup.sl).sup.th missing input point in the missing output set x.sup.sl.sub.[m]; and by eliminating the integrals of w.sup.sl and x.sup.sl.sub.[m], a log edge likelihood function is given, as shown in formula (10): $\begin{array}{cc}\tilde{L}=\log \int p\left(x_{\left[m\right]}^{\mathrm{sl}},x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}.Math.w^{\mathrm{sl}},{\beta }_{\left[\mathrm{obs}\right]}\right)p\left(w^{\mathrm{sl}}.Math.\alpha \right){\mathrm{dx}}_{\left[m\right]}^{\mathrm{sl}}{\mathrm{dw}}^{\mathrm{sl}}=-\frac{1}{2}\left(N^{\mathrm{sl}}\log 2\pi +\log .Math.\tilde{C}.Math.+{\left(x^{\mathrm{sl}}\right)}_{\left[\mathrm{obs}\right]}^T{\tilde{C}}^{-1}{x}_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}\right)& \left(10\right)\end{array}$ wherein a covariance matrix {tilde over (C)}=Φ.sub.x.sub.[obs]A.sup.−1Φ.sub.x.sub.[obs].sup.T+β.sub.[obs].sup.−1I, and I is a unit matrix; and the extreme values α* and β*.sub.[obs] of α and β.sub.[obs] are found by maximizing partial lead {tilde over (L)}, and the iterative formulae are obtained as shown in formula (11) and formula (12): $\begin{array}{cc}{\alpha }_{i^{\mathrm{sl}}}=\frac{{\tilde{\gamma }}_{i^{\mathrm{sl}}}}{{\tilde{\mu }}_{i^{\mathrm{sl}}}^2}& \left(11\right)\end{array}$ $\begin{array}{cc}{\beta }_{\left[\mathrm{obs}\right]}^{-1}=\frac{{.Math.x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}-{\Phi }_{x_{\left[\mathrm{obs}\right]}}\tilde{\mu }.Math.}^2}{N_{\mathrm{yobs}}^{\mathrm{sl}}-{\beta }_{\left[\mathrm{obs}\right]}\mathrm{Tr}\left(\widetilde{.Math.}{\Phi }_{x_{\left[\mathrm{obs}\right]}}^T{\Phi }_{x_{\left[\mathrm{obs}\right]}}\right)}& \left(12\right)\end{array}$ wherein i.sup.sl=1, 2, . . . , N.sup.sl, {tilde over (γ)}.sub.i=1−α.sub.i{tilde over (Σ)}.sub.ii, {tilde over (Σ)}=(A+β.sub.[obs]Φ.sub.x.sub.[obs].sup.TΦ.sub.x.sub.[obs]).sup.−1 and {tilde over (μ)}=β.sub.[obs]{tilde over (Σ)}Φ.sub.x.sub.[obs].sup.Tx.sup.sl.sub.[obs]; and the method of the RVM model comprises the following steps: {circle around (1)} initializing: α=[1/(N.sup.sl).sup.2, 1/(N.sup.sl).sup.2, . . . , 1/(N.sup.sl).sup.2].sup.T, and β.sub.[obs] is initialized as the variance of 10 times x.sup.sl.sub.[obs]; and the related vector set is initialized as a set of all input samples; {circle around (2)} updating mean values and covariance matrices of joint posterior probabilities of w.sup.sl and x.sup.sl.sub.[m]; {circle around (3)} updating α and β.sub.[obs]; {circle around (4)} if α.sub.j.sub.sl>10.sup.8, removing the (j.sup.sl).sup.th sample from the related vector set; {circle around (5)} calculating the mean posterior probability of each element in x.sup.sl.sub.[m], filling in the corresponding position of the missing input point, and further updating the kernel matrix Φ; {circle around (6)} calculating the value of the log edge likelihood function; {circle around (7)} if the change of an element in α relative to the last iteration or the change of the value of the edge likelihood function is less than the given threshold, stopping the iterative process, thus obtaining a final forecasting model; otherwise, executing step {circle around (2)}; step c: obtaining the forecast value of a future period of time based on RVM single-step iterative forecasting; when a new input x.sup.sl, is given, the value of β.sub.[obs] is replaced with β*.sub.[obs]; and the forecast distribution of the corresponding output x.sup.sl, is approximately normal distribution, as shown in formula (13):
p(t.sup.sl.sub.*|t.sup.sl.sub.[obs])≈∫p(t.sup.sl.sub.*|x.sup.sl.sub.*,w,β*.sub.[obs])p(w|t.sup.sl.sub.[obs],α,β*.sub.[obs])dw=N(t.sup.sl.sub.*|μ.sub.*,σ.sub.*.sup.2)  (13) wherein μ.sub.*=ϕ.sub.*.sup.Tμ.sub.[w|obs], and σ.sub.*.sup.2=(β*.sub.[obs]).sup.−1+ϕ.sub.*.sup.T{tilde over (Σ)}ϕ.sub.* and ϕ.sub.*=[1, K(x.sup.sl.sub.*, x.sup.sl.sub.1), K(x.sup.sl.sub.*, x.sup.sl.sub.2), . . . , K(x.sup.sl.sub.*, x.sup.sl.sub.N.sub.sl)].sup.T; μ.sub.[w|obs] is the mean posterior probability of the weight vector w.sup.sl; and then, the forecast mean value of a future period of time is a point, and the next round of forecasting will take the forecast result of the previous round as input to carry out cyclic rolling forecast; step d: building of LSSVM data set based on incomplete data set; after the parameters α and b are solved, using x.sup.sl.sub.*=αK(x.sup.sl, x.sup.sl.sub.*)+b to calculate an estimated value of the forecast quantity; therefore, inputting data containing missing points filled by moving average into the LSSVM model for training, then forecasting the data containing missing points, replacing the original values of missing points with the forecast result, and repeating the above process until the change of the values of missing points is less than the set threshold, thus completing the data updating process; step e: obtaining the forecast value of a future period of time based on Gaussian process (GP) forecasting; after updating the missing points by LSSVM, inputting the updated data into the GP model, wherein GP is a probability model with kernel; and where the weight w.sup.sl and the noise ε.sup.sl in formula (14) obey zero-mean prior distribution:
p(w.sup.sl)=N(w.sup.sl|0,α.sup.−1I)
p(ε.sup.sl)=N(ε.sup.sl|0,β.sup.−1)  (14) wherein α and β are hyperparameters of the weight and the noise distribution, and the output likelihood function is expressed as formula (15):
p(x.sup.sl|Φ,w.sup.sl,β)=N(x.sup.sl|Φw.sup.sl,β.sup.−1I)  (15) according to Bayes formula, deducing the posterior distribution of the weight as shown in formula (16): $\begin{array}{cc}p\left(w^{\mathrm{sl}}.Math.x^{\mathrm{sl}},\Phi ,\alpha ,\beta \right)=\frac{p\left(x^{\mathrm{sl}}.Math.\Phi ,w^{\mathrm{sl}}\right)p\left(w^{\mathrm{sl}}\right)}{p\left(x^{\mathrm{sl}}.Math.\Phi \right)}=N\left(w^{\mathrm{sl}}.Math.m_{N^{\mathrm{sl}}},.Math.\right)& \left(16\right)\end{array}$ wherein m.sub.N=βΣΦ.sup.Tx.sup.sl and Σ=(αI+βΦ.sup.TΦ).sup.−1; and then selecting the Gaussian kernel function and maximizing the likelihood function, thus obtaining optimal hyperparameters; step f: algorithm selection; calculating the mean absolute percentage error (MAPE) of the validation sets obtained by the RVM algorithm and the GP algorithm respectively, selecting the algorithm with the lowest MAPE, and saving and outputting the hyperparameter values and the information of the kernel function optimized by the selected algorithm; and step g: forecasting; after real data are acquired, using the selected algorithm and the parameters thereof for forecasting.

5. The industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein the scheduling method for multiple by-product gas systems is as follows: first, according to the capacity operation states of a blast furnace gas system, a coke oven gas system and a linz-donawitz converter gas system, the balance of each system is judged with system safety indexes as constraints to give scheduling suggestions of each single-medium system; second, with the adjustable range of each single-gas medium and the load capacity of each adjustable user as constraints, the comprehensive adjustable range of multiple gas media is calculated; finally, the adjustable range for optimal scheduling is obtained by integrating the adjustable range of a pipe network and the adjustable range of multiple media; and the calculation method for the adjustable range of each stage is as follows: (1) solving the adjustable range of a single medium wherein the adjustment of the pipe network is judged according to the state of the pipe network using is a two-stage judgment method; wherein the first stage is a short-term period, and the second stage is a long-term period; each stage is divided into three parts, the state of each part is judged, and the overall state of each stage is judged according to the state of each part, so as to judge the adjustment of the pipe network; the system adjustable range is obtained by repeatedly judging the state of the pipe network under each adjustment value; when the capacity state under this adjustment value is judged to require no adjustment, the adjustment value at this moment can be adopted; the adjustment value within the adjustment range is judged repeatedly, finally obtaining the adjustable range; and supposing that the use flow of blast furnace gas, coke oven gas and linz-donawitz converter gas is x.sup.s, y.sup.s and z.sup.s respectively, the adjustable ranges is defined as [x.sub.1.sup.s, x.sub.2.sup.s], [y.sub.1.sup.s, y.sub.2.sup.s] and [z.sub.1.sup.s, z.sub.2.sup.s] respectively; (2) solving the adjustable range of multiple media wherein the use range of each medium within the adjustable range under the multi-media adjusting user is calculated as shown in formula (24):
Total_cal=bfg_cal×x+cog_cal×y+ldg_cal×z  (24) wherein Total_cal is the gross calorific value of the multi-media adjusting user, and bfg_cal, cog_cal and ldg_cal are respectively the calorific value of each gas medium; and formula (25) is obtained by the elimination method: $\begin{array}{cc}{x}^s=\frac{{\mathrm{Total}}_{-}\mathrm{cal}-\left(\mathrm{cog}\_\mathrm{cal}\times y+\mathrm{ldg\_cal}\times z\right)}{\mathrm{bfg\_cal}}& \left(25\right)\end{array}$ the value ranges [x.sub.1.sup.s′, x.sub.2.sup.s′], [y.sub.1.sup.s′, y.sub.2.sup.s′] and [z.sub.1.sup.s′, z.sub.2.sup.s′] of x.sup.s, y.sup.s and z.sup.s in the case of x.sup.s∈[y.sub.1.sup.s, y.sub.2.sup.s], z.sup.s∈[z.sub.1.sup.s, z.sub.2.sup.s] and x.sup.s∈[x.sub.1.sup.s, x.sub.2.sup.s] are solved, that is the adjustable range of each gas medium after multi-media adjustment is considered; and then, the final adjustable ranges [x.sub.1.sup.s″, x.sub.2.sup.s″], [y.sub.1.sup.s″, y.sub.2.sup.s″] and [z.sub.1.sup.s″, z.sub.2.sup.s″] for optimal scheduling are obtained through interval differencing, as shown in formulae (26)-(28) respectively:
[x.sub.1.sup.s″,x.sub.2.sup.s″]=[x.sub.1.sup.s,x.sub.2.sup.s]−[x.sub.1.sup.s′,x.sub.2.sup.s′]  (26)
[y.sub.1.sup.s″,y.sub.2.sup.s″]=[y.sub.1.sup.s,y.sub.2.sup.s]−[y.sub.1.sup.s′,y.sub.2.sup.s′]  (27)
[z.sub.1.sup.s″,z.sub.2.sup.s″]=[z.sub.1.sup.s,z.sub.2.sup.s]−[z.sub.1.sup.s″,z.sub.2.sup.s′]  (28).

6. The industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein a method for collaborative optimization and scheduling of comprehensive energy is as follows: considering that the adjustable range of byproduct gas and the total steam demand within a forecasting period t.sub.0-t.sub.1 have been determined based on the aforementioned forecasting model, the energy consumption of combined heat and power sets and boilers is distributed to minimize the energy cost and maximize the comprehensive energy efficiency, and the corresponding objective functions are respectively shown in formula (29) and formula (30): $\begin{array}{cc}\min J_{\mathrm{co}}=\underset{i^{\mathrm{mf}}}{\overset{3}{.Math.}}{q}_{i^{\mathrm{mf}}}^u{X}_{i^{\mathrm{mf}}}-q_{\mathrm{power}}^g{P}_{\mathrm{total}}^g-q_{\mathrm{steam}}^g{F}^g& \left(29\right)\end{array}$
max J.sub.effc(KH.sub.steamF.sup.g+KH.sub.eleP.sub.total.sup.g)/Σ.sub.i.sub.mf.sub.=1.sup.3KH.sub.i.sub.mfX.sub.i.sub.mf  (30) wherein J.sub.co represents energy cost, g represents generating capacity, and u represents consumption; $\underset{i^{\mathrm{mf}}}{\overset{3}{.Math.}}{q}_{i^{\mathrm{mf}}}^u{X}_{i^{\mathrm{mf}}}$ is the sum of consumption costs of byproduct gas and purchased coal, q.sub.i.sub.mf.sup.u is the unit price of the (i.sup.mf).sup.th medium, and i.sup.mf=1, 2, 3, respectively representing BFG, coke oven gas (COG) and coal, X.sub.i.sub.mf a is the corresponding consumption amounts: q.sub.power.sup.g and q.sub.steam.sup.g respectively represent the unit price of power generation and the unit price of steam generation of a heat and power system, and P.sub.total.sup.g and F.sup.g are respectively the total amount of self-generating electricity and the total amount of steam generation; J.sub.effc is the comprehensive energy efficiency of a heat and power set; KH.sub.steam and KH.sub.ele respectively represent conversion coefficients of total steam consumption and total power generation into standard coal; KH.sub.i.sub.mf represents the conversion coefficient of each energy medium; and the corresponding constraint conditions are shown in formula (31): $\begin{array}{cc}\{\begin{array}{cc}{X}_1=1/\left(1-{\delta }_B\right)\underset{i_1=1}{\overset{6}{.Math.}}{F}_{\mathrm{BFG},i_1}^{\mathrm{BLR}}& \left(c1\right)\\ X_2=1/\left(1-{\delta }_C\right)\underset{i_1=1}{\overset{6}{.Math.}}{F}_{\mathrm{COG},i_1}^{\mathrm{BLR}}& \left(c2\right)\\ X_3=1/\left(1-{\delta }_M\right)\underset{i_2=1}{\overset{2}{.Math.}}{F}_{M,i_2}^{\mathrm{BLR}}& \left(c3\right)\\ D_{S,i_2}^{\mathrm{BLR}}=f_{i_2}\left(F_{\mathrm{BFG},i_2}^{\mathrm{BLR}},F_{\mathrm{COG},i_2}^{\mathrm{BLR}},F_{M,i_2}^{\mathrm{BLR}}\right),\left(i_2=1,2\right)& \left(c4\right)\\ D_{S,j_1}^{\mathrm{BLR}}=f_{j_1}\left(F_{\mathrm{BFG},j_1}^{\mathrm{BLR}},F_{\mathrm{COG},j_1}^{\mathrm{BLR}}\right),\left(j_1=3,4,5,6\right)& \left(c5\right)\\ \underset{i_2=1}{\overset{2}{.Math.}}{D}_{S,i_2}^{\mathrm{BLR}}=\underset{j_2=1}{\overset{2}{.Math.}}{D}_{T,j_2}^{\mathrm{ST}}& \left(c6\right)\\ \underset{i_3=1}{\overset{4}{.Math.}}{D}_{\mathrm{HS},i_3}^{\mathrm{BLR}}=\underset{j_3=3}{\overset{4}{.Math.}}{D}_{T,j_3}^{\mathrm{ST}}& \left(c7\right)\\ D_{T,j_2}^{\mathrm{ST}}=D_{P,j_2}^{\mathrm{ST}}+D_{S,j_2}^{\mathrm{ST}},\left(j_2=1,2\right)& \left(c8\right)\\ P_{\mathrm{total}}^g=\underset{i_3=1}{\overset{4}{.Math.}}{P}_{i_3}^{\mathrm{ST}}& \left(c9\right)\\ F_M^u=\left(1-{\delta }_S\right)\left(\underset{i_2=1}{\overset{2}{.Math.}}{D}_{S,i_2}^{\mathrm{ST}}+\underset{j_1=3}{\overset{6}{.Math.}}{D}_{S,j_1}^{\mathrm{BLR}}+F^r\right)& \left(c10\right)\\ \underset{i_3=3}{\overset{4}{.Math.}}{D}_{T,i_3}^{\mathrm{ST}}=\left(1-{\delta }_S\right)\underset{j_1=3}{\overset{6}{.Math.}}{D}_{\mathrm{HS},j_1}^{\mathrm{BLR}}& \left(c11\right)\\ F^g=1/\left(1-{\delta }_S\right)\left(F_M^u+\underset{i_3=3}{\overset{4}{.Math.}}{D}_{T,i_3}^{\mathrm{ST}}\right)& \left(c12\right)\\ {\mathrm{Pa}}^L<{\mathrm{Pa}}_0+K_{\mathrm{pa}}{F}^g<{\mathrm{Pa}}^H& \left(c13\right)\end{array}& \left(31\right)\end{array}$ wherein the constraint conditions in formula (31) can be described as follows: c1-c3: X.sub.1 and X.sub.2 are respectively the total flow of BFG and the total flow of COG input of all the boilers, where F.sub.FLG,i.sub.1.sup.BLR and F.sub.COG,i.sub.1.sup.BLR are the BFG and COG consumption flow of the i.sub.1th boiler, X.sub.3 is the total consumption of coal, and F.sub.M,i.sub.2.sup.BLR is the coal consumption of the i.sub.2th boiler; δ.sub.B, δ.sub.C and δ.sub.M are respectively the corresponding loss rates; c4-c6: the relationship between the input and the output of each device, and ƒ(⋅) is the corresponding fitting relationship; D.sub.S,j.sub.2.sup.BLR is the steam generation flow of the i.sub.2th boiler of the power plant, and D.sub.S,j.sub.1.sup.DBLR represents the steam generation flow of the j.sub.1th boiler which are supplying steam into the network separately; D.sub.T,j.sub.2.sup.ST is the steam consumption flow of the j.sub.2th 300 MW steam turbine; c7-c8: balance of generation and consumption of high pressure steam and medium pressure steam; D.sub.HS,j.sub.3.sup.BTR denotes the generation flow of the high pressure steam by the four steam supplying boilers and D.sub.T,j.sub.3.sup.ST means the consumption flow of a j.sub.3th 25 MW steam turbine; D.sub.P,j.sub.2.sup.ST is the steam that used to generate power of a j.sub.2th 300 MW steam turbine, and D.sub.S,j.sub.2.sup.ST is the corresponding extraction steam; c9: constraints on self-generating capacity, and P.sub.i.sub.3.sup.ST is the power generation of the i.sub.3th steam turbine; c10-c12: balance constrains on steam pipe network, F.sup.r is the total recovery amount of waste heat, and δ.sub.s is the loss rate of steam; D.sub.S,j.sub.1.sup.BLR is the steam generation flow of the j.sub.1th steam supplying boiler, and F.sub.M.sup.u is the total steam consumption flow of the steam users (without the four steam turbines); and c13: constraints on pressure of steam pipe network; Pa.sup.L and Pa.sup.H are respectively the upper limit and the lower limit of the pressure of a steam pipe network in the south area, and Pa.sub.0 represents the initial state of the pressure; K.sub.pa is the conversion coefficient of flow to pressure; then the swarm intelligence optimization method is adopted for solution, thus obtaining a scheduling scheme with optimal economy.

7. The distributed industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein an energy trend forecast method is as follows: energy generation-consumption forecasting for aperiodic feature working conditions comprising the following specific steps: (1) preparing data; extracting time series data of blast furnace gas (BFG) consumption flow for cold rolling of a period of time from the real-time relational database of an industrial field; (2) acquiring variational mode decomposition (VMD) parameters optimized by the Bayesian optimization method and decomposing and reorganizing the data; for the original series ƒ.sup.ns(t), the problem of constraint variation is shown in formula (17): $\begin{array}{cc}\underset{\left\{u_{k^{\mathrm{ns}}}^{\mathrm{ns}},{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\}}{\min }\left\{\underset{k^{\mathrm{ns}}=1}{\overset{K^{\mathrm{ns}}}{.Math.}}{.Math.{\partial }_t\left[\left(\delta \left(t\right)+\frac{j^{\mathrm{ns}}}{\pi t}\right)u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)\right]e^{-j{\omega }_{k}t}.Math.}_2^2\right\}s.t.\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)=f^{\mathrm{ns}}\left(t\right)& \left(17\right)\end{array}$ wherein K.sup.ns is the number of modes, u.sub.k.sub.ns.sup.ns(t) is the (k.sup.ns).sup.th mode, ω.sub.k.sub.ns.sup.ns is the center frequency, and δ(t) represents Dirac distribution; and to solve the constraint variation problem of formula (17), the quadratic penalty term and the Lagrange multiplier Δ.sup.ns are used to transform the formula (17) into an unconstrained problem, as shown in formula (18): $\begin{array}{cc}L\left(\left\{u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\},\left\{{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\},{\lambda }^{\mathrm{ns}}\right)={\alpha }^{\mathrm{ns}}\underset{k^{\mathrm{ns}}}{.Math.}{.Math.{\partial }_t\left[\left(\delta \left(t\right)+\frac{j^{\mathrm{ns}}}{\pi t}\right)*u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)\right]e^{-j{\omega }_{k}t}.Math.}_2^2+{.Math.f^{\mathrm{ns}}\left(t\right)-\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right).Math.}_2^2+.Math.{\lambda }^{\mathrm{ns}}\left(t\right),f^{\mathrm{ns}}\left(t\right)-\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right).Math.& \left(18\right)\end{array}$ in order to solve the mode u.sub.k.sub.ns.sup.ns and the center frequency ω.sub.k.sub.ns.sup.ns, the problem in the formula (18) is transformed into a minimization problem and solved by an alternating direction method of multipliers, and then the optimal solutions can be expressed as formula (19) and formula (20): $\begin{array}{cc}{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)=\frac{{\hat{f}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)-\underset{i^{\mathrm{ns}}\ne k^{\mathrm{ns}}}{.Math.}{\hat{u}}_{i^{\mathrm{ns}}}\left({\omega }^{\mathrm{ns}}\right)+\left({\hat{\lambda }}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)/2\right)}{1+2{{\alpha }^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}-{\omega }_{K^{\mathrm{ns}}}^{\mathrm{ns}}\right)}^2}& \left(19\right)\end{array}$ $\begin{array}{cc}{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}=\frac{{\int }_0^{\infty }{\omega }^{\mathrm{ns}}{.Math.{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right).Math.}^{2}d{\omega }^{\mathrm{ns}}}{{\int }_0^{\infty }{.Math.{\hat{u}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right).Math.}^{2}d{\omega }^{\mathrm{ns}}}& \left(20\right)\end{array}$ wherein {circumflex over (ƒ)}.sub.ns(ω.sup.ns), û.sub.i.sub.ns.sup.ns(ω.sup.ns), {circumflex over (λ)}.sup.ns(ω.sup.ns) and u.sub.k.sub.ns.sup.ns(ω.sup.ns) represent Fourier transform forms of corresponding variables respectively; and after u.sub.k.sub.ns.sup.ns and ω.sub.k.sub.ns.sup.ns are updated, the Lagrange multiplier is also updated according to formula (21): $\begin{array}{cc}{\hat{\lambda }}_{t+1}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)={\hat{\lambda }}_t^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)+\tau \left({\hat{f}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)-\underset{k^{\mathrm{ns}}}{.Math.}{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)\right)& \left(21\right)\end{array}$ then the above updating process is repeated until the convergence condition of formula (22) is satisfied: $\begin{array}{cc}\underset{k^{\mathrm{ns}}}{.Math.}\frac{{.Math.{\hat{u}}_{k^{\mathrm{ns}}}\left(t+1\right)-{\hat{u}}_{k^{\mathrm{ns}}}\left(t\right).Math.}_2^2}{{.Math.{\hat{u}}_{k^{\mathrm{ns}}}\left(t\right).Math.}_2}<{\epsilon }^{\mathrm{ns}}& \left(22\right)\end{array}$ wherein ε.sup.ns is a threshold satisfying the convergence condition; and a Bayesian parameter optimization process comprises two parts; first, the previous function must be selected to represent the hypothesis about the optimized function; Gauss processes are selected as a probability model for Bayesian parameter optimization; and second, an acquisition parameter for building an acquisition function from model posteriori to determine the next parameter combination to be estimated; wherein a set of hyperparameters is X.sup.ns=x.sup.ns.sub.1, x.sup.ns.sub.2, . . . , x.sup.ns.sub.n.sub.ns, which has a functional relationship with the target to be optimized; for ƒ.sup.ns(X.sup.ns) and GP is a Gaussian process model of input data X.sup.ns and an objective function ƒ.sup.ns, and S.sup.ns is an acquisition function; first, time series data are decomposed by complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) to obtain the number K.sup.ns.sub.k.sub.ns of VMD reserved series; for the VMD method, the decomposition number K.sup.ns=K.sup.ns.sub.d.sub.ns+K.sup.ns.sub.k.sub.ns wherein the number K.sup.ns.sub.d.sub.ns of discarded series and the balance factor α.sup.ns are optimized by the Bayesian optimization method; and the objective function of the Bayesian optimization method is expressed as formula (23):
ƒ.sub.obj.sup.ns(K.sub.d.sub.ns.sup.ns,α.sup.ns)=e.sub.origin.sup.ns+e.sub.model.sup.ns  (23) wherein e.sub.origin.sup.ns represents the error between the series composed of first K.sub.k.sub.ns.sup.ns components after decomposition by the VMD method and the original series, and e.sub.model.sup.ns is the error of VMD reorganized series in model training; and after the optimization parameters K.sup.ns and α.sup.ns of VMD are obtained, the original series is decomposed adding and combining the first K.sub.k.sub.ns.sup.ns series into a new time series, and input/output samples are built by using the time series and then divided into a training set and a validation set; (3) training a bidirectional multi-layer long short-term memory (LSTM) model; wherein the parameters of the bidirectional multi-layer LSTM model are as follows: the input dimension, the number of LSTM layers, the number of hidden nodes of each LSTM layer, the learning rate, the number of samples of a training batch, and the parameter optimization method is grid search with early stopping, wherein the input dimension is optimized first, then the number of LSTM layers and the number of hidden nodes of each LSTM layer are optimized, and finally, the learning rate and the number of samples of a training batch are optimized; and the training steps of the LSTM model are as follows: {circle around (1)} determining the optimization range and granularity of the parameters, the optimization direction is decreasing the number of LSTM layers, the number of hidden nodes of each LSTM layer is decreased, and the objective function is the root-mean-square error; {circle around (2)} recording errors between searches during grid search, and stopping search when the error is less than a threshold ò; {circle around (3)} when optimizing the number of LSTM layers and the number of hidden nodes of each LSTM layer, setting the learning rate and the number of samples of a training batch to the maximum values to shorten the search time, and then using the optimized parameters to train the bidirectional multi-layer LSTM model, thus obtaining a forecasting model; (4) carrying out forecasting based on the latest data acquired from the database; taking the latest 1,000 sample points, performing VMD decomposition and reorganization to reduce noise first, and then taking the first input dimension sample points as the input to forecast points at the next moment; for the aperiodic feature time series containing missing points, the missing data are filled by cubic spline interpolation first, as shown in formula (15) and formula (16), and then a forecasting method for complete aperiodic feature data is adopted to forecast such data.

8. The distributed industrial energy resource operation optimization platform automatically building intelligent models and algorithms according to claim 1, wherein an energy trend forecast method is as follows: forecasting of gas reserves for factor model comprising the following steps: step 1: pre-processing data; extracting the input flow of a plurality of blast furnace gas pipe networks of a period of time and the consumption flow of each gas consumption link from the real-time relational database of an industrial field as input factors of the model, and adding the sum of the capacities at the previous moment to the input factors; and supposing that the sum of capacities at the current moment is the output value; meanwhile, normalizing all the input data and the output data; and selecting the first 70% of samples as training samples and the remaining 30% as validation samples; step 2: building a relevance vector machine (RVM) regressive model; step 3: building a Gaussian process (GP) regressive model; step 4: building an LSSVM model; and step 5: after the output estimated values of the validation sets of the RVM regressive model, the GP regressive model and the LSSVM model are respectively obtained by using the validation sets, comparing the model performance by using the error evaluation index MAPE, so as to obtain the optimal model; wherein the smaller the MAPE value is, the more accurate the forecasting effect of the model is; selecting the optimal model with lowest MAPE, and forecasting the estimated value of a forecast point of a future period of time.

Description

DESCRIPTION OF DRAWINGS

(1) FIG. 1 shows a distributed operation architecture of an industrial energy resource operation optimization platform.

(2) FIG. 2 is an operation architecture diagram of interaction of a background service on an edge side, a human-computer interface and an intelligent model.

(3) FIG. 3 is a diagram showing on-line and off-line interactive modes of a modeling terminal.

(4) FIG. 4 is a structural diagram of an integrated energy system of a large domestic steel enterprise.

(5) FIG. 5 is a data feature map of working condition scenarios of a square wave.

(6) FIG. 6(a) is an envelope diagram of EMD filtering of oxygen flow of a blast furnace.

(7) FIG. 6(b) is an IMF exploded view of EMD filtering of oxygen flow of a blast furnace.

(8) FIG. 7 shows an EMD filtering effect of oxygen flow of a blast furnace.

(9) FIG. 8 is a trend chart of typical aperiodic feature data.

(10) FIG. 9 is a diagram of outlier detection results of aperiodic feature data.

(11) FIG. 10 is a trend chart of typical class period feature data.

(12) FIG. 11 is a diagram of outlier detection results of class period feature data.

(13) FIG. 12 is a trend chart of discrete missing data.

(14) FIG. 13 is a trend chart of continuous missing data.

(15) FIG. 14(a) is a comparison diagram of filling effects in the case of discrete missing with the miss rate of 5%.

(16) FIG. 14(b) is a comparison diagram of filling effects in the case of discrete missing with the miss rate of 15%.

(17) FIG. 14(c) is a comparison diagram of filling effects in the case of continuous missing with the miss rate of 5%.

(18) FIG. 14(d) is a comparison diagram of filling effects in the case of discrete missing with the miss rate of 15%.

(19) FIG. 15 is a comparison diagram of forecast effects of complete data of class period feature working conditions.

(20) FIG. 16(a) is a schematic diagram of a time series containing missing points.

(21) FIG. 16(b) is a schematic diagram of a built input/output training sample set.

(22) FIG. 17(a) is a comparison diagram of forecast results under working conditions with the miss rate of 5%.

(23) FIG. 17(b) is a comparison diagram of forecast results under working conditions with the miss rate of 10%.

(24) FIG. 17(c) is a comparison diagram of forecast results under working conditions with the miss rate of 30%.

(25) FIG. 18 is a comparison diagram of forecast effects of complete data of aperiodic feature working conditions.

(26) FIG. 19 is a comparison diagram of forecast results under working conditions with the miss rate of 15%.

(27) FIG. 20 is a comparison diagram of forecast results of gas reserves of a factor model.

(28) FIG. 21 is a structural diagram of solving an adjustable range of a single gas medium.

(29) FIG. 22(a) shows an adjustable range of the capacity of a BFG system.

(30) FIG. 22(b) shows an adjustable range of the capacity of a COG system.

(31) FIG. 22(c) shows an adjustable range of the capacity of an LDG system.

(32) FIG. 23 is a structural diagram of a heat and power system.

(33) FIG. 24(a) is a trend chart of the capacity of a BFG system after scheduling.

(34) FIG. 24(b) is a trend chart of the capacity of a COG system after scheduling.

(35) FIG. 24(c) is a trend chart of the capacity of an LDG system after scheduling.

DETAILED DESCRIPTION

(36) To better understand the technical route and implementation solution of the present invention, with data pre-processing, generation-consumption-storage trend forecasting and multi-energy flow optimal scheduling processes of an energy system of a large domestic steel enterprise as an example, the algorithm of the modeling terminal and the hardware design of an embedded industrial controller of the background service on the edge side of the platform are described in detail. The energy system comprises multiple energy flows such as steam, gas and power, as shown in FIG. 4. As important secondary energy for steel enterprises, byproduct gas can replace primary energy (natural gas) in many links due to the advantages of high calorific value and low cost. Therefore, accurate forecasting and optimal scheduling of the whole process of generation-consumption-storage are of great significance for enterprises to conserve energy, reduce emissions, reduce costs and increase efficiency.

(37) (1) Modeling terminal: data pre-processing, energy generation-consumption-storage trend forecasting and multi-energy flow optimal scheduling model

(38) The historical data of the capacity, the total generated flow and the total consumption flow of the gas system are subjected to filtering, missing value fill and outlier correction by a data pre-processing module, and trained as the input of a forecasting model. The parameters are updated adaptively, so as to obtain a forecasting model with the corresponding capacity. Furthermore, with gas capacity safety, steam demand and adjustable user load capacity as constraint conditions and with economy as the target, an optimal scheduling model of an integrated energy system can be built by logic splicing of the modules, and a reasonable scheduling proposal can be given. The following is the sub-working condition processing algorithm of each module required to build the optimal scheduling model of an integrated energy system.

(39) 1) Data Pre-Processing

(40) Filtering

(41) Considering that process data such as energy and production actually acquired in industry generally contain noise, it is necessary to filter the data first to reduce the influence of noise on data forecasting and other results. For real-time data (such as square wave data and data containing white noise) of different working condition scenarios, different intelligent methods are integrated into a filtering module to achieve targeted filtering for different feature data and improve the filtering reliability.

(42) For working condition scenarios with square wave features in the production process, the feature data are obvious. The BFG consumption flow of a blast furnace hot blast stove shown in FIG. 5 is a typical square wave feature scenario. The present invention adopts a moving average filtering method. For the time series X.sup.f with M.sup.f samples, the formula is shown in formula (1):

(43) $\begin{array}{cc}\stackrel{\_}{X}=\frac{1}{N^f}\underset{i^f=n^f-N^f+1}{\overset{n^f}{.Math.}}{X}_{i^f}^f,n^f=N^f,N^f+1,\mathrm{.Math.},M^f& \left(1\right)\end{array}$
wherein X is the average value of sampling values at the past N.sup.f (N.sup.f<M.sup.f) moments, and X.sup.f, is a sampling value at the (i.sup.f).sup.th sampling moment. A temporary storage area for N.sup.f pieces of data is correspondingly opened in a storage area. Each newly sampled datum is saved in the temporary storage area, and the oldest datum is removed at the same time, so as to continuously update and iterate sampled data and realize filtering.

(44) For working condition scenarios with continuous production features, the data generally contain white noise. The present invention adopts an empirical mode decomposition (EMD) method to process the acquired data. In the method, selective reconstruction is carried out on the intrinsic mode function (IMF) to realize denoising and provide support for improving forecasting accuracy. The calculation steps are as follows:

(45) Step 1: calculating all extreme points of a signal series. Respectively connecting local maximum points and minimum points through cubic spline curves to form an upper envelope and a lower envelope, and calculating a new signal h.sub.1.sup.1(t) with low frequency signals removed, as shown in formula (2):
h.sub.1.sup.1(t)=x.sup.fc(t)−m.sub.1.sup.fc(t)  (2)
wherein m.sub.1.sup.fc(t) is the average value of the upper envelope and the lower envelope at the moment t, and x.sup.fc(t) is the original acquisition series. During the decomposition process, IMF is required to satisfy two conditions, i.e., on the whole signal, the difference between the number of extreme points and the number of zero crossing points is not more than 1, and the average value of the upper envelope and the lower envelope at any point is 0. If h.sub.1.sup.1(t) is satisfactory, h.sub.1.sup.1(t) can be considered as the first IMF component IMF.sub.1(t) of x.sup.fc(t), or h.sub.1.sup.1(t) can be used as original data for repeated calculation.

(46) Step 2: supposing the (n.sup.fc).sup.th IMF component IMF.sub.n.sub.fc(t). When the margin r.sub.n.sub.fc.sup.fc(t) is less than the preset value and is a monotonic function or a constant, the calculation ends. At this moment, the original signal can be expressed in the form of formula (3):

(47) $\begin{array}{cc}{x}^{\mathrm{fc}}\left(t\right)=\underset{i^{\mathrm{fc}}=1}{\overset{n^{\mathrm{fc}}}{.Math.}}{\mathrm{IMF}}_{i^{\mathrm{fc}}}\left(t\right)+r_{n^{\mathrm{fc}}}^{\mathrm{fc}}\left(t\right)& \left(3\right)\end{array}$
wherein r.sub.n.sub.fc.sup.fc(t) is a signal average trend or a trend item of the average value. According to the continuous mean-square error criterion, an appropriate order k.sup.fc is selected to reconstruct the signal, and finally, the filtered reconstructed signal {tilde over (x)}.sup.fc(t) can be obtained, as shown in formula (4):

(48) $\begin{array}{cc}{\tilde{x}}^{\mathrm{fc}}\left(t\right)=\underset{j^{\mathrm{fc}}=k^{\mathrm{fc}}+1}{\overset{n^{\mathrm{fc}}}{.Math.}}{\mathrm{IMF}}_{j^{\mathrm{fc}}}\left(t\right)& \left(4\right)\end{array}$

(49) The data of oxygen flow of a blast furnace of a steel enterprise are adopted to conduct a contrast experiment. FIG. 6(a) and FIG. 6(b) are respectively an envelope diagram and an IMF exploded view of EMD filtering of oxygen flow of a blast furnace. The original signal is processed and decomposed into ten IMF components and residual components. Each IMF component is shown in FIG. 6(b). FIG. 7 is an actual effect diagram after filtering with EMD in combination with a continuous mean square error. Results show that outliers and noises in the data are filtered out after being filtered by the platform, which can provide a good data basis for the subsequent energy forecasting and optimal scheduling.

(50) Outlier Detection

(51) Due to long transmission distance, complex pipe network and violent fluctuation in generation and consumption of the industrial energy system, data exceptions easily occur under the influence of environment and sensor in the data acquisition process. The present invention builds corresponding outlier detection models respectively for working condition scenarios with class period and aperiodic characteristics of the data.

(52) Typical working conditions with aperiodic characteristics of production data (such as input flow of first-stage COG) are shown in FIG. 8. An adaptive fuzzy C-means (AFCM) clustering method is proposed to detect outliers of such data. The method can effectively identify outliers and improve the clustering quality by introducing adaptive vectors and adaptive indexes to adjust and control clustering. AFCM is an improvement on the basis of FCM, and the objective function of AFCM is expressed by J.sub.AFCM.sup.od, as shown in formula (5):

(53) $\begin{array}{cc}\{\begin{array}{c}{J}_{\mathrm{AFCM}}^{\mathrm{od}}=\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{\left(w_{j^{\mathrm{od}}}^{\mathrm{od}}\right)}^{p^{\mathrm{od}}}\underset{i^{\mathrm{od}}=1}{\overset{c^{\mathrm{od}}}{.Math.}}{\left(u_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^{m^{\mathrm{od}}}{\left(d_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^2\\ s.t.\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{u}_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}>0,\underset{i^{\mathrm{od}}=1}{\overset{c^{\mathrm{od}}}{.Math.}}{u}_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}=1,\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{w}_{j^{\mathrm{od}}}^{\mathrm{od}}=1\end{array}& \left(5\right)\end{array}$
wherein w.sub.ƒ.sub.od.sup.od represents the adaptivity of a sample point x.sub.j.sub.od.sup.od, i.e., the influence of x.sub.j.sub.od.sup.od on the result of the objective function J.sub.AFCM.sup.od. p.sup.od is an adaptive index affecting the clustering quality and the convergence rate, and the value range is generally 1≤|p.sup.od|≤10, wherein p.sup.od≠0. The constraint condition

(54) $\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{w}_{j^{\mathrm{od}}}^{\mathrm{od}}$
reflects the internal relation between sample data points. u.sub.i.sub.od.sub.j.sub.od.sup.od is the subjection degree of the adaptive model, m.sup.od and c.sup.od are respectively a fuzzy coefficient and a clustering number, and d.sub.i.sub.od.sub.j.sub.od.sup.od represents the distance between a sample data point x.sub.j.sub.od.sup.od and a cluster center v.sub.i.sub.od.sup.od. u.sub.i.sub.od.sub.j.sub.od.sup.od, v.sub.i.sub.od.sup.od and w.sub.j.sub.od.sup.od can be obtained respectively by formula (6) to formula (8):

(55) $\begin{array}{cc}{u}_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}=\frac{{\left(d_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^{-\frac{2}{m^{\mathrm{od}}-1}}}{\underset{i^{\mathrm{od}}=1}{\overset{c^{\mathrm{od}}}{.Math.}}{\left(d_{i^{\mathrm{od}}{j}^{\mathrm{od}}}\right)}^{-\frac{2}{m^{\mathrm{od}}-1}}}& \left(6\right)\end{array}$$\begin{array}{cc}{v}_{i^{\mathrm{od}}}^{\mathrm{od}}=\frac{\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{w_{j^{\mathrm{od}}}^{\mathrm{od}}\left(u_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^{m^{\mathrm{od}}}{x}_{j^{\mathrm{od}}}^{\mathrm{od}}}{\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}{w_{j^{\mathrm{od}}}^{\mathrm{od}}\left(u_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^{m^{\mathrm{od}}}}& \left(7\right)\end{array}$$\begin{array}{cc}{w}_{j^{\mathrm{od}}}^{\mathrm{od}}={\left[\frac{{\left[\underset{j^{\mathrm{od}}=1}{\overset{n^{\mathrm{od}}}{.Math.}}\left(\underset{i^{\mathrm{od}}=1}{\overset{c^{\mathrm{od}}}{.Math.}}{\left(u_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^{m^{\mathrm{od}}}{\left(d_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^2\right)\right]}^{1/n^{\mathrm{od}}}}{\underset{i^{\mathrm{od}}=1}{\overset{c^{\mathrm{od}}}{.Math.}}\left(u_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right){m^{\mathrm{od}}\left(d_{i^{\mathrm{od}}{j}^{\mathrm{od}}}^{\mathrm{od}}\right)}^2}\right]}^{1/p^{\mathrm{od}}}& \left(8\right)\end{array}$

(56) The steps of outlier detection based on AFCM are as follows:

(57) Step 1: performing AFCM clustering on a time series to be detected to obtain c.sup.od cluster centers which are defined as v.sub.1.sup.od, . . . , v.sub.c.sub.od.sup.od;

(58) Step 2: calculating the minimum distance between each data point in the time series and each cluster center;

(59) Step 3: determining the threshold of an abnormal distance through a box figure. The distance value which is more than the threshold will be determined as an outlier. First, arranging the data from small to large, and respectively calculating an upper quartile Q.sub.3.sup.od and a lower quartile Q.sub.1.sup.od; then calculating an interquartile range IQR.sup.od=Q.sub.3.sup.od−Q.sub.1.sup.od; and finally, calculating an outlier cut-off point Q.sup.od=Q.sub.3.sup.od+1.5*IQR.sup.od.

(60) The accuracy of the method is judged by a misjudgment ratio η.sub.w.sub.od.sup.od and an omission factor η.sub.ƒ.sub.od.sup.od, of which the formulae are respectively shown in formula (9) and formula (10):

(61) $\begin{array}{cc}{\eta }_{w^{\mathrm{od}}}^{\mathrm{od}}=\frac{n_{w^{\mathrm{od}}}^{\mathrm{od}}}{n_{t^{\mathrm{od}}}^{\mathrm{od}}}& \left(9\right)\end{array}$$\begin{array}{cc}{\eta }_{f^{\mathrm{od}}}^{\mathrm{od}}=\frac{n_{f^{\mathrm{od}}}^{\mathrm{od}}}{n_{w_r^{\mathrm{od}}}^{\mathrm{od}}}& \left(10\right)\end{array}$
wherein n.sub.w.sub.od.sup.od and n.sub.t.sub.od.sup.od are respectively the number of misjudgments and the number of detection results, and n.sub.ƒ.sub.od.sup.od and n.sub.w.sub.t.sub.od.sup.od are respectively the number of omissions and the number of actual outliers.

(62) With the input flow of first-stage COG as an example, 30,000 samples are selected, the outlier detection results of the local outlier factor (LOF) and the K-Means method are shown in FIG. 9, and the accuracy statistics are shown in Table 1. It can be seen from the table that the method of the present invention effectively overcomes the problem that other methods cannot realize effective detection in the case of different fluctuation ranges at different moments and has high accuracy and efficiency for data outlier detection in such scenarios.

(63) TABLE-US-00001 TABLE 2 Outlier Detection Accuracy of Aperiodic Feature Data Percentage of Omission Misjudgment Detection Outliers (%) Method Factor (%) Ratio (%) Time(s) 1 LOF 1.5903 0.8274 61 K-means 1.5988 0.7631 153 The present 0.12258 0.12258 51 invention 2 LOF 1.6196 0.8515 61 K-means 1.6198 0.7531 157 The present 0.13482 0.14543 61 invention 5 LOF 1.6275 0.8579 61 K-means 1.6264 0.7817 161 The present 0.13888 0.16993 69 invention

(64) The energy generation/consumption data of some production processes (such as blast furnace iron making process) are similar in form but not strictly identical in duration within each period. Such working conditions can be considered to have class period characteristics, and typical data are shown in FIG. 10. For such data, an outlier detection method based on dynamic time warping (DTW) is proposed. The method comprises the following specific detection steps:

(65) Step 1: dividing historical data samples into n.sup.od sample series [x.sub.1.sup.od, . . . , x.sub.n.sub.od.sup.od] with similar length to the data to be detected according to the periodic characteristics;

(66) Step 2: calculating the sum SumD.sub.i.sub.od.sup.od of DTW distances between the series X.sub.i.sub.od.sup.od and other n.sup.od−1 series. When SumD.sub.k.sub.od.sup.od=min{SumD.sub.k.sub.od.sup.od, . . . , SumD.sub.n.sub.od.sup.od} for the (k.sup.od).sup.th series, taking the corresponding X.sub.k.sub.od.sup.od as a center sample, which is recorded as Y.sup.od. Letting the time warping function be C.sup.od={c.sup.od(1), c.sup.od(2), . . . c.sup.od(N)}, where N represents the path length, c.sup.od(n.sub.x.sub.i.sup.od)=(i.sup.od, j.sup.od) represents the (n.sub.x.sub.i.sup.od).sup.th matching point pair composed of the (i.sup.od).sup.th feature vector x.sub.i.sup.od of the series X.sup.od and the (j.sup.od).sup.th feature vector Y.sub.j.sub.od of the series Y.sup.od, and the distance d.sup.od(i.sup.od, j.sup.od) between the two is the local matching distance. The DTW algorithm minimizes the sum of distances through local optimization, as shown in formula (11):

(67) $\begin{array}{cc}{D}^{\mathrm{od}}\left(m_1^{\mathrm{od}},m_2^{\mathrm{od}}\right)=\min \left\{\underset{n^{\mathrm{od}}=1}{\overset{N^{\mathrm{od}}}{.Math.}}{d}^{\mathrm{od}}\left(x_{i^{\mathrm{od}}}^{\mathrm{od}},y_{j^{\mathrm{od}}}^{\mathrm{od}}\right)\right\}& \left(11\right)\end{array}$

(68) Step 3: with the series X.sup.od of a sample to be detected as the reference, scaling the center sample Y.sup.od based on DTW to obtain Y.sub.*.sup.od, wherein the scaling process meets formula (12):

(69) $\begin{array}{cc}{Y}_{*}^{\mathrm{od}}\left[j^{\mathrm{od}}\right]=\frac{1}{.Math.\left\{k^{\mathrm{od}}.Math.x_{i^{\mathrm{od}}}^{\mathrm{od}}.fwdarw.y_{j^{\mathrm{od}}}^{\mathrm{od}}\right\}.Math.}\underset{x_{i^{\mathrm{od}}}^{\mathrm{od}}.fwdarw.y_{j^{\mathrm{od}}}^{\mathrm{od}}}{.Math.}{y}_{j^{\mathrm{od}}}^{\mathrm{od}}& \left(12\right)\end{array}$
wherein k.sup.od represents the number of points corresponding to the series Y.sup.od and the series X.sup.od in the DTW process, and the symbol x.sub.i.sub.od.sup.od.fwdarw.y.sub.i.sub.od.sup.od represents that the point x.sub.i.sub.od.sup.od in the series X.sup.od corresponds to the point y.sub.i.sub.od.sub.od in the series Y.sup.od.

(70) Step 4: performing AFCM clustering on the scaled series Y.sub.*.sup.od to obtain a cluster center v.sub.*.sup.od, and calculating the minimum distance between each point in the sample to be detected and all cluster centers;

(71) Step 5: obtaining a threshold for determining an abnormal distance by calculation through a box figure. The distance value which is more than the threshold will be determined as an outlier.

(72) With the BFG input flow of No. 1 blast furnace as an example, 3,000 samples are selected, and the effect of outlier detection is shown in FIG. 11. It can be seen from the figure that the method of the present invention can effectively detect abnormal period segments in the time series by extracting period information from class period data, which can effectively solve the problem of data outlier detection in such scenarios.

(73) Missing Value Fill

(74) Due to bad measurement environment and serious noise interference in the industrial production field, acquisition points are easily missing at some moment due to temporary failure of sensors, which affects the subsequent analysis and computing of the model. Generally, the missing points may be discrete or continuous. In the process of filling missing points, in order to maximize the filling effect and minimize the error, it is necessary to adopt different methods to fill data respectively. In view of the above situation, the present invention proposes a filling method for multiple scenarios (discrete single-point missing and continuous missing). In the method, the basic features are analyzed based on the historical sample data and divided into different working condition scenarios, and different filling methods are respectively proposed for discrete single-point missing and continuous multi-point missing.

(75) FIG. 12 shows typical discrete missing. For discrete single-point missing, filling can be divided into two conditions according to the data standard deviation. Under the working condition where the standard deviation is less than ε, data are filled by means of moving average, as shown in formula (13):

(76) 0$\begin{array}{cc}{F}_{t^{\mathrm{tb}}+1}^{\mathrm{tb}}=\frac{1}{n^{\mathrm{tb}}}\left(y_{t^{\mathrm{tb}}}^{\mathrm{tb}}+y_{t^{\mathrm{tb}}-1}^{\mathrm{tb}}+.Math.+y_{t^{\mathrm{tb}}+1-n^{\mathrm{tb}}}^{\mathrm{tb}}\right)& \left(13\right)\end{array}$
wherein y.sub.i.sub.tb.sup.tb is the value in the original sample series at the moment t.sup.tb, F.sub.t.sub.tb.sub.+1.sup.tb is the fill value at the missing moment t.sup.tb+1, and n.sup.tb is the number of sample points.

(77) Under the working condition where the standard deviation is more than ε, data are filled by means of cubic spline interpolation. The index value of non-missing points in the sample series is {x.sub.1.sup.tb, . . . , x.sub.n.sub.tb.sup.tb}, and the corresponding data are {y.sub.1.sup.tb, . . . , y.sub.n.sub.tb.sup.tb}. Then on {x.sub.j.sup.tb, x.sub.j+1.sup.tb}, an expression S.sub.j.sup.tb(x.sup.tb) meets constraint conditions (14):

(78) $\begin{array}{cc}\{\begin{array}{c}{S}_j^{\mathrm{tb}}\left(x^{\mathrm{tb}}\right)={a_j^{\mathrm{tb}}\left(x^{\mathrm{tb}}\right)}^3+{b_j^{\mathrm{tb}}\left(x^{\mathrm{tb}}\right)}^2+c_j^{\mathrm{tb}}{x}^{\mathrm{tb}}+d_j^{\mathrm{tb}}\\ S_j^{\mathrm{tb}}\left(x_j^{\mathrm{tb}}\right)=y_j^{\mathrm{tb}}\\ S_j^{\mathrm{tb}}\left(x_{j+1}^{\mathrm{tb}}\right)=y_{j+1}^{\mathrm{tb}}\\ S_j^{{\mathrm{tb}}^{\prime }}\left(x_{j+1}^{\mathrm{tb}}\right)=S_{j+1}^{{\mathrm{tb}}^{\prime }}\left(x_{j+1}^{\mathrm{tb}}\right)\\ S_j^{{\mathrm{tb}}^{″}}\left(x_{j+1}^{\mathrm{tb}}\right)=S_{j+1}^{{\mathrm{tb}}^{″}}\left(x_{j+1}^{\mathrm{tb}}\right)\end{array}& \left(14\right)\end{array}$
wherein j.sup.tb=1, 2, . . . , n.sup.tb−1. S.sub.j.sup.tb(x.sup.tb) is a third expression on [x.sub.jtb, x.sub.j+1.sup.tb], a.sub.j.sup.tb, b.sub.j.sup.tb, c.sub.j.sup.tb and d.sub.j.sup.tb are respectively the corresponding power coefficients of the expression, and S.sub.jtb′(x.sub.j+1.sup.tb) and S.sub.j.sup.tb′(x.sub.j+1.sup.tb) are respectively the first derivative and the second derivative of the expression S.sub.j.sup.tb(x.sup.tb). The constraint conditions specify interpolation conditions as well as continuity at non-endpoints and continuity of derivatives. Then, the value at the missing moment can be calculated based on the fitted curve, so as to realize filling under such working conditions.

(79) For continuous missing in data, an autoregressive integrated moving average model (ARIMA Model) is adopted to fill segments. FIG. 13 shows the trend of continuous missing data. First, the sample series is differentiated and processed as a stationary series, so as to determine the order of parameters of the model. The mathematical model of the ARIMA model can be described as formula (15):

(80) $\begin{array}{cc}{\Delta }^{d^{\mathrm{tb}}}{y}_{t^{\mathrm{tb}}}^{\mathrm{tb}}={\theta }_0^{\mathrm{tb}}+\underset{i^{\mathrm{tb}}=1}{\overset{p^{\mathrm{tb}}}{.Math.}}{\varphi }_{i^{\mathrm{tb}}-1}^{\mathrm{tb}}{\Delta }^{d^{\mathrm{tb}}}{y}_{t^{\mathrm{tb}}-1}^{\mathrm{tb}}+\underset{j^{\mathrm{tb}}=1}{\overset{q^{\mathrm{tb}}}{.Math.}}{\theta }_{j^{\mathrm{tb}}}^{\mathrm{tb}}{\epsilon }_{i^{\mathrm{tb}}-1}^{\mathrm{tb}}& \left(15\right)\end{array}$
wherein Δ.sup.d.sup.tby.sub.t.sub.tb.sup.tb represents a stationary series of y.sub.t.sub.tb.sup.tb differentiated for d.sup.tb times, and represents the random error series of zero-mean white noise. ϕ.sub.t.sub.tb.sup.tb(i.sup.tb=1, 2, . . . , p.sup.tb) and θ.sub.j.sub.tb.sup.tb(j.sup.tb=1, 2, . . . , q.sup.tb) are parameters to be estimated of the model. The ARIMA model comprises three parameters: p.sup.tb, d.sup.tb and q.sup.tb. First, two complete data segments before and after the missing segment are respectively taken as training samples according to the positions of continuous missing data segments, and tested for stationarity by means of Augmented Dickey-Fuller (ADF) test. The series, if non-stationary, shall be transformed into a stationary series through (d.sup.tb).sup.th differential operation. Second, the correlation is identified by an autocorrelation function (ACF) and a partial autocorrelation function (PACF). If ACF attenuates exponentially to zero, an autoregressive process exists. If PACF attenuates exponentially to zero, a moving average process exists. The Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are taken as evaluation indexes of the model, respectively as shown in formula (16) and formula (17):
AIC=2k.sup.tb−2 ln(L.sup.tb)  (16)
BIC=k.sup.tb ln(n.sup.tb)−2 ln(L.sup.tb)  (17)
wherein k.sup.th is the number of parameters of the model, n.sup.th is the number of samples, L.sup.tb is a likelihood Function, and ln(⋅) represents a natural logarithm function. According to formula (16) and formula (17), a (p.sup.tb, q.sup.tb) combination minimizing AIC or BIC is searched within the range of [0, 5], i.e., the optimum parameter of the model where p.sup.tb and q.sup.tb respectively represent the two parameters of the autoregressive moving average model adopted in the invention, p.sup.tb is the number of autoregressive terms, and q.sup.tb is the number of moving average terms. The complete data segments before and after the missing segment are respectively taken as training samples, the data of the missing segment are forecast based on the model, and the average value of the two forecast results is taken as the fill value. In this way, the sample information can be fully utilized so as to improve the filling accuracy.

(81) To sum up, the filling method comprises the following specific steps:

(82) Step 1: judging the information features of time series data to be filled, comprising the miss rate of data, the missing situations of missing points, the periodicity of time series and the standard deviation of time series. If all the missing points in the time series are discrete missing, executing step 2; if all the missing points in the time series are continuous missing, executing steps 4-5; and if discrete missing points and continuous missing points both exist in the time series, executing step 2 first, and then executing steps 4-5 repeatedly for each continuous missing data segment.

(83) Step 2: if the standard deviation of the series is less than ε, filling discrete points by means of moving average. Otherwise, executing step 3.

(84) Step 3: carrying out filling by means of cubic spline interpolation.

(85) Step 4: filling data based on the ARIMA model. Selecting complete series before and after the missing data segment as samples to forecast the data of the missing segment respectively, so as to calculate the average value of the two forecast results as the fill value.

(86) Step 5: establishing an ARIMA model to forecast the missing data segment, and taking the average value of the two forecast results as the fill value.

(87) Adopting the mean absolute error (MAE), the mean absolute percentage error (MAPE) and the root mean square error (RMSE) to measure the filling effect, wherein the three indexes are respectively shown in formula (18) to formula (20):

(88) $\begin{array}{cc}\mathrm{MAE}=\frac{1}{T^{\mathrm{tb}}}\underset{t^{\mathrm{tb}}=1}{\overset{T^{\mathrm{tb}}}{.Math.}}.Math.y^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right)-y_d^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right).Math.& \left(18\right)\end{array}$$\begin{array}{cc}\mathrm{MAPE}=\frac{\underset{t^{\mathrm{tb}}=1}{\overset{T^{\mathrm{tb}}}{.Math.}}\frac{.Math.y^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right)-y_d^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right).Math.}{y_d^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right)}}{T^{\mathrm{tb}}}\times 100\%& \left(19\right)\end{array}$$\begin{array}{cc}\mathrm{RMSE}=\sqrt{\frac{1}{T^{\mathrm{tb}}}\underset{t^{\mathrm{tb}}=1}{\overset{T^{\mathrm{tb}}}{.Math.}}{\left(y^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right)-y_d^{\mathrm{tb}}\left(t^{\mathrm{tb}}\right)\right)}^2}& \left(20\right)\end{array}$
wherein T.sup.tb is the calculation length, y.sup.tb(t.sup.tb) is a forecast value, and y.sub.d.sup.tb(t.sup.tb) is a measured value.

(89) With the data of BFG use flow for cold rolling as an example, selecting the threshold ε=20 of the standard deviation according to the process experience. Adopting adjacent, spline interpolation (SI) and moving average (MA) respectively as comparison methods. FIG. 14 is a comparison diagram of the filling effect of the method of the present invention and that of other methods, and the statistics of the filling accuracy are shown in Table 2. It can be known from the analysis of statistic results that under different data missing scenarios, compared with other methods, the method can obtain the minimum error and improve the filling accuracy. Through complementation between advantages and disadvantages of different filling methods, the method is more suitable for complex data missing environments, and presents a good filling effect for different scenarios, which reflects the effectiveness and robustness of the method.

(90) TABLE-US-00002 TABLE 2 Comparison of Accuracy between the Method of the Present Invention and Other Filling Methods Type of Percentage Standard Filling Missing of Missing Deviation Method MAE MAPE(%) RAISE Adjacent 0.0746 0.0740 0.1717 SI 0.0532 0.0526 0.0781 <20 MA 0.0688 0.0671 0.1291 The method of the 0.0422 0.0414 0.0522 present invention  5% Adjacent 0.9942 2.0229 29.2096 SI 1.1334 2.3879 40.4975 >20 MA 1.4474 3.0169 63.5708 The method 0.9823 2.0114 29.1760 present invention Continuous Adjacent 0.4858 0.4887 1.9784 SI 0.3279 0.3283 0.8693 <20 MA 0.4506 0.4529 1.6783 The method of the 0.2588 0.2586 0.5860 present invention 15% Adjacent 3.1986 4.6025 93.4125 SI 3.2516 5.6826 112.4580 >20 MA 3.8434 6.4286 134.4054 The method of the 3.0635 4.4689 84.8957 present invention Adjacent 0.0441 0.2315 0.0366 SI 0.0717 0.3879 0.1166 <20 MA 0.0465 0.2526 0.0505 The method of the 0.0351 0.1898 0.0288 present invention Discrete  5% Adjacent 0.1296 0.1415 0.4110 SI 0.4084 0.5228 0.7718 >20 MA 0.5953 0.6449 7.1282 The method of the 0.1197 0.1295 0.4455 present invention Adjacent 0.0818 0.4482 0.0630 SI 0.0920 0.4999 0.0775 <20 MA 0.1729 0.9418 0.2204 The method of the 0.0704 0.4295 0.0597 present invention 15% Adjacent 0.6404 0.7193 3.6310 SI 1.1977 1.3247 14.4117 >20 MA 1.7011 1.8690 20.1054 The method of the 0.4789 0.5298 2.0686 present invention

(91) 2) Energy Trend Forecast

(92) Energy Generation-Consumption Forecasting for Class Period Feature Working Conditions

(93) In the production process of the steel industry, because some processes (such as iron making and coking) have relatively mature production flows and fixed processes, energy generation/consumption has class period features, such as BFG generation flow and BFG consumption flow of a hot blast stove. In view of energy forecasting problems under such working conditions, the present invention proposes a combination forecasting model based on the echo state network (ESN) and the least squares support vector machine (LSSVM). The model can be respectively used for dynamic modeling and forecasting for data quality at different moments, and with the training accuracy as the evaluation index, the model with higher accuracy is taken as the final forecast result.

(94) ESN is a novel recursive neural network structure, and has dynamic characteristics of nonlinear systems because of a large dynamical reservoir (DR) inside. The reservoir comprises a large number of randomly generated and sparsely connected neurons, contains the operating state of the system and has the short-term memory function, thus presenting better performance in time series forecasting. Network updating is shown in formula (21) and formula (22):
x.sup.ss(k.sup.ss+1)=ƒ.sup.ss(W.sup.inu.sup.ss(k.sup.ss+1)+Wx.sup.ss(k.sup.ss)+W.sup.backy.sup.ss(k.sup.ss))  (21)
y.sup.ss(k.sup.ss+1)=ƒ.sup.out(W.sup.out(u.sup.ss(k.sup.ss+1),x.sup.ss(k.sup.ss+1),y.sup.ss(k.sup.ss)))  (22)
wherein ƒ.sup.ss(⋅) is a DR internal activation function, which is generally a Sigmoid function tanh( ), x.sup.ss(k.sup.ss) and u.sup.ss(k.sup.ss) are respectively the state variable of DR and the input vector of the model at the moment k.sup.ss, and y.sup.ss(k.sup.ss) is the output of the model. W.sup.in is an input weight matrix. W is an internal neuron connection matrix of DR. To enable DR to have the dynamic memory ability, a sparse connectivity of 1%-5% is generally maintained, and the spectral radius is less than 1. W.sup.back is a feedback matrix of output neurons and DR. ƒ.sup.out(⋅) is an input and output unit activation function, which is generally a linear function. W′ is an output weight matrix.

(95) Different from ESN, an LSSVM regressive model considers a set{x.sub.t.sub.ss.sup.ss, y.sub.t.sub.ss.sup.ss} containing l.sup.ss samples, wherein t.sup.ss=1, . . . , l.sup.ss, the input sample x.sub.t.sub.ss.sup.ss∈R.sup.m, and the output sample y.sub.t.sub.ss.sup.ss∈R. The regressive model can be expressed as formula (23):
y.sup.ss=w.sup.ssφ.sup.ss(x.sup.ss)+b.sup.ss  (23)
wherein φ.sup.ss(⋅) is a non-linear mapping function, and w.sup.ss and b.sup.ss are respectively a weight and a bias. The solution of the LSSVM model can be summarized as the optimization problem with constraints, as shown in formula (24):

(96) $\begin{array}{cc}\min J\left\{w^{\mathrm{SS}},b^{\mathrm{SS}},e^{\mathrm{SS}}\right\}=\frac{1}{2}{\left(w^{\mathrm{SS}}\right)}^T\left(w^{\mathrm{SS}}\right)+\frac{{\gamma }^{\mathrm{SS}}}{2}\underset{i^{\mathrm{SS}}=1}{\overset{N^{\mathrm{SS}}}{.Math.}}{\left(e_{i^{\mathrm{SS}}}^{\mathrm{SS}}\right)}^{2}s.t.{y_{i^{\mathrm{SS}}}^{\mathrm{SS}}\left(w^{\mathrm{SS}}\right)}^T\phi \left(x_{i^{\mathrm{SS}}}^{\mathrm{SS}}\right)+b^{\mathrm{SS}}+e_{i^{\mathrm{SS}}}^{\mathrm{SS}}& \left(24\right)\end{array}$
wherein γ.sup.ss is a penalty coefficient, N.sup.ss is the total number of samples, and e.sub.t.sub.ss.sup.ss is the fitting error of the (i.sup.ss).sup.th sample. The Lagrange multiplier is introduced to transform constrained optimization into unconstrained optimization so as to solve the optimization problem by using Karush-Kuhn-Tucker Conditions (KKT Conditions), thus obtaining an equation set, as shown in formula (25):

(97) $\begin{array}{cc}\left[\begin{array}{cc}0& 1^T\\ 1& K^{\mathrm{SS}}+\frac{2}{{\gamma }^{\mathrm{SS}}}I\end{array}\right]\left[\begin{array}{c}{b}^{\mathrm{SS}}\\ {\alpha }^{\mathrm{SS}}\end{array}\right]=\left[\begin{array}{c}0\\ y^{\mathrm{SS}}\end{array}\right]& \left(25\right)\end{array}$
wherein K.sup.ss(x.sub.p.sub.ss.sup.ss, x.sub.q.sub.ss.sup.ss)=φ.sup.ss(x.sub.p.sub.ss.sup.ss).sup.Tφ.sup.ss(x.sub.q.sub.ss.sup.ss) is a kernel function, which is generally a Gaussian kernel function. I=diag (1, . . . , 1) is a N.sup.ss-dimensional unit matrix, and 1=[1, 1, . . . , 1].sup.T. After α.sup.ss and b.sup.ss are solved, the expression of the forecasting model can be obtained as shown in formula (26):

(98) $\begin{array}{cc}{y}^{\mathrm{SS}}=\underset{i^{\mathrm{SS}}=1}{\overset{N^{\mathrm{SS}}}{.Math.}}{\alpha }_{i^{\mathrm{SS}}}^{\mathrm{SS}}{K}^{\mathrm{SS}}\left(x^{\mathrm{SS}},x_{i^{\mathrm{SS}}}^{\mathrm{SS}}\right)+b^{\mathrm{SS}}& \left(26\right)\end{array}$

(99) The Particle Swarm Optimization (PSO) algorithm is an evolutionary computing method based on bionics, which is highly favored due to advantages of high computing speed and the ability to jump from local optimal solutions and commonly used to optimize algorithm parameters. Each particle has two features of velocity and position, of which the iterative formulae are shown in formula (27) and formula (28) respectively:
v.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss+1)=ω.sup.ssv.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss)+c.sub.1.sup.ss.Math.rand(1).Math.(P.sub.t.sub.ss.sub.,d.sub.ss.sup.ss−x.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss))+c.sub.2.sup.ss.Math.rand(2).Math.(P.sub.g.sub.ss.sub.,d.sub.ss.sup.ss−x.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss))  (27)
x.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss+1)=x.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss)+v.sub.t.sub.ss.sub.,d.sub.ss.sup.ss(t.sup.ss+1)  (28)
wherein x.sub.t.sub.ss.sub.,d.sub.ss.sup.ss, v.sub.t.sub.ss.sub.,d.sub.ss.sup.ss, P.sub.t.sub.ss.sub.,d.sub.ss.sup.ss and P.sub.g.sub.ss.sub.,d.sub.ss.sup.ss respectively represent the current position, the current velocity and the historically optimal position of the (i.sup.ss).sup.th particle in the d.sup.ss-dimensional and the historically optimal position of all particles in the whole particle swarm in space. An inertia weight ω.sup.ss is used to control the influence of the velocity generated in the previous iteration on the velocity of the current iteration; c.sub.1.sup.ss and c.sub.2.sup.ss are learning factors, generally c.sub.1.sup.ss=c.sub.2.sup.ss==2, and rand(1) and rand(2) are random numbers distributed within (0,1).

(100) The energy generation-consumption forecasting method for class period feature working conditions comprises the following specific steps:

(101) Step 1: dividing historical sample data into a training set and a validation set, wherein the training set is used for training the model, and the validation set is used for validating the model accuracy;

(102) Step 2: for the ESN method, optimizing the embedded dimensions of samples and the number of nodes in the reservoir by the PSO algorithm;

(103) Step 3: for the LSSVM method, optimizing the embedded dimensions of samples and the intrinsic parameters (penalty coefficient and kernel parameter) of the algorithm by the PSO algorithm;

(104) Step 4: evaluating the two methods based on the training accuracy, and taking the computed result of the method with higher accuracy as the forecast result of the model.

(105) With the forecasting of the BFG generation flow of No. 1 blast furnace of a steel enterprise as an example, the accuracy of the forecasting method of the present invention is compared with that of the existing forecasting methods, and the forecast results are evaluated by MAPE and RMSE. FIG. 15 shows forecast results corresponding to the BP neural network, the ARIMA model and the model of the present invention, and the forecast time is 60 min (1 min/point). Table 3 shows the accuracy statistic result of each method. It can be seen that compared with existing methods, the method of the present invention has higher accuracy.

(106) TABLE-US-00003 TABLE 3 Comparison of Accuracy between Existing Methods and the Model of the Present Invention Method MAPE (%) RMSE (km.sup.3/h) BP 6 36 ARIMA 5.6 35 The method of 4.1 22 the present invention

(107) The above model considers the case where the historical data are complete. In the actual production process, it is difficult to avoid missing of data points as mentioned above. Although the filling method can ensure the integrity of samples, the reliability cannot be evaluated because the missing data points cannot be acquired. Aiming at this situation, the present invention proposes a method for direct modeling without filling missing points. The idea of ensemble learning is adopted to train the models based on the relevance vector machine (RVM) and Gaussian processes (GP) to give forecast results. In this way, the time of selecting the optimal filling method can be saved, and the generalization and robustness of the forecasting model under different working conditions can be improved.

(108) Embodiments of the present invention are further described below. Due to low sensitivity of devices or changes in working conditions, time series data acquired in industrial fields generally have the characteristics of missing points and high noise. Therefore, such incomplete data cannot be directly used in time series forecasting. FIG. 16(a) shows a time series containing missing points, wherein white boxes represent missing data points and gray filled boxes represent observed data points. The existing methods often adopt moving average, adjacent and spline interpolation for filling and then use the filled data for forecasting. On one hand, such solution is heavily dependent on the filling method, different filling methods have a great influence on the result, and several attempts are required for selecting a method, which is relatively complicated. On the other hand, such method does not consider noise, so it is difficult to ensure the forecasting accuracy. Aiming at the above situation, the present invention integrates two machine learning algorithms, and builds a time series forecasting method based on an incomplete data set, which can effectively forecast class period time series containing noise and missing points, and the steps are as follows:

(109) Step 1: initial filling and data set building based on phase space reconstruction. Extracting time series data of a period of time from the real-time relational database of an industrial field. Detecting and recording positions of missing points in the time series, and filling all missing points with zero values. The function of phase space reconstruction is to build a training sample set. Supposing a discrete time series {x.sub.1.sup.sl, x.sub.2.sup.sl, . . . , x.sub.N.sub.sl.sup.sl}, wherein x.sub.i.sub.sl.sup.sl represents the value acquired at the (i.sup.sl).sup.th sampling moment, and i.sup.s1=1, 2, . . . , N.sup.sl. The (j.sup.sl).sup.th phase vector x.sub.j.sub.sl.sup.sl is shown in formula (29):
x.sub.j.sub.sl.sup.sl=(s.sub.j.sub.sl.sup.sl,s.sub.j.sub.sl.sub.+1.sup.sl, . . . ,x.sub.j.sub.sl.sub.+p.sub.sl.sub.-1.sup.sl).sup.T  (29)
wherein p.sup.sl is an embedded dimension of the phase vector, and the matrix form of N.sup.sl input samples is X.sup.sl=[x.sub.1.sup.sl, x.sub.2.sup.sl, . . . , x.sub.N.sub.sl.sup.sl].sup.T. FIG. 15(b) represents input/output sample set built by phase space reconstruction, which corresponds to the time series containing missing points shown in FIG. 16(a). It can be seen that missing points in input and output are coupled, and missing points “missk” in FIG. 16(b) have the same value in both input and output.

(110) Step 2: building of RVM model based on incomplete data set. Building an initial RVM model, i.e., building a functional relationship between input and output variables. Supposing the mapping relationship between the (i.sup.sl).sup.th function value y.sub.i.sub.sl.sup.sl∈R and an input vector x.sub.j.sub.sl.sup.sl as shown in formula (30):

(111) $\begin{array}{cc}{x}_{i^{\mathrm{sl}}}^{\mathrm{sl}}=\underset{j^{\mathrm{sl}}=1}{\overset{N^{\mathrm{sl}}}{.Math.}}{w}_{j^{\mathrm{sl}}}^{\mathrm{sl}}K\left(x_{i^{\mathrm{sl}}}^{\mathrm{sl}},x_{j^{\mathrm{sl}}}^{\mathrm{sl}}\right)+w_0^{\mathrm{sl}}+{\mathrm{.Math.}}_{i^{\mathrm{sl}}}^{\mathrm{sl}}={\left({\varphi }_{i^{\mathrm{sl}}}^{\mathrm{sl}}\right)}^T{w}^{\mathrm{sl}}+{\mathrm{.Math.}}_{i^{\mathrm{sl}}}^{\mathrm{sl}}& \left(30\right)\end{array}$
wherein ϕ.sub.i.sub.sl.sup.sl=[1, K(x.sub.i.sub.sl.sup.sl, x.sub.1.sup.sl, K(x.sub.i.sub.sl.sup.sl, x.sub.2.sup.sl), . . . , K(x.sub.i.sub.sl.sup.sl, x.sub.N.sub.sl.sup.sl)].sup.T, and i.sup.sl=1, 2, . . . , N.sup.sl w.sup.sl=[w.sub.0.sup.sl, w.sub.1.sup.sl, . . . , w.sub.N.sub.sl.sup.sl].sup.T represents a weight vector, and is a kernel function which is a Gaussian kernel function in the present invention. ε.sub.i.sub.sl.sup.sl represents a Gaussian distribution noise term with the mean value of 0 and the variance of σ.sup.2(β≡σ.sup.2). The weight vector w.sup.sl is considered as a continuous random variable. Supposing that each element obeys normal distribution with the mean value of 0 and different variances, as shown in formula (31):

(112) $\begin{array}{cc}P\left(w^{\mathrm{sl}}.Math.\alpha \right)=\underset{i^{\mathrm{sl}}=0}{\overset{N^{\mathrm{sl}}}{.Math.}}N\left(w_{i^{\mathrm{sl}}}^{\mathrm{sl}}.Math.0,{\alpha }_i^{-1}\right)& \left(31\right)\end{array}$
wherein α=[α.sub.0, α.sub.1 . . . , α.sub.N].sup.T, and A=diag(α). Parameters α and β are considered as random variables which obey Gamma distribution. In the present invention, a noise variance vector is written as β=[β.sub.[m], β.sub.[m], . . . , β.sub.[m], β.sub.[obs], . . . , β.sub.[obs], β.sub.[obs]].sup.T, and β.sub.[obs]=σ.sup.−2. In addition, β also can be written as β=[β.sub.[m].sup.T, β.sub.[obs].sup.T].sup.T in the form of diagonal matrices B.sub.[m]=diag(β.sub.[m]), B.sub.[obs]=diag(β.sub.[obz]) and B=diag(β). A training output set is composed of missing points observable points, wherein x.sup.sl=[(s.sub.[m].sup.sl).sup.T, (x[obs].sup.sl).sup.T].sup.T, s.sub.[m].sup.sl and s.sub.[obs].sup.sl respectively represent the vectors of missing points and observable points in the output set, the number of output missing points and that of output observable points are respectively N.sup.sl.sub.ymiss and N.sup.sl.sub.yobs, and N.sup.sl=N.sup.sl.sub.ymiss+N.sup.sl.sub.yobs. In the present invention, the element β.sub.[m] in β.sub.[m] is defined as 10.sup.8 to ensure that the noise variance of the corresponding missing output point is not affected. The likelihood function is shown in formula (32):

(113) $\begin{array}{cc}p\left(x_{\left[m\right]}^{\mathrm{sl}},x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}.Math.w,\beta \right)=N\left(x^{\mathrm{sl}}.Math.{\Phi }^{\mathrm{sl}}{w}^{\mathrm{sl}},B^{-1}\right)=1/\left[{\left(2\pi \right)}^{N^{\mathrm{sl}}/2}{.Math.{\left(B^{\mathrm{sl}}\right)}^{-1}.Math.}^{1/2}\right]\times \exp \left\{-{\frac{1}{2}\left[x^{\mathrm{sl}}-\Phi w^{\mathrm{sl}}\right]}^{T}B\left[x^{\mathrm{sl}}-\Phi w^{\mathrm{sl}}\right]\right\}& \left(32\right)\end{array}$
wherein Φ.sup.d=[ϕ.sup.sl.sub.1, ϕ.sup.sl.sub.2, . . . , ϕ.sup.sl.sub.N.sub.sl].sup.T is a kernel function matrix. Aiming at the optimization problem of a missing input point x.sup.sl.sub.[m], and considering the characteristics of the sample set built by phase ace reconstruction, each missing input point has a corresponding equivalent missing point in the missing output set. Therefore, the present invention replaces the posterior probability of each element in x.sup.sl.sub.[m] with the posterior probability of the corresponding missing output point, fills the posterior average of x.sup.sl.sub.[m] in the corresponding position of the missing input point, and further updates the kernel matrix Φ.sup.sl as shown in formula (33):

(114) 0$\begin{array}{cc}{\mu }_{x_{\left[m\right]},i_{x_{\left[m\right]}}^{\mathrm{sl}}}=\int x_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}p\left(x_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}.Math.x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}},\alpha ,{\beta }_{\left[\mathrm{obs}\right]}\right){\mathrm{dx}}_{\left[m\right],{\mathrm{corr}}_{-}i}^{\mathrm{sl}}& \left(33\right)\end{array}$
wherein i.sub.x.sub.[m].sup.sl=1, 2, . . . , N.sup.sl.sub.x.sub.[m], the set x.sup.sl.sub.[m] contains N.sup.sl.sub.x.sub.[m] missing input points, x.sup.sl.sub.[m],i.sup.sl represents the (i.sup.sl).sup.th element in x.sup.sl.sub.[m], and μ.sub.x.sub.[m].sub.,i.sub.sl represents the mean posterior probability of x.sub.[m],i.sub.sl.sup.sl. x.sup.sl.sub.[m],corr_i represents the corresponding variable of the (i.sup.sl).sup.th missing input point in the missing output set x.sup.sl.sub.[m]. By eliminating the integrals of w.sup.sl and x.sup.sl.sub.[m], the present invention gives a log edge likelihood function as shown in formula (34):

(115) $\begin{array}{cc}\tilde{L}=\log \int p\left(x_{\left[m\right]}^{\mathrm{sl}},x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}.Math.w^{\mathrm{sl}},{\beta }_{\left[\mathrm{obs}\right]}\right)p\left(w^{\mathrm{sl}}.Math.\alpha \right){\mathrm{dx}}_{\left[m\right]}^{\mathrm{sl}}{\mathrm{dw}}^{\mathrm{sl}}=-\frac{1}{2}\left(N^{\mathrm{sl}}\log 2\pi +\log .Math.\tilde{C}.Math.+{\left(x^{\mathrm{sl}}\right)}_{\left[\mathrm{obs}\right]}^T{\tilde{C}}^{-1}{x}_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}\right)& \left(34\right)\end{array}$
wherein a covariance matrix {tilde over (C)}=Φx[obs]A.sup.−1Φ.sub.x.sub.[obs].sup.T+β.sub.[obs].sup.−1I, and I is a unit matrix. The extreme values α* and β*.sub.[obs]of α and β.sub.[obs] are found by maximizing partial lead {tilde over ( )} L, and the iterative formulae are obtained as shown in formula (35) and formula (36):

(116) $\begin{array}{cc}{\alpha }_{i^{\mathrm{sl}}}=\frac{{\tilde{\gamma }}_{i^{\mathrm{sl}}}}{{\tilde{\mu }}_{i^{\mathrm{sl}}}^2}& \left(35\right)\end{array}$$\begin{array}{cc}{\beta }_{\left[\mathrm{obs}\right]}^{-1}=\frac{{.Math.x_{\left[\mathrm{obs}\right]}^{\mathrm{sl}}-{\Phi }_{x_{\left[\mathrm{obs}\right]}}\tilde{\mu }.Math.}^2}{N_{\mathrm{yobs}}^{\mathrm{sl}}-{\beta }_{\left[\mathrm{obs}\right]}\mathrm{Tr}\left(\widetilde{.Math.}{\Phi }_{x_{\left[\mathrm{obs}\right]}}^T{\Phi }_{x_{\left[\mathrm{obs}\right]}}\right)}& \left(36\right)\end{array}$
wherein i.sup.sl=1, 2, . . . , N.sup.sl, {tilde over (γ)}.sub.i=1−α.sub.i{tilde over (Σ)}.sub.ii, {tilde over (Σ)}=(A+β.sub.[obs]Φ.sub.x.sub.[obs].sup.TΦ.sub.x.sub.[obs]).sup.−1 and {tilde over (μ)}=β.sub.[obs]{tilde over (Σ)}Φ.sub.x.sub.[obs].sup.Tx.sup.sl.sub.[obs]. The method comprises the following steps:

(117) {circle around (1)} Initializing: α=[1/(N.sup.sl).sup.2, 1/(N.sup.sl).sup.2, . . . , 1/(N.sup.sl).sup.2].sup.T, and β.sub.[obs] is initialized as the variance of 10 times x.sup.sl.sub.[obs]. The related vector set is initialized as a set of all input samples;

(118) {circle around (2)} Updating mean values and covariance matrices of joint posterior probabilities of w.sup.sl and x.sup.sl[m];

(119) {circle around (3)} Updating α and β.sub.[obs];

(120) {circle around (4)} If α.sub.j.sub.sl>10.sup.8, removing the (j.sup.sl).sup.th sample from the related vector set;

(121) {circle around (5)} Calculating the mean posterior probability of each element in x.sup.sl.sub.[m] filling in the corresponding position of the missing input point, and further updating the kernel matrix Φ;

(122) {circle around (6)} Calculating the value of the log edge likelihood function;

(123) {circle around (7)} If the change of an element in a relative to the last iteration or the change of the value of the edge likelihood function is less than the given threshold, stopping the iterative process, thus obtaining a final forecasting model; otherwise, executing step {circle around (2)}.

(124) Step 3: obtaining the forecast value of a future period of time based on RVM single-step iterative forecasting. When a new input x.sup.sl.sub.* (the input vector contains no missing point) is given, the value of β.sub.[obs] can be replaced with β.sub.[obs]. The forecast distribution of the corresponding output x.sup.sl.sub.* is approximately normal distribution, as shown in formula (37):
p(t.sup.sl.sub.*|t.sup.sl.sub.[obs])≈∫p(t.sup.sl.sub.*|x.sup.sl.sub.*,β*.sub.[obs])p(w|t.sup.sl.sub.[obs],α,β*.sub.[obs])dw=N(t.sup.sl.sub.*|μ.sub.*,σ.sub.*.sup.2)  (37)
wherein μ.sub.*=ϕ.sub.*.sup.Tμ.sub.[w|obs], σ.sub.*.sup.2=(β*.sub.[obs]).sup.−1+ϕ.sub.*.sup.T{tilde over (Σ)}ϕ.sub.* and ϕ.sub.*=[1, K(x.sup.sl.sub.*, x.sup.sl.sub.1), K(x.sup.sl.sub.*, x.sup.sl.sub.2), . . . , K(x.sup.sl.sub.*, x.sup.sl.sub.N.sub.sl)].sup.T. μ.sub.[w|obs] is the mean posterior probability of the weight vector w.sup.sl. Then, the forecast mean value of a future period of time is a point, and the next round of forecasting will take the forecast result of the previous round as input to carry out cyclic rolling forecast.

(125) Step 4: building of LSSVM data set based on incomplete data set. After the parameters α and b are solved, using x.sup.sl.sub.*=αK(x.sup.sl, x.sup.sl.sub.*)+b to calculate an estimated value of the forecast quantity. Therefore, inputting data containing missing points filled by moving average into the LSSVM model for training, then forecasting the data containing missing points, replacing the original values of missing points with the forecast result, and repeating the above process until the change of the values of missing points is less than the set threshold, thus completing the data updating process.

(126) Step 5: obtaining the forecast value of a future period of time based on GP forecasting. After updating the missing points by LSSVM, inputting the updated data into the GP model, wherein GP is a probability model with kernel. Supposing that the weight w.sup.sl and the noise ε.sup.sl in formula (38) obey zero-mean prior distribution:
p(w.sup.sl)=N(w.sup.sl|0,α.sup.−1I)
p(ε.sup.sl)=N(ε.sup.sl|0,β.sup.−1)  (38)
wherein α and β are hyperparameters of the weight and the noise distribution, and the output likelihood function can be expressed as formula (39):
p(x.sup.sl|Φ,w.sup.sl,β)=N(x.sup.sl|Φw.sup.sl,β.sup.−1I)  (39)

(127) According to Bayes formula, deducing the posterior distribution of the weight as shown in formula (40):

(128) $\begin{array}{cc}p\left(w^{\mathrm{sl}}.Math.x^{\mathrm{sl}},\Phi ,\alpha ,\beta \right)=\frac{p\left(x^{\mathrm{sl}}.Math.\Phi ,w^{\mathrm{sl}}\right)p\left(w^{\mathrm{sl}}\right)}{p\left(x^{\mathrm{sl}}.Math.\Phi \right)}=N\left(w^{\mathrm{sl}}.Math.m_{N^{\mathrm{sl}}},.Math.\right)& \left(40\right)\end{array}$
wherein m.sub.N=βΣΦ.sup.Tx.sup.sl and Σ=(αI+βΦ.sup.TΦ).sup.−1. Then selecting the Gaussian kernel function and maximizing the likelihood function, thus obtaining optimal hyperparameters.

(129) Step 6: algorithm selection. Calculating the MAPE of the validation sets obtained by the RVM algorithm and the GP algorithm respectively, selecting the optimal algorithm, and saving and outputting the hyperparameter values and the information of the kernel function optimized by the corresponding algorithm.

(130) Step 7: forecasting. After real data are acquired, using the optimal algorithm and the parameters thereof for forecasting.

(131) FIG. 17 shows comparison diagrams of forecast results based on nearest neighbor interpolation (NNI), spline interpolation (SI) and the present invention when the data of the training set are missing by 5% (FIG. 17(a)), 10% (FIG. 17(b)) and 30% (FIG. 17(c)), and the accuracy statistics are shown in Table 4. The core algorithm of the comparison method uses RVM to forecast the time series. RMSE and MAPE are taken as the evaluation indexes of the forecast results. It can be known from the statistic result in the table that the method of the present invention has higher forecasting accuracy and can directly model and forecast the data containing missing input without separate filling, which improves the computing accuracy and efficiency.

(132) TABLE-US-00004 TABLE 4 Comparison of Forecast Results of Three Methods Percentage of Forecasting Missing Data Method MAPE (%) RMSE 5% NNI 5.85 36.99 SI 5.28 32.12 The present 4.47 28.46 invention 10% NNI 6.86 37.92 SI 6.59 33.60 The present 5.59 32.91 invention 30% NNI 8.37 46.34 SI 7.37 38.00 The present 6.70 36.41 invention

(133) Energy Generation-Consumption Forecasting for Aperiodic Feature Working Conditions

(134) Another common working condition of production is energy consumption entirely based on the production process, but has no obvious periodic feature (such as cold rolling and hot rolling). In view of energy forecasting problems under such aperiodic feature working conditions, the present invention builds a bidirectional multi-layer long short-term memory forecasting model based on adaptive variational modal decomposition for noise reduction. First, complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) is adopted to decompose energy series, so as to calculate the initial parameters of variational mode decomposition (VMD), which can reduce the noise of the original data, smooth the data and reduce the aperiodicity and uncertainty of the data. Then, the Bayesian optimization method is used to obtain the optimal VMD decomposition parameter M, and the decomposition parameters of the VMD method are determined according to the decomposition results, which can search better parameters within a fixed number of iterations and shorten the time of parameter optimization. The VMD method is used to reduce the noise of the series, then the data are divided into a training set and a validation set, and input/output samples are built. Furthermore, bidirectional multi-layer long short-term memory (LSTM) is built for training, which is used to capture deep features in the series and has better generalization effect than single LSTM. Meanwhile, the measure of early stopping is taken for parameter optimization of the forecasting model, the direction of grid search is from a simple model to a complex model, and search is stopped when the error range is less than a certain degree, so as to shorten the time of parameter optimization. With the structure of decomposition for noise reduction and complex network, the present invention can avoid, with a higher probability, the hysteresis phenomenon of traditional models in forecasting of aperiodic data, and has better generalization ability.

(135) The method comprises the following specific implementation steps:

(136) Step 1: preparing data.

(137) Extracting time series data of BFG consumption flow for cold rolling of a period of time from the real-time relational database of an industrial field.

(138) Step 2: acquiring VMD parameters optimized by the Bayesian optimization method and decomposing and reorganizing the data.

(139) Different from EMD and derived methods thereof, the VMD method is improved based on Hibert transform and Wiener filtering. The variational mode decomposition method is a quasi-orthogonal completely non-recursive decomposition method, which applies a variational problem solving process to the field of signal decomposition, transforms recursive decomposition of the EMD method into non-recursive variational decomposition, and highlights the local features of signals by means of iteratively solving the center frequency and bandwidth of each mode function. VMD has better denoising effect, avoids the phenomenon of spectrum aliasing, and can distinguish series with similar frequencies accurately.

(140) For the original series ƒ.sup.ns(t), the problem of constraint variation is shown in formula (41):

(141) $\begin{array}{cc}\underset{\left\{u_{k^{\mathrm{ns}}}^{\mathrm{ns}},{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\}}{\min }\left\{\underset{k^{\mathrm{ns}}=1}{\overset{N^{\mathrm{ns}}}{.Math.}}{.Math.{\partial }_t\left[\left(\delta \left(t\right)+\frac{j^{\mathrm{ns}}}{\pi t}\right)u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)\right]e^{-j{\omega }_{k}t}.Math.}_2^2\right\}s.t.\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)=f^{\mathrm{ns}}\left(t\right)& \left(41\right)\end{array}$
wherein K.sup.ns is the number of modes, u.sub.k.sub.ns.sup.ns(t) is the (k.sup.ns).sup.th mode, ω.sub.k.sub.ns.sup.ns is the center frequency, and g(t) represents Dirac distribution. To solve the above problem, the quadratic penalty term and the Lagrange multiplier λ.sup.ns are adopted to transform the above problem into an unconstrained problem, as shown in formula (42):

(142) $\begin{array}{cc}L\left(\left\{u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\},\left\{{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}\right\},{\lambda }^{\mathrm{ns}}\right)={\alpha }^{\mathrm{ns}}\underset{k^{\mathrm{ns}}}{.Math.}{.Math.{\partial }_t\left[\left(\delta \left(t\right)+\frac{j^{\mathrm{ns}}}{\pi t}\right)*u_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right)\right]e^{-j{\omega }_{k}t}.Math.}_2^2+{.Math.f^{\mathrm{ns}}\left(t\right)-\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right).Math.}_2^2+.Math.{\lambda }^{\mathrm{ns}}\left(t\right),f^{\mathrm{ns}}\left(t\right)-\underset{k^{\mathrm{ns}}}{.Math.}{u}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left(t\right).Math.& \left(42\right)\end{array}$

(143) In order to solve the mode u.sub.k.sub.ns.sup.ns and the center frequency ω.sub.k.sub.ns.sup.ns, the problem in the above formula is transformed into a minimization problem and solved by an alternating direction method of multipliers (ADMM), and then the optimal solutions can be expressed as formula (43) and formula (44):

(144) $\begin{array}{cc}{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)=\frac{{\hat{f}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)-\underset{i^{\mathrm{ns}}\ne k^{\mathrm{ns}}}{.Math.}{\hat{u}}_{i^{\mathrm{ns}}}\left({\omega }^{\mathrm{ns}}\right)+\left({\hat{\lambda }}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)/2\right)}{1+2{{\alpha }^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}-{\omega }_{K^{\mathrm{ns}}}^{\mathrm{ns}}\right)}^2}& \left(43\right)\end{array}$$\begin{array}{cc}{\omega }_{k^{\mathrm{ns}}}^{\mathrm{ns}}=\frac{{\int }_0^{\infty }{\omega }^{\mathrm{ns}}{.Math.{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right).Math.}^{2}d{\omega }^{\mathrm{ns}}}{{\int }_0^{\infty }{.Math.{\hat{u}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right).Math.}^{2}d{\omega }^{\mathrm{ns}}}& \left(44\right)\end{array}$
wherein {circumflex over (ƒ)}.sup.ns(ωns), u.sub.i.sub.ns.sup.ns(ωns), {umlaut over (λ)}.sup.ns(ω.sup.ns) and ü.sub.k.sub.ns.sup.ns(ω.sup.ns) represent Fourier transform forms of corresponding variables respectively. After u.sub.k.sub.ns.sup.ns and ω.sub.k.sub.ns.sup.ns are updated, the Lagrange multiplier is also updated according to formula (45):

(145) $\begin{array}{cc}{\hat{\lambda }}_{t+1}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)={\hat{\lambda }}_t^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)+\tau \left({\hat{f}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)-\underset{k^{\mathrm{ns}}}{.Math.}{\hat{u}}_{k^{\mathrm{ns}}}^{\mathrm{ns}}\left({\omega }^{\mathrm{ns}}\right)\right)& \left(45\right)\end{array}$

(146) Then the above updating process is repeated until the convergence condition of formula (46) is satisfied:

(147) $\begin{array}{cc}\underset{k^{\mathrm{ns}}}{.Math.}\frac{{.Math.{\hat{u}}_{k^{\mathrm{ns}}}\left(t+1\right)-{\hat{u}}_{k^{\mathrm{ns}}}\left(t\right).Math.}_2^2}{{.Math.{\hat{u}}_{k^{\mathrm{ns}}}\left(t\right).Math.}_2}<{\epsilon }^{\mathrm{ns}}& \left(46\right)\end{array}$
wherein ε.sup.ns is a threshold satisfying the convergence condition. The Bayesian parameter optimization process generally comprises two parts: First, the previous function must be selected to represent the hypothesis about the optimized function. With greater flexibility and tractability, Gauss processes are generally selected as a probability model for Bayesian parameter optimization. Second, Bayesian optimization needs an acquisition parameter for building an acquisition function from model posteriori to determine the next parameter combination to be estimated.

(148) Supposing that a set of hyperparameters is X.sup.ns=x.sup.ns.sub.1, x.sup.ns.sub.2, . . . , x.sup.ns.sub.n.sub.ns, which has a functional relationship with the target to be optimized. For ƒ.sup.ns(X.sup.ns), supposing that GP is a Gaussian process model of input data X.sup.ns and an objective function ƒ.sup.ns, and S.sup.ns is an acquisition function. First, time series data are decomposed by CEEMDAN to obtain the number K.sup.ns.sub.k.sub.ns of VMD reserved series. For the VMD method, the decomposition number K.sup.ns=.sup.ns.sub.d.sub.ns+K.sup.ns.sub.k.sub.ns, wherein the number K.sup.ns.sub.d.sub.ns of discarded series and the balance factor α.sup.ns are optimized by the Bayesian optimization method. The objective function of the Bayesian optimization method of the present invention is expressed in formula (47):
ƒ.sub.obj.sup.ns(K.sub.d.sub.ns.sup.ns,α.sup.ns)=e.sub.origin.sup.ns+e.sub.model.sup.ns  (47)
wherein e.sub.origin.sup.ns, represents the error between the series composed of first K.sub.k.sub.ns.sup.ns components after decomposition by the VMD method and the original series, and e.sub.model.sup.ns is the error of VMD reorganized series in model training. If only e.sub.model.sup.ns is used as the objective function, and K.sub.d.sub.ns.sup.ns will be guided to be optimized to a larger direction, which makes the series smoother and easier to forecast, but causes the error from real series to be greater and greater. After the optimization parameters K.sup.ns and α.sup.ns of VMD are obtained, the original series is decomposed, the first K.sub.k.sub.ns.sup.ns series are added and combined into a new time series, and input/output samples are built by using the time series and then divided into a training set and a validation set.

(149) Step 3: training a bidirectional multi-layer LSTM model. The parameters of the bidirectional multi-layer LSTM model are mainly as follows: the input dimension, the number of LSTM layers, the number of hidden nodes of each LSTM layer, the learning rate and the number of samples of a training batch. The parameter optimization method is grid search with early stopping, wherein the input dimension is optimized first, then the number of LSTM layers and the number of hidden nodes of each LSTM layer are optimized, and finally, the learning rate and the number of samples of a training batch are optimized. The specific steps are as follows:

(150) {circle around (1)} Determining the optimization range and granularity of the parameters, the optimization direction is decreasing model complexity, i.e., the number of LSTM layers is decreased, the number of hidden nodes of each LSTM layer is decreased, and the objective function is the root-mean-square error;

(151) {circle around (2)} Recording errors between searches during grid search, and stopping search when the error is less than the threshold ò;

(152) {circle around (3)} When optimizing the number of LSTM layers and the number of hidden nodes of each LSTM layer, setting the learning rate and the number of samples of a training batch to the maximum values to shorten the search time, and then using the optimized parameters to train the bidirectional multi-layer LSTM model, thus obtaining a forecasting model.

(153) Step 4: carrying out forecasting based on the latest data acquired from the database.

(154) As the VMD method requires a certain amount of data to ensure the stability of decomposition results, the present invention takes the latest 1,000 sample points, performs VMD decomposition and reorganization to reduce noise first, and then takes the first input dimension sample points as the input to forecast points at the next moment.

(155) With the data forecasting of BFG consumption flow for cold rolling as an example, FIG. 18 is a comparison diagram of forecast effects of the LSTM method, the LSTM method with decomposition for noise reduction and the present invention, and the accuracy statistics of forecast results are shown in Table 5. It can be seen from the statistic result that the method of the present invention is superior to other methods in the two evaluation indexes: MAPE and RMSE, indicating the effectiveness of the method.

(156) TABLE-US-00005 TABLE 5 Comparison of Forecast Results Method RMSE MAPE LSTM 0.8990 2.624% VMD_LSTM 0.8931 2.300% The present 0.4684 1.356% invention

(157) For the aperiodic feature time series containing missing points, the present invention first fills the missing data by cubic spline interpolation, as shown in formula (15) and formula (16), and then adopts a forecasting method for complete aperiodic feature data to forecast such data.

(158) With the data of BFG consumption flow for hot rolling as an example, FIG. 19 is a comparison diagram of forecast results of a series with the miss rate of 15% in moving average pre-processing and cubic spline interpolation pre-processing. The statistics of forecasting accuracy are shown in Table 6. It can be seen from the table that compared with other methods, the method of the present invention has higher accuracy and can provide effective guidance for practical industrial application.

(159) TABLE-US-00006 TABLE 6 Accuracy Statistics of Forecast Results Method RMSE MAPE Forecasting method based on 2.8140 6.745% moving average The method of the present 0.9875 2.361% invention

(160) Forecasting of Gas Reserves for Factor Model

(161) Since the capacity of a gas holder is taken as an important sign of balance of a gas system, the forecasting accuracy thereof plays an extremely important role. To solve the short-term forecasting problem, the present invention proposes a forecasting method considering the capacity fluctuation factor. First, the experimental data acquired in an industrial field are divided into a training set and a validation set, and a relevance vector machine model, a Gaussian process model and a least squares support vector machine model are established respectively; and an optimal model is selected through the error evaluation index for forecasting. The present invention adopts the idea of ensemble algorithm to build a factor model forecasting method integrating the RVM algorithm, the GP algorithm and the LSSVM algorithm, combines multiple machine learning algorithms, and compares validation errors of the models to finally obtain a stable model with good performance in all aspects, which greatly improves the forecasting accuracy and the algorithm stability.

(162) The present invention comprises the following specific implementation steps:

(163) Step 1: pre-processing data.

(164) Extracting the input flow of 4 blast furnace gas pipe networks of a period of time and the consumption flow of 24 gas consumption links from the real-time relational database of an industrial field as input factors of the model, and adding the sum of the capacities at the previous moment to the input factors; and supposing that the sum of capacities at the current moment is the output value. Meanwhile, normalizing all the input data and the output data. Selecting the first 70% of samples as training samples and the remaining 30% as validation samples.

(165) Step 2: building a relevance vector machine regressive model;

(166) Step 3: building a Gaussian process regressive model;

(167) Step 4: building an LSSVM model;

(168) Step 5: selecting an optimal model, and forecasting the estimated value of a forecast point of a future period of time.

(169) After the output estimated values of the validation sets of the RVM model, the GP model and the LSSVM model are respectively obtained by using the validation sets, the present invention compares the model performance by using the error evaluation index MAPE, so as to obtain the optimal model. The smaller the MAPE value is, the more accurate the forecast effect of the model is.

(170) FIG. 20 is a comparison diagram of forecast effects of blast furnace gas capacity of the error back propagation neural network (BP) method, the support vector machine (SVM) method and the present invention, and the comparison of the forecasting accuracy is shown in Table 7. It can be known from the statistic result that the method of the present invention has higher accuracy and can provide more reliable support for subsequent optimal scheduling.

(171) TABLE-US-00007 TABLE 7 Comparison of Forecast Results of Three Methods Forecasting Method MAPE RMSE BP 0.0079 2.0225 SVM 0.0439 8.1252 The present 0.0060 1.2988 invention

(172) 3) Scheduling of multiple by-product gas systems The balance of supply and demand of each system is analyzed based on the trend forecast result of the generation-consumption-storage unit of each gas system. When energy safety events occur, it is necessary to give a reasonable scheduling interval according to the safety constraints of the gas capacity, so as to ensure the safe operation of the by-product gas system. First, according to the capacity operation states of the blast furnace gas (BFG) system, the coke oven gas (COG) system and the linz-donawitz converter gas (LDG) system, the balance of each system is judged with system safety indexes as constraints to give scheduling suggestions (whether adjustment is needed and the corresponding adjustable range) of each single-medium system. Second, with the adjustable range of each single-gas medium and the load capacity of each adjustable user (such as power plant and boiler) as constraints, the comprehensive adjustable range of multiple gas media is calculated. Finally, the adjustable range for optimal scheduling is obtained by integrating the adjustable range of each pipe network and the adjustable range of multiple media. The calculation method for the adjustable range of each stage is as follows:

(173) {circle around (1)} Solving the adjustable range of a single medium

(174) The calculation process of the adjustable range of a single medium is shown in FIG. 21. According to the balance state of the system, a two-stage judgment method is proposed. The first stage is a short-term period, and the second stage is a long-term period. Each stage is divided into three parts, the state (ultrahigh high limit, ultrahigh limit, etc.) of each part is judged, and the overall state of each stage is judged according to the state of each part, so as to judge the adjustment of the pipe network. The system adjustable range is obtained by circularly judging the state of the pipe network under each adjustment value. When the capacity state under this adjustment value is judged to require no adjustment, the adjustment value at this moment can be adopted. The adjustment value within the adjustment range is judged circularly, finally obtaining the adjustable range. Supposing that the use flow of blast furnace gas, coke oven gas and linz-donawitz converter gas is x.sup.s, y.sup.s and z.sup.s respectively, the adjustable ranges can be defined as [x.sub.1.sup.s, x.sub.2.sup.s], [y.sub.1.sup.s, y.sub.2.sup.s] and [z.sub.1.sup.s, z.sub.2.sup.s] respectively.

(175) {circle around (2)} Solving the adjustable range of multiple media

(176) According to the principle that the total calorific value of gas used by a multi-media adjusting user is constant, the use range of each medium within the adjustable range under the multi-media adjusting user is calculated as shown in formula (48):
Total_cal=bfg_cal×x+cog_cal×y+ldg_cal×z  (48)
wherein Total_cal is the gross calorific value of the multi-media adjusting user, and bfg_cal, cog_cal and ldg_cal are respectively the calorific value of each gas medium. Formula (49) can be obtained by the elimination method:

(177) $\begin{array}{cc}{x}^s=\frac{{\mathrm{Total}}_{-}{\mathrm{cal}}_{-}\left(\mathrm{cog\_}\mathrm{cal}\times y+\mathrm{ldg\_}\mathrm{cal}\times z\right)}{\mathrm{bfg\_}\mathrm{cal}}& \left(49\right)\end{array}$

(178) The value ranges [x.sub.1.sup.s′, x.sub.2.sup.s′], [y.sub.1.sup.s′, y.sub.2.sup.s′] and [z.sub.1.sup.s′, z.sub.2.sup.s′] of x.sup.s, y.sup.s and z.sup.s in the case of y.sup.s∈[y.sub.1.sup.s, y.sub.2.sup.s], z.sup.s∈[z.sub.1.sup.s, z.sub.2.sup.s] and x.sup.s∈[x.sub.1.sup.s, x.sub.2.sup.s] are solved, that is the adjustable range of each gas medium after multi-media adjustment is considered. Then, the final adjustable ranges [x.sub.1.sup.s″, x.sub.2.sup.s″], [y.sub.1.sup.s″, y.sub.2.sup.s″] and [z.sub.1.sup.s″, z.sub.2.sup.s″] for optimal scheduling are obtained through interval differencing, as shown in formulae (50)-(52) respectively:
[x.sub.1.sup.s″,x.sub.2.sup.s″]=[x.sub.1.sup.s,x.sub.2.sup.s]−[x.sub.1.sup.s′,x.sub.2.sup.s′]  (50)
[y.sub.1.sup.s″,y.sub.2.sup.s″]=[y.sub.1.sup.s,y.sub.2.sup.s]−[y.sub.1.sup.s′,y.sub.2.sup.s′]  (51)
[z.sub.1.sup.s″,z.sub.2.sup.s″]=[z.sub.1.sup.s,z.sub.2.sup.s]−[z.sub.1.sup.s′,z.sub.2.sup.s′]  (52)

(179) FIG. 22(a), FIG. 22(b) and FIG. 22(c) respectively show the adjustable range of each gas system calculated by the method of the present invention in certain scheduling, respectively “−30,000-140,000 m.sup.3/h”, “20,000-40,000 m.sup.3/h” and “−20,000-70,000 m.sup.3/h”. It can be seen from the figure that the method of the present invention can give a reasonable adjustment range of each gas system on the premise of meeting the safety limit of the capacity, so as to provide reliable constraint conditions for subsequent optimal scheduling of comprehensive energy.

(180) 4) Collaborative Optimization and Scheduling of Comprehensive Energy

(181) The combined heat and power (CHP) system of the steel enterprise is shown in FIG. 23. High pressure steam is only supplied to two 25 MW generating sets (ST.sub.3, G.sub.3 and ST.sub.4, G.sub.4), and the rest use medium pressure steam, comprising two 300 MW generating sets (ST.sub.1, G.sub.1 and ST.sub.2, G.sub.2), two 135.sup.T gas boilers (BLR.sub.3 and BLR.sub.4), and two 35.sup.T boilers (BLR.sub.5 and BLR.sub.6). Because the comprehensive energy efficiency of each device is different, it is necessary to achieve the optimal distribution of energy consumption of each set and each boiler to improve energy efficiency and reducing costs.

(182) Considering that the adjustable range of byproduct gas (BFG and COG) and the total steam demand within a forecasting period t.sub.0-t.sub.1 have been determined based on the aforementioned forecasting model, the energy consumption of combined heat and power sets and boilers can be optimally distributed to minimize the energy cost (the sum of cost of self-produced energy (such as BFG, COG and steam) and cost of self-generating electricity (coal price)) and maximize the comprehensive energy efficiency, and the corresponding objective functions are respectively shown in formula (53) and formula (54):

(183) 0$\begin{array}{cc}\min J_{\mathrm{co}}=\underset{i^{\mathrm{mf}}}{\overset{3}{.Math.}}{q}_{i^{\mathrm{mf}}}^u{X}_{i^{\mathrm{mf}}}-q_{\mathrm{power}}^g{P}_{\mathrm{total}}^g-q_{\mathrm{steam}}^g{F}^g& \left(53\right)\end{array}$
max J.sub.effc=(KH.sub.steamF.sup.g+KH.sub.eleP.sub.total.sup.g)/Σ.sub.i.sub.mf.sub.=1.sup.3KH.sub.i.sub.mfX.sub.i.sub.mf  (54)

(184) wherein J.sub.co represents energy cost, g represents generating capacity, and u represents consumption. Since the energy cost of an energy center is calculated independently from that of other plants (iron works, steel works, etc.), the cost of the energy center is taken as the target for optimization in this example

(185) $\underset{i^{\mathrm{mf}}}{\overset{3}{.Math.}}{q}_{i^{\mathrm{mf}}}^u{X}_{i^{\mathrm{mf}}}$
is the sum of consumption costs of byproduct gas and purchased coal, and q.sub.i.sub.mf.sup.u, is the unit price of the (i.sup.mf).sup.th medium (i.sup.mf=1, 2, 3, respectively representing BFG (yuan/m.sup.3), COG (yuan/m.sup.3) and coal (yuan/t)). q.sub.power.sup.g and q.sub.steam.sup.g respectively represent the unit price of power generation (yuan/kwh) and the unit price of steam generation (yuan/t, the price of high pressure steam is the same as that of medium pressure steam) of a heat and power system, and P.sub.total.sup.g and F.sup.g are respectively the total amount of self-generating electricity (kwh) and the total amount of steam generation (t). J.sub.eff is the comprehensive energy efficiency of a heat and power set. KH.sub.steam and KH.sub.ele respectively represent conversion coefficients (kg/m.sup.3 and kg/kwh) of total steam consumption and total power generation into standard coal. KH.sub.i.sub.mf represents the conversion coefficient (kg/m.sup.3 and kg/t) of each energy medium (BFG, COG and coal). The corresponding constraint conditions are shown in formula (55):

(186) $\begin{array}{cc}\{\begin{array}{cc}{X}_1=1/\left(1-{\delta }_B\right)\underset{i_1=1}{\overset{6}{.Math.}}{F}_{\mathrm{BFG},i_1}^{\mathrm{BLR}}& \left(c1\right)\\ X_2=1/\left(1-{\delta }_C\right)\underset{i_1=1}{\overset{6}{.Math.}}{F}_{\mathrm{COG},i_1}^{\mathrm{BLR}}& \left(c2\right)\\ X_3=1/\left(1-{\delta }_M\right)\underset{i_2=1}{\overset{2}{.Math.}}{F}_{M,i_2}^{\mathrm{BLR}}& \left(c3\right)\\ D_{S,i_2}^{\mathrm{BLR}}=f_{i_2}\left(F_{\mathrm{BFG},i_2}^{\mathrm{BLR}},F_{\mathrm{COG},i_2}^{\mathrm{BLR}},F_{M,i_2}^{\mathrm{BLR}}\right),\left(i_2=1,2\right)& \left(c4\right)\\ D_{S,j_1}^{\mathrm{BLR}}=f_{j_1}\left(F_{\mathrm{BFG},j_1}^{\mathrm{BLR}},F_{\mathrm{COG},j_1}^{\mathrm{BLR}}\right),\left(j_1=3,4,5,6\right)& \left(c5\right)\\ \underset{i_2=1}{\overset{2}{.Math.}}{D}_{S,i_2}^{\mathrm{BLR}}=\underset{j_2=1}{\overset{2}{.Math.}}{D}_{T,j_2}^{\mathrm{ST}}& \left(c6\right)\\ \underset{i_3=1}{\overset{4}{.Math.}}{D}_{\mathrm{HS},i_3}^{\mathrm{BLR}}=\underset{j_3=3}{\overset{4}{.Math.}}{D}_{T,j_3}^{\mathrm{ST}}& \left(c7\right)\\ D_{T,j_2}^{\mathrm{ST}}=D_{P,j_2}^{\mathrm{ST}}+D_{S,j_2}^{\mathrm{ST}},\left(j_2=1,2\right)& \left(c8\right)\\ P_{\mathrm{total}}^g=\underset{i_3=1}{\overset{4}{.Math.}}{P}_{i_3}^{\mathrm{ST}}& \left(c9\right)\\ F_M^u=\left(1-{\delta }_S\right)\left(\underset{i_2=1}{\overset{2}{.Math.}}{D}_{S,i_2}^{\mathrm{ST}}+\underset{j_1=3}{\overset{6}{.Math.}}{D}_{S,j_1}^{\mathrm{BLR}}+F^r\right)& \left(c10\right)\\ \underset{i_3=1}{\overset{4}{.Math.}}{D}_{T,i_3}^{\mathrm{ST}}=\left(1-{\delta }_S\right)\underset{j_1=3}{\overset{6}{.Math.}}{D}_{\mathrm{HS},j_1}^{\mathrm{BLR}}& \left(c11\right)\\ F^g=1/\left(1-{\delta }_S\right)\left(F_M^u+\underset{i_3=3}{\overset{4}{.Math.}}{D}_{T,i_3}^{\mathrm{ST}}\right)& \left(c12\right)\\ {\mathrm{Pa}}^L<{\mathrm{Pa}}_0+K_{\mathrm{pa}}{F}^g<{\mathrm{Pa}}^H& \left(c13\right)\end{array}& \left(55\right)\end{array}$

(187) The constraint conditions in formula (55) can be described as follows:

(188) c1-c3: X.sub.1 and X.sub.2 are respectively the total flow of BFG and the total flow of COG input by each boiler, and X.sub.3 is the total consumption of coal. δ.sub.B, δ.sub.C and δ.sub.M are respectively the corresponding loss rates;

(189) c4-c6: the relationship between the input and the output of each device, and ƒ(⋅) is the corresponding fitting relationship;

(190) c7-c8: balance of generation and consumption of high pressure steam and medium pressure steam;

(191) c9: constraints on self-generating capacity;

(192) c10-c12: balance constrains on steam pipe network, F.sup.r is the total recovery amount of waste heat, and δ.sub.s is the loss rate of steam;

(193) c13: constraints on pressure of steam pipe network. Pa.sup.L and Pa.sup.H are respectively the upper limit and the lower limit of the pressure of a steam pipe network in the south area, and Pa.sub.0 represents the initial state of the pressure; and K.sub.pa is the conversion coefficient of flow to pressure;

(194) Then the swarm intelligence optimization method is adopted for solution, thus obtaining a scheduling scheme with optimal economy.

(195) With certain scheduling as an example, the scheduling effects of BFG, COG and LDG systems are shown in FIG. 24(a), FIG. 24(b) and FIG. 24(c) respectively, and the statistics are shown in Table 8. On the premise of satisfying the balance of supply and demand of a steam system and the safe operation of a gas system, a scheduling scheme is given and simulated, and compared with the scheduling based on artificial experience. Results show that the present invention has more reasonable calculation of adjustment amount and lower economic cost of single scheduling. It can be seen that the method of the present invention has good effect on scheduling of the integrated energy system of steel enterprises and can provide a reliable basis for energy conservation, emission reduction, cost reduction and efficiency increase for steel enterprises.

(196) TABLE-US-00008 TABLE 8 Index Statistics after Scheduling Compared with Artificial Experience Cost reduced Whether Whether Scheduling within exceeding exceeding Whether Method 1 h (yuan) safety limit mechanical limit to emit Artificial 0 Yes No Yes experience The present 2835.4 No No No invention

(197) (2) Background Service: Computing Trigger Mechanism with Edge-Cloud Cooperation

(198) The background service is divided into two parts: the cloud service is deployed on the cloud server and receives parameter updating instructions triggered by edge side computing service and human-machine interaction. The background service on the edge side reads and serializes required real-time samples of each energy system and auxiliary information such as production plan messages and production signals from the relational database, takes the processed data samples and relevant information such as plans as input of the “*.iail” model, and forecasts the trend of concerned users in a future period of time to invoke an optimal scheduling model to compute scheduling suggestions and display same on the human-machine interface of the human-computer interface. In general, the background service updates the computed results every 3 minutes to ensure timeliness of data and adaptability of the system.

(199) (3) Human-computer interface: display of computed result and human-machine interaction

(200) The human-computer interface reads the results saved in the database and displays same on the human-machine screen with the refresh frequency of 1 minute. In the case of unplanned temporary changes (such as temporary wind break of blast furnace) in working conditions of production, operators can trigger feedback information to the cloud service through the human-machine interface, and the cloud service will trigger the edge side model to realize self-updating of parameters. Meanwhile, interaction with artificial knowledge data such as forecasting accuracy, scheduling results and scheduling rules through a standardized template can also be realized through an interface of the human-computer interface, and the background service will realize self-updating of the rule base according to artificial knowledge. The practicability and the computing accuracy of the platform are further improved by means of blended learning of human-machine integration.