Multi-Objective Optimal Source–Load Interaction Scheduling of Combined Heat and Power Microgrid Considering Stable Supply and Demand

With the development of smart grids, it has become possible to take demand-side resource utilization into account to improve the comprehensive benefits of combined heat and power microgrids (CHP-MGs). In order to improve the benign interaction between the source and the load of the system, the source side decouples the thermoelectric linkage through energy storage devices and improves the system multi-energy supply capacity by introducing various energy flow forms of energy devices. On the demand side, considering the elasticity of electric heating load and the diversity of heating mode, an integrated demand response (IDR) model is established, and a flexible IDR price compensation mechanism is introduced. On this basis, aiming at the optimal stability of supply and demand and the minimum operating cost of the system, a multi-objective optimal operation model of combined heat and power source–load interaction is constructed, taking into account the user satisfaction with energy consumption and the internal equipment load constraints of the system. Finally, an improved multi-objective optimization algorithm is used to solve the model. The analysis of the algorithm shows that the source–load interaction multi-objective optimal scheduling of the cogeneration microgrid considering the stability of supply and demand can effectively improve the stability of supply and demand and the economy of the system.


INTRODUCTION
Combined heat and power microgrid (CHP-MG) based on the concept of multi-energy complementation, energy cascade utilization, and the coordination and optimization of multitype heterogeneous energy subsystems breaks the mode of discrete planning and operation between energy systems and effectively improves the utilization rate of energy (Alomoush, 2019;Blair and Mabee, 2020;Nojavan et al., 2020;Hemmati et al., 2021a;Wang et al., 2021a). In recent years, many scholars have done many pioneering research studies on the modeling and optimal scheduling of the source and load sides of the CHP-MG system (Freeman et al., 2017;Ippolito and Venturini, 2018;Zhao et al., 2019;Ronaszegi et al., 2020;Singh and Kumar, 2020;Jordehi, 2021).
On the source side, Pashaei-Didani et al. (2019) proposed a fuel cell system including a reformer, a hydrogen storage tank, and a fuel cell using natural gas as a raw material. At the same time, the heat production characteristics of hydrogen production from natural gas reforming and the thermoelectric hydrogen coupling characteristics of fuel cells were considered, which effectively reduced the microscopic grid emissions and costs. Hemmati et al. (2021b) improved system operation/recovery and reduced operation/energy costs by optimizing CHP size, location, and equipment operation. In the above research, the optimization target only considers the economic cost and environmental cost and lacks the consideration of the supply-side stability of the system. For example, frequent switching of the start-up state of the equipment or relatively large adjustment of the operating power not only affects the stability and service life of the equipment but also increases the complexity of the transportation, storage, and scheduling of natural gas and other fuels. For source-side optimization scheduling, Liu and Yang (2022) propose a primal-dual based dynamic weight distributed algorithm for multi-objective optimal scheduling of distributed integrated energy system, so as to reduce the operating cost and environmental cost of the system. It is proved that the algorithm has better flexibility, reliability, and adaptability and lower communication burden. Yi et al. (2020) proposed a distributed, neurodynamic-based approach for economic dispatch in an integrated energy system. Compared with other centralized and distributed optimization methods, it is shown that the proposed distributed optimization method has advantages in convergence speed and computational complexity. Anh and Cao (2020) proposed an optimal energy management (OEM) approach that uses smart optimization techniques to achieve optimal hybrid thermoelectric isolation microgrids. Naderipour et al. (2020) proposed the use of particle swarm optimization algorithm to optimize the configuration of the cogeneration system to reduce the operating cost on the basis of considering the maximum allowable capacity. It can be seen that the centralized algorithm is mostly used in a single energy system. For distributed large-scale energy systems, the distributed algorithm has better solving ability than the centralized algorithm.
For the load side, flexible electrical and thermal loads can also be used as potential targets for optimal scheduling of CHP-MG systems. How to quantify and guide the flexible load to participate in each link of the optimal operation of the system and promote the good interaction between source and load has gradually become a research hotspot (Liu et al., 2018;Mohseni et al., 2021;Zhou, 2019;Dranka et al., 2021;Hca et al., 2021;Jiang et al., 2021;Amir et al., 2019) Integrated demand response should adopt price, policy, contract, and other methods to guide users to change their energy consumption habits and optimize the source-load matching relationship of the energy system, so we can obtain greater comprehensive benefits (Aghamohammadloo et al., 2021;Li et al., 2021a;Salehimaleh et al., 2022;Zhang et al., 2019). Lv et al. (2021) have built an IDR multi-time scale optimal scheduling model for power, gas, and heat loads. Zhang et al. (2015) proposed to use a piecewise linear function to characterize the cost of reducing the load and determine the IDR price according to the marginal cost corresponding to various energy consumption in a specific period. Munoz-Delgado et al. (2016) adopted the method of stochastic programming and generated a large number of discrete scenes based on the scene method to represent the uncertainty of DG output and load and convert the uncertainty problem into a deterministic problem. Shao et al. (2020) by incorporating the demand-side response of the energy hub into the integrated power and natural gas system solved the problem that the complex coupling relationship between power and natural gas may reduce the flexibility of system operation. Wc et al. (2021) proposed an optimization technique based on integrated demand response and tolerance of household energy management.
Most of the above studies focus on user load while ignoring the benign interaction between source and load. Load reduction and transfer will greatly affect user satisfaction with energy consumption. However, the conversion and substitution of heating forms (electric heating/gas heating) on the side of the load has little influence on the user's energy experience. And there are few studies to establish a comprehensive evaluation system that considers both the supply and the demand of users. In addition, most of them only focus on user load, ignoring the benign interaction between source and load. For example, renewable energy on the supply side has great volatility and randomness (Yang et al., 2018;Salama et al., 2021), and the negative impact of this factor can be alleviated through IDR. Li et al. (2019) introduced the interaction between energy supplier and customer into system model development. On this basis, considering the characteristics of different time scales of electricity and heat, an event-triggered distributed algorithm is used to optimize scheduling, smoothing real-time load changes and renewable resource fluctuations while maximizing day-ahead social welfare.
To sum up, this paper proposes a multi-objective coordinated optimal scheduling method for cogeneration microgrid source-load that considers the stability of supply and demand. Compared with the existing research, the difficulties and challenges as well as the innovations and contributions made in this paper are as follows: 1) The existing studies on the source side mostly consider economy and environmental protection but lack the consideration of supply stability. In this paper, the CHP-MG model makes use of energy storage equipment to decouple thermoelectric connection and expand the power supply capacity through multi-energy devices. Meanwhile, the supply stability is taken as one of the criteria to evaluate the operation. 2) For the demand response of the load side, most of them adopt the method of load reduction and transfer, which will greatly affect the user's energy satisfaction. Therefore, in addition to the construction of a comprehensive demand response model that includes the electrical load that can be time-shifted and interrupted, electrical loads that can be time-shifted and uninterrupted, the curtailable electrical load, the curtailable heat load, it also considers the integrated demand response model for heat supply mode conversion and introduces a flexible IDR price compensation mechanism. In addition, a comprehensive satisfaction evaluation system was established considering the demand-side user energy deviation degree, IDR compensation price, supply-side heating and power supply flexibility, and other factors. 3) For the problem that most IDR only focuses on user load but ignores the benign interaction between source and load, the fluctuation of renewable energy is smoothed by the IDR of load side. 4) The solution of CHP-MG model has the characteristics of high dimensionality and non-linearity, so an intelligent algorithm is needed to approximate the optimal solution. In addition, compared with the distributed large-scale energy system, the system is of single cogeneration type, and the amount of information and calculation is small, so the centralized optimization algorithm is adopted. In this paper, the improved multi-objective grey wolf optimizer (MOGWO) is used to carry out the source-charge interaction multi-objective optimal scheduling of microgrids, and through example analysis, it is verified that the algorithm can promote the benign interaction between source and load and ensure the stability of supply and demand (SDS) and economy of CHP-MG operation.

GENERAL FRAMEWORK OF THE SYSTEM
The system is roughly divided into two parts, the source side and the load side, and its framework is shown in Figure 1. The source side contains wind energy and CHP-MG systems, which use wind energy, electric energy, and natural gas as sources to provide customers with both electric and thermal energy through the complementary use of multiple energy sources.
In the CHP-MG system, the gas boiler consumes natural gas to generate heat, and the gas turbine consumes natural gas to generate electricity and at the same time generates a large amount of high-temperature flue gas, which is recovered by the heat exchanger to supply heat energy to the load side, thereby realizing the cascade utilization of energy. When the thermal energy generated by a gas turbine and gas boiler is greater than the demand of heating load, the excess thermal energy is stored in the thermal storage equipment. When the thermal energy produced is less than the heating load demand, the thermal energy storage device releases thermal energy. Electricity storage devices are similar to thermal storage devices.
The load side takes the residential load as the core, and the load is divided into two categories: electrical load and thermal load. The electrical load is divided into uncontrollable load and controllable load. Controllable load is further divided into timeshifted and interrupted, time-shifted and uninterrupted. The heat load is divided into the heat supply mode conversion load and the curtailable heat load. Among them, the heat load with convertible heating mode mainly considers the user's behavior of abandoning the heating mode of electric heating and turning to the gas-heat mode of central heating due to the influence of price. The curtailable heat load mainly considers the flexible hot water load.
In order to ensure user experience and prevent the loss of customer groups, this paper establishes a user satisfaction evaluation system. When user satisfaction is lower than the warning value, the scheme should be abandoned. Among them, user satisfaction includes electricity consumption satisfaction and heat consumption satisfaction.

SOURCE-SIDE MULTI-ENERGY COMPLEMENTARY MODEL
In this paper, the CHP-MG contains three energy forms: heat, electricity, and gas. The system energy equipment mainly includes a wind turbine (WT), micro-turbine (MT), and gas boiler (GB); energy storage equipment includes electric energy storage (EES) and thermal energy storage (TES); energy conversion equipment includes heat exchanger (HE) and electrothermal equipment (EE). Its structure is shown in Figure 2. This system has energy interaction with the external power grid (PG).

Wind Turbine
The research shows that the output power of wind turbine (P w ) depends on the wind speed and the rated output power of wind turbine (P r ). In addition, the wind speed obeys Weibull distribution. Therefore, the probability density function (PDF) of P w can be expressed as Here, ε is the scaling factor; k is the shape factor, reflecting the average wind speed or the average output power of the wind turbine in a certain period; and h (v r /v in ) − 1, in which v in and v r are the cut-in and rated wind speeds, respectively. More details of the probabilistic WT model can be found in the study of Yang et al. (2019). In order to reduce the influence of uncertain wind power on system scheduling, this paper takes the mathematical expectation P t WT of P w at each time period as the reference value (Li et al., 2021b).

Micro-Turbine
This paper focuses on the heat and power supply of MT, ignoring the effect of environmental changes on power generation and combustion efficiency. Its heat generation power Q t MT is Here, P t MT is the electric generation power of MT in period t and η MT , η L are the power generation efficiency and heat loss coefficient of MT.
Operating power constraints and climbing constraints are Here, P MT,max /P MT,min and dP MT,max /dP MT,min are the upper and lower limits of MT operating power and climbing rate and I t MT is the status flag bit of MT, which is 1 when running and 0 when stopped.

Heat Exchanger
The waste heat discharged from the MT is converted into usable heat energy by the HE: Here, Q t HE is the heat production power of HE in period t and η MT is the thermal efficiency of HE.

Gas Boiler
The GB generates heat energy by consuming natural gas to supplement the heat load when HE and TES are insufficient. The relationship between output power and input power is Here, Q GB,max /Q GB,min and dQ GB,max /dQ GB,min are the upper and lower limits of GB operating power and climbing rate.

Electrothermal Equipment
Electrothermal equipment consumes electric energy for heat production, such as air conditioners and electric boilers. The relationship that needs to be satisfied is Here, Q t EE is the heat generation power of EE in period t; η EE is the thermal efficiency of EE; and Q EE,max , Q EE,min are the upper and lower limits of EE heat generation power.

Energy Storage Equipment
In this paper, the energy storage device is mainly used to decouple the complex electrical and thermal connection, so that the CHP-MG system can get rid of the traditional mode of "determining electricity by heat" and "determining heat by electricity." On this basis, the system runs in the direction of supply and demand stability and optimal economy.

Electricity Storage Equipment
EES is often used in microgrids to achieve peak-shaving and valley-filling of electrical loads, smooth operation of generation equipment, and reduction of operating costs. The relationship between power and energy storage capacity during the operation of EES is shown below.
Energy storage constraint: Here, E t EES is the storage capacity of EES in period t; P t−1 EES,c is the charging power of EES in period t-1; P t EES,d is the discharge power of EES in period t; η c /η d is the charge/discharge coefficient of EES; and E EES,max /E EES,min is the maximum/minimum capacity of EES.
Power constraint: Here, P EES,c,max /P EES,c,min is the maximum/minimum charging power of EES; P EES,d,max /P EES,d,min is the maximum/minimum discharge power of EES; and I EES,c /I EES,d is the charge/ discharge status flag of EES, which is 0 when stopped and 1 when running. The mutual exclusion constraint needs to be satisfied, that is, Here, dP EES,c,max /dP EES,c,min is the maximum/minimum ramp rate for charging EES and dP EES,d,max /dP EES,d,min is the maximum/minimum ramp rate of EES discharge.

Heat Storage Equipment
An energy storage device is a storage device for energy, not just for electrical energy. Generally, the peak of electricity consumption occurs during the day, and the peak of heat consumption occurs in the morning and evening, resulting in a mismatch between the electric heat output and the electricity load. Therefore, TES is used in this paper to realize the translation of heat time and solve the difference in the time of electric-thermal load. In order to simplify the model, without considering the influence of environmental factors on heat selfloss rate, the TES constraint relationship is similar to that of EES, as shown below. Energy storage constraint: Here, E t TES is the heat storage of TES in period t; Q t−1 TES,c is the heat storage of EES in period t-1; Q t TES,d is the exothermic power of TES in period t; δ c /δ d is the heat storage/release coefficient of TES; δ is the self-loss coefficient of TES; and E TES,max /E TES,min is the maximum/minimum capacity of TES.
The power constraint and climbing constraint satisfied by TES are consistent with those of EES and will not be discussed in this paper.

COMPREHENSIVE DEMAND RESPONSE MODEL OF LOAD SIDE
The IDR behavior of users on the load side of the CHP-MG system is divided into the response of electric load demand and the response of thermal load demand. It is assumed that the user's IDR mode adopts the following strategies: 1) In the case of satisfying user satisfaction, the thermal load can be flexibly responded by adjusting the demand of electric/gas heat load and the reduction of heat load. 2) In the case of satisfying user satisfaction, the electrical load can be flexibly responded by shifting and interrupting and reducing the electrical load.
Among them, the flexible electrical loads, that is, electrical loads that participate in response, can be divided into the timeshiftable and interruptible electric load, the time-shiftable noninterruptible load, and the curtailable electric load. Flexible heat loads, that is, the heat loads that participate in response, are reducible heat loads. In addition, for flexible heat load, this paper considers the user's choice of heating mode.

Time-Shiftable and Interruptible Load
Time-shiftable and interruptible load means that the load's power consumption time can be shifted and interrupted according to system requirements. The working duration and the number of interruptions of different electrical appliances are different. Assuming that the number of devices participating in demand response is n b,int , the charging time of device n is T b,n , the number Frontiers in Energy Research | www.frontiersin.org May 2022 | Volume 10 | Article 901529 5 of interrupts is q b,n , and the transferable range is [t b,n,s , t b,n,e ], the following relationship must be satisfied: Here, q b,n,max is the maximum number of interruptions of device n; t b,n,si /t b,n,ei is the actual time point of the ith start/end of device n; and U b is the set of working time periods after the device n is shifted and interrupted. The electric load L t b0 /L t b expression before/after the translation and interruption is Here, P b,n is the rated power of the device n; L t b0,n is the ideal response of equipment n to the electrical load demand in time period t; and L t b,n is the response electrical load demand of device n after time shift and interruption in time period t.

Time-Shiftable and Non-Interruptible Power Loads
Time-shiftable and non-interruptible power load means that the load's power consumption time can be shifted according to system requirements, and the working duration of different electrical appliances is different. The modeling process of the electrical loads L t a0 , L t a , etc., before/after translation is only a case when the number of interruptions in 4.1.1 is 0, so this section will not repeat them.

Electric Load Can Be Reduced
In the case of satisfying user satisfaction, the electric load can be appropriately reduced, and the electric load curve of the system load side can be smoothed. The correlation is as follows: Here, ΔL t c is the electrical load difference before/after reduction in period t; L t c is the electrical load after reduction in time period t; L t c0 is the electric load participating in the reduction in time period t; ε t c,e is the electric load reduction coefficient after reduction in time period t; ε c,e,max /ε c,e,min is the maximum/ minimum load reduction factor; and I c,e is the sign indicating whether the electric load is reduced, which is 1 when it is reduced.

Response of Heating Mode
The peak period of electrical loads such as air conditioners and electric furnaces is roughly the same as the peak period of grid power supply. Under the stimulation of price, some users gave up the heating method of electric heating and switched to the gas heating method of central heating. The relevant relationship is as follows: Here, ΔH t a is the heat load difference before/after the response by the heating method in the period t; H t a0 is the original gas heat load participating in the response of the heating mode in time period t; H t a /D t a is the gas/electric heating load demand after the response of the heating mode in the period t; ε t a,h is the replacement coefficient in time period t; ε a,h,max /ε a,h,min is the maximum/minimum replacement coefficient; and I a,h is the sign indicating whether the heating mode response occurs to the heat load, which is 1 when the response occurs.

Reduced Heat Load
The main consideration of the heat load that can be reduced is the hot water load. Under the lowest water temperature T w,min acceptable by the user, appropriate reduction of the heat load can smooth the heat load curve on the load side of the system. The correlation is as follows: Here, C w , ρ w are the specific heat capacity and density of water; v t w is the user's hot water demand in time period t; T w0 is the initial water temperature (15°C); T t w is the ideal average temperature of hot water for the user in time period t; H t b,min is the minimum hot water power required by the user in time period t; ΔH t b is the hot water load difference before/after the reduction in period t; H t b is the hot water load reduced in time period t; H t b0 is the hot water load participating in the reduction in time period t; ε t b,h is the hot water load reduction coefficient in period t, whose value is between 0 and 1; and I b,h is a flag indicating whether the hot water load has been reduced, which is 1 when the reduction occurs.

IDR Price Compensation Mechanism
Price compensation encourages users to actively participate in demand response. Unlike the traditional fixed-ratio price compensation, this paper uses a flexible price compensation with the compensation cost f c shown as follows: Here, c t ab,e , c t c,e are the unit compensation price for participating in time shift and interruption and reducing the response electric load in time period t; c t a,h , c t b,h are the unit compensation price involved in adjusting the electric heating load and reducing the heat load response during period t; and λ t ab,e , λ t c,e , λ t a,h , and λ t b,h are auxiliary variables with values between 0 and 1.

SOURCE-LOAD COORDINATION AND INTERACTIVE OPTIMIZATION SCHEDULING
The optimal operation model of CHP-MG with IDR, decoupling the thermoelectric connection through the energy storage device, considering the multi-energy complementary characteristics, jointly formulates the optimal output plan of each coupling device from both sides of the source and the load, so as to improve the economy and supply and demand of the microgrid stability.

Economics
When considering CHP-MG economy F 1 , both the supply side and the demand side need to be considered. The energy supply cost f CG on the supply side includes the purchase cost of input energy and the aging cost of equipment, while the economy f CD on the demand side is mainly reflected in the compensation cost of users participating in IDR: 1) Large grid electricity purchase cost f grid : Here, P t grid and c t buy are the purchase power and unit price of electricity from the microgrid to the large grid in time period t, respectively.
2) EES charge and discharge aging cost f BT : Here, c t EES is the unit aging price of EES in time period t.
3) TES charge and discharge aging cost f BT : Here, c t TES is the unit aging price of TES in time period t. 4) MT natural gas cost f MT : Here, c t gas is the price of natural gas in time period t and LHV gas is the low calorific value of natural gas, which is 9.7 kWh/m 3 . 5) GB natural gas cost f GB : 6) The cost of natural gas for electric heat to gas heat f Ha is Here, λ Ha is the conversion factor of electricity-heat load.

Stability of Supply and Demand
While pursuing economy, CHP-MG system operators should also focus on providing a stable energy supply for users. On the supply side, the stability of MT, GB, and external power grid is mainly considered. The stability of MT and GB output can reduce the frequent regulation of equipment, reduce the failure rate, and prolong the service life. In addition, as the main consumption equipment of natural gas, the MT and GB can also improve the input stability of natural gas and reduce the complexity of transportation and management. Keeping the power interaction with the external grid stable improves the stability and safety of the system while reducing the stress caused by the CHP-MG system to the outside.
On the demand side, it mainly measures the discrete degree of the difference between the user's electricity load and the wind power load (LW). Because wind power load is easily affected by environmental factors, its load curve is highly discrete, which increases the difficulty in hardware and software regulation of microgrid. The traditional load prediction curve is usually "peakcutting and valley-filling" by an energy storage device, but the load prediction curve often represents a larger possible state of wind turbine output while ignoring other output possibilities. In this paper, the mathematical expectation P t WT of P w in time period t is taken as the "equilibrium point" of output possibility of wind turbine in time period t, so that the optimized operation process of the system works near the "equilibrium point" of wind turbine. In addition, the CHP-MG system adjusts the electrical load curve through IDR and smoothes the difference curve between user power load and wind power load.
In this paper, the weighted value of the change rate of the output curve of each load or equipment is used to evaluate the supply and demand stability coefficient F 2 , which is defined as Here, σ MT , σ GB , σ grid , and σ l−W are the standard deviation of MT power generation, GB heat generation power, external grid output power, and LW; σ σ is the penalty function of equilibrium within the stability of supply and demand, whose value is the standard deviation of matrix σ; and a 1 , a 2 ... a 5 are coefficients greater than 0. The relevant parameters must meet the following conditions: Here, σ MT,max , σ GB,max , σ grid,max , σ l−W,max are upper bounds on σ MT , σ GB , σ grid , σ l−W .

Threshold
Whether the CHP-MG can provide a comfortable energy experience for customers is also an important measure indicating whether the microgrid is competitive in the energy market. To quantify this type of auxiliary service, user satisfaction is evaluated in terms of load deviation rate, IDR price compensation factor, and adjustable rate of electricity and heat production. In this paper, user satisfaction is divided into power supply satisfaction f UC,e and heating satisfaction f UC,h , which need to meet f UC,e ≥ UC e, min , f UC,h ≥ UC h, min .
Here, UC e,min /UC h,min is the lower limit of power supply/thermal satisfaction.

Satisfaction of Power Supply
The CHP-MG system characterizes power supply satisfaction (ESS) by the electric load deviation rate, the power generation adjustable rate, and other indicators. Among them, price compensation for IDR can reduce the negative impact on user comfort caused by electrical load deviation. Its power supply satisfaction is defined as 1) Adjustable rate of electricity production θ e : θ e cθ e,do + (1 − c)θ e,up , c ∈ (0, 1) .
Here, θ e,do /θ e,up indicates the power generation can increase/cut the rate, respectively. θ e,do is given as Here, θ t e,do is the amount of power generation that can be reduced in time period t. The main consideration is the amount of electricity that can be reduced by MT and EES on the premise of meeting the upper and lower limit constraints and climbing constraints described above.
θ e,up is the amount of electricity that can be added in the t period, similar to θ t e,do . The main consideration is the amount of electricity that MT and EES can add on the premise of meeting the upper and lower limit constraints and climbing constraints described above.
2) The electrical load deviation rate d e is related to κ t ab,e (electric load deviation rate due to time shift and interruption), κ t c,e (electric load deviation rate due to curtailment), λ t ab,e , λ t c,e . It is defined as The influence of the deviation rate of power generation on user satisfaction is often dominated by the inferior party, so the exponential variable is introduced as the weight in the above formula. If d ab,e or d c,e is obviously too large, its weight will increase sharply.

Satisfaction With Heat Supply
In the CHP-MG system, heat supply satisfaction (HSS) is characterized by the heat load deviation rate and heat production adjustable rate. Among them, price compensation for IDR can reduce the negative impact on user comfort caused by thermal load deviation. Satisfaction with heat supply is defined as 1) Heat production adjustable rate θ h : Here, θ h,do /θ h,up is the curtailment/increment of power generation in period t. The main consideration is the amount of power that can be cut/increased by GB and TES with the need to satisfy the upper and lower constraints and climbing constraints described above.
2) The heat generation deviation rate d h is related to d a,h (the heat load deviation rate due to the selection of the heating method), d a,h (the heat load deviation rate due to the reduction), λ t a,h , and λ t b,h . Its definition is consistent with the heat generation deviation rate and will not be repeated in this paper.

Constraints
The CHP-MG system constructed in this paper includes the demands of both thermal and electric loads and needs to satisfy the electric and thermal power balance and the large grid interaction power constraints in addition to the equipment operation constraints.
1) Electric power balance: Here, P t WT generates power of WT in time period t and P load is the basic electrical load required by users on the load side in time period t.
2) Thermal power balance: Here, P phot is the basic heat load required by users on the load side in time period T.
3) Large grid interactive power constraints: Here, P grid,max /P grid,min and dP grid,max /dP grid,min are the upper and lower limits of the power purchase and ramp rate of the large grid.
6 IMPROVED MULTI-OBJECTIVE GREY WOLF OPTIMIZER

Multi-Objective Grey Wolf Optimizer
The GWO (Mirjalili et al., 2014) is a new intelligent algorithm proposed by Mirialili et al., in 2014. In 2015, on this basis, an MOGWO (Mirjalili et al., 2016) was proposed. The optimization process of MOGWO is divided into the following steps:

1) Social class stratification
In the MOGWO, the objective function value of grey wolf individuals in each iteration process is calculated, the nondominated individuals are determined, and the excellent population Archive is updated. And the α, β, δ wolves are determined by the roulette method, and the rest are ω wolves, so as to realize the classification of the grey wolf population.

2) Surrounded
The position of the wolf represents a problem solved by the algorithm. The behavior of being surrounded by grey wolves during hunting is defined as Here, t is the current iteration number; D is the distance between the individual and the prey; A and C are synergy coefficient vectors; X p (t) represents the position vector of prey; and X(t) represents the position vector of the current grey wolf. A and C are defined as follows: Here, a is the convergence factor. For the specific definition, see Eq. 41. r 1 and r 2 are random variables between 0 and 1. ω wolf surrounds its prey under the leadership of α, β, δ wolves, and its position update formula is as follows: Here, D α , D β , D δ are the distances between the α, β, δ wolves and the current individual grey wolf, respectively; X α (t), X β (t), X δ (t) are the current position vectors of α, β, δ wolves, respectively; and C 1 , C 2 , C 3 are random variables.
Here, X(t + 1) is the updated position vector for ω wolf.

3) Attack prey
When the prey stops moving, the wolf attacks it to complete the hunt. The timing of the grey wolf attacking its prey is controlled by the value of a. As shown in Eq. 41, the size of a decreases linearly with the decrease of the number of iterations. When a drops to 0, that is, the maximum number of iterations, the encirclement is stopped and the prey is attacked, that is, the iteration is stopped, and the result is output: Here, t max is the maximum number of iterations.

Improved Multi-Objective Grey Wolf Optimizer
The optimization problems of CHP-MG systems are non-linear optimization problems with complex constraints and high solution dimensionality. When using the original MOGWO to solve such problems, it is easy to converge prematurely and fall into local optimum. To remedy such problems, Mirjalili et al. (2016) proposed a reduced-order aggregate model based on the balanced truncation approach and Wang et al. (2021b) proposed a reduced-order small-signal closed-loop transfer function model based on Jordan continued-fraction expansion. However, the dimensionality reduction method mentioned above is difficult to apply under the condition that the system input variables are closely related to the target value or threshold value. Therefore, the improvement of the multi-objective grey wolf optimization algorithm in this paper mainly focuses on the search process and form of wolves (Taha and Elattar, 2018;Rui et al., 2020a), as follows: 1) The MOGWO algorithm's insufficient exploration ability affects the global search ability and local convergence ability of the whole algorithm. The size of control parameter a in Eq. 41a affects the overall detection distance of the wolves. In this paper, the linearly reduced a is changed to the cosine function form of Eq. 42. In the early stage, the decline speed of a slowed down, which improved the global search ability of the algorithm. The local convergence of the algorithm can be ensured by slowing down the later descent speed: a 2 · cos k 2πk max (42) 2) The original MOGWO, based on the crowding of the Archive population, uses a roulette wheel to pick wolves α, β, δ (the less crowded wolf is more likely to be chosen). However, when the three leading wolves are all local optimal solutions, the algorithm still has the risk of falling into local optimal solutions. In this paper, the wolves are divided into three groups by the FCM clustering algorithm, the non-dominant individuals in each group are determined, and Archive1~Archive3 are updated. The leader wolf from each group was selected as Wolf α,Wolf β, and Wolf δ by the above roulette method. This method increases the diversity of the population and avoids falling into local optimum.
3) The three groups of wolves (except the leader wolf) are randomly drawn with probability P to select the hunter, and its position is updated by Eq. 37 (displacement update formula of grey wolf algorithm). The cooperative search ability of each group of pursuers is ensured. The other wolves (except the leader wolf) in the unselected groups were used as vigilant, and the target was the leader wolf in each group to ensure the search ability of the vigilant group within the group. The group update was carried out through Eq. 43 (displacement updating formula of particle swarm optimizer) (Kennedy and Eberhart, 1995): Here, w d i is the inertia weight; w max /w min is the preset maximum/minimum inertia coefficient; f d average , f d min are the average fitness and minimum fitness of all particles in the dth iteration; V k i is the velocity of the kth iteration of the current particle; X k i is the position of the kth iteration of the current particle; P best is the best solution for the individual; G best is the best solution for the population; c 1 and c 2 are learning factors; and rand 1 and rand 2 are random numbers between 0 and 1.
The steps of applying the improved multi-objective grey wolf algorithm to CHP-MG system scheduling optimization are as follows: Step 1: Set the parameters and data of the CHP-MG optimization model.
Step 2: Set parameters such as the number of grey wolves, the maximum number of iterations, and the excellent population Archive and initialize the grey wolf population.
Step 3: The grey wolf population (CHP-MG scheduling scheme) is grouped by the FCM clustering algorithm, and the α, β, and δ wolves are selected according to the improved method at (2), that is, the optimal scheme.
Step 4: The remaining wolves in each group except the leader wolf are randomly selected with probability P, that is, some schemes except the optimal scheme in each group of scheduling schemes are selected, and their positions are updated through the displacement update formula of the grey wolf algorithm.
Step 5: For each group of wolves that are not selected except the head wolf in step 4, take the head wolf of each group as the objective function (the optimal scheme of each group) and update their positions in groups through the displacement update formula of particle swarm optimization.
Step 6: Calculate the objective function value of the whole grey wolf population, filter the threshold, determine the nondominated individuals, and update Archive.
Step 7: Judge whether the maximum number of iterations has been reached. If yes, output Archive. It is the Pareto solution set of the optimal scheduling scheme of CHP-MG system. On the contrary, skip to step 3 until the termination condition is met.

Basic Data
In order to verify the effectiveness of the model and algorithm described in this paper, example simulation was carried out by referring to the system micro-source device parameters, energy storage device parameters, energy purchase and sale price, and load data in the examples in the literature (Li et al., 2021b;Yang et al., 2016;Wang et al., 2020;Li et al., 2021c;Li et al., 2021d;Lu et al., 2021;Zhang et al., 2021). The total scheduling time is 24 h, and the unit scheduling time is 1h.
The prediction curves of basic heat load, hot water load, electric load, and response of heating mode are shown in Supplementary Figure 1. The wind turbine output prediction curve is shown in Supplementary Figure 2. The timeshare price of energy (power grid sales price and natural gas price) is shown in Supplementary Figure 3. Refer to Supplementary Table 1 for the parameters of micro-source equipment. See Supplementary Table 2 for the parameters of the energy storage device. Flexible time-shifting noninterruptible electrical load and time-shifting interruptible electrical load mainly consider the time-shifting characteristics and interruptible characteristics of washing machine (WM), dishwasher (DW), rice cooker (RC), dryer (DY), and sweeper (SE). A total of 500 households in this area are divided into five categories: A, B, C, D, and E, according to their energy consumption habits. The operating characteristics of various household appliances are shown in

Scene
Flexible electrical load Flexible heat load

Analysis of Simulation Results
In order to further verify the validity of the model, four operating modes are selected for comparison, as shown in Table 1.
Using MATLAB software to calculate the above four scenarios, in the improved multi-objective grey wolf algorithm, the total population is 100. The Pareto solution set convergence curves of the four scenarios are shown in Figure 3, and the optimization results are shown in Table 2.
It can be seen from the analysis in Figure 3 that the improved MOGWO has better convergence characteristics than the original MOGWO algorithm, and its search accuracy is also high. At the same time, due to the improvement of the update strategy of Pareto solution set, its diversity and distribution characteristics are better.
According to the analysis in Table 2, compared with traditional optimized operation (Case1), considering flexible electrical load and thermal load significantly improves the economy and stability of supply and demand of the system.
Considering the flexible electrical load (Case2), the total cost is only reduced by $37.79 compared to Case1, while the stability of supply and demand is improved by 52.20%. This is because the flexible electrical loads of Case2 can be timeshifted but not interrupted, and the time-shifted and interrupted electrical loads can improve the stability of supply and demand. However, the total load and energy supply form have not been greatly improved, so the cost cannot be effectively reduced. As for the curtailable electrical load, considering the user experience, the total amount of its participation in the demand response is small, and it can only slightly improve the system economy and the stability of supply and demand.
Considering the flexible heat load (Case3), the total cost is reduced by $158.38 compared to Case1, while the supply-demand stability factor is only increased by 9.77%. This is because the response of Case3's flexible heat load heating mode turns electric heating to gas-heat central heating, which changes the energy supply form and reduces costs. For the curtailable electrical load, while adjusting the heat load curve to improving the stability of supply and demand, it avoids unnecessary waste of heat energy, thereby reducing costs. But similar to curtailable electrical loads, the effect is limited.
As for Case4, the advantages of Case2 and Case3 are taken into account because of both flexible electrical load and flexible thermal load, so that the economy and stability of supply and demand are optimal. The total cost was reduced by $195.71, and supply and demand stability improved by 67.99%.

Analysis of Electric Energy Operation Results
Select the scheme with the lowest total cost in Case1~Case4 and analyze the power operation results through comparison.
The power operation results of Case1 are shown in Figure 4A. During 1-6 and 23-24 electricity price valley/gas price peak period, electricity load is mainly borne by wind energy, large grid, and electric energy stored in EES. During the electricity price peak/gas price valley period from 7 to 12 and 19 to 22, the electricity load is mainly borne by gas turbines, wind energy, and large power grids, and electricity storage equipment releases electricity. During the electricity price average/gas price peak period from 13 to 18, the electricity load is mainly borne by gas turbines, wind energy, and large power grids, and the electricity storage equipment stores excess electricity. The output of the gas turbine decreases relative to the peak period of electricity price/ valley gas price.
The power operation results of Case2~Case4 are shown in Figures 4B-D. Compared with Case1, in the periods 7-24, the electrical load borne by the gas turbine gradually increases, and the electrical load borne by the large power grid gradually decreases. During the periods 1-6, the stored electric energy of the power storage equipment decreased significantly, and during the periods 21-24, the released electric energy increased significantly.

Analysis of Thermal Energy Operation Results
Select the scheme with the smallest total cost in Case1~Case4 and analyze the thermal energy operation results by comparison.
The power operation results of Case1 are shown in Figure 5A. During 1-6 and 23-24 electricity price valley/gas price peak period, the heat load is mainly borne by gas boilers and electric heating equipment. During the electricity price peak/ gas price valley period from 7 to 12 and 19 to 22, the heat load is mainly borne by heat exchangers, gas boilers, and electric heating equipment. During the electricity price average/gas price peak period from 13 to 18, the heat load is mainly borne by heat exchangers, gas boilers, and electric heating equipment. Compared with the electricity price peak/gas price valley period, the output of the heat exchanger is reduced.
The thermal operation results of Case2~Case4 are shown in Figures 5B-D. Compared with Case1, in the period of 7-24, the heat load borne by the heat exchanger gradually increases. During the whole period, the heat load borne by the gas boiler increased significantly.

Economic Analysis
In order to further analyze the influence of the improvement mentioned above on system economy, this paper conducts economic analysis by comparing Case1 to Case4.
The scheme with the lowest total cost from Case1 to Case4 was selected for economic analysis by comparison. The optimization results are shown in Table 3.
Combining Table 3 and Figure 6, the total cost of Case2 is reduced by $37.79 compared to that of Case1, of which the natural gas cost of the gas turbine is increased by $85.80 and that of the gas boiler is reduced by $6.57. The electricity purchase cost of the large grid decreased by $153.67; the aging costs of storage device were reduced by $0.03, which was negligible; an additional IDR compensation cost of $36.67 was incurred. This is because the load side of the Case2 system responds to demand by introducing flexible electrical loads to "cut peaks and fill valleys" for loads. It reduces the probability that the sourceside gas turbine cannot supply energy in time due to factors such as ramp constraints, thereby increasing the natural gas cost of the gas turbine and indirectly reducing the power purchase cost of the large power grid. Compared with that of Case1, the total cost of Case3 decreased by 158.36 dollars, of which the natural gas cost of gas turbine and gas boiler increased by 187.36 dollars and 33.42 dollars, respectively, and the power purchase cost of large power grid decreased by 406.64 dollars. This is because the response of the load-side heating mode of the system makes some users change from electric heating to more economical gas-heat central heating. Therefore, the output proportion of gas-fired boiler and heat exchanger of natural gas equipment on the source side increases, which directly or indirectly increases the natural gas cost consumed by a gas turbine and gas-fired boiler and reduces the power purchase cost of large power networks.
Compared with Case1~Case3, Case4 has the best economic indexes due to taking into account the above two advantages. The total cost decreased by $195.79. The cost of gas for the gas turbine and gas boiler increased by $195.71 and $34.18, respectively. The purchase of electricity by the large grid was reduced by $515.80. In addition, compared with Case1, Case4 increased the satisfaction of electricity consumption by 0.0715 and decreased the satisfaction degree of heat consumption by 0.2224. Among them, the improvement of power consumption satisfaction is because the load-side demand response makes users deviate from the power consumption time and the total amount but also improves the power generation adjustable rate. The reduction in heat satisfaction is smaller, similar to electricity satisfaction.

Analysis of Supply and Demand Stability
In order to further analyze the impact of the above improvements on the stability of supply and demand of the system, this paper  conducts an analysis of the stability of supply and demand by comparing Case1~Case4. The scheme with the smallest supply and demand stability coefficient of Case1~Case4 is selected, and the stability of supply and demand is analyzed by comparison. The optimization results are shown in Table 4.
As shown in Table 4, the overall supply and demand stability of Case2 is improved by 52.20% compared to that of Case1. The dispersion degree of LW, gas turbine, gas boiler, and large grid load curves decreased by 61.65%, 14.64%, 0.96%, and 75.51%, respectively. This is because the load side of the system reduces the dispersion degree of load curves of gas turbine, gas boiler, and large grid by introducing flexible electrical loads for demand response. This is because the load side of the system performs demand response by introducing flexible electrical loads. While reducing the discrete degree of load curves of LW, gas boilers, and large power grids, more gas turbines are involved in power supply, which reduces the output stability of natural gas power generation equipment.
Compared with that of Case1, the overall supply and demand stability of Case3 is improved by 9.77%. The dispersion degree of LW and large grid load curves decreased by 22.53% and 10.90%, respectively, while the dispersion degree of gas turbine and gas boiler load curves increased by 50.13% and 22.94%, respectively. Obviously, the response of load side heating mode of the system not only reduces the dispersion of LW and large power grid load curve, but also makes natural gas heating equipment more involved in energy supply, thus improving the stability of its output.
Compared with that of Case1, the total stability of supply and demand of Case4 increased by 67.99%. The dispersion degree of LW, gas boilers, large power grids, and gas turbines decreases by 78.11%, 31.74%, 24.41%, and 88.92%, respectively. This is because the system load side avoids the adverse effects of Case1 and Case2 by introducing flexible electrical load and flexible thermal load to carry out demand response. Among them, the smoothness of LW curve was greatly improved, mainly reflected in the period from 7 to 24, as shown in Figure 7A. The smoothness of the gas turbine output curve has been greatly improved, mainly reflected in the periods 10-12, 15-19, and 1. The output curve is shown in Figure 7B. The smoothness of the output curve of the gas boiler has been improved to a certain extent as a whole, and its output curve is shown in Figure 7C. The smoothness of the output curve of the large power grid has been greatly improved as a whole, and the output curve is shown in Figure 7D.  This paper proposes a multi-objective optimal source-load interaction scheduling of combined heat and power microgrid considering stable supply and demand. After comparing and analyzing the actual results, the conclusions are as follows: 1) The source side uses energy storage equipment to decouple the thermoelectric connection and introduces energy equipment to ensure the heating and power supply capacity of the system. 2) Flexible thermal and electrical loads are introduced on the load side for demand response, and an IDR compensation mechanism is established. It can effectively improve the economy of the system and the stability of supply and demand, in which the total cost is reduced by 11.65% and the stability of supply and demand is increased by 67.99%. 3) The comprehensive satisfaction evaluation system established by considering factors such as user energy deviation, IDR compensation price, supply-side heating, and power supply flexibility effectively guarantees the user's comfort when participating in demand response. 4) The improved multi-objective grey wolf optimization algorithm adopted realizes the multi-objective optimal scheduling of the source-load interaction, promotes the benign interaction between the system source and the load, and ensures the supply and demand stability and economy of the CHP-MG operation.