Multi-Searcher Optimization for the Optimal Energy Dispatch of Combined Heat and Power-Thermal-Wind-Photovoltaic Systems

This paper proposes a novel multi-searcher optimization (MSO) algorithm for the optimal energy dispatch (OED) of combined heat and power-thermal-wind-photovoltaic systems. The available power of wind turbine (WT) units and photovoltaic (PV) units is approximated with the probability density functions of wind speed and solar irradiance, respectively. The chaos theory is used to implement a wide global search, which can effectively avoid a low-quality local optimum for OED. Besides, a double-layer searcher is designed to guarantee fast convergence to a high-quality optimal solution. Finally, three benchmark functions and an energy system with 27 units are used for testing the performance of the MSO compared with nine other frequently used heuristic algorithms. The simulation results demonstrate that the proposed technique not only can solve the highly nonlinear, non-smooth, and non-convex OED problem of an energy system, but can also achieve a superior performance for the convergence speed and the optimum quality.


Introduction
As one of the critical operation tasks of power systems, economic dispatch (ED) is usually employed to minimize the total operating cost via optimally calculating the active power outputs of all the generators to balance the active power demand under various operating constraints [1].However, the conventional ED only considers the electricity energy dispatch of the thermal units, which cannot satisfy the operation requirement of an integrated energy system (IES) with various energy types and sources.
In light of this issue, the ED problem of an IES has attracted extensive investigation.By incorporating customer aggregators' flexible energy demand into the centralized energy dispatch model, a two-level optimization problem [2] was formed in the heat and power IES.In [3], a regional integrated energy system including wind turbines (WT), photovoltaics (PV), gas turbines and battery energy storage was introduced to minimize the operating costs of the system and determine the optimal coordination between the various energy sources.In order to fully exploit the economic and environmental advantages of the system, considering the difficulties of information collection from subareas, a novel decentralized optimal multi-energy flow (OMEF) for large-scale IESs in a carbon trading market was proposed in [4].The improved differential evolution algorithm [5] was designed to obtain the minimum operation cost of the IES, while considering battery lifetime loss.

•
Traditional ED only considers the electricity energy dispatch of power systems, but does not consider the heat energy dispatch of thermal systems.In contrast, the presented OED can achieve electricity and heat energy dispatch simultaneously, while the stochastic characteristics of wind and PV energy are fully taken into account; • Through combination with the chaos theory, the global search ability of the MSO can be effectively enhanced.Besides, a double-layer searcher is designed to accelerate the convergence of MSO, while a high-quality optimal solution can be guaranteed.Compared with nine existing heuristic algorithms, the proposed MSO can search a higher quality dispatch scheme for OED within a shorter computation time.
The rest of the paper is structured as follows: Section 2 details the mathematical model of OED of combined heat and power-thermal-wind-photovoltaic systems.Section 3 gives the detailed operations of the MSO and its application for OED.A discussion and analyses of case studies is given in Section 4. Section 5 concludes the main contributions of this study.

Mathematical Model of OED of Combined Heat and Power-Thermal-Wind-Photovoltaic Systems
The OED framework of combined heat and power-thermal-wind-photovoltaic systems is illustrated in Figure 1.The supply side consists of thermal units, combined heat and power (CHP) plants, heating plants, WT units and PV units.The generated electricity and heat energy are transmitted to the demand side through lines and pipelines.This study aims to minimize the total operating cost of the whole system while satisfying all the constraints, including energy balance constraints, the capacity limits of all energy sources and feasible operating region constraints of CHP units.

Relationship between Wind Speed and Wind Power Output
The power output of a WT unit is mainly determined by the probability distribution parameter of wind speed.In general, the wind speed probability distribution can be described by different models, including Weibull distribution, Rayleigh distribution, log-normal distribution, and normal distribution.Since the two-parameter Weibull model can adjust the parameters to adapt to the periodic variation of wind speed, it is the most frequently used model for wind speed probability distribution, in which the cumulative distribution function can be written as follows [15,16]: where V represents the speed random variable, v is the wind speed, k is the shape factor of the wind speed probability distribution function, which determines the overall shape of the wind speed probability distribution density curve, and c is the scale factor of the wind speed probability distribution function, which reflects the average wind speed of the wind field and can enlarge or reduce the curve, but does not affect the shape of the distribution density curve.
Assuming that the output characteristic of a WT is a simplified piecewise linear function, its power output can be calculated according to the given wind speed [17]: where pwt is the current maximum power points of the WT unit, p r wt is the rated power of the WT unit, vr is the rated wind speed, v is the current wind speed, and vin and vout are the cut-in and cut-out wind speeds, respectively.
The probability distribution of wind speed is converted into the probability distribution of wind power by a linear transformation [18]: Figure 1.Optimal energy dispatch (OED) framework of combined heat and power-thermalwind-photovoltaic systems.

Relationship between Wind Speed and Wind Power Output
The power output of a WT unit is mainly determined by the probability distribution parameter of wind speed.In general, the wind speed probability distribution can be described by different models, including Weibull distribution, Rayleigh distribution, log-normal distribution, and normal distribution.Since the two-parameter Weibull model can adjust the parameters to adapt to the periodic variation of wind speed, it is the most frequently used model for wind speed probability distribution, in which the cumulative distribution function can be written as follows [15,16]: where V represents the speed random variable, v is the wind speed, k is the shape factor of the wind speed probability distribution function, which determines the overall shape of the wind speed probability distribution density curve, and c is the scale factor of the wind speed probability distribution function, which reflects the average wind speed of the wind field and can enlarge or reduce the curve, but does not affect the shape of the distribution density curve.
Assuming that the output characteristic of a WT is a simplified piecewise linear function, its power output can be calculated according to the given wind speed [17]: where p wt is the current maximum power points of the WT unit, p r wt is the rated power of the WT unit, v r is the rated wind speed, v is the current wind speed, and v in and v out are the cut-in and cut-out wind speeds, respectively.
The probability distribution of wind speed is converted into the probability distribution of wind power by a linear transformation [18]: where V and P WT are the wind speed random variable and the wind power random variable, respectively, and T represents a transformation.For the Weibull function, the discrete part of the wind power output can be written as follows: If the output power of the fan is in a continuously changing range, the probability density function of the Weibull model can be described as follows:

Relationship between Solar Irradiance and Solar Power Output
According to the statistics, solar irradiance is subject to Beta distribution over a certain period of time [19], whose probability density function can be written as follows: where r max is maximum solar irradiance, r is the solar irradiance, µ is the mean, and σ is the standard deviation.Generally speaking, the output power of a PV cell is linearly matched with the solar irradiance.If power loss is not considered, and the influence of the ambient temperature is ignored, the output power of a PV unit can be described by [20,21]: where A is the total area of the PV cell and η is the PV cell efficiency.Therefore, the probability density function of the PV cell can be derived as follows: where p pv is the output power of a PV cell and p max pv is the maximum generated power.

Objective Function
In this study, the objective function of OED is the minimization of the total operating cost, which can be calculated as follows: . ., N P ; j = 1, 2, . . ., N c k = 1, 2, . . ., N k ; l = 1, 2, . . ., N l ; m = 1, 2, . . ., N m (15) where C i (P p i ) with i = 1, 2, . . ., N p represents the cost function of the ith thermal unit, C j [P c j , H c j ] with j = 1, 2, . . ., N c represents the cost function of the jth CHP, C k (H h k ) with k = 1, 2, . . ., N k represents the cost function of the kth heat-only unit, C wt,l (P l ) with l = 1, 2, . . ., N l represents the cost function of the lth WT unit, C pv,m (P m ) with m = 1, 2, . . ., N m represents the cost function of the mth PV unit, and N p , N c , N k , N l and N m are the number of conventional units, CHP units, heat-only units, WT units and PV units, respectively.

Thermal Units
In general, the fuel cost of a thermal unit can be expressed by a quadratic function.However, a thermal unit usually has multiple valves for controlling the power output.When the steam admission valves in thermal units are first opened, a sudden increase in losses is observed.This leads to ripples in the cost function, which is known as valve-point loading.Hence, the fuel cost function needs to increase a rectified sinusoidal component by considering this effect.The operating cost of the thermal units can be written as follows [22]: where a i , b i , c i , d i , and f i are cost coefficients of the ith thermal unit, p pmin i represents the minimum power generation limit of the ith thermal unit, and P p i is the electricity energy output of the ith thermal unit.

CHP Units
The CHP cost function is usually expressed as a quadratic polynomial function of the heat and electricity power output of a CHP unit, which appends a coupling coefficient that relates electricity power and heat.Hence, the operating cost of a CHP unit is determined simultaneously by its heat and electricity energy outputs, as follows [23]: where a j , b j , c j , d j , e j and f j are the operating cost coefficients of the jth CHP unit, P c j is the electricity energy output of the jth CHP unit, and H c j is the heat energy output of the jth CHP unit.

Heat-Only Units
Compared with the CHP unit, the heat-only unit can only generate the heat energy.From the result of experiments for the generation of characteristic of heat only units, the operating cost can be described as a quadratic function [24]: where a k , b k and c k are the operating cost coefficients of the kth heat-only unit, and H h k is the heat energy output of the kth heat-only unit.

WT Units
In general, WT units are privately owned and grid operators need to purchase electricity energy from them with a power purchase agreement and a certain payment.The payment consists of three parts, in which the first part is the direct cost of buying wind power.A linear cost function can be assumed for the direct cost of WT, as follows [25]: C d,wt,l (P l ) = d wt,l P wt,l l = 1, 2, . . ., N l (19) where d wt,l is the direct cost coefficients of the lth wind power and P wt,l is the electricity energy output of the lth wind power.
The second part is the penalty cost for underestimation of wind power output due to its high randomness.When the planned wind power is lower than the available wind power, the excess wind energy is wasted and needs to be translated into the penalty cost.The penalty cost for not using all the available wind power will be linearly related to the difference between the available wind power and the planned wind power, as follows [25]: where K ue,wt,l is the underestimated coefficient of the lth wind power, P wt,rate,l is the rated power generation of the lth wind power, and P wt,l is the scheduled power generation of the lth wind power.The third part is the ancillary cost generated by overestimated wind power output.If the scheduled wind power output is higher than the actual available wind power, then the backup power source should be obtained elsewhere to meet the load demand.Similar to the second part, the ancillary cost is also linearly related to the difference between the available wind power and the scheduled wind power, multiplied by the wind power output probability function.Hence, the ancillary cost can be calculated as: C oe,wt,l (P l ) = K oe,wt,l P wt,l 0 (P wt,l − p wt ) f P WT (p wt )dp wt l = 1, 2, . . ., N l (21) where K oe,wt,l is the overestimated coefficient of the lth wind power.
According to the Equations ( 19)-( 21), the total operating cost of a wind plant can be calculated as follows: C wt,l (P l ) = C d,wt,l (P l ) + C oe,wt,l (P l ) + C ue,wt,l (P l ) l = 1, 2, . . ., N l

PV Units
To avoid unnecessary complexity, we assumed the cost function of PV units was similar to WT units.This addition was composed of three parts and calculated as a proposed approach [25]: C ue,pv,m (P m ) = K ue,pv,m P pv,rate,m P pv,m C oe,pv,m (P m ) = K oe,pv,m where d pv,m is the direct cost coefficients of the mth PV power unit, K ue,pv,m is the underestimated coefficient of the mth PV unit, and K oe,pv,m is the overestimated coefficient of the mth PV unit.

Equality Constraints
The energy balance constraints are the most essential requirements for OED, i.e., the electricity and heat energy outputs from the energy suppliers need to be equal to the energy demands, which can be written as follows: where P d and H d represent the total electricity and heat energy demands, respectively.

Inequality Constraints
In order to guarantee a feasible operation, the energy output of each unit should be limited within its lower and upper bounds: where P p i,min and P p i,max are the lower and upper bounds of the ith thermal unit, respectively, H h k,min and H h k,max are the lower and upper bounds of the kth heat-only unit, respectively, P c j,min (H c j ) and P c j,max (H c j ) are the lower and upper bounds of the jth CHP unit, respectively, and H c j,min (P c j ) and H c j,max (P c j ) are the lower and upper bounds of the jth CHP unit, respectively.
It can be found from Equations ( 33) and ( 34) that the heat energy output of the CHP unit will directly influence its lower and upper bounds of the electricity energy output and vice versa.This coupling relationship is usually described by a feasible operating region constraint, i.e., the energy outputs of the CHP unit must be enclosed by the boundary curve ABCDEF, as illustrated in Figure 2. The BC along the boundary line of the area increases the heat energy output of the unit and the electricity energy output decreases.The heat energy output by the CD unit along the line is decremented.
where Pd and Hd represent the total electricity and heat energy demands, respectively.

Inequality Constraints
In order to guarantee a feasible operation, the energy output of each unit should be limited within its lower and upper bounds:  It can be found from Equations ( 33) and ( 34) that the heat energy output of the CHP unit will directly influence its lower and upper bounds of the electricity energy output and vice versa.This coupling relationship is usually described by a feasible operating region constraint, i.e., the energy outputs of the CHP unit must be enclosed by the boundary curve ABCDEF, as illustrated in Figure 2

Multi-Searcher Optimization
Most of the existing metaheuristic algorithms are inspired by nature phenomena, and the proposed MSO is developed from their advantages and the chaos theory.Overall, the MSO consists of a two-layer searcher, in which the first-layer searcher is responsible for a wide global search and the second-layer searcher is implemented for a deep local search.

Chaos Theory
In many practical systems, chaos is a ubiquitous non-linear phenomenon.Chaotic motion can reach every state of a certain scale, according to its regularity and ergodicity, which is better than a simple random search.Therefore, for many optimization algorithms, chaotic theory is frequently used to introduce random parameters in the initialization solution or the replacement algorithm to enhance ergodicity and accelerate global optimal convergence [26].Due to the non-repetitive nature of chaos, the chaotic mapping-based optimization algorithm has more comprehensive search ability than the original one.This feature can make the optimization algorithm effectively escape a low-quality local optimum and avoid the premature convergence.
Here, the logistic mapping formula is used to generate the chaotic variables with chaotic state during the initialization process, thereby, the chaotic iterative formula can be written as follows: where Z t is the tth chaotic vector and κ is a positive real parameter, i.e., the controllable variables of the optimization.

Double Layer Searcher
As shown in Figure 3, in the process of initialization, a series of searchers are randomly distributed within the search area by using Equation (36), which is the so-called global searcher (GS) in the MSO.Surrounded by each global searcher and centered around this searcher within a specified radius, some other searchers are randomly distributed, known as local searchers (LS), which are subject to the global searcher.These LS aims to perceive the surrounding dynamic change process in the vicinity of each GS and then perform a local search.When the searching process of the GS is finished, the LS will initiate a subsequent search.This feature can effectively improve the convergence speed and avoid a low-quality local optimum.The search radius of the GS depends on the distance from the global best (Gbest) and the level of each global searcher.The level will be equal to 0 if a global searcher is finding the minimum cost, and vice versa.Hence, the search radius of the GS is calculated as follows: where rmax and rmin are the pre-set maximum and minimum radii, kGS represents the level of GS, nGS is the number of global searchers, Gbest is described as the global optimal solution, and dist (GS,Gbest) denotes the distance between the global searcher and the global best.

Random Walk Rule
As shown in Figure 4  The search radius of the GS depends on the distance from the global best (Gbest) and the level of each global searcher.The level will be equal to 0 if a global searcher is finding the minimum cost, and vice versa.Hence, the search radius of the GS is calculated as follows: where r max and r min are the pre-set maximum and minimum radii, k GS represents the level of GS, n GS is the number of global searchers, Gbest is described as the global optimal solution, and dist (GS,Gbest) denotes the distance between the global searcher and the global best.

Random Walk Rule
As shown in Figure 4, after searching using the LS, each GS begins to move to a new position by a random motion technique called "random walk".Movement in the random walk is done from a point toward the target position, which leads to the GS being located in the area around the target position.The search radius of the GS depends on the distance from the global best (Gbest) and the level of each global searcher.The level will be equal to 0 if a global searcher is finding the minimum cost, and vice versa.Hence, the search radius of the GS is calculated as follows: where rmax and rmin are the pre-set maximum and minimum radii, kGS represents the level of GS, nGS is the number of global searchers, Gbest is described as the global optimal solution, and dist (GS,Gbest) denotes the distance between the global searcher and the global best.

Random Walk Rule
As shown in Figure 4  To randomly move from one position to another, the searcher should consider two variables, including the distance and direction angle.The distance is determined by a coefficient of the distance between two points, and the direction angle is a random value.The distance between the global searcher and the global best (known as d) is found by using Equation (37).λ is a random value with a distribution of uniform probability in a range from 0 to 1, and θ is a random value with a distribution of uniform probability in a range from -π/4 to π/4.When all global searchers move according to the random walk rule, the new position can be obtained, and the obtained solution will be compared with the current global optimal solution.If the solution quality is higher, the global optimal solution will be replaced.

Constraint Processing
In this paper, all equality constraints (Equations ( 27)-( 28)) are processed by the penalty function method to avoid the "over limit" phenomenon.The equality constraint optimization is transformed into a problem without an equality constraint, and the fitness function is constructed by using a penalty coefficient M, which is as follows: Moreover, the penalty factor M is used for guaranteeing a feasible solution for the proposed OED.In general, a larger M can effectively avoid an infeasible solution, but easily results in a low optimization accuracy.Hence, the penalty factor M should be carefully set via a proper balance between the constraint violation and the optimization accuracy.Through trial-and-error, it was set to be 10 7 for OED.Obviously, when the equality constraint is satisfied, the penalty function is equal to 0, and the value of the fitness function is the cost function.

Execution Procedure
In order to solve a specific optimization problem, the number of optimization variables, constraints, and the objective function should be given to the proposed MSO.In fact, the number of optimization variables will determine the parameters setting of the MSO, while the fitness function is completely dependent on both of the objective function and constraints.Therefore, the proposed MSO can be used for addressing different benchmark test functions or the proposed OED.
At last, the overall execution procedure of the MSO for OED is given in the Figure 5.

•
Step 1: Input system parameters and unit parameters, including the operating cost coefficients of the thermal units, CHP units and heat-only units, the feasible operating region of the CHP units, and the main parameters of renewable energy resources such as wind speed and solar irradiance.Set the algorithm parameters, including the number of GS, the number of LS, the maximum search radius, minimum search radius, and the maximum iteration number.

•
Step 2: Initialize the double-layer searcher based on Equations ( 35)-(36) in chaos theory; and set the current iteration number T = 1.

•
Step 3: Determine the global searcher level and searcher radius with Equation (37).Particularly, the level 0 is assigned to the GS corresponding to the least cost, and the highest level is assigned to the GS corresponding to the highest cost.

•
Step 4: Implement the global searcher and local searcher separately, where the local searcher is surrounded by each global searcher and centered around the searcher within a specified radius.

•
Step 5: Calculate the objective function with Equation (15), which is equal to the total operating cost.

•
Step 6: Evaluate the objective function using the penalty function method from Equation (38), based on the given M.

•
Step 7: Determine the distance and direction based on random walk theory and guide global searcher migration, in which both λ and θ are the random values with a distribution of uniform probability.

•
Step 8: Judge the termination of MSO.If the termination condition is met, then output the economic scheduling result, else turn to step 3 and enter the next iteration.
Through trial-and-error, the main parameters of different algorithms are given in Table 1.All the simulations were undertaken in MATLAB R2016b on a personal computer with an Intel(R) Core TM i7 CPU at 3.4 GHz with 16 GB of RAM.
Through trial-and-error, the main parameters of different algorithms are given in Table 1.All the simulations were undertaken in MATLAB R2016b on a personal computer with an Intel(R) Core TM i7 CPU at 3.4 GHz with 16 GB of RAM.In order to assure a fair comparison, all the parameters of each algorithm have been carefully set based on two rules, as follows: (1) Common parameters: The maximum iteration number and the population size are the common parameters of all of the algorithms, which should be set to the same values for each algorithm.Generally speaking, both a larger maximum iteration number and a larger population size can effectively obtain a higher quality optimum, however, this consumes more computation time.Hence, these two parameters were set to achieve a proper balance between the optimum quality and the computation time.Through trial-and-error, they were set to be 300 and 150 for all the algorithms, respectively.
(2) Independent parameters: Each algorithm has its own parameters, e.g., the crossover probability of GA and the maximum velocity of PSO.These parameters can be determined according to the optimum quality with the given maximum iteration number and population size.Through trial-and-error, the optimal parameters are given in Table 1.

Benchmark Test Function
In order to verify the validity of the algorithm, three different types of functions were selected for numerical experiments.The number of iteration steps was set to 100.
The Sphere function is a simpler unimodal function that examines the convergence speed of the algorithm, where the mathematical model can be written as follows: The global advantage is x i = 0, f (x) = 0.The comparison of the fitness convergence curves in the unimodal function is shown in the Figure 6.It can be seen that the MSO can converge to the global optimum with a fastest convergence rate compared to the other nine algorithms.
Appl.Sci.2019, 9, x FOR PEER REVIEW 13 of 22 The Sphere function is a simpler unimodal function that examines the convergence speed of the algorithm, where the mathematical model can be written as follows: The global advantage is xi = 0, f(x) = 0.The comparison of the fitness convergence curves in the unimodal function is shown in the Figure 6.It can be seen that the MSO can converge to the global optimum with a fastest convergence rate compared to the other nine algorithms.The step function is a discontinuous step function that can be used to verify the validity of the algorithm, which can be described as follows: The global best advantage is xi = 0, f(x) = 0.The comparison results of the fitness convergence curves of the discontinuous step function are shown in the Figure 7. Similarly, the proposed MSO also had the fastest convergence rate among all of the algorithms, while a global optimum could be guaranteed.The step function is a discontinuous step function that can be used to verify the validity of the algorithm, which can be described as follows: The global best advantage is x i = 0, f (x) = 0.The comparison results of the fitness convergence curves of the discontinuous step function are shown in the Figure 7. Similarly, the proposed MSO also had the fastest convergence rate among all of the algorithms, while a global optimum could be guaranteed.
The Rastrigin function is a multimodal function that can be used to test the algorithm's ability to jump out of local optimum, which can be described as follows: The global best advantage is x i = 0, f (x) = 0.The comparison results of the fitness convergence curves of the discontinuous step function are shown in the Figure 8.It can be seen that the global optimization ability of the MSO was more obvious, and other comparison algorithms were trapped in a low-quality local optimum.
The global best advantage is xi = 0, f(x) = 0.The comparison results of the fitness convergence curves of the discontinuous step function are shown in the Figure 7. Similarly, the proposed MSO also had the fastest convergence rate among all of the algorithms, while a global optimum could be guaranteed.The Rastrigin function is a multimodal function that can be used to test the algorithm's ability to jump out of local optimum, which can be described as follows:  In the figures (like Figure 6), the MSO seems to have an advantage immediately in iteration 0 or 1.The main reason for this is that the MSO adopts chaotic mapping for the population initialization, which can reach every state of a certain scale due to its regularity and ergodicity.On the other hand, the other algorithms only adopt a simple random search-based initialization, such that they easily lead to a low optimization efficiency during the initialization process.Hence, the MSO can effectively obtain a higher quality initial best solution compared to other algorithms.

Simulation Model
The test system included 13 thermal units with the valve point effect, 4 CHP units, 4 heat-only units, 3 WT units and 3 PV units.The proposed OED was implemented as a day-ahead or hour-ahead dispatch, thereby the demand profile was determined by the day-ahead or hour-ahead load forecasting according to the historical demand data [33][34][35].The total electricity energy demand was 3837.136MW and total heat energy demand was 615.372MWth.
The main parameters of the thermal units, CHP units and heat-only units are provided in Table 2, Table 3 and Table 4, respectively.The feasible operating region of the CHP unit is shown in Figure 9.The main parameters of renewable energy resources are given in Table 5.All the parameters of the energy sources are steady values, which are usually acquired from the energy suppliers when they connect into the integrated energy system.In this paper, these parameters were referred from the published papers [18,[36][37][38][39][40].The capacities of the WT units in these papers were 130 MW, 94 MW and 94 MW, respectively.Moreover, the capacity of the PV unit was 150 MW, separately.The maximum iteration number of each algorithm was set to 300.In the figures (like Figure 6), the MSO seems to have an advantage immediately in iteration 0 or 1.The main reason for this is that the MSO adopts chaotic mapping for the population initialization, which can reach every state of a certain scale due to its regularity and ergodicity.On the other hand, the other algorithms only adopt a simple random search-based initialization, such that they easily lead to a low optimization efficiency during the initialization process.Hence, the MSO can effectively obtain a higher quality initial best solution compared to other algorithms.

Simulation Model
The test system included 13 thermal units with the valve point effect, 4 CHP units, 4 heat-only units, 3 WT units and 3 PV units.The proposed OED was implemented as a day-ahead or hour-ahead dispatch, thereby the demand profile was determined by the day-ahead or hour-ahead load forecasting according to the historical demand data [33][34][35].The total electricity energy demand was 3837.136MW and total heat energy demand was 615.372MWth.
The main parameters of the thermal units, CHP units and heat-only units are provided in Tables 2-4, respectively.The feasible operating region of the CHP unit is shown in Figure 9.The main parameters of renewable energy resources are given in Table 5.All the parameters of the energy sources are steady values, which are usually acquired from the energy suppliers when they connect into the integrated energy system.In this paper, these parameters were referred from the published papers [18,[36][37][38][39][40].The capacities of the WT units in these papers were 130 MW, 94 MW and 94 MW, respectively.Moreover, the capacity of the PV unit was 150 MW, separately.The maximum iteration number of each algorithm was set to 300.

Discussion and Analysis
Figure 10 shows the convergence of total operating costs obtained by different algorithms and the detailed optimal dispatch strategies of all the energy suppliers obtained by different algorithms are listed in Table 6.It can be seen that the MSO obtains the lowest total operating cost, which can verify that the MSO can effectively jump out of a low-quality local optimum.The MSO started to converge after 114 iterations.

Discussion and Analysis
Figure 10 shows the convergence of total operating costs obtained by different algorithms and the detailed optimal dispatch strategies of all the energy suppliers obtained by different algorithms are listed in Table 6.It can be seen that the MSO obtains the lowest total operating cost, which can verify that the MSO can effectively jump out of a low-quality local optimum.The MSO started to converge after 114 iterations.Figure 11 represents the box and whisker plots of various algorithms that were ran ten times.From the comparison of the box and whisker diagrams, it can be found that the convergence stability of the MSO was not as good as GA and GWO, but the range of variation was still acceptable, and the fitness function was lower than that of GA and GWO.
Figure 11 represents the box and whisker plots of various algorithms that were ran ten times.From the comparison of the box and whisker diagrams, it can be found that the convergence stability of the MSO was not as good as GA and GWO, but the range of variation was still acceptable, and the fitness function was lower than that of GA and GWO.In order to verify the calculation speed of the MSO, Figure 12 shows the average execution time comparison bar graph, obtained by running the various algorithms 10 times.It was found that the execution time of MSO was the shortest among all the algorithms, which was about 4.96 times faster than SA.

Conclusions
In this paper, a novel MSO is proposed for the OED of combined heat and power-thermal-wind-photovoltaic systems, which has the following advantages: In order to verify the calculation speed of the MSO, Figure 12 shows the average execution time comparison bar graph, obtained by running the various algorithms 10 times.It was found that the execution time of MSO was the shortest among all the algorithms, which was about 4.96 times faster than SA.
Figure 11 represents the box and whisker plots of various algorithms that were ran ten times.From the comparison of the box and whisker diagrams, it can be found that the convergence stability of the MSO was not as good as GA and but the range of variation was still acceptable, and the fitness function was lower than that of GA and GWO.In order to verify the calculation speed of the MSO, Figure 12 shows the average execution time comparison bar graph, obtained by running the various algorithms 10 times.It was found that the execution time of MSO was the shortest among all the algorithms, which was about 4.96 times faster than SA.

Conclusions
In this paper, a novel MSO is proposed for the OED of combined heat and power-thermal-wind-photovoltaic systems, which has the following advantages:

Conclusions
In this paper, a novel MSO is proposed for the OED of combined heat and power-thermal-windphotovoltaic systems, which has the following advantages: 1.
The introduction of chaos theory can effectively maintain the non-repetition of the initial population and the diversity of the searcher, which can improve the global search ability of MSO.

2.
The double-layer searcher can effectively make the proposed MSO quickly converge to the best advantage, and avoid the defects that cause the traditional heuristic algorithm to easily fall into a low-quality local optimum.

3.
The superior performance of the MSO was verified by the benchmark test function and a specific engineering optimization of OED, compared with various heuristic algorithms.

Figure 1 .
Figure 1.Optimal energy dispatch (OED) framework of combined heat and power-thermal-wind-photovoltaic systems.
the lower and upper bounds of the jth CHP unit, respectively, and H c j,min (P c j ) and H c j,max (P c j ) are the lower and upper bounds of the jth CHP unit, respectively.
. The BC along the boundary line of the area increases the heat energy output of the unit and the electricity energy output decreases.The heat energy output by the CD unit along the line is decremented.

Figure 2 .
Figure 2. Feasible operation region of the combined heat and power (CHP) units.

Figure 2 .
Figure 2. Feasible operation region of the combined heat and power (CHP) units.

Figure 4 .
Figure 4. Random walk rule of the global searcher (GS).

Figure 4 .
Figure 4. Random walk rule of the global searcher (GS).

Figure 4 .
Figure 4. Random walk rule of the global searcher (GS).

Figure 6 .
Figure 6.Comparison of convergence curves of ten algorithms in single peak test function.

Figure 6 .
Figure 6.Comparison of convergence curves of ten algorithms in single peak test function.

Figure 7 .
Figure 7.Comparison of convergence curves of ten algorithms in a discontinuous step test function.Figure 7. Comparison of convergence curves of ten algorithms in a discontinuous step test function.

Figure 7 .
Figure 7.Comparison of convergence curves of ten algorithms in a discontinuous step test function.Figure 7. Comparison of convergence curves of ten algorithms in a discontinuous step test function.
advantage is xi = 0, f(x) = 0.The comparison results of the fitness convergence curves of the discontinuous step function are shown in the Figure 8.It can be seen that the global optimization ability of the MSO was more obvious, and other comparison algorithms were trapped in a low-quality local optimum.

Figure 8 .
Figure 8.Comparison of convergence curves of ten algorithms in multi-peak test function.

Figure 8 .
Figure 8.Comparison of convergence curves of ten algorithms in multi-peak test function.

Figure 9 .
Figure 9. Feasible operation region of the CHP unit.

Figure 9 .
Figure 9. Feasible operation region of the CHP unit.

Figure 10 .
Figure 10.Convergence process of a total operating costs obtained by different algorithms.

Figure 10 .
Figure 10.Convergence process of a total operating costs obtained by different algorithms.

Figure 11 .
Figure 11.Box and whisker plots of total operating costs obtained by ten algorithms in 10 runs.

Figure 12 .
Figure 12.Average execution time obtained by ten algorithms in 10 runs.

Figure 11 .
Figure 11.Box and whisker plots of total operating costs obtained by ten algorithms in 10 runs.

Figure 11 .
Figure 11.Box and whisker plots of total operating costs obtained by ten algorithms in 10 runs.

Figure 12 .
Figure 12.Average execution time obtained by ten algorithms in 10 runs.

Figure 12 .
Figure 12.Average execution time obtained by ten algorithms in 10 runs. , p i,max are the lower and upper bounds of the ith thermal unit, respectively, H h k,min and H h k,max are the lower and upper bounds of the kth heat-only unit, respectively, P

Table 1 .
Summary of the most important parameters used for different algorithms.

Table 2 .
Parameters of thermal units.

Table 3 .
Parameters of CHP units.

Table 4 .
Parameters of heat-only units.

Table 5 .
Parameters of renewable energy resources.

Table 5 .
Parameters of renewable energy resources.

Table 6 .
Comparative results of optimal solutions obtained by different algorithms.

Table 6 .
Comparative results of optimal solutions obtained by different algorithms.
C k (H h k ) Cost function of the kth heat-only unit C wt,l (P l ) Cost function of the lth WT unit C pv,m (P m ) Cost function of the mth PV unit p pmin i Minimum power generation limit of the ith thermal unit P Heat energy output of the kth heat-only unit d wt,l Direct cost coefficients of the lth wind power P wt,l Electricity energy output of the lth wind power K ue,wt,l Underestimated coefficient of the lth wind power P wt,rate,l Rated power generation of the lth wind power P wt,l Scheduled power generation of the lth wind power K oe,wt,l Overestimated coefficient of the lth wind power d pv,m Direct cost coefficients of the mth PV power unit K ue,pv,m Underestimated coefficient of the mth PV unit K oe,pv,m Overestimated coefficient of the mth PV unit P d , H d Total electricity and heat energy demands, respectively P Lower and upper bounds of the ith thermal unit, respectively H h k,min , H h k,maxLower and upper bounds of the kth heat-only unit, respectively P c j,min (H c j ), P c j,max (H c j ) Lower and upper bounds of the jth CHP unit, respectively H c j,min (P c j ), H c j,max (P c j ) Lower and upper bounds of the jth CHP unit, respectivelyx min , x max Lower and upper limits of the optimization variable, respectively r max , r min Pre-set maximum and minimum radius k GS Random value with a distribution of uniform probability in the range of 0-1