Uncertain Interval Forecasting for Combined Electricity-Heat-Cooling-Gas Loads in the Integrated Energy System Based on Multi-Task Learning and Multi-Kernel Extreme Learning Machine

: The accurate prediction of electricity-heat-cooling-gas loads on the demand side in the integrated energy system (IES) can provide signiﬁcant reference for multiple energy planning and stable operation of the IES. This paper combines the multi-task learning (MTL) method, the Bootstrap method, the improved Salp Swarm Algorithm (ISSA) and the multi-kernel extreme learning machine (MKELM) method to establish the uncertain interval prediction model of electricity-heat-cooling-gas loads. The ISSA introduces the dynamic inertia weight and chaotic local searching mechanism into the basic SSA to improve the searching speed and avoid falling into local optimum. The MKELM model is established by combining the RBF kernel function and the Poly kernel function to integrate the superior learning ability and generalization ability of the two functions. Based on the established model, weather, calendar information, social–economic factors, and historical load are selected as the input variables. Through empirical analysis and comparison discussion, we can obtain: (1) the prediction results of workday are better than those on holiday. (2) The Bootstrap-ISSA-MKELM based on the MTL method has superior performance than that based on the STL method. (3) Through comparing discussion, we discover the established uncertain interval prediction model has the superior performance in combined electricity-heat-cooling-gas loads prediction.


Introduction
With the rapid development of China's economy and the rapid growth of energy demand, energy development is facing great challenges from resources and the environment. It is a new direction of energy development and reform to construct a clean, low-carbon, safe, and highly efficient energy system to maximize the development and utilization of renewable energy and to improve the efficiency of energy utilization. In this context, the integrated energy system (IES) has become an important direction of the future transformation of China's energy system. The development of the IES with multi-energy and mutual aid centered on electric energy is conducive to further improving the energy efficiency of users and optimizing the energy structure. The accurate forecasting of renewable energy generation on the supply side and multi-energy loads on the demand side is the foundation of the construction and stable operation of the IES.

Research Motivation
With the exploitation and utilization of renewable energy and the rapid development of economy and society, the form of energy supply and demand in China has become more For the electric energy subsystem, it is the most significant subsystem in the IES, interacting with other energy subsystems frequently. The electric energy mainly originates from distributed renewable energy generating units, the electric power network, combined cooling-heat-power (CCHP) units, and energy storage batteries, etc. The consumption of electric energy primarily includes energy storage batteries, electric hydrogen generation system, electric boiler and heat pump, electric refrigerator, electric vehicle charging pile, users electric load, etc. The coupling and transformation between electric energy and other energy sources contain converting electric energy into gas energy through power to gas (P2G) converter, converting electric energy to heat energy through electric heating boiler and heat pump, and converting electric energy to cooling energy through electric refrigerator.
For the heat energy subsystem, the heat energy primarily originates from the exterior heating network, heat storage tank, CCHP units, combined heat and power (CHP) units, electric boiler, heat pump, etc. The consumption of heat energy primarily contains users heat load, absorption refrigerator, heat storage tank, etc. The heat energy subsystem can generate electricity employing waste heat energy produced by the CCHP units and distribute the electricity to the electric energy subsystem. The heat energy subsystem can also use the gas energy generated by the CCHP units to produce electricity and heat energy to satisfy the heat demand.
For the cooling energy subsystem, the cooling energy mainly comes from the CCHP units, cooling energy storage equipment, electric refrigerator, absorption refrigerator, etc. The consumption of cooling energy includes air conditioning, refrigerators, and other refrigeration equipment. Generally, the cooling energy subsystem generates cooling energy through the electric refrigerator employing the electric energy generated by the For the heat energy subsystem, the heat energy primarily originates from the exterior heating network, heat storage tank, CCHP units, combined heat and power (CHP) units, electric boiler, heat pump, etc. The consumption of heat energy primarily contains users heat load, absorption refrigerator, heat storage tank, etc. The heat energy subsystem can generate electricity employing waste heat energy produced by the CCHP units and distribute the electricity to the electric energy subsystem. The heat energy subsystem can also use the gas energy generated by the CCHP units to produce electricity and heat energy to satisfy the heat demand.
For the cooling energy subsystem, the cooling energy mainly comes from the CCHP units, cooling energy storage equipment, electric refrigerator, absorption refrigerator, etc. The consumption of cooling energy includes air conditioning, refrigerators, and other refrigeration equipment. Generally, the cooling energy subsystem generates cooling energy through the electric refrigerator employing the electric energy generated by the electric energy subsystem. It can also use the heat energy from the CCHP units to generate cooling energy through the absorption refrigerator.
For the gas energy subsystem, the source of gas energy mainly includes an exterior gas distribution network, electric hydrogen generation system, gas storage tank, etc. The consumption of gas mainly contains CHP units, CCHP units, a gas storage tank, gas filling station, and users' gas load. The coupling and conversion between gas energy and other energy sources includes converting gas energy into electric energy, heat energy and cooling energy through CHP units and CCHP units. The gas storage tank maintains the balance between the supply and demand of gas by charging and discharging gas energy.
As previous description, there exist strong coupling relations among the electric energy subsystem, the heat energy subsystem, the cooling energy subsystem, and the gas energy subsystem, which makes the IES coordinated operation to a maximum degree. Therefore, in order to predict the combined electricity-heat-cooling-gas loads on the demand side of the IES accurately, the complicated coupling relations are need to be taken into account.

Literature Review
Two factors need to be considered when predicting the combined loads on the demand side of the IES. One is the selection of input variables for the combined load prediction model, and the other is the methods employed to forecast the combined load. As is shown in previous literature, the input variables selected for load forecasting model of a single energy type are greatly different. For example, in terms of the electricity load prediction, historical load data, weather condition, calendar information, and electricity price are considered by Yang et al. [8] to propose a probability density prediction methodology to forecast the electric power load on the basis of the quantile regression of the Gauss process. Guo et al. [9] established a probability density prediction approach based on deep learning mechanism taking atmospheric pollutants factors (e.g., CO, SO2, PM2,5, etc.), date related factors (e.g., seasonal, monthly, daily data), weather related factors (e.g., daily temperatures, rainfall amount, etc.), and economic factors (e.g., economic structure influencing residential living standards) into account to predict electricity consumption. Kim et al. [10] took the weather and holiday related variables as the critical factors to predict peak load demand. Recently, a few studies employed feature selection methods to select the informative input variables for load forecasting model to improve the prediction accuracy. For example, autocorrelation function is employed by Yang et al. [11] to select the optimal input feature to improve the half-hour electricity load prediction accuracy. In terms of the heat load prediction, the weather condition and historical load data are considered by Geysen et al. [12]. Suryanarayana et al. [13] selected the last 7 days of heat load and the next 24 h temperature prediction results as the heat load forecasting input variables. Nigitz et al. [14] considered the ambient temperature for each hour of the day as the critical influencing factor of heat load forecasting. In terms of the cooling load prediction, Zhou et al. [15] considered the historical load data and weather conditions as the input variables for cooling load prediction. Ahmad et al. [16] considered climate related factors as the input variables for the cooling load prediction model. In terms of the gas load prediction, temperature, seasonal index, price of natural gas and oil, exchange rate, as well as the population and gross domestic product (GDP) of the Istanbul province are considered by Beyca et al. [17] to input into the established gas consumption prediction model. Additionally, Shi et al. [18], taking the historical load data, weather condition, economic data, and calendar information into account, established a combined electricity-heat-gas load forecasting model. Tan et al. [19] established a combined electricityheat-cooling-gas load prediction model on the basis of multi-task learning and the least square support vector machine, considering weather-related factors, calendar information, primary economic factors, and technical conditions.
Previous literatures researching on the load forecasting mainly focused on single energy load. However, with the rapid development of computer technology and the construction of the integrated energy systems combined with multi-energy types, the multitask learning (MTL) method gradually replaces the single task learning (STL) method in the application of load forecasting [48]. It is difficult for the traditional STL method to obtain sufficient information of data distribution through training, which results in ignoring the complicated inter-relation among various factors, and hence decreasing the learning effect [49]. Considering about this, the multi-task learning (MTL) method is proposed, which can exactly explore the inter-relations among different tasks [50,51]. In load forecasting applications, Nunes et al. [52] compared the prediction results obtained by the artificial neural networks (ANN) on the basis of the STL and the MTL methods, and then the comparison results demonstrated that the ANN based on the MTL method has a superior performance than that based on the STL method. Shireen et al. [53] combined the MTL method with time series regression and the results showed that the prediction accuracy has been improved effectively. Tan et al. [19] combined the MTL method with the least square support vector machine to apply in electricity-heat-cooling-gas load forecasting, and through comparing with other models, it can be found that the established model is the most time-saving and can obtain the forecasting results with the highest accuracy. In the IES, energy conversion elements (such as CCHP units, electric boilers, electric refrigerators, energy storage, etc.) are used to realize the coupling of electricity, heat, cooling and gas, which will result in strong coupling relationships existing among electricity, heat, cooling and gas sub-systems. In the IESs, a large amount of energy conversion and sharing information of electricity, heat, cooling, and gas exists in the data, and the characteristics of various energy conversion are complicated. Considering this, the MTL method is employed to establish the electricity-heat-cooling-gas loads forecasting model, which can effectively employ the complicated shared information of various energy conversion in the IES, so as to improve the multi-energy load forecasting accuracy.
The forecasting techniques, including traditional statistical methods and AI-based methods, generally belong to the deterministic prediction model. However, the electricityheat-cooling-gas load on the demand side in the IES is a complicated dynamic system affected by weather related factors, calendar information, social-economic factors, and technical conditions with nonlinear, high-dimensional, and uncertain characteristics. Considering the various factors influencing the electricity-heat-cooling-gas loads on the demand side in the IES and the errors that exist in the process of statistics and monitoring, there are uncertainties in the prediction of the combined electricity-heat-cooling-gas loads on the demand side in the IES. Contemporarily, the existing researches mainly focused on single points prediction, which cannot characterize the influence of uncertainties on the combined electricity-heat-cooling-gas loads prediction. Therefore, it is necessary to consider the influence of uncertainties when forecasting the combined electricity-heat-cooling-gas loads on the demand side of the IES, and extend the single point prediction to the uncertainty interval prediction, so as to improve the reliability of the prediction results. In recent years, interval prediction theory is a common method to research on the prediction problem with uncertainties. The prediction results are expressed in the form of interval to represent the uncertainty characteristics of the prediction problem. Compared with the single point prediction method, the interval prediction method can obtain the change range of the predicted object in a period of time in the future by means of interval expression, and quantify the reliability of the prediction model according to the interval width and coverage probability, so as to make the prediction results of the model with high reliability [54]. At present, the commonly used interval prediction methods mainly include mean variance estimation method, upper and lower limit estimation method, delta method, Bayesian method and Bootstrap method [55][56][57][58][59][60][61][62]. The Bootstrap method can approximate the sample distribution characteristics without subjective assumption, but only through repeated sampling, and then can effectively construct the prediction interval, which has the advantages of high accuracy and reliability [61]. Therefore, in this paper, the Bootstrap method is used to establish the electricity-heat-cooling-gas load uncertain interval prediction model.
As previous introduction, extreme learning machine (ELM) method optimized by various intelligent optimization approaches is widely applied in load forecasting. ELM method is a kind of machine learning method established based on a feed-forward neuron network, which is suitable for supervised learning and unsupervised learning problems. The connection weights between input layer and hidden layer and the threshold value of hidden layer neurons are generated randomly. In the training process, the only optimal solution can be obtained by setting the number of hidden layer neurons [63,64]. Compared with the traditional Back-Propagation neural network algorithm, the ELM method has fast learning speed and good generalization performance. However, due to the randomness of ELM feature mapping and its sensitivity to noise, the robustness and generalization ability of the ELM method are reduced. Based on this, Kernel Extreme Learning Machine (KELM) is proposed by introducing the kernel method and regularization coefficient into the ELM method, the influence of random mapping and sample noise on the model can be decreased, and the prediction accuracy of the model can be improved [65,66]. However, the KELM method only considers one kernel form in the prediction, which is difficult to apply to the prediction of the combined electricity-heat-cooling-gas load on the demand side of the IES with multiple data characteristic samples. Therefore, it is necessary to use multi-kernel form to establish a multi-energy load forecasting model. Simultaneously, when employing the multi-kernel extreme learning machine (MKELM) model, it is necessary to determine the relevant parameters in the model. The model parameters set by subjective knowledge always lead to low prediction accuracy. Therefore, this paper uses a new bionic intelligent optimization algorithm named Salp Swarm Algorithm (SSA) to optimize the parameters of the MKELM model. In order to improve the convergence performance of SSA and make it better to solve the high-dimensional complex parameters optimization problem in the optimization process, an improved SSA (ISSA) is proposed in this paper, which introduces the dynamic inertia weight and chaotic local searching mechanism into the basic SSA to improve the searching speed and avoid falling into local optimum.
Above all, this paper established an uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads in the IES based on the MTL method, the Bootstrap method, the MKELM method, and the ISSA.

Our Contributions
Contemporarily, researches on load forecasting of various energy types are relatively mature. With the continuous development of energy reform and energy transformation, the construction of the IES has received great attention, and the accurate prediction of energy on supply and demand sides of the IES is the foundation of the construction of the IES and can ensure the safe and stable operation of the IES. However, it is difficult to find the studies focused on forecasting the combined load based on the coupling relations among integrated energy subsystems. At the aim of filling this gap and exploring the coupling relations deeply, this research establishes a hybrid load prediction model for the IES on the basis of the MTL method, the Bootstrap method, the MKELM, and the ISSA. The contributions of this research are as follows: (1) The complicated coupling relations among the various integrated energy subsystems are untangled systematically, and the significant input variables of the established hybrid load prediction model are determined. (2) The uncertain interval forecasting model of the combined electricity-heat-cooling-gas loads in the IES is established based on the MTL method, the Bootstrap method, the MKELM, and the ISSA. (3) Aiming at solving the problems of local optimum, low convergence speed, and low precision of the basic SSA, the improved SSA (ISSA) is proposed by introducing the inertia weight mechanism and chaotic local searching strategy into the basic SSA. (4) The MKELM model is proposed by combining the RBF kernel function and the Poly kernel function, so as to integrate the superior learning ability and the generalization ability of the two kernel functions.

Organization of the Research
The rest parts of the research are organized as follows. The methodology of the MTL method, the Bootstrap method, the MKELM, and the ISSA as well as the calculation process of the established uncertain interval forecasting model are introduced in Section 2. The numerical analysis is conducted in Section 3. The comparison discussion on the forecasting results obtained by the MTL method and STL method as well as other optimization algorithms and kernel functions of the ELM are provided in Section 4. Section 5 summarizes the primary conclusions of this research.

Materials and Methods
An uncertain interval forecasting model for the combined electricity-heat -cooling-gas loads on the demand side of the IES is established in this paper based on the MTL method, the Bootstrap method, the MKELM, and the ISSA. The established model is not only a combination of several methods, but also makes great improvement for the MKELM and the original SSA.

The Multi-Task Learning Method
The MTL method is a classical inductive transfer mechanism method in the machine learning field. It employs the information hidden in multiple related tasks to improve the generalization ability of algorithms. The MTL method improves the performance of the prediction model by using the sharing mechanism to train multiple tasks in parallel. When studying a problem, it can learn and acquire knowledge of other related problems by using shared information [67]. The information sharing mechanism of the MTL method is shown in Figure 2.

Organization of the Research
The rest parts of the research are organized as follows. The methodology of the method, the Bootstrap method, the MKELM, and the ISSA as well as the calcu process of the established uncertain interval forecasting model are introduced in S 2. The numerical analysis is conducted in Section 3. The comparison discussion o forecasting results obtained by the MTL method and STL method as well as othe mization algorithms and kernel functions of the ELM are provided in Section 4. Sec summarizes the primary conclusions of this research.

Materials and Methods
An uncertain interval forecasting model for the combined electricit -cooling-gas loads on the demand side of the IES is established in this paper based MTL method, the Bootstrap method, the MKELM, and the ISSA. The established mo not only a combination of several methods, but also makes great improvement f MKELM and the original SSA.

The Multi-Task Learning Method
The MTL method is a classical inductive transfer mechanism method in the ma learning field. It employs the information hidden in multiple related tasks to impro generalization ability of algorithms. The MTL method improves the performance prediction model by using the sharing mechanism to train multiple tasks in pa When studying a problem, it can learn and acquire knowledge of other related pro by using shared information [67]. The information sharing mechanism of the method is shown in Figure 2.  In the IES, various energy subsystems, such as electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem, are coupled by energy flow through different energy conversion equipment, including CCHP units, electric boiler, gas boiler, electric refrigerator, absorption refrigerator, etc. According to the description in Section 1.2, there exist a strong coupling relation among different energy subsystems. A large amount of energy conversion sharing information exists in the data, which are difficult to be extracted and summarized by traditional manual feature extraction methods. Therefore, the MTL method is combined into the multi-energy load forecasting model, which can effectively employ the complicated sharing information of electricity, heat, cooling, gas, and other energy conversion, help machine learning and bionic intelligent optimization algorithms to train better abstract characteristics, and improve the forecasting precision of multi-energy loads in the IES [19].
Suppose there are T target tasks, and the data of each task originates from the generalized function space X × Y, X ∈ R n , Y ∈ R. Each task contains m training sample sets {(x 1t , y 1t ), (x 2t , y 2t ), · · · , (x mt , y mt )}. Each task has a sample distribution P t , and each task is different but related.
In the MTL method, the goal of the multi-tasks is to learn is the training function of T tasks, and y it is the training sample.
In order to balance the correlation and constraint among different tasks, it is necessary to consider the learning optimization goals among different tasks to obtain better learning effects [68]. There are differences among tasks, but the association among tasks can be realized by sharing information path.

The Bootstrap Method
The Bootstrap method is an important statistical method in non-parametric statistics to estimate the variance of statistics and then to evaluate the interval. It is a statistical method based on repeated sampling. The prediction interval established by this method can not only accurately represent the uncertainty degree of the prediction object, but also analyze the change range of the time series of the prediction object in a period of time in the future, which can effectively make up for the limitations of the single point prediction method [69]. In the combined electricity-heat-cooling-gas loads forecasting, considering the characteristics, various uncertainties will affect the forecasting accuracy of the prediction model, while combining the Bootstrap method into the multi-energy load forecasting model to predict the interval of multi-energy load can improve the prediction precision and make the forecasting results better applied to engineering practice.
Supposing the sample of the combined electricity-heat-cooling-gas load prediction is x i is the set of influencing factors affecting the combined electricity-heatcooling-gas load, t i is the actual values of the electricity-heat-cooling-gas load, therefore, where f (x i ) is the forecasting values of the electricity-heat-cooling-gas load, f (x i ) is the actual regression values of the electricity-heat-cooling-gas load forecasting model, and ε(x i ) is the data noise [70].
In the Equation (2), t i − f (x i ) is the total error between the actual value and the forecasting value of the electricity-heat-cooling-gas load, and f (x i ) − f (x i ) is the prediction model error of the electricity-heat-cooling-gas load. The total error consists of prediction model error and data noise. Assuming that all errors are unbiased, and then the variables can be estimated according to the variance.
The Bootstrap method can obtain M sub-samples by repeatedly sampling the electricityheat-cooling-gas load samples, and then train the electricity-heat-cooling-gas load forecasting model, and estimate the model error variance according to Equation (3).
where σ 2 ∆ (x) is the error variance of the forecasting model, representing the model uncertainties caused by human subjective factors such as the structure and the parameters setting of the forecasting model. f (x) is the expected values of the output of M prediction models, and f k (x) is the prediction value of the output of the k-th model.
In the Bootstrap method, the data noise variance represents the uncertainties of the electricity-heat-cooling-gas loads samples caused by random noise, and can be estimated according to Equation (4) where σ 2 ε (x) is the data noise variance, and E represents the expected value. At the aim of effectively predicting the data noise variance σ 2 ε (x), we need to construct the square residual sequence of the prediction model, and combine it into a new sample set , and then train the prediction model according to the new sample set. Among them, r i can be obtained through Equation (5).
When training the new sample set, in order to maximize the probability of noise variance appearing in the electricity-heat-cooling-gas load samples, the maximum likelihood method is used to train the new model, so as to ensure the prediction confidence level. The objective function is After the training of the forecasting model, the forecasting interval of the electricityheat-cooling-gas load with significance level α can be constructed as follows: where 1 − α/2 is the quantile in the standard Gaussian distribution [71].

The Multi-Kernel Extreme Learning Machine
The MKELM method transfers the low dimensional data into Hilbert space to reduce the randomness and complexity of the model network through regularization theory and kernel mapping based on the extreme learning machine (ELM) and kernel extreme learning machine (KELM), which can effectively improve the prediction accuracy and generalization ability of the model.
The ELM is a new machine learning method proposed by Professor Huang of Nanyang Polytechnic University in Singapore [72]. Its structure is shown in Figure 3. The ELM is an effective single hidden layer feedforward network, which only needs to set the number of hidden layer nodes in the network structure without adjusting the network input weights and bias values. Compared with neural network and other methods, the ELM method has superior prediction performance and better generalization ability in classification and regression prediction.
Given N training samples, where (x i , t i ) is the input value of the electricity-heatcooling-gas loads samples, x in ] is the expected output value of the electricity-heat-cooling-gas loads samples, and the amount of hidden nodes of the ELM model is L, then the output should be where β is the output weight from hidden node to output node, G represents the motivation function, a i and b i are the parameters of i-th hidden layer, and h x j is the output vector of the hidden layer.
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 33 number of hidden layer nodes in the network structure without adjusting the network input weights and bias values. Compared with neural network and other methods, the ELM method has superior prediction performance and better generalization ability in classification and regression prediction.
is the expected output value of the electricity-heat-cooling-gas loads samples, and the amount of hidden nodes of the ELM model is L, then the output should be where β is the output weight from hidden node to output node, G represents the mo- Among Equation (9), Furthermore, it can be expressed by matrix expression of linear system as Among Equation (9), In the ELM model, if the motivation function is infinitely differentiable, the bias value and input weight of the structure do not need to be changed, and then the ELM training problem can be transformed into solving β. The least square method is usually used to determine the output weight of the ELM model: where H + is the Moore Penrose generalized inverse of the hidden layer output matrix H [73].
When solving H + , the sample data set may have multi-collinearity problem, which may lead to non-singularity, and then have a certain impact on the prediction results. Considering this, the parameter 1/C is introduced into diagonal matrix HH T to make the characteristic root of HH T deviate from zero. Therefore, the weights β of the ELM model can be solved through Equation (12) [74].
The ELM model can effectively avoid the problems of slow calculation speed and overfitting of a traditional neural network, since it does not need to train a large number of model parameters in the training process. However, there also exist disadvantages in the ELM model. For example, the process of selecting the number of hidden nodes is complicated and time-consuming, which needs a lot of experiments or subjective experience to select the optimal number of neurons resulting in the decline of prediction performance and generalization ability of the model.
Aiming at solving this demerit, Professor Huang proposed the KELM model, which can effectively avoid the selection of hidden layer. By transforming the output vector of hidden layer into kernel function mapping of support vector machine, the output random fluctuation is reduced, and the generalization ability and robustness of the model are improved [72].
The definition of the kernel matrix in the KELM model is where Ω ELM is the kernel matrix of the KELM model and K(x i , x j ) is the kernel function. Then, Equation (12) can be transformed into Therefore, the output expression of the KELM model can be written as [74] f Generally, the KELM model provides various selections of kernel functions, such as Poly kernel, Gaussian RBF kernel, etc. The selection of the kernel function has a great impact on the forecasting results of the combined electricity-heat-cooling-gas loads. The single kernel model is usually difficult to adapt to the samples with multiple data characteristics. Hence, to overcome the shortcomings and limitations of single KELM, this paper combines different kernel functions and makes full use of the advantages of different kernel functions in learning and generalization, and constructs a multi-kernel extreme learning machine model.
According to Equation (15), the selection of kernel function has a significant influence on the prediction results of the model. Kernel functions are divided into global kernel (such as RBF kernel function) and local kernel (such as poly kernel function), which are where δ is the kernel parameter of the RBF kernel function and d is the kernel parameter of the Poly kernel function [75,76]. The data characteristic distribution of RBF kernel function and Poly kernel function when x = 0 and x = 0.5 are provided in Figures 4 and 5. From Figure 4, it can be found that the RBF kernel function can make the near data affect the kernel value, and has better local learning ability, but it has the disadvantage of poor generalization performance. The Poly kernel function can identify the influence of remote data on kernel value, and its global generalization ability is strong. However, when the data is concentrated to the central point, the influence of remote data on kernel value is significantly reduced, and its learning ability becomes weak. Therefore, it can be seen that the identification ability of different kernel functions for sample data features is quite different. As the integrated embodiment of multiple influencing factors, the electricity-heat-cooling-gas load is a dynamic evolution process of multiple data features including trend characteristics, periodic characteristics, volatility characteristics, and chaotic characteristics. Thus, it is difficult to use a single kernel function to make a high-precision regression prediction for the electricity-heat-cooling-gas load with multiple data distribution characteristics. Considering this, this paper combines the RBF kernel function and the Poly kernel function to construct a hybrid kernel function, so as to integrate the superior learning ability and generalization ability of the two kernel functions. Based on the Mercer theory, the kernel matrix satisfies semi-positive definiteness, which means the linear combination of multiple single kernel functions is still a kernel function. Based on this, the hybrid kernel function is constructed as follows where ρ is the weight coefficient of   The MKELM model can be constructed by introducing the Equation (18) into the Equation (15). The parameters need to be determined including regularization coefficient C, RBF kernel parameter δ , Poly kernel parameter d , and weight coefficient ρ . In this   The MKELM model can be constructed by introducing the Equation (18) into the Equation (15). The parameters need to be determined including regularization coefficient C, RBF kernel parameter δ , Poly kernel parameter d , and weight coefficient ρ . In this  The MKELM model can be constructed by introducing the Equation (18) into the Equation (15). The parameters need to be determined including regularization coefficient C, RBF kernel parameter δ, Poly kernel parameter d, and weight coefficient ρ. In this paper, the improved SSA bionic intelligent optimization algorithm is used to optimize these four parameters.

The Improved SSA
The Salp Swarm Algorithm (SSA) is a new heuristic optimization algorithm proposed by an Australian scholar named Seyedali Mirjalili in 2017 [77]. It is proposed by imitating the behavior of the chain group of the salp in the process of looking for food. A single salp and a salp chain are shown in Figure 6.
1, 9, x FOR PEER REVIEW 14 of 33 paper, the improved SSA bionic intelligent optimization algorithm is used to optimize these four parameters.

The Improved SSA
The Salp Swarm Algorithm (SSA) is a new heuristic optimization algorithm proposed by an Australian scholar named Seyedali Mirjalili in 2017 [77]. It is proposed by imitating the behavior of the chain group of the salp in the process of looking for food. A single salp and a salp chain are shown in Figure 6. In the SSA, the initial population is divided into leader and follower. The leader leads the Salp chain. The specific steps of the SSA are as follows.
(1) Parameters setting The position of the initial population in the SSA is set randomly, and the position matrix of the salp swarm S is as follows 11 12 1 where ij s represents the value of the j-th variable of the i-th salp, 1, 2, , ij s can be calculated by Equation (20) with random distribution: In the SSA, the initial population is divided into leader and follower. The leader leads the Salp chain. The specific steps of the SSA are as follows.
(1) Parameters setting Five parameters need to be set in the SSA, which are the initial population number Agents_no, the amount of variables dim, the maximum number of iteration Max_iteration, the lower bound of variables lb = [lb 1 , lb 2 , · · ·], and the upper bound of variables ub = [ub 1 , ub 2 , · · ·].
(2) Population initialization The position of the initial population in the SSA is set randomly, and the position matrix of the salp swarm S is as follows where s ij represents the value of the j-th variable of the i-th salp, i = 1, 2, · · · , n, j = 1, 2, · · · , d. s ij can be calculated by Equation (20) with random distribution: where S(i, j) demonstrates the value of the i-th row and j-th column in the matrix s ij , rand(i, j) implies the random matrix, and ub(i) and lb(i) mean the upper bound and the lower bound of the i-th salp, respectively. (

3) Fitness function construction
The fitness function is written as follows: In the matrix OS, the position of the salp with the best fitness value is designated as the variable F as the source food, which is tracked by the salp chain.

(4) Iteration
In the SSA, all the salp will change their positions to perform global search and avoid local optimization. Leaders update the location of food sources as shown below: where x 1 j represents the position of the first salp (leader) in the j-th dimension, F j is the position of the food source, ub j and lb j demonstrate the upper bound and lower bound, respectively, and c 1 , c 2 , and c 3 are random constants. c 2 and c 3 are generated randomly in the interval [0, 1], which determine the iteration direction of the next position of the j-th dimension and step length. c 1 is a critical parameter in the SSA, and the definition is as below: where l is the current number of iterations and L is the maximum number of iterations. The SSA uses Newton's law of motion to update the follower's position as where x i j is the position of the i-th follower in the j-th dimension, t means the time, v 0 represents the initial velocity and i ≥ 2.
In the SSA, all steps are executed in the iteration process until reaching the end standard of the SSA iteration, except initialization.
In order to further improve the optimization performance of the SSA and make it better for solving the optimization problem of high-dimensional complex parameters in the MKELM prediction model, this paper introduces dynamic inertia weight and chaotic local searching mechanism into the basic SSA intelligent optimization algorithm, which named improved SSA (ISSA) in this paper, so as to improve the optimization searching speed of the SSA and avoid falling into local optimum easily.
It can be seen from Equation (22) that the leader updates the position according to the information position of the parent generation in the iteration process, which mainly focuses on the area near the existing individuals.
As a result, the search scope of the algorithm is limited, leading to the oscillation phenomenon in the later optimization stage, which will reduce the convergence accuracy and convergence speed of the SSA. In order to solve this problem, inertia weight mechanism is introduced into the SSA, which can control the balance between global searching and local searching of the whole algorithm by adjusting its own value, so as to strengthen the convergence ability of the algorithm. Therefore, the inertia weight mechanism is introduced on the basis of the Equation (22), and the updated equation is obtained as below: where w is the inertia weight, which can control the influence of parent position information on the new population, so as to adjust the searching trend and expand the searching scope. Generally speaking, using a larger weight in the early stage of the iteration will help to quickly locate the range of the optimal solution, and the convergence speed of the algorithm can be accelerated by reducing the weight in the later stage of the iteration. It is difficult to make optimal adjustment according to the population optimization and iteration by using fixed inertia weight. Therefore, the dynamic inertia weight with nonlinear time-varying updating is employed in this paper to improve the basic SSA.
where t is the number of iteration steps of the current population.
With the increase of iteration steps, the inertia weight between individuals decreases from the initial value, which can avoid the oscillation phenomenon of the SSA in the later stage of optimization and accelerate the population convergence speed of the SSA.
Generally speaking, when the SSA enters the later stage of the optimization iteration, the population tends to approach the local optimum. Chaotic motion has the advantages of regularity and ergodicity, which can traverse all states in the specified range without repetition. Therefore, in order to further enhance the optimization performance of the SSA, the chaotic local searching strategy is introduced into the SSA to reduce the shortcoming of easily falling into local optimum and improve the searching efficiency of the SSA in the solution space. The specific process of introducing chaotic local searching strategy into the SSA is as follows.
(1) The individual s ij is generated into chaotic variables as where ϕ j i is the chaotic variable, s min i and s max i are the minimum value and maximum value of the j-th variable.
(2) Generating chaotic sequence by logic mapping as (3) The chaotic sequence is transformed into the searching space of the target solution to generate a new position, and compared with the group optimal individual. The better value is selected to update, so as to complete the local searching process. This procedure can be written as The improved SSA has been enhanced in global searching and local searching, hence it has higher convergence accuracy and better optimization performance. The improved SSA can also help the MKELM prediction model to explore the optimal parameters values to improve the accuracy in the combined electricity-heat-cooling-gas loads prediction.

The Established Uncertain Interval Forecasting Model for the Combined Electricity-Heat-Cooling-Gas Loads in the IES
The established uncertain interval forecasting model of the combined electricity-heatcooling-gas loads on the demand side in the IES contains the Bootstrap method, the MTL method, and the MKELM optimized by the improved SSA. The role of each method is illustrated in Figure 7. The specific modeling steps of the uncertain interval forecasting model are described as below.
The improved SSA has been enhanced in global searching and local searching, hence it has higher convergence accuracy and better optimization performance. The improved SSA can also help the MKELM prediction model to explore the optimal parameters values to improve the accuracy in the combined electricity-heat-cooling-gas loads prediction.

The Established Uncertain Interval Forecasting Model for the Combined Electricity-Heat-Cooling-Gas Loads in the IES
The established uncertain interval forecasting model of the combined electricity-heat-cooling-gas loads on the demand side in the IES contains the Bootstrap method, the MTL method, and the MKELM optimized by the improved SSA. The role of each method is illustrated in Figure 7. The specific modeling steps of the uncertain interval forecasting model are described as below. (1) The selection of model variables and the preprocessing of sample data According to the characteristics of the combined electricity-heat-cooling-gas loads forecasting, the input variables of the forecasting model are selected, which are the primary factors influencing the demand of the electricity-heat-cooling-gas loads, including weather, calendar information, social-economic factors, and historical load. Considering that each variable has different units, the sample data of the selected variables are normalized as the Equation (30), so that they can be analyzed in the same dimension.
where ˆi x is the standardized value, i x is the sample value, and min(x) and max(x) are the minimum value and maximum value, respectively.
(2) Sample division The sample data is classified into training sample data and testing sample data to construct the training set and testing set for the combined electricity-heat-cooling-gas loads forecasting model. Additionally, then the Bootstrap method is employed to construct M pseudo sample sets by repeated sampling. (1) The selection of model variables and the preprocessing of sample data According to the characteristics of the combined electricity-heat-cooling-gas loads forecasting, the input variables of the forecasting model are selected, which are the primary factors influencing the demand of the electricity-heat-cooling-gas loads, including weather, calendar information, social-economic factors, and historical load. Considering that each variable has different units, the sample data of the selected variables are normalized as the Equation (30), so that they can be analyzed in the same dimension.
wherex i is the standardized value, x i is the sample value, and min(x) and max(x) are the minimum value and maximum value, respectively.
(2) Sample division The sample data is classified into training sample data and testing sample data to construct the training set and testing set for the combined electricity-heat-cooling-gas loads forecasting model. Additionally, then the Bootstrap method is employed to construct M pseudo sample sets by repeated sampling.

(3) Multi-tasks construction
There exist strong coupling relations among electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem in the IES as description in Section 2. The load demand of each energy subsystem is related and influenced each other, which is suitable to employ the MTL method. Therefore, the MTL method is incorporated into the uncertain interval forecasting model for the electricity-heat-cooling-gas loads. Hence, four single tasks can be constructed according to the load demand characteristics of electricity, heat, cooling, and gas, which can be fed into the model as a module in the prediction model.

(4) Parameters setting of the ISSA and the optimal values for the parameters in the MKELM
The parameters of the ISSA algorithm are set, and the regularization coefficient C, the RBF kernel parameter δ, the poly kernel parameter d and the weight coefficient ρ in the MKELM model are optimized by using the ISSA, and the optimal values of these four parameters are finally determined. The root mean square error (RMSE) of the prediction results is utilized as the fitness function of the improved SSA, which can be calculated as Equation (31).
Additionally, the group position of the ISSA is updated according to Equation (25). The chaotic local searching is carried out for the population optimal solution in a single iteration employing Equations (27)- (29). When the termination condition is reached, the optimization of the ISSA is completed, and then the optimal values of four parameters in the MKELM model can be obtained.  The performance of the proposed uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads is evaluated by two indicators, which are the prediction interval coverage probability (PICP) and the prediction interval normalized averaged width (PINAW). The PICP represents the probability that the forecast interval covers the actual electricity-heat-cooling-gas loads. The larger the PICP, the higher the reliability of the interval obtained by the prediction model is, which can better avoid the electricity-heat-cooling-gas loads from exceeding the upper and lower limits of the prediction interval. The PINAW demonstrates the uncertainty degree of the prediction interval of the model. The narrower the width, the lower the uncertainty of the prediction interval and the better the accuracy of the model. The calculation process of the PICP and PINAW are shown in Equations (32) and (33). (33) where N is the number of testing samples for the electricity-heat-cooling-gas loads prediction, R is the range of the target values defined as the difference between its minimum and maximum values, U(x i ) and L(x i ) are the upper and lower bounds of the i-th prediction interval of the electricity-heat-cooling-gas loads, and b i = 1, if the actual load value is contained in the prediction interval, otherwise, b i = 0. The calculation process of the established uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads in the IES is illustrated in Figure 8.
The calculation process of the established uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads in the IES is illustrated in Figure 8.

Results
Aiming at verifying the feasibility and effectiveness of the established uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads in the IES, the industrial park in northern China is selected as the research object for empirical analysis. The internal energy supply and demand structure of this industrial park is illustrated in Figure 1. The IES of the selected industrial park is composed of electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem, including CCHP unit, Photovoltaic unit, electric boiler, gas boiler, electric refrigerator, absorption refrigerator, other energy generation and conversion equipment, as well as electric energy storage, heat storage and other energy storage equipment to improve the operation flexibility and operation economy of the IES.

Results
Aiming at verifying the feasibility and effectiveness of the established uncertain interval forecasting model for the combined electricity-heat-cooling-gas loads in the IES, the industrial park in northern China is selected as the research object for empirical analysis. The internal energy supply and demand structure of this industrial park is illustrated in Figure 1. The IES of the selected industrial park is composed of electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem, including CCHP unit, Photovoltaic unit, electric boiler, gas boiler, electric refrigerator, absorption refrigerator, other energy generation and conversion equipment, as well as electric energy storage, heat storage and other energy storage equipment to improve the operation flexibility and operation economy of the IES.

Variables Selection and Data Preprocess
In the application of the established uncertain interval forecasting model for the electricity-heat-cooling-gas loads in the IES, four factors, namely, weather, calendar information, social-economic factors, and historical load, are considered as the input variables of the forecasting model. The specific variables of four factors are listed in Figure 9. As can be seen from Figure 9, the weather factor contains the temperature index and humidity index. The calendar information factor includes workday and holiday (weekends included). The social-economic factor is the critical factor of the combined load prediction of the IES which is consisted of electricity price, heat energy price, gas price, and per capita GDP.
The historical load is made up of historical load data of each energy subsystem, including historical electric load, historical heat load, historical cooling load, and historical gas load.
In the application of the established uncertain interval forecasting model for the electricity-heat-cooling-gas loads in the IES, four factors, namely, weather, calendar information, social-economic factors, and historical load, are considered as the input variables of the forecasting model. The specific variables of four factors are listed in Figure 9. As can be seen from Figure 9, the weather factor contains the temperature index and humidity index. The calendar information factor includes workday and holiday (weekends included). The social-economic factor is the critical factor of the combined load prediction of the IES which is consisted of electricity price, heat energy price, gas price, and per capita GDP. The historical load is made up of historical load data of each energy subsystem, including historical electric load, historical heat load, historical cooling load, and historical gas load. The hourly data of different variables from January 2018 to June 2018 are selected as the training sample, and the hourly data of different variables of July 2018 is treated as the test sample for the established uncertain interval forecasting model of the industrial park's combined electricity-heat-cooling-gas loads. The data are collected from the Meteorological Bureau and Energy statistics of the province where the industrial park located and through site investigation. Additionally, Matlab R2014a is employed for programming.

Forecasting Results Analysis
Before forecasting the electricity-heat-cooling-gas loads in the IES of the industrial park, it is necessary to determine the parameters of the ISSA. Specifically, Agents_no = 100, dim = 4, Max_iteration = 200, and the lower and upper bounds of variables are 0 and 10. This paper selects a workday (2018.07.12) and a holiday (2018.07.22) as the testing objects. The ISSA is employed to optimize the parameters of the MKELM, including regularization coefficient C, RBF kernel parameter δ , Poly kernel parameter d , and weight coefficient ρ , and the optimization results are shown in Table 1.  Figure 10 shows the iterative optimization convergence curve of the ISSA and the basic SSA in finding the optimal parameters values of the MKELM model. It can be seen that the ISSA converges in steps 50 and 41 on workday and holiday, respectively, while the basic SSA converges in steps 69 and 54, respectively. Moreover, the fitness function values of the ISSA in workday and holiday are 0.054 and 0.038, respectively, which are The hourly data of different variables from January 2018 to June 2018 are selected as the training sample, and the hourly data of different variables of July 2018 is treated as the test sample for the established uncertain interval forecasting model of the industrial park's combined electricity-heat-cooling-gas loads. The data are collected from the Meteorological Bureau and Energy statistics of the province where the industrial park located and through site investigation. Additionally, Matlab R2014a is employed for programming.

Forecasting Results Analysis
Before forecasting the electricity-heat-cooling-gas loads in the IES of the industrial park, it is necessary to determine the parameters of the ISSA. Specifically, Agents_no = 100, dim = 4, Max_iteration = 200, and the lower and upper bounds of variables are 0 and 10. This paper selects a workday (2018.07.12) and a holiday (2018.07.22) as the testing objects. The ISSA is employed to optimize the parameters of the MKELM, including regularization coefficient C, RBF kernel parameter δ, Poly kernel parameter d, and weight coefficient ρ, and the optimization results are shown in Table 1.  Figure 10 shows the iterative optimization convergence curve of the ISSA and the basic SSA in finding the optimal parameters values of the MKELM model. It can be seen that the ISSA converges in steps 50 and 41 on workday and holiday, respectively, while the basic SSA converges in steps 69 and 54, respectively. Moreover, the fitness function values of the ISSA in workday and holiday are 0.054 and 0.038, respectively, which are better than those of the basic SSA (0.079 and 0.047, respectively). Therefore, compared with the basic SSA, the ISSA has higher convergence accuracy and convergence speed. The main reason is that the ISSA can effectively solve the oscillation problem of the basic SSA in the iterative convergence process by introducing the dynamic inertia weight mechanism, thus improving the stability and searching efficiency of the algorithm. At the same time, the chaotic local searching mechanism is introduced into the basic SSA, which can reduce the risk of the SSA falling into local optimum, and improve the optimization performance of the SSA.
Furthermore, according to the established uncertain interval forecasting model of the combined electricity-heat-cooling-gas loads in the IES, the forecasting intervals of various loads on workday and holiday at 90% confidence level can be obtained, as shown in Figures 11 and 12. with the basic SSA, the ISSA has higher convergence accuracy and convergence speed. The main reason is that the ISSA can effectively solve the oscillation problem of the basic SSA in the iterative convergence process by introducing the dynamic inertia weight mechanism, thus improving the stability and searching efficiency of the algorithm. At the same time, the chaotic local searching mechanism is introduced into the basic SSA, which can reduce the risk of the SSA falling into local optimum, and improve the optimization performance of the SSA. Furthermore, according to the established uncertain interval forecasting model of the combined electricity-heat-cooling-gas loads in the IES, the forecasting intervals of various loads on workday and holiday at 90% confidence level can be obtained, as shown in Figures 11 and 12. (a) Electric load (b) Heat load same time, the chaotic local searching mechanism is introduced into the basic SSA, which can reduce the risk of the SSA falling into local optimum, and improve the optimization performance of the SSA. Furthermore, according to the established uncertain interval forecasting model of the combined electricity-heat-cooling-gas loads in the IES, the forecasting intervals of various loads on workday and holiday at 90% confidence level can be obtained, as shown in Figures 11 and 12.  (c) Cooling load (d) Gas load According to the above prediction results, it can be deduced that: (1) Electric, heat, cooling, and gas loads have peak, flat, and valley periods, but the number of peak, flat, and valley is different. The prediction results of the uncertain interval forecasting model established in this paper can better reflect the fluctuation According to the above prediction results, it can be deduced that: (1) Electric, heat, cooling, and gas loads have peak, flat, and valley periods, but the number of peak, flat, and valley is different. The prediction results of the uncertain interval forecasting model established in this paper can better reflect the fluctuation characteristics of various loads. However, generally, the prediction interval is relatively large in the peak and valley periods, which indicates that the load is often difficult to predict at the time point with large fluctuation. (2) The electric load, heat load, and gas load fluctuate greatly, but the prediction interval of electric load is large. The PINAW values of electric load prediction on a workday (July 12) and holiday (July

Discussion
In order to verify the effectiveness and forecasting accuracy of the established uncertain interval prediction model, including the MTL method, the Bootstrap method, the MKELM method, and the ISSA, for the combined electricity-heat-cooling-gas loads in the IES, the STL method as well as some other optimization algorithms and other kernel functions of the ELM are selected to comparatively discuss about the prediction accuracy.

Comparison Discussion on the Forecasting Results Obtained by the MTL Method and STL Method
There exist strong coupling relations among electric, heat, cooling, and gas loads in the IES, therefore, the MTL method is suitable for the uncertain interval prediction multi-energy loads. According to the prediction results in the previous section, it can be seen that the multi-load uncertain interval prediction model containing the MTL method has superior prediction performance. In order to further verify the effectiveness and feasibility of the MTL method, the STL method is employed to replace the MTL method in the established uncertain interval prediction model, hence we can comparatively analyze the prediction results obtained by the STL method and the MTL method. In order to ensure the comparability of the prediction results, the same input sample data, output sample data, and the model parameter values are used for the two kinds of prediction models with different task types. The difference is that in the multi-load uncertain interval prediction model of the IES with the STL method, electric load, heat load, cooling load, and gas load will be input into the model for prediction, respectively, while in the model with the MTL method, electric load, heat load, cooling load, and gas load will be input into the model at the same time, and the uncertain interval prediction results of the four kinds of loads can be obtained by one time.
The prediction results of the multi-load uncertain interval prediction model of the IES with the STL method on workday (July 12) and holiday (July 22) are shown in Figures 13 and 14.
The prediction performance of the multi-load uncertain interval prediction model for the IES of the industrial park with the STL method and the MTL method is shown in Table 2. It can be found that the PICP values of the model with the STL method and the MTL method are all 100% coverage of electric load, heat load, and cooling load on the workday and holiday. The PICP values of the model with the MTL method for gas load are both 100% on the workday and holiday, while the PICP values of the model with the STL method for gas load on the workday (July 12) and holiday (July 22) are 95.83% and 93.75%, respectively, which indicates that two results points and three results points in the prediction results of the model with the STL method do not fall into the prediction range, respectively. According to PINAW values, the prediction range width of the model with the MTL method on the workday and holiday are both narrower than that of the model with the STL method on workday and holiday, which demonstrates that the uncertainty of the prediction interval obtained by the MTL method is lower and the forecasting accuracy of the MTL method is higher than the STL method. This is primarily because the MTL mechanism can comprehensively consider the complicated coupling relations among electric, heat, cooling, and gas loads in the IES, and make full use of the information of electric, heat, cooling, and gas loads. When one load is predicted, the information contained other related loads will be utilized, so as to improve the prediction accuracy of combined electricity-heat-cooling-gas loads.   The prediction performance of the multi-load uncertain interval prediction model for the IES of the industrial park with the STL method and the MTL method is shown in Table 2. It can be found that the PICP values of the model with the STL method and the MTL method are all 100% coverage of electric load, heat load, and cooling load on the workday and holiday. The PICP values of the model with the MTL method for gas load are both 100% on the workday and holiday, while the PICP values of the model with the STL method for gas load on the workday (July 12) and holiday (July 22) are 95.83% and 93.75%, respectively, which indicates that two results points and three results points in the prediction results of the model with the STL method do not fall into the prediction range, respectively. According to PINAW values, the prediction range width of the model with the MTL method on the workday and holiday are both narrower than that of the model with the STL method on workday and holiday, which demonstrates that the uncertainty of the prediction interval obtained by the MTL method is lower and the forecasting accuracy of the MTL method is higher than the STL method. This is primarily because the MTL mechanism can comprehensively consider the complicated coupling relations among electric, heat, cooling, and gas loads in the IES, and make full use of the information of electric, heat, cooling, and gas loads. When one load is predicted, the information contained other related loads will be utilized, so as to improve the prediction accuracy of combined electricity-heat-cooling-gas loads.   Additionally, we also compare the calculation time between the model with the MTL method and the model with the STL method, as illustrated in Figure 15. When forecasting the electricity-heat-cooling-gas loads, the model with the MTL method takes 795 s, while the prediction model with the STL method takes 514 s, 479 s, 415 s, and 497 s to forecast electric load, heat load, cooling load, and gas load, respectively. The total time consuming of the model with the STL method is 1905 s. It can be seen that the prediction model with the MTL method takes much less time than the prediction model with the STL method. The primary reason is that the MTL method only trains one time to forecast the loads of electricity, heat, cooling, and gas, while the STL method needs to train four times. Therefore, the MTL mechanism effectively shortens the training time of the model, and has strong practical engineering value and significance.
MTL method and the model with the STL method, as illustrated in Figure 15. When forecasting the electricity-heat-cooling-gas loads, the model with the MTL method takes 795 s, while the prediction model with the STL method takes 514 s, 479 s, 415 s, and 497 s to forecast electric load, heat load, cooling load, and gas load, respectively. The total time consuming of the model with the STL method is 1905 s. It can be seen that the prediction model with the MTL method takes much less time than the prediction model with the STL method. The primary reason is that the MTL method only trains one time to forecast the loads of electricity, heat, cooling, and gas, while the STL method needs to train four times. Therefore, the MTL mechanism effectively shortens the training time of the model, and has strong practical engineering value and significance.  The PICP and PINAW values of the proposed model and five comparison models are shown in Figure 16 taking the electric load, heat load, cooling load, and gas load on a workday and holiday as an example. It can be seen that according to the PICP index, the interval coverage rate of the proposed model in this paper integrating the MTL method, The PICP and PINAW values of the proposed model and five comparison models are shown in Figure 16 taking the electric load, heat load, cooling load, and gas load on a workday and holiday as an example. It can be seen that according to the PICP index, the interval coverage rate of the proposed model in this paper integrating the MTL method, the Bootstrap method, the ISSA, and the MKELM method is 100%, which can satisfy the requirement of 90% confidence level, followed by MTL-Bootstrap-SSA-MKELM, MTL-Bootstrap-ISSA-rELM, MTL-Bootstrap-ISSA-pELM, MTL-Bootstrap-SSA-ELM, and MTL-Bootstrap-ELM. This result demonstrates that the prediction model with the MKELM method and the ISSA has the highest interval coverage rate. From the PINAW index, the interval width of the proposed model in this paper integrating the MTL method, the Bootstrap method, the ISSA, and the MKELM method is the narrowest, which are 1.42 and 4.02 on the workday and holiday, respectively, followed by MTL-Bootstrap-SSA-MKELM, MTL-Bootstrap-ISSA-rELM, MTL-Bootstrap-ISSA-pELM, MTL-Bootstrap-SSA-ELM, and MTL-Bootstrap-ELM. This comparison result indicates that the interval width of the prediction model optimized by the ISSA is narrower than that optimized by the basic SSA and not optimized by the SSA, and the interval width of the MKELM model is narrower than that of (Poly or RBF) single kernel ELM model. This indicates that the proposed ISSA method and the MKELM method in this paper are effective and feasible in the uncertain interval prediction of the combined electricity-heat-cooling-gas loads in the IES.

Conclusions
Based on the analysis of the complex coupling relationship among the electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem, it can be discovered that there exists strong coupling relationships among various sub-energy systems. Therefore, it is necessary to forecast the combined electrici- Above all, the established uncertain interval prediction model MTL-Bootstrap-ISSA -MKELM in this paper has the superior performance in combined electricity-heat -coolinggas loads uncertain interval prediction, of which the interval prediction effect is superior than other five comparison models.

Conclusions
Based on the analysis of the complex coupling relationship among the electric energy subsystem, heat energy subsystem, cooling energy subsystem, and gas energy subsystem, it can be discovered that there exists strong coupling relationships among various sub-energy systems. Therefore, it is necessary to forecast the combined electricity-heat-cooling-gas loads instead of predicting a single energy load. Additionally, the accurate prediction of combined electricity-heat-cooling-gas loads on the demand side can provide important reference and decision-making basis for multiple energy planning and optimal scheduling, demand side management, as well as safe and stable operation of multi-energy systems in the IES of the industrial park. Hence, establishing an uncertain interval forecasting model with high accuracy for combined electricity-heat-cooling-gas loads in the IES is of great significance. This paper combines the MTL method, the Bootstrap method, the ISSA and the MKELM method to establish the uncertain interval prediction model of combined electricity-heat-cooling-gas loads on the demand side of the IES in the industrial park. The ISSA proposed in this paper introduces the dynamic inertia weight and chaotic local searching mechanism into the basic SSA optimization algorithm, so as to improve the searching speed of the basic SSA and avoid falling into local optimum easily. The MKELM model is also established in this paper by combining the RBF kernel function and the Poly kernel function, so as to integrate the superior learning ability and generalization ability of the two kernel functions. Based on the established uncertain interval prediction model, taking a typical industrial park in a northern province of China as an example, according to the natural, technological, and social conditions in the region, the model test and verification are carried out for the load data of electricity, heat, cooling, and gas in the IES of the industrial park.
When forecasting the electricity-heat-cooling-gas loads in the IES of the industrial park employing the established uncertain interval prediction model, four factors including weather, calendar information, social-economic factors, and historical load are selected as the input variables of the prediction model. The results of the empirical analysis and comparative analysis are as follows: (1) The prediction results of the electricity-heat-cooling-gas loads of the IES on the workday and holiday forecasted by the established uncertain interval prediction model integrating the MTL method, the Bootstrap method, the ISSA, and the MKELM method are in good accordance with the actual loads. Additionally, the prediction results on the workday are better than those on the holiday, although the fluctuation of electric, heat, cooling, and gas loads are varied on the workday and holiday. (2) Through comparing the prediction accuracy of the Bootstrap-ISSA-MKELM based on the MTL method and the Bootstrap-ISSA-MKELM based on the STL method, it can be found that the Bootstrap-ISSA-MKELM based on the MTL method has superior performance than that based on the STL method. This is primarily because the MTL mechanism can comprehensively consider the complicated coupling relations among electric, heat, cooling, and gas loads in the IES, and make full use of the information of multi-energy loads. When one load is predicted, the information contained other related loads will be utilized, so that the prediction accuracy of electric, heat, cooling, and gas loads has been greatly improved. Additionally, the MTL mechanism can effectively shorten the calculation time. This is mainly because that the uncertain interval prediction results of the four kinds of loads can be obtained by one time via the model based on the MTL method as the multi-energy loads data can be input into the model at one time. However, the STL mechanism needs to train each energy load one by one, which makes the cumulative training time becomes longer.