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

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 benefts 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 benefts 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 efectiveness of the model and the improved algorithm.


Introduction
Te 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][4][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. Te 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.
Te 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][9][10][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 CCHPtype 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. Tere exist multiple agents in the integrated energy system with their demands. Tere may also be conficts between the diferent agents. Te 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. Tey 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 modifed particle swarm optimization algorithm with negative knowledge to resolve the mixedmodel 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 frefy 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 verifed 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 sufers 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 CCHPtype PIES. At frst, 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. Te standard COA and the improved one were utilized to obtain the proft of the energy system operator, thus enhancing the revenue of the operator.
Te 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. Te implementation method is introduced in Section 4. Case analysis and conclusions are presented in Section 5 and Section 6, respectively.

The PIES Multiagent Model
Te 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. Te schematic diagram of the Stackelberg game of the PIES is shown in Figure 1.

Te Description of the System.
Te 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 ofered 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 ofer of the ESO.
Te ESO connects the source and the load and acts as a bridge for energy interaction between the EU and the EP. Te 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. Te ESO has a more fexible pricing mechanism than that of the conventional grid, encouraging consumers to participate in demand response. Te ESO sufers 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.
Te PIES structure is shown in Figure 2, where the natural gas fowing 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. Te 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. Te thermal storage tank (TST) and batteries store the excess energy.
where I sell ESO represents the revenue of the operator from the sold energy; C cost ESO represents the operating cost of the operator; C buy ESO , C grid ESO , and C penal ESO 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; p buy e and p buy h represent the prices which the ESO pays for electricity and heat; p sell e and p sell h denote the prices obtained by selling electricity and heat by the ESO, respectively; p buy g and p sell g represent the time-of-use electricity price and feed-in tarifs, respectively; P es and P hs denote the electrical and heat power sold by the ESO to the ESP; P eu and P hu represent the electrical and heat power sold by the ESO to the EU, respectively; P ep and P hp denote the electrical and heat power purchased by ESO from the EP, respectively; P ey and P hy represent the electrical and heat power purchased by the ESO from the the ESP, respectively; P h 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 ofers of the ESO to ensure that the buying and selling prices are within the market price range for the electrical and thermal energies. Tese constraints are given by where p hmax and p hmin indicate the upper and lower limits of the heating price, respectively.

Model of the Energy Producer.
Te EP maximizes its revenue by adjusting the output of each device in the CCHP system after the ESO gives the energy ofer, and its proft is the diference 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; p gas denotes the price of the natural gas; V g 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; v C and v N denote the carbon dioxide and nitrogen oxide volumes of emissions. Te 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.
Te power relationship for the CCHP system is described by where η MT denotes the power generation efciency of the gas turbine and η RE represents the heat production efciency of the waste heat boiler. Te balance of the electrical power is expressed by where 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 fxed electric load and shifted electrical load, respectively.
In the actual operation, the following constraints should be satisfed to fulfll the constraints on the power output and creep rate of the gas turbine and the gas boiler at time t: 4 International Transactions on Electrical Energy Systems 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. Te energy storage device should satisfy some constraints, such as the thermal storage tank constraint, described as where H tst (t) represents the real-time quantity of heat storage; η h ch and η h dis denote the charging and discharging efciencies of the thermal energy; H min tst and H max tst represent the minimum and maximum heat storage capacity, respectively; and η h denotes the energy self-loss rate of the thermal storage tank. Te energy storage equipment should ensure the consistency of the reserves at the beginning and the end of the daily cycle. Te 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 frst, probability scenarios are generated. Te 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.
Te expression for the output of the photovoltaic power system is represented by where 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 coefcient; 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. Te probability density function is defned by 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 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. Te output of the wind at time t can be expressed as a function of the intermittent wind speed as 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 P rated represents the rated power of the wind turbine. Te natural incoming wind ] usually follows the Weibull probability distribution. Te probability density function is defned by where φ denotes the shape parameter and θ represents the scale parameter.

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. Te method integrates the time series characteristics of a linear trend, seasonal variation, and stochastic fuctuation and combines with the exponential smoothing method to have better forecasting capability. Te basic idea is to decompose the time series with linear trend, seasonal variation, and stochastic fuctuation and combine them with the exponential smoothing method to International Transactions on Electrical Energy Systems estimate the long-term trend, increment of trend, and seasonal fuctuation respectively. A forecasting model is then built and extrapolated to the forecast values.
Te multiplicative Holt-Winter model consists of three smoothing equations and one forecasting equation. Te equation system of the parameters is defned by Te forecasting equation is defned by where 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; α, β, c∈[0,1] denote the smoothing parameter; and k represents the number of moments to be predicted. Te EU adjusts the energy use plan according to the ESO ofer, and the objective function is the diference between the utility of the user and the cost of energy purchase which is defned 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 v e , u e , v h , and u h represent the preference constants for the widely used quadratic utility function [28][29][30]. d e and d h denote the electrical and thermal load demands of the EU. C buy EU denotes the cost of energy purchase by the user as follows: Te demand response of the EU is divided into shifted electrical load l te and curtailable thermal load l rh , while the fxed load of the EU is divided into fxed electrical load l fe and fxed 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.

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

Basic Elements of the Game.
According to the above section, the ESO, the EP, and the EU are all independent agents in the PIES. Te 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, infuencing the price setting of the ESO. Te decisions among the agents are sequential and afect 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]. Te master-slave game model is given by Te participant set, strategy set, and payof function as the three elements of the game model described in (20) can be described as

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. Te set of strategies (L e * , L u * , and L p * ) is the equilibrium solution of the Stackelberg game under the following constraints as follows: 6 International Transactions on Electrical Energy Systems 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 proft by adjusting its set of strategies individually [13]. Te following conditions should be satisfed to ensure the existence of a Stackelberg equilibrium [32]: (1) Te 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 satisfes 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 . Tus, 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 satisfed.
F ESO , F EP , and F EU can be calculated from equations (1), (3), and (15) respectively. Terefore, it can be concluded that the four objective functions are continuous. Accordingly, conditions (3) and (4) are also satisfed.
In summary, the existence of the equilibrium solution of the Stackelberg game is proved.

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.

Te Coyote Optimization Algorithm.
Te main diference between the COA and the most intelligent algorithms is that coyotes are divided into several groups, and the internal social infuences are considered [25]. Te 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. 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.

Coyote Growth in the Group.
Te adaptive ability of coyotes can be evaluated according to the objective function defned 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: Te alpha coyote and other ones naturally infuence 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: Te 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 new soc p,t c � soc p,t c + r 1 δ α + r 2 δ t , δ α � alpha − X rc1 , 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.

International Transactions on Electrical Energy Systems
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. Te 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 frst 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.

Migration of Coyotes.
In the COA, coyotes migrate between groups with probability P e . At frst, 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 diferent groups swap their positions with a probability P e , which is calculated as

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 sufers from the defects such as quickly falling into the local optimum, low convergence rate, and insufcient exploration ability for a high-dimensionality and complex model. Terefore, this manuscript adopts a novel growth approach.

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]. Te 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][35][36]. Te 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 fnding the global optimum. Terefore, 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 efect of the unevenly distributed initial population within the search for an optimum, and makes it easier to escape from a local minimum.
Te tent mapping is given as Te Bernoulli shift transformation is applied to the tent mapping as follows: where mod denotes a function to fnd 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.
Te equation for the initial population generated by the tent mapping is defned by soc p,t c,j � lb j + z n ub j − lb j . (32)

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 signifcant impact on the pack, and if the worst coyote is strengthened, the best and the median coyotes inevitably grow more optimally. Te 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 afect the other coyotes in the group. Terefore, 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. Te beetle antennae search (BAS) is a heuristic algorithm that simulates the search of a beetle for food [37]. Te 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 to 8 International Transactions on Electrical Energy Systems (35) where 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 ftness values of the left and right whiskers, respectively. Te convergence speed of the coyote algorithm and its global search ability can be improved by strengthening the growth of coyotes through BAS [38].

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 signifcantly infuence 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. Terefore, the number of groups is set to 20, with fve 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 efect 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 efect of the global solution and enhances the global search ability. Terefore, 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. Te specifc distribution is shown in Figure 3.
Te fowchart of the BCOA is shown in Figure 4.

Ofspring Generation Employing Genetic Crossover
Operations. Te 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 fnding 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 ftness 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. Te frst period fuctuates slightly around a CR of 0.5, with higher diversity and enhanced exploration ability. Te later period jumps signifcantly around a CR of 0.5, producing new solutions dominated by one of the operations, with reduced diversity and enhanced exploitation ability. Te calculation is as follows: where MaxDT denotes the maximum number of iterations.

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. Terefore, this paper combines the improved coyote algorithm with quadratic programming to propose a distributed solution method. Te 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 proft 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 tarif 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

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

International Transactions on Electrical Energy Systems
Tese 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 p-values are all less than 0.05, and the null hypothesis is rejected for all functions. Tere exists an obvious diference between the four algorithms.

Basic Data. A PIES in Northern
China is employed as a practical example to evaluate the optimal dispatching strategy of the mentioned PIES model. Te general overview of this PIES is shown in Figure 5. Te trading patterns are   Functions [− 100, 100] 0 similar in winter and summer, while only the output of the new energy is somewhat diferent. Tis 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. Te feed-in tarif of the grid is 0.35 ¥/kWh, and the heating price ranges between 0.15 and 0.5 ¥/kWh. Consider that the customer preference coefcients 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.

Algorithm Comparison.
To verify the BCOA optimization for the constructed model, the improved BCOA algorithm is compared with the conventional COA and the improved diferential evolutionary algorithm (ADE) to optimize the revenue objective function of the ESO. Figure 7 depicts the comparison results. Te population of coyotes is adjusted as 100, while the number of coyote groups is chosen as 20 with fve 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 fnds 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 fnally fnds the highest optimal daily revenue, and the optimal solution is superior to the other algorithms. Figure 8 depicts the comparisons of the convergence processes of the ESO, the EP, and the EU benefts. Te benefts 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 diferent convergence trends, verifying the    existence of a game process among agents. Moreover, the ESO benefts show an overall upward trend, while the lowerlevel benefts tend to decline, refecting 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 proft. Te ESO and the EP benefts 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. Te 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 tarif to satisfy the constraint. Te 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.  International Transactions on Electrical Energy Systems Figure 10 shows the operation of the energy storage device. Te electric and thermal loads before and after the demand response are compared in Figure 11. Te EU reduces its energy costs and improves the overall beneft by adjusting the use time of electric cars and washing machines and other equipment. After adopting the demand response strategy, the electric load fuctuation 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. Te 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. Te 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.

Analysis of the Game Results.
To verify the rationality and efectiveness of the proposed optimization strategy, two diferent scenarios are  (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.
Te leader is the ESO, and the followers are only producers and energy users.
Te above scenarios are adopted to calculate the corresponding payofs 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 benefts of each subject under two scenarios is considered separately: (1) Te 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 loadside 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. Te revenue of the producer decreases because the load side cuts some of the load.

Conclusions
Te 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. Tis improves the convergence speed, the global search ability, and the solution stability of the algorithm. Tese 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.
Te 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. Te transferable electric load is shifted to the valley through the demand response of the consumers, improving the consumer surplus and reducing the load fuctuation.
Te 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: Termal storage tank

Data Availability
Te data that support the fndings of this study are available on request from the corresponding author. Te data are not publicly available due to privacy or ethical restrictions.