A bi-level co-expansion planning of integrated electric and heating system considering demand response

1 Introduction

The widespread deployment of combined heat and power (CHP) units augments the interconnection between the power system and the heating system (Lin et al., 2020; Khatibi et al., 2021). In this context, the integrated electric and heating system (IEHS) has gained massive attention in recent years. IEHS is an important part of the Energy Internet (Long et al., 2022), and it has received extensive attention in both industry and academia.

Extensive studies on the optimal operation of IEHS can be found in the literature. Li et al. (2016) proposed a combined heat and power dispatch (CHPD) method that exploits the flexibility of the district heating system (DHS) for better wind accommodation in the power system. And Xue et al. (2020) developed a heterogeneous decomposition algorithm to tackle the multi-agent problem in CHPD. In addition to economic dispatch, other operation problems such as unit commitment (Anand et al., 2019) and optimal energy flow (Yao et al., 2021) are also studied.

The aforementioned studies are based on given energy facilities which may not provide sufficient operation flexibility once the load increases (Fu et al., 2020). On this account, the co-expansion planning (CEP) of IEHS is another research priority. Li et al. (2021) proposed a CEP method for hybrid concentrating solar power and CHP plant in IEHS. Cheng et al. (2019) developed a CEP model aimed at minimizing cost and emissions. Martinez Cesena et al. (2016) proposed a CEP model for IEHS to address long-term price uncertainty. Cao et al. (2020) proposed a data-driven method to solve the CEP with uncertainty.

However, the aforementioned studies usually formulate the CEP problem as a single-layer model, which is unreasonable. In industry practice, the investment and operation are determined by the generation company and system operator, respectively (Pineda and Morales, 2019). The generation company determines the investment strategy, which is submitted to the system operator. And the system operator develops a cost-effective operation strategy, which regulates the unit operation. Thus, it is urgent to develop a CEP framework to capture this feature.

Considering the impact of the operation strategies on investment decisions, fully exploiting flexibility resources in IEHS operation will result in more beneficial planning strategies. Demand response (DR) is an important way to mobilize flexible resources on the demand side (Huang et al., 2019; Gjorgievski et al., 2021), which has been studied extensively in quantitative evaluation (D’hulst et al., 2015), economic dispatch (Alipour et al., 2019; Majidi et al., 2019) and unit commitment (Mansour-Saatloo et al., 2020). In contrast, few works focus on the performance of DR on CEP problems. DR could optimize the load curves by intentionally modifying the energy consumption patterns of users to facilitate IEHS in balancing supply and demand (Yin et al., 2016). Thus, DR can be regarded as a potential tool to postpone or even cancel unnecessary facilities planning.

Accordingly, the focus of this paper is to develop a bi-level co-expansion planning (BLCEP) model with DR constraints for IEHS. The main contributions of our work are as follows:

1) We develop a BLCEP framework of IEHS in this paper. The investment of power and heat sources are optimized in the upper level, while the operation problem is optimized in the lower level. In this way, the optimal results can be obtained through the game of the upper-level problem and the lower-level problem.

2) DR is introduced into the BLCEP model, and its cost is considered in the objective function. The electric and heat load curve is optimized by DR, and the operation flexibility of the demand side is exploited. In this way, unnecessary investment in IEHS facilities is prevented.

3) The proposed BLCEP model is transformed into a single-layer model by replacing the lower-level operation model with corresponding Karush-Kuhn-Tucker (KKT) conditions. Therefore, the proposed model can be solved directly with commercial solvers such as CPLEX and Gurobi.

The remainder of this paper is as follows. In Section 2, we proposed a BLCEP framework of IEHS. Based on the framework, a BLCEP model with DR constraints is given in Section 3. In Section 4, the proposed BLCEP model is transformed into a single-level mixed integer linear optimization problem. In Section 5, numerical simulations of the P6H8 system are performed to demonstrate the effectiveness of the proposed model. And the conclusion is given in Section 6.

2 Framework

A typical IEHS is comprised of an electric power system (EPS) and DHSs. In this regard, the expansion planning for IEHS includes planning for both power and heat sources. Figure 1 depicts the BLCEP framework of IEHS. An investment model is presented in the upper level, which is utilized to determine the optimal investment strategy for energy facilities. Based on the investment decision given by the upper level, the IEHS operation problem is optimized to minimize the operation costs at the lower level and feedback to the upper level.

FIGURE 1

FIGURE 1. Framework of BLCEP of IEHS problem.

Noting that the investment model in the upper level and the operation model in the lower level are mutually influenced. The less investment in expansion, the stricter operation conditions, and the larger operation cost. Similarly, lower operation costs require more relaxed operation conditions, i.e., more investment in energy facilities. Hence, the optimal solution of bi-level programming can be regarded as the gaming between the investment level and operation level (Zeng et al., 2017; Wang et al., 2021).

3 Bi-level co-expansion planning of integrated electric and heating system model formulation

Based on the framework given in Section 2, the BLCEP model is formulated as follows.

3.1 Investment model in the upper level

Denote $x_{g,y}^G$, $x_{w,y}^{WD}$, $x_{e,y}^{EB}$ as binary variables representing the installation status of generators, wind farms and electric boilers, respectively. The investment model in the upper level aims to minimize the total cost $C^{UL}$ as

$\min C^{UL}={\displaystyle \sum }_{y=1}^Y\frac{1}{{\left(1+dr\right)}^{y-1}}\left(C_I\left(y\right)+C_O\left(y\right)\right)(1)$

The first term indicates the investment cost of energy facilities which consist of generators, wind farms and electric boilers, i.e.,

$\begin{array}{l}{C}_I\left(y\right)={\displaystyle \sum }_{g\in {\mathcal{G}}_{CT}^C\cup {\mathcal{G}}_{CHP}^C}\left(c_g^G\cdot \left(x_{g,y}^G-x_{g,y-1}^G\right)\right)\\ +{\displaystyle \sum }_{w\in {\mathcal{W}}^C}\left(c_w^{WD}\cdot \left(x_{w,y}^{WD}-x_{w,y-1}^{WD}\right)\right)\\ +{\displaystyle \sum }_{e\in {\mathcal{E}}^C}\left(c_e^{EB}\cdot \left(x_{e,y}^{EB}-x_{e,y-1}^{EB}\right)\right)\end{array}(2)$

where dr indicates the discount rate, and it generally takes a value between 6% and 8% (Mu et al., 2020). The second term of (1) indicates operation cost in $y$ year, i.e.,

$C_O\left(y\right)={\displaystyle \sum }_{d=1}^{N{D}_y}{\displaystyle \sum }_{t=1}^{{\mathcal{T}}_d}\left(C_{CHP}\left(p_{g,t}^G,h_{g,t}^G\right)+C_{CT}\left(p_{g,t}^G\right)+C_{WD}\left(p_{w,t}^{WD}\right)+C_{DR}\left(\Delta p_{b,t}^{LD},\Delta h_{i,t}^{LD}\right)\right)(3)$

The upper-level model is subjected to investment constraints. 1) duplicate investment of energy facilities is prohibited (4)–(6). 2) annual investment cost is limited by budget (7).

${\displaystyle \sum }_{y=1}^Y{x}_{g,y}^G=1,x_{g,y-1}^G\le x_{g,y}^G,\forall g\in {\mathcal{G}}_{CT}^C\cup {\mathcal{G}}_{CHP}^C,(4)$

${\displaystyle \sum }_{y=1}^Y{x}_{w,y}^{WD}=1,x_{w,y-1}^{WD}\le x_{w,y}^{WD},\forall w\in {\mathcal{W}}^C,(5)$

${\displaystyle \sum }_{y=1}^Y{x}_{e,y}^{EB}=1,x_{e,y-1}^{EB}\le x_{e,y}^{EB},\forall e\in {\mathcal{E}}^C,(6)$

$C_I\left(y\right)\le {\overline{C}}_I,\forall y\in \left\{\mathrm{1,2},\cdots ,Y\right\},(7)$

where ${\mathcal{G}}_{CT}^C$ is the set of candidate conventional thermal units, ${\mathcal{G}}_{CHP}^C$ is the set of candidate CHP units, ${\mathcal{W}}^C$ is the set of candidate wind farms, and ${\mathcal{E}}^C$ is the set of candidate electric boilers. And ${\overline{C}}_I$ indicates the maximum annual investment.

3.2 Operation model in the lower level

The objective function $C^{LL}$ in the lower level aims to minimize the total operation cost consisting of generator fuel cost wind curtailment cost and DR cost as

$\begin{array}{l}\min C^{LL}={\displaystyle \sum }_{y=1}^Y\frac{1}{{\left(1+dr\right)}^{y-1}}({\displaystyle \sum }_{d=1}^{N{D}_y}{\displaystyle \sum }_{t=1}^{{\mathcal{T}}_d}{C}_{CHP}\left(p_{g,t}^G,h_{g,t}^G\right)+{\displaystyle \sum }_{d=1}^{N{D}_y}{\displaystyle \sum }_{t=1}^{{\mathcal{T}}_d}{C}_{CT}\left(p_{g,t}^G\right)\\ +{\displaystyle \sum }_{d=1}^{N{D}_y}{\displaystyle \sum }_{t=1}^{{\mathcal{T}}_d}{C}_{WD}\left(p_{w,t}^{WD}\right)+{\displaystyle \sum }_{d=1}^{N{D}_y}{\displaystyle \sum }_{t=1}^{{\mathcal{T}}_d}{C}_{DR}\left(\Delta p_{b,t}^{LD},\Delta h_{i,t}^{LD}\right))\end{array}(8)$

where $C_{CHP}\left(\cdot \right)$, $C_{CT}\left(\cdot \right)$, $C_{WD}\left(\cdot \right)$indicate the cost function of CHP units, conventional thermal units, and wind farms, respectively. $p_{g,t}^G$, $h_{g,t}^G$, $p_{w,t}^{WD}$ represent the output of generators and wind farms. $C_{DR}\left(\cdot \right)$is the cost function to compensate users for participating in DR as

$C_{DR}\left(\Delta p_{b,t}^{LD},\Delta h_{i,t}^{LD}\right)={\lambda }_{ELD}\cdot \left|\Delta p_{b,t}^{LD}\right|+{\lambda }_{HLD}\cdot \left|\Delta h_{i,t}^{LD}\right|(9)$

where $\Delta p_{b,t}^{LD}$ and $\Delta h_{i,t}^{LD}$ are load shifting by DR. ${\lambda }_{ELD}$ and ${\lambda }_{HLD}$ are the cost coefficient of power and heat load shifting, respectively.

In addition to the objective, operation constraints of DHS and EPS would have to be considered in the model as follows.

3.2.1 District heating system operation constraints

Generally, heat is generated by heating facilities such as CHP units, electric boilers in heat stations, i.e.,

${\displaystyle \sum }_{g\in {\kappa }_i^{CHP}}{h}_{g,t}^G+{\displaystyle \sum }_{e\in {\kappa }_i^{EB}}{h}_{e,t}^{EB}=h_{i,t}^{HS},\forall i\in {\mathcal{N}}^{HS},t\in {\mathcal{T}}_d,(10)$

where $h_{e,t}^{EB}$ indicates heat generated in electric boiler e, and $h_{i,t}^{HS}$ indicates heat generated in heat station connecting node i. ${\kappa }_i^{CHP}$ and ${\kappa }_i^{EB}$ denote the set of CHP unit and electric boiler connecting heat node i, respectively. The output of CHP units is limited by their operation feasible region ${\mathcal{R}}_g^G$, i.e.,

$\left(p_{g,t}^G,h_{g,t}^G\right)\in {\mathcal{R}}_g^G,\forall g\in {\mathcal{G}}_{CHP},t\in {\mathcal{T}}_d.(11)$

As for electric boilers, they can be formulated as

$h_{e,t}^{EB}={\eta }_e^{EB}\cdot p_{e,t}^{EB}\forall e\in \mathcal{E},t\in {\mathcal{T}}_d,(12)$

$0\le p_{e,t}^{EB}\le {\overline{P}}_e^{EB},\forall e\in \mathcal{E}\mathcal{\backslash }{\mathcal{E}}^C,t\in {\mathcal{T}}_d,(13)$

$0\le p_{e,t}^{EB}\le x_{e,y}^{EB}{\overline{P}}_e^{EB},\forall e\in {\mathcal{E}}^C,t\in {\mathcal{T}}_d,(14)$

where $\mathcal{E}$ is the set of all electric boilers. ${\eta }_e^{EB}$ is the efficiency of electricity to heat. And ${\overline{P}}_e^{EB}$ is maximum consumption of electric boilers e. The heat generated in heat stations will be delivered to head loads $h_{i,t}^{LD}$ through hot water, i.e.,

$h_{i,t}^{HS}=c\cdot m_{i,t}^{HS}\cdot \left({\tau }_{i,t}^{NS}-{\tau }_{i,t}^{NR}\right),\forall i\in {\mathcal{N}}^{HS},t\in {\mathcal{T}}_d,(15)$

$h_{i,t}^{LD}=c\cdot m_{i,t}^{LD}\cdot \left({\tau }_{i,t}^{NS}-{\tau }_{i,t}^{NR}\right),\forall i\in {\mathcal{N}}^{LD},t\in {\mathcal{T}}_d,(16)$

where $c$ indicate the density of water. $m_{i,t}^{HS}$ and $m_{j,t}^{LD}$ is the mass flow rate in heat station and heat load connecting with node i. ${\tau }_{i,t}^{NS}$ and ${\tau }_{i,t}^{NR}$ denote temperature at heat node i of the supply and return pipelines. To guarantee the heating quality, the supply temperature in heat stations and return temperature in heat loads should be limited:

${\underset{¯}{\tau }}_i^{NS}\le {\tau }_{i,t}^{NS}\le {\overline{\tau }}_i^{NS},\forall i\in {\mathcal{N}}^{HS},t\in {\mathcal{T}}_d,(17)$

${\underset{¯}{\tau }}_i^{NR}\le {\tau }_{i,t}^{NR}\le {\overline{\tau }}_i^{NR},\forall i\in {\mathcal{N}}^{LD},t\in {\mathcal{T}}_d.(18)$

Let ${\mathcal{S}}_i^{p+}$ and ${\mathcal{S}}_i^{p-}$ denote pipelines starting and ending of heat node i, respectively. Similarly, ${\mathcal{S}}_i^{LD}$ and ${\mathcal{S}}_i^{HS}$ denote the heat loads and heat stations connected with heat node i, respectively. According to the law of continuity, the total mass flow rate injected into the heat node is 0, i.e.,

$\{\begin{array}{l}{\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p+}}{m}_{p,t}^S+{\displaystyle \sum }_{j\in {\mathcal{S}}_i^{LD}}{m}_{j,t}^{LD}={\displaystyle \sum }_{j\in {\mathcal{S}}_i^{HS}}{m}_{j,t}^{HS}+{\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p-}}{m}_{p,t}^S,\forall i\in \mathcal{N},t\in {\mathcal{T}}_d,\\ {\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p+}}{m}_{p,t}^R+{\displaystyle \sum }_{j\in {\mathcal{S}}_i^{LD}}{m}_{j,t}^{LD}={\displaystyle \sum }_{j\in {\mathcal{S}}_i^{HS}}{m}_{j,t}^{HS}+{\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p-}}{m}_{p,t}^R,\forall i\in \mathcal{N},t\in {\mathcal{T}}_d,\end{array}(19)$

where $m_{p,t}^S$ and $m_{p,t}^R$ indicate mass flow rate in supply and return pipelines, respectively. $m_{j,t}^{LD}$ and $m_{j,t}^{HS}$ indicate the mass flow rate in the load and heat station. Traditionally, the DHS operates in quality regulation mode, which means constant mass flow rate and variable heating temperature regulation strategy (Wang et al., 2019). Based on this, the mass flow rate is fixed in this paper to linearize the DHS model (Shao et al., 2020).

According to the energy conservation, the temperatures of mass flow rate would mix when mass flow into the same heat node, i.e.,

$\{\begin{array}{l}{\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p+}}\left({\tau }_{p,t}^{PS,out}\cdot m_{p,t}^S\right)={\tau }_{i,t}^{NS}\cdot {\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p+}}{m}_{p,t}^S,\forall i\in \mathcal{N},t\in {\mathcal{T}}_d,\\ {\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p-}}\left({\tau }_{p,t}^{\mathit{PR},out}\cdot m_{p,t}^R\right)={\tau }_{i,t}^{NR}\cdot {\displaystyle \sum }_{p\in {\mathcal{S}}_i^{p-}}{m}_{p,t}^R,\forall i\in \mathcal{N},t\in {\mathcal{T}}_d,\end{array}(20)$

where ${\tau }_{p,t}^{PS,out}$ and ${\tau }_{p,t}^{\mathit{PR},out}$ denote the temperature at the inlet of the supply and return pipeline p, respectively. Besides, the temperature at the heat node is the same as those at the inlet of the supply and return pipeline, i.e.,

$\{\begin{array}{l}{\tau }_{p,t}^{PS,in}={\tau }_{i,t}^{NS},\forall p\in {\mathcal{S}}_i^{p+},i\in \mathcal{N},t\in {\mathcal{T}}_d,\\ {\tau }_{p,t}^{\mathit{PR},in}={\tau }_{i,t}^{NR},\forall p\in {\mathcal{S}}_i^{p-},i\in \mathcal{N},t\in {\mathcal{T}}_d,\end{array}(21)$

Considering heating dissipation during heating transfer, the outlet temperature of mass flow rate in pipelines is lower than the inlet temperature.

$\{\begin{array}{l}{\tau }_{p,t}^{PS,out}={\tau }_{p,t}^{am}+\left({\tau }_{p,t}^{PS,in}-{\tau }_{p,t}^{am}\right)\cdot e^{-\frac{{\lambda }_p{L}_p}{m_{p,t}^S{A}_p\rho c}},\forall p\in {\mathcal{S}}_i^p,t\in {\mathcal{T}}_d,\\ {\tau }_{p,t}^{\mathit{PR},out}={\tau }_{p,t}^{am}+\left({\tau }_{p,t}^{\mathit{PR},in}-{\tau }_{p,t}^{am}\right)\cdot e^{-\frac{{\lambda }_p{L}_p}{m_{p,t}^R{A}_p\rho c}},\forall p\in {\mathcal{S}}_i^p,t\in {\mathcal{T}}_d,\end{array}(22)$

where $\rho$ indicate specific heat capacity. $A_p$, ${\lambda }_p$, and $L_p$ indicate the cross-sectional area, heat conductivity coefficient, and length of pipeline p, respectively.

3.1.2 Electric power system operation constraints

Let ${\mathcal{S}}_b^{CT}$ denote the set of conventional thermal units connecting bus b, and the set of all conventional thermal units is denoted by ${\mathcal{G}}_{CT}:={\cup }_{b\in \mathcal{B}}{\mathcal{S}}_b^{CT}$. Similarly, the set of CHP units connecting bus b, all CHP units, wind farms connecting bus b, and all wind farms are defined by ${\mathcal{S}}_b^{CHP}$, ${\mathcal{G}}_{CHP}$, ${\mathcal{S}}_b^{WD}$ and $\mathcal{W}$, respectively. By introducing binary investment decision variables into the direct flow model, the operation constraints on the power system are given as follows, which consist of power flow balance (23), transmission capacity (24), generator output (25)–(28), ramp up/down (29)–(32), and the spinning reserve limit (33).

${\displaystyle \sum }_{g\in {\mathcal{G}}_{TU}\cup {\mathcal{G}}_{CHP}}{p}_{g,t}^G+{\displaystyle \sum }_{w\in \mathcal{W}}{p}_{w,t}^{WD}={\displaystyle \sum }_{e\in \mathcal{E}}{p}_{e,t}^{EB}+{\displaystyle \sum }_{b\in \mathcal{B}}{D}_{b,t},\forall t\in {\mathcal{T}}_d,(23)$

$\left|{\displaystyle \sum }_{b\in \mathcal{B}}S{F}_{b,l}\cdot \left({\displaystyle \sum }_{g\in {\mathcal{S}}_b^{TU}\cup {\mathcal{S}}_b^{CHP}}{p}_{g,t}^G+{\displaystyle \sum }_{w\in {\mathcal{S}}_b^{WD}}{p}_{w,t}^{WD}-{\displaystyle \sum }_{e\in {\mathcal{S}}_b^{EB}}{p}_{e,t}^{EB}-D_{b,t}\right)\right|\le F_l,\forall l\in \mathcal{L},t\in {\mathcal{T}}_d,(24)$

$0\le p_{w,t}^{WD}\le {\overline{P}}_w^{WD},\forall w\in \mathcal{W}\mathcal{\backslash }{\mathcal{W}}^C,t\in {\mathcal{T}}_d,(25)$

$0\le p_{w,t}^{WD}\le x_{w,a,y}^{WD}{\overline{P}}_w^{WD},\forall w\in {\mathcal{W}}^C,t\in {\mathcal{T}}_d,(26)$

${\underset{¯}{P}}_g^G\le p_{g,t}^G\le {\overline{P}}_g^G,\forall g\in \left({\mathcal{G}}_{CT}\backslash {\mathcal{G}}_{CT}^C\right)\cup \left({\mathcal{G}}_{CHP}\backslash {\mathcal{G}}_{CHP}^C\right),t\in {\mathcal{T}}_d,(27)$

$x_{g,y}^G{\underset{¯}{P}}_g^G\le p_{g,t}^G\le x_{g,y}^G{\overline{P}}_g^G,\forall g\in {\mathcal{G}}_{CT}^C\cup {\mathcal{G}}_{CHP}^C,t\in {\mathcal{T}}_d,(28)$

$0\le r{u}_{g,t}\le R{U}_g,0\le r{d}_{g,t}\le R{D}_g,\forall g\in {\mathcal{G}}_{CT}\cup {\mathcal{G}}_{CHP},t\in {\mathcal{T}}_d,(29)$

$r{u}_{g,t}\le {\overline{P}}_g^G-p_{g,t}^G,r{d}_{g,t}\le p_{g,t}^G-{\underset{¯}{P}}_g^G,\forall g\in {\mathcal{G}}_{CT}\cup {\mathcal{G}}_{CHP},t\in {\mathcal{T}}_d,(30)$

$r{u}_{g,t}\le x_{g,y}^G{\overline{P}}_g^G-p_{g,t}^G,r{d}_{g,t}\le p_{g,t}^G-x_{g,y}^G{\underset{¯}{P}}_g^G,\forall g\in {\mathcal{G}}_{CT}^C\cup {\mathcal{G}}_{CHP}^C,t\in {\mathcal{T}}_d,(31)$

$-R{D}_g\cdot \Delta t\le p_{g,t}^G-p_{g,t-1}^G\le R{U}_g\cdot \Delta t,\forall g\in {\mathcal{G}}_{CT}\cup {\mathcal{G}}_{CHP},t\in {\mathcal{T}}_d,(32)$

${\displaystyle \sum }_{g\in {\mathcal{G}}_{TU}\cup {\mathcal{G}}_{CHP}}r{u}_{g,t}\ge SRU,{\displaystyle \sum }_{g\in {\mathcal{G}}_{TU}\cup {\mathcal{G}}_{CHP}}r{d}_{g,t}\ge SRD,\forall t\in {\mathcal{T}}_d,(33)$

where $D_{b,t}$ refers to the actual value of power load at bus b, and $S{F}_{b,l}$ is shift factor of bus b to line l. ${\overline{P}}_w$ is the predicted output of wind farms. ${\underset{¯}{P}}_g^G$ and ${\overline{P}}_g^G$ are the lower and upper limits of the power output of generator g, respectively. $R{U}_g$ and $R{D}_g$ are upward and downward ramping capacities of generator g, respectively. $r{u}_{g,t}$ and $r{d}_{g,t}$ are the upward and downward spinning reserve of generator g. $SRU$ and $SRD$ are EPS upward and downward spinning reserve capacities requirement.

3.2.3 Demand response constraints

The demand response is considered in this paper. Denote ${\dot{D}}_{b,t}$ and ${\dot{h}}_{i,t}^{LD}$ as the predicted value of power and heat load. Thus, the actual value of power and heat load as

$D_{b,t}={\dot{D}}_{b,t}+\Delta p_{b,t}^{LD},h_{i,t}^{LD}={\dot{h}}_{i,t}^{LD}+\Delta h_{i,t}^{LD},\forall b\in \mathcal{B},i\in {\mathcal{N}}^{LD},t\in {\mathcal{T}}_d.(34)$

To guarantee user satisfaction, the total amount of load in the subperiod $\mathcal{T}$ is constant, and the demand change is limited.

${\displaystyle \sum }_{t\in \mathcal{T}}\Delta p_{b,t}^{LD}=0,{\displaystyle \sum }_{t\in \mathcal{T}}\Delta h_{i,t}^{LD}=0,\forall b\in \mathcal{B},i\in {\mathcal{N}}^{LD},(35)$

$\left|\Delta p_{b,t}^{LD}\right|\le {ϵ}_{LD}\cdot {\dot{D}}_{b,t},\left|\Delta h_{i,t}^{LD}\right|\le {ϵ}_{HD}\cdot {\dot{h}}_{i,t}^{LD},\forall b\in \mathcal{B},i\in {\mathcal{N}}^{LD},t\in {\mathcal{T}}_d,(36)$

where ${ϵ}_{LD}$ and ${ϵ}_{HD}$ are power and heat load change rate, respectively.

4 Solving strategy

For brevity, the BLCEP model proposed in Section 3 is rewritten in matrix form, as follows

$\begin{array}{l}\underset{\mathit{x}\in \left\{\mathrm{0,1}\right\}}{\min }{\mathit{a}}^{\mathbf{T}}\mathit{x}+{\mathit{b}}^{\mathbf{T}}\mathit{y}\\ s.t.{\mathit{c}}^{\mathbf{T}}\mathit{x}\le \mathit{d},\\ \underset{\mathit{y}\in \mathbf{R}}{\min }{\mathit{f}}^{\mathbf{T}}\mathit{y}\\ s.t.{\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}\le \mathit{e}:\mathit{\lambda }.\end{array}(37)$

where $\mathit{x}$ represents binary investment variables in upper level, including $x_{g,y}^G$, $x_{w,y}^{WD}$ and $x_{e,y}^{EB}$. $\mathit{y}$ represents continuous operation variables in lower level, including $p_{g,t}^G$, $p_{w,t}^{WD}$, $\Delta p_{b,t}^{LD}$, $h_{g,t}^G$, $h_{e,t}^{EB}$, $h_{i,t}^{HS}\Delta h_{i,t}^{LD}$, ${\tau }_{i,t}^{NS}$, ${\tau }_{i,t}^{NR}$, ${\tau }_{p,t}^{PS,out}$ and ${\tau }_{p,t}^{\mathit{PR},out}$. ${\mathit{a}}^{\mathbf{T}}\mathit{x}+{\mathit{b}}^{\mathbf{T}}\mathit{y}$ represents (1). ${\mathit{c}}^{\mathbf{T}}\mathit{x}\le \mathit{d}$ represents (4)–(7). ${\mathit{f}}^{\mathbf{T}}\mathit{y}$ represents (8). ${\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}\le \mathit{e}$ represents (10)–(36). $\mathit{\lambda }$ is a vector of dual variables of lower-level constraint ${\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}\le \mathit{e}$.

Since the lower-level operation problem with fixed investment decision is linear, it can be replaced by corresponding KKT conditions as follows:

$\begin{array}{l}\underset{\mathit{x}\in \left\{\mathrm{0,1}\right\},\mathit{y}\in \mathbf{R}}{\min }{\mathit{a}}^{\mathbf{T}}\mathit{x}+{\mathit{b}}^{\mathbf{T}}\mathit{y}\\ s.t.{\mathit{a}}^{\mathbf{T}}\mathit{x}\le \mathit{d},\\ {\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}\le \mathit{e},\\ {\mathit{f}}^{\mathbf{T}}+\mathit{\lambda }{\mathit{s}}^{\mathbf{T}}=0,\\ \mathit{\lambda }\ge 0,\\ \mathit{\lambda }\left({\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}-\mathit{e}\right)=0.\end{array}(38)$

In (38), the non-linear constraint $\mathit{\lambda }\left({\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}-\mathit{e}\right)=0$ can be handled using the Fortuny-Amat transformation as

$\{\begin{array}{l}\mathit{\lambda }\le \mathit{\mu }\cdot {\rm M},\\ \mathit{e}-{\mathit{r}}^{\mathbf{T}}\mathit{x}-{\mathit{s}}^{\mathbf{T}}\mathit{y}\le \left(1-\mathit{\mu }\right)\cdot \mathit{M},\\ \mathit{\mu }\in \left\{\mathrm{0,1}\right\}.\end{array}(39)$

where $\mathit{\mu }$ a vector of binary auxiliary variable, and ${\rm M}$ is a large enough parameter. Replace$\mathit{\lambda }\left({\mathit{r}}^{\mathbf{T}}\mathit{x}+{\mathit{s}}^{\mathbf{T}}\mathit{y}-\mathit{e}\right)=0$ in (38) with (39), the proposed BLCEP model is transformed into a single-level mixed integer linear optimization problem that can be solved directly by commercial solvers.

5 Case study

In this section, a modified P6H8 system is utilized to testify the proposed model. The modified P6H8 system is shown in Figure 2, which comprised a six-buses electric power system and an eight-nodes district heating system. The planning horizon is 10 years, while each year is divided into three typical days (transition season, summer, and winter) with hourly timesteps. The candidate facilities data is presented in Table 1. The annual growth rate of the power load is 2.5%, and that of the heat load is 4%. The discount rate is 8%. The load change rate ${ϵ}_{LD}$ and ${ϵ}_{HD}$ are 15%.

FIGURE 2

FIGURE 2. A modified P6H8 system.

TABLE 1

TABLE 1. Candidate facilities data.

To compare and analyze the proposed model, two cases are employed as follows:

Case 1: Determine the planning strategy of IEHS using the BLCEP method without DR.

Case 2: Determine the planning strategy of IEHS using the BLCEP method with DR constraints developed in this paper.

All the test is performed on 11th Gen Intel(R) Core (TM) i7-1165G7 at 2.80GHz CPU, 16GB RAM system. The proposed BLCEP model is coded using Matlab R2020a with Gurobi 9.5.0 Solver.

5.1 Case 1

In Case 1, the IEHS is planned and operated in a coordinated mode without demand response. The installation year of candidate facilities is presented in Table 2, and the co-planning cost is presented in Table 3. During the whole planning horizon, five facilities (i.e., G1, G2, C1, W1, and E1) would have to be installed to satisfy increased loads at a total cost of $ 1.37×10⁸.

TABLE 2

TABLE 2. Installation year of candidate facilities.

TABLE 3

TABLE 3. Cost comparison between Case 1 and 2.

5.2 Case 2

We consider a BLCEP method that considers demand response in Case 2. As shown in Figure 3, the peak power loads in the first summer and winter are transferred to valley periods by DR to facilitate the supply-demand balance of power. Due to the tight coupling of power and heat, valley heat loads can be shifted to peak periods in summer for dispatching CHP units to supply peak power loads. Hence, G2 is canceled, as shown in Table 2. As shown in Figure 4. the peak heat loads in the third winter are transferred to valley periods by DR to facilitate the supply-demand balance of heat, so that C1 is postponed.

FIGURE 3

FIGURE 3. Power loads curves in the first year.

FIGURE 4

FIGURE 4. Heat loads curves in the third year.

As shown in Table 2, Case 2 cancels G2 installation and postpones C1 installation in comparison to Case 1. Hence, the investment cost could be reduced by $ 0.43 × 10⁸ even if G1 is installed ahead. Considering the shifted loads would like to be supplied by cost-effective units, the generator fuel cost is reduced by $ 0.13 × 10⁸. Besides, wind curtailment costs are reduced by $ 2.02 × 10⁶ by exploiting the flexibility of loads to consume wind power. As a result, the total cost is reduced by $ 0.44 × 10⁸.

6 Conclusion

The proposed BLCEP framework achieves the game optimization of investment and operation, which can be utilized in CEP in other systems. DR is introduced in the BLCEP model and analyzed in the case study. It is found that DR would help the generation company with decision support to develop a more beneficial planning strategy.

The BLCEP model proposed in this paper is oriented to deterministic scenarios. However, intermittent renewable power outputs and variable load demands will bring challenges to IEHS economic and safe operation (Fu, 2022), ultimately affecting planning strategies. In our future works, the uncertainty of renewable power and load demands will be considered in our model.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://docs.google.com/spreadsheets/d/13TcFGtE5cLdSJiCrPx0PsfU3gLrKAqE4/edit?usp=sharing&ouid=115924522171410221124&rtpof=true&sd=true.

Author contributions

YD: Conceptualization, Methodology, Formal Analysis, Writing-Original Draft YX: Conceptualization, Supervision, Writing-Reviewing and Editing LL: Writing- Validation, Reviewing and Editing CY: Visualization JZ: Validation.

Conflict of interest

CY was employed by Guodian Nanjing Automation Co., Ltd.

The remaining 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

Alipour, M., Zare, K., Seyedi, H., and Jalali, M. (2019). Real-time price-based demand response model for combined heat and power systems. Energy 168, 1119–1127. doi:10.1016/j.energy.2018.11.150

CrossRef Full Text | Google Scholar

Anand, H., Narang, N., and Dhillon, J. S. (2019). Multi-objective combined heat and power unit commitment using particle swarm optimization. Energy 172, 794–807. doi:10.1016/j.energy.2019.01.155

CrossRef Full Text | Google Scholar

Cao, Y., Wei, W., Wang, J., Mei, S., Shafie-khah, M., and Catalao, J. P. S. (2020). Capacity planning of energy hub in multi-carrier energy networks: A data-driven robust stochastic programming approach. IEEE Trans. Sustain. Energy 11, 3–14. doi:10.1109/TSTE.2018.2878230

CrossRef Full Text | Google Scholar

Cheng, Y., Zhang, N., Lu, Z., and Kang, C. (2019). Planning multiple energy systems toward low-carbon society: A decentralized approach. IEEE Trans. Smart Grid 10, 4859–4869. doi:10.1109/TSG.2018.2870323

CrossRef Full Text | Google Scholar

D’hulst, R., Labeeuw, W., Beusen, B., Claessens, S., Deconinck, G., and Vanthournout, K. (2015). Demand response flexibility and flexibility potential of residential smart appliances: Experiences from large pilot test in Belgium. Appl. Energy 155, 79–90. doi:10.1016/j.apenergy.2015.05.101

CrossRef Full Text | Google Scholar

Fu, X., Guo, Q., and Sun, H. (2020). Statistical machine learning model for stochastic optimal planning of distribution networks considering a dynamic correlation and dimension reduction. IEEE Trans. Smart Grid 11, 2904–2917. doi:10.1109/TSG.2020.2974021

CrossRef Full Text | Google Scholar

Fu, X. (2022). Statistical machine learning model for capacitor planning considering uncertainties in photovoltaic power. Prot. Control Mod. Power Syst. 7, 5. doi:10.1186/s41601-022-00228-z

CrossRef Full Text | Google Scholar

Gjorgievski, V. Z., Markovska, N., Abazi, A., and Duić, N. (2021). The potential of power-to-heat demand response to improve the flexibility of the energy system: An empirical review. Renew. Sustain. Energy Rev. 138, 110489. doi:10.1016/j.rser.2020.110489

CrossRef Full Text | Google Scholar

Huang, W., Zhang, N., Kang, C., Li, M., and Huo, M. (2019). From demand response to integrated demand response: Review and prospect of research and application. Prot. Control Mod. Power Syst. 4, 12. doi:10.1186/s41601-019-0126-4

CrossRef Full Text | Google Scholar

Khatibi, M., Bendtsen, J. D., Stoustrup, J., and Molbak, T. (2021). Exploiting power-to-heat assets in district heating networks to regulate electric power network. IEEE Trans. Smart Grid 12, 2048–2059. doi:10.1109/TSG.2020.3044348

CrossRef Full Text | Google Scholar

Li, X., Wu, X., Gui, D., Hua, Y., and Guo, P. (2021). Power system planning based on CSP-CHP system to integrate variable renewable energy. Energy 232, 121064. doi:10.1016/j.energy.2021.121064

CrossRef Full Text | Google Scholar

Li, Z., Wu, W., Shahidehpour, M., Wang, J., and Zhang, B. (2016). Combined heat and power dispatch considering pipeline energy storage of district heating network. IEEE Trans. Sustain. Energy 7, 12–22. doi:10.1109/TSTE.2015.2467383

CrossRef Full Text | Google Scholar

Lin, C., Wu, W., Wang, B., Shahidehpour, M., and Zhang, B. (2020). Joint commitment of generation units and heat exchange stations for combined heat and power systems. IEEE Trans. Sustain. Energy 11, 1118–1127. doi:10.1109/TSTE.2019.2917603

CrossRef Full Text | Google Scholar

Long, H., Fu, X., Kong, W., Chen, H., Zhou, Y., and Yang, F. (2022). Key technologies and applications of rural energy internet in China. Inf. Process. Agric. doi:10.1016/j.inpa.2022.03.001

CrossRef Full Text | Google Scholar

Majidi, M., Mohammadi-Ivatloo, B., and Anvari-Moghaddam, A. (2019). Optimal robust operation of combined heat and power systems with demand response programs. Appl. Therm. Eng. 149, 1359–1369. doi:10.1016/j.applthermaleng.2018.12.088

CrossRef Full Text | Google Scholar

Mansour-Saatloo, A., Agabalaye-Rahvar, M., Mirzaei, M. A., Mohammadi-Ivatloo, B., Abapour, M., and Zare, K. (2020). Robust scheduling of hydrogen based smart micro energy hub with integrated demand response. J. Clean. Prod. 267, 122041. doi:10.1016/j.jclepro.2020.122041

CrossRef Full Text | Google Scholar

Martinez Cesena, E. A., Capuder, T., and Mancarella, P. (2016). Flexible distributed multienergy generation system expansion planning under uncertainty. IEEE Trans. Smart Grid 7, 348–357. doi:10.1109/TSG.2015.2411392

CrossRef Full Text | Google Scholar

Mu, Y., Chen, W., Yu, X., Jia, H., Hou, K., Wang, C., et al. (2020). A double-layer planning method for integrated community energy systems with varying energy conversion efficiencies. Appl. Energy 279, 115700. doi:10.1016/j.apenergy.2020.115700

CrossRef Full Text | Google Scholar

Pineda, S., and Morales, J. M. (2019). Solving linear bilevel problems using big-Ms: Not all that glitters is gold. IEEE Trans. Power Syst. 34, 2469–2471. doi:10.1109/TPWRS.2019.2892607

CrossRef Full Text | Google Scholar

Shao, C., Shahidehpour, M., and Ding, Y. (2020). Market-based integrated generation expansion planning of electric power system and district heating systems. IEEE Trans. Sustain. Energy 11, 2483–2493. doi:10.1109/TSTE.2019.2962756

CrossRef Full Text | Google Scholar

Wang, D., Zhi, Y., Jia, H., Hou, K., Zhang, S., Du, W., et al. (2019). Optimal scheduling strategy of district integrated heat and power system with wind power and multiple energy stations considering thermal inertia of buildings under different heating regulation modes. Appl. Energy 240, 341–358. doi:10.1016/j.apenergy.2019.01.199

CrossRef Full Text | Google Scholar

Wang, L., Zhu, Z., Jiang, C., and Li, Z. (2021). Bi-level robust optimization for distribution system with multiple microgrids considering uncertainty distribution locational marginal price. IEEE Trans. Smart Grid 12, 1104–1117. doi:10.1109/TSG.2020.3037556

CrossRef Full Text | Google Scholar

Xue, Y., Li, Z., Lin, C., Guo, Q., and Sun, H. (2020). Coordinated dispatch of integrated electric and district heating systems using heterogeneous decomposition. IEEE Trans. Sustain. Energy 11, 1495–1507. doi:10.1109/TSTE.2019.2929183

CrossRef Full Text | Google Scholar

Yao, S., Gu, W., Lu, S., Zhou, S., Wu, Z., Pan, G., et al. (2021). Dynamic optimal energy flow in the heat and electricity integrated energy system. IEEE Trans. Sustain. Energy 12, 179–190. doi:10.1109/TSTE.2020.2988682

CrossRef Full Text | Google Scholar

Yin, R., Kara, E. C., Li, Y., DeForest, N., Wang, K., Yong, T., et al. (2016). Quantifying flexibility of commercial and residential loads for demand response using setpoint changes. Appl. Energy 177, 149–164. doi:10.1016/j.apenergy.2016.05.090

CrossRef Full Text | Google Scholar

Zeng, Q., Zhang, B., Fang, J., and Chen, Z. (2017). A bi-level programming for multistage co-expansion planning of the integrated gas and electricity system. Appl. Energy 200, 192–203. doi:10.1016/j.apenergy.2017.05.022

CrossRef Full Text | Google Scholar