Integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings

Abstract

Disclosed is an integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings, comprising the following steps. Step 10: respectively establish a district heating network model considering transmission delay and heat loss and a building model considering thermal storage capacity. Step 20: establish an integrated energy system optimization model consisting of a combined cooling, heat and power system model, the district heating network model and the building model. Step 30: solve the integrated energy system optimization model to obtain an optimal scheduling plan, control outputs of a gas turbine and a gas boiler per hour according to the optimal scheduling plan, and purchase electricity from a power grid and a wind power. According to the method, both the district heating network and buildings are included in a scheduling scope, so that the load adjustment with multiple degrees of freedom can be achieved.

Description

CROSS-REFERENCE TO RELATED APPLICATION

This application is a 371 of international application of PCT application serial no. PCT/CN2018/074412, filed on Jan. 29, 2018, which claims the priority benefit of China application no. 201710019950.2, filed on Jan. 11, 2017. The entirety of each of the above mentioned patent applications is hereby incorporated by reference herein and made a part of this specification

TECHNICAL FIELD

The present invention belongs to the field of combined heat and power scheduling of integrated district energy systems, and more particularly, relates to an integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings.

BACKGROUND

In recent years, wind power, as a renewable energy, has developed rapidly due to the mature technology, good economic performance, and high energy efficiency. By the end of 2015, the total installed capacity of wind power in the world reaches 423 GW, and the newly added installed capacity is 63 GW, wherein China accounts for 30.5 GW. The wind power in China is mainly developed in the three northern regions, which are rich in wind power resources and have great demand for heat load in winter. However, with the rapid development of the wind power, the wind power absorption is limited due to the strong coupling relationship between heat output and power output of the combined heat and power units. About 70% of the heat load is supplied by centralized combined heat and power (CHP) units in Jilin. The combined heat and power unit is generally operated in a mode of “thermal load following” in winter. This operation mode greatly limits the power output of the CHP unit. The heat load is high and the power load is low during the nighttime, so that the power supply exceeds the demand which causes serious wind abandonment during this period. According to the data from the National Energy Administration, in the first half of 2016, the national average wind power utilization hours are 917 h, the abandoned wind volume is 323 TWh, and the average wind power absorption rate is 21%. In short, the wind power absorption has become a key problem for the sustainable development of the wind power industry. In order to solve the problem of wind power absorption, electric power personnel have conducted a lot of researches, such as battery, electric boiler, etc. Considering the temporal and spatial relationship between the wind power and the heat load, that is, the regions and periods that have much wind power resources are generally the regions and periods with large demand for heat load, from the perspective of overall energy consumption, the district heating system can be used to provide more space for the wind power absorption.

SUMMARY

Technical Problem

The technical problem to be solved by the present invention is to provide an integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings. The method not only can utilize the transmission delay of the district heating system to coordinate the heat supply and demand balance, but also can utilize the thermal storage capacity of buildings to change heat load temporal distribution, thus improving an operational flexibility of a combined heat and power system, effectively improving a problem of wind power abandonment, and improving an overall economy of the system.

Technical Solution

In order to solve the technical problem above, the embodiments of the present invention provide an integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings, and the method comprises the following steps of:

step 10) respectively establishing a district heating network model considering transmission delay and heat loss and a building model considering thermal storage capacity;

step 20) establishing an integrated energy system optimization model consisting of a combined cooling, heat and power system model, the district heating network model, and the building model; and step 30) solving the integrated energy system optimization model to obtain an optimal scheduling plan, controlling outputs of a gas turbine and a gas boiler per hour according to the optimal scheduling plan, and purchasing electricity from a power grid and a wind power.

As a preferred embodiment, the establishing a district heating network model in the step 10) comprises:

step 101) establishing a district heating network pipeline model, which specifically comprises steps 1011) to 1015):

step 1011) establishing a nodal flow equilibrium equation, as shown in Equations (1) and (2):

$\begin{array}{cc}\sum _{k\in S_{ps,i}^e}{q}_{ps,k,t}=\sum _{k\in S_{ps,i}^s}{q}_{ps,k,t}\forall i\in S_{ns},t\in S_t& \mathrm{Equation}\left(1\right)\\ \sum _{k\in S_{\mathrm{pr},i}^e}{q}_{\mathrm{pr},k,t}=\sum _{k\in S_{\mathrm{pr},i}^s}{q}_{\mathrm{pr},k,t}\forall i\in S_{\mathrm{nr}},t\in S_t& \mathrm{Equation}\left(2\right)\end{array}$

wherein q_ps,k,t represents a water flow of a supply pipeline k at a time t in unit of kg/h; q_pr,k,t represents a water flow of a return pipeline k at the time t in unit of kg/h; S_ps,i^e represents a set of supply pipelines ended at a node i; S_pr,i^e represents a set of return pipelines ended at the node i, S_ps,i^s, S_ps,i^e represents a set of supply pipelines started at the node i, S_pr,i^s represents a set of return pipelines started at the node i, S_ns represents a set of supply pipeline nodes, S_nr represents a set of return pipeline nodes, and S_t represents a set of scheduling time periods;

step 1012) establishing a pipeline pressure loss equation, as shown in Equations (3) to (5):

$\begin{array}{cc}\Delta p_{ps,k,t}={\mu }_p·q_{ps,k,t}^2\forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(3\right)\\ \Delta p_{\mathrm{pr},k,t}={\mu }_p·q_{\mathrm{pr},k,t}^2\forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(4\right)\\ \sum _{k\in S_{ps}}\Delta p_{ps,k,t}+\sum _{k\in S_{\mathrm{pr}}}\Delta p_{ps,k,t}=\sum _{i\in S_{\mathrm{pu}}}\Delta p_{\mathrm{pu},i,t}\forall t\in S_t& \mathrm{Equation}\left(5\right)\end{array}$

wherein Δp_ps,k,t represents a pressure loss of the supply pipeline k at the time t in unit of m; μ_p represents a pressure loss factor, S_ps represents a set of supply pipelines, Δp_pr,k,t represents a pressure loss of the return pipeline k at the time t in unit of m; S_pr represents a set of return pipelines, Δp_pu,i,t represents a pressure provided by a water pump i at the time t, and S_pu represents a set of water pumps in a pipeline;

step 1013) establishing a temperature-flow-heat equation, as shown in Equations (6) and (7):Q_ps,k,tⁱⁿ=q_ps,k,t·T_ps,k,tⁱⁿ/λ∀t∈S_t,k∈S_ps Q_ps,k,t^out=q_ps,k,t·T_ps,k,t^out/λ∀t∈S_t,k∈S_ps  Equation (6)Q_pr,k,tⁱⁿ=q_pr,k,t·T_pr,k,tⁱⁿ/λ∀t∈S_t,k∈S_p Q_pr,k,t^out=q_pr,k,t·T_pr,k,t^out/λ∀t∈S_t,k∈S_p  Equation (7)

wherein Q_ps,k,tⁱⁿ represents an inlet heat power of the supply pipeline k at the time t in unit of kW; C represents a specific heat capacity of water; T_ps,k,tⁱⁿ represents an inlet temperature of the supply pipeline k at the time t in unit of ° C.; λ represents a unit conversion factor; Q_ps,k,t^out represents an outlet heat power of the supply pipeline k at the time t in unit of kW; T_ps,k,t^out represents an outlet temperature of the supply pipeline k at the time t in unit of ° C.; Q_pr,k,tⁱⁿ represents an inlet heat power of the return pipeline k at the time t in unit of kW; T_pr,k,tⁱⁿ represents an inlet temperature of the return pipeline k at the time t in unit of ° C.; Q_pr,k,t^out represents an outlet heat power of the return pipeline k at the time t in unit of kW; and T_pr,k,t^out represents an outlet temperature of the return pipeline k at the time t in unit of ° C.;

step 1014) establishing a temperature fusion equation: according to the first law of thermodynamics, if the water flow of each pipeline ended at the node i forms a stable temperature field after fusion at the node i, then inlet temperatures of the pipelines started at the node i are all equal and equal to a node temperature, as shown in Equations (8) to (11):

$\begin{array}{cc}\sum _{k\in S_{ps,i}^e}{T}_{ps,k,t}^{\mathrm{out}}·q_{ps,k,t}=T_{ns,i,t}·\sum _{k\in S_{ps,i}^e}{q}_{ps,k,t}\forall i\in S_{ns},t\in S_t& \mathrm{Equation}\left(8\right)\\ \sum _{k\in S_{\mathrm{pr},i}^e}{T}_{\mathrm{pr},k,t}^{\mathrm{out}}·q_{\mathrm{pr},k,t}=T_{\mathrm{nr},i,t}·\sum _{k\in S_{\mathrm{pr},i}^e}{q}_{\mathrm{pr},k,t}\forall i\in S_{\mathrm{nr}},t\in S_t& \mathrm{Equation}\left(9\right)\\ T_{ns,i,t}=T_{ps,k,t}^{in}\forall i\in S_{ns},t\in S_t,k\in S_{ps,i}^s& \mathrm{Equation}\left(10\right)\\ T_{\mathrm{nr},i,t}=T_{\mathrm{pr},k,t}^{in}\forall i\in S_{\mathrm{nr}},t\in S_t,k\in S_{\mathrm{pr},i}^s& \mathrm{Equation}\left(11\right)\end{array}$

wherein T_ns,i,t represents a temperature of the node i of the supply pipeline at the time t in unit of ° C.; T_nr,i,t and represents a temperature of the node i of the return pipeline at the time t in unit of ° C.; and

step 1015) establishing a district heating network transmission delay equation:

calculating a water flow rate of hot water in the pipeline, as shown in Equations (12) and (13):

$\begin{array}{cc}{v}_{ps,k,t}=\frac{q_{ps,k,t}}{\rho {\pi \left(d_k/2\right)}^2}/\lambda \forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(12\right)\\ v_{\mathrm{pr},k,t}=\frac{q_{\mathrm{pr},k,t}}{\rho {\pi \left(d_k/2\right)}^2}/\lambda \forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(13\right)\end{array}$

wherein ν_ps,k,t represents a water flow rate of hot water in the supply pipeline k at the time t in unit of m/s; ρ represents a density of hot water; d_k represents an inner diameter of the pipeline k in unit of m; and ν_pr,k,t represents a water flow rate of hot water in the return pipeline k at the time t in unit of m/s;

constraint conditions of the water flow rate of the hot water satisfying Equations (14) and (15):ν_ps,k,t^min≤ν_ps,k,t≤ν_ps,k,t^max∀k∈S_ps,t∈S_t  Equation (14)ν_pr,k,t^min≤ν_pr,k,t≤ν_pr,k,t^max∀k∈S_pr,t∈S_t  Equation (15)

wherein ν_ps,k,t^min represents a lower limit of the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_ps,k,t represents the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_ps,k,t^max represents an upper limit of the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_pr,k,t represents a lower limit of the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s; ν_pr,k,t represents the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s; and ν_pr,k,t^max represents an upper limit of the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s;

calculating a transmission time of hot water in the pipeline, as shown in Equations (16) and (17):

$\begin{array}{cc}{\tau }_{ps,k,t}=\sum _{j\in S_{ps,k}}\frac{l_j}{v_{ps,j,t}}/\lambda \forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(16\right)\\ {\tau }_{\mathrm{pr},k,t}=\sum _{j\in S_{\mathrm{pr},k}}\frac{l_j}{v_{\mathrm{pr},j,t}}/\lambda \forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(17\right)\end{array}$

wherein τ_ps,k,t represents a transmission time of the supply pipeline k at the time t in unit of h; l_j represents a length of a pipeline j in unit of m; ν_ps,j,t represents a water flow rate of hot water in a supply pipeline j at the time t in unit of m/s; S_ps,k represents a set of pipelines of hot water flowing from a heat source to the supply pipeline k; τ_pr,k,t represents a transmission time of the return pipeline k at the time t in unit of h; S_pr,k represents a set of pipelines of hot water flowing from the heat source to the return pipeline k; and ν_pr,j,t represents a water flow rate of hot water in a return pipeline j at the time t in unit of m/s;

rounding actual transmission times calculated by Equations (16) and (17), as shown in Equations (18) and (19):τ_ps,k,t^sp=round(τ_ps,k,t/Δt)∀k∈S_ps,t∈S_t  Equation (18)τ_pr,k,t^sp=round(τ_pr,k,t/Δt)∀k∈S_pr,t∈S_t  Equation (19)

wherein τ_ps,k,t^sp represents a transmission period of the supply pipeline k at the time t in unit of h; τ_pr,k,t^sp represents a transmission time period of the return pipeline k at the time t in unit of h; and Δt represents a scheduling time scale in unit of h;

after considering a transmission delay and a transmission heat loss of the heating network, inlet and outlet temperatures of the pipeline satisfying constraints shown in Equations (20) and (21):

$\begin{array}{cc}{Q}_{ps,k_1,t}^{in}-Q_{ps,k_2,t+{\tau }_{ps,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}=\left(1-{\mu }_{\mathrm{hn}}\sum _{j\in S_{ps,k_2}}{l}_j\right)·Q_{ps,k_1,t}^{in}\forall t\in S_t,k_1\in S_{ps,\mathrm{hs},}{k}_2\in S_{ps}& \mathrm{Equation}\left(20\right)\\ Q_{\mathrm{pr},k_1,t}^{in}-Q_{\mathrm{pr},k_2,t+{\tau }_{\mathrm{pr},k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}=\left(1-{\mu }_{\mathrm{hn}}\sum _{j\in S_{\mathrm{pr},k_2}}{l}_j\right)·Q_{\mathrm{pr},k_1,t}^{in}\forall t\in S_t,k_1\in S_{\mathrm{pr},m},k_2\in S_{\mathrm{pr}}& \mathrm{Equation}\left(21\right)\end{array}$

wherein Q_ps,k1,tⁱⁿ represents an inlet heat power of a supply pipeline k₁ at the time t in unit of kW;

$Q_{ps,k_2,t+{\tau }_{ps,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}$ represents an outlet heat power of a supply pipeline k₂ at a time t+τ_ps,k2,t^sp in unit of kW; μ_hn represents a heat loss rate of the heating network; S_ps,k2 represents a set of pipelines between the heat source and the supply pipeline k₂; S_ps,hs represents a set of supply pipelines connected with the heat source; Q_pr,k1,tⁱⁿ represents an inlet heat power of a return pipeline k₁ at the time t in unit of kW;

$Q_{pr,k_2,t+{\tau }_{pr,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}$ represents an outlet heat power of a return pipeline k₂ at a time t+τ_ps,k2,t^sp in unit of kW; S_pr,k2 represents a set of pipelines between the heat source and the return pipeline k₂; S_pr,m represents a set of return pipelines connected with a heat exchanger m; τ_ps,k2,t^sp represents a delay period of hot water flowing from the heat source to the supply pipeline k₂ at the time t; τ_pr,k2,t^sp represents a delay period of hot water flowing from the return pipeline k₂ to the heat source at the time t; and l_j represents a length of the pipeline j in unit of m; and

step 102) establishing a heat exchanger model:

in the heating network, coupling the heat source with a primary heat supply network by a primary heat exchanger, the model being shown in Equations (22) and (23):Q_ps,k1,tⁱⁿ−Q_pr,k2,t^out=η_ex,1·(Q_gt,t+Q_gb,t)∀t∈S_t,k₁∈S_ps,hs,k₂∈S_pr,hs  Equation (22)q_ps,k1,t=q_pr,k2,t∀t∈S_t,k₁∈S_ps,hs,k₂∈S_pr,hs  Equation (23)

wherein Q_ps,k1,tⁱⁿ represents the inlet heat power of the supply pipeline k₁ at the time t in unit of kW; Q_pr,k2,t^out represents the outlet heat power of the return pipeline k₂ at the time t in unit of kW; Q_gt,t represents heat output of the gas turbine at the time t in unit of kW; Q_gb,t represents heat output of the gas boiler at the time t in unit of kW; η_ex,1 represents a heat exchange efficiency of the primary heat exchanger; q_ps,k1,t represents a water flow of hot water in a supply pipeline k₁ at the time t in unit of kg/h; q_pr,k2,t represents a water flow of hot water in a return pipeline k₂ at the time t in unit of kg/h; and S_pr,hs represents a set of return pipelines connected with the heat source;

• • coupling the primary heat supply network with a secondary heat supply network by a secondary heat exchanger, the model being shown in Equations (24) and (25):

$\begin{array}{cc}{Q}_{ps,k_1,t}^{\mathrm{out}}-Q_{\mathrm{pr},k_2,t}^{in}=\sum _{n\in S_{\mathrm{ra},m}}{Q}_{\mathrm{ra},n,t}/{\eta }_{\mathrm{ex},2}\forall t\in S_t,k_1\in S_{ps,m},k_2\in S_{\mathrm{pr},m}& \mathrm{Equation}\left(24\right)\\ q_{ps,k_1,t}=q_{\mathrm{pr},k_2,t}\forall t\in S_t,k_1\in S_{ps,m},k_2\in S_{\mathrm{pr},m}& \mathrm{Equation}\left(25\right)\end{array}$

wherein Q_ps,k1t^out represents the inlet heat power of the supply pipeline k₁ at the time t in unit of kW; Q_pr,k2,tⁱⁿ represents the outlet heat power of the return pipeline k₂ at the time t in unit of kW; Q_ra,n,t represents a heat dissipation power of a user radiator n at the time t in unit of kW; η_ex,2 represents a heat exchange efficiency of the secondary heat exchanger; S_ps,m represents a set of supply pipelines connected with a secondary heat exchanger m; and S_pr,m represents a set of return pipelines connected with the secondary heat exchanger m.

As a preferred embodiment, Δt=0.5 h, λ=3600, C=4.168 kJ/(kg·° C.), and ρ=960 kg/m³.

As a preferred embodiment, the establishing a building model in the step 10) comprises:

step 111) establishing an indoor temperature change model, as shown in Equations (26) and (27):

$\mathrm{Equation}\left(26\right)$ $\{\begin{array}{c}{T}_{n,t+1}^{in}=\left(T_{n,t}^{\mathrm{out}}+Q_{\mathrm{ra},n,t}/{\eta }_{\mathrm{air}}\right)·\left(1-e^{-\Delta t/T_c}\right)+T_{n,t}^{in}·e^{-\Delta t/T_c}{U}_{\mathrm{she},m,t}=1\\ T_{n,t+1}^{in}=T_{n,t}^{\mathrm{out}}·\left(1-e^{-\Delta t/T_c}\right)+T_{n,t}^{in}·e^{-\Delta t/T_c}{U}_{\mathrm{she},m,t}=0\end{array}{T}_t^{\min }\le T_{n,t}^{in}\le T_t^{\max }\mathrm{Equation}\left(27\right)$

wherein T_n,t+1ⁱⁿ represents an indoor temperature of a building n at a time t+1 in unit of ° C.; T_n,t^out represents an outdoor temperature of the building n at a time t in unit of ° C.; η_air represents a thermal conductivity of air in unit of kW/° C.; T_c represents a scheduling period; T_n,tⁱⁿ represents an indoor temperature of the building n at the time t in unit of ° C.; U_she,m,t represents a switching state of a secondary heat exchanger m at the time t, i.e., U_she,m,t=1 represents that the secondary heat exchanger m is switched on at the time t, and U_she,m,t=0 represents that the secondary heat exchanger m is switched off at the time t; T_t^min represents a lower limit of an indoor temperature at the time t in unit of ° C.; and T_t^max represents an upper limit of the indoor temperature at the time t in unit of ° C.; and

step 112) calculating a heat supply index, as shown in Equations (28) and (29):

$\begin{array}{cc}{Q}_{\mathrm{res},n,t}^d=K_n^a·A_n·{10}^{-3}& \mathrm{Equation}\left(28\right)\\ \sum _t{Q}_{\mathrm{ra},n,t}=\sum _t{Q}_{\mathrm{res},n,t}^d& \mathrm{Equation}\left(29\right)\end{array}$

wherein Q_res,n,t^d represents a design heat load of the building n at the time t in unit of kW; K_n^a represents an area thermal index of the building n in unit of W/m²; and A_n represents an area of the building n in unit of m².

As a preferred embodiment, the step 20) comprises:

step 201) establishing an objective function, as shown in Equation (30):

$\begin{array}{cc}\min C_{\mathrm{total}}=C_e+C_g+C_{\mathrm{om}}+C_{\mathrm{wt}}{C}_e=\sum _t\left(K_{\mathrm{gd},e,t}·P_{\mathrm{gd},t}+K_{\mathrm{wt},e,t}·P_{\mathrm{wt},t}\right)·\Delta t{C}_g=\sum _t\frac{K_g}{H_{\mathrm{ng}}}·\left(Q_{\mathrm{gb},t}/{\eta }_{\mathrm{gb},h}+P_{\mathrm{gt},t}/{\eta }_{\mathrm{gt},e}\right)·\Delta t{C}_{\mathrm{om}}=\sum _t\left(K_{\mathrm{gt},\mathrm{om}}·P_{\mathrm{gt},t}+K_{\mathrm{gb},\mathrm{om}}·Q_{\mathrm{gb},t}\right)·\Delta t{C}_{\mathrm{wt}}=\sum _t\delta ·\left(P_{\mathrm{wt},t}^{\mathrm{pre}}-P_{\mathrm{wt},t}\right)·\Delta t& \mathrm{Equation}\left(30\right)\end{array}$

wherein C_total represents a total cost for daily operation in unit of ¥ (RMB); C_e represents an electricity purchasing cost for daily operation in unit of ¥; C_g represents a gas purchasing cost for daily operation in unit of ¥; C_om represents a maintenance cost for daily operation in unit of ¥; C_wt represents a wind power abandonment penalty for daily operation in unit of ¥; P_gd,t represents an amount of electricity purchased from the power grid at a time t in unit of kW; K_gd,e,t represents an electricity price for purchasing electricity from the power grid at the time t in unit of ¥/kWh; P_wt,t represents an amount of electricity purchased from the wind power at the time t in unit of kW; K_wt,e,t represents an electricity price for purchasing electricity from the wind power at the time t in unit of ¥/kWh; η_gb,h represents a heating efficiency of the gas boiler; P_gt,t represents an output of the gas turbine at the time t in unit of kW; η_gt,e represents a power generation efficiency of the gas turbine; K_g represents a unit price of gas in unit of ¥/m³; H_ng is a heat value of gas in unit of kWh/m³; K_gt,om, represents an operation and maintenance cost of the gas turbine in unit of ¥/kWh; K_gb,om represents an operation and maintenance cost of the gas boiler in unit of ¥/kWh; δ represents a wind power penalty cost in unit of ¥/kWh; and P_wt,t^pre represents a predicted wind power output at the time t in unit of kW; and

step 202) establishing constraint conditions, which comprises steps 2021) to 2024):

step 2021) establishing an electric power equilibrium constraint, as shown in Equation (31):

$\begin{array}{cc}{P}_{\mathrm{gt},t}+P_{\mathrm{gd},t}+P_{\mathrm{wt},t}=\sum _{m\in S_{\mathrm{she}}}\sum _{n\in S_{\mathrm{ra},m}}{P}_{n,t}& \mathrm{Equation}\left(31\right)\end{array}$

wherein S_she represents a set of secondary heat exchangers; S_ra,m represents a set of user radiators connected with a secondary heat exchanger m; and P_n,t represents an electrical load of a building n in unit of kW;

step 2022) establishing a gas turbine operation constraint, as shown in Equations (32) to (34):Q_gt,t=(1−η_gt,e−η_gt,loss)η_hr,h/η_gt,e·P_gt,t  Equation (32)P_gt,t^min≤P_gt,t≤P_gt,t^max  Equation (33)−P_gt,dw^max≤P_gt,t−P_gt,t−1≤P_gt,up^max  Equation (34)

wherein η_gt,loss represents a loss rate of the gas turbine; η_hr,h represents a recovery efficiency of a heat recover device; P_gt,t^max represents an upper limit of an operating power of the gas turbine at the time t in unit of kW; P_gt,t^min represents a lower limit of the operating power of the gas turbine at the time t in unit of kW; P_gt,dw^max represents an upper limit of a ramping down power of the gas turbine in unit of kW; P_gt,up^max represents an upper limit of a ramping up power of the gas turbine in unit of kW; and P_gt,t−1 represents an output of the gas turbine at a time t−1 in unit of kW;

step 2023) establishing a minimum start-stop time constraint: comprising a minimum operation time constraint of the gas turbine shown in Equation (35), a stop time constraint of the gas turbine shown in Equation (36), a minimum operation time constraint of the secondary heat exchanger shown in Equation (37), and a stop time constraint of the secondary heat exchanger shown in Equation (38):

$\begin{array}{cc}\{\begin{array}{c}{\tau }_{\mathrm{gt},t}^{\mathrm{on}}=\left({\tau }_{\mathrm{gt},t-1}^{\mathrm{on}}+U_{\mathrm{gt},t}·\Delta t\right)·U_{\mathrm{gt},t}\\ {\tau }_{\mathrm{gt},t}^{\mathrm{off}}=\left({\tau }_{\mathrm{gt},t-1}^{\mathrm{off}}+\left(1-U_{\mathrm{gt},t}\right)·\Delta t\right)·\left(1-U_{\mathrm{gt},t}\right)\end{array}& \mathrm{Equation}\left(35\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{gt},t}^{\mathrm{on}}\ge {\tau }_{\mathrm{gt}}^{\mathrm{on},\min }\\ {\tau }_{\mathrm{gt},t}^{\mathrm{off}}\ge {\tau }_{\mathrm{gt}}^{\mathrm{off},\min }\end{array}& \mathrm{Equation}\left(36\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{she},m,t}^{\mathrm{on}}=\left({\tau }_{\mathrm{she},m,t-1}^{\mathrm{on}}+U_{\mathrm{she},m,t}·\Delta t\right)·U_{\mathrm{she},m,t}\\ {\tau }_{\mathrm{she},m,t}^{\mathrm{off}}=\left({\tau }_{\mathrm{she},m,t-1}^{\mathrm{off}}+\left(1-U_{\mathrm{she},m,t}\right)·\Delta t\right)·\left(1-U_{\mathrm{she},m,t}\right)\end{array}& \mathrm{Equation}\left(37\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{she},m,t}^{\mathrm{on}}\ge {\tau }_{\mathrm{she},m}^{\mathrm{on},\min }\\ {\tau }_{\mathrm{she},m,t}^{\mathrm{off}}\ge {\tau }_{\mathrm{she},m}^{\mathrm{off},\min }\end{array}& \mathrm{Equation}\left(38\right)\end{array}$

wherein η_gt,t^on represents a continuous start-up time of the gas turbine at the time t in unit of h; η_gt,t−1^on represents a continuous start-up time of the gas turbine at the time t−1 in unit of h; U_gt,t represents an operating state of the gas turbine at the time t, U_gt,t=1 represents that the gas turbine is operated at the time t, and U_gt,t=0 represents that the gas turbine is shut down at the time t; τ_gt,t^off represents a continuous shutdown time of the gas turbine at the time t in unit of h; τ_gt,t−1^off represents a continuous shutdown time of the gas turbine at the time t−1 in unit of h; τ_gt^on,min represents a lower limit of the continuous start-up time of the gas turbine in unit of h; τ_gt^off,min represents a lower limit of the continuous shutdown time of the gas turbine in unit of h; τ_she,m,t^on represents a continuous start-up time of the secondary heat exchanger at the time t in unit of h; τ_she,m,t−1^on represents a continuous start-up time of the secondary heat exchanger at the time t−1 in unit of h; U_she,m,t represents an on-off state of the secondary heat exchanger m at the time t, U_she,m,t=1 represents that the secondary heat exchanger m is switched on at the time t, and U_she,m,t=0 represents that the secondary heat exchanger m is switched off at the time t; τ_she,m,t^off represents a continuous stop time of the secondary heat exchanger at the time t in unit of h; τ_she,m,t−1^off represents a continuous stop time of the secondary heat exchanger at the time t−1 in unit of h; τ_she,m^on,min represents a lower limit of the continuous start-up time of the secondary heat exchanger in unit of h; and τ_she,m^off,min represents a lower limit of the continuous stop time of the secondary heat exchanger in unit of h;

step 2024) establishing a tie-line power constraint, as shown in Equation (39):

$\begin{array}{cc}\{\begin{array}{c}{P}_{\mathrm{gd}}^{\min }\le P_{\mathrm{gd},t}\le P_{\mathrm{gd}}^{\max }\\ 0\le P_{\mathrm{wt},t}\le P_{\mathrm{wt},t}^{\mathrm{pre}}\end{array}& \mathrm{Equation}\left(39\right)\end{array}$

wherein P_gd^min represents a lower limit of purchasing electricity from the power grid in unit of kW; and P_gd^max represents an upper limit of purchasing electricity from the power grid in unit of kW.

As a preferred embodiment, in the step 30), heating network parameters are substituted into Equations (12) and (13) to obtain a water flow rate of each pipeline section; the heating network parameters and the obtained water flow rates are substituted into Equations (16) and (17) to obtain a specific delay of each pipeline section; the specific delay is substituted into Equations (18) and (19) to obtain a delay period of each pipeline section; finally, the delay period of each pipeline section, the heating network parameters and system parameters are substituted into an integrated energy system optimization model to obtain the optimal scheduling plan; the outputs of the gas turbine and the gas boiler are controlled according to the optimal scheduling plan, and the electricity is purchased from the power grid and the wind power.

Beneficial Effects

Compared with the prior art, the embodiments of the present invention have the following advantages: according to the integrated energy system optimization method considering the thermal inertia of the district heating network and buildings provided by the embodiments of the present invention, a complete district heating network model is firstly established, comprising the nodal flow equilibrium, the pressure loss equation, the node temperature fusion equation, the transmission delay equation, and other operation constraints. Secondly, it is more practical to take the buildings as a thermal energy storage unit, and it is not necessary to additionally install thermal energy storage and other equipment, which can effectively improve the economic performance. Finally, a complete operational model involving all parts of source-network-load in the integrated energy system is established. The model can realize multiple-degree-of-freedom scheduling and increase the system operational flexibility. The optimization model can change the output of devices not only by the district heating network but also by the buildings. For example, the demand-side response method is the collaborative optimization of the source and the load. However, the model in this invention is the source-network-load collaborative optimization, and has larger adjustment range and the higher system operational flexibility. The operation optimization model of the integrated energy system can greatly improve the wind power absorption rate and has the better economic performance. In view of the thermal storage properties of the buildings, the heat load distribution can be changed to increase the heat load supply during the daytime when the output of the wind power is small, so as to store the heat in the buildings, and release the heat from the building during the nighttime when the output of the wind power is large, so as to reduce the output of the CHP unit and increase the wind power absorption. In view of the time delay of the district heating network, the heat supply and demand can be balanced in a longer time scale, and the unit output is directly changed without changing the heat load distribution, thus realizing a larger time scale and a larger capacity of output adjustment.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a structure diagram of an integrated energy system according to an embodiment of the present invention.

FIG. 2 is a structure diagram of a district heating network according to an embodiment of the present invention.

FIG. 3 is a diagram of a district heating network transmission delay according to an embodiment of the present invention.

FIG. 4 is a diagram of a building model according to an embodiment of the present invention.

FIG. 5 is a structure diagram of a model according to an embodiment of the present invention.

FIG. 6 is a distribution diagram of a primary heat supply network according to an embodiment of the present invention.

DETAILED DESCRIPTION

In order to make the object, technical solution and advantages of the present invention clearer, the present invention is further described in detail below with reference to the drawings and the embodiments. It should be understood that the detailed embodiments described here are only intended to explain the present invention, rather than to limit the present invention.

Taking a combined cooling, heat and power system as an example, a structure of an integrated energy system is shown in FIG. 1. It is assumed that a district energy supply agent manages the operation of a CHP unit to satisfy heat and power load of users in the district. Both a heat supply system and a power supply system comprise three parts of source-network-load. The heat energy produced by the CHP is injected into a primary heating network through a primary heat exchanger, and then is injected into each secondary heating network through a secondary heating exchanger, and hot water in the secondary heat exchanger supplies heat to a building through each radiator. Insufficient heat load is supplemented by a gas boiler. Power transmitted by the CHP is injected into a 110 kV power transmission network through a transformer, then is injected into a 10 kV distribution network through a distribution transformer, and is finally delivered to the user. Insufficient power can be supplemented by buying electricity from a main network or a wind power.

An integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings according to an embodiment of the present invention comprises the following steps of:

step 10) respectively establishing a district heating network model considering transmission delay and a building model considering thermal storage capacity;

step 20) establishing an integrated energy system optimization model consisting of a combined cooling, heat and power system model, the district heating network model, and the building model; and

step 30) solving the integrated energy system optimization model to obtain an optimal scheduling plan, controlling outputs of a gas turbine and a gas boiler per hour according to the optimal scheduling plan, and purchasing electricity from a power grid and a wind power.

In the embodiment above, the establishing a district heating network model in the step 10) comprises:

step 101) a district heating network pipeline model is established, which specifically comprises steps 1011) to 1015):

step 1011) a nodal flow equilibrium equation is established: as a structure diagram of a district heating network pipeline shown in FIG. 2, and according to the Kirchhoff's law, a sum of flows to a certain node is equal to a sum of flows out of the node, so that the equilibrium equations shown in Equations (1) and (2) are respectively applied to the water supply and return pipelines:

$\begin{array}{cc}\sum _{k\in S_{ps,i}^e}{q}_{ps,k,t}=\sum _{k\in S_{ps,i}^s}{q}_{ps,k,t}\forall i\in S_{ns},t\in S_t& \mathrm{Equation}\left(1\right)\\ \sum _{k\in S_{\mathrm{pr},i}^e}{q}_{\mathrm{pr},k,t}=\sum _{k\in S_{\mathrm{pr},i}^s}{q}_{\mathrm{pr},k,t}\forall i\in S_{\mathrm{nr}},t\in S_t& \mathrm{Equation}\left(2\right)\end{array}$

wherein q_ps,k,t represents a water flow of a supply pipeline k at a time t in unit of kg/h; q_pr,k,t represents a water flow of a return pipeline k at the time t in unit of kg/h; S_ps,i^e represents a set of supply pipelines ended at a node i; S_pr,i^e represents a set of return pipelines ended at the node i, S_ps,i^s S_ps,i^e represents a set of supply pipelines started at the node i, S_pr,i^s represents a set of return pipelines started at the node i, S_ns represents a set of supply pipeline nodes, S_nr represents a set of return pipeline nodes, and S_t represents a set of scheduling time periods;

step 1012) a pipeline pressure loss equation is established: a pressure loss of a pipeline is proportional to a square of a water flow in the pipeline, as shown in Equations (3) and (4); and according to the Kirchhoff's law, a sum of pipeline pressure drops is equal to a sum of pressures provided by water pumps, as shown in Equation (5):

$\begin{array}{cc}\Delta p_{ps,k,t}={\mu }_p·q_{ps,k,t}^2\forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(3\right)\\ \Delta p_{\mathrm{pr},k,t}={\mu }_p·q_{\mathrm{pr},k,t}^2\forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(4\right)\\ \sum _{k\in S_{ps}}\Delta p_{ps,k,t}+\sum _{k\in S_{\mathrm{pr}}}\Delta p_{ps,k,t}=\sum _{i\in S_{\mathrm{pu}}}\Delta p_{\mathrm{pu},i,t}\forall t\in S_t& \mathrm{Equation}\left(5\right)\end{array}$

wherein Δp_ps,k,t represents a pressure loss of the supply pipeline k at the time t in unit of m; μ_p represents a pressure loss factor, S_ps represents a set of supply pipelines, Δp_pr,k,t represents a pressure loss of the return pipeline k at the time t in unit of m; S_pr represents a set of return pipelines, Δp_pu,i,t represents a pressure provided by a water pump i at the time t, and S_pu represents a set of water pumps in a pipeline;

step 1013) due to heat loss in the heating network, an inlet temperature of a pipeline is different from an outlet temperature of the pipeline. Therefore, a pipeline has two temperature variables, two heat variables and one water flow variable. A temperature-flow-heat equation is established, as shown in Equations (6) and (7):Q_ps,k,tⁱⁿ=q_ps,k,t·T_ps,k,tⁱⁿ/λ∀t∈S_t,k∈S_ps Q_ps,k,t^out=q_ps,k,t·T_ps,k,t^out/λ∀t∈S_t,k∈S_ps  Equation (6)Q_pr,k,tⁱⁿ=q_pr,k,t·T_pr,k,tⁱⁿ/λ∀t∈S_t,k∈S_p Q_pr,k,t^out=q_pr,k,t·T_pr,k,t^out/λ∀t∈S_t,k∈S_p  Equation (7)

wherein Q_ps,k,tⁱⁿ represents an inlet heat power of the supply pipeline k at the time t in unit of kW; C represents a specific heat capacity of water, and preferably, C=4.168 kJ/(kg·° C.); T_ps,k,tⁱⁿ represents an inlet temperature of the supply pipeline k at the time t in unit of ° C.; λ represents a unit conversion factor, and in the embodiment, λ is preferably 3600; Q_ps,k,t^out represents an outlet heat power of the supply pipeline k at the time t in unit of kW; T_ps,k,t^out represents an outlet temperature of the supply pipeline k at the time t in unit of ° C.; Q_pr,k,tⁱⁿ represents an inlet heat power of the return pipeline k at the time t in unit of kW; T_pr,k,tⁱⁿ represents an inlet temperature of the return pipeline k at the time t in unit of ° C.; Q_pr,k,t^out represents an outlet heat power of the return pipeline k at the time t in unit of kW; and T_pr,k,t^out represents an outlet temperature of the return pipeline k at the time t in unit of ° C.;

step 1014) a temperature fusion equation is established: as a structure diagram of a district heating network pipeline shown in FIG. 2, and according to the first law of thermodynamics, the heat flowing into a certain node is equal to the heat flowing out of the node. If the water flow of each pipeline ended at the node i forms a stable temperature field after fusion at the node i, then inlet temperatures of the pipelines started at the node i are all equal and equal to a node temperature, as shown in Equations (8) to (11):

$\begin{array}{cc}\sum _{k\in S_{ps,i}^e}{T}_{ps,k,t}^{\mathrm{out}}·q_{ps,k,t}=T_{ns,i,t}·\sum _{k\in S_{ps,i}^e}{q}_{ps,k,t}\forall i\in S_{ns},t\in S_t& \mathrm{Equation}\left(8\right)\\ \sum _{k\in S_{\mathrm{pr},i}^e}{T}_{\mathrm{pr},k,t}^{\mathrm{out}}·q_{\mathrm{pr},k,t}=T_{\mathrm{nr},i,t}·\sum _{k\in S_{\mathrm{pr},i}^e}{q}_{\mathrm{pr},k,t}\forall i\in S_{\mathrm{nr}},t\in S_t& \mathrm{Equation}\left(9\right)\\ T_{ns,i,t}=T_{ps,k,t}^{in}\forall i\in S_{ns},t\in S_t,k\in S_{ps,i}^s& \mathrm{Equation}\left(10\right)\\ T_{\mathrm{nr},i,t}=T_{\mathrm{pr},k,t}^{in}\forall i\in S_{\mathrm{nr}},t\in S_t,k\in S_{\mathrm{pr},i}^s& \mathrm{Equation}\left(11\right)\end{array}$

wherein T_ns,i,t represents a temperature of the node i of the supply pipeline at the time t in unit of ° C.; and T_nr,i,t represents a temperature of the node i of the return pipeline at the time t in unit of ° C.; and

step 1015) a district heating network transmission delay equation is established:

as a diagram of a district heating network transmission delay shown in FIG. 3, a water flow rate of hot water in the pipeline is calculated, and the water flow rate of the hot water in the pipeline is proportional to the water flow of the pipeline, as shown in Equations (12) and (13):

$\begin{array}{cc}{v}_{ps,k,t}=\frac{q_{ps,k,t}}{\rho {\pi \left(d_k/2\right)}^2}/\lambda \forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(12\right)\\ v_{\mathrm{pr},k,t}=\frac{q_{\mathrm{pr},k,t}}{\rho {\pi \left(d_k/2\right)}^2}/\lambda \forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(13\right)\end{array}$

wherein ν_ps,k,t represents a water flow rate of hot water in the supply pipeline k at the time t in unit of m/s; ρ represents a density of hot water, and in the embodiment, ρ=960 kg/m³; d_k represents an inner diameter of the pipeline k in unit of m; and ν_pr,k,t represents a water flow rate of hot water in the return pipeline k at the time t in unit of m/s;

if the water flows too fast, the pipeline can be in an unstable hydraulic condition, and if the water flows too slow, a heating effect can be affected. Constraint conditions of the water flow rate of the hot water are established, satisfying Equations (14) and (15):ν_ps,k,t^min≤ν_ps,k,t≤ν_ps,k,t^max∀k∈S_ps,t∈S_t  Equation (14)ν_pr,k,t^min≤ν_pr,k,t≤ν_pr,k,t^max∀k∈S_pr,t∈S_t  Equation (15)

wherein ν_ps,k,t^min represents a lower limit of the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_ps,k,t represents the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_ps,k,t^max represents an upper limit of the water flow rate of the hot water in the supply pipeline k at the time t in unit of m/s; ν_pr,k,t represents a lower limit of the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s; ν_pr,k,t represents the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s; and ν_pr,k,t^max represents an upper limit of the water flow rate of the hot water in the return pipeline k at the time t in unit of m/s;

a transmission time of hot water in the pipeline is calculated, as shown in Equations (16) and (17):

$\begin{array}{cc}{\tau }_{ps,k,t}=\sum _{j\in S_{ps,k}}\frac{l_j}{v_{ps,j,t}}/\lambda \forall k\in S_{ps},t\in S_t& \mathrm{Equation}\left(16\right)\\ {\tau }_{\mathrm{pr},k,t}=\sum _{j\in S_{\mathrm{pr},k}}\frac{l_j}{v_{\mathrm{pr},j,t}}/\lambda \forall k\in S_{\mathrm{pr}},t\in S_t& \mathrm{Equation}\left(17\right)\end{array}$

wherein τ_ps,k,t represents a transmission time of the supply pipeline k at the time t in unit of h; l_j represents a length of a pipeline j in unit of m; ν_ps,j,t represents a water flow rate of hot water in a supply pipeline j at the time t in unit of m/s; S_ps,k represents a set of pipelines of hot water flowing from a heat source to the supply pipeline k; τ_pr,k,t represents a transmission time of the return pipeline k at the time t in unit of h; S_pr,k represents a set of pipelines of hot water flowing from the heat source to the return pipeline k; and ν_pr,j,t represents a water flow rate of hot water in a return pipeline j at the time t in unit of m/s;

since actual transmission times are calculated by Equations (16) and (17), and a scheduling command is executed in unit of integer time period in the scheduling optimization model, the actual transmission times calculated by Equations (16) and (17) are rounded, as shown in Equations (18) and (19):τ_ps,k,t^sp=round(τ_ps,k,t/Δt)∀k∈S_ps,t∈S_t  Equation (18)τ_pr,k,t^sp=round(τ_pr,k,t/Δt)∀k∈S_pr,t∈S_t  Equation (19)

wherein τ_ps,k,t^sp represents a transmission period of the supply pipeline k at the time t in unit of h; τ_pr,k,t^sp represents a transmission time period of the return pipeline k at the time t in unit of h; and Δt represents a scheduling time scale in unit of h; and preferably, Δt=0.5.

As a diagram of a district heating network transmission delay shown in FIG. 4, after considering a transmission delay and a transmission heat loss of the heating network, inlet and outlet temperatures of the pipeline satisfy constraints shown in Equations (20) and (21):

$\begin{array}{cc}{Q}_{ps,k_1,t}^{in}-Q_{ps,k_2,t+{\tau }_{ps,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}=\left(1-{\mu }_{\mathrm{hn}}\sum _{j\in S_{ps,k_2}}{l}_j\right)·Q_{ps,k_1,t}^{in}\forall t\in S_t,k_1\in S_{ps,\mathrm{hs}},k_2\in S_{ps}& \mathrm{Equation}\left(20\right)\\ Q_{\mathrm{pr},k_1,t}^{in}-Q_{\mathrm{pr},k_2,t+{\tau }_{\mathrm{pr},k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}=\left(1-{\mu }_{\mathrm{hn}}\sum _{j\in S_{\mathrm{pr},k_2}}{l}_j\right)·Q_{\mathrm{pr},k_1,t}^{in}\forall t\in S_t,k_1\in S_{\mathrm{pr},m},k_2\in S_{\mathrm{pr}}& \mathrm{Equation}\left(21\right)\end{array}$

wherein Q_ps,k1,tⁱⁿ represents an inlet heat power of a supply pipeline k₁ at the time t in unit of kW;

$Q_{ps,k_2,t+{\tau }_{ps,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}$ represents an outlet heat power of a supply pipeline k₂ at a time t+τ_ps,k2,t^sp in unit of kW; μ_hn represents a heat loss rate of the heating network; S_ps,k2 represents a set of pipelines between the heat source and the supply pipeline k₂; S_ps,hs represents a set of supply pipelines connected with the heat source; Q_pr,k1,tⁱⁿ represents an inlet heat power of a return pipeline k₁ at the time t in unit of kW;

$Q_{pr,k_2,t+{\tau }_{pr,k_2,t}^{\mathrm{sp}}}^{\mathrm{out}}$ represents an outlet heat power of a return pipeline k₂ at a time t+τ_ps,k2,t^sp in unit of kW; S_pr,k2 represents a set of pipelines between the heat source and the return pipeline k₂; S_pr,m represents a set of return pipelines connected with a heat exchanger m; τ_ps,k2,t^sp represents a delay period of hot water flowing from the heat source to the supply pipeline k₂ at the time t; τ_pr,k2,t^sp represents a delay period of hot water flowing from the return pipeline k₂ to the heat source at the time t; and l_j represents a length of the pipeline j in unit of m; and

step 102) a heat exchanger model is established:

in the heating network, the heat generated by the CHP is injected into a primary heat supply network through coupling the heat source with the primary heat supply network by a primary heat exchanger, and the model is shown in Equation (22):Q_ps,k1,tⁱⁿ−Q_pr,k2,t^out=η_ex,1·(Q_gt,t+Q_gb,t)∀t∈S_t,k₁∈S_ps,hs,k₂∈S_pr,hs  Equation (22)

the water flow equilibrium of the water supply and return pipelines in the primary heat exchanger is ensured, as shown in Equation (23):q_ps,k1,t=q_pr,k2,t∀t∈S_t,k₁∈S_ps,hs,k₂∈S_pr,hs  Equation (23)

wherein Q_ps,k1,tⁱⁿ represents the inlet heat power of the supply pipeline k₁ at the time t in unit of kW; Q_pr,k2,t^out represents the outlet heat power of the return pipeline k₂ at the time t in unit of kW; Q_gt,t represents heat output of the gas turbine at the time t in unit of kW; Q_gb,t represents heat output of the gas boiler at the time t in unit of kW; η_ex,1 represents a heat exchange efficiency of the primary heat exchanger; q_ps,k1,t represents a water flow of hot water in a supply pipeline k₁ at the time t in unit of kg/h; q_pr,k2,t represents a water flow of hot water in a return pipeline k₂ at the time t in unit of kg/h; and S_pr,hs represents a set of return pipelines connected with the heat source;

the primary heat supply network is coupled with a secondary heat supply network by a secondary heat exchanger, and the model is shown in Equation (24):

$\begin{array}{cc}{Q}_{ps,k_1,t}^{\mathrm{out}}-Q_{\mathrm{pr},k_2,t}^{in}=\sum _{n\in S_{\mathrm{ra},m}}{Q}_{\mathrm{ra},n,t}/{\eta }_{\mathrm{ex},2}\forall t\in S_t,k_1\in S_{ps,m},k_2\in S_{\mathrm{pr},m}& \mathrm{Equation}\left(24\right)\end{array}$

the water flow equilibrium of the water supply and return pipelines in the secondary heat exchanger is ensured, as shown in Equation (25):q_ps,k1,t=q_pr,k2,t∀t∈S_t,k₁∈S_ps,m,k₂∈S_pr,m  Equation (25)

wherein Q_pr,k1,t^out represents the inlet heat power of the supply pipeline k₁ at the time t in unit of kW; Q_pr,k2,tⁱⁿ represents the outlet heat power of the return pipeline k₂ at the time t in unit of kW; Q_ra,n,t represents a heat dissipation power of a user radiator n at the time t in unit of kW; η_ex,2 represents a heat exchange efficiency of the secondary heat exchanger; S_ps,m represents a set of supply pipelines connected with a secondary heat exchanger m; and S_pr,m represents a set of return pipelines connected with the secondary heat exchanger m.

In the embodiment above, the establishing a building model in the step 10) comprises:

step 111) an indoor temperature change model is established.

Considering that residential heating is a main component of the heat load and has huge adjustable potential, the heat load in the embodiment is the heat load of the residential heating. As a diagram of a building model shown in FIG. 4, the heat load of residents in a small district can be adjusted by adjusting the switching of the secondary heat exchanger. Assuming that N_m users are arranged under the secondary heat exchanger m, an indoor temperature change of a user n can be expressed as Equation (26), and in order to ensure a heating comfort level of residents, the temperature is required to satisfy Equation (27):

Equations (26) and (27) are as follows:

$\mathrm{Equation}\left(26\right)$ $\{\begin{array}{c}{T}_{n,t+1}^{in}=\left(T_{n,t}^{\mathrm{out}}+Q_{\mathrm{ra},n,t}/{\eta }_{\mathrm{air}}\right)·\left(1-e^{-\Delta t/T_c}\right)+T_{n,t}^{in}·e^{-\Delta t/T_c}{U}_{\mathrm{she},m,t}=1\\ T_{n,t+1}^{in}=T_{n,t}^{\mathrm{out}}·\left(1-e^{-\Delta t/T_c}\right)+T_{n,t}^{in}·e^{-\Delta t/T_c}{U}_{\mathrm{she},m,t}=0\end{array}{T}_t^{\min }\le T_{n,t}^{in}\le T_t^{\max }\mathrm{Equation}\left(27\right)$

wherein T_n,t+1ⁱⁿ represents an indoor temperature of a building n at a time t+1 in unit of ° C.; T_n,t^out represents an outdoor temperature of the building n at a time t in unit of ° C.; η_air represents a thermal conductivity of air in unit of kW/° C.; T_c represents a scheduling period, and preferably, T_c=24 h; T_n,tⁱⁿ represents an indoor temperature of the building n at the time t in unit of t; U_she,m,t represents a switching state of a secondary heat exchanger m at the time t, i.e., U_she,m,t=1 represents that the secondary heat exchanger m is switched on at the time t, and U_she,m,t=0 represents that the secondary heat exchanger m is switched off at the time t; T_t^min represents a lower limit of an indoor temperature at the time t in unit of ° C.; and T_t^max represents an upper limit of the indoor temperature at the time t in unit of ° C.; and

step 112) a heat supply index is calculated: Equation (28) is an algorithm of an area thermal index of a design heat load of residential heating; in order to ensure a residential heating quality, the total heat supply is required to be equal to the total design heat load while adjusting the heat load of residents, as shown in Equation (29):

Equations (28) and (29) are as follows:

$\begin{array}{cc}{Q}_{\mathrm{res},n,t}^d=K_n^a·A_n·{10}^{-3}& \mathrm{Equation}\left(28\right)\\ \sum _t{Q}_{\mathrm{ra},n,t}=\sum _t{Q}_{\mathrm{res},n,t}^d& \mathrm{Equation}\left(29\right)\end{array}$

wherein Q_res,n,t^d represents a design heat load of the building n at the time t in unit of kW; K_n^a represents an area thermal index of the building n in unit of W/m²; and A_n represents an area of the building n in unit of m².

In the embodiment above, the step 20) comprises the following steps. As a structure diagram of an optimization model shown in FIG. 5, a district energy supply agent establishes an optimization model to solve a day-ahead scheduling plan of the unit and a power supply plan to achieve optimal economy. Since a transmission delay variable in the district heating network model above is a time-related variable and needs to be superimposed with a time variable, the transmission delay variable cannot be solved by general commercial optimization software, and it is considered that in actual district heating network operation, a quality adjustment mode is generally adopted, i.e., a temperature of the hot water is adjusted to ensure a fixed flow, and the fixed water flow is generally a design flow. Therefore, when the design water flow of the pipeline is determined, the transmission delay of each section of the pipeline is also determined, which can be used as a parameter input model, and the model can be solved by the general commercial software. The design water flow of the pipeline can be known by inquiring a design code (GB50019-2003) for a heating and ventilation air conditioner.

In step 201), an objective function is established, and the optimization model aims at economy, comprising an electricity purchasing cost, an operation and maintenance cost, a gas purchasing cost and a wind power abandonment penalty, as shown in Equation (30):

$\begin{array}{cc}\min C_{\mathrm{total}}=C_e+C_g+C_{\mathrm{om}}+C_{\mathrm{wt}}{C}_e=\sum _t\left(K_{\mathrm{gd},e,t}·P_{\mathrm{gd},t}+K_{\mathrm{wt},e,t}·P_{\mathrm{wt},t}\right)·\Delta t{C}_g=\sum _t\frac{K_g}{H_{\mathrm{ng}}}·\left(Q_{\mathrm{gb},t}/{\eta }_{\mathrm{gb},h}+P_{\mathrm{gt},t}/{\eta }_{\mathrm{gt},e}\right)·\Delta t{C}_{\mathrm{om}}=\sum _t\left(K_{\mathrm{gt},\mathrm{om}}·P_{\mathrm{gt},t}+K_{\mathrm{gb},\mathrm{om}}·Q_{\mathrm{gb},t}\right)·\Delta t{C}_{\mathrm{wt}}=\sum _t\delta ·\left(P_{\mathrm{wt},t}^{\mathrm{pre}}-P_{\mathrm{wt},t}\right)·\Delta t& \mathrm{Equation}\left(30\right)\end{array}$

wherein C_total represents a total cost for daily operation in unit of ¥; C_e represents an electricity purchasing cost for daily operation in unit of ¥; C_g represents a gas purchasing cost for daily operation in unit of ¥; C_om, represents a maintenance cost for daily operation in unit of ¥; C_wt represents a wind power abandonment penalty for daily operation in unit of ¥; P_gd,t represents an amount of electricity purchased from the power grid at a time t in unit of kW; K_gd,e,t represents an electricity price for purchasing electricity from the power grid at the time t in unit of ¥/kWh; P_wt,t represents an amount of electricity purchased from the wind power at the time t in unit of kW; K_wt,e,t represents an electricity price for purchasing electricity from the wind power at the time t in unit of ¥/kWh; η_gb,h represents a heating efficiency of the gas boiler; P_gt,t represents an output of the gas turbine at the time t in unit of kW; η_gt,e represents a power generation efficiency of the gas turbine; K_g represents a unit price of gas in unit of ¥/m³; H_ng is a heat value of gas in unit of kWh/m³; K_gt,om represents an operation and maintenance cost of the gas turbine in unit of ¥/kWh; K_gb,om represents an operation and maintenance cost of the gas boiler in unit of ¥/kWh; δ represents a wind power penalty cost in unit of ¥/kWh; and P_wt,t^pre represents a predicted wind power output at the time t in unit of kW.

In step 202), constraint conditions are established, which comprises steps 2021) to 2024).

In step 2021), an electric power equilibrium constraint is established, as shown in Equation (31):

$\begin{array}{cc}{P}_{\mathrm{gt},t}+P_{\mathrm{gd},t}+P_{\mathrm{wt},t}=\sum _{m\in S_{\mathrm{she}}}\sum _{n\in S_{\mathrm{ra},m}}{P}_{n,t}& \mathrm{Equation}\left(31\right)\end{array}$

wherein S_she represents a set of secondary heat exchangers; S_ra,m represents a set of user radiators connected with a secondary heat exchanger m; and P_n,t represents an electrical load of a building n in unit of kW;

step 2022) a gas turbine operation constraint is established, as shown in Equations (32) to (34):Q_gt,t=(1−η_gt,e−η_gt,loss)η_hr,h/η_gt,e·P_gt,t  Equation (32)P_gt,t^min≤P_gt,t≤P_gt,t^max  Equation (33)

wherein η_gt,loss represents a loss rate of the gas turbine; η_hr,h represents a recovery efficiency of a heat recover device; P_gt,t^max represents an upper limit of an operating power of the gas turbine at the time t in unit of kW; P_gt,t^min represents a lower limit of the operating power of the gas turbine at the time t in unit of kW; P_gt,dw^max represents an upper limit of a ramping down power of the gas turbine in unit of kW; P_gt,up^max represents an upper limit of a ramping up power of the gas turbine in unit of kW; and P_gt,t−1 represents an output of the gas turbine at a time t−1 in unit of kW.

In step 2023), a minimum start-stop time constraint is established: in order to prevent mechanical loss on the gas turbine and the secondary heat exchanger caused by frequent starting and stopping, it is necessary to limit a minimum operation and a stop time of the gas turbine established, and a minimum operation and a stop time of the secondary heat exchanger, specifically comprising a minimum operation time constraint of the gas turbine shown in Equation (35), a stop time constraint of the gas turbine shown in Equation (36), a minimum operation time constraint of the secondary heat exchanger shown in Equation (37), and a stop time constraint of the secondary heat exchanger shown in Equation (38):

$\begin{array}{cc}\{\begin{array}{c}{\tau }_{\mathrm{gt},t}^{\mathrm{on}}=\left({\tau }_{\mathrm{gt},t-1}^{\mathrm{on}}+U_{\mathrm{gt},t}·\Delta t\right)·U_{\mathrm{gt},t}\\ {\tau }_{\mathrm{gt},t}^{\mathrm{off}}=\left({\tau }_{\mathrm{gt},t-1}^{\mathrm{off}}+\left(1-U_{\mathrm{gt},t}\right)·\Delta t\right)·\left(1-U_{\mathrm{gt},t}\right)\end{array}& \mathrm{Equation}\left(35\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{gt},t}^{\mathrm{on}}\ge {\tau }_{\mathrm{gt}}^{\mathrm{on},\min }\\ {\tau }_{\mathrm{gt},t}^{\mathrm{off}}\ge {\tau }_{\mathrm{gt}}^{\mathrm{off},\min }\end{array}& \mathrm{Equation}\left(36\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{she},m,t}^{\mathrm{on}}=\left({\tau }_{\mathrm{she},m,t-1}^{\mathrm{on}}+U_{\mathrm{she},m,t}·\Delta t\right)·U_{\mathrm{she},m,t}\\ {\tau }_{\mathrm{she},m,t}^{\mathrm{off}}=\left({\tau }_{\mathrm{she},m,t-1}^{\mathrm{off}}+\left(1-U_{\mathrm{she},m,t}\right)·\Delta t\right)·\left(1-U_{\mathrm{she},m,t}\right)\end{array}& \mathrm{Equation}\left(37\right)\\ \{\begin{array}{c}{\tau }_{\mathrm{she},m,t}^{\mathrm{on}}\ge {\tau }_{\mathrm{she},m}^{\mathrm{on},\min }\\ {\tau }_{\mathrm{she},m,t}^{\mathrm{off}}\ge {\tau }_{\mathrm{she},m}^{\mathrm{off},\min }\end{array}& \mathrm{Equation}\left(38\right)\end{array}$

wherein τ_gt,t^on represents a continuous start-up time of the gas turbine at the time t in unit of h; τ_gt,t−1^on represents a continuous start-up time of the gas turbine at the time t−1 in unit of h; U_gt,t represents an operating state of the gas turbine at the time t, U_gt,t=1 represents that the gas turbine is operated at the time t, and U_gt,t=0 represents that the gas turbine is shut down at the time t; τ_gt,t^off represents a continuous shutdown time of the gas turbine at the time t in unit of h; τ_gt,t−1^off represents a continuous shutdown time of the gas turbine at the time t−1 in unit of h; τ_gt^on,min represents a lower limit of the continuous start-up time of the gas turbine in unit of h; τ_gt^off,min represents a lower limit of the continuous shutdown time of the gas turbine in unit of h; τ_she,m,t^on represents a continuous start-up time of the secondary heat exchanger at the time t in unit of h; τ_she,m,t−1^on represents a continuous start-up time of the secondary heat exchanger at the time t−1 in unit of h; U_she,m,t represents an on-off state of the secondary heat exchanger m at the time t, U_she,m,t=1 represents that the secondary heat exchanger m is switched on at the time t, and U_she,m,t=0 represents that the secondary heat exchanger m is switched off at the time t; τ_she,m,t^off represents a continuous stop time of the secondary heat exchanger at the time t in unit of h; τ_she,m,t−1^off represents a continuous stop time of the secondary heat exchanger at the time t−1 in unit of h; τ_she,m^on,min represents a lower limit of the continuous start-up time of the secondary heat exchanger in unit of h; and τ_she,m^off,min represents a lower limit of the continuous stop time of the secondary heat exchanger in unit of h.

In step 2024), a tie-line power constraint is established, as shown in Equation (39):

$\begin{array}{cc}\{\begin{array}{c}{P}_{\mathrm{gd}}^{\min }\le P_{\mathrm{gd},t}\le P_{\mathrm{gd}}^{\max }\\ 0\le P_{\mathrm{wt},t}\le P_{\mathrm{wt},t}^{\mathrm{pre}}\end{array}& \mathrm{Equation}\left(39\right)\end{array}$

wherein P_gd^min represents a lower limit of purchasing electricity from the power grid in unit of kW; and P_gd^max represents an upper limit of purchasing electricity from the power grid in unit of kW.

In the embodiment above, in the step 30), heating network parameters are substituted into Equations (12) and (13) to obtain a water flow rate of each pipeline section; the heating network parameters and the obtained water flow rates are substituted into Equations (16) and (17) to obtain a specific delay of each pipeline section; the specific delay is substituted into Equations (18) and (19) to obtain a delay period of each pipeline section; finally, the delay period of each pipeline section, the heating network parameters and system parameters are substituted into an integrated energy system optimization model to obtain the optimal scheduling plan; the outputs of the gas turbine and the gas boiler are controlled according to the optimal scheduling plan, and the electricity is purchased from the power grid and the wind power.

In the embodiment, the objective function and the constraint of the operation optimization model are both linear, so that the optimization operation model of the integrated energy system established according to the embodiment of the present invention is a typical mixed integer linear programming model. In the embodiment, the variables to be optimized comprise a day-ahead output plan of the unit, an electricity purchase plan, a temperature of a district heating network operating pipeline, and a control state of the secondary heat exchanger.

In the embodiment, the district heating network model containing the constraints of the nodal flow equilibrium, the node temperature fusion, the transmission delay, the transmission heat loss and the like is firstly established, and the operation optimization model of the integrated energy system is secondly established by taking a building as a heat storage unit and combining the district heating network model. In the embodiment of the present invention, the heating network and the user are included in scheduling, thermoelectric coupling is utilized, and the wind power absorption is promoted from the perspective of a thermodynamic system, so that the wind power absorption can be greatly improved, and the system operation cost can be effectively reduced. In view of the thermal storage capacity of the building, the heat load distribution is changed to change the output distribution of the CHP unit to realize more wind power absorption during the nighttime; and in the operation optimization model, in view of a district heating network delay, the output of the CHP unit and the heat load of the user are staggered to form a leading supply delay on the time scale, thus improving the wind power absorption and the system operation economy.

An embodiment is illustrated as follows. Taking an actual heating district in Jilin as an example, layout of the primary heating network in the district is shown in FIG. 6, with 50 pipelines, 24 nodes and 26 secondary heat exchangers in total. In addition to the district energy supply agent (DESA), an independent heat source is provided in the district for standby during peak heating, and V1/V2 is a peak-load regulating stop valve. A total installed capacity of wind power in Jilin is 5000 MW, which is absorbed by many districts, and in order to facilitate the research on the wind power absorption capacity of a single district, this case sets a 5 MW virtual wind turbine to be connected through a 10 kV distribution grid. The gas price in this district is 2.3 V/m³, the electricity price is 5.25 ¥/kWh, and the wind power penalty is 0.2 ¥/kWh. For the convenience of comparative analysis: Case 1 is set as a basic case, and the transmission delay of the district heating network and the thermal energy storage properties of the buildings are not considered in this case; only the heat storage properties of the house are considered in Case 2; and Case 3 is the model proposed in the embodiment of the present invention. Operational results of the cases are shown in Table 1.

	TABLE 1
	Transmission			Wind power absorption
	delay of	Thermal	Economy (*10{circumflex over ( )}4¥)	(MWh)
	district heating	storage of		Electricity		Wind	Total	Absorption	Absorption
Case	network	house	Gas	purchasing	Maintenance	abandonment	costs	amount	rate
1	x	x	11.20	1.60	1.73	1.13	15.66	30.51	51.6%
2	x	✓	10.87	2.67	1.40	0.33	15.27	51.40	86.9%
3	✓	✓	9.87	2.93	1.33	0.07	14.20	56.93	96.2%

1) Economy Analysis

It can be seen in Table 1 that the daily operation cost of Case 2 is saved by ¥ 3900 compared with Case 1, with a saving rate of 2.5%. In Case 2, the thermal storage of buildings is considered, i.e., a thermal storage device is introduced. By readjusting the heat load distribution, more wind power is absorbed, which increases the electricity purchasing cost, but simultaneously reduces the wind abandonment penalty. In addition, due to the introduction of the thermal storage device, the output is increased during the daytime when the power load is high, and the heat is stored in the buildings, while the output is reduced during the nighttime when the power load is low, and the heat stored in the buildings is released, so that the unit output is optimized, the gas and operational maintenance costs are reduced, and the economy is improved. In Case 3, the daily operation cost is saved by ¥ 10700, with a saving rate of 7.0% compared with Case 2, and is saved by ¥14600, with a saving rate of 9.3%, compared with Case 1. Thus, it can be seen that, considering the transmission delay of the district heating network, a larger space can be provided for the wind power absorption.

2) Wind Power Absorption Analysis

A wind power absorption amount is 30.51 MWh in Case 1, with an absorption rate of 51.6%. A wind power absorption amount is 51.04 MWh in Case 2, with a wind power absorption rate increased by 35.3% compared with Case 1. A wind power absorption amount is 56.93 MWh in Case 3, with an absorption rate of 96.2%, which is increased by 44.6% compared with Case 1 and is increased by 9.3% compared with Case 2. During the daytime (07:00-21:00), the power load is high but the heat load is low, and a fan output is low, and the wind power during this period can be absorbed in the three cases. During the nighttime (21:00-24:00) and the early morning (00:00-07:00), the power load is low and the fan output is large, so that serious wind power abandonment occurs during this period in Case 1; and in Case 2, due to the consideration of the heat storage properties of the building and the district heating network delay, the wind power absorption during this period is greatly improved. The foregoing is merely the preferred embodiments of the present invention, and it should be noted that those of ordinary skills in the art may further make a plurality of improvements and decorations without departing from the principle of the present invention, and these improvements and decorations shall also fall within the protection scope of the present invention.

Claims

1. An integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings, the integrated energy system comprises a heat supply system and a power supply system, the method comprising the following steps: step 10: respectively establishing a district heating network model considering transmission delay and heat loss and a building model considering thermal storage capacity;step 20: establishing an integrated energy system optimization model consisting of a combined cooling, heat and power system model, the district heating network model, and the building model; andstep 30: solving the integrated energy system optimization model to obtain an optimal scheduling plan for operating the integrated energy system, controlling outputs of a gas turbine and a gas boiler of the heat supply system per hour according to the optimal scheduling plan, and purchasing electricity from a power grid and a wind power, for the power supply system according to the optimal scheduling plan,wherein the establishing a district heating network model in the step 10 comprises:step 101: establishing a district heating network pipeline model, which specifically comprises steps 1011 to 1015:step 1011: establishing a nodal flow equilibrium equation, as shown in Equations (1) and (2):

2. The integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings according to claim 1, wherein Δt=0.5 h, λ=3600, C=4.168 kJ/(kg·° C.), and ρ=960 kg/m3.

3. The integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings according to claim 1, wherein the establishing a building model in the step 10 comprises: step 111: establishing an indoor temperature change model, as shown in Equations (26) and (27):

4. The integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings according to claim 1, wherein the step 20 comprises: step 201: establishing an objective function, as shown in Equation (30):

5. The integrated energy system operational optimization method considering thermal inertia of district heating networks and buildings according to claim 1, wherein in the step 30, heating network parameters are substituted into Equations (12) and (13) to obtain a water flow rate of each pipeline section; the heating network parameters and the obtained water flow rates are substituted into Equations (16) and (17) to obtain a specific delay of each pipeline section; the specific delay is substituted into Equations (18) and (19) to obtain a delay period of each pipeline section; finally, the delay period of each pipeline section, the heating network parameters and system parameters are substituted into an integrated energy system optimization model to obtain the optimal scheduling plan; the outputs of the gas turbine and the gas boiler are controlled according to the optimal scheduling plan, and the electricity is purchased from the power grid and the wind power.