# Optimization Operation of the Park-Level Integrated Energy System Based on the Improved Coyote Optimization Algorithm

Authors: Jiyong Li; Yasai Wang; Yifei Zheng; Chengye Liu; Xiaoshuai Xu; Ran Chen
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

Conventional optimization methods cannot fully satisfy the interests of multiparticipants and protect the privacy of participants in the integrated energy system that observe changes in the energy market structure. To allocate the benefits among the stakeholders in the integrated energy system and improve renewable energy accommodation, the manuscript proposes an optimal dispatching strategy for a park-level integrated energy system employing the Stackelberg game. Firstly, the benefits and cost models of each stakeholder of the integrated energy system are constructed by considering the integrated demand response and the uncertainty of renewable energy output. A master-slave game model that contains the energy system operator, energy producer, and energy users is then established, and the existence of the Stackelberg equilibrium is demonstrated. Furthermore, a distributed algorithm is proposed to resolve the game model by combining an improved coyote optimization algorithm with quadratic programming. Due to the shortcomings of the conventional coyote optimization algorithm, such as slow convergence rate and quickly falling into local optimum, a beetle antennae search is utilized to strengthen the optimal and the worst coyotes and to improve the convergence speed, global search ability, and optimization accuracy of the standard coyote algorithm. Finally, an industrial park in Northern China is adopted as an illustration to evaluate the effectiveness of the model and the improved algorithm.


---

## Body

## 1. Introduction

The over-exploitation of fossil fuels like coal and oil has led to severe problems, such as environmental pollution and global warming [1, 2]. Increasing the accommodation capacity of cleaner energy such as photovoltaics is crucial to the reductions in emissions of greenhouse gases and the sustainable development of human society [3–5]. However, clean energy combined with conventional energy networks can result in a large amount of abandonment of wind power and photovoltaics, which is a great waste of resources and does not meet the current energy demand well [6]. In this context, the energy Internet with multienergy coupling and also the combination of information and energy technologies will become the mainstream mode of energy supply, and several national and international energy organizations attempt to promote the development of the energy Internet [7]. Since the park-level integrated energy system (PIES) serves as a basis of the energy Internet, its optimization draws a lot of interest in current energy research. The PIES integrates the combined cooling heating power (CCHP), energy storage devices, renewable energy devices, and other energy conversion equipment to realize energy gradient utilization and a high-proportion renewable energy accommodation by coordinating and optimizing energy generation, transmission, storage, and consumption.The authors have done a lot of research on the PIES and have constructed various energy systems based on the PIES and resolved the constructed system models employing feasible strategies. Various studies have been performed on combined cooling, heating, and power-type integrated energy systems such as Refs. [8–11]. Xu et al. [9] established an optimization model of combined cooling, heating, and power-type multimicrogrid by considering the interaction power of microgrids. Wang et al. [10] constructed a CCHP-type integrated energy system by combining solar and compressed air energy storage and employed a nondominated sorting genetic algorithm to design a multiobjective optimal operation strategy. Ju et al. [11] constructed a hybrid energy system of the CHP and renewable energy driven by distributed energy and constructed a multiobjective operation optimization and evaluation model.However, the mentioned research has employed a centralized optimization approach, which could not appropriately describe the interaction between electricity prices and loads and protect the privacy and security of the participants. There exist multiple agents in the integrated energy system with their demands. There may also be conflicts between the different agents. The game theory has been widely utilized in the optimization operations of integrated energy systems to satisfy the interests of multiple agents. At the same time, the game process can also well describe the interaction between electricity price and load. In Ref. [12], the energy hub model was utilized to construct an optimization model of the coalition game, considering the integrated demand response. A multileader multifollower game model was constructed in Ref. [13] to resolve the equilibrium strategies between multiple distributed energy stations and multiple energy consumers. Ma et al. [14] constructed a master-slave transaction model consisting of operators and PV prosumers based on the approach called distributed energy management. Wang et al. [15] considered a scenario with a variable called the integrated price of electricity of energy retailers and heat and established a Stackelberg game optimization framework, including the integrated energy retailer, energy supplier, and load aggregator for cooperative optimization of multiple agents. However, the aforementioned studies only cover the source, grid, and load segments. They do not consider the impact of energy storage devices on the transaction and the uncertainty of both photovoltaic and wind power outputs.In addition, many feasible schemes are proposed to resolve the PIES optimal scheduling model. Hou et al. [16] employed an improved particle swarm algorithm to optimize an integrated electro-thermal-hydrogen energy system by considering the uncertainty of the renewable energy output. Binary particle swarm algorithms and particle swarm algorithms were adopted in Ref. [17] to resolve the multiobjective unit commitment problem caused by binary and real variables in the model, respectively. Delice et al. [18] proposed a new modified particle swarm optimization algorithm with negative knowledge to resolve the mixed-model two-sided assembly line balancing problem. An algorithm [19] was proposed that combines an adaptive gravitational search algorithm (AGSA) with pattern search (PS) called AGSA-PS. Koessler and Almomani [20] tested three methods of hybridizing the PSO and the PS to enhance the global minima and robustness. A hybrid optimization method [21], the GA-SQP, in which the genetic algorithm (GA) is a stochastic method was combined with the sequential quadratic programming (SQP) method, which was a deterministic method. An algorithm [22] was proposed that consists of the combination of both genetic algorithm (GA) and the particle swarm optimization (PSO). A new metaheuristic method [23] was suggested that was inspired by the life cycle of water striders. A new approach to the firefly algorithm [24] was suggested that was based on opposition-based learning (OBFA) to enhance the global search ability of the original algorithm. As a novel swarm intelligence optimization algorithm, the coyote optimization algorithm (COA) was proposed [25], and testing functions verified its superiority over the other algorithms. However, due to the large-scale, nonlinear, and nonconvex characteristics of the optimal scheduling model of the PIES, the COA still suffers from slow convergence speed, low convergence accuracy, and other shortcomings in the optimization process of the actual PIES.Facing the above issues, the manuscript established an operational optimization strategy for multiple competing agents under the master-slave game framework for a CCHP-type PIES. At first, the probability scenarios method was introduced to describe the wind speed, solar radiation, and load uncertainty, while the scenario set of the wind, photovoltaics, and load was derived by utilizing the scenario reduction techniques. Next, the revenue model of each agent of PIES was established under the master-slave game framework, and the existence of the equilibrium solutions in the transaction game is proved. Besides, a distributed solution algorithm was proposed to protect the privacy of each agent, and the improved COA was presented. The standard COA and the improved one were utilized to obtain the profit of the energy system operator, thus enhancing the revenue of the operator.The rest of this article is organized as follows. Section2 introduces the structure of the PIES system and the energy transaction process and establishes the mathematical model of each agent. Section 3 constructs the Stackelberg game model and proves the existence of the Stackelberg equilibrium for the constructed game model. The implementation method is introduced in Section 4. Case analysis and conclusions are presented in Section 5 and Section 6, respectively.


## 2. The PIES Multiagent Model

The proposed PIES consists of an energy system operator (ESO), energy user (EU), and energy producer (EP) to form a complete system with the operator that functions as a bridge for transactions between the multiple agents. The schematic diagram of the Stackelberg game of the PIES is shown in Figure1.Figure 1
The schematic diagram of the Stackelberg game.


### 2.1. The Description of the System

The energy producer with the CCHP system is the primary energy source in the integrated energy system providing electrical and heat energy to the operator and determining the output of each device based on the energy offered given by the operator. As the load side of the system, the energy users buy both electrical and heat energy from the operator and adjust the electrical and thermal consumptions based on the offer of the ESO.The ESO connects the source and the load and acts as a bridge for energy interaction between the EU and the EP. The ESO sets energy prices based on supply and demand, thus maximizing revenue by buying at lower prices and selling electric and thermal energy at higher prices. The ESO has a more flexible pricing mechanism than that of the conventional grid, encouraging consumers to participate in demand response. The ESO suffers from inadequate energy supply, purchasing power from the grid at a higher price when the sum of the electric power supplied from the EP is still less than the electric demand of the EU and compensating for the EU when the ESO cannot meet the thermal load.The PIES structure is shown in Figure2, where the natural gas flowing into the system is converted into both electrical and thermal energies by the energy production equipment and into various energy supplies to the load side through the energy conversion equipment. The gas turbine generates both electrical and thermal energies, where electrical energy is supplied to the electrical load along with electricity generated by wind and photovoltaics. In contrast, thermal energy is recovered through the waste heat boiler and then aggregated with thermal energy produced by the boiler, a part of which is supplied to the thermal load through the heat exchanger, and the rest is supplied to the cold load through the chiller. The thermal storage tank (TST) and batteries store the excess energy.Figure 2
The schematic diagram of the PIES.


### 2.2. Model of the Energy System Operator

The ESO, as the upper level, maximizes its revenue by considering the supply and demand for electricity and heat energy and market information to determine the price of buying and selling energy, buying the energy from the EP and the ESP, and selling it to the EU and the ESP. The ESO aims to maximize its profit by(1)maxFeso=∑t=124IESOsell−CESOcostΔt,IESOsell=peselltPeut+Pest+phselltPhut+Phst,CESOcost=CESObuy+CESOgrid+CESOpenal,CESObuy=pebuytPept+Peyt+phbuytPhpt+Phyt,CESOgrid=pgbuytPgbuyt−pgselltPgsellt,CESOpenal=Phtξht,where IESOsell represents the revenue of the operator from the sold energy; CESOcost represents the operating cost of the operator; CESObuy, CESOgrid, and CESOpenal stand for the energy purchasing costs of the ESO from the EP, the electricity purchasing cost from the grid, and the penalty cost of the interrupting heat supply, respectively; pebuy and phbuy represent the prices which the ESO pays for electricity and heat; pesell and phsell denote the prices obtained by selling electricity and heat by the ESO, respectively; pgbuy and pgsell represent the time-of-use electricity price and feed-in tariffs, respectively; Pes and Phs denote the electrical and heat power sold by the ESO to the ESP; Peu and Phu represent the electrical and heat power sold by the ESO to the EU, respectively; Pep and Php denote the electrical and heat power purchased by ESO from the EP, respectively; Pey and Phy represent the electrical and heat power purchased by the ESO from the the ESP, respectively; Ph and ξh denote the amount of heat load loss and the price of heat interruption penalty, respectively.To ensure that entities will not skip the ESO to trade or trade with the outside world directly, some constraints should be implied on the offers of the ESO to ensure that the buying and selling prices are within the market price range for the electrical and thermal energies. These constraints are given by(2)pgsellt<pesellt<pgbuyt,pgsellt<pebuyt<pgbuyt,phmint<phsellt<phmaxt,phmint<phbuyt<phmaxt,where phmax and phmin indicate the upper and lower limits of the heating price, respectively.


### 2.3. Model of the Energy Producer

The EP maximizes its revenue by adjusting the output of each device in the CCHP system after the ESO gives the energy offer, and its profit is the difference between the revenue from the energy salesIEP, the generation cost Cgen, and the emission cost Cemi as follows:(3)maxFEP=∑t=124IEP−Cgen−CemiΔt,IEP=pebuytPept+phbuytPhpt,Cgen=pgasVgt+∑i=1IJiPit,Cemi=PemiFgentvC+vN,where Cgen represents the generation cost of the CCHP system; pgas denotes the price of the natural gas; Vg represents the volume of the natural gas purchased; Ji denotes the maintenance cost factor of the ith CCHP device; Pi denotes the output power of the ith device; I denote the total number of devices in the CCHP; Pemi represents the gas emission factor; Fgen denotes the heat energy of the consumed gas; vC and vN denote the carbon dioxide and nitrogen oxide volumes of emissions.The electrical and heat energy sold by the EP can be provided jointly by the CCHP devices as follows:(4)Pept=Ppvt+Pwtt+PMTt,Phpt=PGBt+PREt,where Ppv and Pwt represent the power generated by photovoltaics (PV) and wind turbine (WT), respectively; PMT represents the power generated by the gas turbine; PGB and PRE denote the thermal output powers of the gas boiler and the waste heat boiler, respectively.The power relationship for the CCHP system is described by(5)PRE=1−ηMTηMTPMTηRE,where ηMT denotes the power generation efficiency of the gas turbine and ηRE represents the heat production efficiency of the waste heat boiler.The balance of the electrical power is expressed by(6)Ppvt+Pwtt+PMTt+Pnett+Pey=lfet+ltet,where Pnet denotes the amount of power trading between the CHP system and the grid, while positive denotes electricity sales and negative represents electricity purchase; Pey represents the power of the energy storage, while positive denotes discharging, negative represents charging; and lfe and lte represent the fixed electric load and shifted electrical load, respectively.In the actual operation, the following constraints should be satisfied to fulfill the constraints on the power output and creep rate of the gas turbine and the gas boiler at timet:(7)0≤PMTt≤PMT,max,0≤PGBt≤PGB,max,PMT,down≤PMTt−PMTt−1≤PMT,up,PGB,down≤PGBt−PGBt−1≤PGB,up,where PMT, max and PGB,max denote the rated capacity of the gas turbine and the gas boiler, respectively; PMT,up and PGB,up represent the upper limits of the ramp rate of the gas turbine and the boiler; and PMT,down and PGB,down denote the lower limits of the ramp rate of the gas turbine and the boiler, respectively.The energy storage device should satisfy some constraints, such as the thermal storage tank constraint, described as(8)Htstt=Htstt−11−ηh+ηchhPhyt−Phstηdish,Htstmin≤Htstt≤Htstmax,Pchmin≤Phst−Phst−1≤Pchmax,Pdismin≤Phyt−Phyt−1≤Pdismax,where Htst (t) represents the real-time quantity of heat storage; η h ch and η h dis denote the charging and discharging efficiencies of the thermal energy; Htstmin and Htstmax represent the minimum and maximum heat storage capacity, respectively; and ηh denotes the energy self-loss rate of the thermal storage tank. The energy storage equipment should ensure the consistency of the reserves at the beginning and the end of the daily cycle.The method of a probability-based scenario is adopted in the manuscript to describe the stochastic volatility of the wind and the photovoltaic output considering the substantial uncertainty in solar radiation and wind speed [26]. At first, probability scenarios are generated. The uncertainty sampling of the wind speed, photovoltaics, and load is then completed employing Latin hypercube, while the Weibull and Beta distributions are utilized for the wind speed and the light intensity, respectively. Various sampling scenarios are generated after the sampling is completed. Finally, scenario reduction techniques are employed to complete the scenario reduction and obtain the probability of occurrence for each scenario.The expression for the output of the photovoltaic power system is represented by(9)PPV,t=PSTCGtGSTC1+kTt−TSTC,Tt=Tairt+0.01381+0.031Tairt1−0.042VWGt,Tairt=12Tmax−Tmin+Tmax−Tminsin2πt−tp24,where PPV,t represents the output power of the PV power system at time t; PSTC denotes the maximum output power of solar panel under the standard test; GSTC represents light intensity under the standard test, which is 1000 W/m2; k denotes the temperature coefficient; T(t) represents the actual temperature of the solar panel at time t; TSTC denotes the solar panel temperature under the standard test, which is 25°C; Tair represents the outdoor temperature at time t; VW denotes the wind speed; Tmax denotes the maximum daily temperature; Tmin represents the minimum daily temperature; and tp represents the average daily temperature.Light intensity usually follows a beta distribution. The probability density function is defined by(10)fGt=Γα+βΓαΓβGtGmaxα−11−GtGmaxβ−1,where G(t) denotes the light intensity at time t; Γ represents the gamma probability density function; α and β denote the shape factor of the beta probability density function; and Gmax denotes the maximum daily light intensity.When the installed capacity is determined, the maximum value of wind power output at each moment is determined by the actual conditions such as weather and environment. The output of the wind at timet can be expressed as a function of the intermittent wind speed as(11)PWPP,t=0,0≤νt≤νin,νt>νout,νt−νinνrated−νinPrated,νin≤νt<νrated,Prated,νrated≤νt≤νout,where νt denotes the real-time wind speed; νin represents the cut-in wind speed; νout represents the cut-out wind speed; νrated denotes the rated wind speed; and Prated represents the rated power of the wind turbine.The natural incoming windν usually follows the Weibull probability distribution. The probability density function is defined by(12)fwν=φθνθφ−1e−ν/θφ,where φ denotes the shape parameter and θ represents the scale parameter.


### 2.4. Model of Energy Users

In a real integrated energy system, there exists a stochastic nature of real-time load changes with time and season. In this paper, a short-term load forecasting method based on the robust Holt–Winter model is used to forecast the load. The method integrates the time series characteristics of a linear trend, seasonal variation, and stochastic fluctuation and combines with the exponential smoothing method to have better forecasting capability.The basic idea is to decompose the time series with linear trend, seasonal variation, and stochastic fluctuation and combine them with the exponential smoothing method to estimate the long-term trend, increment of trend, and seasonal fluctuation respectively. A forecasting model is then built and extrapolated to the forecast values.The multiplicative Holt–Winter model consists of three smoothing equations and one forecasting equation. The equation system of the parameters is defined by(13)at=αytSt−1+1−αat−1+bt−1,bt=βat−at−1+1−βbt−1,St=γytat+1−γSt−1.The forecasting equation is defined by(14)y^t+k=at+kbtSt−l+k,where yt represents the observed value of the time series at time t; at denotes the stable component of the load data at time t; bt represents the tendency trend component of the load data at time t; St denotes the seasonal component of the load data at time t; l represents the seasonal length; α, β, γ∈[0,1] denote the smoothing parameter; and k represents the number of moments to be predicted.The EU adjusts the energy use plan according to the ESO offer, and the objective function is the difference between the utility of the user and the cost of energy purchase which is defined by(15)maxFEU=UEU−CEUbuy.UEU denotes the utility function of the user, which is the sum of the satisfaction gained by the user from the consumed electricity and heat energy [27], which is expressed as follows:(16)UEU=vedet−ue2det2+vhdht−uh2dht2,where ve, ue, vh, and uh represent the preference constants for the widely used quadratic utility function [28–30]. de and dh denote the electrical and thermal load demands of the EU.CEUbuy denotes the cost of energy purchase by the user as follows:(17)CEUbuy=peselltPeut+phselltPhut.The demand response of the EU is divided into shifted electrical loadlte and curtailable thermal load lrh, while the fixed load of the EU is divided into fixed electrical load lfe and fixed thermal load lfh, which is formulated as(18)det=lfet+ltet,dht=lfht−lrht.lte and lrh should satisfy the following constraints as follows:(19)0≤ltet≤ltemax,0≤lrht≤lrhmax,where ltemax and lrhmax present the upper limits of the response volume set to 20% of the total load.


## 2.1. The Description of the System

The energy producer with the CCHP system is the primary energy source in the integrated energy system providing electrical and heat energy to the operator and determining the output of each device based on the energy offered given by the operator. As the load side of the system, the energy users buy both electrical and heat energy from the operator and adjust the electrical and thermal consumptions based on the offer of the ESO.The ESO connects the source and the load and acts as a bridge for energy interaction between the EU and the EP. The ESO sets energy prices based on supply and demand, thus maximizing revenue by buying at lower prices and selling electric and thermal energy at higher prices. The ESO has a more flexible pricing mechanism than that of the conventional grid, encouraging consumers to participate in demand response. The ESO suffers from inadequate energy supply, purchasing power from the grid at a higher price when the sum of the electric power supplied from the EP is still less than the electric demand of the EU and compensating for the EU when the ESO cannot meet the thermal load.The PIES structure is shown in Figure2, where the natural gas flowing into the system is converted into both electrical and thermal energies by the energy production equipment and into various energy supplies to the load side through the energy conversion equipment. The gas turbine generates both electrical and thermal energies, where electrical energy is supplied to the electrical load along with electricity generated by wind and photovoltaics. In contrast, thermal energy is recovered through the waste heat boiler and then aggregated with thermal energy produced by the boiler, a part of which is supplied to the thermal load through the heat exchanger, and the rest is supplied to the cold load through the chiller. The thermal storage tank (TST) and batteries store the excess energy.Figure 2
The schematic diagram of the PIES.


## 2.2. Model of the Energy System Operator

The ESO, as the upper level, maximizes its revenue by considering the supply and demand for electricity and heat energy and market information to determine the price of buying and selling energy, buying the energy from the EP and the ESP, and selling it to the EU and the ESP. The ESO aims to maximize its profit by(1)maxFeso=∑t=124IESOsell−CESOcostΔt,IESOsell=peselltPeut+Pest+phselltPhut+Phst,CESOcost=CESObuy+CESOgrid+CESOpenal,CESObuy=pebuytPept+Peyt+phbuytPhpt+Phyt,CESOgrid=pgbuytPgbuyt−pgselltPgsellt,CESOpenal=Phtξht,where IESOsell represents the revenue of the operator from the sold energy; CESOcost represents the operating cost of the operator; CESObuy, CESOgrid, and CESOpenal stand for the energy purchasing costs of the ESO from the EP, the electricity purchasing cost from the grid, and the penalty cost of the interrupting heat supply, respectively; pebuy and phbuy represent the prices which the ESO pays for electricity and heat; pesell and phsell denote the prices obtained by selling electricity and heat by the ESO, respectively; pgbuy and pgsell represent the time-of-use electricity price and feed-in tariffs, respectively; Pes and Phs denote the electrical and heat power sold by the ESO to the ESP; Peu and Phu represent the electrical and heat power sold by the ESO to the EU, respectively; Pep and Php denote the electrical and heat power purchased by ESO from the EP, respectively; Pey and Phy represent the electrical and heat power purchased by the ESO from the the ESP, respectively; Ph and ξh denote the amount of heat load loss and the price of heat interruption penalty, respectively.To ensure that entities will not skip the ESO to trade or trade with the outside world directly, some constraints should be implied on the offers of the ESO to ensure that the buying and selling prices are within the market price range for the electrical and thermal energies. These constraints are given by(2)pgsellt<pesellt<pgbuyt,pgsellt<pebuyt<pgbuyt,phmint<phsellt<phmaxt,phmint<phbuyt<phmaxt,where phmax and phmin indicate the upper and lower limits of the heating price, respectively.


## 2.3. Model of the Energy Producer

The EP maximizes its revenue by adjusting the output of each device in the CCHP system after the ESO gives the energy offer, and its profit is the difference between the revenue from the energy salesIEP, the generation cost Cgen, and the emission cost Cemi as follows:(3)maxFEP=∑t=124IEP−Cgen−CemiΔt,IEP=pebuytPept+phbuytPhpt,Cgen=pgasVgt+∑i=1IJiPit,Cemi=PemiFgentvC+vN,where Cgen represents the generation cost of the CCHP system; pgas denotes the price of the natural gas; Vg represents the volume of the natural gas purchased; Ji denotes the maintenance cost factor of the ith CCHP device; Pi denotes the output power of the ith device; I denote the total number of devices in the CCHP; Pemi represents the gas emission factor; Fgen denotes the heat energy of the consumed gas; vC and vN denote the carbon dioxide and nitrogen oxide volumes of emissions.The electrical and heat energy sold by the EP can be provided jointly by the CCHP devices as follows:(4)Pept=Ppvt+Pwtt+PMTt,Phpt=PGBt+PREt,where Ppv and Pwt represent the power generated by photovoltaics (PV) and wind turbine (WT), respectively; PMT represents the power generated by the gas turbine; PGB and PRE denote the thermal output powers of the gas boiler and the waste heat boiler, respectively.The power relationship for the CCHP system is described by(5)PRE=1−ηMTηMTPMTηRE,where ηMT denotes the power generation efficiency of the gas turbine and ηRE represents the heat production efficiency of the waste heat boiler.The balance of the electrical power is expressed by(6)Ppvt+Pwtt+PMTt+Pnett+Pey=lfet+ltet,where Pnet denotes the amount of power trading between the CHP system and the grid, while positive denotes electricity sales and negative represents electricity purchase; Pey represents the power of the energy storage, while positive denotes discharging, negative represents charging; and lfe and lte represent the fixed electric load and shifted electrical load, respectively.In the actual operation, the following constraints should be satisfied to fulfill the constraints on the power output and creep rate of the gas turbine and the gas boiler at timet:(7)0≤PMTt≤PMT,max,0≤PGBt≤PGB,max,PMT,down≤PMTt−PMTt−1≤PMT,up,PGB,down≤PGBt−PGBt−1≤PGB,up,where PMT, max and PGB,max denote the rated capacity of the gas turbine and the gas boiler, respectively; PMT,up and PGB,up represent the upper limits of the ramp rate of the gas turbine and the boiler; and PMT,down and PGB,down denote the lower limits of the ramp rate of the gas turbine and the boiler, respectively.The energy storage device should satisfy some constraints, such as the thermal storage tank constraint, described as(8)Htstt=Htstt−11−ηh+ηchhPhyt−Phstηdish,Htstmin≤Htstt≤Htstmax,Pchmin≤Phst−Phst−1≤Pchmax,Pdismin≤Phyt−Phyt−1≤Pdismax,where Htst (t) represents the real-time quantity of heat storage; η h ch and η h dis denote the charging and discharging efficiencies of the thermal energy; Htstmin and Htstmax represent the minimum and maximum heat storage capacity, respectively; and ηh denotes the energy self-loss rate of the thermal storage tank. The energy storage equipment should ensure the consistency of the reserves at the beginning and the end of the daily cycle.The method of a probability-based scenario is adopted in the manuscript to describe the stochastic volatility of the wind and the photovoltaic output considering the substantial uncertainty in solar radiation and wind speed [26]. At first, probability scenarios are generated. The uncertainty sampling of the wind speed, photovoltaics, and load is then completed employing Latin hypercube, while the Weibull and Beta distributions are utilized for the wind speed and the light intensity, respectively. Various sampling scenarios are generated after the sampling is completed. Finally, scenario reduction techniques are employed to complete the scenario reduction and obtain the probability of occurrence for each scenario.The expression for the output of the photovoltaic power system is represented by(9)PPV,t=PSTCGtGSTC1+kTt−TSTC,Tt=Tairt+0.01381+0.031Tairt1−0.042VWGt,Tairt=12Tmax−Tmin+Tmax−Tminsin2πt−tp24,where PPV,t represents the output power of the PV power system at time t; PSTC denotes the maximum output power of solar panel under the standard test; GSTC represents light intensity under the standard test, which is 1000 W/m2; k denotes the temperature coefficient; T(t) represents the actual temperature of the solar panel at time t; TSTC denotes the solar panel temperature under the standard test, which is 25°C; Tair represents the outdoor temperature at time t; VW denotes the wind speed; Tmax denotes the maximum daily temperature; Tmin represents the minimum daily temperature; and tp represents the average daily temperature.Light intensity usually follows a beta distribution. The probability density function is defined by(10)fGt=Γα+βΓαΓβGtGmaxα−11−GtGmaxβ−1,where G(t) denotes the light intensity at time t; Γ represents the gamma probability density function; α and β denote the shape factor of the beta probability density function; and Gmax denotes the maximum daily light intensity.When the installed capacity is determined, the maximum value of wind power output at each moment is determined by the actual conditions such as weather and environment. The output of the wind at timet can be expressed as a function of the intermittent wind speed as(11)PWPP,t=0,0≤νt≤νin,νt>νout,νt−νinνrated−νinPrated,νin≤νt<νrated,Prated,νrated≤νt≤νout,where νt denotes the real-time wind speed; νin represents the cut-in wind speed; νout represents the cut-out wind speed; νrated denotes the rated wind speed; and Prated represents the rated power of the wind turbine.The natural incoming windν usually follows the Weibull probability distribution. The probability density function is defined by(12)fwν=φθνθφ−1e−ν/θφ,where φ denotes the shape parameter and θ represents the scale parameter.


## 2.4. Model of Energy Users

In a real integrated energy system, there exists a stochastic nature of real-time load changes with time and season. In this paper, a short-term load forecasting method based on the robust Holt–Winter model is used to forecast the load. The method integrates the time series characteristics of a linear trend, seasonal variation, and stochastic fluctuation and combines with the exponential smoothing method to have better forecasting capability.The basic idea is to decompose the time series with linear trend, seasonal variation, and stochastic fluctuation and combine them with the exponential smoothing method to estimate the long-term trend, increment of trend, and seasonal fluctuation respectively. A forecasting model is then built and extrapolated to the forecast values.The multiplicative Holt–Winter model consists of three smoothing equations and one forecasting equation. The equation system of the parameters is defined by(13)at=αytSt−1+1−αat−1+bt−1,bt=βat−at−1+1−βbt−1,St=γytat+1−γSt−1.The forecasting equation is defined by(14)y^t+k=at+kbtSt−l+k,where yt represents the observed value of the time series at time t; at denotes the stable component of the load data at time t; bt represents the tendency trend component of the load data at time t; St denotes the seasonal component of the load data at time t; l represents the seasonal length; α, β, γ∈[0,1] denote the smoothing parameter; and k represents the number of moments to be predicted.The EU adjusts the energy use plan according to the ESO offer, and the objective function is the difference between the utility of the user and the cost of energy purchase which is defined by(15)maxFEU=UEU−CEUbuy.UEU denotes the utility function of the user, which is the sum of the satisfaction gained by the user from the consumed electricity and heat energy [27], which is expressed as follows:(16)UEU=vedet−ue2det2+vhdht−uh2dht2,where ve, ue, vh, and uh represent the preference constants for the widely used quadratic utility function [28–30]. de and dh denote the electrical and thermal load demands of the EU.CEUbuy denotes the cost of energy purchase by the user as follows:(17)CEUbuy=peselltPeut+phselltPhut.The demand response of the EU is divided into shifted electrical loadlte and curtailable thermal load lrh, while the fixed load of the EU is divided into fixed electrical load lfe and fixed thermal load lfh, which is formulated as(18)det=lfet+ltet,dht=lfht−lrht.lte and lrh should satisfy the following constraints as follows:(19)0≤ltet≤ltemax,0≤lrht≤lrhmax,where ltemax and lrhmax present the upper limits of the response volume set to 20% of the total load.


## 3. Stackelberg Game Model

This section establishes a Stackelberg game model to analyze the trading of multiple energies. Besides, the existence of the game equilibrium is proved.


### 3.1. Basic Elements of the Game

According to the above section, the ESO, the EP, and the EU are all independent agents in the PIES. The ESO sets the price strategy to obtain maximum revenue, while the EP and EU make optimal adjustments according to the price signals of the ESO, influencing the price setting of the ESO. The decisions among the agents are sequential and affect each other. Since the ESO is a manager with the decision priority, the game of three agents can be constructed as a Stackelberg game [31].The master-slave game model is given by(20)Y=E∪P∪U;L;F.The participant set, strategy set, and payoff function as the three elements of the game model described in (20) can be described as(1)
The game participants include energy system operatorE, energy user U, and energy producer P(2)
The strategies set is defined asL = {Le, Lu, Lp}, where the strategy of the ESO, Le, represents the offer of energy, expressed as Le=(p buy e, p buy h, p sell e, p sell h); the strategy of the EP, Lp, denotes the output of gas turbines and boilers at each period, denoted as Lp = (PMT, PGB); the strategy of the user, Lu, represents the amount of demand response at each moment, denoted as Lu = (lte, lrh)(3)
The payoff setF = {FESO, FEP, FEU} is the objective function of the four agents, calculated by equations (1), (3), and (15)


### 3.2. Stackelberg Game Equilibrium

When all the followers respond optimally to the prior strategy of the leader and the leader accepts this response, the two-level game reaches Stackelberg equilibrium [13]. Let Le∗ be the equilibrium strategy of the ESO and Lu∗ and Lp∗ be the optimal response strategies of the EU and EP, respectively. The set of strategies (Le∗, Lu∗, and Lp∗) is the equilibrium solution of the Stackelberg game under the following constraints as follows:(21)FesoLe∗,Lu∗,Lp∗≥FesoLe−i∗,Lei,Lu∗,Lp∗,FEULe∗,Lu∗,Lp∗≥FEULe∗,Lu,Lp∗,FEPLe∗,Lu∗,Lp∗≥FEPLe∗,Lu∗,LP,∀LC∈Le,where L∗e(−i) stands for the other strategies except Lei. When the strategies of each agent in the game are equilibrium solutions, no agent on either side can increase its profit by adjusting its set of strategies individually [13].The following conditions should be satisfied to ensure the existence of a Stackelberg equilibrium [32]:(1)
The strategy setsLe, Lu, and Lp are nonempty, bounded, and convex sets in the Euclidean space(2)
FEP and FEU represent quasi-concave functions concerning Lp and Lu, respectively(3)
FESO denotes a continuous function of Le, Lu, Ls, and Lp(4)
FEP and FEU denote continuous functions of Lp, and Lu, respectivelyAccording to the above PIES model, the strategy of the ESO and the strategy of the EP should satisfy constraints equations (2) and (7), respectively, while the strategy of the EU should satisfy constraints (18) and (19). Accordingly, the set of strategies of each agent satisfies condition (1).It can be concluded from equations (15)–(17) that the utility function of the EU is convex, and the other terms are linear or constant functions about lte and lrh. Thus, FEU represents a convex function of Lu. Similarly, the objective functions of EP are convex. Moreover, since the convex function must be a quasi-concave function, condition (2) is also satisfied.FESO, FEP, and FEU can be calculated from equations (1), (3), and (15) respectively. Therefore, it can be concluded that the four objective functions are continuous. Accordingly, conditions (3) and (4) are also satisfied.In summary, the existence of the equilibrium solution of the Stackelberg game is proved.


## 3.1. Basic Elements of the Game

According to the above section, the ESO, the EP, and the EU are all independent agents in the PIES. The ESO sets the price strategy to obtain maximum revenue, while the EP and EU make optimal adjustments according to the price signals of the ESO, influencing the price setting of the ESO. The decisions among the agents are sequential and affect each other. Since the ESO is a manager with the decision priority, the game of three agents can be constructed as a Stackelberg game [31].The master-slave game model is given by(20)Y=E∪P∪U;L;F.The participant set, strategy set, and payoff function as the three elements of the game model described in (20) can be described as(1)
The game participants include energy system operatorE, energy user U, and energy producer P(2)
The strategies set is defined asL = {Le, Lu, Lp}, where the strategy of the ESO, Le, represents the offer of energy, expressed as Le=(p buy e, p buy h, p sell e, p sell h); the strategy of the EP, Lp, denotes the output of gas turbines and boilers at each period, denoted as Lp = (PMT, PGB); the strategy of the user, Lu, represents the amount of demand response at each moment, denoted as Lu = (lte, lrh)(3)
The payoff setF = {FESO, FEP, FEU} is the objective function of the four agents, calculated by equations (1), (3), and (15)


## 3.2. Stackelberg Game Equilibrium

When all the followers respond optimally to the prior strategy of the leader and the leader accepts this response, the two-level game reaches Stackelberg equilibrium [13]. Let Le∗ be the equilibrium strategy of the ESO and Lu∗ and Lp∗ be the optimal response strategies of the EU and EP, respectively. The set of strategies (Le∗, Lu∗, and Lp∗) is the equilibrium solution of the Stackelberg game under the following constraints as follows:(21)FesoLe∗,Lu∗,Lp∗≥FesoLe−i∗,Lei,Lu∗,Lp∗,FEULe∗,Lu∗,Lp∗≥FEULe∗,Lu,Lp∗,FEPLe∗,Lu∗,Lp∗≥FEPLe∗,Lu∗,LP,∀LC∈Le,where L∗e(−i) stands for the other strategies except Lei. When the strategies of each agent in the game are equilibrium solutions, no agent on either side can increase its profit by adjusting its set of strategies individually [13].The following conditions should be satisfied to ensure the existence of a Stackelberg equilibrium [32]:(1)
The strategy setsLe, Lu, and Lp are nonempty, bounded, and convex sets in the Euclidean space(2)
FEP and FEU represent quasi-concave functions concerning Lp and Lu, respectively(3)
FESO denotes a continuous function of Le, Lu, Ls, and Lp(4)
FEP and FEU denote continuous functions of Lp, and Lu, respectivelyAccording to the above PIES model, the strategy of the ESO and the strategy of the EP should satisfy constraints equations (2) and (7), respectively, while the strategy of the EU should satisfy constraints (18) and (19). Accordingly, the set of strategies of each agent satisfies condition (1).It can be concluded from equations (15)–(17) that the utility function of the EU is convex, and the other terms are linear or constant functions about lte and lrh. Thus, FEU represents a convex function of Lu. Similarly, the objective functions of EP are convex. Moreover, since the convex function must be a quasi-concave function, condition (2) is also satisfied.FESO, FEP, and FEU can be calculated from equations (1), (3), and (15) respectively. Therefore, it can be concluded that the four objective functions are continuous. Accordingly, conditions (3) and (4) are also satisfied.In summary, the existence of the equilibrium solution of the Stackelberg game is proved.


## 4. The Solution of the Stackelberg Game Model Based on the Improved Coyote Optimization Algorithm

A large-scale nonlinear programming problem should be resolved to optimize the objective function of the leader ESO, while the improved coyote optimization algorithm (beetle antennae search coyote optimization algorithm, the BCOA) can be employed to reduce the solution complexity. As for the lower level of the model, since the objective function contains quadratic terms, quadratic programming can resolve the problem.


### 4.1. The Coyote Optimization Algorithm

The main difference between the COA and the most intelligent algorithms is that coyotes are divided into several groups, and the internal social influences are considered [25]. The COA should only set some control parameters, including the number of coyotes Np, the number of coyotes in each group Nc, the population of coyotes N, and the maximum number of iterations MaxIter.


#### 4.1.1. Initialization of Coyote Populations

The coyote populationN consists of Np groups of coyotes with Nc coyotes in each group, while the coyotes can be initialized in the search space [lb, ub] through the following equation:(22)socc,jp,t=lbj+rjubj−lbj,where c = [1, 2, ...Nc], p = [1, 2, ...Np], and j = [1, 2, ..., D], while D is the dimension of the optimization problem, and rj is a random number generated by a uniform probability distribution within [0, 1]. socp,t c,j is randomly initialized for the jth dimension of the cth coyote in the pth group, and lbj and ubj denote the lower and upper bounds of the jth dimension of the coyote, respectively.


#### 4.1.2. Coyote Growth in the Group

The adaptive ability of coyotes can be evaluated according to the objective function defined by(23)fitcp,t=fsoccp,t.Naturally, the alpha coyote is the best socially adapted coyote. In the COA, it corresponds to the best (minimum or maximum) objective function value, as follows:(24)alphap,t=soccp,t|argc=1,2,…,Ncminfsoccp,t.The alpha coyote and other ones naturally influence the social behaviors of coyotes. In the COA, the median coyotect p,t j is employed to represent the cultural tendencies of each group of coyotes as follows:(25)ctjp,t=mediansocc,jp,t.The social status of coyotes can be updated according to the alpha coyoteδα and median coyote δt, while the growth of coyotes within the group can be described by(26)new_soccp,t=soccp,t+r1δα+r2δt,δα=alpha−Xrc1,δt=ct−Xrc2.In (26), new_soccp,t represents the coyote social condition after the update; both random numbers r1 and r2 are within the probability uniform distribution [0,1], and Xrc1and Xrc2 denote any two coyotes in the group that are not equal to c.


#### 4.1.3. Birth and Death of Coyotes

The birth and death of coyotes can increase the diversity of a population. Young pups are born according to the following equation:(27)pupjp,t=socrb1,jp,t,rj<Psorj=j1,socrb2,jp,t,rj≥Ps+Pαorj=j2,Rj,otherwise.In (27), rj denotes a random number generated with a uniform probability distribution within [0, 1], j1 and j2 represent two randomly generated dimensions, and Rj denotes a random number generated randomly within the jth dimensional decision variable. The probabilities of the scatter Ps and the association Pα are calculated as follows:(28)Ps=1D,Pα=1−Ps2,where D denotes the dimension. After the coyote pup is born, the newly produced pup is first evaluated for social adaptability and then compared to the group’s worst adaptable and oldest coyote. If the pup adapts better, the oldest coyote dies while the pup is retained, and the age of the pup is set to 0; otherwise, the pup dies.


#### 4.1.4. Migration of Coyotes

In the COA, coyotes migrate between groups with probabilityPe. At first, a coyote is randomly assigned to a group. However, as the coyote grows, the group expels and forces it to migrate. Two random coyotes from different groups swap their positions with a probability Pe, which is calculated as(29)Pe=0.005Nc2.


### 4.2. Improving the Algorithm of the Coyote Optimization

In the COA, the alpha and median coyotes in the group lead to the growth of the whole group. Although the mechanism of coyote birth and death can jump out of the local optimum, the conventional coyote algorithm still suffers from the defects such as quickly falling into the local optimum, low convergence rate, and insufficient exploration ability for a high-dimensionality and complex model. Therefore, this manuscript adopts a novel growth approach.


#### 4.2.1. Initializing Populations Based on Tent Chaotic Sequences

Chaos is a characteristic of nonlinear dynamical systems with bounded unstable dynamical behavior, ergodicity, and nonperiodic behavior [33]. The idea of utilizing chaotic sequences instead of random sequences has been employed in the optimization theory in the literature due to the advantages of chaos like randomness and ergodicity [34–36]. The population is initialized in the standard coyote algorithm utilizing the rand function, while the sample is not distributed uniformly and the distribution of individuals has some extreme values, reducing the probability of finding the global optimum. Therefore, a uniformly distributed chaotic sequence generated by tent mapping is utilized in the manuscript in the initialization stage to enhance the traversal of the population, to improve the exploration ability, to reduce the adverse effect of the unevenly distributed initial population within the search for an optimum, and makes it easier to escape from a local minimum.The tent mapping is given as(30)zi+1=2zi,0≤z≤0.5,21−zi,0.5<z≤1.The Bernoulli shift transformation is applied to the tent mapping as follows:(31)zi+1=2zimod1,where mod denotes a function to find the remainder of the division. Equation (31) is employed to obtain the chaotic variables Zn and apply them to the solution space of the real problem.The equation for the initial population generated by the tent mapping is defined by(32)socc,jp,t=lbj+znubj−lbj.


#### 4.2.2. Leading the Growth of the Best and the Worst Coyotes with the Beetle Antennae Search Strategy

If only the best coyote leads the growth of the coyotes in the group, the algorithm can easily fall into local optimality, while the birth and death mechanisms of a coyote determine the global search ability limitations. According to the idea of the barrel principle, the worst coyote in the group also has a significant impact on the pack, and if the worst coyote is strengthened, the best and the median coyotes inevitably grow more optimally. The best coyote leads the other coyotes in the group to grow, while the guidance of the best coyote can improve the convergence accuracy in the local search process. However, the best coyote falling into local optimum can also affect the other coyotes in the group. Therefore, to further optimize the growth of the best and worst coyotes in the group, the improved coyote algorithm (beetle antennae search coyote optimization algorithm, the BCOA) is proposed to strengthen the growth of the best and the worst coyotes based on the beetle antennae search strategy.The beetle antennae search (BAS) is a heuristic algorithm that simulates the search of a beetle for food [37]. The beetle determines the concentration of food in the left and right directions by its left and right whiskers and moves towards the direction of a higher concentration as follows:(33)xr=x+d0dir,xl=x−d0dir,where dir denotes a random direction vector; Xr and Xl denote a position located in the right and left search regions, respectively; d0 represents the length of the antenna, which should be long enough to cover the appropriate initial search region in the iteration to escape from the local optimum in the initial state and then gradually decay over time.(34)d0t=0.95d0t−1+0.01.By analogizing coyotes to beetles based on comparing the adaptations of left and right and moving towards a better direction, the best and the worst coyotes grow according to(35)xdt=xdt−1+dirδtsignfxr−fxl,where xd denotes the social condition of the best and worst coyotes in the group after completing the update; δt represents the step size of the search; sign denotes the sign function, and f(xl) and f(xr) denote the fitness values of the left and right whiskers, respectively.The convergence speed of the coyote algorithm and its global search ability can be improved by strengthening the growth of coyotes through BAS [38].


#### 4.2.3. Dynamic Adjustment of the Number of Coyote Groups

In the COA, the number of coyote groups and the number of coyotes in each group can significantly influence the performance of the algorithm. Assuming that the total number of coyotes is given, the higher the number of groups, the smaller the number of each group with the smaller growth space, but the search ability for the global optimal solution is enhanced. Conversely, the higher the number of each group, the stronger the local search ability. Therefore, the number of groups is set to 20, with five coyotes per group in the early iterations of the algorithm and 10 with ten coyotes per group in the later iterations.By setting in this way, in the late search period, the number of groups is high, which enhances the positive feedback effect of the global solution and the local search ability; in the early search period, the number of groups is low, which weakens the positive feedback effect of the global solution and enhances the global search ability. Therefore, dynamically adjusting the number of coyotes in a group parameter not only improves the operability but also can better balance the exploration and exploitation ability. In addition, random grouping after dynamically adjusting the parameters eliminates the process of coyote repulsion and admission by the group and improves the operability. The specific distribution is shown in Figure3.Figure 3
The allocation of the number of groups Np and the number of coyotes in the group Nc.The flowchart of the BCOA is shown in Figure4.Figure 4
The flowchart of the BCOA.


#### 4.2.4. Offspring Generation Employing Genetic Crossover Operations

The conventional birth of the coyote is shown in (27). To prevent falling into the local optimum, the genetic crossover strategy is introduced to increase the diversity of the population further and expand the search space, thus enhancing the probability of finding the global optimum.Crossover variation is applied to the random dimension of the two parental coyotes, taking the one with better adaptability as the newborn pup and then the fitness value of a new pup is calculated. If the pup is worse than all the older coyotes, it dies; otherwise, the coyote with the oldest and worst adaptability is replaced.Cross-probability CR setting. The first period fluctuates slightly around a CR of 0.5, with higher diversity and enhanced exploration ability. The later period jumps significantly around a CR of 0.5, producing new solutions dominated by one of the operations, with reduced diversity and enhanced exploitation ability. The calculation is as follows:(36)CR=0.5×sin2π×0.25×t+π×tMaxDT+1,where MaxDT denotes the maximum number of iterations.


### 4.3. Solution Method of the Improved Coyote Algorithm Combined with Quadratic Programming

Although the conventional centralized solution method exposes much information about each agent, such as the objective function and the equipment information, each agent cannot divulge its trade secrets to its competitors in the actual electricity market. Therefore, this paper combines the improved coyote algorithm with quadratic programming to propose a distributed solution method. The steps of the distributed solution algorithm are shown as follows:(1)
Initialize the coyote population with the chaotic mapping and send the electricity and heat prices determined by the ESO to the lower level(2)
Employ quadratic programming based on pricing signals from the ESO to resolve the objective functions of the EP and the EU and send the energy trading scheme back to the ESO(3)
Calculate the ESO profit based on the feedback power of the lower level(4)
Update the prices as equations (26) and (35), then replace the best-adapted price with the objective function of the ESO, perform the selection operation, and take the optimal tariff as the internal electricity price for the next iteration(5)
If the game has reached equilibrium, output the result; otherwise, go to the next iteration


### 4.4. The Validation of the BCOA Algorithm

To verify the performance of the BCOA that is run on four standard test functions with the PSO, DE, and GWO algorithms, the results of the four algorithms are then compared. The expressions and variable ranges of these four functions are shown in Table1.Table 1
The benchmark functions.
FunctionsDimRangefminF1x=∑i=1nxi230[−100, 100]0F2x=∑i=1nxi+∏i=130xi30[−10, −10]0F3x=∑i=1n∑j=1ixj230[−100, 100]0F4x=maxxi,1≤i≤n30[−100, 100]0All algorithms were run 30 times independently on the standard test function, and then the mean and standard deviation of the obtained optimal solutions were recorded. In all experiments, the population size and the maximum number of iterations are set to 40 and 200, respectively.These four functions were resolved by employing the MATLAB software employing the PSO, the DE, and the GWO algorithms and compared with the results produced by the BCOA algorithm as shown in Table2.Table 2
The statistical results of the benchmark functions.
FunctionsBCOAPSODEGWOp valuesF1Ave1.43E −471.08E +031.89E +013.72E −315.6E −12STD4.85E −477.56E +024.83E +018.64E −31F2Ave5.43E −358.46E +015.83E −021.18E −183.2E −10STD7.73E −352.21E +012.21E −019.52E −19F3Ave4.68E −528.84E +038.42E +021.40E −042.8E −11STD1.57E −553.60E +034.59E +022.36E −04F4Ave2.73E −282.08E +013.36E +012.83E −068.4E −6STD4.52E −292.95E +006.98E +003.96E −06Table2 depicts that the BCOA proposed in paper has the best comprehensive performance. In addition, the p-values are all less than 0.05, and the null hypothesis is rejected for all functions. There exists an obvious difference between the four algorithms.


## 4.1. The Coyote Optimization Algorithm

The main difference between the COA and the most intelligent algorithms is that coyotes are divided into several groups, and the internal social influences are considered [25]. The COA should only set some control parameters, including the number of coyotes Np, the number of coyotes in each group Nc, the population of coyotes N, and the maximum number of iterations MaxIter.


### 4.1.1. Initialization of Coyote Populations

The coyote populationN consists of Np groups of coyotes with Nc coyotes in each group, while the coyotes can be initialized in the search space [lb, ub] through the following equation:(22)socc,jp,t=lbj+rjubj−lbj,where c = [1, 2, ...Nc], p = [1, 2, ...Np], and j = [1, 2, ..., D], while D is the dimension of the optimization problem, and rj is a random number generated by a uniform probability distribution within [0, 1]. socp,t c,j is randomly initialized for the jth dimension of the cth coyote in the pth group, and lbj and ubj denote the lower and upper bounds of the jth dimension of the coyote, respectively.


### 4.1.2. Coyote Growth in the Group

The adaptive ability of coyotes can be evaluated according to the objective function defined by(23)fitcp,t=fsoccp,t.Naturally, the alpha coyote is the best socially adapted coyote. In the COA, it corresponds to the best (minimum or maximum) objective function value, as follows:(24)alphap,t=soccp,t|argc=1,2,…,Ncminfsoccp,t.The alpha coyote and other ones naturally influence the social behaviors of coyotes. In the COA, the median coyotect p,t j is employed to represent the cultural tendencies of each group of coyotes as follows:(25)ctjp,t=mediansocc,jp,t.The social status of coyotes can be updated according to the alpha coyoteδα and median coyote δt, while the growth of coyotes within the group can be described by(26)new_soccp,t=soccp,t+r1δα+r2δt,δα=alpha−Xrc1,δt=ct−Xrc2.In (26), new_soccp,t represents the coyote social condition after the update; both random numbers r1 and r2 are within the probability uniform distribution [0,1], and Xrc1and Xrc2 denote any two coyotes in the group that are not equal to c.


### 4.1.3. Birth and Death of Coyotes

The birth and death of coyotes can increase the diversity of a population. Young pups are born according to the following equation:(27)pupjp,t=socrb1,jp,t,rj<Psorj=j1,socrb2,jp,t,rj≥Ps+Pαorj=j2,Rj,otherwise.In (27), rj denotes a random number generated with a uniform probability distribution within [0, 1], j1 and j2 represent two randomly generated dimensions, and Rj denotes a random number generated randomly within the jth dimensional decision variable. The probabilities of the scatter Ps and the association Pα are calculated as follows:(28)Ps=1D,Pα=1−Ps2,where D denotes the dimension. After the coyote pup is born, the newly produced pup is first evaluated for social adaptability and then compared to the group’s worst adaptable and oldest coyote. If the pup adapts better, the oldest coyote dies while the pup is retained, and the age of the pup is set to 0; otherwise, the pup dies.


### 4.1.4. Migration of Coyotes

In the COA, coyotes migrate between groups with probabilityPe. At first, a coyote is randomly assigned to a group. However, as the coyote grows, the group expels and forces it to migrate. Two random coyotes from different groups swap their positions with a probability Pe, which is calculated as(29)Pe=0.005Nc2.


## 4.1.1. Initialization of Coyote Populations

The coyote populationN consists of Np groups of coyotes with Nc coyotes in each group, while the coyotes can be initialized in the search space [lb, ub] through the following equation:(22)socc,jp,t=lbj+rjubj−lbj,where c = [1, 2, ...Nc], p = [1, 2, ...Np], and j = [1, 2, ..., D], while D is the dimension of the optimization problem, and rj is a random number generated by a uniform probability distribution within [0, 1]. socp,t c,j is randomly initialized for the jth dimension of the cth coyote in the pth group, and lbj and ubj denote the lower and upper bounds of the jth dimension of the coyote, respectively.


## 4.1.2. Coyote Growth in the Group

The adaptive ability of coyotes can be evaluated according to the objective function defined by(23)fitcp,t=fsoccp,t.Naturally, the alpha coyote is the best socially adapted coyote. In the COA, it corresponds to the best (minimum or maximum) objective function value, as follows:(24)alphap,t=soccp,t|argc=1,2,…,Ncminfsoccp,t.The alpha coyote and other ones naturally influence the social behaviors of coyotes. In the COA, the median coyotect p,t j is employed to represent the cultural tendencies of each group of coyotes as follows:(25)ctjp,t=mediansocc,jp,t.The social status of coyotes can be updated according to the alpha coyoteδα and median coyote δt, while the growth of coyotes within the group can be described by(26)new_soccp,t=soccp,t+r1δα+r2δt,δα=alpha−Xrc1,δt=ct−Xrc2.In (26), new_soccp,t represents the coyote social condition after the update; both random numbers r1 and r2 are within the probability uniform distribution [0,1], and Xrc1and Xrc2 denote any two coyotes in the group that are not equal to c.


## 4.1.3. Birth and Death of Coyotes

The birth and death of coyotes can increase the diversity of a population. Young pups are born according to the following equation:(27)pupjp,t=socrb1,jp,t,rj<Psorj=j1,socrb2,jp,t,rj≥Ps+Pαorj=j2,Rj,otherwise.In (27), rj denotes a random number generated with a uniform probability distribution within [0, 1], j1 and j2 represent two randomly generated dimensions, and Rj denotes a random number generated randomly within the jth dimensional decision variable. The probabilities of the scatter Ps and the association Pα are calculated as follows:(28)Ps=1D,Pα=1−Ps2,where D denotes the dimension. After the coyote pup is born, the newly produced pup is first evaluated for social adaptability and then compared to the group’s worst adaptable and oldest coyote. If the pup adapts better, the oldest coyote dies while the pup is retained, and the age of the pup is set to 0; otherwise, the pup dies.


## 4.1.4. Migration of Coyotes

In the COA, coyotes migrate between groups with probabilityPe. At first, a coyote is randomly assigned to a group. However, as the coyote grows, the group expels and forces it to migrate. Two random coyotes from different groups swap their positions with a probability Pe, which is calculated as(29)Pe=0.005Nc2.


## 4.2. Improving the Algorithm of the Coyote Optimization

In the COA, the alpha and median coyotes in the group lead to the growth of the whole group. Although the mechanism of coyote birth and death can jump out of the local optimum, the conventional coyote algorithm still suffers from the defects such as quickly falling into the local optimum, low convergence rate, and insufficient exploration ability for a high-dimensionality and complex model. Therefore, this manuscript adopts a novel growth approach.


### 4.2.1. Initializing Populations Based on Tent Chaotic Sequences

Chaos is a characteristic of nonlinear dynamical systems with bounded unstable dynamical behavior, ergodicity, and nonperiodic behavior [33]. The idea of utilizing chaotic sequences instead of random sequences has been employed in the optimization theory in the literature due to the advantages of chaos like randomness and ergodicity [34–36]. The population is initialized in the standard coyote algorithm utilizing the rand function, while the sample is not distributed uniformly and the distribution of individuals has some extreme values, reducing the probability of finding the global optimum. Therefore, a uniformly distributed chaotic sequence generated by tent mapping is utilized in the manuscript in the initialization stage to enhance the traversal of the population, to improve the exploration ability, to reduce the adverse effect of the unevenly distributed initial population within the search for an optimum, and makes it easier to escape from a local minimum.The tent mapping is given as(30)zi+1=2zi,0≤z≤0.5,21−zi,0.5<z≤1.The Bernoulli shift transformation is applied to the tent mapping as follows:(31)zi+1=2zimod1,where mod denotes a function to find the remainder of the division. Equation (31) is employed to obtain the chaotic variables Zn and apply them to the solution space of the real problem.The equation for the initial population generated by the tent mapping is defined by(32)socc,jp,t=lbj+znubj−lbj.


### 4.2.2. Leading the Growth of the Best and the Worst Coyotes with the Beetle Antennae Search Strategy

If only the best coyote leads the growth of the coyotes in the group, the algorithm can easily fall into local optimality, while the birth and death mechanisms of a coyote determine the global search ability limitations. According to the idea of the barrel principle, the worst coyote in the group also has a significant impact on the pack, and if the worst coyote is strengthened, the best and the median coyotes inevitably grow more optimally. The best coyote leads the other coyotes in the group to grow, while the guidance of the best coyote can improve the convergence accuracy in the local search process. However, the best coyote falling into local optimum can also affect the other coyotes in the group. Therefore, to further optimize the growth of the best and worst coyotes in the group, the improved coyote algorithm (beetle antennae search coyote optimization algorithm, the BCOA) is proposed to strengthen the growth of the best and the worst coyotes based on the beetle antennae search strategy.The beetle antennae search (BAS) is a heuristic algorithm that simulates the search of a beetle for food [37]. The beetle determines the concentration of food in the left and right directions by its left and right whiskers and moves towards the direction of a higher concentration as follows:(33)xr=x+d0dir,xl=x−d0dir,where dir denotes a random direction vector; Xr and Xl denote a position located in the right and left search regions, respectively; d0 represents the length of the antenna, which should be long enough to cover the appropriate initial search region in the iteration to escape from the local optimum in the initial state and then gradually decay over time.(34)d0t=0.95d0t−1+0.01.By analogizing coyotes to beetles based on comparing the adaptations of left and right and moving towards a better direction, the best and the worst coyotes grow according to(35)xdt=xdt−1+dirδtsignfxr−fxl,where xd denotes the social condition of the best and worst coyotes in the group after completing the update; δt represents the step size of the search; sign denotes the sign function, and f(xl) and f(xr) denote the fitness values of the left and right whiskers, respectively.The convergence speed of the coyote algorithm and its global search ability can be improved by strengthening the growth of coyotes through BAS [38].


### 4.2.3. Dynamic Adjustment of the Number of Coyote Groups

In the COA, the number of coyote groups and the number of coyotes in each group can significantly influence the performance of the algorithm. Assuming that the total number of coyotes is given, the higher the number of groups, the smaller the number of each group with the smaller growth space, but the search ability for the global optimal solution is enhanced. Conversely, the higher the number of each group, the stronger the local search ability. Therefore, the number of groups is set to 20, with five coyotes per group in the early iterations of the algorithm and 10 with ten coyotes per group in the later iterations.By setting in this way, in the late search period, the number of groups is high, which enhances the positive feedback effect of the global solution and the local search ability; in the early search period, the number of groups is low, which weakens the positive feedback effect of the global solution and enhances the global search ability. Therefore, dynamically adjusting the number of coyotes in a group parameter not only improves the operability but also can better balance the exploration and exploitation ability. In addition, random grouping after dynamically adjusting the parameters eliminates the process of coyote repulsion and admission by the group and improves the operability. The specific distribution is shown in Figure3.Figure 3
The allocation of the number of groups Np and the number of coyotes in the group Nc.The flowchart of the BCOA is shown in Figure4.Figure 4
The flowchart of the BCOA.


### 4.2.4. Offspring Generation Employing Genetic Crossover Operations

The conventional birth of the coyote is shown in (27). To prevent falling into the local optimum, the genetic crossover strategy is introduced to increase the diversity of the population further and expand the search space, thus enhancing the probability of finding the global optimum.Crossover variation is applied to the random dimension of the two parental coyotes, taking the one with better adaptability as the newborn pup and then the fitness value of a new pup is calculated. If the pup is worse than all the older coyotes, it dies; otherwise, the coyote with the oldest and worst adaptability is replaced.Cross-probability CR setting. The first period fluctuates slightly around a CR of 0.5, with higher diversity and enhanced exploration ability. The later period jumps significantly around a CR of 0.5, producing new solutions dominated by one of the operations, with reduced diversity and enhanced exploitation ability. The calculation is as follows:(36)CR=0.5×sin2π×0.25×t+π×tMaxDT+1,where MaxDT denotes the maximum number of iterations.


## 4.2.1. Initializing Populations Based on Tent Chaotic Sequences

Chaos is a characteristic of nonlinear dynamical systems with bounded unstable dynamical behavior, ergodicity, and nonperiodic behavior [33]. The idea of utilizing chaotic sequences instead of random sequences has been employed in the optimization theory in the literature due to the advantages of chaos like randomness and ergodicity [34–36]. The population is initialized in the standard coyote algorithm utilizing the rand function, while the sample is not distributed uniformly and the distribution of individuals has some extreme values, reducing the probability of finding the global optimum. Therefore, a uniformly distributed chaotic sequence generated by tent mapping is utilized in the manuscript in the initialization stage to enhance the traversal of the population, to improve the exploration ability, to reduce the adverse effect of the unevenly distributed initial population within the search for an optimum, and makes it easier to escape from a local minimum.The tent mapping is given as(30)zi+1=2zi,0≤z≤0.5,21−zi,0.5<z≤1.The Bernoulli shift transformation is applied to the tent mapping as follows:(31)zi+1=2zimod1,where mod denotes a function to find the remainder of the division. Equation (31) is employed to obtain the chaotic variables Zn and apply them to the solution space of the real problem.The equation for the initial population generated by the tent mapping is defined by(32)socc,jp,t=lbj+znubj−lbj.


## 4.2.2. Leading the Growth of the Best and the Worst Coyotes with the Beetle Antennae Search Strategy

If only the best coyote leads the growth of the coyotes in the group, the algorithm can easily fall into local optimality, while the birth and death mechanisms of a coyote determine the global search ability limitations. According to the idea of the barrel principle, the worst coyote in the group also has a significant impact on the pack, and if the worst coyote is strengthened, the best and the median coyotes inevitably grow more optimally. The best coyote leads the other coyotes in the group to grow, while the guidance of the best coyote can improve the convergence accuracy in the local search process. However, the best coyote falling into local optimum can also affect the other coyotes in the group. Therefore, to further optimize the growth of the best and worst coyotes in the group, the improved coyote algorithm (beetle antennae search coyote optimization algorithm, the BCOA) is proposed to strengthen the growth of the best and the worst coyotes based on the beetle antennae search strategy.The beetle antennae search (BAS) is a heuristic algorithm that simulates the search of a beetle for food [37]. The beetle determines the concentration of food in the left and right directions by its left and right whiskers and moves towards the direction of a higher concentration as follows:(33)xr=x+d0dir,xl=x−d0dir,where dir denotes a random direction vector; Xr and Xl denote a position located in the right and left search regions, respectively; d0 represents the length of the antenna, which should be long enough to cover the appropriate initial search region in the iteration to escape from the local optimum in the initial state and then gradually decay over time.(34)d0t=0.95d0t−1+0.01.By analogizing coyotes to beetles based on comparing the adaptations of left and right and moving towards a better direction, the best and the worst coyotes grow according to(35)xdt=xdt−1+dirδtsignfxr−fxl,where xd denotes the social condition of the best and worst coyotes in the group after completing the update; δt represents the step size of the search; sign denotes the sign function, and f(xl) and f(xr) denote the fitness values of the left and right whiskers, respectively.The convergence speed of the coyote algorithm and its global search ability can be improved by strengthening the growth of coyotes through BAS [38].


## 4.2.3. Dynamic Adjustment of the Number of Coyote Groups

In the COA, the number of coyote groups and the number of coyotes in each group can significantly influence the performance of the algorithm. Assuming that the total number of coyotes is given, the higher the number of groups, the smaller the number of each group with the smaller growth space, but the search ability for the global optimal solution is enhanced. Conversely, the higher the number of each group, the stronger the local search ability. Therefore, the number of groups is set to 20, with five coyotes per group in the early iterations of the algorithm and 10 with ten coyotes per group in the later iterations.By setting in this way, in the late search period, the number of groups is high, which enhances the positive feedback effect of the global solution and the local search ability; in the early search period, the number of groups is low, which weakens the positive feedback effect of the global solution and enhances the global search ability. Therefore, dynamically adjusting the number of coyotes in a group parameter not only improves the operability but also can better balance the exploration and exploitation ability. In addition, random grouping after dynamically adjusting the parameters eliminates the process of coyote repulsion and admission by the group and improves the operability. The specific distribution is shown in Figure3.Figure 3
The allocation of the number of groups Np and the number of coyotes in the group Nc.The flowchart of the BCOA is shown in Figure4.Figure 4
The flowchart of the BCOA.


## 4.2.4. Offspring Generation Employing Genetic Crossover Operations

The conventional birth of the coyote is shown in (27). To prevent falling into the local optimum, the genetic crossover strategy is introduced to increase the diversity of the population further and expand the search space, thus enhancing the probability of finding the global optimum.Crossover variation is applied to the random dimension of the two parental coyotes, taking the one with better adaptability as the newborn pup and then the fitness value of a new pup is calculated. If the pup is worse than all the older coyotes, it dies; otherwise, the coyote with the oldest and worst adaptability is replaced.Cross-probability CR setting. The first period fluctuates slightly around a CR of 0.5, with higher diversity and enhanced exploration ability. The later period jumps significantly around a CR of 0.5, producing new solutions dominated by one of the operations, with reduced diversity and enhanced exploitation ability. The calculation is as follows:(36)CR=0.5×sin2π×0.25×t+π×tMaxDT+1,where MaxDT denotes the maximum number of iterations.


## 4.3. Solution Method of the Improved Coyote Algorithm Combined with Quadratic Programming

Although the conventional centralized solution method exposes much information about each agent, such as the objective function and the equipment information, each agent cannot divulge its trade secrets to its competitors in the actual electricity market. Therefore, this paper combines the improved coyote algorithm with quadratic programming to propose a distributed solution method. The steps of the distributed solution algorithm are shown as follows:(1)
Initialize the coyote population with the chaotic mapping and send the electricity and heat prices determined by the ESO to the lower level(2)
Employ quadratic programming based on pricing signals from the ESO to resolve the objective functions of the EP and the EU and send the energy trading scheme back to the ESO(3)
Calculate the ESO profit based on the feedback power of the lower level(4)
Update the prices as equations (26) and (35), then replace the best-adapted price with the objective function of the ESO, perform the selection operation, and take the optimal tariff as the internal electricity price for the next iteration(5)
If the game has reached equilibrium, output the result; otherwise, go to the next iteration


## 4.4. The Validation of the BCOA Algorithm

To verify the performance of the BCOA that is run on four standard test functions with the PSO, DE, and GWO algorithms, the results of the four algorithms are then compared. The expressions and variable ranges of these four functions are shown in Table1.Table 1
The benchmark functions.
FunctionsDimRangefminF1x=∑i=1nxi230[−100, 100]0F2x=∑i=1nxi+∏i=130xi30[−10, −10]0F3x=∑i=1n∑j=1ixj230[−100, 100]0F4x=maxxi,1≤i≤n30[−100, 100]0All algorithms were run 30 times independently on the standard test function, and then the mean and standard deviation of the obtained optimal solutions were recorded. In all experiments, the population size and the maximum number of iterations are set to 40 and 200, respectively.These four functions were resolved by employing the MATLAB software employing the PSO, the DE, and the GWO algorithms and compared with the results produced by the BCOA algorithm as shown in Table2.Table 2
The statistical results of the benchmark functions.
FunctionsBCOAPSODEGWOp valuesF1Ave1.43E −471.08E +031.89E +013.72E −315.6E −12STD4.85E −477.56E +024.83E +018.64E −31F2Ave5.43E −358.46E +015.83E −021.18E −183.2E −10STD7.73E −352.21E +012.21E −019.52E −19F3Ave4.68E −528.84E +038.42E +021.40E −042.8E −11STD1.57E −553.60E +034.59E +022.36E −04F4Ave2.73E −282.08E +013.36E +012.83E −068.4E −6STD4.52E −292.95E +006.98E +003.96E −06Table2 depicts that the BCOA proposed in paper has the best comprehensive performance. In addition, the p-values are all less than 0.05, and the null hypothesis is rejected for all functions. There exists an obvious difference between the four algorithms.


## 5. Case Study

### 5.1. Basic Data

A PIES in Northern China is employed as a practical example to evaluate the optimal dispatching strategy of the mentioned PIES model. The general overview of this PIES is shown in Figure5. The trading patterns are similar in winter and summer, while only the output of the new energy is somewhat different. This article only analyses the energy transaction process and the operating state of the system in winter, considering its length limitation. Table 3 shows the time-of-use prices of the grid and the gas price. The feed-in tariff of the grid is 0.35 ¥/kWh, and the heating price ranges between 0.15 and 0.5 ¥/kWh. Consider that the customer preference coefficients for electrical and heat energies are as presented in Ref. [13]. Figure 6 shows the prediction curves of renewable energy and load power for a typical winter day based on the probability scenarios method, while the device parameters of each system are presented in Table 4.Figure 5
The general overview of the case.Table 3
Electricity and heat prices.
—Times (h)Price(¥/kWh)Grid price1 : 00–9:000.49 : 00–12 : 00, 16 : 00–20 : 00, 23 : 00–1:000.812 : 00–16 : 00, 20 : 00–23 : 001.25Gas priceThe whole day0.35Figure 6
The curve of the power prediction.Table 4
The operation parameters for each type of equipment.
Equipment parametersNumerical valuesMaximum charging and discharging powers of the battery (kW)200, 200Maximum charging and discharging powers of the TST (kW)200, 150Heat storage and release efficiencies of the TST0.98, 0.98Charging and discharging efficiencies of the battery0.95, 0.95Energy self-loss rate of the TST0.02The capacity of the GB (kW)700The capacity of the GT (kW)400Upper and lower ramp rate limits of the GT (kW)220, −220Upper and lower ramp rate limits of the GB (kW)400, −220The efficiency of the GT0.35The efficiency of the GB0.9The efficiency of the WHB0.83


### 5.2. Algorithm Comparison

To verify the BCOA optimization for the constructed model, the improved BCOA algorithm is compared with the conventional COA and the improved differential evolutionary algorithm (ADE) to optimize the revenue objective function of the ESO. Figure7 depicts the comparison results.Figure 7
The comparison of convergence curves of different algorithms.The population of coyotes is adjusted as 100, while the number of coyote groups is chosen as 20 with five coyotes in each group, and the maximum number of iterations is selected as 100. Table5 summarizes that the improved BCOA algorithm converges at 26 iterations, while the conventional COA and the ADE algorithms converge at 39 and 55 iterations, respectively, demonstrating the superiority of the BCOA over the other algorithms concerning the convergence speed. Given the optimal value, the BCOA finds the revenue of the ESO as 10781 ¥, the revenue of the COA as 9503 ¥, and the revenue of the ADE as 8803 ¥. According to the obtained results, the BCOA finally finds the highest optimal daily revenue, and the optimal solution is superior to the other algorithms.Table 5
The optimization results of different algorithms.
AlgorithmsESO revenue (¥)Number of iterationsADE880355COA950339BCOA1078126


### 5.3. Analysis of the Game Results

Figure8 depicts the comparisons of the convergence processes of the ESO, the EP, and the EU benefits. The benefits of each agent converge at the 26th iteration, demonstrating the fast convergence rate of the proposed solution method and the existence of the game equilibrium. As the iteration number increases, the game agents have different convergence trends, verifying the existence of a game process among agents. Moreover, the ESO benefits show an overall upward trend, while the lower-level benefits tend to decline, reflecting the leadership position of the ESO in the Stackelberg game. After the 26th iteration, the game equilibrium reaches between the leader and the followers, and each agent cannot adjust its strategy individually to obtain a higher profit. The ESO and the EP benefits are 10781, and 5511 ¥, respectively, while the consumer surplus of the EU is 12621 ¥.Figure 8
The convergence curves of income of each player.Figure9 shows the price of the ESO for energy transactions after the convergence of the iterative process. The green and blue lines in Figure 9(a) indicate the upper and lower bounds of the heat prices of the park within which the ESO should set more competitive heat prices for the followers. In Figure 9(b), the electricity price set by the ESO should be between the time-of-use electricity prices of the grid and the feed-in tariff to satisfy the constraint. The peaks of electricity selling prices are at 12 : 00 and 21 : 00 because these are the peak hours for the electricity consumption of the customers and the peaks of the PV and the WT outputs. On the other hand, higher selling prices can promote renewable energy accommodation.Figure 9
Energy prices of the ESO. (a) The curve of heat price. (b) The curve of electricity price.
(a)(b)Figure10 shows the operation of the energy storage device. The electric and thermal loads before and after the demand response are compared in Figure 11. The EU reduces its energy costs and improves the overall benefit by adjusting the use time of electric cars and washing machines and other equipment. After adopting the demand response strategy, the electric load fluctuation is smoothed out, which plays the role of peak load shifting, reducing the energy consumption burden of the system and many hidden problems of the grid and the PIES. The heat load has been cut-in all hours, while the cut is less in the hours with lower heat load to ensure the comfort of the EU.Figure 10
The details of the energy storage output.Figure 11
The electric and heat load curves before and after optimization.Figures12 and 13 show the electric and thermal energy outputs of each device in the CCHP system at the equilibrium state. In Figure 12, all the electricity from the WT and the PV is sold to the ESO, sending the surplus energy to energy storage devices to improve renewable energy accommodation. At night, the EU takes advantage of the lower price of electricity to charge some devices like electric vehicles, while the ESP also buys electricity at a lower price. During the peak period of electrical load, the supply of the EP cannot meet the electrical load, while the ESP and the grid supplements the shortfall. The thermal load is at its peak between 9 : 00–11 : 00 and 22 : 00 and requires a portion of the thermal energy supplied by storage devices in addition to the thermal energy supplied by the GB and the GT.Figure 12
Electricity output optimal results of the PIES.Figure 13
The heat output optimal results of the PIES.To verify the rationality and effectiveness of the proposed optimization strategy, two different scenarios are selected, and the returns of each agent are calculated for comparison.(1)
Scenario 1: the proposed optimization strategy is utilized to construct the PIES model based on the Stackelberg game.(2)
Scenario 2: no energy storage devices are considered. The leader is the ESO, and the followers are only producers and energy users.The above scenarios are adopted to calculate the corresponding payoffs of the game participants, as shown in Table6.Table 6
Income comparison of each player under different scenarios.
ScenariosBenefits (¥)ESOEPEU11078159641262129996337011391A comparison between Scenarios 1 and 2 reveals that the revenue of the ESO, the EP, and the EU can be improved after introducing the energy storage, with the revenue of the ESO increasing by 785 ¥. Energy storage devices also help the consumption of renewable energy sources, reducing the rate of wind and light curtailment and improving the revenue of the overall system.To verify the validity of the model proposed in this paper regarding the load side, the following comparison of the benefits of each subject under two scenarios is considered separately:(1)
The optimization of the load side is not considered, and the data obtained from the forecast is used for the load-side power(2)
Considering the level ability of the electric load and the curtailing of the thermal load, the load-side power utilizes the optimized data in Figure11From the comparison results of the two scenarios in Table7, when the load-side adjustability is considered, the energy cost decreases from 22875 to 20254, and the utility function increases from 32106 to 32875. Overall, the load-side objective function increases from 9231 to 12621. It shows that the model proposed in this paper can reduce energy costs and improve the economy of energy use while ensuring the comfort of energy use. The revenue of the producer decreases because the load side cuts some of the load.Table 7
The comparison of income of each player under different scenarios.
ScenariosBenefits (¥)ESOEPEU1968168399231210781596412621


## 5.1. Basic Data

A PIES in Northern China is employed as a practical example to evaluate the optimal dispatching strategy of the mentioned PIES model. The general overview of this PIES is shown in Figure5. The trading patterns are similar in winter and summer, while only the output of the new energy is somewhat different. This article only analyses the energy transaction process and the operating state of the system in winter, considering its length limitation. Table 3 shows the time-of-use prices of the grid and the gas price. The feed-in tariff of the grid is 0.35 ¥/kWh, and the heating price ranges between 0.15 and 0.5 ¥/kWh. Consider that the customer preference coefficients for electrical and heat energies are as presented in Ref. [13]. Figure 6 shows the prediction curves of renewable energy and load power for a typical winter day based on the probability scenarios method, while the device parameters of each system are presented in Table 4.Figure 5
The general overview of the case.Table 3
Electricity and heat prices.
—Times (h)Price(¥/kWh)Grid price1 : 00–9:000.49 : 00–12 : 00, 16 : 00–20 : 00, 23 : 00–1:000.812 : 00–16 : 00, 20 : 00–23 : 001.25Gas priceThe whole day0.35Figure 6
The curve of the power prediction.Table 4
The operation parameters for each type of equipment.
Equipment parametersNumerical valuesMaximum charging and discharging powers of the battery (kW)200, 200Maximum charging and discharging powers of the TST (kW)200, 150Heat storage and release efficiencies of the TST0.98, 0.98Charging and discharging efficiencies of the battery0.95, 0.95Energy self-loss rate of the TST0.02The capacity of the GB (kW)700The capacity of the GT (kW)400Upper and lower ramp rate limits of the GT (kW)220, −220Upper and lower ramp rate limits of the GB (kW)400, −220The efficiency of the GT0.35The efficiency of the GB0.9The efficiency of the WHB0.83


## 5.2. Algorithm Comparison

To verify the BCOA optimization for the constructed model, the improved BCOA algorithm is compared with the conventional COA and the improved differential evolutionary algorithm (ADE) to optimize the revenue objective function of the ESO. Figure7 depicts the comparison results.Figure 7
The comparison of convergence curves of different algorithms.The population of coyotes is adjusted as 100, while the number of coyote groups is chosen as 20 with five coyotes in each group, and the maximum number of iterations is selected as 100. Table5 summarizes that the improved BCOA algorithm converges at 26 iterations, while the conventional COA and the ADE algorithms converge at 39 and 55 iterations, respectively, demonstrating the superiority of the BCOA over the other algorithms concerning the convergence speed. Given the optimal value, the BCOA finds the revenue of the ESO as 10781 ¥, the revenue of the COA as 9503 ¥, and the revenue of the ADE as 8803 ¥. According to the obtained results, the BCOA finally finds the highest optimal daily revenue, and the optimal solution is superior to the other algorithms.Table 5
The optimization results of different algorithms.
AlgorithmsESO revenue (¥)Number of iterationsADE880355COA950339BCOA1078126


## 5.3. Analysis of the Game Results

Figure8 depicts the comparisons of the convergence processes of the ESO, the EP, and the EU benefits. The benefits of each agent converge at the 26th iteration, demonstrating the fast convergence rate of the proposed solution method and the existence of the game equilibrium. As the iteration number increases, the game agents have different convergence trends, verifying the existence of a game process among agents. Moreover, the ESO benefits show an overall upward trend, while the lower-level benefits tend to decline, reflecting the leadership position of the ESO in the Stackelberg game. After the 26th iteration, the game equilibrium reaches between the leader and the followers, and each agent cannot adjust its strategy individually to obtain a higher profit. The ESO and the EP benefits are 10781, and 5511 ¥, respectively, while the consumer surplus of the EU is 12621 ¥.Figure 8
The convergence curves of income of each player.Figure9 shows the price of the ESO for energy transactions after the convergence of the iterative process. The green and blue lines in Figure 9(a) indicate the upper and lower bounds of the heat prices of the park within which the ESO should set more competitive heat prices for the followers. In Figure 9(b), the electricity price set by the ESO should be between the time-of-use electricity prices of the grid and the feed-in tariff to satisfy the constraint. The peaks of electricity selling prices are at 12 : 00 and 21 : 00 because these are the peak hours for the electricity consumption of the customers and the peaks of the PV and the WT outputs. On the other hand, higher selling prices can promote renewable energy accommodation.Figure 9
Energy prices of the ESO. (a) The curve of heat price. (b) The curve of electricity price.
(a)(b)Figure10 shows the operation of the energy storage device. The electric and thermal loads before and after the demand response are compared in Figure 11. The EU reduces its energy costs and improves the overall benefit by adjusting the use time of electric cars and washing machines and other equipment. After adopting the demand response strategy, the electric load fluctuation is smoothed out, which plays the role of peak load shifting, reducing the energy consumption burden of the system and many hidden problems of the grid and the PIES. The heat load has been cut-in all hours, while the cut is less in the hours with lower heat load to ensure the comfort of the EU.Figure 10
The details of the energy storage output.Figure 11
The electric and heat load curves before and after optimization.Figures12 and 13 show the electric and thermal energy outputs of each device in the CCHP system at the equilibrium state. In Figure 12, all the electricity from the WT and the PV is sold to the ESO, sending the surplus energy to energy storage devices to improve renewable energy accommodation. At night, the EU takes advantage of the lower price of electricity to charge some devices like electric vehicles, while the ESP also buys electricity at a lower price. During the peak period of electrical load, the supply of the EP cannot meet the electrical load, while the ESP and the grid supplements the shortfall. The thermal load is at its peak between 9 : 00–11 : 00 and 22 : 00 and requires a portion of the thermal energy supplied by storage devices in addition to the thermal energy supplied by the GB and the GT.Figure 12
Electricity output optimal results of the PIES.Figure 13
The heat output optimal results of the PIES.To verify the rationality and effectiveness of the proposed optimization strategy, two different scenarios are selected, and the returns of each agent are calculated for comparison.(1)
Scenario 1: the proposed optimization strategy is utilized to construct the PIES model based on the Stackelberg game.(2)
Scenario 2: no energy storage devices are considered. The leader is the ESO, and the followers are only producers and energy users.The above scenarios are adopted to calculate the corresponding payoffs of the game participants, as shown in Table6.Table 6
Income comparison of each player under different scenarios.
ScenariosBenefits (¥)ESOEPEU11078159641262129996337011391A comparison between Scenarios 1 and 2 reveals that the revenue of the ESO, the EP, and the EU can be improved after introducing the energy storage, with the revenue of the ESO increasing by 785 ¥. Energy storage devices also help the consumption of renewable energy sources, reducing the rate of wind and light curtailment and improving the revenue of the overall system.To verify the validity of the model proposed in this paper regarding the load side, the following comparison of the benefits of each subject under two scenarios is considered separately:(1)
The optimization of the load side is not considered, and the data obtained from the forecast is used for the load-side power(2)
Considering the level ability of the electric load and the curtailing of the thermal load, the load-side power utilizes the optimized data in Figure11From the comparison results of the two scenarios in Table7, when the load-side adjustability is considered, the energy cost decreases from 22875 to 20254, and the utility function increases from 32106 to 32875. Overall, the load-side objective function increases from 9231 to 12621. It shows that the model proposed in this paper can reduce energy costs and improve the economy of energy use while ensuring the comfort of energy use. The revenue of the producer decreases because the load side cuts some of the load.Table 7
The comparison of income of each player under different scenarios.
ScenariosBenefits (¥)ESOEPEU1968168399231210781596412621


## 6. Conclusions

The manuscript suggests an optimal dispatch model based on the Stackelberg game, with the ESO as the upper-level leader and the EP and the EU as the lower-level followers, considering the privacy, economy, and stability of the agents. Each agent pursues the highest return under stable operation and formulates its respective trading strategies to reach Stackelberg equilibrium after several games.An implementation method is proposed to determine the game equilibrium, which can protect the privacy of each agent. To resolve the problems of uneven initial population distribution and quickly fall into local optimum in the conventional COA, some improvements are introduced, such as chaotic mapping of initial populations, utilizing the beetle antennae search strategy to strengthen the best and the worst coyotes, dynamic grouping, and genetic crossover. This improves the convergence speed, the global search ability, and the solution stability of the algorithm. These improvements increase the daily revenue of the ESO by 11.8%, while the BCOA-based optimized system increases the revenue while ensuring the stable operation of the system.The algorithm analysis shows that the proposed game model enhances the revenue of each agent and reduces the pressure of energy consumption at the maximum load by introducing the energy storage provider. The transferable electric load is shifted to the valley through the demand response of the consumers, improving the consumer surplus and reducing the load fluctuation.The optimal operation method of the proposed PIES is employed to obtain the optimal equilibrium strategies for each agent in the game process, which can be considered as a reference value for market decisions. To enhance system integrity, future research will investigate the impact of the inclusion of other agents into the game model and the inclusion of more energy conversion and storage devices.



---

*Source: 6204918-2022-10-31.xml*