An Interval Optimization-Based Approach for Electric–Heat–Gas Coupled Energy System Planning Considering the Correlation between Uncertainties

Abstract

In this paper, a novel methodological framework for energy hub (EH) planning, considering the correlation between renewable energy source (RES) and demand response (DR) uncertainties, is proposed. Unlike other existing works, our study explicitly considers the potential correlation between the uncertainty of integrated energy system operations (i.e., wind speed, light intensity, and demand response). Firstly, an EH single-objective interval optimization model is established, which aims at minimizing investment and operation costs. The model fully considers the correlation between various uncertain parameters. Secondly, the correlation between uncertainties is dealt with by the interval models of multidimensional parallelism and affine coordinate transformation, which are transformed into a deterministic optimization problem by the interval order relationship and probability algorithm, and then solved by a genetic algorithm. Finally, an experimental case is analyzed, and the results show that the research method in this paper has good engineering practicability. At the same time, different correlations among uncertainties have different influences on integrated energy system planning. Correlation and influence are positively correlated.

1. Introduction

In recent years, with increasingly prominent environmental problems and the gradual depletion of energy, it has become an irresistible trend to make full use of RES to generate electricity. In the past, fossil energy was the main source of electricity generation. Compared with the traditional electric energy utilization mode, the integrated energy system realizes flexible conversion and unified management of different energy sources, with the help of advanced conversion and storage equipment. It provides a new solution for carbon emission control and economic construction [1].

In order to achieve the optimal economic effect, equipment selection and installed capacity of the system should be scientifically planned. Therefore, there are many studies on integrated energy system planning. The concept of an energy hub is put forward. It realizes the standardized modeling of an integrated energy system and lays an important foundation for planning research [2]. An energy system composed of electric, gas, cold and heat loads, and combined heat and power, boiler and cooling unit, can achieve the best energy management [3]. An integrated energy system combining wind and solar could provide a more sustainable way to increase the utilization of renewable energy [4]. Rational use of an environmental energy power supply and cogeneration can reduce the loss of energy and damage to the environment [5]. At the same time, economic optimization can be achieved. Considering the equipment selection and capacity planning of the integrated energy system, coupled with electricity, heat and gas [6], and the comprehensive energy optimization model, an operating cost is established [7,8]. According to the above literature, these studies all ignore the impact of parameter uncertainty on planning. The randomness and volatility of renewable energy source output and the uncertainty of demand response caused by users’ energy usage habits have a direct impact on integrated energy system planning. Therefore, in order to make up for the above deficiencies, the new type of multi-energy integrated system should take into account demand response to electricity and thermal energy, as well as the uncertainties of RES, and optimize the system by modeling energy conversion, storage, production, and consumption [9]. The input energy of the hub includes electricity and natural gas, and the output demand includes heat and electricity. The impact of randomness on wind power generation and hub power demand is evaluated [10]. Considering the uncertainty of the system, the energy hub will change the optimal operating point to meet consumers’ demand for electricity and heat [11]. The influence of uncertainties related to wind power generation and load forecasting on energy hub dispatching can be dealt with by adopting a stochastic programming method [12]. A stochastic programming method is adopted to solve different uncertainties, and a sustainable framework for optimizing the design of energy hubs is established [13].

Research on the uncertainty of renewable energy and demand response has been relatively sufficient, such as stochastic theory and fuzzy theory, but few studies consider the correlation between parametric uncertainties on the planning and operation of integrated energy systems uncertainties. However, because the uncertainties of RES output will directly affect the electricity selling price of the system operator, the energy use of consumers will be affected. Although the above viewpoints have been discussed in previous studies, they have not been fully considered in energy system planning. It is necessary to fully consider the high correlation between the uncertainty of DR and RES. Only by considering this more realistic feature can we obtain a more scientific and persuasive planning result than exists in current research.

In order to solve the above problems, a single objective interval optimization planning method for integrated energy systems, which considers the correlation between uncertainties, is proposed.

The main innovation points of this paper are as follows:

• An interval optimization approach for electric–heat–gas coupled energy system planning is proposed. The optimization of system economic benefits can be achieved by considering the influence of the uncertainty correlation between renewable energy sources and demand response.
• A new method was used to solve the correlation of the uncertainties of renewable energy sources and DR. The correlation between uncertainties is dealt with by a multidimensional parallelepiped interval model and affine coordinate transformation.

The rest of this article is as follows. The second section introduces the EH planning framework, considering uncertainty. The third section introduces the uncertainty modeling that was carried out for wind generator (WG), photovoltaic (PV) and load. The correlation model between uncertainties is also provided. The fourth section introduces the EH programming model of multi-objective interval optimization by considering the correlation between uncertainties, and proposes the specific objective function and constraint conditions. The fifth section introduces the corresponding solution strategy. The sixth section introduces the example study. The seventh section introduces the summary.

2. EH Planning Framework Considering Uncertainty Correlation

In this paper, a typical park-level EH containing electricity, heat and gas is taken as an example to carry out research. The basic structure is shown in Figure 1. As shown in Figure 1, the input end of the EH system is connected to the external power grid and the natural gas network respectively, and the output end is connected to users on the demand side. During normal operation, the EH uses all kinds of internal energy conversion equipment (including combined heat and power (CHP), gas boiler (GB) and electric boiler (EB)) to convert primary input energy into the electrical energy and heat energy required by the system, and it cooperates with wind, light and other RES at the same time so as to realize the energy supply to end users. In addition, in order to improve EH operation efficiency, it is assumed that the system is also equipped with electric energy storage (ES), thermal energy storage (TS) and an intelligent energy meter [14]. According to the demand response characteristics, all appliances of users are divided into uncontrollable load, time shiftable load, and energy fungible load. Among them, uncontrollable load is the load that cannot respond to demand, such as lighting facilities, televisions, and so on. The adjustment of time shiftable load energy usage time will not affect user comfort, so users can flexibly adjust the usage time of time shiftable load. Such loads include washing machines, electric cars, etc. In reality, there are two types of energy fungible load. Firstly, there are household appliances that can switch between multiple energy types (natural gas or electricity), including electric-gas mixed air-conditionings, and residential kitchen equipment, such as electric-gas mixed stoves or ovens. Secondly, some appliances, such as gas stoves and microwave ovens, consume different types of energy but achieve the same function (heating and cooking). Therefore, they can replace each other and are equivalent to the energy fungible load as a whole. Similar examples are wall-hung boilers and electric heaters.

In this paper, the correlation between the parameter uncertainties of the system is considered. By adopting an interval method to model the uncertainty, including the WG, PV output uncertainty and the demand response uncertainty, the decision maker can optimize the type and the installed capacity of the equipment in EH, and realize the economic operation of the system. Therefore, this paper proposes a single-objective interval optimization (SOIO) model for EH programming which considers the uncertainty correlation. See Figure 2.

3. Uncertainty and Correlation Modeling

3.1. Uncertainty

The accurate modeling of parameter uncertainty is the foundation of this paper. In this section, the WG and PV output uncertainty and load uncertainty considering the demand response will be modeled.

3.1.1. Wind Generator

The output power of a wind generator is mainly affected by wind speed at its position. In nature, because wind speed is random and intermittent, the output of a wind generator is an uncertain variable. For this reason, the load factor ${\tilde{k}}_t^{\mathrm{WG}}$ is defined in this paper to represent the ratio between the actual generating power of the wind generator and its rated capacity under wind speed at time period $t$. It is an interval variable describing the variation of the WG load rate at random wind speed, ${\tilde{k}}_t^{\mathrm{WG}}=[{\underset{\_}{k}}_t^{\mathrm{WG}},{\overline{k}}_t^{\mathrm{WG}}]$. Then, wind power output can be expressed as:

$${\tilde{P}}_t^{\mathrm{WG}}={\tilde{k}}_t^{\mathrm{WG}}{M}^{\mathrm{WG}}$$

(1)

where

• ${\tilde{P}}_t^{\mathrm{WG}}$, generating power of WG;
• ${\tilde{k}}_t^{\mathrm{WG}}$, WG load factor;
• $M^{\mathrm{WG}}$, installed capacity of WG.

3.1.2. Photovoltaic Power Generation

The intensity of solar light and ambient temperature are key factors affecting the output power of photovoltaic power generation [15]. Considering the uncertainty of temperature and light in nature, similar to the WG model, the operating characteristics of photovoltaic power generation can be expressed by Formula (2). ${\tilde{k}}_t^{\mathrm{PV}}=[{\underset{\_}{k}}_t^{\mathrm{PV}},{\overline{k}}_t^{\mathrm{PV}}]$ is an interval variable representing the random variation of PV load rate.

$${\tilde{P}}_t^{\mathrm{PV}}={\tilde{k}}_t^{\mathrm{PV}}{M}^{\mathrm{PV}}$$

(2)

where

• ${\tilde{P}}_t^{\mathrm{PV}}$, generating power of PV;
• ${\tilde{k}}_t^{\mathrm{PV}}$, PV load factor;
• $M^{\mathrm{PV}}$, installed capacity of PV.

3.1.3. Load Modeling

Social human factors are often taken into account in load modeling [16]. In the RTP-based DR project, on the premise of maintaining the same comfort level, the user adjusts energy use demand voluntarily, according to the dynamic price signal provided by the system operator, to achieve positive interaction with the system. For example, consumers take the initiative to reduce unnecessary electricity consumption during periods of higher electricity prices. From the perspectives of demand nature and adjustable potential, EH loads can be roughly divided into three categories [17,18,19]: uncontrolled load (UL), time shiftable load (TSL), and energy fungible load (EFL). The above types of loads are modeled as follows.

• Uncontrollable Load

UL is a load of high energy importance and unresponsive to RTP. Common ULs in EH include lighting, heating, etc. Its operation characteristics can be expressed in terms (3) and (4).

$$P_t^{\mathrm{UL}}={\gamma }_t^{\mathrm{UL}-\mathrm{e}}{P}_t^{\mathrm{DE},0}$$

(3)

$$H_t^{\mathrm{UL}}={\gamma }_t^{\mathrm{UL}-\mathrm{h}}{H}_t^{\mathrm{DE},0}$$

(4)

where

• $P_t^{\mathrm{UL}}-H_t^{\mathrm{UL}}$, electrical–thermal load demand of UL under RTP;
• ${\gamma }_t^{\mathrm{UL}-\mathrm{e}}-{\gamma }_t^{\mathrm{UL}-\mathrm{h}}$, proportion of UL in total electricity–heat load demand;
• $P_t^{\mathrm{DE},0}-H_t^{\mathrm{DE},0}$, the electric–thermal load value under the base price.

2.

Time Shiftable Load

TSL refers to a load refers to the total energy consumption is fixed, but can flexibly change energy use time. Common TSLs include water heaters, washing machines, electric cars, etc.

In the RTP-based DR project, the user adjusts the running time of TSL according to the dynamic electricity price obtained. In actual operation, as users have different living habits and behavioral preferences, the parameters ${\tilde{\epsilon }}_t^{\mathrm{TSL}}$ and ${\tilde{\epsilon }}_{tt\prime }^{\mathrm{TSL}}$ can be defined as uncertain variables. In this paper, the interval method is adopted to model the uncertainty. So, ${\tilde{\epsilon }}_t^{\mathrm{TSL}}=[{\underset{\_}{\epsilon }}_t^{\mathrm{TSL}},{\overline{\epsilon }}_t^{\mathrm{TSL}}]$, ${\tilde{\epsilon }}_{tt\prime }^{\mathrm{TSL}}=[{\underset{\_}{\epsilon }}_{tt\prime }^{\mathrm{TSL}},{\overline{\epsilon }}_{tt\prime }^{\mathrm{TSL}}]$. The load transfer process can be modeled with price elasticity [20] as:

$$\begin{array}{r}{\tilde{P}}_t^{\mathrm{TSL}}={\gamma }_t^{\mathrm{TSL}}{P}_t^{\mathrm{DE},0}\left[1+\frac{{\tilde{\epsilon }}_t^{\mathrm{TSL}}\left({\rho }_t^{\mathrm{e}}-{\rho }_t^{\mathrm{e},0}\right)}{{\rho }_t^{\mathrm{e},0}}\right]+\\ \left[{\displaystyle \sum _{t\prime \in T,t\prime \ne t}\frac{{\gamma }_{t\prime }^{\mathrm{TSL}}{P}_{t\prime }^{\mathrm{DE},0}{\tilde{\epsilon }}_{tt\prime }^{\mathrm{TSL}}\left({\rho }_{t\prime }^{\mathrm{e}}-{\rho }_{t\prime }^{\mathrm{e},0}\right)}{{\rho }_{t\prime }^{\mathrm{e},0}}}\right]\end{array}$$

(5)

where

• ${\tilde{P}}_t^{\mathrm{TSL}}$, electrical load demand of TSL under RTP;
• ${\gamma }_t^{\mathrm{TSL}}$, proportion of TSL in the total electricity load demand of EH;
• ${\rho }_t^{\mathrm{e},0}-{\rho }_t^{\mathrm{e}}$, base electricity price-RTP electricity price;
• ${\tilde{\epsilon }}_t^{\mathrm{TSL}}$, TSL price self-elasticity coefficient;
• ${\tilde{\epsilon }}_{tt\prime }^{\mathrm{TSL}}$, TSL price cross-elasticity coefficient.

3.

Energy Fungible Load

EFL refers to a load that has a fixed energy time but can flexibly choose the form of energy (electricity or gas) according to needs. In EH, typical EFLs include air conditioning equipment based on electricity–gas mixed refrigeration, resident kitchen equipment, etc. In the RTP-based DR project, the user can choose the EFL form according to the real-time price of different energy in each time period, without losing comfort level. The process can be modeled by alternatives in microeconomics [21] as:

$${\tilde{P}}_t^{\mathrm{EFL}}={\gamma }_t^{\mathrm{EFL}-\mathrm{e}}{P}_t^{\mathrm{DE},0}\left[1+\frac{{\tilde{\epsilon }}_t^{\mathrm{EFL}}\left({\rho }_t^{\mathrm{e}}-{\rho }^{\mathrm{h}}\right)}{{\rho }^{\mathrm{h}}}\right]$$

(6)

$${\tilde{H}}_t^{\mathrm{EFL}}={\gamma }_t^{\mathrm{EFL}-\mathrm{h}}{H}_t^{\mathrm{DE},0}-{\delta }^{\mathrm{EFL}}\left({\tilde{P}}_t^{\mathrm{EFL}}-{\gamma }_t^{\mathrm{EFL}-\mathrm{e}}{P}_t^{\mathrm{DE},0}\right)$$

(7)

where

• ${\tilde{P}}_t^{\mathrm{EFL}}-{\tilde{H}}_t^{\mathrm{EFL}}$, electrical–thermal load demand of EFL under RTP;
• ${\gamma }_t^{\mathrm{EFL}-\mathrm{e}}-{\gamma }_t^{\mathrm{EFL}-\mathrm{h}}$, proportion of EFL in the total electricity–heat load demand of EH;
• ${\tilde{\epsilon }}_t^{\mathrm{EFL}}$, substitute price elasticity coefficient;
• ${\delta }^{\mathrm{EFL}}$, electric–thermal conversion efficiency of EFL;
• ${\rho }^{\mathrm{h}}$, hot selling price for EH to end users.

4.

Total Load Demand of System includes DR

In Formulas (3)–(7), the total EH load in each period under the DR mode can be calculated as follows:

$${\tilde{P}}_t^{\mathrm{DE}}=P_t^{\mathrm{UL}}+{\tilde{P}}_t^{\mathrm{TSL}}+{\tilde{P}}_t^{\mathrm{EFL}}$$

(8)

$${\tilde{H}}_t^{\mathrm{DE}}=H_t^{\mathrm{UL}}+{\tilde{H}}_t^{\mathrm{EFL}}$$

(9)

where

• ${\tilde{P}}_t^{\mathrm{DE}}-{\tilde{H}}_t^{\mathrm{DE}}$, electrical-thermal load value of the system under RTP.

3.2. Correlation

As the uncertainty of RES output will directly affect the electricity selling price of the system operator, and thus affect the user’s energy usage, there is a correlation between uncertainties. An integrated energy system considering the uncertainty and correlation of RES and DR can effectively reduce the influence of fluctuations of external factors on EH, and thus increase the overall resilience of the energy hub [22,23,24,25]. In this paper, a multidimensional parallelepiped interval model is used to deal with the correlation between uncertainties [26]. For a q-dimension problem of general form, the uncertainty field is, correspondingly, a multidimensional parallelism. For each uncertain variable, its edge interval is shown in Equation (10), and the correlation coefficient ρ_ij is used to describe the correlation between two interval variables, as shown in Equation (11).

$$U_{\mathrm{i}}^{\mathrm{I}}=[U_i^L,U_i^R], i=1,2\dots q$$

(10)

$${\rho }_{ij}=\frac{U_j^w}{U_i^w\tan {\theta }_{ij}}$$

(11)

where, θ_ij is the correlation angle between any two variables, U_i and U_j; ρ_ij range in [−1, 1]. When ρ_ij = 1, the variables U_i and U_j are completely positively correlated. When ρ_ij = 0, variables U_i and U_j are independent of each other. When ρ_ij = −1, the variables U_i and U_j are completely negatively correlated.

4. EH Programming Model with Multi-Objective Interval Optimization Considering the Correlation between Uncertainties

In this paper, a single objective optimization model considering parameter correlation is constructed based on the multidimensional parallelepiped interval model [26]; the general expression is as follows:

$$\{\begin{array}{l}\underset{\mathit{X}}{\min } f(\mathit{X},\mathit{U})\\ \mathsf{s}.\mathsf{t}.h_i(\mathit{X},\mathit{U})\le {\tilde{b}}_i=[b_i^{\min },b_i^{\max }],i=1,2,\dots ,l,X \in {\Omega }^{\mathrm{n}}\\ \mathit{U}\in \Omega ({\mathsf{U}}^I,\rho )\end{array}$$

(12)

where, f(X,U) represents the vector consisting of the objective function; h_i(.) is inequality constraint condition; ${\tilde{b}}_i$ is the uncertain constraint for allowable value interval $i$; and $I$ represents the number of constraints in the model. In addition, X and U, respectively, represent the vector space composed of optimization variables and interval variables. U is a q-dimensional uncertainty vector, whose uncertainty domain $\Omega$ is a geometrically multidimensional parallelohedron, which is determined by edge interval vector U^I and correlation coefficient vector $\rho$.

4.1. Objective Function

As mentioned above, this paper takes the economy of the system as the target and constructs the EH interval optimization programming model, as shown in Equation (13). The physical meaning of Equation (16) is: for EH investors, the cost of implementing DR project $C^{\mathrm{DR}}$ is equal to the difference between the total revenue before and after the introduction of DR (generated by providing power supply and heating services to users). If the total operating income of EH investors after the implementation of DR is less than the income before the implementation of DR, the change is in the form of “cost”. On the contrary, it is manifested as a “benefit” form (namely, negative cost). In order to reduce the complexity of the research, it is assumed that the electricity–gas purchase price of EH from the external energy market is established based on long-term contracts and controlled by the government, so the impact of the above price uncertainty can be ignored in the EH planning model. Taking into account the price characteristics of different energies, the wholesale market price of EH’s purchasing of electricity from the external power network is set as a real-time price, and the price of gas purchased from the external gas network is set as a fixed price.

The physical meaning of Equation (17) is: “carbon emission tax [27]”, which can be used to transform the carbon emission of the environmental indicator system into a cost economic indicator to realize the unification of environmental and economic attributes in the measurement.

$$\min f_1=C^{\mathrm{inv}}+C^{\mathrm{opt}}+C^{\mathrm{DR}}+C_{\mathrm{t}}^{ctp}$$

(13)

$$C^{\mathrm{inv}}={\displaystyle \sum _{i\in I}\frac{r{(1+r)}^{y_i}}{{(1+r)}^{y_i}-1}}{c}_i^{\mathrm{inv}}{M}_i$$

(14)

$$C^{\mathrm{opt}}={\displaystyle \sum _{i\in I}{c}_i^{\mathrm{mai}}{M}_i}+\tau {\displaystyle \sum _{t\in T}\Delta t(c_t^{\mathrm{ele}}{P}_t^{\mathrm{ele}}+c_t^{\mathrm{gas}}{G}_t^{\mathrm{gas}})}$$

(15)

$$C^{\mathrm{DR}}=\tau {\displaystyle \sum _{t\in T}\Delta t\left\{\begin{array}{l}[{\rho }_t^{\mathrm{e},0}{P}_t^{\mathrm{DE},0}({\gamma }_t^{\mathrm{TSL}}+{\gamma }_t^{\mathrm{EFL}-\mathrm{e}})+{\rho }^{\mathrm{h}}{\gamma }_t^{\mathrm{EFL}-\mathrm{h}}{H}_t^{\mathrm{DE},0}]-\\ [{\rho }_t^{\mathrm{e}}({\tilde{P}}_t^{\mathrm{TSL}}+{\tilde{P}}_t^{\mathrm{EFL}})+{\rho }^{\mathrm{h}}{\tilde{H}}_t^{\mathrm{EFL}}]\end{array}\right\}}$$

(16)

$$C_{\mathrm{t}}^{ctp}=B_{\lambda }^{ct}(e{g}^{gsp}{P}_t^{gsp}{r}_t+e{g}^{gt}{\displaystyle \sum _{g\in {\Omega }_G^{gt}}{P}_{g,t}^{gt}{r}_t})$$

(17)

where

• $C^{ \mathrm{inv}}$, annual value of EH total investment cost;
• $C^{ \mathrm{opt}}$, annual operating cost;
• $C^{\mathrm{DR}}$, cost of DR project;
• $C_{\mathrm{t}}^{ctp}$, carbon emission cost;
• $c_t{}^{\mathrm{inv}}$, investment cost of equipment per unit capacity;
• $c_i^{\mathrm{mai}}$, annual fixed maintenance cost of equipment;
• $c_t{}^{\mathrm{ele}}-c_t{}^{\mathrm{gas}}$, purchase electricity–gas price for EH;
• $e{g}^{gsp}$, the comprehensive carbon emission intensity of the power grid;
• $e{g}^{gt}$, the carbon emission per unit capacity of the micro gas turbine;
• $B_{\lambda }^{ct}$, carbon tax rate applied in the λ year;
• $r$, discount rate, which is set as 8% in this paper;
• $\Delta t$, duration of a single period, which is set as 1 h in this paper;
• $\tau$, number of days in a year;
• $y$, service life of the equipment;
• $T$, EH running period;
• ${\gamma }_t^{\mathrm{EFL}-\mathrm{e}}-{\gamma }_t^{\mathrm{EFL}-\mathrm{h}}$, proportion of EFL in the total electricity–heat load demand of EH.

4.2. Constraints

4.2.1. System Configuration Constraint

The system configuration constraint requires that the maximum installed capacity of all kinds of equipment in EH should not exceed a certain limit to meet the limitation of physical space

$$0\le M_i\le M_i^{\max }, \forall i\in I$$

(18)

where

• $M$, installed capacity for equipment;
• $M^{\max }$, upper limit installed capacity for equipment.

4.2.2. Pricing Constraint

In order to avoid the negative impact of excessive RTP fluctuations on the interests of users, the price level should be controlled within a certain range. On the other hand, in order to ensure that the demand side is willing to participate in the DR project, users’ expected total energy use cost under RTP should be lower than their cost expenditure under the base electricity price mode.

$${\rho }^{\mathrm{e},\min }\le {\rho }_t^{\mathrm{e}}\le {\rho }^{\mathrm{e},\max },\forall t\in T$$

(19)

$$\begin{array}{c}{\displaystyle \sum _{t\in T}\left[{\rho }_t^{\mathrm{e},0}{P}_t^{\mathrm{DE},0}({\gamma }_t^{\mathrm{TSL}}+{\gamma }_t^{\mathrm{EFL}-\mathrm{e}})+{\rho }^{\mathrm{h}}{\gamma }_t^{\mathrm{EFL}-\mathrm{h}}{H}_t^{\mathrm{DE},0}\right]}\ge \\ {\displaystyle \sum _{t\in T}\left[{\rho }_t^{\mathrm{e}}({\tilde{P}}_t^{\mathrm{TSL}}+{\tilde{P}}_t^{\mathrm{EFL}})+{\rho }^{\mathrm{h}}{\tilde{H}}_t^{\mathrm{EFL}}\right]}\end{array}$$

(20)

where

• ${\rho }^{\mathrm{e},\max }-{\rho }^{\mathrm{e},\min }$, the upper or lower limits of RTP allowable fluctuation range

4.2.3. Security Constraint

System security constraints mainly include EH’s internal electricity–heat–gas power balance constraints (21)–(23) and EH’s electricity–gas purchase quantity constraints (24) and (25) from an external energy network:

$$P_t^{\mathrm{ele}}+P_t^{\mathrm{CHP}}+{\tilde{P}}_t^{\mathrm{WG}}+{\tilde{P}}_t^{\mathrm{PV}}+P_t^{\mathrm{ES}-\mathrm{dch}}={\tilde{P}}_t^{\mathrm{DE}}+P_t^{\mathrm{EB}}+P_t^{\mathrm{ES}-\mathrm{ch}}, \forall t\in T$$

(21)

$$H_t^{\mathrm{CHP}}+H_t^{\mathrm{GB}}+H_t^{\mathrm{EB}}+H_t^{\mathrm{TS}-\mathrm{dch}}={\tilde{H}}_t^{\mathrm{DE}}+H_t^{\mathrm{TS}-\mathrm{ch}},\forall t\in T$$

(22)

$$G_t^{\mathrm{gas}}=G_t^{\mathrm{CHP}}+G_t^{\mathrm{GB}} , \forall t\in T$$

(23)

$$0\le P_t^{\mathrm{ele}}\le P^{\mathrm{ele}-\max },\forall t\in T$$

(24)

$$0\le G_t^{\mathrm{gas}}\le G^{\mathrm{gas}-\max },\forall t\in T$$

(25)

where

• $P^{\mathrm{CHP}}$, electrical power produced by CHP;
• $P^{\mathrm{EB}}$, power consumed by EB;
• $G^{\mathrm{CHP}}-G^{\mathrm{GB}}$, gas consumption of CHP-GB;
• $G^{\mathrm{gas}}$, gas purchase quantity from the external market;
• $G^{\mathrm{gas}-\max }$, the upper limit of EH’s gas purchase quantity;
• $P^{\mathrm{ES}-\mathrm{ch}}-P^{\mathrm{ES}-\mathrm{dch}}$, charge–discharge power of ES;
• $H^{\mathrm{TS}-\mathrm{ch}}-H^{\mathrm{TS}-\mathrm{dch}}$, charge–discharge power of TS;
• $H^{\mathrm{CHP}}-H^{\mathrm{EB}}-H^{\mathrm{GB}}$, heat generation power of CHP-EB-GB;
• $P^{\mathrm{ele}}$, purchasing power from the external market;
• $P^{\mathrm{ele}-\max }$, the upper limit of EH’s purchasing power.

4.2.4. Constraints on Operating Characteristics of Components

• Cogeneration Unit

The generating power of a cogeneration unit is a primary function of the amount of natural gas it consumes and is less than its configured capacity. The thermoelectric power of a cogeneration unit is limited by its thermoelectric ratio.

$$P_t^{\mathrm{CHP}}={\eta }^{\mathrm{e}}{G}_t^{\mathrm{CHP}}{v}^{\mathrm{gas}}$$

(26)

$$0\le P_t^{\mathrm{CHP}}\le M^{\mathrm{CHP}}$$

(27)

$$H_t^{\mathrm{CHP}}=P_t^{\mathrm{CHP}}{\beta }^{\mathrm{CHP}}$$

(28)

where

• $v^{\mathrm{gas}}$, calorific value of natural gas;
• $M^{\mathrm{CHP}}$, installed capacity for CHP;
• ${\beta }^{\mathrm{CHP}}$, CHP heat-electric ratio;
• ${\eta }^{\mathrm{e}}$, CHP generation efficiency.

2.

Electric–Gas Boilers

Electric–gas boilers use electricity and natural gas as energy sources and convert them into heat. For the above equipment, the relationship between heat generation power and power–gas consumption can be expressed by Formulas (29)–(30). In addition, the maximum heat generation power is limited by installed capacity, as shown in inequalities (31)–(32).

$$H_t^{\mathrm{EB}}={\eta }^{\mathrm{EB}}{P}_t^{\mathrm{EB}}$$

(29)

$$H_t^{\mathrm{GB}}={\eta }^{\mathrm{GB}}{G}_t^{ \mathrm{GB}}{v}^{\mathrm{gas}}$$

(30)

$$0\le H_t^{\mathrm{EB}}\le M^{\mathrm{EB}}$$

(31)

$$0\le H_t^{\mathrm{GB}}\le M^{\mathrm{GB}}$$

(32)

where

• $M^{\mathrm{EB}}-M^{\mathrm{GB}}$, installed capacity of electric-gas boiler;
• ${\eta }^{\mathrm{EB}}-{\eta }^{\mathrm{GB}}$, thermal efficiency of electric-gas boiler.

3.

Electrical–Thermal Energy Storage Equipment

In EH operation, the relationship between the energy storage state of the electrical–thermal energy storage device and the change of charge–discharge-thermal power can be expressed by Equations (33) and (34). In order to prevent excessive charging and releasing of the equipment, SOC constraints should be applied, as shown in Equations (35) and (36). In addition, the charging and discharging power constraints of the energy storage equipment are shown in Equations (37) and (38), in which charge and discharge electricity–heat cannot be carried out at the same time. Finally, in order to ensure the sustainable operation of the energy storage equipment, it is required that the energy storage state at the beginning and the end should remain unchanged within a scheduling cycle, as shown in Equations (39) and (40).

$$E_t^{\mathrm{ES}}=\left(1-{\lambda }^{\mathrm{ES}}\right)E_{t-1}^{\mathrm{ES}}+\Delta t\left(P_t^{\mathrm{ES}-\mathrm{ch}}{\eta }^{\mathrm{ES}-\mathrm{ch}}-\frac{P_t^{\mathrm{ES}-\mathrm{dch}}}{{\eta }^{\mathrm{ES}-\mathrm{dch}}}\right)$$

(33)

$$E_t^{\mathrm{TS}}=\left(1-{\lambda }^{\mathrm{TS}}\right)E_{t-1}^{\mathrm{TS}}+\Delta t\left(H_t^{\mathrm{TS}-\mathrm{ch}}{\eta }^{\mathrm{TS}-\mathrm{ch}}-\frac{H_t^{\mathrm{TS}-\mathrm{dch}}}{{\eta }^{\mathrm{TS}-\mathrm{dch}}}\right)$$

(34)

$${\mu }^{\mathrm{ES}-\min }{M}^{\mathrm{ESE}}\le E_t^{\mathrm{ES}}\le {\mu }^{\mathrm{ES}-\max }{M}^{\mathrm{ESE}}$$

(35)

$${\mu }^{\mathrm{TS}-\min }{M}^{\mathrm{TSE}}\le E_t^{\mathrm{TS}}\le {\mu }^{\mathrm{TS}-\max }{M}^{\mathrm{TSE}}$$

(36)

$$0\le P_t^{\mathrm{ES}-\mathrm{ch}}\le {\alpha }_t^{\mathrm{ES}-\mathrm{ch}}{\xi }^{\mathrm{ES}}{M}^{\mathrm{ESE}}$$

(37)

$$0\le H_t^{\mathrm{TS}-\mathrm{ch}}\le {\alpha }_t^{\mathrm{TS}-\mathrm{ch}}{\xi }^{\mathrm{TS}}{M}^{\mathrm{TSE}}$$

(38)

$$0\le P_t^{\mathrm{ES}-\mathrm{dch}}\le {\alpha }_t^{\mathrm{ES}-\mathrm{dch}}{\xi }^{\mathrm{ES}}{M}^{\mathrm{ESE}}$$

(39)

$$0\le H_t^{\mathrm{TS}-\mathrm{dch}}\le {\alpha }_t^{\mathrm{TS}-\mathrm{dch}}{\xi }^{\mathrm{TS}}{M}^{\mathrm{TSE}}$$

(40)

where

• $M^{\mathrm{ESE}}-M^{\mathrm{TSE}}$, installed capacity of ES-TS;
• ${\xi }^{\mathrm{ES}}-{\xi }^{\mathrm{TS}}$, power–capacity ratio of ES and TS;
• ${\lambda }^{\mathrm{ES}}-{\lambda }^{\mathrm{TS}}$, self-discharge–heat rate of ES-TS;
• $E^{\mathrm{ES}}-E^{\mathrm{TS}}$, energy storage state of ES-TS;
• ${\alpha }^{\mathrm{ES}-\mathrm{ch}}-{\alpha }^{\mathrm{ES}-\mathrm{dch}}$, 0–1 variable in charge–discharge state of ES;
• ${\alpha }^{\mathrm{TS}-\mathrm{ch}}-{\alpha }^{\mathrm{TS}-\mathrm{dch}}$, 0–1 variable in charge–discharge state of TS;
• ${\eta }^{\mathrm{ES}-\mathrm{ch}}-{\eta }^{\mathrm{ES}-\mathrm{dch}}$, charge–discharge efficiency of ES;
• ${\eta }^{\mathrm{TS}-\mathrm{ch}}-{\eta }^{\mathrm{TS}-\mathrm{dch}}$, charge–discharge efficiency of TS;
• ${\mu }^{\mathrm{ES}-\min }-{\mu }^{\mathrm{ES}-\max }$, minimum–maximum energy storage coefficient of ES;
• ${\mu }^{\mathrm{TS}-\min }-{\mu }^{\mathrm{TS}-\max }$, minimum–maximum energy storage coefficient of TS.

5. Solution Method

In order to solve the above model effectively, affine coordinate transformation is first adopted in this paper to transform an uncertain interval optimization problem with parameter correlation into a typical uncertain interval optimization problem. Secondly, the objective function and constraint conditions with interval variables are treated according to the interval order relationship and probability degree, respectively, and converted into a conventional deterministic single-objective optimization problem. On this basis, the genetic algorithm (GA) is further used to solve the problem.

5.1. Correlation Processing Based on Affine Coordinate Transformation

First, the theory of affine coordinate system is adopted to transform the uncertainty domain $\Omega$ of variable U in the original coordinate system into variable P in the affine coordinate system, and the relationship between variables $U_i(i=1,2,\dots ,\mathsf{q})$ and $P_i(i=1,2,\dots ,\mathsf{q})$ can be obtained:

$$\left[\begin{array}{l}{U}_1\\ U_2\\ ⋮\\ U_q\end{array}\right]=A\left[\begin{array}{l}{P}_1\\ P_2\\ ⋮\\ P_q\end{array}\right]+\left[\begin{array}{l}{U}_1^c\\ U_2^c\\ ⋮\\ U_q^c\end{array}\right]$$

(41)

where $U^c$ represents the midpoint vector of the interval, $U_i^c=\left(U_i^L+U_i^R\right)/2$, and matrix A represents the transformation matrix. The elements of matrix A can be obtained through the following formula:

$$a_{ij}=\{\begin{array}{ll}0,& j>i\\ \frac{\cos {\theta }_k-{\displaystyle \sum _{m=1}^{j-1}{a}_{im}{a}_{jm}}}{a_{jj}},& j<i\\ \sqrt{1-{\displaystyle \sum _{l=1}^{j-1}{a}_{il}^2}}& j=i\end{array}$$

(42)

After the above transformation, a general interval optimization problem can be obtained.

$$\{\begin{array}{l}\underset{\mathit{X}}{\min } F(\mathit{X},\mathsf{P})\\ {s}.{t}.H_i(\mathit{X},\mathsf{P})\le {\tilde{b}}_i=[b_i^{\min },b_i^{\max }],i=1,2,\dots ,l,\mathit{X}\in {\Omega }^{\mathrm{n}}\\ \mathrm{P}\in {\mathrm{P}}^I=\left[{{P}}^L,{{P}}^R\right],P_i\in P_i^I=\left[P_i^L,P_i^R\right],i=1,2,\dots ,q\end{array}$$

(43)

5.2. Deterministic Transformation of Interval Optimization Model

In this section, the objective function and the uncertainty of the constraint condition are transformed, respectively.

5.2.1. Objective Function Transformation

For the objective function $f(\mathit{X},\mathsf{P})$, the influence caused by uncertain variable P at decision variable X can be expressed as interval number $\left[\underset{\_}{f}(\mathit{X}),\overline{f}(\mathit{X})\right]$. Where $\underset{\_}{f}(\mathit{X})$ and $\overline{f}(\mathit{X})$, respectively, represent the lower limit and upper limit of the value fluctuation of the objective function through the interval analysis method [26], as shown in Equation (44).

$$\{\begin{array}{l}\underset{\_}{f}(\mathit{X})=\underset{\mathit{U}}{\min }f(\mathit{X},\mathsf{P})\\ \overline{f}(\mathit{X})=\underset{\mathit{U}}{\max }f(\mathit{X},\mathsf{P})\end{array}$$

(44)

In order to quantitatively judge the quality of the target interval and find the optimal decision variable, the interval order relation is used in this paper [26] to process the objective function in Equation (43) and convert it into an equivalent deterministic target consisting of the midpoint $f_{}^{\mathrm{m}}(\mathit{X})$ and radius values $f_{}^{\mathrm{w}}(\mathit{X})$ of the interval:

$$\underset{\mathit{X}}{\min }f(\mathit{X},\mathsf{P})=\underset{\mathit{X}}{\min }\left[\langle f_{}^{\mathsf{m}}(\mathit{X}),f_{}^{\mathsf{w}}(\mathit{X})\rangle \right]$$

(45)

Among them, $f^{\mathsf{m}}( \mathit{X} )=[\underset{\_}{f}( \mathit{X} )+\overline{f}( \mathit{X} )] /2$ and $f^{\mathsf{w}}( \mathit{X} )=[\overline{f}( \mathit{X} )-\underset{\_}{f}( \mathit{X} )] / 2$ respectively, reflect the expected benefits of the planning scheme and its sensitivity to the impact of uncertainties.

In practical application, because decision makers have different preferences for investment risk and return, in order to flexibly meet the above requirements, this paper adopts the linear weighted summing method to integrate $f^{\mathsf{m}}(\mathit{X})$ and $f^{\mathsf{w}}(\mathit{X})$ in Equation (45), and, finally, obtains the optimization objective function in standard form as follows

$$\underset{\mathit{X}}{\min }{\phi }_{}{f}_{}^{\mathrm{m}}(\mathit{X})+(1-{\phi }_{})f_{}^{\mathrm{w}}(\mathit{X})$$

(46)

where $\phi$ is the weight coefficient reflecting decision makers’ expectations of each optimization goal and volatility preference.

5.2.2. Constraint Transformation

For the constraint h_i(X,P) in the model, the possible values caused by uncertain variable P at the decision variable X can be expressed as interval numbers $[\underset{\_}{h_i}(\mathit{X}),\overline{h_i}(\mathit{X})]$. In this paper, the interval probability is used to transform it [28]. Compared with the interval order relation, the interval probability degree uses the quantized probability degree function to judge whether the relevant interval satisfies the given relation constraint. Therefore, its mathematical meaning and objectivity are stronger, so it is more suitable for the transformation of the constraint conditions containing interval numbers.

According to the above methods, the interval constraint condition h_i(X,P) in Equation (42) can be transformed into the following deterministic form:

$$\psi (h_i(\mathit{X})\le {\tilde{b}}_i)\ge {\lambda }_i, i=1,2,\dots ,l$$

(47)

where Ψ(.) is the interval probability representation function. In practical engineering, many methods, such as fuzzy set, grey number and probability, can be used to establish the concrete form of Ψ(.), and the judgment of interval relation is realized accordingly. Without a loss of generality, this paper uses the probabilistic proportional analysis method to calculate the above function values [28].

Here, $h_i(\mathit{X})$ is the possible interval corresponding to the uncertainty constraint i at X, that is $h_i(\mathit{X})=[{\underset{\_}{h}}_i(\mathit{X}),{\overline{h}}_i(\mathit{X})]$, which can be determined by the following formula:

$$\{\begin{array}{l}{\underset{\_}{h}}_i(\mathit{X})=\underset{\mathit{U}}{\min } h_i(\mathit{X},\mathsf{P})\\ {\overline{h}}_i(\mathit{X})=\underset{\mathit{U}}{\max } h_i(\mathit{X},\mathsf{P})\end{array}$$

(48)

In addition, ${\lambda }_i$ in Formula (47) represents the possible degree limit for the constraint i, whose size determines the feasible domain of the optimization variable X in the model. A larger choice ${\lambda }_i$ indicates that the decision maker has a stricter requirement for this constraint condition.

Through the above steps, the single-objective interval optimization model established in this paper is transformed into a conventional deterministic single-objective optimization problem, which can be solved by using the classical single-objective solution algorithm.

5.3. Algorithm Flow

In order to solve the above model efficiently, a solution method combining a genetic algorithm [29] and interval analysis method is adopted in this paper, and it is based on affine coordinate transformation. Its basic process is shown in Figure 3.

The main steps of the above algorithm are described as follows:

(1)

Parameter initialization. Read the input data of the system and set relevant parameters.

(2)

Create the generated population. The initial population (i.e., the initial solution set X) is generated by random functions.

(3)

According to Equations (41)–(43), affine coordinate transformation is performed on the uncertainty domain $\Omega$ of variable U to obtain the uncertainty domain of variable P in the new coordinate system.

(4)

For various group individuals, interval analysis method is used to calculate the upper and lower bounds of each objective function and constraint conditions.

(5)

According to Equations (44)–(47), the midpoint and radius of the objective function, as well as the possibility degree of all constraints, are calculated to achieve the deterministic transformation of the original model.

(6)

Determine individual fitness and determine whether it conforms to the optimization criteria. If so, output the result. If not, proceed to the next step.

(7)

Cross selection produces new individuals, and fitness calculation is carried out to select individuals with high fitness.

(8)

Generate variation to get a new generation of population, and return to Step 6.

6. Results and Discussion

6.1. Parameter Setting

In order to verify the effectiveness of the method proposed in this paper, a typical EH structure is used for simulation analysis [30], as shown in Figure 1. The technical and economic parameters of various hub components are presented [31,32,33,34] in

Table 1. Parameters of EH equipment.

Device Type	Technical Parameters	Economic Parameters
Cogeneration	$y^{CHP}$ = 20 years ${\eta }^e$ = 0.35 ${\beta }^{CHP}$ = 1.3	$c^{CHP-inv}$ = 4500 yuan/kW $c^{CHP-mai}$ = 45 yuan/kW
Wind generator	$y^{WG}$ = 25 years	$c^{WG-inv}$ = 7000 yuan/kW $c^{WG-mai}$ = 90 yuan/kW
Photovoltaic	$y^{PV}$ = 25 years	$c^{PV-inv}$ = 12,000 yuan/kW $c^{PV-mai}$ = 240 yuan/kW
Electric boiler	$y^{EB}$ = 15 years ${\eta }^{EB}$ = 0.95	$c^{EB-inv}$ = 1000 yuan/kW $c^{EB-mai}$ = 40 yuan/kW
Gas boiler	$y^{GB}$ = 15 years ${\eta }^{GB}$ = 0.86	$c^{GB-inv}$ = 800 yuan/kW $c^{GB-mai}$ = 32 yuan/kW
Electrical energy storage	$y^{ES}$ = 10 years ${\eta }^{ES-ch}$ = 0.95 ${\eta }^{ES-dch}$ = 0.95 ${\mu }^{ES-\min }$ = 0.2 ${\mu }^{ES-\max }$ = 0.9 ${\lambda }^{ES}$ = 0.001 ${\xi }^{ES}$ = 0.2	$c^{ES-inv}$ = 1800 yuan/kWh $c^{ES-mai}$ = 18 yuan/kWh
Thermal energy storage	$y^{TS}$ = 20 years ${\eta }^{TS-ch}$ = 0.88 ${\eta }^{TS-dch}$ = 0.88 ${\mu }^{TS-\min }$ = 0.1 ${\mu }^{TS-\max }$ = 0.9 ${\lambda }^{TS}$ = 0.01 ${\xi }^{TS}$ = 0.2	$c^{TS-inv}$ = 200 yuan/kWh $c^{TS-mai}$ = 2 yuan/kWh

. According to the EH load level on the user side, and combined with the simulation results obtained by adjusting the upper limit value of the equipment in the early stage, the maximum configurable capacity of CHP, renewable distributed power supply and electric–thermal boiler were set as 2000 kW, 1000 kW and 1000 kW, respectively. In addition, considering the power capacity characteristics of energy storage devices in the current market, the upper limit of energy storage in EH was set as 400 kW/2000 kWh. In the actual project, the above constraint values can be set flexibly according to the demand of EH investors and specific environmental conditions. In addition, the electricity price EH purchases from external markets is shown in Figure 4, and the gas purchase price is 3.45 yuan/m³ [31]. We showed that the base electricity price of electricity sold without DR is always 0.8 yuan/kWh [32], the selling heat price is 0.7 yuan/kWh [33], carbon emission tax rate is 65.86 yuan/t, and CO₂ emission intensity in the grid is 0.92 kg/kWh [34].

The total power–heat demand curve of the system and the time-division variation of various loads are shown in Figure 5 and

Table 2. Proportion of different types of loads in EH.

Period	${\mathit{\gamma }}_{\mathit{t}}^{\mathbf{UL}-\mathbf{e}}$	${\mathit{\gamma }}_{\mathit{t}}^{\mathbf{TSL}}$	${\mathit{\gamma }}_{\mathit{t}}^{\mathbf{EFL}-\mathbf{e}}$	${\mathit{\gamma }}_{\mathit{t}}^{\mathbf{UL}-\mathbf{h}}$	${\mathit{\gamma }}_{\mathit{t}}^{\mathbf{EFL}-\mathbf{h}}$
22:00–7:00	80%	10%	10%	30%	70%
7:00–8:00 11:00–18:00	20%	60%	20%	50%	50%
8:00–11:00 18:00–22:00	40%	20%	40%	40%	60%

.

The sunrise force prediction curves of wind power and PV are shown in Figure 6 and Figure 7. According to the prediction accuracy of the existing commercial weather software, the prediction errors of wind power and photovoltaic output are set as ±20% and ±15% of their predicted values, respectively. The gray bands in Figure 6 and Figure 7 indicate the fluctuation range of output due to prediction errors. In addition, the price elasticity coefficient of EFL-TSL in the system is shown in

Table 3. Price elasticity of EFL-TSL demands.

Modulus of Elasticity	${\tilde{\mathit{\epsilon }}}_{\mathit{t}}^{\mathbf{TSL}}$	${\tilde{\mathit{\epsilon }}}_{\mathit{t}{\mathit{t}}^{\prime }}^{\mathbf{TSL}}$	${\tilde{\mathit{\epsilon }}}_{\mathit{t}}^{\mathbf{EFL}}$
Band	[−2, −1.45]	[0.06, 0.08]	[−1.6, −1]

, and the fluctuation range is assumed to be ±15%.

6.2. Analysis of Calculation Result

6.2.1. The Influence of Correlation

In order to reveal the impact of correlation on EH operation, this section makes a comparative analysis of EH planning schemes under the two conditions of whether correlation exists or not and a consideration of the size of the correlation coefficient.

Based on the parameter setting in Section 5.1, EH planning and its cost–benefit situation are considered in the case of whether there is correlation among light intensity, wind speed and load uncertainty. The specific calculation results are shown in the

Table 4. Optimization results with and without correlation.

Equipment	The Installed Capacity
	With Correlation	No Correlation
CHP (KW)	1209	1286
EB (KW)	956	974
GB (KW)	972	998
WG (KW)	892	859
PV (KW)	443	395
ES (kWh)	890	902
TS (kWh)	813	832
Investment cost (ten thousand yuan)	197.7	160.2
Operating cost (ten thousand yuan)	474.6	534.1
DR cost (ten thousand yuan)	43.4	45.5
Carbon emission cost (ten thousand yuan)	94.1	98.7
Total cost (ten thousand yuan)	809.8	838.5

.

Table 4 compares the optimization results in the two cases. Compared with the planning scheme without considering the correlation, the installation capacity of WG and PV in the planning scheme, considering the correlation, is greatly increased. Although the investment cost is slightly higher, the total economic cost and environmental benefits have better performance. This means that the effect of correlation does benefit the system.

Specifically, the operating cost of the planning scheme considering the correlation was reduced by 11%, the total cost was reduced by 3.4%, and the carbon emission cost was reduced by 4.7%. This is mainly because, taking into account the correlation, the uncertainty of RES will have a greater impact on the demand response of users, so that WG and PV replace some of the electricity purchased from the grid. This makes full use of RES resources in EH and effectively reduces carbon emissions. At the same time, it avoids purchasing electricity from external market during peak hours, thus reducing the total operating cost of the system and offsetting the impact of increased investment cost.

Figure 8, Figure 9, Figure 10 and Figure 11 shows the results of the operation with or without correlation. Figure 8 and Figure 9 show the operating results of the component without considering the correlation. Figure 10 shows the optimal real-time electricity price of the system when the correlation is considered. Figure 11 shows the change of users’ load demands in EH when the correlation is considered.

As can be seen from the simulation results, when correlation is not considered, WG and PV installation capacity is small due to low matching between demand response and RES output uncertainty. The efficiency of RES is low. The economic benefit of the system is poor. As can be seen from Figure 10, the electricity optimal real-time price changes after considering the correlation. Among them, peak electricity prices mainly appear between 18:00 and 21:00, while trough electricity prices mainly appear between 1:00 and 6:00. Figure 11 shows that the demand side can be adjusted flexibly according to the operating state of the system. When RES output is large, electricity selling price is low, and electricity load demand is large at this time. Therefore, there is consistency between users’ daily load curve and RES output under the consideration of correlation. At the same time, during off-peak hours, relatively low electricity prices lead consumers to shift their heat demand to electric heating and reduce the heat load demand, thus reducing the heat generated by CHP and imported heat from external markets.

We considered the correlation between light intensity, wind speed and the uncertainty of DR, and analyzed the change of each optimization index under the premise of different correlation. In order to compare the influence of different correlation coefficients on the planning results, we set four correlation values.

Table 5. Optimization results under different correlation coefficient.

Equipment	The Installed Capacity
	ρij = 0.2	ρij = 0.4	ρij = 0.6	ρij = 0.8
CHP (KW)	1239	1209	1192	1178
EB (KW)	963	956	949	942
GB (KW)	982	972	969	965
WG (KW)	874	892	912	927
PV (KW)	421	443	459	474
ES (kWh)	894	890	886	882
TS (kWh)	819	813	807	801
Investment cost (ten thousand yuan)	187.4	197.7	204.9	210.6
Operating cost (ten thousand yuan)	505.2	474.6	450.1	431.1
DR cost (ten thousand yuan)	44.6	43.4	40.1	38.1
Carbon emission cost (ten thousand yuan)	96.0	94.1	90.8	88.3
Total cost (ten thousand yuan)	833.2	809.8	785.9	768.1

compares the optimization results in four cases. Obviously, the stronger the correlation, the larger the installation capacity of the WG and PV, the smaller the operating cost and total cost, and the better the economic benefit.

Specifically, when the maximum correlation coefficient is 0.8, the operating cost, total cost and carbon emission cost can be reduced by 19%, 8.4% and 10.5%, respectively. When the correlation coefficient is at least 0.2, the operating cost, total cost and carbon emission cost can be reduced by 5.4%, 0.6% and 2.7%, respectively. It can be seen from Figure 12 that economic benefit and environmental benefit are positively correlated with the correlation. This is mainly because the large correlation coefficient indicates that users’ demand response is more sensitive to RES uncertainty, and the daily load demand of users is more consistent with RES output. Therefore, when RES output is large, TSL and EFL will be promoted to move to this time. It can not only reduce the electricity purchase cost during peak load hours, but also increase the consumption of RES. When the correlation coefficient is large, the system demand for RES power generation increases, which will greatly increase the installation capacity of WG and PV. Therefore, in EH planning, considering the correlation with the installation capacity of RES is helpful in obtaining better expected benefits.

6.2.2. Sensitivity Analysis

In the above study, carbon emission tax is carried out at a fixed price. However, this price may vary considerably with technological innovation, policy reform, and other externalities [35]. In order to elucidate the influence of its change on the system, we have conducted a sensitivity analysis. The optimization results under different carbon emission tax rates (from 66 yuan/t to 264 yuan/t) are shown in Figure 13 and Figure 14. Figure 13 shows the variation of installed capacity of WG and PV under different carbon emission tax rates. With the increase in carbon tax, RES installation capacity increased. However, after rising to twice the initial price, it tended to saturate. This is because higher carbon tax leads EH to favor the use of RES generation to avoid higher carbon tax costs. However, restricted by the consumption of RES and economic benefits, the installation capacity of RES tends to be stable after reaching a certain carbon tax price. At the same time, higher RES output will reduce electricity selling price, thus stimulating the demand response of users and reducing CO₂ emissions. It can be inferred from this curve that the too narrow price range in RTP may lead to inefficient demand response, which is not conducive to a reduction in CO₂ emissions. This curve could also be used by the government for determining the pricing regulation schemes that properly reflect the impacts of prices on emission reductions.

As can be seen from Figure 14, when the carbon emission tax rate increases, the trend of total investment cost is not constant. As RES installation capacity increases, operating costs will decrease, which will offset the impact of higher investment costs. This is because reasonable pricing of environmental costs can help stimulate market enthusiasm for investment in RES power generation. It also increases the willingness of grid operators to use RES, thus increasing the installed capacity of WG and PV. And at the same time, it will force grid optimization to reduce the operating cost of the system. When the carbon tax is higher than 132 yuan/t, the high installation price of PV-WG and the high carbon tax will lead to the increase of the total cost. The carbon tax is 132 yuan/t, corresponding to the valley value of the cost.

7. Conclusions

In this paper, a single objective interval optimization planning method for integrated energy systems is proposed, which considers the correlation between uncertainties. The model in this paper takes the non-direct-controlled DR project based on real-time price as an example, and considers the influence of the correlation of the three uncertainties of light intensity, wind speed and DR so as to realize the optimization of the economic benefits of EH. According to the characteristics of the model: Firstly, the correlation between uncertainties is dealt with by the multidimensional parallelepiped interval model and affine coordinate transformation; then, by using the interval order relation and the possibility degree method, it is transformed into a deterministic single objective optimization problem. Compared with traditional robust optimization focusing on the optimal value in the worst case and stochastic optimization based on scenarios with high probability of occurrence, interval optimization considers both the best case and the worst case. Therefore, interval optimization can flexibly consider the optimality and robustness of the EH plan according to the decision maker’s subjective preference for risk. By flexibly adjusting the weights for optimality and robustness, the planning scheme obtained by interval optimization is generally more economical than robust optimization, and more robust than stochastic optimization. As such, it has good engineering practical value.

Based on the experimental cases, this paper draws the following conclusions:

• The correlation between renewable energy sources and demand response uncertainty will affect the operating results of EH. Due to the strong correlation, the load demand response can better fit the RES output curve. So, compared with the traditional planning method without considering the correlation, the model based on this paper can obtain higher economic benefits, and the economic environmental benefits are positively correlated with the correlation. Therefore, in EH planning, considering the correlation between RES output uncertainty and demand response uncertainty can achieve better returns.
• Carbon emission tax directly affects the optimal planning result of the system. Higher carbon tax helps to promote the installation and output of renewable energy sources and reduce CO₂ emissions. However, when the price of carbon tax is higher, or doubled, the promotion effect tends to flatten out.