# Research on Cooperative Optimization Control Strategy of Multitype Energy Sources

Authors: Renliang Liu; Yong Chen; Zhe Tan; Chao Liu; Kun Yang; Liang Pei
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

With the increasing consumption of energy, the efficiency of energy utilization urgently needs to be further improved. At present, the coordinated optimization control of the integrated energy system of multiple types of cooling, heating, and power equipment is an important way to improve the comprehensive energy efficiency of the regional power grid and reduce the operating cost of the power grid. Aiming at this scenario, this paper establishes a fine energy storage model by analyzing the uncertainty of wind power output and considering the influence of low temperature and other conditions on the energy storage device in the energy storage side. On the load side, the influence of comprehensive demand response of electricity and heat on system operation is analyzed, and a combined cold, heat, and electricity supply system including renewable energy and energy storage device is established. Aiming at the optimal total cost of the integrated energy cooling and heating triple power system collaborative optimization control, the crow search algorithm is used to iteratively optimize the configuration model of the triple power system. The four models are considered in the article include groundless heat pump and energy storage, excluding joint demand response, including ground source heat pump and traditional energy storage model, excluding joint demand response, including ground source heat pump and fine energy storage model, excluding joint demand response, including ground source heat pump and fine energy storage model, taking into account the joint demand response. The simulation results show that the cooperative optimization control strategy of the combined cooling, heating, and power system with renewable energy and fine energy storage device model can enhance the system’s schedulable space, improve the comprehensive energy utilization efficiency, and have considerable economic benefits.


---

## Body

## 1. Introduction

Energy, as the most basic driving force for the growth and development of the national economy, is the basis for human survival, and it is also the current research focus of the international community. In the next few years, the production and utilization of energy will inevitably accelerate the consumption of fossil fuels, thus causing inevitable damage to our environment [1]. In recent years, in order to cope with the challenges of limited fossil fuels and global environmental issues, distributed power generation technology has developed rapidly around the world [2]. At the same time, the integrated energy system has developed rapidly at home and abroad due to its characteristics such as the ability to absorb new energy. In this regard, because the combined cooling heating and power (CCHP) technology can not only meet the demand for refrigeration but also meet the demand for cogeneration, and the economy is getting higher and higher. Therefore, it is more and more favored by people and is the main direction of future energy technology development.The integrated energy system has the characteristics of multienergy flow, strong coupling, and strong uncertainty. Therefore, an efficient coordinated control strategy is very important to improve the performance of the cooling, heating, and power integrated energy system. In terms of control strategy, a model predictive control MPC controller based on random scenes combined with a nonlinear Modelica model provides a richer architectural description and captures architectural dynamics more accurately than linear models [3]. Literature [4] proposed a multitime-scale energy management framework for smart photovoltaic (PV) systems, which can calculate optimal plans for battery operation, power purchase, and application usage. A smart photovoltaic system is a local energy community that includes multiple buildings and homes equipped with photovoltaic panels and batteries. In terms of model establishment, literature [5] developed a two-stage stochastic mixed integer linear programming model to determine the key operating decisions of the microgrid system when the power demand is uncertain. Its purpose is to minimize the overall system cost of the building under investigation within the specified time frame. Literature [6] always takes into account the uncertainties related to wind power output, power, heating, and cooling loads and proposes a hierarchical two-stage robust optimal scheduling model for multiple CCHP microgrid systems, which helps to minimize the operating cost of CCHP microgrid and transmission line failure. An integrated community energy system (ICES) scheduling model is proposed [7], which includes a variety of renewable energy and energy conversion devices, which can be used to develop an integrated energy system with a high penetration rate of renewable energy in the future. A new model and optimized dispatching method for combined heat and power (P2G) and carbon capture system (CCS), which solves the carbon source required by P2G and the carbon emissions of combined heat and power through optimized dispatch in the integrated energy system (IES) question [8]. Literature [9] proposes a two-stage robust scheduling model of MES that considers integrated thermal and power demand response (DR) to drive robust operation decision-making. A new type of CCHP system [10] combining ORC and ST, using solar energy and waste heat as the ORC driving heat source, and a collaborative optimization method that considers system configuration optimization and operation strategy, thereby achieving the improvement of system energy efficiency and economy.In the above-mentioned literature review, in terms of model establishment, some use nonlinear models for predictive control, and some use two-stage hybrid linear programming. Although many scholars have proposed a variety of models for optimal scheduling of the CCHP system. It can have a certain optimization effect in a specific environment, but the considerations are not very comprehensive. For example, in terms of wind power consumption, wind power cannot be absorbed well, and the uncertainty of wind power has not been well considered. Based on this, this paper establishes a refined energy storage model, considers the consumption of renewable energy, and proposes an optimized control strategy to solve these problems.The contents of this paper are structured as follows: In Section2, the CCHP system is introduced, and then the energy storage device model is built, including the battery energy storage model and the heat storage model. In Section 3, the demand response of electric load and heat load is studied, and the user’s electric load demand response is constructed based on the principle of price-based demand response. At the same time, this section mainly studies the user’s heat load demand response by considering the heating load. In Section 4, a collaborative optimization model is established. First, the uncertain new problem of wind power is transformed into a certainty problem. Then the objective function is established with the lowest operating cost of the integrated energy system. Finally, it is solved iteratively through the crow search algorithm. In Section 5, case analysis is used to verify the feasibility of the algorithm. In the last section, it is a summary of this article. The overall research idea is shown in Figure 1.Figure 1
Flow chart of research method.


## 2. System Description

### 2.1. CCHP System

The seasonal climate in southern China is changeable, with abundant renewable energy such as geothermal energy and offshore wind energy. There are geographical locations and technical conditions for the development of CCHP and ground source heat pumps. The IES structure established in the article is shown in Figure2. It is composed of photovoltaic (PV), wind turbine (WT), energy storage (ES) equipment, CCHP unit, and ground source heat pump (GSHP). The heat load demand is provided by ground source heat pumps and CCHP units. Energy storage equipment includes electric energy storage (EES) and heat storage tank (HST). The electric load demand is jointly supplied by the upper-level grid, wind turbines, and CCHP units. The CCHP unit operates in the system with “heat-fixed electricity” mode.Figure 2
Structure of the IES.The CCHP system can provide multiple types of energy at the same time. It uses the high-grade energy of natural gas combustion to produce electricity and uses a bromine cooler to recover waste heat for supplying cold and heat loads. A typical CCHP system consists of a gas turbine (GT) and a bromine cooler (BC), and its mathematical model is as follows:(1)PGTt=QGTt1−βGT−βlβGT,PBC,ht=QGTtδBCλBC,h,where PGTt, QGTt, PBC,ht represent generating power and exhaust gas residual heat of GT and the heating power of BC at time t, respectively; βGT, βl, δBC, λBC,h are power generation efficiency of GT and heat loss coefficient of GT, flue gas recovery of BC, and heating coefficient of BC, respectively.GSHP converts shallow energy into high quality energy by inputting part of high-quality electricity. Compared with the previous electric boiler and other equipment, it can output four to seven times the amount of cold or heat of the electricity consumed. By the end of 2019, China’s total installed capacity of ground source heat pumps had reached 20,000 MW, and the heating and cooling floor space of hydrothermal ground source heat pumps exceeded 150 million square meters. The mathematical model of its heating capacity can be expressed as follows:(2)Php,ht=ηop,hPhpt,where Php,ht, Phpt represent the heating power and power consumption of GSHP at time t, respectively. ηop,h is the heating efficiency of GSHP.


### 2.2. Energy Storage Device Modeling

As an important part of IES, energy storage devices can store excess energy and provide energy when users’ energy needs cannot be met. And it can effectively suppress the uncertainty of new energy output and has the effect of cutting peaks and filling valleys for system operation. The traditional energy charging and discharging mathematical model can be written as follows:(3)EESt=1−μEESt−1+UPchtηchΔt−VPdistηdisΔt,(4)U+V≤1,U∈0,1,V∈0,1,where EESt is the energy storage capacity of energy storage unit at time t; Pcht, Pdist, respectively, represent the charging and discharging rates of the energy storage unit at time t; μ, ηch, ηdis are the energy loss coefficient of energy storage unit, the conversion efficiency of energy charging, discharging, respectively. U, V are 0 and 1 variables. The energy storage unit charges when U = 1 and discharges when V = 1.


#### 2.2.1. The Model of Battery Energy Storage

The traditional energy storage model ignores the influence of temperature and other factors on the battery. The actual operating capacityEEESn is related to the ambient temperature. The mathematical expression is shown as follows:(5)EEESn=EEESN1+δaTat−TN,where EEESN is the rated capacity of EES in standard state, and the temperature TN in standard state is 25°C; Tat is the ambient temperature at t, δa is the capacity temperature coefficient, whose value is 0.6%In addition, the charging and discharging rate of EES is related to the service life. Too high charging and discharging rate will accelerate EES aging. Therefore, its recharging and discharging power cannot be too large, which can be expressed as follows:(6)PEES,ch≤achEEESntΔt,PEES,dis≤adisEEESntΔt,where PEES,ch, PEES,dis, respectively, represents the charge and discharge power of EES at time t. ach and adis are the charge and discharge coefficients of the battery. In order to ensure the safety of charge and discharge and the life of the battery, they are generally between 0.1 and 0.4.


#### 2.2.2. The Model of Thermal Storage Tank

The fine storage model of heat storage tank in this paper includes atmospheric heat storage tank, circulating pump and heat exchanger, and other additional equipment. Taking its analysis of heating mode in winter as an example, its principle is shown in Figure3. The HST exchanges heat with the user side only through the heat exchanger.Figure 3
Connection diagram of heat storage tank and heat supply network.When HST is storing and releasing heat energy, due to the different internal temperatures, water will be layered, with high-temperature hot water in the upper layer, low-temperature hot water in the lower layer, and oblique-temperature layer in the middle. When heating, hot water flows out from the upper pipe, passes through the heat exchanger, and then enters the tank from the lower pipe, while the inclined temperature layer moves up. On the contrary, low-temperature hot water flows out of the lower pipe during heat storage and flows in from the upper pipe after being heated by the heat exchanger, while the thermocline moved down. The HST mathematical model is shown below:(7)ENHST=ρcVTh−Tl,where ENHST is the rated heat storage when heating; ρ and c are, respectively, the density and specific heat capacity of water; V is the volume of HST water storage, which is a fixed value in this paper. Th and Tl are, respectively, the highest and lowest temperature for HST to maintain the water temperature in the tank.In order to satisfy the optimal operation of HST model in long time scale, the influence of liquid and temperature in the tank is ignored. If the water temperature of the hot water layer nozzle and the low temperature layer nozzle is maintained atTh and Tl, and the heat charge and release power of HST are determined by controlling the water flow rate It, then the mathematical model of HST can be expressed as follows:(8)EtHST=ρcVTh−Tl,Tl≤TtHST≤Th,PtHST=ρcItHSTTh−Tl,ΔEtHST=ρcItHSTTh−Tl−ΦtHST,loss,where EtHST, TtHST, PtHST are the energy storage, water temperature in the tank, and heat charge and release power of HST at time t, respectively, and ΦtHST,loss represents the heat loss of HST at time t.The heat loss of HST is mainly composed of two parts: one is the heat loss caused by the heat transfer of low temperature hot water and high temperature hot water in the thermocline layer and the heat transfer in the inner wall of the tank. It mainly shows its influence on HST energy storage efficiency. The charging and releasing efficiency of HST can be expressed byEchHST and ηdisHST. The second is the heat loss from the outer surface of the HST tank to the surrounding air. For the second heat loss Φthst,loss, the mathematical expression is as follows:(9)ΦtHST,loss=kSTtHST−Tta,S=4πaENHSTρcTh−Tl,where Tta represents the ambient temperature of HST at time t; K is the convection heat coefficient that comprehensively considers the influence of outer wall, insulation layer, and ambient temperature on the heat of HST working medium; S is the surface area of HST; and a is the height-diameter ratio of HST.According to the heat loss characteristics of the HST above, the energy storage state of the HST can be represented byEtHST/ENHST, or TtHST−T1/Th−Tl. Then equation (3) can be rewritten to obtain the fine energy storage model of HST:(10)EtHST=Et−1HST+ηchHSTPtHST,ch−PtHST,disηdisHSTΔt−kSEt−1HSTTh−TlENHST+Tl−TtaΔt.


### 2.3. Comprehensive Analysis of Demand Response of Electric Heating

#### 2.3.1. Demand Response of Electrical Load

In this paper, the PBDR (price-based demand response) modeling is based on the price-based demand response principle. PBDR is adopted to guide the demand side electricity consumption, so as to reduce the peak and valley load difference. The demand elasticity of user electrical load is expressed by the following formula:(11)Nt1,t2=ΔPt1P0t1D0t2ΔD0t2Nt1,t2≥0,t1≠t2,Nt1,t2≤0,t1=t2,where P0t1 represents the electric load at time t1 before taking the demand response, D0t2 represents electricity price at time t2 before taking the demand response, ΔPt1 represents the electric load at time t1 after taking the demand response, ΔDt2 represents electricity price at time t2 after taking the demand response, when t1 = t2, Nt1,t2 represent the self-elastic coefficient; when t1 ≠ t2, Nt1,t2 represent cross elasticity coefficient.After taking PBDR, the electrical load change expression can be shown as follows:(12)ΔP1P01ΔP2P02⋮ΔP24P024=N11…N124N21…N224⋮…⋮N241…N2424ΔD1D01ΔD2D02⋮ΔD24D024.After taking PBDR, the electric load at timet is(13)Pload,et=P0t1+ΔPtP0t.


#### 2.3.2. Heat Load Demand Response

Residential heat load includes heating load and domestic hot water. Considering that the heat generated by hot water accounts for a small proportion of the total heat load, and it is a rigid load with no elasticity. In order to reflect the impact of thermal load elasticity on the absorption of renewable energy in the system, only the residential heating load with elasticity is considered in this paper.The relationship between the indoor temperature variation of the heating building, the thermal power required by the building, and the current ambient temperature is as follows:(14)Pload,ht=M1RTint+1−e−Δt/τTint1−e−Δt/τ−Toutt,τ=RC,where Pload,ht is the heating power of the building at time t; Tint, Toutt are, respectively, the indoor and outdoor temperatures of the building at time t; R is the equivalent thermal resistance of the building, in unit of °C/kW; M1 is the amount of heating for residents; and C is the indoor air heat capacity, in unit of (kW·h)/°C.For residents, their perception of temperature comfort has a certain degree of ambiguity. When the heating temperature changes within a certain range, the comfort experience of residents will not be affected. The indoor temperature meets the following constraints:(15)Tin,min≤Tint≤Tin,max,where Tin,max, Tin,min, respectively, represent the maximum and minimum temperature acceptable to indoor residents, in unit of °C.


### 2.4. Wind Power Uncertainty

Wind power output has volatility and randomness, which usually results in a certain error between the actual output of wind power and the predicted output [11]. At present, in the research on the uncertainty of wind power, most literatures use Gaussian distribution function to fit wind power output, and the fitted value of normal distribution differs greatly from the actual value in certain specific periods. Therefore, the General distribution is used to describe the uncertainty of wind power output, and its function form is as follows:(16)fx=αβe−αx−γ1+e−αx−γβ+1,Fx=1+e−αx−γβ,F−1y=γ−1αlny−1/β−1,where f(x), F(x), F−1(y) are, respectively, the probability density function of the general distribution, the cumulative distribution function, and its inverse function; α, β, and γ are all shape parameters, where α > 0, β > 0, −∞ < γ < +∞; y represents the cumulative probability.In this paper, typical scenarios of wind power output forecasting errors are considered. By reducing it, the uncertainty of wind power output is transformed into a deterministic problem for solving. The prediction error ΔPwt,t of wind power output follows the general distribution with shape parameters α, β, and γ. The Latin hypercube sampling method is used to simulate wind power fluctuation scene sets, generating M initial scene sets. Then sampling synchronous back substitution method is used to reduce the generated scene. Every time a scene is subtracted, the remaining scene and its corresponding probability are updated, so that the scene after the reduction is closest to the initial scene probability. After repeated back-generation, the final prediction error scenario set of N wind power generating units are obtained, which is added to the predicted value to obtain the actual wind power output in the corresponding period.


### 2.5. Collaborative Optimization Modeling

At present, the synergistic scheduling potential of integrated energy use has not been well explored. There are still a lot of problems, including serious wind and light abandoning phenomenon, not well considered TOU price, and so on. There is no better model for optimal scheduling. In order to make better use of energy coordination, this paper considers many factors to model the optimization of integrated energy. The goal is to minimize the total operating cost of the system by adjusting the cost of all parts, including gas purchase cost, conventional unit power generation cost, and new energy unit power generation cost, within a certain cycle range, considering the constraints of load demand and so on. Detailed description is as follows.


#### 2.5.1. Objective Function

The article aims to optimize the total cost of IES operation, and the establishment of the fine energy storage model IES source storage and load coordination optimization scheduling objective function is as follows:(17)minF=∑t=1TCfut+Cwt+Cet+Cgridt+Cwtt,where F is the total economic cost of system operation; T is the scheduling period; and Cfut, Cwt, Cet, Cgridt, and Cwttd correspond to the gas purchase costs during t, unit operating costs, pollutant treatment costs, electricity interaction costs with the superior grid, and wind curtailment costs.Gas purchase costs can be expressed as follows:(18)Cfut=PGTtβGT⋅ΔtLCH4DCH4,where PGT, βGT, DCH4, and LCH4 correspond to the power of gas turbine at t time (kW), gas utilization efficiency, natural gas price (yuan/m3), and low calorific value (kW·h/m3). Δt is the intervals between the two periods (for example, if the period is 24 hours and the interval is adjusted once every hour, Δt is 1 h).Unit operation costs can be expressed as follows:(19)Cwt=∑i=1KPitWi+∑j=1N1PjtWj+∑n=1N2PES,ntWnΔt,where Pit, Pjt, PES,nt correspond to the power (kW) of controllable generator set i, new energy generator set j, and energy storage equipment n at time t. Wi, Wj, and Wn are the operating cost per unit time of controllable generator set i, new energy generator unit j, and energy storage equipment n (yuan/h). Where Pit includes the power generation power of wind turbines PWTt and the power generation power of photovoltaic generators PPVt. Among them, PES,nt includes electric energy storage PEES,nt and thermal energy storage PHST,nt. Where K, N1, and N2 correspond to the types of controllable units, new energy units, and energy storage equipment.Pollutant treatment costs can be expressed as follows:(20)Cet=∑l=1LθlElPitΔt,where Pit, El, θl, and L correspond to the power (kW) of the generator set at time t, the emission coefficient (/kW·h), the treatment cost (yuan) of pollutant l, and the type of pollutant.Electrical energy interaction costs with the upper grid can be expressed as follows:(21)Cgridt=PgbtmaxPgridt,0Δt+Pgstmax−Pgridt,0Δt,where Pgbt and Pgst correspond to the electricity purchase and sale price t (yuan/kW·h) of the grid at time t and Pgridt corresponds to the power of the grid connection line at time t, greater than 0 meaning electricity purchase, and less than 0 meaning electricity sale.The cost of wind power discards can be expressed as follows:(22)Cwtt=PwttDwtΔt,where Pwtt and Dwt correspond to the discarded power (kW) from wind power at time t and the cost (yuan/kW·h) of discards per unit of wind power generated.


#### 2.5.2. Restrictions

Operation needs to be within reasonable limits. In order to adjust IES safely and efficiently, the system and unit must meet the following constraints.(1)
Power balance constraint(23)Pload,et+Phpt=PGTt+PWTt+Pgridt+PPVt+PEEStPload,ht=PBC,ht+Php,htPHSTt,wherePload,et, Pload,ht correspond to the system electrical load power and thermal load power. The other variables have been described above.(2)
Gas purchase power constraint(24)0≤QGT,buy≤QGT,buymax,QGT,buy=PGTtΔtβGTLCH4,wherePGTt, βGT, QGT,buy, and QGT,buymax correspond to the gas turbine power generation, the gas turbine power generation efficiency, the amount of gas purchased, and the maximum amount of gas purchased.(3)
Operation constraints of controllable generator sets(25)Pi,min≤Pit≤Pi,max,−λi,downΔt≤Pit−Pit−1≤λi,upΔt,wherePi,max, Pi,min correspond to the upper and lower limits of the output of controllable unit i. λi,up, λi,down correspond to the upper and lower limits of the ramp rate of controllable unit i.(4)
Wind power constraint(26)0≤PWT≤PWTmax,wherePWT is the power generation output of the fan, and PWTmax the maximum predicted output of wind power generation.(5)
Photovoltaic generation constraint(27)0≤PPV≤PPVmax,wherePPV PPV is the photovoltaic power generation output, and PPVmax is the maximum predicted output of photovoltaic power generation.(6)
Interacting power constraints with the upper grid(28)Pgrid,min≤Pgridt≤Pgrid,max,wherePgrid,min is the minimum values of the interactive power of the power system and Pgrid,max is the maximum values of the interactive power of the system.(7)
Energy storage equipment constraints(29)PES,min≤PESt≤PES,max,Emin≤EESt≤Emax,EES0=EES24,where PES,max, PES,min are the maximum and minimum values of the charging and discharging power of the energy storage device within the constraints.(1) Crow Search Algorithm. Crow Search Algorithm (CSA) is a swarm intelligence search algorithm proposed to simulate crow’s action of hiding food, which has the advantages of fewer iterations and better searching ability. Assume that the number of crows is I (the number of decision variables in this paper), yi and n are the position of crow i after the nth search (the position in this paper refers to the value of decision variables, including PGT, Pi, PWT, PPV, Pgrid, and PES). At this time, the best hiding location it finds is called the crow’s memory location, denoted as Mi,n (the optimal value of the decision variable at this time), where i=1,2,…,I, n=1,2,…,nmax, nmax indicates the maximum number of searches. The main steps are as follows: Step 1: Initializeyi,1, Mi,1=yi,1. Give each decision variable an initial value.Step 2: Search for a new location. Each bird in the colony randomly selects another crow as the target for the search location. Assume that in iterationn, crow i follows Crow z, then there are two situations:(1)
Crowz does not find himself being followed by crow iCrowi will be close to the memory position MZ,n of crow z, and the position of crow i will be updated through equation (30).(30)yi.n+1=yi,t+Ci×fl×Mz,n−yi,n,whereyi,n+1, fl, Ci, Mz,n correspond to the position of crow i after (n + 1)th search, the step length, the uniformly distributed random number in the interval [0.1], and the best position to find at this time.(2)
Crowz finds himself being followed by crow iCrowi cannot get close to Crow z’s memory location and can only find a random new location to update its position.According to the description of the above two cases, the update formula of crowi’s position is summarized as follows:(31)yi,n+1=yi,n+1=yi,n+Ci×fl×Mz,n−yi,n,Cz≥APi,n,generate a random position,otherwise,whereAPi,t is the probability that crow can find being followed by crow i when nth search.Step 3: Update the memory location. The crows in the group update their positions each time. Determine whether to replace the memory location through the evaluation functionF (in this article, the objective function is used to represent the evaluation function of the system operating cost). The specific update formula is as follows:(32)Mi,n+1=yi,n+1,Fyi,n+1≤FMi,n,Mi,n,otherwise.Step 4: Repeat Step 2 to Step 3 to search fornmax times or the evaluation function does not meet the conditions to proceed to the next step.Step 5: Output optimization results. The optimal memory location in the crow group is the final optimization result of the algorithm.The above are the solving steps of crow algorithm, combined with the objective function of collaborative optimization. The model solving process is shown in Figure4.Figure 4
CSA flow chart for solving the optimization problem of the IES.


## 2.1. CCHP System

The seasonal climate in southern China is changeable, with abundant renewable energy such as geothermal energy and offshore wind energy. There are geographical locations and technical conditions for the development of CCHP and ground source heat pumps. The IES structure established in the article is shown in Figure2. It is composed of photovoltaic (PV), wind turbine (WT), energy storage (ES) equipment, CCHP unit, and ground source heat pump (GSHP). The heat load demand is provided by ground source heat pumps and CCHP units. Energy storage equipment includes electric energy storage (EES) and heat storage tank (HST). The electric load demand is jointly supplied by the upper-level grid, wind turbines, and CCHP units. The CCHP unit operates in the system with “heat-fixed electricity” mode.Figure 2
Structure of the IES.The CCHP system can provide multiple types of energy at the same time. It uses the high-grade energy of natural gas combustion to produce electricity and uses a bromine cooler to recover waste heat for supplying cold and heat loads. A typical CCHP system consists of a gas turbine (GT) and a bromine cooler (BC), and its mathematical model is as follows:(1)PGTt=QGTt1−βGT−βlβGT,PBC,ht=QGTtδBCλBC,h,where PGTt, QGTt, PBC,ht represent generating power and exhaust gas residual heat of GT and the heating power of BC at time t, respectively; βGT, βl, δBC, λBC,h are power generation efficiency of GT and heat loss coefficient of GT, flue gas recovery of BC, and heating coefficient of BC, respectively.GSHP converts shallow energy into high quality energy by inputting part of high-quality electricity. Compared with the previous electric boiler and other equipment, it can output four to seven times the amount of cold or heat of the electricity consumed. By the end of 2019, China’s total installed capacity of ground source heat pumps had reached 20,000 MW, and the heating and cooling floor space of hydrothermal ground source heat pumps exceeded 150 million square meters. The mathematical model of its heating capacity can be expressed as follows:(2)Php,ht=ηop,hPhpt,where Php,ht, Phpt represent the heating power and power consumption of GSHP at time t, respectively. ηop,h is the heating efficiency of GSHP.


## 2.2. Energy Storage Device Modeling

As an important part of IES, energy storage devices can store excess energy and provide energy when users’ energy needs cannot be met. And it can effectively suppress the uncertainty of new energy output and has the effect of cutting peaks and filling valleys for system operation. The traditional energy charging and discharging mathematical model can be written as follows:(3)EESt=1−μEESt−1+UPchtηchΔt−VPdistηdisΔt,(4)U+V≤1,U∈0,1,V∈0,1,where EESt is the energy storage capacity of energy storage unit at time t; Pcht, Pdist, respectively, represent the charging and discharging rates of the energy storage unit at time t; μ, ηch, ηdis are the energy loss coefficient of energy storage unit, the conversion efficiency of energy charging, discharging, respectively. U, V are 0 and 1 variables. The energy storage unit charges when U = 1 and discharges when V = 1.


### 2.2.1. The Model of Battery Energy Storage

The traditional energy storage model ignores the influence of temperature and other factors on the battery. The actual operating capacityEEESn is related to the ambient temperature. The mathematical expression is shown as follows:(5)EEESn=EEESN1+δaTat−TN,where EEESN is the rated capacity of EES in standard state, and the temperature TN in standard state is 25°C; Tat is the ambient temperature at t, δa is the capacity temperature coefficient, whose value is 0.6%In addition, the charging and discharging rate of EES is related to the service life. Too high charging and discharging rate will accelerate EES aging. Therefore, its recharging and discharging power cannot be too large, which can be expressed as follows:(6)PEES,ch≤achEEESntΔt,PEES,dis≤adisEEESntΔt,where PEES,ch, PEES,dis, respectively, represents the charge and discharge power of EES at time t. ach and adis are the charge and discharge coefficients of the battery. In order to ensure the safety of charge and discharge and the life of the battery, they are generally between 0.1 and 0.4.


### 2.2.2. The Model of Thermal Storage Tank

The fine storage model of heat storage tank in this paper includes atmospheric heat storage tank, circulating pump and heat exchanger, and other additional equipment. Taking its analysis of heating mode in winter as an example, its principle is shown in Figure3. The HST exchanges heat with the user side only through the heat exchanger.Figure 3
Connection diagram of heat storage tank and heat supply network.When HST is storing and releasing heat energy, due to the different internal temperatures, water will be layered, with high-temperature hot water in the upper layer, low-temperature hot water in the lower layer, and oblique-temperature layer in the middle. When heating, hot water flows out from the upper pipe, passes through the heat exchanger, and then enters the tank from the lower pipe, while the inclined temperature layer moves up. On the contrary, low-temperature hot water flows out of the lower pipe during heat storage and flows in from the upper pipe after being heated by the heat exchanger, while the thermocline moved down. The HST mathematical model is shown below:(7)ENHST=ρcVTh−Tl,where ENHST is the rated heat storage when heating; ρ and c are, respectively, the density and specific heat capacity of water; V is the volume of HST water storage, which is a fixed value in this paper. Th and Tl are, respectively, the highest and lowest temperature for HST to maintain the water temperature in the tank.In order to satisfy the optimal operation of HST model in long time scale, the influence of liquid and temperature in the tank is ignored. If the water temperature of the hot water layer nozzle and the low temperature layer nozzle is maintained atTh and Tl, and the heat charge and release power of HST are determined by controlling the water flow rate It, then the mathematical model of HST can be expressed as follows:(8)EtHST=ρcVTh−Tl,Tl≤TtHST≤Th,PtHST=ρcItHSTTh−Tl,ΔEtHST=ρcItHSTTh−Tl−ΦtHST,loss,where EtHST, TtHST, PtHST are the energy storage, water temperature in the tank, and heat charge and release power of HST at time t, respectively, and ΦtHST,loss represents the heat loss of HST at time t.The heat loss of HST is mainly composed of two parts: one is the heat loss caused by the heat transfer of low temperature hot water and high temperature hot water in the thermocline layer and the heat transfer in the inner wall of the tank. It mainly shows its influence on HST energy storage efficiency. The charging and releasing efficiency of HST can be expressed byEchHST and ηdisHST. The second is the heat loss from the outer surface of the HST tank to the surrounding air. For the second heat loss Φthst,loss, the mathematical expression is as follows:(9)ΦtHST,loss=kSTtHST−Tta,S=4πaENHSTρcTh−Tl,where Tta represents the ambient temperature of HST at time t; K is the convection heat coefficient that comprehensively considers the influence of outer wall, insulation layer, and ambient temperature on the heat of HST working medium; S is the surface area of HST; and a is the height-diameter ratio of HST.According to the heat loss characteristics of the HST above, the energy storage state of the HST can be represented byEtHST/ENHST, or TtHST−T1/Th−Tl. Then equation (3) can be rewritten to obtain the fine energy storage model of HST:(10)EtHST=Et−1HST+ηchHSTPtHST,ch−PtHST,disηdisHSTΔt−kSEt−1HSTTh−TlENHST+Tl−TtaΔt.


## 2.2.1. The Model of Battery Energy Storage

The traditional energy storage model ignores the influence of temperature and other factors on the battery. The actual operating capacityEEESn is related to the ambient temperature. The mathematical expression is shown as follows:(5)EEESn=EEESN1+δaTat−TN,where EEESN is the rated capacity of EES in standard state, and the temperature TN in standard state is 25°C; Tat is the ambient temperature at t, δa is the capacity temperature coefficient, whose value is 0.6%In addition, the charging and discharging rate of EES is related to the service life. Too high charging and discharging rate will accelerate EES aging. Therefore, its recharging and discharging power cannot be too large, which can be expressed as follows:(6)PEES,ch≤achEEESntΔt,PEES,dis≤adisEEESntΔt,where PEES,ch, PEES,dis, respectively, represents the charge and discharge power of EES at time t. ach and adis are the charge and discharge coefficients of the battery. In order to ensure the safety of charge and discharge and the life of the battery, they are generally between 0.1 and 0.4.


## 2.2.2. The Model of Thermal Storage Tank

The fine storage model of heat storage tank in this paper includes atmospheric heat storage tank, circulating pump and heat exchanger, and other additional equipment. Taking its analysis of heating mode in winter as an example, its principle is shown in Figure3. The HST exchanges heat with the user side only through the heat exchanger.Figure 3
Connection diagram of heat storage tank and heat supply network.When HST is storing and releasing heat energy, due to the different internal temperatures, water will be layered, with high-temperature hot water in the upper layer, low-temperature hot water in the lower layer, and oblique-temperature layer in the middle. When heating, hot water flows out from the upper pipe, passes through the heat exchanger, and then enters the tank from the lower pipe, while the inclined temperature layer moves up. On the contrary, low-temperature hot water flows out of the lower pipe during heat storage and flows in from the upper pipe after being heated by the heat exchanger, while the thermocline moved down. The HST mathematical model is shown below:(7)ENHST=ρcVTh−Tl,where ENHST is the rated heat storage when heating; ρ and c are, respectively, the density and specific heat capacity of water; V is the volume of HST water storage, which is a fixed value in this paper. Th and Tl are, respectively, the highest and lowest temperature for HST to maintain the water temperature in the tank.In order to satisfy the optimal operation of HST model in long time scale, the influence of liquid and temperature in the tank is ignored. If the water temperature of the hot water layer nozzle and the low temperature layer nozzle is maintained atTh and Tl, and the heat charge and release power of HST are determined by controlling the water flow rate It, then the mathematical model of HST can be expressed as follows:(8)EtHST=ρcVTh−Tl,Tl≤TtHST≤Th,PtHST=ρcItHSTTh−Tl,ΔEtHST=ρcItHSTTh−Tl−ΦtHST,loss,where EtHST, TtHST, PtHST are the energy storage, water temperature in the tank, and heat charge and release power of HST at time t, respectively, and ΦtHST,loss represents the heat loss of HST at time t.The heat loss of HST is mainly composed of two parts: one is the heat loss caused by the heat transfer of low temperature hot water and high temperature hot water in the thermocline layer and the heat transfer in the inner wall of the tank. It mainly shows its influence on HST energy storage efficiency. The charging and releasing efficiency of HST can be expressed byEchHST and ηdisHST. The second is the heat loss from the outer surface of the HST tank to the surrounding air. For the second heat loss Φthst,loss, the mathematical expression is as follows:(9)ΦtHST,loss=kSTtHST−Tta,S=4πaENHSTρcTh−Tl,where Tta represents the ambient temperature of HST at time t; K is the convection heat coefficient that comprehensively considers the influence of outer wall, insulation layer, and ambient temperature on the heat of HST working medium; S is the surface area of HST; and a is the height-diameter ratio of HST.According to the heat loss characteristics of the HST above, the energy storage state of the HST can be represented byEtHST/ENHST, or TtHST−T1/Th−Tl. Then equation (3) can be rewritten to obtain the fine energy storage model of HST:(10)EtHST=Et−1HST+ηchHSTPtHST,ch−PtHST,disηdisHSTΔt−kSEt−1HSTTh−TlENHST+Tl−TtaΔt.


## 2.3. Comprehensive Analysis of Demand Response of Electric Heating

### 2.3.1. Demand Response of Electrical Load

In this paper, the PBDR (price-based demand response) modeling is based on the price-based demand response principle. PBDR is adopted to guide the demand side electricity consumption, so as to reduce the peak and valley load difference. The demand elasticity of user electrical load is expressed by the following formula:(11)Nt1,t2=ΔPt1P0t1D0t2ΔD0t2Nt1,t2≥0,t1≠t2,Nt1,t2≤0,t1=t2,where P0t1 represents the electric load at time t1 before taking the demand response, D0t2 represents electricity price at time t2 before taking the demand response, ΔPt1 represents the electric load at time t1 after taking the demand response, ΔDt2 represents electricity price at time t2 after taking the demand response, when t1 = t2, Nt1,t2 represent the self-elastic coefficient; when t1 ≠ t2, Nt1,t2 represent cross elasticity coefficient.After taking PBDR, the electrical load change expression can be shown as follows:(12)ΔP1P01ΔP2P02⋮ΔP24P024=N11…N124N21…N224⋮…⋮N241…N2424ΔD1D01ΔD2D02⋮ΔD24D024.After taking PBDR, the electric load at timet is(13)Pload,et=P0t1+ΔPtP0t.


### 2.3.2. Heat Load Demand Response

Residential heat load includes heating load and domestic hot water. Considering that the heat generated by hot water accounts for a small proportion of the total heat load, and it is a rigid load with no elasticity. In order to reflect the impact of thermal load elasticity on the absorption of renewable energy in the system, only the residential heating load with elasticity is considered in this paper.The relationship between the indoor temperature variation of the heating building, the thermal power required by the building, and the current ambient temperature is as follows:(14)Pload,ht=M1RTint+1−e−Δt/τTint1−e−Δt/τ−Toutt,τ=RC,where Pload,ht is the heating power of the building at time t; Tint, Toutt are, respectively, the indoor and outdoor temperatures of the building at time t; R is the equivalent thermal resistance of the building, in unit of °C/kW; M1 is the amount of heating for residents; and C is the indoor air heat capacity, in unit of (kW·h)/°C.For residents, their perception of temperature comfort has a certain degree of ambiguity. When the heating temperature changes within a certain range, the comfort experience of residents will not be affected. The indoor temperature meets the following constraints:(15)Tin,min≤Tint≤Tin,max,where Tin,max, Tin,min, respectively, represent the maximum and minimum temperature acceptable to indoor residents, in unit of °C.


## 2.3.1. Demand Response of Electrical Load

In this paper, the PBDR (price-based demand response) modeling is based on the price-based demand response principle. PBDR is adopted to guide the demand side electricity consumption, so as to reduce the peak and valley load difference. The demand elasticity of user electrical load is expressed by the following formula:(11)Nt1,t2=ΔPt1P0t1D0t2ΔD0t2Nt1,t2≥0,t1≠t2,Nt1,t2≤0,t1=t2,where P0t1 represents the electric load at time t1 before taking the demand response, D0t2 represents electricity price at time t2 before taking the demand response, ΔPt1 represents the electric load at time t1 after taking the demand response, ΔDt2 represents electricity price at time t2 after taking the demand response, when t1 = t2, Nt1,t2 represent the self-elastic coefficient; when t1 ≠ t2, Nt1,t2 represent cross elasticity coefficient.After taking PBDR, the electrical load change expression can be shown as follows:(12)ΔP1P01ΔP2P02⋮ΔP24P024=N11…N124N21…N224⋮…⋮N241…N2424ΔD1D01ΔD2D02⋮ΔD24D024.After taking PBDR, the electric load at timet is(13)Pload,et=P0t1+ΔPtP0t.


## 2.3.2. Heat Load Demand Response

Residential heat load includes heating load and domestic hot water. Considering that the heat generated by hot water accounts for a small proportion of the total heat load, and it is a rigid load with no elasticity. In order to reflect the impact of thermal load elasticity on the absorption of renewable energy in the system, only the residential heating load with elasticity is considered in this paper.The relationship between the indoor temperature variation of the heating building, the thermal power required by the building, and the current ambient temperature is as follows:(14)Pload,ht=M1RTint+1−e−Δt/τTint1−e−Δt/τ−Toutt,τ=RC,where Pload,ht is the heating power of the building at time t; Tint, Toutt are, respectively, the indoor and outdoor temperatures of the building at time t; R is the equivalent thermal resistance of the building, in unit of °C/kW; M1 is the amount of heating for residents; and C is the indoor air heat capacity, in unit of (kW·h)/°C.For residents, their perception of temperature comfort has a certain degree of ambiguity. When the heating temperature changes within a certain range, the comfort experience of residents will not be affected. The indoor temperature meets the following constraints:(15)Tin,min≤Tint≤Tin,max,where Tin,max, Tin,min, respectively, represent the maximum and minimum temperature acceptable to indoor residents, in unit of °C.


## 2.4. Wind Power Uncertainty

Wind power output has volatility and randomness, which usually results in a certain error between the actual output of wind power and the predicted output [11]. At present, in the research on the uncertainty of wind power, most literatures use Gaussian distribution function to fit wind power output, and the fitted value of normal distribution differs greatly from the actual value in certain specific periods. Therefore, the General distribution is used to describe the uncertainty of wind power output, and its function form is as follows:(16)fx=αβe−αx−γ1+e−αx−γβ+1,Fx=1+e−αx−γβ,F−1y=γ−1αlny−1/β−1,where f(x), F(x), F−1(y) are, respectively, the probability density function of the general distribution, the cumulative distribution function, and its inverse function; α, β, and γ are all shape parameters, where α > 0, β > 0, −∞ < γ < +∞; y represents the cumulative probability.In this paper, typical scenarios of wind power output forecasting errors are considered. By reducing it, the uncertainty of wind power output is transformed into a deterministic problem for solving. The prediction error ΔPwt,t of wind power output follows the general distribution with shape parameters α, β, and γ. The Latin hypercube sampling method is used to simulate wind power fluctuation scene sets, generating M initial scene sets. Then sampling synchronous back substitution method is used to reduce the generated scene. Every time a scene is subtracted, the remaining scene and its corresponding probability are updated, so that the scene after the reduction is closest to the initial scene probability. After repeated back-generation, the final prediction error scenario set of N wind power generating units are obtained, which is added to the predicted value to obtain the actual wind power output in the corresponding period.


## 2.5. Collaborative Optimization Modeling

At present, the synergistic scheduling potential of integrated energy use has not been well explored. There are still a lot of problems, including serious wind and light abandoning phenomenon, not well considered TOU price, and so on. There is no better model for optimal scheduling. In order to make better use of energy coordination, this paper considers many factors to model the optimization of integrated energy. The goal is to minimize the total operating cost of the system by adjusting the cost of all parts, including gas purchase cost, conventional unit power generation cost, and new energy unit power generation cost, within a certain cycle range, considering the constraints of load demand and so on. Detailed description is as follows.


### 2.5.1. Objective Function

The article aims to optimize the total cost of IES operation, and the establishment of the fine energy storage model IES source storage and load coordination optimization scheduling objective function is as follows:(17)minF=∑t=1TCfut+Cwt+Cet+Cgridt+Cwtt,where F is the total economic cost of system operation; T is the scheduling period; and Cfut, Cwt, Cet, Cgridt, and Cwttd correspond to the gas purchase costs during t, unit operating costs, pollutant treatment costs, electricity interaction costs with the superior grid, and wind curtailment costs.Gas purchase costs can be expressed as follows:(18)Cfut=PGTtβGT⋅ΔtLCH4DCH4,where PGT, βGT, DCH4, and LCH4 correspond to the power of gas turbine at t time (kW), gas utilization efficiency, natural gas price (yuan/m3), and low calorific value (kW·h/m3). Δt is the intervals between the two periods (for example, if the period is 24 hours and the interval is adjusted once every hour, Δt is 1 h).Unit operation costs can be expressed as follows:(19)Cwt=∑i=1KPitWi+∑j=1N1PjtWj+∑n=1N2PES,ntWnΔt,where Pit, Pjt, PES,nt correspond to the power (kW) of controllable generator set i, new energy generator set j, and energy storage equipment n at time t. Wi, Wj, and Wn are the operating cost per unit time of controllable generator set i, new energy generator unit j, and energy storage equipment n (yuan/h). Where Pit includes the power generation power of wind turbines PWTt and the power generation power of photovoltaic generators PPVt. Among them, PES,nt includes electric energy storage PEES,nt and thermal energy storage PHST,nt. Where K, N1, and N2 correspond to the types of controllable units, new energy units, and energy storage equipment.Pollutant treatment costs can be expressed as follows:(20)Cet=∑l=1LθlElPitΔt,where Pit, El, θl, and L correspond to the power (kW) of the generator set at time t, the emission coefficient (/kW·h), the treatment cost (yuan) of pollutant l, and the type of pollutant.Electrical energy interaction costs with the upper grid can be expressed as follows:(21)Cgridt=PgbtmaxPgridt,0Δt+Pgstmax−Pgridt,0Δt,where Pgbt and Pgst correspond to the electricity purchase and sale price t (yuan/kW·h) of the grid at time t and Pgridt corresponds to the power of the grid connection line at time t, greater than 0 meaning electricity purchase, and less than 0 meaning electricity sale.The cost of wind power discards can be expressed as follows:(22)Cwtt=PwttDwtΔt,where Pwtt and Dwt correspond to the discarded power (kW) from wind power at time t and the cost (yuan/kW·h) of discards per unit of wind power generated.


### 2.5.2. Restrictions

Operation needs to be within reasonable limits. In order to adjust IES safely and efficiently, the system and unit must meet the following constraints.(1)
Power balance constraint(23)Pload,et+Phpt=PGTt+PWTt+Pgridt+PPVt+PEEStPload,ht=PBC,ht+Php,htPHSTt,wherePload,et, Pload,ht correspond to the system electrical load power and thermal load power. The other variables have been described above.(2)
Gas purchase power constraint(24)0≤QGT,buy≤QGT,buymax,QGT,buy=PGTtΔtβGTLCH4,wherePGTt, βGT, QGT,buy, and QGT,buymax correspond to the gas turbine power generation, the gas turbine power generation efficiency, the amount of gas purchased, and the maximum amount of gas purchased.(3)
Operation constraints of controllable generator sets(25)Pi,min≤Pit≤Pi,max,−λi,downΔt≤Pit−Pit−1≤λi,upΔt,wherePi,max, Pi,min correspond to the upper and lower limits of the output of controllable unit i. λi,up, λi,down correspond to the upper and lower limits of the ramp rate of controllable unit i.(4)
Wind power constraint(26)0≤PWT≤PWTmax,wherePWT is the power generation output of the fan, and PWTmax the maximum predicted output of wind power generation.(5)
Photovoltaic generation constraint(27)0≤PPV≤PPVmax,wherePPV PPV is the photovoltaic power generation output, and PPVmax is the maximum predicted output of photovoltaic power generation.(6)
Interacting power constraints with the upper grid(28)Pgrid,min≤Pgridt≤Pgrid,max,wherePgrid,min is the minimum values of the interactive power of the power system and Pgrid,max is the maximum values of the interactive power of the system.(7)
Energy storage equipment constraints(29)PES,min≤PESt≤PES,max,Emin≤EESt≤Emax,EES0=EES24,where PES,max, PES,min are the maximum and minimum values of the charging and discharging power of the energy storage device within the constraints.(1) Crow Search Algorithm. Crow Search Algorithm (CSA) is a swarm intelligence search algorithm proposed to simulate crow’s action of hiding food, which has the advantages of fewer iterations and better searching ability. Assume that the number of crows is I (the number of decision variables in this paper), yi and n are the position of crow i after the nth search (the position in this paper refers to the value of decision variables, including PGT, Pi, PWT, PPV, Pgrid, and PES). At this time, the best hiding location it finds is called the crow’s memory location, denoted as Mi,n (the optimal value of the decision variable at this time), where i=1,2,…,I, n=1,2,…,nmax, nmax indicates the maximum number of searches. The main steps are as follows: Step 1: Initializeyi,1, Mi,1=yi,1. Give each decision variable an initial value.Step 2: Search for a new location. Each bird in the colony randomly selects another crow as the target for the search location. Assume that in iterationn, crow i follows Crow z, then there are two situations:(1)
Crowz does not find himself being followed by crow iCrowi will be close to the memory position MZ,n of crow z, and the position of crow i will be updated through equation (30).(30)yi.n+1=yi,t+Ci×fl×Mz,n−yi,n,whereyi,n+1, fl, Ci, Mz,n correspond to the position of crow i after (n + 1)th search, the step length, the uniformly distributed random number in the interval [0.1], and the best position to find at this time.(2)
Crowz finds himself being followed by crow iCrowi cannot get close to Crow z’s memory location and can only find a random new location to update its position.According to the description of the above two cases, the update formula of crowi’s position is summarized as follows:(31)yi,n+1=yi,n+1=yi,n+Ci×fl×Mz,n−yi,n,Cz≥APi,n,generate a random position,otherwise,whereAPi,t is the probability that crow can find being followed by crow i when nth search.Step 3: Update the memory location. The crows in the group update their positions each time. Determine whether to replace the memory location through the evaluation functionF (in this article, the objective function is used to represent the evaluation function of the system operating cost). The specific update formula is as follows:(32)Mi,n+1=yi,n+1,Fyi,n+1≤FMi,n,Mi,n,otherwise.Step 4: Repeat Step 2 to Step 3 to search fornmax times or the evaluation function does not meet the conditions to proceed to the next step.Step 5: Output optimization results. The optimal memory location in the crow group is the final optimization result of the algorithm.The above are the solving steps of crow algorithm, combined with the objective function of collaborative optimization. The model solving process is shown in Figure4.Figure 4
CSA flow chart for solving the optimization problem of the IES.


## 2.5.1. Objective Function

The article aims to optimize the total cost of IES operation, and the establishment of the fine energy storage model IES source storage and load coordination optimization scheduling objective function is as follows:(17)minF=∑t=1TCfut+Cwt+Cet+Cgridt+Cwtt,where F is the total economic cost of system operation; T is the scheduling period; and Cfut, Cwt, Cet, Cgridt, and Cwttd correspond to the gas purchase costs during t, unit operating costs, pollutant treatment costs, electricity interaction costs with the superior grid, and wind curtailment costs.Gas purchase costs can be expressed as follows:(18)Cfut=PGTtβGT⋅ΔtLCH4DCH4,where PGT, βGT, DCH4, and LCH4 correspond to the power of gas turbine at t time (kW), gas utilization efficiency, natural gas price (yuan/m3), and low calorific value (kW·h/m3). Δt is the intervals between the two periods (for example, if the period is 24 hours and the interval is adjusted once every hour, Δt is 1 h).Unit operation costs can be expressed as follows:(19)Cwt=∑i=1KPitWi+∑j=1N1PjtWj+∑n=1N2PES,ntWnΔt,where Pit, Pjt, PES,nt correspond to the power (kW) of controllable generator set i, new energy generator set j, and energy storage equipment n at time t. Wi, Wj, and Wn are the operating cost per unit time of controllable generator set i, new energy generator unit j, and energy storage equipment n (yuan/h). Where Pit includes the power generation power of wind turbines PWTt and the power generation power of photovoltaic generators PPVt. Among them, PES,nt includes electric energy storage PEES,nt and thermal energy storage PHST,nt. Where K, N1, and N2 correspond to the types of controllable units, new energy units, and energy storage equipment.Pollutant treatment costs can be expressed as follows:(20)Cet=∑l=1LθlElPitΔt,where Pit, El, θl, and L correspond to the power (kW) of the generator set at time t, the emission coefficient (/kW·h), the treatment cost (yuan) of pollutant l, and the type of pollutant.Electrical energy interaction costs with the upper grid can be expressed as follows:(21)Cgridt=PgbtmaxPgridt,0Δt+Pgstmax−Pgridt,0Δt,where Pgbt and Pgst correspond to the electricity purchase and sale price t (yuan/kW·h) of the grid at time t and Pgridt corresponds to the power of the grid connection line at time t, greater than 0 meaning electricity purchase, and less than 0 meaning electricity sale.The cost of wind power discards can be expressed as follows:(22)Cwtt=PwttDwtΔt,where Pwtt and Dwt correspond to the discarded power (kW) from wind power at time t and the cost (yuan/kW·h) of discards per unit of wind power generated.


## 2.5.2. Restrictions

Operation needs to be within reasonable limits. In order to adjust IES safely and efficiently, the system and unit must meet the following constraints.(1)
Power balance constraint(23)Pload,et+Phpt=PGTt+PWTt+Pgridt+PPVt+PEEStPload,ht=PBC,ht+Php,htPHSTt,wherePload,et, Pload,ht correspond to the system electrical load power and thermal load power. The other variables have been described above.(2)
Gas purchase power constraint(24)0≤QGT,buy≤QGT,buymax,QGT,buy=PGTtΔtβGTLCH4,wherePGTt, βGT, QGT,buy, and QGT,buymax correspond to the gas turbine power generation, the gas turbine power generation efficiency, the amount of gas purchased, and the maximum amount of gas purchased.(3)
Operation constraints of controllable generator sets(25)Pi,min≤Pit≤Pi,max,−λi,downΔt≤Pit−Pit−1≤λi,upΔt,wherePi,max, Pi,min correspond to the upper and lower limits of the output of controllable unit i. λi,up, λi,down correspond to the upper and lower limits of the ramp rate of controllable unit i.(4)
Wind power constraint(26)0≤PWT≤PWTmax,wherePWT is the power generation output of the fan, and PWTmax the maximum predicted output of wind power generation.(5)
Photovoltaic generation constraint(27)0≤PPV≤PPVmax,wherePPV PPV is the photovoltaic power generation output, and PPVmax is the maximum predicted output of photovoltaic power generation.(6)
Interacting power constraints with the upper grid(28)Pgrid,min≤Pgridt≤Pgrid,max,wherePgrid,min is the minimum values of the interactive power of the power system and Pgrid,max is the maximum values of the interactive power of the system.(7)
Energy storage equipment constraints(29)PES,min≤PESt≤PES,max,Emin≤EESt≤Emax,EES0=EES24,where PES,max, PES,min are the maximum and minimum values of the charging and discharging power of the energy storage device within the constraints.(1) Crow Search Algorithm. Crow Search Algorithm (CSA) is a swarm intelligence search algorithm proposed to simulate crow’s action of hiding food, which has the advantages of fewer iterations and better searching ability. Assume that the number of crows is I (the number of decision variables in this paper), yi and n are the position of crow i after the nth search (the position in this paper refers to the value of decision variables, including PGT, Pi, PWT, PPV, Pgrid, and PES). At this time, the best hiding location it finds is called the crow’s memory location, denoted as Mi,n (the optimal value of the decision variable at this time), where i=1,2,…,I, n=1,2,…,nmax, nmax indicates the maximum number of searches. The main steps are as follows: Step 1: Initializeyi,1, Mi,1=yi,1. Give each decision variable an initial value.Step 2: Search for a new location. Each bird in the colony randomly selects another crow as the target for the search location. Assume that in iterationn, crow i follows Crow z, then there are two situations:(1)
Crowz does not find himself being followed by crow iCrowi will be close to the memory position MZ,n of crow z, and the position of crow i will be updated through equation (30).(30)yi.n+1=yi,t+Ci×fl×Mz,n−yi,n,whereyi,n+1, fl, Ci, Mz,n correspond to the position of crow i after (n + 1)th search, the step length, the uniformly distributed random number in the interval [0.1], and the best position to find at this time.(2)
Crowz finds himself being followed by crow iCrowi cannot get close to Crow z’s memory location and can only find a random new location to update its position.According to the description of the above two cases, the update formula of crowi’s position is summarized as follows:(31)yi,n+1=yi,n+1=yi,n+Ci×fl×Mz,n−yi,n,Cz≥APi,n,generate a random position,otherwise,whereAPi,t is the probability that crow can find being followed by crow i when nth search.Step 3: Update the memory location. The crows in the group update their positions each time. Determine whether to replace the memory location through the evaluation functionF (in this article, the objective function is used to represent the evaluation function of the system operating cost). The specific update formula is as follows:(32)Mi,n+1=yi,n+1,Fyi,n+1≤FMi,n,Mi,n,otherwise.Step 4: Repeat Step 2 to Step 3 to search fornmax times or the evaluation function does not meet the conditions to proceed to the next step.Step 5: Output optimization results. The optimal memory location in the crow group is the final optimization result of the algorithm.The above are the solving steps of crow algorithm, combined with the objective function of collaborative optimization. The model solving process is shown in Figure4.Figure 4
CSA flow chart for solving the optimization problem of the IES.


## 3. Example Analysis

### 3.1. Scene Introduction

The purpose of case verification in this section is to verify the effectiveness and economy of the model in wind power consumption and system optimization. A simulation example is carried out based on the data in Figures5 and 6. The corresponding forecast data of wind turbine, photovoltaic, and electric-heat load are shown in Figure 5. The peak-valley time-of-use electricity price of the system and the typical daily outdoor temperature in winter are shown in Figures 6 and 7.Figure 5
Curve of electric and heat load and forecast output of wind turbine, photovoltaic.Figure 6
Electricity price.Figure 7
Ambient temperature in winter.Natural gasDCH4 and LCH4 take 2.54 yuan/m3 and 9.7 kW·h/m3, respectively. Gas turbine heat dissipation loss coefficient βl = 0.15; bromine cooler flue gas recovery rate λBC=0.85; and heating power λBC,h=1.2. The thermal heating user M1 = 400, the equivalent thermal resistance of the building R = 16, the price demand response self-elasticity coefficient Nt1,t1 is −0.2, and the mutual elasticity coefficient Nt1,t2 is 0.03. Human body comfort temperature Tin,max is 24°C and Tin,min is 22°C. Water density ρ = 971 kg/m3, height-to-diameter ratio a=1, specific heat capacity c=1.167×10−3 kW·h/(kg·K), the commutation coefficient k = 2.668 W/(m2·K), and the heat capacity of the air C = 0.525.The data of each unit of the IES are shown in Table1. Energy storage device data are shown in Table 2. The discharge and treatment cost data of various pollutants are shown in Table 3.Table 1
Operation data of each unit of the system.
ParameterGTGSHPPWPVGridPower limit/kW50030400100100Power lower limit/kW500000Upper limit of climbing rate/(kW/min)64———Lower limit of climbing rate/(kW/min)53———Efficiency0.243———Operation and maintenance unit price (yuan/kW·h)0.0530.0260.0290.025—Table 2
Parameters of energy storage device.
ParameterEESHSTParameterEESHSTCapacity/(kW·h)200300Initial energy storage state0.20.2Charge and discharge rate0.90.88Maximum energy storage state0.90.9Self-consumption rate0.0010.01Minimum energy storage state0.20Operation and maintenance unit price (yuan/kW·h)0.0510.045Maximum charge and discharge power/kW5050Table 3
Pollutant discharge and treatment costs.
TypesSO2NOXCO2Gas turbine emission standards (g/kWh)0.0234.795170.16Emission standard for purchased electricity (g/kWh)6.42.32696Treatment costs of various pollutants (yuan/t)100019509.75


### 3.2. Optimal Operation Analysis

In order to verify the effectiveness of the IES source-storage-load coordination optimization model that takes into account the fine energy storage model in system operation scheduling and wind power consumption. The following four simulation scenarios are set up in the article for verification, of which scenario 4 is the detailed analysis model in the paper. Scenario 1, nonground source heat pump and energy storage, excluding IDR; scenario 2, including ground source heat pump and traditional energy storage model, excluding IDR; scenario 3, including ground source heat pump and fine energy storage model, excluding IDR; scenario 4, including ground source heat pump and fine energy storage model, taking IDR into account.Figure8 shows the electrical and thermal load demand curves of the system considering IDR. Compared with before the response, the peak-valley difference of electrical and thermal loads decreases to different degrees. Under the demand response mechanism, part of the peak load is transferred to the valley period, which makes the load curve smoother after coupling and makes the system have better effect of peak cutting and valley filling. It also improves the economy of IES operation.Figure 8
Electric and heat load curves before and after response.Taking the ground-source heat pump, the fine energy storage model, and the analysis of scenario 4 considering IDR as an example, the optimized output of each device in the system is shown in Figures9 and 10.Figure 9
Optimization result of power supply of the system.Figure 10
Results of system heating optimization.After considering IDR, part of the electricity and heat load during the peak period of the system is transferred to the valley period, which improves the wind power consumption rate and operating economy. At 1 : 00–7:00 and 23 : 00–24 : 00, the electricity price is in the lowest period. At this point, the power generation cost of the CCHP unit in the system is greater than the electricity price of the grid, so the power generation power of the CCHP unit “fixed by heat” can be reduced by reducing the indoor temperature to improve the wind power consumption capacity of the system. Ground source heat pump at this stage of the heating cost is relatively low, preferentially used for heating.At 7 : 00–8:00 and 11 : 00–18 : 00, the electricity price is moderate. At this time, the power generation cost of the system CCHP unit is still greater than the grid electricity price. When the electrical load demand is small, reduce the indoor temperature to reduce the system power generation cost; in the peak period of electric load demand, the room temperature is raised and the power of CCHP unit is increased, so as to reduce the cost of IES buying electricity from external grid in the whole dispatching cycle. At this time, the thermal power of the system is first provided by the CCHP unit. When there is a shortage of thermal power, the ground source heat pump is started for heating.At 8 : 00–11 : 00 and 18 : 00–23 : 00, the electricity price is at the highest value, and the system power generation cost is lower than the grid electricity price. From 18 : 00 to 23 : 00, the electricity/heat load demand is relatively high. By increasing the room temperature, the power generation of the CCHP unit “using heat” is maximized, and the system power purchase cost is reduced. In this period, the thermal power is completely provided by the CCHP unit, the excess thermal energy of the system is stored in the HST, and the ground source heat pump does not work.During the scheduling timeT, EES is driven by the peak-valley price mechanism to store energy during the peak period of electricity price and release energy during the peak period of electricity price and load, so as to reduce the operating cost of IES and cut peaks and fill valleys. The heat storage tank stores the excess heat energy in the system and provides energy when the heat load demand cannot be met; this reduces the operating cost of the system and further validates the economics of the model.


### 3.3. Analysis of Example Results

The operating cost and wind power consumption rate of IES are shown in Tables4 and 5. It can be seen from Table 4 that Scenarios 2 and 3 are compared with Scenario 1, in terms of abandonment cost, gas purchase cost, and total cost, which have significantly reduced, which proves the economy of the ground source heat pump and energy storage device for the system. Scenario 3, which considers the refined energy storage model, has slightly higher operating costs than Scenario 2. It is because the fine energy storage model considers the influence of environmental temperature and other constraints, which increases the operating cost. However, the model accurately describes the equipment output, making it more practical and schedulable, and its increased cost is within an acceptable range. In Scenario 4, compared with Scenario 1, the cost of gas purchase was reduced by 25.65%, the cost of pollutant emission was reduced by 20.82%, the cost of abandonment was reduced by 100%, and the total cost was reduced by 10.12%.Table 4
Operation cost of each scenario.
Cost category value (yuan/day)Scene 1Scene 2Scene 3Scene 4Purchase cost of gas57894347.604382.404362.10Operation and maintenance costs466.67447.37440.22434.42Pollutant treatment cost6755.0455.5653.05Power purchase cost454.05858.15872.62609.47Electricity sales revenue1005286.90295.51252.32Wind curtailment cost21.064.224.890Total cost5792.785425.485460.185206.72Table 5
Total amount and rate of wind power absorbed by the system.
ParameterScene 1Scene 2Scene 3Scene 4Total wind power consumption (kW)4592.555294.385255.065471Wind power consumption rate83.94%96.77%96.05%100%It can be seen from Table5 that after considering ground source heat pumps and energy storage equipment, the system’s wind power consumption rate has increased significantly, but electricity generated by the wind is still discarded at night; Scenario 4, after comprehensive consideration of IDR, ground source heat pumps, and energy storage equipment, the system. The wind power consumption rate reaches 100%, which solves the phenomenon of wind abandonment caused by the system’s “heat-fixed power” operation mode. Through the analysis, the results verify that the IES system with the fine energy storage model in the paper has the ability to completely absorb wind power and significant economic benefits. It provides corresponding solutions for the future wind power consumption.


## 3.1. Scene Introduction

The purpose of case verification in this section is to verify the effectiveness and economy of the model in wind power consumption and system optimization. A simulation example is carried out based on the data in Figures5 and 6. The corresponding forecast data of wind turbine, photovoltaic, and electric-heat load are shown in Figure 5. The peak-valley time-of-use electricity price of the system and the typical daily outdoor temperature in winter are shown in Figures 6 and 7.Figure 5
Curve of electric and heat load and forecast output of wind turbine, photovoltaic.Figure 6
Electricity price.Figure 7
Ambient temperature in winter.Natural gasDCH4 and LCH4 take 2.54 yuan/m3 and 9.7 kW·h/m3, respectively. Gas turbine heat dissipation loss coefficient βl = 0.15; bromine cooler flue gas recovery rate λBC=0.85; and heating power λBC,h=1.2. The thermal heating user M1 = 400, the equivalent thermal resistance of the building R = 16, the price demand response self-elasticity coefficient Nt1,t1 is −0.2, and the mutual elasticity coefficient Nt1,t2 is 0.03. Human body comfort temperature Tin,max is 24°C and Tin,min is 22°C. Water density ρ = 971 kg/m3, height-to-diameter ratio a=1, specific heat capacity c=1.167×10−3 kW·h/(kg·K), the commutation coefficient k = 2.668 W/(m2·K), and the heat capacity of the air C = 0.525.The data of each unit of the IES are shown in Table1. Energy storage device data are shown in Table 2. The discharge and treatment cost data of various pollutants are shown in Table 3.Table 1
Operation data of each unit of the system.
ParameterGTGSHPPWPVGridPower limit/kW50030400100100Power lower limit/kW500000Upper limit of climbing rate/(kW/min)64———Lower limit of climbing rate/(kW/min)53———Efficiency0.243———Operation and maintenance unit price (yuan/kW·h)0.0530.0260.0290.025—Table 2
Parameters of energy storage device.
ParameterEESHSTParameterEESHSTCapacity/(kW·h)200300Initial energy storage state0.20.2Charge and discharge rate0.90.88Maximum energy storage state0.90.9Self-consumption rate0.0010.01Minimum energy storage state0.20Operation and maintenance unit price (yuan/kW·h)0.0510.045Maximum charge and discharge power/kW5050Table 3
Pollutant discharge and treatment costs.
TypesSO2NOXCO2Gas turbine emission standards (g/kWh)0.0234.795170.16Emission standard for purchased electricity (g/kWh)6.42.32696Treatment costs of various pollutants (yuan/t)100019509.75


## 3.2. Optimal Operation Analysis

In order to verify the effectiveness of the IES source-storage-load coordination optimization model that takes into account the fine energy storage model in system operation scheduling and wind power consumption. The following four simulation scenarios are set up in the article for verification, of which scenario 4 is the detailed analysis model in the paper. Scenario 1, nonground source heat pump and energy storage, excluding IDR; scenario 2, including ground source heat pump and traditional energy storage model, excluding IDR; scenario 3, including ground source heat pump and fine energy storage model, excluding IDR; scenario 4, including ground source heat pump and fine energy storage model, taking IDR into account.Figure8 shows the electrical and thermal load demand curves of the system considering IDR. Compared with before the response, the peak-valley difference of electrical and thermal loads decreases to different degrees. Under the demand response mechanism, part of the peak load is transferred to the valley period, which makes the load curve smoother after coupling and makes the system have better effect of peak cutting and valley filling. It also improves the economy of IES operation.Figure 8
Electric and heat load curves before and after response.Taking the ground-source heat pump, the fine energy storage model, and the analysis of scenario 4 considering IDR as an example, the optimized output of each device in the system is shown in Figures9 and 10.Figure 9
Optimization result of power supply of the system.Figure 10
Results of system heating optimization.After considering IDR, part of the electricity and heat load during the peak period of the system is transferred to the valley period, which improves the wind power consumption rate and operating economy. At 1 : 00–7:00 and 23 : 00–24 : 00, the electricity price is in the lowest period. At this point, the power generation cost of the CCHP unit in the system is greater than the electricity price of the grid, so the power generation power of the CCHP unit “fixed by heat” can be reduced by reducing the indoor temperature to improve the wind power consumption capacity of the system. Ground source heat pump at this stage of the heating cost is relatively low, preferentially used for heating.At 7 : 00–8:00 and 11 : 00–18 : 00, the electricity price is moderate. At this time, the power generation cost of the system CCHP unit is still greater than the grid electricity price. When the electrical load demand is small, reduce the indoor temperature to reduce the system power generation cost; in the peak period of electric load demand, the room temperature is raised and the power of CCHP unit is increased, so as to reduce the cost of IES buying electricity from external grid in the whole dispatching cycle. At this time, the thermal power of the system is first provided by the CCHP unit. When there is a shortage of thermal power, the ground source heat pump is started for heating.At 8 : 00–11 : 00 and 18 : 00–23 : 00, the electricity price is at the highest value, and the system power generation cost is lower than the grid electricity price. From 18 : 00 to 23 : 00, the electricity/heat load demand is relatively high. By increasing the room temperature, the power generation of the CCHP unit “using heat” is maximized, and the system power purchase cost is reduced. In this period, the thermal power is completely provided by the CCHP unit, the excess thermal energy of the system is stored in the HST, and the ground source heat pump does not work.During the scheduling timeT, EES is driven by the peak-valley price mechanism to store energy during the peak period of electricity price and release energy during the peak period of electricity price and load, so as to reduce the operating cost of IES and cut peaks and fill valleys. The heat storage tank stores the excess heat energy in the system and provides energy when the heat load demand cannot be met; this reduces the operating cost of the system and further validates the economics of the model.


## 3.3. Analysis of Example Results

The operating cost and wind power consumption rate of IES are shown in Tables4 and 5. It can be seen from Table 4 that Scenarios 2 and 3 are compared with Scenario 1, in terms of abandonment cost, gas purchase cost, and total cost, which have significantly reduced, which proves the economy of the ground source heat pump and energy storage device for the system. Scenario 3, which considers the refined energy storage model, has slightly higher operating costs than Scenario 2. It is because the fine energy storage model considers the influence of environmental temperature and other constraints, which increases the operating cost. However, the model accurately describes the equipment output, making it more practical and schedulable, and its increased cost is within an acceptable range. In Scenario 4, compared with Scenario 1, the cost of gas purchase was reduced by 25.65%, the cost of pollutant emission was reduced by 20.82%, the cost of abandonment was reduced by 100%, and the total cost was reduced by 10.12%.Table 4
Operation cost of each scenario.
Cost category value (yuan/day)Scene 1Scene 2Scene 3Scene 4Purchase cost of gas57894347.604382.404362.10Operation and maintenance costs466.67447.37440.22434.42Pollutant treatment cost6755.0455.5653.05Power purchase cost454.05858.15872.62609.47Electricity sales revenue1005286.90295.51252.32Wind curtailment cost21.064.224.890Total cost5792.785425.485460.185206.72Table 5
Total amount and rate of wind power absorbed by the system.
ParameterScene 1Scene 2Scene 3Scene 4Total wind power consumption (kW)4592.555294.385255.065471Wind power consumption rate83.94%96.77%96.05%100%It can be seen from Table5 that after considering ground source heat pumps and energy storage equipment, the system’s wind power consumption rate has increased significantly, but electricity generated by the wind is still discarded at night; Scenario 4, after comprehensive consideration of IDR, ground source heat pumps, and energy storage equipment, the system. The wind power consumption rate reaches 100%, which solves the phenomenon of wind abandonment caused by the system’s “heat-fixed power” operation mode. Through the analysis, the results verify that the IES system with the fine energy storage model in the paper has the ability to completely absorb wind power and significant economic benefits. It provides corresponding solutions for the future wind power consumption.


## 4. Conclusion

In this paper, a cooperative optimal control strategy of cooling, heating, and power multitypes of energy is proposed. A refined energy storage model was proposed, which considered the uncertainty of wind power output and the influence of low temperature on the energy storage device. Through the simulation experiment on the model of cold, hot, and electric power supply system, the following conclusions are drawn: in the process of energy cooperative optimization control, the comprehensive consideration of electrical and thermal flexible load demand response can effectively reduce the load peak-valley difference, improve the wind power consumption capacity, and make the system have low economic benefits. At the same time, considering environmental temperature, capacity, and other influencing factors, the output mode is more in line with the reality, and the system optimization scheduling is more reliable. Finally, it is verified that the IES system with the fine energy storage model can absorb the wind power well and has certain economic benefits in dealing with the comprehensive energy optimization scheduling.



---

*Source: 1207430-2021-12-07.xml*