Optimal Economic Dispatch for an Industrial Park With Consideration of an Elastic Energy Cloud Model With Integrated Demand Response Uncertainty

Under China’s vigorous development of integrated energy services, the Integrated Energy Service Agency (IESA) is responsible for purchasing energy from external markets and selling energy to multi-energy users (MEUs). Currently, an increase in the various forms of energy in industrial parks has caused great uncertainty for MEUs participating in an integrated demand response (IDR) but has also provided an opportunity for industrial parks to optimize energy conservation. Therefore, determining how to build an elastic energy cloud model with IDR uncertainty and add the uncertainty and randomness of the cloud model to the optimal scheduling of an industrial park integrated energy system is a key problem. In this paper, an optimal economic dispatch model of an industrial park is proposed and considers the uncertain elastic energy of IDR. In this model, the IESA is responsible for the reasonable scheduling of equipment for optimal operation, the establishment of integrated energy retail prices for MEUs, and the goal of maximizing the net income of the IESA. First, the functional relationship among the self-elastic coefficient, retail energy prices, and IDR variation is considered. A cloud model of the self-elastic coefficient is constructed to indirectly represent the multiple uncertainties of the elastic energy in the industrial park. Second, this paper compares and analyzes the economic benefits and IDR potential of the industrial park by considering only single power users in different intervals and the selection of cloud drop elements of MEUs in all intervals. Finally, a new scene random sampling method based on interval contributions (SRS-IC) is employed to solve the optimization model, and a typical example is used to demonstrate that the model and method can guarantee the overall economy of the industrial park, improve the computational efficiency, and explore the IDR potential of MEUs.

B. PARAMETERS f Net Optimization objective function ( ) T Length of the dispatch period (h) R Total Net revenue from IESA sales ( ) C Om Total operating cost of integrated energy system ( ) C Buy Total transmission power cost purchased by the industrial park from external power distribution, gas, cold and heat networks ( ) R Rev,sal,total Total revenue from IESA sales ( ) R Rev,cot,total Total cost and benefit of P2G and GSHP ( ) R E_DL/GL/HL/CL Total retail electricity/gas/heat/cold energy revenues of the IESA ( ) R E_P2G Cost and benefit generated by power-to-natural gas equipment ( ) R E_GSHP Cost and benefit of refrigeration and heat generation of GSHP ( ) C PV/WT/MT/GB/ES_Om Operating cost of PV/WT/MT/ GB/ES ( ) P t Users of shiftable gas/cold energy participating in the output power increase/decrease of the IDR program (kW)

I. INTRODUCTION
Under the background of China's energy market reform and the gradual development of integrated energy services on the demand side, the Integrated Energy Service Agency (IESA) can develop more energy sales packages for multi-energy users (MEUs) [1], [2]. With the increase in multi-type energy varieties containing MEUs, when an industrial park meets the requirements of MEUs participating in an integrated demand response (IDR), the optimized energy-saving effect of the demand-side integrated energy system (IES) is not significant, especially for multi-type energy varieties in industrial parks that have multiple uncertainty characteristics [3], [4]. Currently, the forms of energy use of integrated energy terminal MEUs are divided into rigid energy and elastic energy. Rigid energy is not subject to changes in some incentives, while elastic energy is the part of energy or time that can be adjusted by the stimulation of some factors [5], [6]. Therefore, this paper establishes an elastic energy cloud model that comprehensively considers the uncertainty of MEU participation in IDR and then designs a set of economic scheduling strategies that consider the uncertainty and randomness of the elastic energy cloud model. This strategy not only guarantees overall economic benefits for the industrial park but also maximizes the IDR potential of MEUs. Many studies have focused on uncertainty models, strategies, solutions and demand response (DR) mechanisms. Zhao et al. [7] proposed a multistage robust optimization (RO) model for unit commitment that considered the wind power output and DR uncertainties, and the DR was modeled as an uncertain price-elastic demand curve. Wei et al. [8] presented a two-stage two-level decision model for a smart grid retailer that considered the DR and market price uncertainty. Yi et al. [9] established a robust interval model that considered the uncertainties of renewable energy power generation output and DR. Then, a robust counterpart transformation method combined with a multi-objective genetic algorithm was applied to solve the model. Yang et al. [10] proposed a two-layer nested model for a hydro-wind complementary system that accurately considered the characteristics of intermittency, the fluctuation characteristics of wind power, and forecasts. Li and Jayaweera [11] proposed a series of linear prediction models for load management that considered different assumptions on the stationarity of customer load data and then further designed a DR scheduling scheme based on utility cost minimization with different customer clustering sizes. Ju et al. [12] proposed a multi-objective flexible risk aversion model for power-to-gas (P2G)-based virtual power plant (VPP) operation that considered the uncertainty risk described by a conditional value-at-risk method and robust theory. Pan et al. [13] proposed a synchronously decentralized adaptive robust programming model for an IES that considered the uncertainties of photovoltaic (PV) energy and flexible IDR loads. The above studies have contributed significantly to the uncertainty modeling of demand-side IESs based on their respective approaches.
Some uncertain factors have also been considered in the operation of demand-side IESs, such as economic dispatch strategies and solution methods. For example, Majidi and Zare [14] proposed a model for both the price and power demand uncertainty of renewable units, and a stochastic programming method combined with an interval optimization method was employed to solve the model. Chen et al. [15] presented a two-stage distributed robust hydrothermal-wind economic dispatch strategy, which had the advantages of traditional stochastic optimization (SO) and adjustable RO methods, and a delayed constraint generation algorithm was applied to solve the strategy. Zheng et al. [16] proposed a comprehensive system model for a coupled natural gas and electricity network that considered the impact of the power output uncertainty of wind farms. Then, an interval algorithm was employed to solve an uncertainty analysis of the coupled system and compared with a Monte Carlo (MC) method. Bai et al. [17] proposed an interval optimization-based coordinated operation strategy for a gas-electric IES that considered DR and wind power uncertainty. Bai et al. [18] adopted a SO model to deal with the uncertainty generated by a large number of wind power forecast scenarios. Then, a scenario reduction algorithm was applied to solve the model. Furthermore, Zhang et al. [9] proposed a multi-objective robust scheduling method that considers the uncertainty of the predicted values of input variables, such as the load demand and renewable energy power output. Then, a simultaneous backward reduction (SBR) method was applied to reduce the scenario that meets the requirements. Liu et al. [19] proposed a stochastic RO framework to solve the distributionally robust unit commitment (UC) problem, and the wind power moment information of uncertainty parameters could be obtained from the statistics of historical data.
Moreover, some researchers have studied the factors influencing the DR uncertainty in the operation of demand-side IESs, such as the DR mechanism. For example, Niu et al. [20] established an optimization model comprising a VPP model with an uncertain DR. Sun et al. [21] described the response mechanisms of different types of demand-side resources in detail. Then, they proposed a day-ahead scheduling model based on fuzzy chance constrained programming, and the user response behavior under time-of-use (TOU) was modeled with uncertainty. Wu et al. [22] proposed a multi-time period flexible uncertain DR optimization model that considered the uncertainty impacts of the DR on the optimal operation. Ju et al. [10] proposed a bi-level stochastic scheduling model for a VPP that considered the uncertainties of the wind power plant and PV output, price-based demand response and incentive-based demand response. Then, a robust SO theory was introduced to solve the model. Wang et al. [23] proposed an incentive-based DR unit commitment model that considered different types of users according to their DR characteristics. Wang et al. [24] established a fuzzy load response model under the TOU electricity price that considered the uncertainty of the price DR, wind power output and system load. Gu et al. [25] proposed a bi-level optimal low-carbon economic dispatch model for an industrial park that considered multi-energy price incentives and the IDR of MEUs. In [26], the cloud model theory was applied to the feature mining of wind power prediction training samples. In [27], a cloud model was applied to address large group decision-making problems based on the cloud model in a linguistic environment, and the model was proven a powerful tool for expressing the fuzziness and randomness of the linguistic information.
Previous studies have mainly focused on the uncertainty analysis of the wind power output and electric load DR and have less comprehensively involved the analysis of the IESA's IDR uncertainty optimization scheduling strategy for electricity, heat, cold and gas MEUs. First, there are three methods to deal with the objective function: RO, distributed RO and random programming. Second, almost all fuzzy day-ahead optimal scheduling models with DR uncertainty are solved by transforming the expectation of the uncertainty model and the fuzzy opportunity constraint into a clear equivalence model according to uncertain programming theory. In addition, some previous scholars have affected the operation of industrial parks by adjusting the flexible load in the IDR scheme corresponding to the multi-energy pricing strategy. It can be found that with an increase in the various types of energy in an industrial park, it is difficult to quantify an existing mathematical model of uncertainty, the existing uncertainty model has difficulty measuring randomness and fuzziness at the same time, and it is difficult to describe the uncertainty characteristics of elastic energy. Therefore, it is necessary to use a cloud model to represent the multiple uncertainty characteristics of elastic energy in an industrial park, to tap the potential energy use of MEUs, to improve the economic and energy conservation benefits of the industrial park, and to optimize an economic scheduling scheme of uncertainty and randomness for MEUs to participate in IDR.
Motivated by the abovementioned research gap, this paper proposes an uncertainty optimal economic dispatch strategy for an industrial park with consideration of the elastic energy cloud model. This optimization model includes the IESA setting retail energy prices based on external market prices, taking into account the IESA's retail energy prices and net/total revenue constraints, with the optimal net income of the IESA as the overall goal. In an actual IDR process, the elastic coefficient matrix can represent the functional relationship between the rate of change in the retail energy price and the rate of change in IDR. Considering the uncertainty characteristics of MEUs, such as diversification, the cloud model can be used to represent the multiple uncertainties of elastic energy in the industrial park, and then an improved random sampling method is used to cut the scene to obtain the overall optimal operation results of the industrial park. The contributions of this paper are as follows: 1) An uncertain optimal economic dispatch strategy for an industrial park considering an elastic energy cloud model is introduced, and this strategy considers the elastic energy varieties for MEUs, including electricity, heat, cold, and gas, and has multiple uncertainty characteristics, such as the time transfer and spatial distribution. The strategy can also help the IESA develop a reasonable integrated energy retail price optimization strategy and improve the IDR potential of MEUs. 2) An elastic energy cloud model with three eigenvalue parameters (expection, entropy, and hyper entropy) is established, and this model can represent the uncertainty of the self-elasticity coefficient, the IESA's integrated energy retail price, and MEU participation in IDR. 3) A new concept of interval contribution in the elastic energy cloud model is proposed, which can help analyze the contribution of the selection of different interval elements in the model and the contribution to the overall economic benefits and IDR potential of the industrial park. 4) A new scene random sampling method based on interval contribution (SRS-IC) is employed, which adopts the SBR method to improve an MC random sampling method. This method can effectively improve the system operation efficiency and reduce the error of uncertain optimization results.
The rest of this paper is organized as follows. Part 2 introduces cloud model theory and the elastic energy cloud model. Part 3 introduces the structure and optimization model of the demand-side IES. Part 4 presents the optimization solution strategy for the proposed model in detail. In part 5, case studies and simulation results are presented. Part 6 presents the summary and conclusions. In part 7, future work is presented.

II. ESTABLISHMENT OF THE ELASTIC ENERGY CLOUD MODEL A. CLOUD MODEL THEORY
Cloud model theory can realize the uncertainty conversion between qualitative concepts and their numerical representations [26]. The cloud model adopted in this paper is a normal cloud with the most general adaptability. As shown in figure 1, the abscissa and ordinate represent expection (Ex) and its determinacy (µ) respectively. The cloud model can adopt three eigenvalue parameters Ex, entropy (En), and hyper entropy (He), to reflect the distribution characteristics of the concept as a whole [26]. Among these parameters, Ex is the expectation of cloud droplets belonging to a qualitative concept in the universe; En is the degree of uncertainty of a qualitative concept and is measured by both the fuzziness and randomness of the qualitative concept; and He is the degree of uncertainty of the entropy En (refer to definitions 1 and 2 in Appendix). The cloud model can use a membership function to describe the fuzziness of the uncertainty concept, the randomness of the membership function and cloud droplets. The cloud model perfectly combines fuzziness and randomness. At present, the core algorithm of the cloud model is divided into a forward Gaussian cloud generator (FGCG) and a backward Gaussian cloud generator (BGCG) [29], [31]. The specific algorithm flow is described in (a) and (b) as follows:

1) FGCG ALGORITHM GENERATE CLOUD DROPLETS
Input: Three eigenvalue parameters, Ex, En, and He, and the number of cloud droplets N .
(1) Calculate the average value of x ϕ : (2) Calculate the first-order absolute center moment of the sample: (3) Calculate the second-order center moment of the sample: (4) Calculate the entropy: (5) Solve the hyper entropy: Inspired by the above cloud model idea, this paper first considers the difficulty of obtaining historical data samples of the actual energy demand of MEUs and the integrated energy retail price of the IESA, and then three historical eigenvalue parameters of the cloud model are considered to describe the uncertain distribution law of the self-elastic coefficient. Among these parameters, Ex of the self-elastic coefficient is the most representative point of the qualitative concept. En is a random measure of qualitative concepts. On one hand, it reflects the degree of dispersion of qualitative concepts, and on the other hand, it determines the cloud droplet determination acceptable to concepts in the domain space. He reflects the acceptance of the qualitative concept. The smaller is the value of He, the more concentrated are the cloud droplets, while the larger the value of He is, the more discrete the cloud droplets are. At present, the FGCG and BGCG algorithms in the cloud model can effectively convert between the qualitative concepts of fuzzy variables and clouds. Additionally, a large number of cloud titration values in the cloud model will also cause difficulties in data processing in the uncertain optimal scheduling process of industrial parks.

B. ELASTIC ENERGY CLOUD MODEL
In the actual optimal dispatching process of an industrial park, the actual energy demand of MEUs is often uncertain, and the integrated energy retail price of the IESA is bound to be uncertain. Therefore, the three historical numerical eigenvalues of the cloud model can represent the uncertainty distribution of the self-elastic coefficient. These values can also represent the multiple uncertainty characteristics of changes in elastic energy and changes in the retail energy price of the IESA in the industrial park, which are described as follows: whereε t dl /ε t hl /ε t cl /ε t gl represent the self-elastic coefficients of electricity/heat/cold/gas MEUs at time t represented by the cloud model, respectively;˜ c t dl,prc /˜ c t hl,prc /˜ c t cl,prc /˜ c t gl,prc represent the changes in the retail prices of electricity/heat/cold/gas of the IESA after the implementation of IDR at time t by using the cloud model, respectively; and˜ P t dl,elc,aft /˜ Q t hl,elc,aft /˜ Q t cl,elc,aft /˜ Q t gl,elc,aft represent the actual changes in elastic electrical/heat/cold/gas loads after the implementation of IDR at time t by the cloud model, respectively.
According to formulas (A-3) ∼ (A-6) in Appendix, formula (9) can be simplified into the following form (10): Finally, the same formulas (A-7) ∼ (A-17) in Appendix are derived to obtain a piecewise linear function of the total return of the IESA and the corresponding constraint conditions under uncertain optimal scheduling.
In figure 2, the abscissa represents the expected value of cloud droplets, and the ordinate represents the degree of certainty. For example, figure 2 (a) contains four cloud models with 3 specific numerical eigenvalues. Cloud models 1, 2, 3, and 4 adopt (Ex 1 , En 1 , He 1 ), (Ex 2 , En 2 , He 2 ), (Ex 3 , En 3 , He 3 ), and (Ex 4 , En 4 , He 4 ), respectively, to represent the overall characteristics of their concepts. When the point in a cloud model that best represents its qualitative concept is closer to Ex, the cloud droplet density is higher. In the opposite case, the droplets become sparser. In addition, when the Ex value is different, the cloud image of the cloud model will be shifted to the left and right along the abscissa. When the En value is different, the larger is the value of En, the wider is the value range of a cloud droplet. In contrast, the smaller is the value of En, the smaller is the value range of a cloud droplet. When the He value is different, the larger the value of He, the more discrete are the cloud droplets. For example, compared with the other cloud models, cloud model 3 more easily demonstrates this trend.
When there is only a single form of energy use in an industrial, the physical meaning of the traditional uncertain optimization (such as RO, distributed RO and random programming) model can also be expressed. However, when VOLUME 9, 2021 there are the four energy forms of electricity, heat, cold and gas in an industrial park, it is necessary to construct a four-dimensional elastic energy cloud model to characterize the multiple uncertainty characteristics of multi-type energy varieties containing MEU members. Therefore, from the perspective of model construction, the construction of uncertainty multidimensional cloud model is essentially different from the traditional uncertainty optimization modeling. In addition, according to the basic theory of multidimensional cloud models [26], an n-dimensional normal cloud model can be combined mathematically with n one-dimensional normal clouds. This approach can not only can guarantee the mutual independence of cloud droplet arrays in the cloud model but can also describe the complex qualitative concept of an industrial park containing n types of random energy combinations.
(3) Calculate the determinacy: is a cloud droplet in the domain. (5) Repeat steps 1 to 4 until N cloud droplets are formed. According to the above solution steps, the calculation of the four-dimensional normal cloud model is complex, and it is impossible to effectively describe the uncertainty characteristics of energy user types in the four dimensions of electricity, heat, cold and gas in a spatial coordinate system. Therefore, the dimensionality reduction of the multi-dimensional cloud model can effectively simplify the uncertainty model in this paper. For example, figure 2 (b) shows one cloud model of a spatial distribution. Combined with analysis of the results of the calculation examples, the time (T ), expection (Ex), and determinacy (µ) are used to form four three-dimensional coordinate systems. Figure 3 shows an optimal scheduling structure with information exchange between the IESA and MEUs of an industrial park. Figure 4 shows a busbar terminal structure [25]. The overall structure of the IES includes the IESA and MEUs. The IESA buys energy from external energy markets and sets different retail energy prices for MEUs. The IESA is responsible for the unified and optimized management of the operation of comprehensive energy supply, energy storage and energy conversion equipment in the industrial park. The integrated energy supply equipment includes distributed energy such as WTs, PVs, MTs, RECs and GBs. The energy storage equipment includes ES_BT, ES_HC, ES_CD, and ES_NG; the integrated energy supply and energy conversion devices include P2G; and the integrated energy hub conversion devices comprise an EC, EH, HX, AC and GSHP. The form of energy use of integrated energy terminal MEUs are divided into rigid energy (basic energy) and elastic energy. Rigid energy is not subject to some excitation changes and mainly includes basic (rigid) electrical/gas/cold/heat loads (BDL/BGL/BCL/BHL). Elastic energy is the part that can be used to adjust the energy or time with some factors and mainly includes elastic electrical/gas/cold/heat loads (EDL/EGL/ ECL/EHL). In figures 3 and 4, the symbols used for each equipment type in the industrial park are defined according to the parameter nomenclature.

A. OPTIMIZATION MODEL 1) OBJECTIVE FUNCTION
The optimization in this paper aims to maximize the net income of the IESA, that is, to maximize the difference between the total income of the IESA and the total operating cost of the system. The total retail energy revenue of the  IESA, f Net , mainly includes the sales of electricity, heat, cold, and natural gas and mainly comes from the IDR benefits of rigid energy and elastic energy users in the industrial park. In addition, the production costs and benefits mainly come from P2G and GSHP energy conversion. The total operating cost of the system mainly includes the operating cost of the integrated power supply equipment (WTs, PVs, MTs, GBs, and ES) and the cost of power transmission purchased by the industrial park from the distribution/heat/cold/gas pipeline network. Therefore, the optimization objective function f Net is shown as follows: where R Total represents the total revenue of the IESA; C Om represents the total operating cost of each price of energy supply equipment and energy storage equipment in the industrial park; and C Buy represents the total power transmission cost of each industrial park purchased from external distribution, heat, cold, and natural gas networks.
R Total = R Rev,sal,total + R Rev,cot,total (13) Here, the first term in (13) indicates the total retail sales revenue of the IESA, and the second term in (13) indicates the costs and benefits of the power-to-gas and ground-source heat pump units in the industrial park. (14) Here, the first through fourth terms in (14) indicate the total revenue of the IESA from retail electricity, heat, cold, and gas energy, respectively.
Here, the first term in (15) represents the production costs and benefits of using the P2G device to convert power to gas energy. The second term in (15) represents the production costs and benefits of using the GSHP for heating and cooling generation.
Here, c t dl,prc,aft /c t hl,prc,aft /c t cl,prc,aft /c t gl,prc,aft indicate the real-time retail cost prices of electricity/heat/cold/natural gas after the IESA implements IDR at time t ( /kWh), respectively. P t dl,ld,aft /Q t hl,ld,aft /Q t cl,ld,aft /Q t gl,ld,aft in (16) represent the real-time electricity/heat/cold/gas load after the IESA implements IDR at time t (kW), respectively.
Here, c t p2g,cos t indicates the benefit-cost coefficient of the P2G device converting electric power to natural gas at time t ( /kWh). Q t P2G_EX in (17) means the power from gas converted through P2G equipment at time t (kW). c t co,cos t /c t he,cos t in (17) indicate the revenue and cost coefficients of the GSHP unit generated by refrigeration and thermal energy at time t ( /kWh). Q t HP_C /Q t HP_H in (17) mean the cooling and heating power of the GSHP unit at time t (kW).
Here, the first through fifth terms in (18) indicate the operating costs of the WTs, PVs, MTs, GBs, and ES, respectively. The details are as follows: Here, c wt , c pv , c mt , and c gb are the operating cost coefficients of the WTs, PVs, MTs, and GBs ( /kWh), respectively. c bt , c hc , c cd , and c ng in (19) are the operating cost coefficients of the integrated energy storage devices (batteries, heat, cold, and natural gas tanks) ( /kWh), respectively. P t WT , P t PV , P i,t MT , and Q j,t GB in (19) mean the output power of the WTs, PVs, MTs, and GBs (kW), respectively. P t (19) mean the charging or discharging power of the batteries, heat, cold, and natural gas storage tanks in the industrial park (kW), respectively.
The IESA considered in this paper is not only used for retail energy but is also used for the electricity, heat, cold, and natural gas energy purchased from external markets. The details are as follows: Here, c t dl,buy , c t hl,buy , c t cl,buy , and c t gl,buy represent the power cost coefficients of purchases by the industrial park from power distribution, heat, cold, and natural gas networks ( /kWh), respectively. P t G,Buy , Q t H,Buy , Q t C,Buy , and Q t NG,Buy in (20) represent the power purchased by the industrial park from power distribution, heat, cold, and natural gas networks (kW), respectively.

2) ESTABLISHMENT OF THE INTEGRATED ENERGY RETAIL PRICE AND IDR MODEL
The balanced constraints of this section comprise the following four parts: (21) ∼ (24). In the industrial park, the IESA is responsible for setting the real-time retail energy prices of electricity, heat, cold and gas, while the MEUs change their energy demand in real time according to the price signal. Considering the price-type IDR model, the elastic coefficient can describe the functional relationship between the change rate of a retail energy price and the change rate of the IDR [24]. To facilitate data processing, only the self-elastic coefficient is considered in this model, which is described as follows: where λ dl /λ hl /λ cl /λ gl are the proportions of initial elastic energy in the electricity/heat/cold/gas load predicted before IDR, respectively. P t (21)  in (21) are the changes in the retail prices of electricity/heat/cold/gas after IDR at time t ( /kWh), respectively. In combination with formula (21), in (16), the real-time retail energy cost prices after IDR and the above retail energy price changes should satisfy the following constraints: In addition, the real-time electricity, heat, cold, and natural gas elastic load variation and the reference values of the initial elastic load variation after the implementation of IDR should satisfy the following constraint conditions:   (23) mean the actual initial elastic load predictions before IDR at time t (kW).
Combined with formulas (21) and (23), in (16), the total loads of electricity, heat, cold, and gas in the industrial park in real time after the implementation of DR should meet the following power balance constraints: By combining (22) and (24), formula (16) can be converted into the following form: where formula (25) contains the product of continuous variables. To facilitate the programming calculation, the above formula can introduce the auxiliary variables, κ t dl,prc,,ref , Then, based on the further derivation of (A-1) ∼ (A-7) in Appendix, the quadratic function expression of the total revenue of MEUs after the implementation of IDR can be finally VOLUME 9, 2021 obtained as follows: where the MEUs have the product constraint of continuous variables in the quadratic function expression of the total return after IDR implementation. Therefore, a piecewise linear function processing method [28] can be adopted, and the Big M method [25] can be used to transform the nonlinear model into a linearized model for solving. Please refer to (A-8) ∼ (A-17) in Appendix for the detailed formula derivation.

3) BALANCED OPTIMIZATION CONSTRAINTS
The balanced constraints in this part comprise the following five formulas: (27) ∼ (31).
(1) The energy flow balances of electricity, heat, cold, natural gas, and auxiliary flue gas buses are as follows: where the first constraint in (27) is the electrical power balance, the second is the heat power balance, the third is the cold power balance, the fourth is the natural gas power balance, and the fifth is the auxiliary flue gas power balance. Please refer to figure 4 for specific symbol details.
(2) The constraints of the integrated energy storage devices are as follows: where the first to fourth constraints in (28) are the relationships between changes in the storage energy of batteries, heat, cold, and natural gas tanks and changes in the energy storage charge or discharge power efficiency, power, and selfdischarge efficiency and duration, respectively.
(3) The energy balance constraints of the power supply devices are as follows: where the first constraint in (29) is the thermoelectric energy balance of MTs and the thermoelectric ratio of MTs. The second constraint in (29) is the gas-to-electric energy balance of MTs. The third through fifth constraints in (29) are the electric-to-gas balance of P2G, gas-to-thermal energy balance of GBs, and heat-to-thermal energy balance of RECs, respectively. (4) The energy conversion constraints of the devices are as follows: where the first to fifth constraints in (30) are the energy conversions of the EC, EH, AC, REC, HX, and GSHP, respectively.
The equality constraints of shiftable electrical, heat, cold, and gas load variation are as follows: Formula (24) describes the load power balance constraint of MEUs considering that the elastic load curve should follow its energy continuity, timing and overall properties after IDR implementation. Therefore, the specific constraints are as follows: where the first and second constraints in the first row and the first and second constraints in the second row represent the increases or decreases in the total shiftable electrical, heat, cold, and gas loads and should be equal to the sum of the load changes during period T .
(1) The constraints of each type of power supply device are as follows: where the first to fifth constraints in (32) are the output power limitations of the MTs, WTs, PVs, RECs, and GBs, respectively. P i,t MT,max /P i,t MT,min , P t WT,max /P t WT,min , P t PV,max /P t PV,min , Q t REC,max /Q t REC,min , and Q t GB,max /Q t GB,min are the allowable maximum and minimum output power of MTs, WTs, PVs, RECs, and GBs at time t (kW), respectively. R i MT,up /R i MT,dn are the limitations of the climbing force constraint coefficient of MTs (kW/min), respectively.
(2) The constraints of the integrated energy conversion devices are as follows: where Q t GSHP,max /Q t GSHP,min , Q t P2G,max /Q t P2G,min , Q t EC_EX,max / Q t EC_EX,min , Q t EH_EX,max /Q t EH_EX,max , Q t EH_EX,max /Q t EH_EX,max and Q t HX_EX,max /Q t HX_EX,min are the allowable maximum/minimum output power limitations of the GSHP, P2G, EC, EH, AC, and HX at time t (kW), respectively.
(3) The power transmission constraints of the electrical, heat, cold, and gas networks are as follows: where P t G,Buy,max /Q t H,Buy,max /Q t C,Buy,max /Q t NG,Buy,max are the maximum electrical/heat/cold/natural gas power that the IESA buys from the external distribution/heat/cold/gas networks at time t (kW), respectively.
(4) The capacity constraints of the integrated energy storage devices are as follows: where the first and second constraints in (35) represent the charging and discharging power of the battery storage device, respectively. The third constraint in (35) represents the minimum and maximum capacity constraints of the device. The fourth and fifth constraints in (35) represent the charge and discharge power of the heat/cold tank storage device, respectively. The sixth constraint in (35) indicates the minimum and maximum capacity constraints of the device. The seventh and eighth constraints in (35) represent the charging/discharging power of the air storage tank device, respectively. The ninth constraint in (35) means that the minimum and maximum capacity constraints of the device. ξ t bt/ng/cd/hc,c /ξ t bt/ng/cd/hc,d are binary 0-1 variables. (5) The inequality constraints of shiftable electrical, heat, cold, and gas load variation are as follows: (6) The constraints of the changes in the energy retail price and elastic load are as follows: where the first through fourth constraints in (37) are the upper and lower limits of the variation in the retail prices of electricity, heat, cold and gas, respectively. The fifth through eighth constraints in (37) are the upper and lower limits of the total elastic load after the electricity, heat, cold, and gas loads participate in IDR, respectively. (7) The constraints of the rates of change in retail energy prices are as follows: To guarantee the rationality of the retail energy price of the IESA, the constraint relationship between the rate of change in the retail energy price and the rate of change in the IDR should be comprehensively considered in combination with the actual IDR income model proposed in this paper. Therefore, the rates of change in the IESA's retail energy prices should obey the following constraints: where the first through fourth constraints in (38) are the upper and lower limits of the real-time rates of change in the retail energy prices set by the IESA, respectively. (8) The constraints of the IESA's net/total energy sales income are as follows: Considering the multiple uncertainty characteristics of elastic energy represented by the cloud model in the industrial park, it is easy for the cloud droplet distribution in the cloud model to have a high or low degree of dispersion. As a result, the IESA's energy sales revenue under the uncertain optimal scheduling scenario is unreasonable. Therefore, in the uncertain optimization scenario, the IESA needs to set a reasonable constraint range of net energy sales/total income based on the deterministic optimization results as follows: where the first through fourth constraints in (39) are the upper and lower limits of the IESA's net income from the sales of electricity, heat, cold, and gas under uncertain optimal dispatching, respectively. The fifth constraint in (39) is the upper and lower limits of the IESA's total energy sales revenue under uncertain optimal scheduling and is proposed on the basis of the IESA's net income from energy sales.

IV. OPTIMIZATION SOLUTION STRATEGY A. SIMULTANEOUS BACKWARD REDUCTION (SBR) METHOD
To simplify the random sampling and combination scenario of cloud droplet data samples in the elastic energy cloud model of the industrial park, the computational efficiency of the industrial park in the uncertain optimal dispatching process is improved. Based on the interval contribution of the electric load cloud model, this paper proposes a new SRS-IC method and adopts the SBR method [9] to improve a Monte Carlo random sampling (MC-RS) method, and this approach not only avoids infeasibility of the optimization results but also reflects the economics. This method first needs to reduce the sample random combination scene size of the self-elastic coefficient in the cloud model and only keeps a few representative sample scene sets to replace the original scene sets with a large amount of data to ensure the minimum probability distance between the two scene sets.
The specific steps of the SBR method are as follows: Step 1: Suppose that the set of initial scenes is D and that the number of initial scenes is Ns, £ s (s = 1, 2, . . . , N s ). The probability of each initial scenario is P r (s). The initial scene set is ω, and the set of cut scenes is β. In the actual calculation, ω and P r (s) are initially set to empty collections.
Step 2: During the lth iteration, calculate the Euclidean distance (ED: min D l,s (£ l , £ s )) between the probability of £ l of each scene and the nearest scene of £ s , and then select the probability P r (s ) of scenario £ s occurring, with l ∈ R, s ∈ R and l = s . Then, according to (39), the scene £ s with the shortest distance from scene £ l is calculated as follows: Step 3: According to formula (40), determine the minimum probability distance D min between the scene to be deleted and the reserved scene as follows: Step 4: After each scene is cut, the number of remaining scenes needs to be modified, and the probability of the scene being cut is added to the nearest scene, R = R − { }, J = J + { }, and P r (s ) = P r (s ) + P r ( ), to ensure that the sum of the probabilities of the newly generated scene is 1.
Step 5: Repeat (Step 2) ∼ (Step 4) until the number of remaining scenes reaches the desired number of scenes.

B. FLOW CHART OF THE UNCERTAINTY OPTIMIZATION ECONOMIC DISPATCHING STRATEGY
The uncertainty optimal economic dispatching strategy for an industrial park proposed in this paper considering the elastic energy cloud model is shown in figure 5.

V. CASE STUDIES A. CASE DESCRIPTION
This section establishes an optimal economic dispatch strategy for an industrial park considering the IDR uncertainty elastic energy cloud model. The model and algorithm in this paper are written in MATLAB 2018a and run on a computer with an Intel (R) Core (TM) i7-9700 CPU with 3.00 GHz main frequency and 32 GB memory. The correctness and validity of the proposed method are verified by the following case scenarios. In this paper, a typical case is used for analysis and verification-a discrete manufacturing industrial park. It is assumed that the industrial park has four types of MEUselectricity, heat, cold and gas energy-and that the total power of electricity, heat, cold and gas energy in the example is approximately equal to 4095, 4470, 4335 and 3024 kW, respectively [25]. Figure 6 shows the external and initial energy reference prices of the retail electricity, heat, cold and natural gas unit costs for MEUs established by the IESA. It is assumed that the initial retail electricity, heat, cold and gas prices are 0.55, 0.5, 0.5 and 0.75 /kWh, respectively. The typical load prediction curve of MEUs on a certain day is shown in figure 7. Please refer to Section 2 for descriptions of the main equipment types in the industrial park and refer to figures 3 and 4 for equipment symbols and names. In this case, the daily electric load curve of the discrete manufacturing industrial park shows a ''double peak, one low peak and two concaves'' pattern, the daily heat load curve shows a ''single peak'' pattern, the daily cold load curve shows a ''single peak, two low peaks'' pattern, and the daily gas load curve shows a ''double peak, two concaves'' pattern. The relevant equipment parameters of each industrial park are shown in Table 1. The unit of active power in Table 1 is kW, and the unit of the cost coefficient is /kWh.

1) SCENARIO DESCRIPTION
In this paper, the following 7 different scenarios are set for analysis of the uncertainty optimal economic dispatch strategy of the discrete manufacturing industrial park. Scenario 1 is a certain optimal scheduling strategy, while scenarios 2∼7 are uncertain optimal scheduling strategies, and the selection of all self-elastic coefficient cloud drop elements in VOLUME 9, 2021 FIGURE 6. Flow chart of the uncertainty optimization economic dispatching strategy for an industrial park considering the elastic energy cloud model. the cloud model is considered. Scenario 2 considers only the elastic load cloud model, while scenarios 3∼7 comprehensively consider the elastic energy cloud model of electricity, heat, cold and gas users.
In scenarios 3∼5, the traditional MC-RS method is used for the random sampling of cloud droplets with the self-elastic coefficient of MEUs. Scenarios 6 and 7 based on scenario 5 adopt the improved SBR method to improve the MC-RS method and set the desired scenario reduction numbers. The corresponding scenario settings are shown in Table 2, and the detailed descriptions of scenario construction are as follows: Description (a): The self-elastic coefficients of electricity, heat, cold and gas energy users adopt constant values, and ε t dl /ε t hl /ε t cl /ε t gl = −0.

2) UNCERTAINTY OPTIMIZATION RESULTS ANALYSIS
(1) Uncertainty optimization strategy analysis with the electric load cloud model (a) Schematic diagram of the electric load cloud model based on interval contributions According to the ''3En rule'', the self-elastic coefficient cloud droplet elements (hereinafter referred to as ''elements'') in the cloud model mainly fall within the interval [Ex-3En, Ex + 3En], and the contribution of those elements outside that interval to the qualitative concept can be ignored [29]. In Scenario 2, only the self-elastic coefficient of power users is represented by the digital feature of the cloud. Based on the above ''3En rule'', in view of the randomness and fuzziness of the selection of elements in different intervals in the electric load cloud model, the contribution of the selection of elements in 10 different intervals to the overall economic benefit of the industrial park and the potential of MEUs to participate in IDR were analyzed. The specific schematic diagram is shown in figure 8.
In figure 8, all elements in Interval 1 are referred to as ''backbone elements'', which contribute 50% to the qualitative concept as defined in Appendix [29]. All elements in Interval 2 are called ''skeleton elements less than 0'', and all elements in Interval 3 are called ''skeleton elements greater than 0.'' The elements in Interval 4 are called ''all basic elements'', which contribute 68.26% to the qualitative concept [29]. All elements in Interval 5 are called ''base elements less than 0'', and all elements in Interval 6 are called base elements greater than 0. All elements in Interval 7 are called ''peripheral elements less than 0'', and all elements in Interval 8 are called ''peripheral elements greater than 0''. All elements in the above two parts contribute 27.18% to the qualitative concept [29]. All elements in Interval 9 are called ''weak peripheral elements less than 0'', and all elements in Interval 10 are called ''weak peripheral elements greater than 0''. All elements in the above two parts contribute 4.3% to the qualitative concept.
(b) Uncertainty analysis of revenues, costs, and system running times Table 3 shows the IESA total/net revenues, system operating costs, MEUs' energy costs and system running times in scenario 2.
In Table 3, the corresponding IESA net income of each interval is consistent with the fluctuation trend of the total income. Combined with figure 13, when the element values are in Intervals 2 and 3, the IESA net income corresponding to ''all the backbone elements are less than 0'' is higher than that corresponding to ''all the backbone elements are greater than 0''. When the element values are in Intervals 5 and 6, the net gain of the IESA for ''all basic elements are less than 0'' is higher than that for ''all basic elements are greater than 0''. When the element values are in Intervals 7 and 8, the net income of the IESA corresponding to ''all peripheral elements are less than 0'' is higher than that corresponding to ''all peripheral elements are greater than 0''. When the element values are in Intervals 9 and 10, the IESA net income corresponding to ''all the weak peripheral elements are less than 0'' is higher than that corresponding to ''all the weak peripheral elements are greater than 0''. Therefore, when the values of the interval elements are all less than 0, this provides a positive contribution to the net income growth of the IESA; otherwise, this provides a negative contribution to the VOLUME 9, 2021 FIGURE 9. Uncertain optimal prices of electrical, heat, cold and natural gas energy determined by the IESA in scenario 2.
net income growth of the IESA. The comprehensive energy cost of MEUs showed no significant fluctuation. In addition, there is a large difference between the total return of Interval 6 and the total return of the other intervals, which may be because only the electric load cloud model is considered here or because the positive and negative elements of these two intervals have a certain randomness and fuzziness in the random sampling process.
(c) Retail energy pricing uncertainty analysis Figure 9 shows the uncertainty optimization strategy of the retail energy prices developed by the IESA for MEUs when only the elastic electric load cloud model is considered in scenario 2.
In figure. 9, When only the cloud model is considered to represent the uncertain characteristics of electricity users, the electricity sales strategy formulated by the IESA also has uncertain characteristics, which will indirectly lead to changes in the heat, cold and gas energy prices of the IESA in the overall optimization process. Under the condition of uncertain optimization operation, the IESA's energy sales price will change greatly with the selection of elements in different intervals. Obviously, the price curve of electricity sold in Interval 1 is between those in Intervals 2 and 3. When all values of elements are greater than or less than 0, the price of electricity sold in the IESA is too low or too high, while the prices of heat, cold and gas sold are consistent with the price of electricity sold. Similarly, Intervals 4, 5, and 6 have the same pricing pattern as Intervals 1, 2, and 3. By comparing the sales price curves of Intervals 7 and 8 and Intervals 9 and 10, it can also be found that when the element values extracted from each interval are both less than 0 or greater than 0, the sales price strategy formulated by the IESA is more reasonable.
(d) IDR uncertainty analysis Figure 10 shows MEUs participating in the uncertainty optimization strategy of the IDR when considering the elastic electric load cloud model of the industrial park in scenario 2.
In figure 10, Under the condition of certain optimization operation, MEUs use less energy during the peak period of the IESA energy sales price and try to use more energy during the trough of the energy purchasing price or normal period. The load of thermal energy users increases during the peak period, and the effect of the thermal load DR does not reach the desired effect of peak load clipping. This may be due to the actual energy demand decision at 11:00-16:00, or it may be due to the consideration of only the elastic load cloud model. In the case of uncertainty optimization, the IDR potential of MEUs will also be different when the values of the cloud droplets of the self-elasticity coefficient are located in different ranges. In Intervals 1, 2, and 3, the element values are the IDR potential > Interval 1 > Interval 3 within Interval 2. The distribution of the IDR potential in Intervals 4, 5, and 6 is basically consistent with the above scenario. Furthermore, through the analysis of the IDR potential uncertainty TABLE 3. Comparison of the revenues, costs and system running times from intervals one to ten in scenario 2. optimization results of Intervals 4,5,6,7,8,9 and 10, when the values of all elements in the cloud model are less than or greater than 0, there is an uncertainty optimization interval for the IDR potential of MEUs. However, if the values of the cloud droplet elements in the interval are all less than 0, the influence on the IDR potential of MEUs is positive, and vice versa.
(2) Uncertainty optimization strategy analysis with the elastic energy cloud model (a) Uncertainty analysis of revenues, costs, and system running times Table 4 shows the IESA total/net revenues, system operating costs, MEUs' energy costs and system running times for scenarios 1 to 7 in the park. In Table 4, scenario 1 shows the certain optimal operation results, in which the total/net benefit of the IESA is the minimum, and the integrated energy cost of MEUs is the highest. Scenario 2 considers only the elastic electric load cloud model. In this scenario, the total IESA/net income is higher than that in scenario 1, and the integrated energy cost of MEUs is lower than that in scenario 1. Scenarios 3∼7 comprehensively consider the elastic energy cloud models of electricity, heat, cold and gas. From the perspective of the total/net income of the IESA, it is found that the total/net income of the IESA is higher than that of scenario 2 under scenarios 3∼7, while the integrated energy cost of MEUs is approximately equal to or lower than that of scenario 2. The optimization result of scenario 3 is the closest to the real solution, in which the total/net income of the IESA is 11338.3/5987.5 ( ), and the integrated energy cost of MEUs is 5350.8 ( ). After each scenario is further reduced to the expected value, the total/net income of the IESA and the comprehensive energy cost of MEUs under scenario 7 are the closest to those under scenario 3, and the error rate is only approximately 0.9%, saving time about 1634.85s.
(b) Retail energy pricing uncertainty analysis Figure 11 shows MEUs participating in the uncertainty optimization strategy of the IDR when considering the elastic energy load cloud model of the industrial park for scenarios 1 to 7.
Under the conditions of scenario 1, the IESA's retail energy price fluctuates significantly. Under the circumstance of uncertainty optimization operation, the IESA's electricity, heat and gas revenue periods in scenario 2 are basically the same as those in scenario 1, in which the electricity sales price tends to be stable, the heat and gas sales prices continue to increase at the peak hour, and the cold sales prices transfer from the original prices from 05:00-08:00, 15:00 and 19:00-23:00 to those from 05:00-08:00, and the retail price changes most obviously. When the electrical, heat, cold and gas elastic energy cloud models are comprehensively considered in scenarios 3∼7, the traditional MC-RS method is adopted in scenarios 3, 4, and 5 to conduct the random sampling of cloud droplets with a self-elastic coefficient in the cloud model. Scenarios 6 and 7 are scenario reductions for the sample combination of scenario 5. When the expected scenario reduction is 200, the IESA's sales price fluctuation is almost consistent with that of Scenario 3, and the system operation time is the shortest.
(c) IDR uncertainty analysis Figure 12 shows MEUs participating in the uncertainty optimization strategy of IDR when considering the elastic energy cloud model of the industrial park in scenarios 1 to 7. Figure 12 combined with figure. 10 and Table 4 shows that the IDR potential obtained through the optimization of scenarios 3∼7 is significantly better than that through  scenarios 1 and 2. The traditional MC-RS method is adopted in scenarios 3, 4 and 5 for random combination sampling of elements in the cloud model. Scenario 3 has the largest number of samples, and the optimization result for this scenario can be considered to be the closest to the real solution. In scenario 7, the SBR method is adopted to reduce the random combination sampling scenario and greatly reduces the system running time. In addition, the integrated energy use of MEUs in scenario 7 is very similar to that in scenario 3. Therefore, this conclusion also verifies the effectiveness of the uncertainty optimization strategy proposed in this paper.
In fact, the potential numerical value of MEUs participating in IDR can also be represented by the elastic energy cloud model at multiple time scales. For example, the elastic energy cloud model obtained for scenario 7 corresponding to figure 12 is presented in figure 15. In addition, figure 15 also shows the actual load IDR potential characteristics of MEUs.

VI. SUMMARY AND CONCLUSION
This paper proposes an optimal economic dispatch strategy for an industrial park with consideration of an IDR uncertainty elastic energy cloud model under the background of China's energy market reform. Based on the interval contribution of the electric load cloud model, a new SRS-IC method is proposed to solve this model. Moreover, a typical industrial park is used to verify the effectiveness of the method. The conclusions of this paper are summarized as follows: a) An elastic energy cloud model with electricity, heat, cold and gas MEUs is proposed to describe the randomness and fuzziness of elastic energy in the industrial park.
b) A new concept of an elastic electrical load cloud model based on interval contribution is proposed, and this model can effectively analyze the overall economic benefit of an industrial park and the contribution of MEUs to the IDR potential.
c) Compared with the traditional MC-RS method, a new SRS-IC method is proposed, which can not only save about 1634s of system running time but can also ensure that the error rate of the sampling results is about 0.9%.
d) The multiple uncertainty characteristics of the elastic energy, are fully considered and combined with the above SRS-IC method to propose a new uncertainty optimal dispatch strategy for an industrial park.
e) The above uncertainty optimization strategy can improve the cooperation environment for the IESA and MUEs, and the economic benefits of the IESA, and can also optimize the energy use structure of the industrial park and further explore the IDR potential of MEUs.

VII. FUTURE WORK
Currently, an elastic energy cloud model that can represent both randomness and fuzziness has been established, and a correlation analysis of multiple energy cloud models has been carried out. Therefore, determining how to comprehensively evaluate the value of elastic energy in an industrial park from the perspectives of technology, economy and benefit will be considered. Moreover, based on the value evaluation of elastic energy, determining how to further design the integrated demand management mechanism and the operation mechanism of a multi-energy market will be a focus in the future.

APPENDIX
The quadratic function matching method is as follows: The vertex coordinates can be expressed as: In formula (14), the total revenue of electric power users after participating in DR is taken as an example, and the formula can be derived as follows: where, by combining (22), (23) and (14), the functional relationship between the retail electricity price change and the initial reference price is obtained as follows: where, by combining (22), (23) and (14), the functional relationship between the actual elastic electrical load after participating in the demand response, the predicted value of the initial elastic electrical load, the change in the retail electricity price, the initial reference price, and the self-elasticity coefficient can be obtained as follows: Substituting (A-3) and (A-6) into (A-1), we can further obtain the quadratic function expression of the total return after the participation of power users in DR. The details are as follows: The total revenue of power users after the implementation of DR can be expressed indirectly as follows:  Set x = κ t dl,prc,ref , a = ε t dl λ dl , b = 1 + ε t dl λ dl , c = 1, ε t dl < 0, and λ dl = 0.65.  Table 3.  Table 4.
By referring to the quadratic function substitution method, the vertex coordinates can be expressed as follows: where: According to the distribution state of cloud drops with the self-elastic coefficient, it is necessary to reasonably analyze the approximate range of the x-coordinate of the vertex. The details are as follows: − Suppose the general expression of quadratic piecewise linear function f (x) is: x min = κ t dl,prc,ref ,min , x max = κ t dl,prc,ref ,max , where ς 1 , ς 2 , µ 1 , and µ 2 can be expressed as follows: where ϕ = κ t dl,prc,ref β. Then, add the constraints of ϕ as follows: where M is a relatively large constant. Note 1: Similarly, for the convenience of description, this paper adopts the subscripts dl, hl, cl, and gl to represent the piecewise linear functions of electricity, heat, cold and gas, respectively.
Note 2: Under the analysis of the uncertainty optimization scenario, the fuzzy parameters in each parameter vector can be expressed by the superscript ∼ in combination with formula (38).
Definition 1 [26], [30]: Let U be the universe of discourse and C be a qualitative concept related to U . If x ∈ U is a random instantiation of concept C and µ (x) ∈ [0, 1] is the certainty degree of x belonging to C, which satisfies µ (x) : U ∈ [0, 1] , ∀x ∈ U then, a one-to-many mapping can be defined as a normal cloud, and (x, µ (x)) can be called a cloud drop.
Definition 2 [26], [30]: Let U be the universe of discourse and C be a qualitative concept related to U . C contains three numerical characteristics: expection Ex, entropy En, and hyper entropy He. If x ∈ U is a random instantiation of concept C and µ (x) ∈ [0, 1] is the certainty degree of x belonging to C, which satisfies µ (x) : U ∈ [0, 1] , ∀x ∈ U then, ε = R N (Ex, |y|), y = R N (En, He), and the membership function satisfies the following form: Then, the distribution of the random variable X composed of all cloud drops is called the normal cloud model. See Figures 13-15.