Optimal configuration and operation of the regional integrated energy system considering carbon emission and integrated demand response

Under the “carbon peaking and carbon neutrality” development strategy, in order to suppress load fluctuations and promote renewable energy consumptions in the regional integrated energy system involving concentrating solar power stations, a double-layer optimization model based on the improved non-dominated sorting genetic algorithm-II (NSGA-II) and mixed integer linear programming (MILP) is proposed. The upper layer completes the capacity configuration process based on multiple objectives to minimize the annual planning cost and the net emission of pollutants. The lower layer is designed to minimize the annual operating cost and optimize the output of the devices and the load curves through the participation of the integrated demand response process for flexible loads and the whole process of carbon emission including carbon capture, carbon utilization, and carbon trading mechanisms to obtain the best operating plan. The final results indicate that the participation of concentrating solar power stations can improve the level of coordinated optimization, and the improved NSGA-II is stronger than the conventional one in convergence ability. Besides, considering the whole process of carbon emission and integrated demand response is capable of decreasing the annual operating cost and net carbon emission to improve the economy and environmental protection of the system significantly.

1 Introduction

In recent years, environmental problems have become more serious, so the exploration of clean energy will become an inevitable trend in the future development (Fan et al., 2021). The regional integrated energy system (RIES) breaks the barrier between energy planning and operation, and its internal multi-energy coupling equipment can realize energy gradient utilization (Chen et al., 2022; Wang B. et al., 2022), which plays a huge part in realizing the goal of promoting economic and environmental benefits.

In the field of RIES low-carbon operation, relevant studies mostly focus on the regulatory means and economic mechanisms such as carbon capture, carbon storage, and carbon trading. In the work of Dong et al. (2022), a carbon capture model combined with power-to-gas and gas-fired units was installed, which effectively minimized the carbon cost of the total system through the improved energy hub formulation. In the work of Zhang D. et al. (2021), a carbon storage device with pollutant treatments and carbon capture systems (CCSs) was regarded as the structure of an evolutionary IES, which significantly improved renewable energy penetration under different specifications. In the work of Yan et al. (2023), a seasonal-stepped carbon trading mechanism was introduced, and the impacts of economics on optimal dispatch are also considered comprehensively. However, few studies have been conducted on cost-effective carbon utilization processes compared with costly carbon storage processes.

Research on the participation of flexible loads in integrated demand response (IDR) is generally classified into price type, incentive type, and substitution type according to their guiding modes. Yang et al. (2020) divided the IDR into price-based and alternative parts with the process of rolling optimization and finally showed that the aforementioned method can restrain the fluctuation of the loads. Wang et al. (2021) put forward an uncertain model of the demand response through the energy coupling matrices to investigate the impact of price incentives under different scenarios in the RIES to realize the improvement in the load profile. Shao et al. (2019) promoted the IDR to smart buildings and then developed a real-time exchange market with the feasible region method to adjust energy consumption behaviors. Zhang et al. (2022) proposed a multi-objective model considering two-dimensional demand responses with spatiotemporal coupling characteristics to obtain the control strategy among different benefit weights. However, most research on IDR focuses on electricity and heat at present, while cool energy and gas energy are rarely considered.

In addition, with the steady development of concentrating solar power (CSP) in renewable energy, photothermal power generation has gradually attracted wide attention for its advantages of good controllability and high adjustability. Zhao et al. (2019) analyzed the influences of wind power uncertainty on optimal dispatch and put forward a stochastic model for the combination of the CSP stations and wind farms according to the simulation results. In the work of Zhang G. et al. (2021), CSP stations were introduced as cogenerated units and a multi-dispatch method for the RIES was proposed. The results show that the participation of CSP stations can reduce operating costs. On the premise of considering operating cost, Jiang et al. (2020) built an exchange model between CSP stations and energy markets to participate in DR programs, which significantly improved the energy operating efficiency through the price elasticity matrix of the electrical and the heat loads. However, the existing literature has generally ignored the potential of CSP stations operating in conjunction with the aforementioned CCS.

In view of the problems mentioned previously, the main contributions of this paper are shown as follows:

1) An energy hub (EH) with the participation of the CSP station and the whole process of carbon emission including carbon capture, carbon utilization, and carbon trading is established by considering the power-to-gas equipment.

2) The analysis includes flexible loads such as electricity, heat, cooling, and gas and successively elaborates them for the uncontrollable loads, transferable loads, curtailable loads, and fungible loads, which are refined to reflect the IDR.

3) A double-layer optimization model of improved NSGA-II and MILP is constructed. The upper layer takes the annual planning cost and net pollutant discharge as the target for device selection and capacity configuration through the improved NSGA-II, while the lower layer regards the annual operating cost as the target to optimize the output of each device through the Cplex solver.

This paper is organized as follows. In Section 2, the basic structure and the operating principle of the RIES are introduced. Section 3 is focused on the expression of the proposed CSP station, the whole process of carbon emission, and the IDR. In Section 4, a double-layer model considering the solving methods is developed to realize the optimized process. Case studies are conducted in Section 5, in which the proposed model is simulated under different scenarios. At last, the conclusions are given in Section 6.

2 Basic structure of the RIES

This paper focuses on the RIES shown in Figure 1. The system includes the wind turbine (WT), photovoltaic (PV), CSP station, electric boiler (EB), electric chiller (EC), waste heat boiler (WHB), absorption chiller (AC), gas turbine (GT), gas boiler (GB), carbon capture system (CCS), power to gas (P2G), and storage battery (SB). The energy input sources of the RIES are electricity and natural gas, and the loads include electricity, heat, cooling, and gas. The CSP station can be regarded as cogenerated equipment. Similarly, GT consumes natural gas to generate heat, which can be recovered by WHB. For the CCS, the CO₂ captured is mainly from coal-fired plants in the power grid and gas-fired units (GT and GB) in the RIES. The power consumption of P2G can be supplied by abandoned wind and light, thus realizing the absorptive process of renewable energy. The entered electricity and natural gas of the EH can be purchased from power grid companies and natural gas companies, respectively, and then reasonably distributed to various energy conversion equipment and user-side loads.

FIGURE 1

FIGURE 1. Structure diagram of the RIES.

3 The whole process of carbon emission and integrated demand response in the RIES

3.1 Conventional device model of the RIES

The conventional devices of the RIES mainly include gas-fired units, EB, EC, AC, and WHB. The models are shown as follows:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}}={\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{e}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{h}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{e}}={\mathit{\eta }}_{\mathbf{G}\mathbf{T},\mathbf{e}}{\mathit{H}}_{\mathbf{g}}{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{h}}={\mathit{\eta }}_{\mathbf{G}\mathbf{T},\mathbf{h}}{\mathit{H}}_{\mathbf{g}}{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{B},\mathbf{h}}={\mathit{\eta }}_{\mathbf{G}\mathbf{B},\mathbf{h}}{\mathit{H}}_{\mathbf{g}}{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{B}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{h}}={\mathit{\eta }}_{\mathbf{E}\mathbf{B},\mathbf{h}}{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{B}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{C},\mathbf{c}}={\mathit{\eta }}_{\mathbf{E}\mathbf{C},\mathbf{c}}{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{C}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{A}\mathbf{C},\mathbf{c}}={\mathit{\eta }}_{\mathbf{A}\mathbf{C},\mathbf{c}}{\mathit{P}}_{\mathit{t}}^{\mathbf{h},\mathbf{A}\mathbf{C}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{H}\mathbf{B},\mathbf{h}}={\mathit{\eta }}_{\mathbf{W}\mathbf{H}\mathbf{B},\mathbf{h}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{h}},\end{array}\end{array}\end{array}\end{array}\right.(1)$

where $P_t^{\mathrm{G}\mathrm{T}}$, $P_t^{\mathrm{G}\mathrm{T},\mathrm{e}}$, and $P_t^{\mathrm{G}\mathrm{T},\mathrm{h}}$ are the total output, electrical output, and thermal output of GT at time $t$, respectively; $P_t^{\mathrm{G}\mathrm{B},\mathrm{h}}$, $P_t^{\mathrm{E}\mathrm{B},\mathrm{h}}$, $P_t^{\mathrm{E}\mathrm{C},\mathrm{c}}$, $P_t^{\mathrm{A}\mathrm{C},\mathrm{c}}$, and $P_t^{\mathrm{W}\mathrm{H}\mathrm{B},\mathrm{h}}$ are the output of GB, EB, EC, AC, and WHB at time $t$, respectively; $Q_t^{\mathrm{G}\mathrm{T}}$ and $Q_t^{\mathrm{G}\mathrm{B}}$ are the amount of natural gas consumed by GT and GB at time $t$, respectively; $P_t^{\mathrm{e},\mathrm{E}\mathrm{B}}$, $P_t^{\mathrm{e},\mathrm{E}\mathrm{C}}$, and $P_t^{\mathrm{h},\mathrm{A}\mathrm{C}}$ are the input power of EB, EC, and AC at time $t$, respectively; ${\eta }_{\mathrm{G}\mathrm{T},\mathrm{e}}$ and ${\eta }_{\mathrm{G}\mathrm{T},\mathrm{h}}$ are the electrical efficiency and thermal efficiency of GT, respectively; ${\eta }_{\mathrm{G}\mathrm{B},\mathrm{h}}$, ${\eta }_{\mathrm{E}\mathrm{B},\mathrm{h}}$, ${\eta }_{\mathrm{E}\mathrm{C},\mathrm{c}}$, ${\eta }_{\mathrm{A}\mathrm{C},\mathrm{c}}$, and ${\eta }_{\mathrm{W}\mathrm{H}\mathrm{B},\mathrm{h}}$ are the corresponding energy conversion efficiency of GB, EB, EC, AC, and WHB, respectively; and $H_{\mathrm{g}}$ is the calorific value of the natural gas.

3.2 CSP thermoelectric conversion model

As an emerging form of power generation in recent years, the CSP station is mainly divided into tower type, trough type, disk type, and linear Fresnel type, among which the tower type has been widely used in engineering practice for its advantages of strong economy and good performance (Gorman et al., 2021). In this paper, the heat storage tank and EB are aggregated to model the internal and external energy transfer relationships of the tower-type CSP station. The typical structure is shown in Figure 2.

FIGURE 2

FIGURE 2. EB–CSP framework.

The heat energy collected by the heat collector in the optical field can be stored in the heat storage tank through the heat transfer mediums, and it can also be used to generate electricity through the heat cycle. The expression of photothermal conversion of the heat collector is shown as follows:

${\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{h}}={\mathit{\eta }}_{\mathbf{s}-\mathbf{h}}{\mathit{S}}_{\mathit{S}\mathbf{F}}{\mathit{D}}_{\mathit{t}},(2)$

where $P_t^{\mathrm{S}\mathrm{F},\mathrm{h}}$ is the thermal power obtained by the heat collector at time $t$; ${\eta }_{\mathrm{s}-\mathrm{h}}$ is the photothermal conversion efficiency; $S_{\mathrm{S}\mathrm{F}}$ is the area of the heliostat field; and $D_t$ is the direct solar radiation index (DNI) at time $t$.

The heat storage tank can store the converted heat energy and can also release the heat energy to meet the generation demand or directly supply the heat to the load side. The model of the heat storage tank is shown in the following formula:

${\mathit{S}}_{\mathit{t}}^{\mathbf{H}\mathbf{Q}}={\mathit{S}}_{\mathit{t}-\mathbf{1}}^{\mathbf{H}\mathbf{Q}}\left(\mathbf{1}-{\mathit{\rho }}_{\mathbf{h}}\right)+{\mathit{\eta }}_{\mathbf{c}\mathbf{h}\mathbf{a}}^{\mathbf{H}\mathbf{Q}}{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}-\frac{{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}}{{\mathit{\eta }}_{\mathbf{d}\mathbf{i}\mathbf{s}}^{\mathbf{H}\mathbf{Q}}},(3)$

where $S_t^{\mathrm{H}\mathrm{Q}}$ and $S_{t-1}^{\mathrm{H}\mathrm{Q}}$ are the thermal storage of the heat storage tank at time $t$ and time $t-1$, respectively; ${\rho }_{\mathrm{h}}$ is the energy self-loss coefficient; ${\eta }_{\mathrm{c}\mathrm{h}\mathrm{a}}^{\mathrm{H}\mathrm{Q}}$ and ${\eta }_{\mathrm{d}\mathrm{i}\mathrm{s}}^{\mathrm{H}\mathrm{Q}}$ are the charging and releasing efficiency of the heat storage tank, respectively; and $P_t^{\mathrm{T}\mathrm{S},\mathrm{c}}$ and $P_t^{\mathrm{T}\mathrm{S},\mathrm{f}}$ are the charging and releasing power of the heat storage tank at time $t$, respectively.

To keep the CSP station in stable operation after the participation of EB, its internal heat cycle must meet the following relations:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{h}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{i}\mathbf{n}}={\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{r}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{e}\mathbf{f}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{H}\mathbf{L}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{a}\mathbf{b}\mathbf{n},\mathbf{h}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{i}\mathbf{n}}={\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{h}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{l}\mathbf{o}\mathbf{a}\mathbf{d}},\end{array}\right.(4)$

where $P_t^{\mathrm{S}\mathrm{F},\mathrm{r}}$ and $P_t^{\mathrm{T}\mathrm{S},\mathrm{e}\mathrm{f}}$ are the thermal power directly generated by the heat collector and the heat storage tank at time $t$, respectively; $P_t^{\mathrm{E}\mathrm{B},\mathrm{i}\mathrm{n}}$ and $P_t^{\mathrm{E}\mathrm{B},\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}}$ are the thermal power provided by EB to the heat storage tank and the heat load side, respectively; and $P_t^{\mathrm{H}\mathrm{L}}$ and $P_t^{\mathrm{a}\mathrm{b}\mathrm{n},\mathrm{h}}$ are the thermal power provided to the heat load side and heat energy loss at time $t$, respectively.

The generation power of the CSP station mainly comes from the heating power of the heat collector and the heat storage tank:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}}={\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{d}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{d}\mathbf{f}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{d}}={\mathit{\eta }}_{\mathbf{d}}{\mathit{P}}_{\mathit{t}}^{\mathbf{S}\mathbf{F},\mathbf{r}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{d}\mathbf{f}}=\left(\mathbf{1}-{\mathit{\eta }}_{\mathbf{f}}\right){\mathit{\eta }}_{\mathbf{d}}{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{e}\mathbf{f}},\end{array}\right.(5)$

where $P_t^{\mathrm{S}\mathrm{F},\mathrm{d}}$ and $P_t^{\mathrm{T}\mathrm{S},\mathrm{d}\mathrm{f}}$ are the heat of the collector and the storage tank that provides for power generation, respectively; ${\eta }_{\mathrm{d}}$ and ${\eta }_{\mathrm{f}}$ are the thermoelectric conversion efficiency of the heat collector and released loss rate of the heat storage tank, respectively.

The heat energy provided to the load side is expressed as follows:

${\mathit{P}}_{\mathit{t}}^{\mathbf{H}\mathbf{L}}={\mathit{\eta }}_{\mathbf{h}}\left({\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{e}\mathbf{f}}\right),(6)$

where ${\eta }_{\mathrm{h}}$ is the heat transfer efficiency of the heat storage tank.

3.3 The whole process of carbon emission including P2G collaborative operation

3.3.1 Operating characteristics of P2G

In this paper, P2G can absorb the power of abandoned wind and light, which is used to generate natural gas. The energy consumption of P2G is shown in the following equation:

${\mathit{P}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G}}={\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{A}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{A}},(7)$

where $P_t^{\mathrm{W}\mathrm{A}}$ and $P_t^{\mathrm{V}\mathrm{A}}$ are the power of abandoned wind and light at time $t$, respectively.

The amount of CO₂ consumed in P2G can be formulated as follows:

${\mathit{Q}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}},\mathbf{P}\mathbf{2}\mathbf{G}}={\mathit{\alpha }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}{\mathit{\eta }}_{\mathbf{P}\mathbf{2}\mathbf{G}}{\mathit{P}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G}},(8)$

where ${\alpha }_{\mathrm{C}{\mathrm{O}}_2}$ is the CO₂ consumption per unit of power; ${\eta }_{\mathrm{P}2\mathrm{G}}$ is the conversion efficiency of P2G.

The natural gas produced by P2G can be calculated as follows:

${\mathit{Q}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G},\mathbf{C}{\mathbf{H}}_{\mathbf{4}}}=\frac{\mathbf{3.6}{\mathit{\eta }}_{\mathbf{P}\mathbf{2}\mathbf{G}}{\mathit{P}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G}}}{{\mathit{H}}_{\mathbf{g}}}.(9)$

3.3.2 Principles of carbon capture and carbon utilization

The CCS mainly includes the absorption tower, regeneration tower, compressor, and other structural units. The absorber uses a specific solution to absorb CO₂ from the flue gas and transfers it to the regenerator, where it is heated and separated. Then, CO₂ is compressed in the compressor for transporting. Therefore, the energy consumption of carbon capture and gas treatment generated by the three links mentioned previously are the main sources of the total energy consumption in the CCS (Yan et al., 2017), whose expression is shown as follows:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}\mathbf{S}}={\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{A}\mathbf{P}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{A}\mathbf{P}}={\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}\mathbf{S},\mathbf{r}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}\mathbf{S},\mathbf{f}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}={\mathit{\lambda }}_{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}\left({\mathit{Q}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}}+{\mathit{Q}}_{\mathit{t}}^{\mathbf{S}\mathbf{T}}\right),\end{array}\right.(10)$

where $P_t^{\mathrm{C}\mathrm{A}\mathrm{P}}$ and $P_t^{\mathrm{D}\mathrm{E}\mathrm{A}\mathrm{L}}$ are carbon capture energy consumption and gas treatment energy consumption of the CCS, respectively; $P_t^{\mathrm{C}\mathrm{C}\mathrm{S},\mathrm{r}}$ and $P_t^{\mathrm{C}\mathrm{C}\mathrm{S},\mathrm{f}}$ are the operating energy consumption and fixed energy consumption of carbon capture, respectively; ${\lambda }_{\mathrm{D}\mathrm{E}\mathrm{A}\mathrm{L}}$ is the unit energy consumption coefficient of flue gas treatment; and $Q_t^{\mathrm{C}\mathrm{S}}$ and $Q_t^{\mathrm{S}\mathrm{T}}$ are the flue gas treatment provided by carbon source units and the flue gas storage tank at time $t$, respectively. The operating energy consumption and fixed energy consumption of carbon capture are, respectively, satisfied.

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}\mathbf{S},\mathbf{r}}={\mathit{\lambda }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}{\mathit{E}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}\mathbf{S},\mathbf{f}}=\mathbf{0.1}{\mathit{\lambda }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}{\mathit{E}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}},\end{array}\right.(11)$

where ${\lambda }_{\mathrm{C}{\mathrm{O}}_2}$ is the operating energy consumption coefficient of unit CO₂ treated by the CCS; $E_t^{\mathrm{C}{\mathrm{O}}_2}$ is CO₂ that is captured.

The capacity of the flue gas storage tank is shown as follows:

${\mathit{Q}}_{\mathit{t}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}={\mathit{Q}}_{\mathit{t}-\mathbf{1}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}+{\mathit{Q}}_{\mathit{t}}^{\mathbf{I}\mathbf{N}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{S}\mathbf{T}},(12)$

where $Q_t^{DEAL}$ and $Q_{t-1}^{DEAL}$ are the capacities of the flue gas storage tank at time $t$ and $t-1$, respectively.

In this paper, a combined operation strategy of wind power–photovoltaic–CSP–carbon capture based on the participation of the CCS and new energy units is proposed. The output of new energy units is partly used for carbon capture, partly used for flue gas treatment, and the rest is transported to the power grid. The energy consumption process of the carbon capture is presented in Figure 3 in the following section, and the energy consumption process of the flue gas treatment is similarly omitted.

FIGURE 3

FIGURE 3. Flow chart of carbon capture energy.

The expression of joint operation is shown as follows:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{A}\mathbf{P}}={\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{C}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{C}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}={\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{D}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{W}}={\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{C}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{A}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{V}}={\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{C}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{A}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}}={\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{C}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{D}},\end{array}\end{array}\end{array}\right.(13)$

where $P_t^{\mathrm{W}}$ and $P_t^{\mathrm{V}}$ are the predicted output power of WT and PV at time $t$, respectively; $P_t^{\mathrm{W}\mathrm{C}}$, $P_t^{\mathrm{V}\mathrm{C}}$, and $P_t^{\mathrm{C}\mathrm{C}}$ are the energy consumption of carbon capture provided by WT, PV, and CSP station at time $t$, respectively; $P_t^{\mathrm{W}\mathrm{D}}$, $P_t^{\mathrm{V}\mathrm{D}}$, and $P_t^{\mathrm{C}\mathrm{D}}$ are the energy consumption of flue gas treatment provided by WT, PV, and CSP station at time $t$, respectively; and $P_t^{\mathrm{W}\mathrm{N}}$, $P_t^{\mathrm{V}\mathrm{N}}$, and $P_t^{\mathrm{C}\mathrm{N}}$ are the net output of power generation (on-grid power) provided by WT, PV, and CSP station at time $t$, respectively.

Carbon utilization refers to sending captured CO₂ into P2G to participate in the synthesis of CH₄ as its raw material. This process can decrease the carbon storage cost and increase the operating flexibility of P2G. The carbon utilization process satisfies the following relations:

${\mathit{Q}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}},\mathbf{P}\mathbf{2}\mathbf{G}}=\left(\mathbf{1}-{\mathit{\eta }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right){\mathit{E}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{F}\mathbf{C}},(14)$

where ${\eta }_{\mathrm{C}{\mathrm{O}}_2}$ is the heat loss rate of CO₂; $Q_t^{\mathrm{F}\mathrm{C}}$ is the stored CO₂ at time $t$.

3.3.3 Carbon trading mechanism

The carbon trading mechanism regards the carbon emission as a commodity and controls it through the trading of carbon emission rights between producers and markets. If the actual carbon emission is lower than the allocated, the surplus quotas can be sold to carbon trading markets. Otherwise, the corresponding quotas need to be purchased additionally (Wang X. et al., 2022).

1) Quota models of carbon emission

Carbon emission quota is the amount of emission allowance allocated by the regulatory authorities to each carbon source within the RIES, which varies according to the type of equipment (Chen et al., 2021). In this paper, there are two main carbon sources, namely, coal-fired power plants in the power grid (superior purchasing power) and gas-fired units of the system (GT and GB). Then, the quota models of carbon emission can be expressed as follows:

$\left\{\begin{array}{c}{\mathit{E}}_{\mathbf{q}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}={\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}+{\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{T}}+{\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{B}},\\ {\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}={\mathit{\mu }}_{\mathbf{e}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}},\\ \begin{array}{c}{\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{T}}={\mathit{\mu }}_{\mathbf{g}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}},\\ {\mathit{E}}_{\mathbf{q}}^{\mathbf{G}\mathbf{B}}={\mathit{\mu }}_{\mathbf{g}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{B},\mathbf{h}},\end{array}\end{array}\right.(15)$

where $E_{\mathrm{q}}^{\mathrm{G}\mathrm{R}\mathrm{I}\mathrm{D}}$, $E_{\mathrm{q}}^{\mathrm{G}\mathrm{T}}$, and $E_{\mathrm{q}}^{\mathrm{G}\mathrm{B}}$ are carbon emission quotas of power purchase, GT, and GB, respectively; ${\mu }_{\mathrm{e}}$ and ${\mu }_{\mathrm{g}}$ are the baseline credits for carbon emission per unit of power consumption and per unit of gas consumption for coal-fired and gas-fired units, respectively; $P_t^{\mathrm{G}\mathrm{R}\mathrm{I}\mathrm{D}}$ is the purchased power at time $t$; and $T$ is the operating cycle, which values as 24h.

2) Practical models of carbon emission

Since the CCS can absorb CO₂, the actual model of carbon emission after considering it can be expressed as follows:

$\left\{\begin{array}{c}{\mathit{E}}_{\mathbf{a}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}={\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}+{\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{T}}+{\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{B}}-{\mathit{E}}_{\mathbf{a}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}},\\ {\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}={\mathit{\delta }}_{\mathbf{e}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}},\\ \begin{array}{c}{\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{T}}={\mathit{\delta }}_{\mathbf{g}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}},\\ {\mathit{E}}_{\mathbf{a}}^{\mathbf{G}\mathbf{B}}={\mathit{\delta }}_{\mathbf{g}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{B},\mathbf{h}},\\ {\mathit{E}}_{\mathit{a}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{E}}_{\mathit{t}}^{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}},\end{array}\end{array}\right.(16)$

where $E_{\mathrm{a}}^{\mathrm{G}\mathrm{R}\mathrm{I}\mathrm{D}}$, $E_{\mathrm{a}}^{\mathrm{G}\mathrm{T}}$, and $E_{\mathrm{a}}^{\mathrm{G}\mathrm{B}}$ are the practical carbon emission of power purchase, GT, and GB, respectively; $E_{\mathrm{a}}^{\mathrm{C}{\mathrm{O}}_2}$ is the total amount of CO₂ that is captured; and ${\delta }_{\mathrm{e}}$ and ${\delta }_{\mathrm{g}}$ are the carbon emission intensity of coal-fired units and gas-fired units, respectively.

3) Ladder-type carbon trading

After getting the quotas of carbon emission and the practical model through the process previously, the transaction volume of carbon emission rights that participate in trading markets can be formulated as follows:

${\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}={\mathit{E}}_{\mathbf{a}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}-{\mathit{E}}_{\mathbf{q}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}.(17)$

Compared with the traditional carbon trading mechanism, the ladder-type carbon trading mechanism has more strict controls over carbon emissions. The principle is dividing carbon emissions into multiple intervals through the stepped pricing method. With the increase in carbon emissions, the transaction cost within the corresponding interval will increase (Zhang et al., 2016; Fu et al., 2022; Li et al., 2022). Accordingly, the cost of ladder-type carbon trading can be expressed as follows:

${\mathit{F}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}=\left\{\begin{array}{c}{\mathit{B}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\le \mathit{l},\\ {\mathit{B}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\left[{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\left(\mathbf{1}+{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right)-\mathit{l}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right]\mathit{l}\le {\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\le \mathbf{2}\mathit{l},\\ {\mathit{B}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\left[{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\left(\mathbf{1}+\mathbf{2}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right)-\mathbf{3}\mathit{l}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right]\mathbf{2}\mathit{l}\le {\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\le \mathbf{3}\mathit{l}\\ {\mathit{B}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\left[{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\left(\mathbf{1}+\mathbf{3}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right)-\mathbf{6}\mathit{l}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right]\mathbf{3}\mathit{l}\le {\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\le \mathbf{4}\mathit{l},\\ {\mathit{B}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\left[{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\left(\mathbf{1}+\mathbf{4}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right)-\mathbf{10}\mathit{l}{\mathit{\gamma }}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\right]{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\ge \mathbf{4}\mathit{l},\end{array}\right.,(18)$

where $B_{\mathrm{C}{\mathrm{O}}_2}$ is the benchmark price of the carbon trading; ${\gamma }_{\mathrm{C}{\mathrm{O}}_2}$ is the growth rate of the carbon tax price; and $l$ is the interval length of carbon emission.

3.4 Integrated demand response of flexible loads considering compensation cost

IDR refers to the process in which users reasonably adjust energy-use modes and participate in energy interaction according to different energy prices or incentive mechanisms to optimize the load curves (Shao et al., 2021). Therefore, it is necessary to classify the loads first in the analysis of such problems. In addition to uncontrollable loads (invariable loads) in previous studies, the loads are divided into variable loads such as transferable loads, curtailable loads, and fungible loads in this paper. In addition, compared to conventional studies that only consider electrical and heat loads, this paper incorporates all types of loads involving electricity, heat, cooling, and gas into the IDR. Then, according to the characteristics of the variable loads mentioned previously, the process of IDR involved in this paper is able to divide it into the price type and the substitutable type.

3.4.1 Price-based demand response

Since different types of loads have variant sensitivities to the same price signal, this paper divides the loads involved in the price-based demand response into curtailable loads (CLs) and transferable loads (TLs) and then builds the models of them in turn.

1) Properties and modeling for CL

CL uses the price-demand elasticity matrix to describe its characteristics. In allusion to the t line and the j column element $e_{t,j}$ in the price-based elastic matrix $E\left(t,j\right)$, the elastic coefficient of the load at time t to energy price at time j is defined as follows:

${\mathit{e}}_{\mathit{t},\mathit{j}}=\frac{\Delta {\mathit{L}}_{\mathit{L},\mathit{t}}^{\mathit{i}}/{\mathit{L}}_{\mathit{L},\mathit{t}}^{\mathit{i}}}{\Delta {\mathit{\kappa }}_{\mathit{j}}/{\mathit{\kappa }}_{\mathit{j}}^{\mathbf{0}}},(19)$

where $\Delta L_{L,t}^i$ and $L_{L,t}^i$ are the variable amount and initial amount of the load for class $i$ participating in the price-based demand response at time $t$, respectively; $\Delta {\kappa }_j$ and ${\kappa }_j^0$ are the variable amounts of energy prices and initial energy prices at time $j$ after the IDR, respectively.

Then, the CL variation of class $i$ at time $t$ after the IDR can be presented as follows:

$\Delta {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathit{i}}={\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathit{i}\mathbf{0}}\left[{\displaystyle \sum _{\mathit{j}=\mathbf{1}}^{\mathit{T}}}{\mathit{E}}_{\mathbf{C}\mathbf{L}}^{\mathit{i}}\left(\mathit{t},\mathit{j}\right)\frac{{\mathit{\kappa }}_{\mathit{j}}-{\mathit{\kappa }}_{\mathit{j}}^{\mathbf{0}}}{{\mathit{\kappa }}_{\mathit{j}}^{\mathbf{0}}}\right],(20)$

where $L_{\mathrm{C}\mathrm{L},t}^{i0}$ is the amount of initial CL for class $i$ at time $t$; $E_{\mathrm{C}\mathrm{L}}^i\left(t,j\right)$ is the price-demand elastic matrix of CL for class $i$; and ${\kappa }_j$ is the corresponding load energy price of class $i$ at time $j$.

2) Properties and modeling for TL

TL can realize flexible controls of working hours and power at different periods under the total load unchanged premise according to the energy prices of users’ own demand response. Taking time-sharing energy prices and incentive measures as signals, users can be guided to transfer the peak loads to the normal or trough period (Azzam et al., 2023). Similarly, after expressing the characteristics of IDR in the price-demand elasticity matrix, the TL variation in class $i$ at time $t$ can be formulated as follows:

$\Delta {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathit{i}}={\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathit{i}\mathbf{0}}\left[{\displaystyle \sum _{\mathit{j}=\mathbf{1}}^{\mathit{T}}}{\mathit{E}}_{\mathbf{T}\mathbf{L}}^{\mathit{i}}\left(\mathit{t},\mathit{j}\right)\frac{{\mathit{\kappa }}_{\mathit{j}}-{\mathit{\kappa }}_{\mathit{j}}^{\mathbf{0}}}{{\mathit{\kappa }}_{\mathit{j}}^{\mathbf{0}}}\right],(21)$

where $L_{\mathrm{T}\mathrm{L},t}^{i0}$ is the amount of initial TL for class $i$ at time $t$; $E_{\mathrm{T}\mathrm{L}}^i\left(t,j\right)$ is the price-demand elastic matrix of TL for class $i$.

3.4.2 Substitution-type demand response

For a certain type of the heat load that can be directly supplied by heat or electricity (electricity to heat), the electrical energy can be consumed in periods of low electricity prices, while the heat energy can be directly consumed in periods of high electricity prices to meet different needs so as to realize the mutual substitution between the electrical and heat energy. The fungible load (FL) model can be expressed as follows:

$\left\{\begin{array}{c}\Delta {\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{e}}={\mathit{\epsilon }}_{\mathbf{e},\mathbf{h}}\cdot \Delta {\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{h}},\\ {\mathit{\epsilon }}_{\mathbf{e},\mathbf{h}}=\frac{{\mathit{\nu }}_{\mathbf{e}}{\mathit{\varphi }}_{\mathbf{e}}}{{\mathit{\nu }}_{\mathbf{h}}{\mathit{\varphi }}_{\mathbf{h}}},\end{array}\right.(22)$

where $\Delta L_{\mathrm{F}\mathrm{L},t}^{\mathrm{e}}$ and $\Delta L_{\mathrm{F}\mathrm{L},t}^{\mathrm{h}}$ are the fungible electricity load and heat load, respectively; ${\epsilon }_{\mathrm{e},\mathrm{h}}$ is the electric-heating substitution coefficient, which varies with the time; ${\nu }_{\mathrm{e}}$ and ${\nu }_{\mathrm{h}}$ are the unit calorific values of electrical energy and heat energy, respectively; and ${\varphi }_{\mathrm{e}}$ and ${\varphi }_{\mathrm{h}}$ are the energy utilization rates of electricity and heat, respectively.

Through the IDR process mentioned previously, the following conditions of load balance can be obtained as follows:

$\left\{\begin{array}{c}{\mathit{L}}_{\mathbf{e}}^{\mathbf{D}\mathbf{R}}={\mathit{L}}_{\mathbf{e}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{e}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{e}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{e}},\\ {\mathit{L}}_{\mathbf{h}}^{\mathbf{D}\mathbf{R}}={\mathit{L}}_{\mathbf{h}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{h}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{h}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{h}},\\ \begin{array}{c}{\mathit{L}}_{\mathbf{c}}^{\mathbf{D}\mathbf{R}}={\mathit{L}}_{\mathbf{c}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{c}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{c}},\\ {\mathit{L}}_{\mathbf{g}}^{\mathbf{D}\mathbf{R}}={\mathit{L}}_{\mathbf{g}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{g}}+\mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{g}},\end{array}\end{array}\right.(23)$

where $L_{\mathrm{e}}$, $L_{\mathrm{h}}$, and $L_{\mathrm{c}}$ are the initial loads of electricity, heat, and cooling, respectively; $L_{\mathrm{g}}$ is the amount of natural gas required by the initial gas load; $L_{\mathrm{e}}^{\mathrm{D}\mathrm{R}}$, $L_{\mathrm{h}}^{\mathrm{D}\mathrm{R}}$, and $L_{\mathrm{c}}^{\mathrm{D}\mathrm{R}}$ are the loads of electricity, heat, and cooling after participating in the IDR, respectively; and $L_{\mathrm{g}}^{\mathrm{D}\mathrm{R}}$ is the amount of natural gas required by the gas load after participating in the IDR. Then, the load compensation cost can be calculated as follows:

$\left\{\begin{array}{c}{\mathit{F}}_{\mathbf{c}\mathbf{u}\mathbf{t}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}\left({\mathit{c}}_{\mathbf{c}\mathbf{u}\mathbf{t}}^{\mathbf{e}}\cdot \Delta {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{e}}+{\mathit{c}}_{\mathbf{c}\mathbf{u}\mathbf{t}}^{\mathbf{h}}\cdot \Delta {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{h}}+{\mathit{c}}_{\mathbf{c}\mathbf{u}\mathbf{t}}^{\mathbf{c}}\cdot \Delta {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{c}}+{\mathit{c}}_{\mathbf{c}\mathbf{u}\mathbf{t}}^{\mathbf{g}}\cdot \Delta {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathbf{g}}\right),\\ {\mathit{F}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}\left({\mathit{c}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}^{\mathbf{e}}\cdot \Delta {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{e}}+{\mathit{c}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}^{\mathbf{h}}\cdot \Delta {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{h}}+{\mathit{c}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}^{\mathbf{c}}\cdot \Delta {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{c}}+{\mathit{c}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}^{\mathbf{g}}\cdot \Delta {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathbf{g}}\right),\\ {\mathit{F}}_{\mathbf{s}\mathbf{u}\mathbf{b}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}\left({\mathit{c}}_{\mathbf{s}\mathbf{u}\mathbf{b}}^{\mathbf{e}}\cdot \Delta {\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{e}}+{\mathit{c}}_{\mathbf{s}\mathbf{u}\mathbf{b}}^{\mathbf{h}}\cdot \Delta {\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathbf{h}}\right),\end{array}\right.(24)$

where $F_{\mathrm{c}\mathrm{u}\mathrm{t}}$, $F_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$, and $F_{\mathrm{s}\mathrm{u}\mathrm{b}}$ are the compensation costs that can be reduced, transferred, and replaced, respectively; $c_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\mathrm{e}}$, $c_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\mathrm{h}}$, $c_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\mathrm{c}}$, and $c_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\mathrm{g}}$ are the compensation prices of electrical, heat, cool, and gas loads that can be reduced per unit of power, respectively; $c_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{e}}$, $c_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{h}}$, $c_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{c}}$, and $c_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}^{\mathrm{g}}$ are the compensation prices of electrical, heat, cool, and gas loads that can be transferred per unit of power, respectively; and $c_{\mathrm{s}\mathrm{u}\mathrm{b}}^{\mathrm{e}}$ and $c_{\mathrm{s}\mathrm{u}\mathrm{b}}^{\mathrm{h}}$ are the compensation prices of electrical and heat load that can be replaced per unit of power, respectively.

4 Double-layer optimization model of the RIES

Based on the RIES with the participation of CSP stations as shown in Figure 1, this paper establishes a double-layer optimization model that considers the whole process of carbon emission and the IDR to demonstrate the innovations. The upper layer randomly generates the planned capacities of the devices and transmits them to the lower layer. The lower layer constrains the output of the devices and transmits the results back to the upper layer. Finally, the upper layer revises the capacity configuration of each device again. The iterative process between the upper layer and the lower layer leads to the most optimal configuration and the lowest annual operating cost of the RIES.

4.1 Models of upper programming

The upper model aims at minimizing the annual planning cost and the annual net pollutant emission of the RIES, and the decision variables are the installed capacities of different devices. The mathematical formulas are shown as follows:

$\left\{\begin{array}{c}\mathbf{min}{\mathit{F}}_{\mathbf{p}\mathbf{c}}={\mathit{F}}_{\mathbf{i}\mathbf{n}\mathbf{v}}+{\mathit{F}}_{\mathbf{r}\mathbf{c}}\\ \mathbf{min}{\mathit{E}}_{\mathbf{e}\mathbf{n}\mathbf{v}}=\mathbf{365}{\displaystyle \sum }_{\mathit{s}\in \left\{\mathit{s}\mathit{u}\mathit{m},\mathit{w}\mathit{i}\mathit{n},\mathit{t}\mathit{r}\mathit{a}\right\}}{\mathit{\theta }}_{\mathit{s}}\left({\mathit{E}}_{\mathbf{r}}^{\mathbf{S}{\mathbf{O}}_{\mathbf{2}}}+{\mathit{E}}_{\mathbf{r}}^{\mathbf{N}{\mathbf{O}}_{\mathbf{x}}}+{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}\right)\end{array}\right.(25)$

where $F_{\mathrm{p}\mathrm{c}}$ and $E_{\mathrm{e}\mathrm{n}\mathrm{v}}$ are the annual planning cost and the annual net emission of various pollutants during the operating life of the RIES, respectively; $F_{\mathrm{i}\mathrm{n}\mathrm{v}}$ and $F_{\mathrm{r}\mathrm{c}}$ are the annual investment cost and the annual operating cost, respectively; and $E_{\mathrm{r}}^{\mathrm{S}{\mathrm{O}}_2}$ and $E_{\mathrm{r}}^{\mathrm{N}{\mathrm{O}}_{\mathrm{x}}}$ are the net emissions of SO₂ and NO_x, respectively.

The annual operating cost $F_{\mathrm{r}\mathrm{c}}$ is introduced separately as the objective function of the lower model, while other variables can be expressed as follows:

$\left\{\begin{array}{c}{\mathit{F}}_{\mathbf{i}\mathbf{n}\mathbf{v}}={\displaystyle \sum }_{\mathit{k}}{\mathit{c}}_{\mathit{k}}{\mathit{S}}_{\mathbf{max}}^{\mathit{k}}\\ {\mathit{E}}_{\mathbf{r}}^{\mathbf{S}{\mathbf{O}}_{\mathbf{2}}}={\mathit{\xi }}_{\mathbf{e}}^{\mathbf{S}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}\\ {\mathit{E}}_{\mathbf{r}}^{\mathbf{N}{\mathbf{O}}_{\mathbf{x}}}={\mathit{\xi }}_{\mathbf{e}}^{\mathbf{N}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}\end{array}\right.(26)$

where $k$ indicates different types of the devices; $c_k$ is the installed investment cost of per unit capacity for device $k$; $S_{\max }^k$ is the installed capacity of device $k$; and ${\xi }_{\mathrm{e}}^{\mathrm{S}}$ and ${\xi }_{\mathrm{e}}^{\mathrm{N}}$ are the SO₂ and the NO_x emission intensity of coal-fired units in the power grid, respectively.

4.2 Models of lower operating

The lower layer intends to minimize the annual operating cost of the RIES, and the decision variable is the output of each device. The mathematical model is expressed in the following equation:

$\left\{\begin{array}{c}\mathbf{min}{\mathit{F}}_{\mathbf{r}\mathbf{c}}=\mathbf{365}{\displaystyle \sum _{\mathit{m}=\mathbf{1}}^{\mathit{N}}}\frac{\mathit{m}}{{\left(\mathbf{1}+\mathbf{r}\right)}^{\mathit{m}}}\bullet {\displaystyle \sum }_{\mathit{s}\in \left\{\mathit{s}\mathit{u}\mathit{m},\mathit{w}\mathit{i}\mathit{n},\mathit{t}\mathit{r}\mathit{a}\right\}}{\mathit{\theta }}_{\mathit{s}}\\ \left({\mathit{F}}_{\mathbf{e}}+{\mathit{F}}_{\mathbf{g}\mathbf{a}\mathbf{s}}+{\mathit{F}}_{\mathbf{o}\mathbf{c}}+{\mathit{F}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}+{\mathit{F}}_{\mathbf{d}\mathbf{r}}+{\mathit{F}}_{\mathbf{f}\mathbf{c}}\right),\\ {\mathit{F}}_{\mathbf{e}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{C}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}},\\ \begin{array}{c}{\mathit{F}}_{\mathbf{g}\mathbf{a}\mathbf{s}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{C}}_{\mathit{t}}^{\mathbf{G}\mathbf{A}\mathbf{S}}{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{A}\mathbf{S}},\\ {\mathit{F}}_{\mathbf{o}\mathbf{c}}={\mathit{c}}_{\mathbf{D}\mathbf{P}}^{\mathit{k}}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{P}}_{\mathit{t}}^{\mathit{k}},\\ \begin{array}{c}{\mathit{F}}_{\mathbf{d}\mathbf{r}}={\mathit{F}}_{\mathbf{c}\mathbf{u}\mathbf{t}}+{\mathit{F}}_{\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{s}}+{\mathit{F}}_{\mathbf{s}\mathbf{u}\mathbf{b}},\\ {\mathit{F}}_{\mathbf{f}\mathbf{c}}={\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}{\mathit{k}}_{\mathbf{f}\mathbf{c}}{\mathit{Q}}_{\mathit{t}}^{\mathbf{F}\mathbf{C}}.\end{array}\end{array}\end{array}\right.(27)$

This paper adopts the ladder-type carbon trading method, and the corresponding cost $F_{{\mathrm{C}\mathrm{O}}_2}$ is shown in Eq. 18.

Here, $N$ is the planning year; $\mathrm{r}$ is the discounted rate; and ${\theta }_s$ is the proportion of season $s$ in the year, and three typical seasons involving summer (sum), winter (win), and transitional season (tra) are considered in this paper; $F_{\mathrm{e}}$, $F_{\mathrm{g}\mathrm{a}\mathrm{s}}$, $F_{\mathrm{o}\mathrm{c}}$, $F_{\mathrm{d}\mathrm{r}}$, and $F_{\mathrm{f}\mathrm{c}}$ are the costs of electrical purchase, gas purchase, operation or maintenance for each device, IDR compensation, and carbon storage in the RIES, respectively; $C_t^{\mathrm{G}\mathrm{R}\mathrm{I}\mathrm{D}}$ and $C_t^{\mathrm{G}\mathrm{A}\mathrm{S}}$ are the unit price of electrical purchase and gas purchase at time $t$, respectively; $Q_t^{\mathrm{G}\mathrm{A}\mathrm{S}}$ is the flow of gas purchase at time $t$; $c_{\mathrm{D}\mathrm{P}}^k$ is the unit maintenance cost of device $k$; $P_t^k$ is the output of device $k$ at time $t$; and $k_{\mathrm{f}\mathrm{c}}$ is the unit carbon storage price.

4.3 Constraint conditions

4.3.1 Constraints of the CSP station

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}+{\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}\le \mathbf{1},\\ {\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}\cdot \mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}\cdot {\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}},\\ \begin{array}{c}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}\cdot \mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}\cdot {\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{c}}\cdot {\mathit{P}}_{\mathit{t}}^{\mathbf{T}\mathbf{S},\mathbf{f}}=\mathbf{0},\\ \begin{array}{c}\mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}}\le {\mathit{S}}_{\mathit{t}}^{\mathbf{H}\mathbf{Q}}\le {\mathit{S}}_{\mathbf{max}}^{\mathbf{H}\mathbf{Q}},\\ {\mathit{S}}_{\mathbf{0}}^{\mathbf{H}\mathbf{Q}}={\mathit{S}}_{\mathbf{24}}^{\mathbf{H}\mathbf{Q}},\\ \begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}}\le {\mathit{P}}_{\mathbf{max}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}},\\ \left|{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}}-{\mathit{P}}_{\mathit{t}-\mathbf{1}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}}\right|\le \mathbf{\Delta }{\mathit{P}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}},\end{array}\end{array}\end{array}\end{array}\right.(28)$

where ${\mu }_t^{\mathrm{T}\mathrm{S},\mathrm{c}}$ and ${\mu }_t^{\mathrm{T}\mathrm{S},\mathrm{f}}$ are variables of $0-1$, which represent the charging and releasing state of the heat storage tank at time $t$, respectively; $S_{\max }^{\mathrm{H}\mathrm{Q}}$ is the capacity of the heat storage tank; $S_0^{\mathrm{H}\mathrm{Q}}$ and $S_{24}^{\mathrm{H}\mathrm{Q}}$ are the starting and ending values of heat storage during the day, respectively; and $P_{\max }^{\mathrm{C}\mathrm{S}\mathrm{P},\mathrm{e}}$ and $\Delta P^{\mathrm{C}\mathrm{S}\mathrm{P},\mathrm{e}}$ are the upper power limit and the climbing rate of the CSP station, respectively.

4.3.2 Constraints of the CCS

$\left\{\begin{array}{c}{\mathit{E}}_{\mathbf{min}}^{{\mathbf{C}\mathbf{O}}_{\mathbf{2}}}\le {\mathit{E}}_{\mathit{t}}^{{\mathbf{C}\mathbf{O}}_{\mathbf{2}}}\le {\mathit{E}}_{\mathbf{max}}^{{\mathbf{C}\mathbf{O}}_{\mathbf{2}}},\\ {\mathit{Q}}_{\mathbf{min}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}\le {\mathit{Q}}_{\mathit{t}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}}\le {\mathit{Q}}_{\mathbf{max}}^{\mathbf{D}\mathbf{E}\mathbf{A}\mathbf{L}},\\ \begin{array}{c}\mathbf{0}\le {\mathit{Q}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}}\le {\mathit{Q}}_{\mathbf{max}}^{\mathbf{C}\mathbf{S}},\\ \mathbf{0}\le {\mathit{Q}}_{\mathit{t}}^{\mathbf{S}\mathbf{T}}\le {\mathit{Q}}_{\mathbf{max}}^{\mathbf{S}\mathbf{T}},\\ \mathbf{0}\le {\mathit{Q}}_{\mathit{t}}^{\mathbf{I}\mathbf{N}}\le {\mathit{Q}}_{\mathbf{max}}^{\mathbf{I}\mathbf{N}},\end{array}\end{array}\right.(29)$

where $E_{\max }^{{\mathrm{C}\mathrm{O}}_2}$ and $E_{\min }^{{\mathrm{C}\mathrm{O}}_2}$ are the upper and lower limits of CO₂ captured by the CCS, respectively; $Q_{\max }^{\mathrm{D}\mathrm{E}\mathrm{A}\mathrm{L}}$ and $Q_{\min }^{\mathrm{D}\mathrm{E}\mathrm{A}\mathrm{L}}$ are the upper and lower limits of the storage in the flue gas storage tank, respectively; $Q_{\max }^{\mathrm{C}\mathrm{S}}$ and $Q_{\max }^{\mathrm{S}\mathrm{T}}$ are the upper limits of the flue gas treatment provided by the carbon source unit and flue gas storage tank, respectively; and $Q_{\max }^{\mathrm{I}\mathrm{N}}$ is the upper limit of the flue gas that flows into the flue gas storage tank.

4.3.3 Constraints of P2G

$\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G}}\le {\mathit{P}}_{\mathbf{max}}^{\mathbf{P}\mathbf{2}\mathbf{G}},(30)$

where $P_{\max }^{\mathrm{P}2\mathrm{G}}$ is the maximum energy consumption of P2G.

4.3.4 Constraints of new energy units and purchasing energy

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{N}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{W}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{N}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{V}},\\ \begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{N}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{S}\mathbf{P},\mathbf{e}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}\le {\mathit{P}}_{\mathbf{max}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}},\\ \mathbf{0}\le {\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{A}\mathbf{S}}\le {\mathit{Q}}_{\mathbf{max}}^{\mathbf{G}\mathbf{A}\mathbf{S}},\end{array}\end{array}\right.(31)$

where $P_{\max }^{\mathrm{G}\mathrm{R}\mathrm{I}\mathrm{D}}$ and $Q_{\max }^{\mathrm{G}\mathrm{A}\mathrm{S}}$ are the maximum purchased power and gas, respectively.

4.3.5 Constraints of electrical energy storage devices

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}+{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le \mathbf{1},\\ {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\cdot \mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{e}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\cdot {\mathit{S}}_{\mathbf{max}}^{\mathbf{e}},\\ \begin{array}{c}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\cdot \mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{e}}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\cdot {\mathit{S}}_{\mathbf{max}}^{\mathbf{e}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\cdot {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}=\mathbf{0},\\ \begin{array}{c}{\mathit{S}}_{\mathit{t}}^{\mathbf{e}}={\mathit{S}}_{\mathit{t}-\mathbf{1}}^{\mathbf{e}}+\left({\mathit{\eta }}_{\mathbf{e}}^{\mathbf{c}\mathbf{h}\mathbf{a}}{\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}-\frac{{\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}}{{\mathit{\eta }}_{\mathbf{e}}^{\mathbf{d}\mathbf{i}\mathbf{s}}}\right)\Delta \mathit{t},\\ \mathbf{0.1}{\mathit{S}}_{\mathbf{max}}^{\mathbf{e}}\le {\mathit{S}}_{\mathit{t}}^{\mathbf{e}}\le {\mathit{S}}_{\mathbf{max}}^{\mathbf{e}},\\ {\mathit{S}}_{\mathbf{0}}^{\mathbf{e}}={\mathit{S}}_{\mathbf{24}}^{\mathbf{e}},\end{array}\end{array}\end{array}\right.(32)$

where ${\mu }_t^{\mathrm{c}\mathrm{h}\mathrm{a}}$ and ${\mu }_t^{\mathrm{d}\mathrm{i}\mathrm{s}}$ are variables of $0-1$, which represent the charging and releasing state of the storage battery at time $t$, respectively; $S_{\max }^{\mathrm{e}}$ is the capacity of the storage battery; $P_t^{\mathrm{c}\mathrm{h}\mathrm{a}}$ and $P_t^{\mathrm{d}\mathrm{i}\mathrm{s}}$ refer to the charging and releasing power of the storage battery, respectively; $S_t^{\mathrm{e}}$ and $S_{t-1}^{\mathrm{e}}$ are the charge capacities of the storage battery at time $t$ and time $t-1$, respectively; ${\eta }_{\mathrm{e}}^{\mathrm{c}\mathrm{h}\mathrm{a}}$ and ${\eta }_{\mathrm{e}}^{\mathrm{d}\mathrm{i}\mathrm{s}}$ are the charging and releasing efficiency of the storage battery, respectively; $\Delta t$ is the unit operating period, which values as 1; and $S_0^{\mathrm{e}}$ and $S_{24}^{\mathrm{e}}$ are the starting and ending values of the storage battery during the day, respectively.

4.3.6 Constraints of other devices

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathit{k}}\le {\mathit{S}}_{\mathbf{max}}^{\mathit{k}},\\ \left|{\mathit{P}}_{\mathit{t}}^{\mathit{k}}-{\mathit{P}}_{\mathit{t}-\mathbf{1}}^{\mathit{k}}\right|\le \mathbf{\Delta }{\mathit{P}}^{\mathit{k}},\end{array}\right.(33)$

where $S_{\max }^k$ and $\Delta P^k$ are the capacity and climbing rate of device $k$, respectively.

4.3.7 Constraints of integrated demand response

1) Curtailable loads:

${\mathit{L}}_{\mathbf{C}\mathbf{L},\mathbf{min}}^{\mathit{i}}\le \mathbf{\Delta }{\mathit{L}}_{\mathbf{C}\mathbf{L},\mathit{t}}^{\mathit{i}}\le {\mathit{L}}_{\mathbf{C}\mathbf{L},\mathbf{max}}^{\mathit{i}},(34)$

where $L_{\mathrm{C}\mathrm{L},\max }^i$ and $L_{\mathrm{C}\mathrm{L},\min }^i$ are the upper and lower limits of the actual reduction for CL $i$ at time $t$, respectively.

2) Transferable loads:

$\left\{\begin{array}{c}{\displaystyle \sum _{\mathit{t}=\mathbf{1}}^{\mathit{T}}}\mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathit{i}}=\mathbf{0},\\ {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathbf{min}}^{\mathit{i}}\le \mathbf{\Delta }{\mathit{L}}_{\mathbf{T}\mathbf{L},\mathit{t}}^{\mathit{i}}\le {\mathit{L}}_{\mathbf{T}\mathbf{L},\mathbf{max}}^{\mathit{i}},\end{array}\right.(35)$

where $L_{\mathrm{T}\mathrm{L},\max }^i$ and $L_{\mathrm{T}\mathrm{L},\min }^i$ are the upper and lower limits of the actual transfer for TL $i$ at time $t$, respectively.

3) Fungible loads:

$\mathbf{0}\le \mathbf{\Delta }{\mathit{L}}_{\mathbf{F}\mathbf{L},\mathit{t}}^{\mathit{i}}\le {\mathit{L}}_{\mathbf{F}\mathbf{L},\mathbf{max}}^{\mathit{i}},(36)$

where $L_{\mathrm{F}\mathrm{L},\max }^i$ is the maximum substitution power of FL $i$ at time $t$.

4.3.8 Balance of the loads

1) Load balancing before participating in IDR:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{e}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{B}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{C}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}={\mathit{L}}_{\mathbf{e}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{l}\mathbf{o}\mathbf{a}\mathbf{d}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{B},\mathbf{h}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{H}\mathbf{B},\mathbf{h}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{H}\mathbf{L}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{h},\mathbf{A}\mathbf{C}}={\mathit{L}}_{\mathbf{h}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{C},\mathbf{c}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{A}\mathbf{C},\mathbf{c}}={\mathit{L}}_{\mathbf{c}},\\ {\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{A}\mathbf{S}}+{\mathit{Q}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G},\mathbf{C}{\mathbf{H}}_{\mathbf{4}}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{B}}={\mathit{L}}_{\mathbf{g}}.\end{array}\end{array}\right.(37)$

2) Load balancing after participating in IDR:

$\left\{\begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{R}\mathbf{I}\mathbf{D}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{V}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{C}\mathbf{N}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{T},\mathbf{e}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{B}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{e},\mathbf{E}\mathbf{C}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}={\mathit{L}}_{\mathbf{e}}^{\mathbf{D}\mathbf{R}},\\ {\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{B},\mathbf{l}\mathbf{o}\mathbf{a}\mathbf{d}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{G}\mathbf{B},\mathbf{h}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{W}\mathbf{H}\mathbf{B},\mathbf{h}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{H}\mathbf{L}}-{\mathit{P}}_{\mathit{t}}^{\mathbf{h},\mathbf{A}\mathbf{C}}={\mathit{L}}_{\mathbf{h}}^{\mathbf{D}\mathbf{R}},\\ \begin{array}{c}{\mathit{P}}_{\mathit{t}}^{\mathbf{E}\mathbf{C},\mathbf{c}}+{\mathit{P}}_{\mathit{t}}^{\mathbf{A}\mathbf{C},\mathbf{c}}={\mathit{L}}_{\mathbf{c}}^{\mathbf{D}\mathbf{R}},\\ {\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{A}\mathbf{S}}+{\mathit{Q}}_{\mathit{t}}^{\mathbf{P}\mathbf{2}\mathbf{G},\mathbf{C}{\mathbf{H}}_{\mathbf{4}}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{T}}-{\mathit{Q}}_{\mathit{t}}^{\mathbf{G}\mathbf{B}}={\mathit{L}}_{\mathbf{g}}^{\mathbf{D}\mathbf{R}}.\end{array}\end{array}\right.(38)$

4.4 Model solving

For the double-layer optimization model constructed previously, the improved NSGA-II is employed for the upper layer. As for the lower layer, nonlinear problems are transformed into linear problems according to the MILP method, and then, the Cplex solver is called for the corresponding calculation.

4.4.1 Multi-objective problems of the upper layer solved by the improved NSGA-II

This paper uses the NSGA-II for capacity configuration with different research objects and proposes an improved NSGA-II to compare with the conventional NSGA-II. The NSGA-II treats each sub-objective in a high-dimensional multi-objective optimization problem equally without introducing weights, which can avoid the influence of local optimal solutions on capacity configuration. The conventional NSGA-II is mainly composed of selection, crossover, mutation, and non-dominated sorting, among which the mutation process is usually realized by simulated binary mutation operators (SBMO) that leads to low population diversity and search efficiency. Therefore, this paper perfected the mutation process of the conventional algorithm by using adaptive mixed mutation operators. The main principle is to mix the simulated binary mutation operator and the normal distributed mutation operator (NDMO) in an adaptive way to determine the proportion of the two operators in different periods within the algorithm. The model is shown as follows:

$\left\{\begin{array}{c}{\mathit{z}}_{\mathbf{1},\mathit{j}}=\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}-{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\left[\left(\mathbf{1}+\mathit{\beta }\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\mathit{\beta }\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right]+\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\\ \left[\left(\mathbf{1}+\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right],\mathit{\lambda }\le \mathbf{0.5},\\ \begin{array}{c}{\mathit{z}}_{\mathbf{1},\mathit{j}}=\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}-{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\left[\left(\mathbf{1}+\mathit{\beta }\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\mathit{\beta }\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right]+\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\\ \left[\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right],\mathit{\lambda }>\mathbf{0.5},\\ \begin{array}{c}{\mathit{z}}_{\mathbf{2},\mathit{j}}=\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}-{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\left[\left(\mathbf{1}-\mathit{\beta }\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}+\mathit{\beta }\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right]+\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\\ \left[\left(\mathbf{1}+\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right],\mathit{\lambda }\le \mathbf{0.5},\\ \begin{array}{c}{\mathit{z}}_{\mathbf{2},\mathit{j}}=\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}-{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\left[\left(\mathbf{1}-\mathit{\beta }\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}+\mathit{\beta }\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right]+\frac{{\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}}}{\mathbf{2}\cdot {\mathit{k}}_{\mathbf{G}\mathbf{e}\mathbf{n}\mathbf{s}}}\\ \left[\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{1},\mathit{j}}+\left(\mathbf{1}-\left|\mathit{N}\left(\mathbf{0,1}\right)\right|\right){\mathit{p}}_{\mathbf{2},\mathit{j}}\right],\mathit{\lambda }>\mathbf{0.5},\end{array}\end{array}\end{array}\end{array}\right.(39)$

where $z_{1,j}$ and $z_{2,j}$ are the percentages of SBMO and NDMO, respectively; $k_{\mathrm{G}\mathrm{e}\mathrm{n}}$ and $k_{\mathrm{G}\mathrm{e}\mathrm{n}\mathrm{s}}$ are the number of the current population iterations and the maximum population iterations, respectively; $\beta$ and $\left|N\left(\mathrm{0,1}\right)\right|$ are random variables of simulated binary and normal distribution, respectively; $p_{1,j}$ and $p_{2,j}$ are the mutation probabilities of SBMO and NDMO, respectively; and $\lambda$ is the period of mutation.

As can be known from the formula mentioned previously, in the early stages of the algorithm, the proportion of SBMO should be higher to expand the search limit. In addition, in the later stages of the algorithm, the proportion of NDMO should be higher to promote search accuracy. By matching the search range and the search accuracy, the accuracy of the calculations can be significantly improved.

4.4.2 Transformation and solution of the lower nonlinear model

There are some nonlinear terms in the lower layer that needed to be transformed. For instance, the nonlinear terms in the ladder-type carbon trading model are solved by introducing segmenting points and auxiliary variables. As for nonlinear terms within the constraints, the Big-M method is used to deal with them, and then, the original nonlinear constraints are equitably transformed into mixed integer linear constraints by introducing several 0–1 variables. After the aforementioned process, the nonlinear problems of the lower layer are transformed into linear problems, which can be solved by the Cplex solver. The specific MILP process is shown as follows:

1) Big-M method for nonlinear constraint problems

The method is introduced by taking two subproblems of the storage battery cannot charge or release at the same moment and the charging or releasing power constraint as examples. The conversion methods of other nonlinear constraints are similar and will not be described here.

The general mathematical expression that the storage battery cannot charge or release at the same moment can be written as follows:

$\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\perp {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\ge \mathbf{0},(40)$

where “$0\le a\perp b\ge 0$”means $a\ge 0,b\ge 0$ and $ab=0$.

By using the Big-M method and introducing $0-1$ variables, the aforementioned equation can be equivalently transformed into

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\le \mathbf{M}\cdot {\mathit{\mu }}_{\mathit{t}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le \mathbf{M}\cdot \left(\mathbf{1}-{\mathit{\mu }}_{\mathit{t}}\right),\end{array}\right.(41)$

where ${\mu }_t$ is the binary variable of $0-1$; $\mathrm{M}$ is a large constant.

The general expression of charging or releasing power constraint can be presented as follows:

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\le {\mathit{P}}_{\mathbf{max}}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le {\mathit{P}}_{\mathbf{max}}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}},\\ \begin{array}{c}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}+{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le \mathbf{1},\\ {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\in \left\{\mathbf{0,1}\right\},{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\in \left\{\mathbf{0,1}\right\},\end{array}\end{array}\right.(42)$

where $P_{\max }$ is the maximum charging and releasing power of the storage battery.

According to the Big-M method, the aforementioned equation can be equivalently transformed into

$\left\{\begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\le {\mathit{P}}_{\mathbf{max}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\mathbf{M},\\ \begin{array}{c}\mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le {\mathit{P}}_{\mathbf{max}},\\ \mathbf{0}\le {\mathit{P}}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le {\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\mathbf{M},\\ \begin{array}{c}{\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}+{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\le \mathbf{1},\\ {\mathit{\mu }}_{\mathit{t}}^{\mathbf{c}\mathbf{h}\mathbf{a}}\in \left\{\mathbf{0,1}\right\},{\mathit{\mu }}_{\mathit{t}}^{\mathbf{d}\mathbf{i}\mathbf{s}}\in \left\{\mathbf{0,1}\right\}.\end{array}\end{array}\end{array}\right.(43)$

2) Linearization of the ladder-type carbon trading

The model of ladder-type carbon trading is detailed in Eq. 18 previously, and the concrete implementation of its linearization is shown as follows.

The aforementioned formula is a piecewise function of five sections, so six piecewise points $w_1,w_2,\cdots ,w_6$, six continuous auxiliary variables $u_1,u_2,\cdots ,u_6$, and five binary auxiliary variables $v_1,v_2,\cdots ,v_5$ are added to satisfy the following expression:

$\left\{\begin{array}{c}{\mathit{u}}_{\mathbf{1}}+{\mathit{u}}_{\mathbf{2}}+\cdots +{\mathit{u}}_{\mathbf{6}}=\mathbf{1},\\ {\mathit{v}}_{\mathbf{1}}+{\mathit{v}}_{\mathbf{2}}+\cdots +{\mathit{v}}_{\mathbf{5}}=\mathbf{1},\\ \begin{array}{c}{\mathit{u}}_{\mathbf{1}}\ge \mathbf{0},{\mathit{u}}_{\mathbf{2}}\ge \mathbf{0},\cdots ,{\mathit{u}}_{\mathbf{6}}\ge \mathbf{0},\\ {\mathit{u}}_{\mathbf{1}}\le {\mathit{v}}_{\mathbf{1}},{\mathit{u}}_{\mathbf{2}}\le {\mathit{v}}_{\mathbf{1}}+{\mathit{v}}_{\mathbf{2}},{\mathit{u}}_{\mathbf{3}}\le {\mathit{v}}_{\mathbf{2}}+{\mathit{v}}_{\mathbf{3}},\\ {\mathit{u}}_{\mathbf{4}}\le {\mathit{v}}_{\mathbf{3}}+{\mathit{v}}_{\mathbf{4}},{\mathit{u}}_{\mathbf{5}}\le {\mathit{v}}_{\mathbf{4}}+{\mathit{v}}_{\mathbf{5}},{\mathit{u}}_{\mathbf{6}}\le {\mathit{v}}_{\mathbf{5}}.\end{array}\end{array}\right.(44)$

Then, the aforementioned equations can be transformed into the following linear expression:

$\left\{\begin{array}{c}{\mathit{E}}_{\mathbf{r}}^{\mathbf{R}\mathbf{I}\mathbf{E}\mathbf{S}}={\displaystyle \sum _{\mathit{n}=\mathbf{1}}^{\mathbf{6}}}{\mathit{u}}_{\mathit{n}}{\mathit{w}}_{\mathit{n}},\\ {\mathit{F}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}={\displaystyle \sum _{\mathit{n}=\mathbf{1}}^{\mathbf{6}}}{\mathit{u}}_{\mathit{n}}{\mathit{F}}_{\mathbf{C}{\mathbf{O}}_{\mathbf{2}}}\left({\mathit{w}}_{\mathit{n}}\right).\end{array}\right.(45)$

Taking all the previous factors into consideration, the problems of MILP that are covered in this paper can be finally converted as follows:

$\left\{\begin{array}{c}\mathbf{min}\mathit{c}\mathit{x}\\ \mathbf{s}.\mathbf{t}.\mathit{A}\mathit{x}\left\{\begin{array}{c}\ge \\ =\\ \le \end{array}\right\}=\mathit{b},\\ \begin{array}{c}{\mathit{x}}_{\mathbf{min}}\le {\mathit{x}}_{\mathit{p}}\le {\mathit{x}}_{\mathbf{max}}\mathit{p}\in \mathbf{I},\\ {\mathit{x}}_{\mathit{q}}\in \left\{\mathbf{0,1}\right\}\mathit{q}\in \mathbf{J},\end{array}\end{array}\right.(46)$

where $cx$ is the objective function, which stands for the annual operating cost of the RIES; $A$ and $b$ are the coefficient matrix and its corresponding value of subfunctions, respectively; $x_p$ is the continuous variable, which represents the upper and lower limits of constraint conditions; and $x_q$ is the integer variable, which values as 0 or 1.

4.4.3 Solving procedure

Combined with the double-layer model, the solving procedure in this paper can be described as follows:

1) We input the efficiency and the unit cost of each device, as well as the related load data, forecast output power of WT and PV, time-of-use electrical price and gas price, etc.

2) We initialize the improved NSGA-II. The basic parameters of the algorithm in the upper programming model should be reasonably set. In this paper, the population number is 50, and the maximum number of iterations is 500. Then initial values are assigned to the decision variables and the number of iterations to generate a random initial population and start the iteration process.

3) Based on the two optimization objectives of each population, the fitness (objective function value) can be calculated, respectively, and the corresponding capacity configuration of the device is substituted into the lower operating model as the constraint condition of running for each device.

4) The MILP method transforms the nonlinear model into a linear one in the lower layer, and the Cplex solver is called to calculate for it. Then, the lower layer returns the calculation results of the whole carbon emission process and the output power of the devices under IDR to the upper layer as the constraints.

5) According to the constraints returned by the lower model, the upper model completes the non-dominated sorting process and calculates the crowding degree of the objective functions.

6) The adaptive mixed mutation operators are used to complete the mutation process of the NSGA-II and combine it with the crossover process to merge the populations. Then, the non-dominated sorting process and the calculation of the crowding degree are performed again to decrease the error.

7) The elite retention strategy is introduced to generate new populations and update their positions through the tournament selection mechanism after the competitive process is completed.

8) We determine whether the termination condition is just satisfied. If the current number of iterations does not reach the maximum one, we skip back to step (3) to continue iterating to continuously improve the accuracy of the calculation. Otherwise, the configuration capacity and the planned operation scheme of each device can be output.

The specific solving process is presented in Figure 4.

FIGURE 4

FIGURE 4. Flow chart of solving the double-layer optimization model based on the improved NSGA-II and MILP method.

5 Example analysis

5.1 Comparative analysis under different cases

For further study impacts of the whole carbon emission process and the IDR on the operation and capacity configuration in the system, this paper makes some improvements on the example shown in the work of Wei et al. (2022) and Zeng et al. (2023), and the following four cases are built for comparative analysis:

Case I: The whole process of carbon emission and the IDR involving flexible loads are not considered.

Case II: The whole process of carbon emission is not considered, and the IDR involving flexible loads is considered.

Case III: The whole process of carbon emission is considered, and the IDR involving flexible loads is not considered.

Case IV: The whole process of carbon emission and the IDR involving flexible loads are both considered.

According to the four cases constructed previously, the outcome is enumerated in Table 1.

TABLE 1

TABLE 1. Cost results of comparative cases within the year (unit: ten thousand Yuan/year).

As shown in Table 1, compared with case I, the costs of carbon trading and annual operating in case III are reduced by 81.61% and 2.84%, respectively. This is because case III considers the whole process of carbon emission and sends the captured CO₂ into P2G for carbon utilization, which effectively reduces the net emission of CO₂. Moreover, introducing the ladder-type carbon trading mechanism makes the initial carbon quotas offset part of the carbon trading cost, thus reducing the annual operating cost. In addition, compared with case I, the cost of purchasing electricity and gas in case II is reduced by 73.43% and 16.47%, respectively. This is because the process of IDR can significantly cut down the peak-time flexible loads and scale up the valley-time flexible loads, especially the change in the electrical load, which allows the energy purchasing method to be more economical and selective. Compared with case III, the cost of carbon storage in case IV is lower. This is because the electrical load is decreased during peak periods after considering the IDR, which leads to the output power reduction of new energy units. In addition, the reduced part is converted to abandoned wind and light, which increases the amount of CO₂ consumed by P2G. Compared with case II, the compensatory cost of IDR in case IV is slightly lower. This is because the peak-to-valley differences of the electrical load are reduced by increasing the on-grid power of new energy units under the combined operation demand, which decreases the compensatory cost of IDR. In addition, the annual operating cost, electricity and gas purchasing cost, device operation and maintenance cost, and carbon trading cost in case IV are all smaller than those in case III, because considering the IDR in the whole process of carbon emission can transfer parts of flexible loads in high-price periods to low-price periods and lower their energy consumption. Moreover, the response process of the mutual substitution involving electricity, heat, and other energy sources also significantly reduces the energy purchasing cost, thus making the operating modes more reasonable, and the economy and the environmental protection tend to be more coordinated in the system.

In addition, compared to case IV, cases I, II, and III have the lowest cost for their respective parts, but the other costs are at a higher level and the overall performance is poor. Case IV has generally lower costs for each subcomponent and has the best overall performance with annual operating cost improvement rates of 3.25%, 2.37%, and 0.42%.

Meanwhile, in an effort to validate the availability and the superiority of the algorithm that is proposed in this paper, quantitative analysis is added on the basis of qualitative analysis to compare the improved NSGA-II with the conventional one. Assuming that the number of populations is 50, the maximum number of iterations is 500, and running each standard test function 30 times, the relevant parameters of test functions (Tao et al., 2019) and the final results including the standard deviation (SD), average fitness (AF), optimal value (OV), worst value (WV), and convergence curve can be obtained as shown in Table 2 and Figure 5.

TABLE 2

TABLE 2. Relevant parameters of test functions and the calculation results.

FIGURE 5

FIGURE 5. Average convergence curve of the test. (A) Results of the function Sphere, and (B) results of the function Rastrigrin.

As is shown in Table 2, both the conventional algorithm and improved NSGA-II can complete the search process, but the improved NSGA-II exhibits higher search accuracy and more divergent data, which indicates the better population diversity. In addition, as shown in Figure 5, the improved NSGA-II is superior to the conventional one in terms of convergence speed, convergence accuracy, global search ability, and optimized stability.

Table 3 shows the configuration of each device obtained by the two algorithms in case IV. As shown in Table 3, capacity configuration results under the improved NSGA-II are generally lower than those under the conventional one, especially for new energy units and storage batteries. This shows that the improved NSGA-II is able to reduce the annual investment cost effectively and then further decrease the annual planning cost of the system, which can diminish approximately 4.03 million Yuan. In addition, the net emission of pollutants is also cut down so that environmental protection can be improved. Main parameters of various types of the devices are listed in Table 4 below.

TABLE 3

TABLE 3. Configuration results of different algorithms.

TABLE 4

TABLE 4. Main parameters of various types of the devices.

5.2 Analysis of optimized comparison

This paper takes a typical day in winter to analyze the operating results of the entire RIES. The operating period is 24 h, and the unit operating period is 1 h. The output forecast of renewable energy on a typical winter day is shown in Figure 6. In addition, the corresponding time-of-use price curves are shown in Figure 7.

FIGURE 6

FIGURE 6. Output forecast of renewable energy on a typical winter day in the RIES.

FIGURE 7

FIGURE 7. Time-of-use price curves.

The situation of carbon emission in different cases is shown in Figure 8, and the flexible load curves such as electricity, heat, cooling, and gas under the IDR as well as their detailed composition in case IV are presented in Figure 9.

FIGURE 8

FIGURE 8. Net emission and capture volume of CO₂.

FIGURE 9

FIGURE 9. Detailed composition of flexible loads. (A) Electrical load, (B) heat load, (C) cool load, and (D) gas load.

As shown in Figure 8, compared with case I and case II, net carbon emissions in case III and case IV are significantly diminished after considering the whole process of carbon emission. Compared with case III, case IV has lower carbon emission and higher carbon capture, which can reduce carbon trading costs to enhance environmental benefits. This is because the IDR of peak cutting and valley filling can reduce some of the purchased power and make the output of new energy units and GTs increase slightly; as a result, the energy consumption of the CCS also increases.

As presented in Figure 9, compared with case III, the CLs, TLs, and FLs of each period in case IV are dynamically changed after considering the IDR, and the maximum peak-to-valley differences of each flexible load are reduced by 36.32%, 35.75%, 26.31%, and 11.42%, thus smoothing the load curves and realizing the process of peak cutting and valley filling. For instance, according to the electrical load in Figure 9A, the variation of CL is reflected in peak periods; the variation of TL is larger in peak and valley periods and smaller in normal periods; and the variation of FL is positive in valley periods and negative in peak and normal periods, which can better reflect the peak and valley characteristics of electricity prices. In addition, the peak-to-valley differences of the flexible electrical load curve in Figure 9A are effectively reduced to make the cost of power purchase continuously decrease. In addition, it is easily known from Figure 9A that the CL of electricity is reflected in peak periods (10:00–11:00 and 17:00–20:00) when the electrical load and price are high, and the reduced load would be converted into incentive subsidies to users. The TL of electricity is shifted from the high-electricity-price periods (10:00–11:00 and 17:00–20:00) to the low-electricity-price periods (23:00–08:00 and 12:00–16:00), which realizes the process of transferring electrical load from peak periods to the valley called peak cutting and valley filling to promote the economy of using the electricity. The FL of electricity converts part of the electrical load into heat load in high electricity price periods to alleviate the mismatch and imbalance between the electrical load demand and the supply capacity. In addition, during the periods of low electricity prices, part of the heat load is transformed into the electrical load for central heating.

The analysis method of the flexible heat load shown in Figure 9B is similar to that of the electrical load. As for the flexible cool and gas loads shown in Figures 9C, D, the system smoothes their load curves by cutting or transferring to realize the process of stabilizing the load fluctuation, which makes the output of the corresponding device more reasonable and perfect.

The aforementioned analysis shows that the IDR under case IV can achieve peak cutting and valley filling to improve economic benefits. Meanwhile, the peak-to-valley differences of each load are reduced by 430.41, 70.60, 33.22 kW, and 29.08 m³ compared with case III, which makes the whole system more flexible in terms of energy purchase so that the environmental benefits are significantly increased.

5.3 Analysis of operating results

Figure 10 shows the optimized operation scheme of storage devices involving the storage battery and heat storage tank of the CSP station. Since the energy storage devices can be used as a part of the flexible load, they have abilities to participate in the whole peak-cutting and valley-filling process. For the storage battery, during the valley period of electrical consumption from 23: 00 to 07: 00, the electrical load and price are both lower, so the storage battery charges at this time to cope with the subsequent peak load. During the peak periods of 10: 00–11: 00 and 17: 00–20: 00, the electrical load and price are both higher, so the storage battery discharges to relieve the pressure in the grid. For the CSP heat storage tank, the main thermal source is the heat-conducting medium of the heat absorber in the optical field, so it will be affected by the operation mode of the CSP station in most periods. For example, from 09: 00 to 12: 00, the heat storage tank can provide thermal energy to the load side. In addition, during the period of 13: 00–16: 00, the light is abundant so that the heat storage tank collects heat from the optical field to meet the demand of thermal use in the evening (17: 00–20: 00). During the period of 23: 00–08: 00, the partial output of the electric boiler will flow into the CSP heat storage tank to realize the heat storage backup process.

FIGURE 10

FIGURE 10. Charge/discharge power and the capacity of storage devices. (A) Storage battery and (B) CSP heat storage tank.

The input of the two energy storage equipment makes the operation of the system more reasonable, and the maximum variation range of power storage and heat storage can reach 770 and 911 kW, respectively, which enables the flexibility of energy storage increase significantly to reduce the cost of energy purchase.

Figure 11 shows the optimized operation results of the four energy sources on the typical winter day, which are analyzed as follows:

1) Balance of electricity

As shown in Figure 11A, during the period of 23: 00–05: 00, the electrical demand is mainly satisfied through WT and GT to reduce the costs of electricity generation and device operation. Since the electrical price in this period is at the trough value, the flexible heat load will be converted into electrical load through the IDR, and the system will purchase a small amount of electricity for the storage battery. During the periods of 06: 00–09: 00 and 12: 00–16: 00, the electricity is afforded by WT, PV, CSP station, and GT in priority, and the insufficient part will be supplemented through electrical purchase. During the periods of 10: 00–11: 00 and 17: 00–22: 00, the electrical demands of users are relatively larger. At this time, except for the output of new energy units and GTs, the electrical shortage can be provided by storage battery and electrical purchase, in turn, to maintain the balance. In addition, the participation of flexible electrical load in IDR during this period can further reduce and transfer the peak electricity consumption, thus significantly decreasing the power impact of electrical devices on the grid and their required configuration capacities.

2) Balance of heat

As shown in Figure 11B, during the period of 22: 00–08: 00, the heat is mainly afforded by GB and WHB while the demand is larger than at any other time, so the total amount and the peak value of heat energy can be reduced by cutting and transferring part of the flexible heat load. In addition, due to the low electricity price at this time, a small amount of the heat that is output by EB will be transported to the load side as well as the residual output of EB will flow into the heat storage tank of the CSP station to meet the lower capacity limit. During the period of 09: 00–21: 00, the heat is mainly provided by WHB, and part of the flexible electrical load will be converted into heat load through the IDR. During the period from 09: 00 to 12: 00, as the light gradually increases, the heat storage tank of the CSP station releases the heat to the load side, which reduces the output of GB to 0. During the period from 13: 00 to 16: 00, WHB takes most of the heat load, and the heat storage tank of the CSP station uses this period to charge as the reserve for the next period. During the period from 17: 00 to 21: 00, the gas prices are high so that the output of GB is decreased, and at the same time, the output of electrical generation in the CSP station is continuously reduced; as a result, the heat storage tank of the CSP station releases the heat to the load side. If the demand cannot be satisfied, EB will provide part of the required heat.

3) Balance of cooling

As shown in Figure 11C, the output of EC completes the supply of cool load in the valley and normal time of electrical price, such as the periods of 21: 00–09: 00 and 12: 00–16: 00. During the periods of high electricity price such as 10: 00–11: 00 and 17: 00–19: 00, the output of EC is less than before and AC bears all the cool load, which is helpful for peak cutting and valley filling. Since there is little demand for cool energy in winter, the cool mode of the system is relatively flexible. As a result, the flexible cooling load only needs to be partially reduced and transferred at the peak.

4) Balance of gas

As shown in Figure 11D, the gas load and the consumption of gas-fired units are satisfied through gas purchase during the periods of 22: 00–08: 00 and 13: 00–16: 00. During the period from 22: 00 to 05: 00, although the gas price is at the trough value, the consumption of gas is generally small, and as a result, the gas purchase is still relatively small. During the period from 06: 00 to 08: 00, the gas purchase is growing and GT becomes the main source of gas consumption and continuously increases its output to satisfy the demands of electricity and heat. During the period from 13: 00 to 16: 00, the consumption of the gas-fired units is maintained within a certain range due to the influence of the climbing rate. At this time, the gas load fluctuates slightly but the amplitude is small; therefore, the corresponding gas purchase is basically at a stable level. During the periods of 09: 00–12: 00 and 17: 00–21: 00, constraints on carbon capture and flue gas treatment make regulating abilities of the system decline, which leads to the abandonment of wind and light. At this time, P2G makes efforts to absorb the abandoned power of wind and light, which can synthesize natural gas to alleviate the imbalance between gas purchase and gas demands. Moreover, due to the high gas price in this period, the flexible gas load will be reduced and transferred accordingly.

FIGURE 11

FIGURE 11. Optimal balance scheduling results on a typical winter day. (A) Electrical load, (B) heat load, (C) cool load, and (D) gas load.

6 Conclusion

In this paper, an operating model of RIES that considers the whole process of carbon emission and IDR is constructed. Meanwhile, a double-layer optimal configuration method based on the improved NSGA-II and MILP is proposed with the participation of the CSP station. The following conclusions can be drawn by setting up four cases:

1) The proposed double-layer optimal configuration model can reasonably optimize the capacities and the output of the devices to obtain the optimal operation scheme. Meanwhile, the participation of the CSP station can improve the coordinated optimization abilities of the RIES.

2) Compared to not considering carbon emissions, introducing the whole process of carbon emission can significantly reduce the net carbon emission, in which the combined operation strategy of wind power–photovoltaic–CSP–carbon capture can coordinate the output of each device and suppress the fluctuation of renewable energy, the process of carbon utilization decreases the storage cost and promotes P2G to absorb the abandoned power of wind and light, and the ladder-type carbon trading mechanism can strictly control the cost of carbon emission to enhance the economy and environmental protection.

3) Compared to not considering IDR, classifying flexible loads including electricity, heat, cooling, and gas into curtailable, transferable, or fungible loads for participating in the IDR can effectively smooth the load curves, reduce the peak-to-valley differences, and realize multi-energy complementarities. Meanwhile, the compensatory mechanism of the IDR can rationally adjust the energy consumption strategies to optimize the energy structure and improve energy efficiency.

4) Compared with the conventional NSGA-II, the improved NSGA-II can obtain the optimal non-dominated solution sets in the upper layer, and its capacity configuration and operation scheme have lower annual planning cost and net pollutant discharge. Meanwhile, the optimized accuracy of the improved NSGA-II is higher, and the convergent speed and ability reach the optimal level.

At present, this paper mainly studies the day-ahead configuration and operation of the single-region RIES. In the follow-up work, the whole process of carbon emission and IDR will be further extended to the energy network with multi-RIES interconnection by considering the effects of source-load uncertainty on capacity configuration and the step-by-step refinement of time scales to realize the decoupling, coordinated optimization, and stable operation of different energies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

XZ: methodology, algorithm finding, and experimental platform provision. JW: programming, simulation, and writing original draft preparation. YZ: algorithm refinement, editing, supervision, and review. HL: raw data provision, algorithm testing, and article refinement. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (No. 52167014); Small and Medium-sized Enterprise Innovation Fund of Gansu Province (No. 22CX3JA002); and Natural Science Foundation Project of Gansu Province (No. 21JR7RA211).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Azzam, S. M., Elshabrawy, T., and Ashour, M. (2023). A bi-level framework for supply and demand side energy management in an islanded microgrid. IEEE Trans. Ind. Inf. 191, 220–231. doi:10.1109/tii.2022.3144154

CrossRef Full Text | Google Scholar

Chen, B., Wu, W., Guo, Q., and Sun, H. (2022). An efficient optimal energy flow model for integrated energy systems based on energy circuit modeling in the frequency domain. Appl. Energy 326, 119923. doi:10.1016/j.apenergy.2022.119923

CrossRef Full Text | Google Scholar

Chen, W., Ma, Y., and Bai, C. (2021). The impact of carbon emission quota allocation regulations on the investment of low-carbon technology in electric power industry under peak-valley price policy. IEEE Trans. Eng. Manage. early access, 1–18. doi:10.1109/tem.2021.3121002

CrossRef Full Text | Google Scholar

Dong, W., Lu, Z., He, L., Geng, L., Guo, X., and Zhang, J. (2022). Low-carbon optimal planning of an integrated energy station considering combined power-to-gas and gas-fired units equipped with carbon capture systems. Int. J. Electr. Power & Energy Syst. 138, 107966. doi:10.1016/j.ijepes.2022.107966

CrossRef Full Text | Google Scholar

Fan, H., Yu, Z., Xia, S., and Li, X. (2021). Review on coordinated planning of source-network-load-storage for integrated energy systems. Front. Energy Res. 9, 641158. doi:10.3389/fenrg.2021.641158

CrossRef Full Text | Google Scholar

Fu, X., Zeng, G., Zhu, X., Zhao, J., Huang, B., and Liu, J. (2022). Optimal scheduling strategy of grid-connected microgrid with ladder-type carbon trading based on Stackelberg game. Front. Energy Res. 10, 961341. doi:10.3389/fenrg.2022.961341

CrossRef Full Text | Google Scholar

Gorman, B. T., Lanzarini-Lopes, M., Johnson, N. G., Miller, J. E., and Stechel, E. B. (2021). Techno-economic analysis of a concentrating solar power plant using redox-active metal oxides as heat transfer fluid and storage media. Front. Energy Res. 9, 734288. doi:10.3389/fenrg.2021.734288

CrossRef Full Text | Google Scholar

Jiang, P., Dong, J., and Huang, H. (2020). Optimal integrated demand response scheduling in regional integrated energy system with concentrating solar power. Appl. Therm. Eng. 166, 114754. doi:10.1016/j.applthermaleng.2019.114754

CrossRef Full Text | Google Scholar

Li, Y., Bu, F., Gao, J., and Li, G. (2022). Optimal dispatch of low-carbon integrated energy system considering nuclear heating and carbon trading. J. Clean. Prod. 378, 134540. doi:10.1016/j.jclepro.2022.134540

CrossRef Full Text | Google Scholar

Shao, C., Ding, Y., Siano, P., and Lin, Z. (2019). A framework for incorporating demand response of smart buildings into the integrated heat and electricity energy system. IEEE Trans. Ind. Electron. 662, 1465–1475. doi:10.1109/tie.2017.2784393

CrossRef Full Text | Google Scholar

Shao, C., Ding, Y., Siano, P., and Song, Y. (2021). Optimal scheduling of the integrated electricity and natural gas systems considering the integrated demand response of energy hubs. IEEE Syst. J. 153, 4545–4553. doi:10.1109/jsyst.2020.3020063

CrossRef Full Text | Google Scholar

Tao, J., Xu, W., and Li, Y. (2019). Optimal operation of integrated energy system combined cooling heating and power based on multi-objective algorithm. Sci. Tech. Engrg. 19(33), 200–205. doi:10.3969/j.issn.1671-1815.2019.33.029

CrossRef Full Text | Google Scholar

Wang, B., Sun, H., and Song, X. (2022). Optimal dispatching modeling of regional power-heat-gas interconnection based on multi-type load adjustability. Front. Energy Res. 10, 931890. doi:10.3389/fenrg.2022.931890

CrossRef Full Text | Google Scholar

Wang, L., Hou, C., Ye, B., Wang, X., Yin, C., and Cong, H. (2021). Optimal operation analysis of integrated community energy system considering the uncertainty of demand response. IEEE Trans. Power Syst. 364, 3681–3691. doi:10.1109/tpwrs.2021.3051720

CrossRef Full Text | Google Scholar

Wang, X., Xin, A., Chen, X., Fang, L., Jia, Q., Ma, L., et al. (2022). Optimal dispatching of ladder-type carbon trading in integrated energy system with advanced adiabatic compressed air energy storage. Front. Energy Res. 10, 933786. doi:10.3389/fenrg.2022.933786

CrossRef Full Text | Google Scholar

Wei, Z., Ma, X., Guo, Y., Wei, P., Lu, B., and Zhang, H. (2022). Optimized operation of integrated energy system considering demand response under carbon trading mechanism. Electr. Power Constr. 431, 1–9. doi:10.12204/j.issn.1000-7229.2022.01.001

CrossRef Full Text | Google Scholar

Yan, M., Li, Y., Chen, G., Zhang, L., Mao, Y., and Ma, C. (2017). A novel flue gas pre-treatment system of post-combustion CO₂ capture in coal-fired power plant. Chem. Eng. Res. Des. 128, 331–341. doi:10.1016/j.cherd.2017.10.005

CrossRef Full Text | Google Scholar

Yan, N., Ma, G., Li, X., and Guerrero, J. M. M. (2023). Low-carbon economic dispatch method for integrated energy system considering seasonal carbon flow dynamic balance. IEEE Trans. Sustain. Energy 141, 576–586. doi:10.1109/tste.2022.3220797

CrossRef Full Text | Google Scholar

Yang, H., Li, M., Jiang, Z., and Zhang, P. (2020). Multi-time scale optimal scheduling of regional integrated energy systems considering integrated demand response. IEEE Access 8, 5080–5090. doi:10.1109/access.2019.2963463

CrossRef Full Text | Google Scholar

Zeng, X., Zhang, J., and Wang, X. (2023). Optimal configuration of regional integrated energy system after taking into account multiple uncertainties and the participation of concentrating solar power stations. High. Volt. Eng. 491, 353–363. doi:10.13336/j.1003-6520.hve.20211326

CrossRef Full Text | Google Scholar

Zhang, D., Yun, Y., Wang, X., He, J., and Dong, H. (2021). Economic dispatch of integrated electricity-heat-gas energy system considering generalized energy storage and concentrating solar power plant. Autom. Electr. Power Syst. 4519, 33–42. doi:10.7500/AEPS20210220002

CrossRef Full Text | Google Scholar

Zhang, D., Zhu, H., Zhang, H., Goh, H. H., Liu, H., and Wu, T. (2022). Multi-objective optimization for smart integrated energy system considering demand responses and dynamic prices. IEEE Trans. Smart Grid 132, 1100–1112. doi:10.1109/tsg.2021.3128547

CrossRef Full Text | Google Scholar

Zhang, G., Xie, P., Huang, S., Chen, Z., Du, M., Tang, N., et al. (2021). Modeling and optimization of integrated energy system for renewable power penetration considering carbon and pollutant reduction systems. Front. Energy Res. 9, 767277. doi:10.3389/fenrg.2021.767277

CrossRef Full Text | Google Scholar

Zhang, N., Hu, Z., Dai, D., Dang, S., Yao, M., and Zhou, Y. (2016). Unit commitment model in smart grid environment considering carbon emissions trading. IEEE Trans. Smart Grid 71, 420–427. doi:10.1109/tsg.2015.2401337

CrossRef Full Text | Google Scholar

Zhao, S., Fang, Y., and Wei, Z. (2019). Stochastic optimal dispatch of integrating concentrating solar power plants with wind farms. Int. J. Electr. Power & Energy Syst. 109, 575–583. doi:10.1016/j.ijepes.2019.01.043

CrossRef Full Text | Google Scholar