#### 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.

#### 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. Section 2 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 Figure 1.

##### 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 Figure 2, 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.

##### 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 bywhere represents the revenue of the operator from the sold energy; represents the operating cost of the operator; , , and 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; and represent the prices which the ESO pays for electricity and heat; and denote the prices obtained by selling electricity and heat by the ESO, respectively; and represent the time-of-use electricity price and feed-in tariffs, respectively; and denote the electrical and heat power sold by the ESO to the ESP; and represent the electrical and heat power sold by the ESO to the EU, respectively; and denote the electrical and heat power purchased by ESO from the EP, respectively; and represent the electrical and heat power purchased by the ESO from the the ESP, respectively; 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 bywhere *p*_{hmax} and *p*_{hmin} 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 sales *I*_{EP}, the generation cost *C*_{gen}, and the emission cost *C*_{emi} as follows:where *C*_{gen} represents the generation cost of the CCHP system; denotes the price of the natural gas; represents the volume of the natural gas purchased; *J*_{i} denotes the maintenance cost factor of the *i*^{th} CCHP device; *P*_{i} denotes the output power of the *i*^{th} device; *I* denote the total number of devices in the CCHP; *P*_{emi} represents the gas emission factor; *F*_{gen} denotes the heat energy of the consumed gas; and 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:where *P*_{pv} and *P*_{wt} represent the power generated by photovoltaics (PV) and wind turbine (WT), respectively; *P*_{MT} represents the power generated by the gas turbine; *P*_{GB} and *P*_{RE} denote the thermal output powers of the gas boiler and the waste heat boiler, respectively.

The power relationship for the CCHP system is described bywhere *η*_{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 bywhere *P*_{net} denotes the amount of power trading between the CHP system and the grid, while positive denotes electricity sales and negative represents electricity purchase; *P*_{ey} represents the power of the energy storage, while positive denotes discharging, negative represents charging; and *l*_{fe} and *l*_{te} 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 time *t*:where *P*_{MT, max} and *P*_{GB,max} denote the rated capacity of the gas turbine and the gas boiler, respectively; *P*_{MT,up} and *P*_{GB,up} represent the upper limits of the ramp rate of the gas turbine and the boiler; and *P*_{MT,down} and *P*_{GB,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 aswhere *H*_{tst} (*t*) represents the real-time quantity of heat storage; *η h ch* and *η h dis* denote the charging and discharging efficiencies of the thermal energy; and 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 bywhere *P*_{PV,t} represents the output power of the PV power system at time *t*; *P*_{STC} denotes the maximum output power of solar panel under the standard test; *G*_{STC} represents light intensity under the standard test, which is 1000 W/m^{2}; *k* denotes the temperature coefficient; *T(t)* represents the actual temperature of the solar panel at time *t*; *T*_{STC} denotes the solar panel temperature under the standard test, which is 25°C; *T*_{air} represents the outdoor temperature at time *t*; *V*_{W} denotes the wind speed; *T*^{max} denotes the maximum daily temperature; *T*^{min} represents the minimum daily temperature; and *t*_{p} represents the average daily temperature.

Light intensity usually follows a beta distribution. The probability density function is defined bywhere *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 *G*^{max} 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 time *t* can be expressed as a function of the intermittent wind speed aswhere *ν*_{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 *P*_{rated} represents the rated power of the wind turbine.

The natural incoming wind *ν* usually follows the Weibull probability distribution. The probability density function is defined bywhere *φ* 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

The forecasting equation is defined bywhere *y*_{t} represents the observed value of the time series at time *t*; *a*_{t} denotes the stable component of the load data at time *t*; *b*_{t} represents the tendency trend component of the load data at time *t*; *S*_{t} 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

*U*_{EU} 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:where , *u*_{e}, , and *u*_{h} represent the preference constants for the widely used quadratic utility function [28–30]. *d*_{e} and *d*_{h} denote the electrical and thermal load demands of the EU.

denotes the cost of energy purchase by the user as follows:

The demand response of the EU is divided into shifted electrical load *l*_{te} and curtailable thermal load *l*_{rh}, while the fixed load of the EU is divided into fixed electrical load *l*_{fe} and fixed thermal load *l*_{fh}, which is formulated as

*l *_{te} and *l*_{rh} should satisfy the following constraints as follows:where *l*_{temax} and *l*_{rhmax} 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

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 operator *E*, energy user U, and energy producer (2)The strategies set is defined as *L* = {*L*_{e}, *L*_{u}, *L*_{p}}, where the strategy of the ESO, *L*_{e}, represents the offer of energy, expressed as *L*_{e}=(*p* buy *e*, *p* buy *h*, *p* sell *e*, *p* sell *h*); the strategy of the EP, *L*_{p}, denotes the output of gas turbines and boilers at each period, denoted as *L*_{p} = (*P*_{MT}, *P*_{GB}); the strategy of the user, *L*_{u}, represents the amount of demand response at each moment, denoted as *L*_{u} _{=} (*l*_{te}, *l*_{rh})(3)The payoff set *F* = {*F*_{ESO}, *F*_{EP}, *F*_{EU}} 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 *L*_{e} be the equilibrium strategy of the ESO and *L*_{u} and *L*_{p} be the optimal response strategies of the EU and EP, respectively. The set of strategies (*L*_{e}, *L*_{u}, and *L*_{p}) is the equilibrium solution of the Stackelberg game under the following constraints as follows:where *L**e(*−*i)* stands for the other strategies except *L*_{ei}. 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 sets *L*_{e}, *L*_{u}, and *L*_{p} are nonempty, bounded, and convex sets in the Euclidean space(2)*F*_{EP} and *F*_{EU} represent quasi-concave functions concerning *L*_{p} and *L*_{u,} respectively(3)*F*_{ESO} denotes a continuous function of *L*_{e}, *L*_{u}, *L*_{s}, and *L*_{p}(4)*F*_{EP} and *F*_{EU} denote continuous functions of *L*_{p}, and *L*_{u}, respectively

According 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 *l*_{te} and *l*_{rh}. Thus, *F*_{EU} represents a convex function of *L*_{u}. Similarly, the objective functions of EP are convex. Moreover, since the convex function must be a quasi-concave function, condition (2) is also satisfied.

*F *_{ESO}, *F*_{EP}, and *F*_{EU} 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 *N*_{p}, the number of coyotes in each group *N*_{c}, the population of coyotes *N*, and the maximum number of iterations *MaxIter*.

###### 4.1.1. Initialization of Coyote Populations

The coyote population *N* consists of *N*_{p} groups of coyotes with *N*_{c} coyotes in each group, while the coyotes can be initialized in the search space [lb, ub] through the following equation:where *c* = [1, 2, ...*N*_{c}], *p* = [1, 2, ...*N*_{p}], and *j* = [1, 2, ..., *D*], while *D* is the dimension of the optimization problem, and *r*_{j} is a random number generated by a uniform probability distribution within [0, 1]. *socp,t c,j* is randomly initialized for the *j*^{th} dimension of the *c*^{th} coyote in the *p*^{th} group, and *lb*_{j} and *ub*_{j} denote the lower and upper bounds of the *j*^{th} 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

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:

The alpha coyote and other ones naturally influence the social behaviors of coyotes. In the COA, the median coyote *ct p,t j* is employed to represent the cultural tendencies of each group of coyotes as follows:

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

In (26), new_soc_{c}^{p,t} represents the coyote social condition after the update; both random numbers *r*_{1} and *r*_{2} are within the probability uniform distribution [0,1], and *X*_{rc1}*and X*_{rc2} 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:

In (27), *r*_{j} denotes a random number generated with a uniform probability distribution within [0, 1], *j*_{1} and *j*_{2} represent two randomly generated dimensions, and *R*_{j} denotes a random number generated randomly within the *j*^{th} dimensional decision variable. The probabilities of the scatter *P*_{s} and the association *P*_{α} are calculated as follows: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 probability *P*_{e}. 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 *P*_{e}, which is calculated as

##### 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

The Bernoulli shift transformation is applied to the tent mapping as follows:where mod denotes a function to find the remainder of the division. Equation (31) is employed to obtain the chaotic variables *Z*_{n} 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

###### 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:where **dir** denotes a random direction vector; *X*_{r} and *X*_{l} denote a position located in the right and left search regions, respectively; *d*_{0} 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.

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 towhere *x*_{d} 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*(*x*_{l}) and *f*(*x*_{r}) 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 Figure 3.

The flowchart of the BCOA is shown in Figure 4.

###### 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:where Max*DT* 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 Table 1.

All 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 Table 2.

Table 2 depicts that the BCOA proposed in paper has the best comprehensive performance. In addition, the -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 Figure 5. 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.

##### 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. Figure 7 depicts the comparison results.

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. Table 5 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.

##### 5.3. Analysis of the Game Results

Figure 8 depicts the comparisons of the convergence processes of the ESO, the EP, and the EU benefits. The benefits of each agent converge at the 26^{th} 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 26^{th} 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 9 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.

**(a)**

**(b)**

Figure 10 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.

Figures 12 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.

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 Table 6.

A 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 Figure 11

From the comparison results of the two scenarios in Table 7, 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.

#### 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.

#### Abbreviations

PIES: | Park-level integrated energy system |

CCHP: | Combined cooling heating and power |

DR: | Demand response |

TOU: | Time of use |

COA: | Coyote optimization algorithm |

BCOA: | Beetle antennae search coyote optimization algorithm |

BAS: | Beetle antennae search |

EU: | Energy users |

ESP: | Energy storage provider |

EP: | Energy producer |

ESO: | Energy system operator |

PV: | Photovoltaics |

GB: | Gas boiler |

WHB: | Waste heat boiler |

GT: | Gas turbine |

TST: | Thermal storage tank |

WT: | Wind turbine |

: | Price |

: | Power |

C: | Cost |

I: | Income |

ƞ: | Efficiency |

l: | Load |

F: | Objective function |

L: | Strategy set. |

#### Data Availability

The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.