Optimal Non-Convex Combined Heat and Power Economic Dispatch via Improved Artiﬁcial Bee Colony Algorithm

: It is well accepted that combined heat and power (CHP) generation can increase the efﬁciency of power and heat generation at the same time. With the increasing penetration of CHPs, determination of economic dispatch of power and heat becomes more complex and challenging. The CHP economic dispatch (CHPED) problem is a challenging optimization problem due to non-linearity and non-convexity in both objective function and constraints. Hence, in this paper a novel meta-heuristic algorithm, namely improved artiﬁcial bee colony (IABC) algorithm is proposed to solve the CHPED problem. The valve-point effects, power losses as well as the feasible operation region of CHP units are taken into account in the proposed CHPED problem model and the optimal dispatch of power/heat outputs of CHP units is determined via the proposed IABC algorithm. The proposed algorithm is applied on three test systems, in which two of them are large-scale CHPED benchmarks. The obtained results and comprehensive comparison with available methods, demonstrate the superiority of the proposed algorithm for dealing with non-convex and constrained CHPED problem.


Motivation and Problem Statement
Heat is considered to be a byproduct of power generation in conventional power generation systems and when it is not fully used that results in lower efficiency. Co-generation systems or combined heat and power (CHP) generation systems use the heat from a power plant and send it around to interested consumers. Thus, co-generation plants can produce both heat and electricity with better energy efficiency and fuel usage [1]. In recent years, CHP systems have attracted more attention due to their higher efficiency (up to 85%), network loss reduction, and rapid return of investment [2,3] compared to conventional systems. The complexity of the economic dispatch problem will be increased by including the CHP systems. Hence, it is necessary to propose appropriate solution procedure to obtain optimal schedules for both heat and power.

Contributions
Many literature works are listed in the previous section that concentrated on the solution of the CHPED problem. Given that the CHPED problem is a non-convex and nonlinear problem of optimization, there is no mathematical or metaheuristic algorithm that can guarantee the optimal global solution to these problems. Because of high economic saving potential of better algorithms, this paper focuses on solution methodology of CHPED problem. In this study, a related problem in the literature [12,15,[20][21][22][23][24] is used to compare the results obtained with the methods previously applied to CHPED problem. In this paper a method based on improved artificial bee colony (IABC) algorithm is proposed to solve the CHPED problem. ABC is a heuristic optimization technique, which is based on the intelligent search behavior of honey bee swarm. It provides a population-based search procedure in which individuals (which called foods positions) are modified by the artificial bees, which their aim is to discover the places of food sources with the highest nectar amount [32]. The main contributions of this work can be summarized as follows: 1.
Proposing an improved version of artificial bee colony algorithm for dealing with non-convex optimization problems.

2.
Studying the effectiveness and performance of the proposed algorithm using normal and large-scale test systems and benchmark functions.

3.
Implementation of the proposed algorithm on CHPED problem with different sizes and characteristics. 4.
Compared with available methods in the literature, achieving feasible and better results for large-scale CHPED test systems.

Paper Organization
The rest of this paper is organized as follows. Section 2 provides the mathematical formulation of the CHPED problem considering valve-point effects, transmission losses and regional heat dispatch. Section 3 describes the proposed IABC algorithm. Section 5 gives the step by step procedure of the proposed IABC algorithm for solving the CHPED problem. Several case studies are presented in Section 6. Finally, conclusions are given in Section 7.

Chp Economic Dispatch Problem Formulation
The considered co-generation system in this study consists of power-only units, heat-only units and CHP units. The objective of the CHPED problem is to minimize total cost of serving the heat and power demands. The total cost can be stated as sum of the costs of generating heat and power as follows [20].
The cost functions of the different units can be expressed using the following quadratic functions.
The quadratic function approximation (2) is widely used in the literature for modelling the cost function of power-only units [17]. Usually, an absolute sinusoidal term is added to the quadratic cost function for modeling the valve-point effects [33] as follows.
Therefore, the CHPED problem becomes non-convex. The objective function (1), should be minimized subject to the following technical constraints [19]: where (6) models the power production and consumption balance. The power transmission system loss is calculated using B− matrix coefficients using (7). The heat production and demand balance is modeled using (8). The capacity limits of the power-only units and heat-only units are bounded using (9) and (10), respectively. The production limits of heat and power generation of CHP units are modeled using (11) and (12). It is observed from these equations that the upper and lower limits of power generation of CHP units are functions of produced heat (or vice versa). This heat-power dual dependency is presented using feasible operation region (FOR) for a specific CHP. The FOR of CHP units represents either a convex region or non-convex region as described in [34]. In the case of non-convex region, the CHPED problem becomes more complicated, due to non-convexity in both the objective function and the constraints. Typical FORs of CHP units is presented in Figures 1 and 2. As it can be observed, Figure 1 represents a convex region while Figure 2 represents a non-convex region [35].

Original Abc Algorithm
This algorithm is based on particle swarm intelligence that is inspired by the behavior of honey bees finding food. Bee colony algorithm first was proposed by Karaboga [32]. In this algorithm there are three categories of bees, i.e., employed, onlooker and scout bees. The population of employed bees and onlooker bees are equal (i.e., half of the colony).
Employed bees are responsible for exploiting the nectar sources and providing the waiting bees (onlooker bees) in the hive with information about the nature of the locations of the food source which they exploit. Onlooker bees wait in the hive and decide to exploit a food source based on knowledge exchanged by the bees they are working. Scouts either search the area randomly to find a new food source according to their internal motivation or based on potential external clues [36]. The process of finding food source in honey bee colony can be divided into three parts [32]:

1.
Employed bees discover food sources and determine the quality of nectar and share its location with others bees.

2.
Onlooker bees decide based on the quality of the food sources found by employed bees and follow the location of food sources of employed bees.

3.
If the food source of an employed bee is abandoned, it becomes scout bee and discover new food source randomly.
In ABC algorithm, each multi-dimensional particle (or food source) is shown as follows [32].
Thus, the i-th particle is shown as follows.
That bee employed is associated with a single place of food source. Therefore, the number of places of food supply is equal to the number of bees employed. In each iteration of ABC algorithm, employed bees discover food sources as follows [32].
In the above equation k is a random integer that it selected from the set {1, ..., SN}. After production of new solutions in each iteration, the fitness (nectar) function is calculated from the following expression [32].
where f (X i ) is the objective function value for X i to be minimized. After calculation of objective function fitness, if Fit(X new The onlooker bees select the employed bees location based on the fitness value of their corresponding food sources. For this purpose, the possibility of choosing the food source location is calculated as follows. As the nectar quantity of food sources (fitness of solutions) increases, so does the number of onlookers visiting them, which facilitates convergence to the optimal solution [36].
For each iteration a random real number is generated for each source within the range [0,1]. If the probability (Pi in (17) associated with this source is greater than this random number, the onlooker bee modifies the location of this source of food by using (15). After the food source is evaluated from (16), if the fitness value is improved, then the onlooker bee replaces the old food source location by the new one, otherwise it keeps the old location.
If after a certain number of iterations, employed bee's food source location does not improved, the food source location is abandoned and this location is replaced with a random new location by the scout bee from: where x max j and x min j are upper and lower bounds for j-th decision variable x j , respectively.

Improved Abc (Iabc) Algorithm
The ABC algorithm has been implemented successfully in various optimization problems such as in hydroelectric generation estimation [37] and parameter estimation of solar cells [38]. However, it still attracts the attention of many researchers to improve its performance. Most of these methods modify Equation (15). For example, Kraboga proposed a new search equation for employed bees as follows [36].
In the above equation MR is modification rate which is equal to 0.8, R i,j is a uniformly distributed random number in the interval [0, 1]. Also, in [39] Gao and Liu proposed a new search equation as follows.
x new i,j = x best j + rand(−1, 1) × (x old r 1 ,j − x old r 2 ,j ) (20) where r 1 and r 2 are mutually different random integers selected from the set {1, ..., SN}. x best j is the individual x j corresponds to the particle with the best fitness in the current population.
The results reported using the above modifications indicate that both of the above search rules are very effective approaches in the optimization problems solved by ABC algorithm [39].
In this paper, a hybrid search technique is proposed which combines the above search formulas as follows.
The main distinguishing features of the proposed IABC algorithm are as follows: 1.
As it is observed from Equation (21), as the difference between x r 1 ,j and x r 2 ,j decreases, the disturbance of position x i,j decreases. Therefore, the length of step is adaptively reduced by approaching to an optimal solution, and hence the algorithm converges to the optimal solution. 2.
It is observed from Equation (18) that the algorithm automatically jumps form local optimal or even non-optimal points, since scout bees are generated when no progress made in the search for a specified food source (or solution).

3.
Onlooker bees capability included in this algorithm enables comparison of the behavior of all food sources (or solutions) simultaneously. In other words, it is observed from Equation (17) that if a specified solution (or food source) i has a small P i , then it is a good solution, and hence it is not updated by onlooker bees. Otherwise, it is replaced with new position by onlooker bees.

Investigation of Iabc Algorithm on the Benchmark Functions
To investigate the performance of proposed IABC algorithm, two studies are conducted here as follows.

Study-I: Investigations on Six Benchmark Functions
In this study, six well-known benchmark functions which have different characteristics are examined. These functions which have been employed in [39,40] for evaluation of ABC algorithm and its variants, are given in Table 1.
The proposed IABC algorithm is applied to the above benchmark functions. Similar to the settings of GABC algorithm given in [40], the following settings are considered for the proposed IABC algorithm evaluation: Population size (or SN) is 80, maximum iterations number (Iter max ) is 5000 and the number of trial runs is 30.
For the above settings, the mean value and standard deviation of the results are presented in Table 2. The obtained results are compared with MABC [39], ABC [40] and GABC [40]. It can be observed from Table 2 that the proposed IABC algorithm converges to better results in comparison with ABC, MABC and GABC algorithms, in terms of the mean and standard deviation of the results.

Name Formula D (Problem Dimension) Search Space Global Minimum
Schaffer [39]

Study-Ii: Investigations on Large-Scale Benchmark Functions
To examine the capability of proposed IABC algorithm for solution of large-scale optimization problems, it is implemented on the 300 variables version of the benchmark functions f 4 -f 6 [41]. Table 3 compares the mean value obtained by the proposed IABC algorithm with GSO [41], GA [41], PSO [41], EP [41], ES [41], ABC [36], and MABC [39]. It is evidently observed from this table that IABC algorithm outperforms the above existing approaches, since the obtained solution by IABC is very close to the global optimal solution in all considered benchmark functions. Besides, in Table 4 the performance of IABC algorithm is compared with the basic ABC [36] and MABC [39] in terms of the best and worst obtained solutions for the above three benchmarks. It is observed from this table that the proposed IABC approach gives smaller values for both the best and worst solutions. especially for f 4 and f 5 these values are both zero, which means the algorithm always converges to the global optimal point.

Implementation of Iabc on the Chped Problem
To solve the CHPED problem by the IABC algorithm, the following steps are performed.

1.
Step-1: The first stage in IABC algorithm is initialization of the employed bees. Every food source location is a candidate solution of CHPED problem. The position of each food source (X i ) is a vector of all real power and heat outputs of the units as presented in the following.
i , X i ] The initial population of employed bees for power-only and heat-only units are determined from (26) and (27), respectively. Also, the population of employed bees for CHP units is determined from (28) and (29), respectively. 2.
Step-2: By setting Iter = 1 (where Iter is the iteration number of the algorithm), discover new food source locations by employed bees using Equation (21).

3.
Step-3: In this step the objective function value for the population of bees are calculated at the current iteration. Since the CHPED is a constrained optimization problem, it is converted to an unconstrained problem using penalty coefficient (λ). λ is assumed to be 10000 for all test systems studied in the following section. Hence, the objective function will be as follows.
Step-4: Fitness of i-th food source is calculated from Equation (16). If the new food source fitness is better than the old, then the old food source is replaced with the new location (obtained in Step-2, and B Scout,i = 0 (B Scout,i is a counter that determines limit value for converting i-th employed bee to scout bee), otherwise old location is preserved and B Scout,i = B Scout,i + 1.

5.
Step-5: At this step onlooker bees select food source of employed bees by using the roulette wheel criterion given in (17). Based on the value of P i for each food source, the onlooker bees modify the selected locations of employed bees by using (21), as follows: If P i > rand(0, 1), then the fitness of new food source is calculated from (16). If it is better than old fitness value, the old food source is replaced with new location and B Scout,i = 0, otherwise the old food location is kept and B Scout,i = B Scout,i + 1.

6.
Step-6: After the completion of the food source update process for employed and onlooker bees, if B Scout,i > Limit, then that food source is abandoned, and the employed bee is converted to a scout bee. The scout bee selects its new food source randomly in the space via (18).

7.
Step-7: Check the stopping criterion. If the algorithm converged, then go to Step-8, else Iter = Iter + 1 and go to Step-2 and repeat the above procedure. In this paper, the stopping criterion is reaching to the maximum number of iterations in each run. In other words, if Iter = Iter max , then the algorithm stopped. 8.
Step-8: Stop. To clarify the optimization process for energy engineers, the implemented method is presented in Figure 3.

Case Studies
In this section, the effectiveness and validity of the proposed method is evaluated by implementing it on three different test systems. The numerical study is performed using MATLAB 7.5 software on a PC with an Intel Core i7, 2.93 GHz CPU and 8 GB of RAM. The obtained results using the proposed IABC are compared with the reported results in the literature. The parameters used in the algorithm for different case studies are presented in Table 5.  [23]. Data of test system I is provided in Table 6. By investigating the published papers, it was found that there were three different loss matrix data for this system. Hence, we have solved the problem for this system using the available data in three cases.
The optimal dispatches of the units in this case are provided in Table 7. The obtained results using the proposed IABC algorithm are compared with the results of line-up competition algorithm (LCA) [42], teaching learning-based optimization (TLBO) [42], oppositional teaching learning-based optimization (TLBO) [42], conventional PSO (CPSO) [15] and time varying acceleration coefficients PSO (TVAC-PSO) [15] in Table 7. The distribution of the total cost for 50 independent runs is presented in Figure 4. It is inferred from this figure that in 33 runs the obtained value for cost is less than the mean value, which means the algorithm is capable of attaining a solution better that than the mean value in 66% of trails. The minimum, average and maximum values of the obtained costs for these runs are also provided in Table 7. The convergence characteristics of the proposed method for this case is depicted in Figure 5. It is observed that the proposed IABC algorithm converges quickly in early iterations and hence, the number of maximum runs can be decreased to save the solution time.

Case II
In this case, the coefficients of the network loss matrix are assumed to be as follows [44]. The optimal solution results using the proposed algorithm are presented in Table 8 and is compared with the real coded genetic algorithm (RCGA) [44] and Grey wolf optimization algorithm (GWO) [45]. It can be observed that the average of the obtained costs using the proposed algorithm for 50 independent runs is lower than the minimum cost obtained using RCGA.
Using the above data, the problem is solved using the proposed algorithm and the results are presented in Table 9. The obtained results are compared with particle swarm optimization (PSO) [22], real-coded genetic algorithm (RCGA) [22], evolutionary programming (EP) [22], artificial immune system (AIS) [21], bee colony optimization (BCO) [22] and differential evolution (DE) [23] in Table 9. It can be found that the proposed algorithm outperforms all of the previously proposed algorithms in the literature in less time.

Test System Ii (24-Unit System)
This system is one of the large benchmarks for the CHPED problem, which is proposed in [15]. Data of cost functions of the second test system is presented in Table 10. The valve-point effects are considered in this test system. This test system consists of 13 power only units, 6 CHP units, and 5 heat-only units. The total electricity demand is 2350 MW and the total heat demand is 1250 MWth. The power only units' data are based on the 13-unit challenging standard economic dispatch test system [46]. Data of the test system can be accessed from [15]. Table 11 shows the obtained results using the proposed algorithm for this system. The obtained results in this table are compared with the recent algorithms such as TLBO [43], OTLBO [43], CPSO [15], TVAC-PSO [15], arithmetic crossover harmony search (ACHS) [47], group search optimization method (GSO) [48], improved GSO (IGSO) [48] and Grey wolf optimization algorithm (GWO) [45]. As it can be observed from this table, the proposed IABC algorithm obtains a better solution compared to other reported algorithms in the literature. It should be mentioned that the minimum, average, and maximum cost of 50 independent runs are also presented in Table 11. Distribution of total costs in these trail runs for this test system is depicted in Figure 6. It is inferred from this figure that in 27 runs the obtained cost by the proposed IABC algorithm is less than the average value of 50 trials. This means that the algorithm is able to find a solution better that than the mean value in 54% of trails. The convergence characteristics of the proposed IABC algorithm for 24-unit test system is provided in Figure 7. As it can be observed from this figure the proposed algorithm is converged to the optimal solution in earlier iterations.   Figure 6. Distribution of total costs for 50 independent runs for 24-unit test system.

Test System Iii (48-Unit System)
This system consists of 48 units which includes 26 power-only, 12 CHP, and 10 heat-only units. The data of this system is given in [15]. The total electricity demand is 4700 MW and the total heat demand is 2500 MWth. Table 12 summarizes the optimal heat and power dispatches using the proposed method. The obtained results using the proposed algorithm are compared with the recent algorithms such as TLBO [43], OTLBO [43], CPSO [15], TVAC-PSO [15], group search optimization method (GSO) [48] and improved GSO (IGSO) [48].
It should be mentioned that the total costs presented in this table are directly quoted from the corresponding references. Some of the algorithms have obtained a lower cost, but with the expense of violating some of the technical constraints. Table 13 presents the feasibility analysis of the different algorithms. As it can be observed from this table, the solutions obtained by TLBO, OTLBO, GSO, and IGSO are not feasible. The CHP units which violate their corresponding feasible operation regions are indicate in Table 13 by × symbol. For example, for the TLBO method the CHP unit 38 violates its feasible region and the total power and heat mismatch is −69.9975. The minimum, mean and maximum obtained total costs are presented for 50 independent runs. Distribution of total costs for 50 independent runs for this test system is presented in Figure 8. It is observed from this figure that in 28 runs the obtained solution is better than the average solution, and if we ignore solution 34, the remaining solutions are also very close to the average solution. This shows good diversity of the obtained solutions by this algorithm. The convergence characteristics of the proposed algorithm for this system is depicted in Figure 9. As it can be observed from this figure the proposed algorithm is converged to the optimal solution in earlier iterations, similar to previous test cases, which confirms the capability of the algorithm in dealing with large scale CHPED problems. To calculate the annual cost-saving, the studied hourly load is considered to be the average load during a year. Therefore, in the 48-unit test system, as an example, the saving will be 6,082,861 ($/year) if the result of the proposed algorithm is compared with the result of the best available feasible algorithm (TVAC-PSO). This comparison shows that the proposed optimization algorithm reduces the system operation cost considerably.

Conclusions
A new approach based on the improved artificial bee colony (IABC) algorithm is proposed in this paper for an efficient solution of CHPED problem. Different characteristics and constraints such as valve-point effect, power losses, feasible operation region of CHP units, and capacity limits of units are taken into account in the formulation. The effectiveness of the proposed IABC algorithm is verified using standard benchmark functions and statistical analysis. It is found that the proposed algorithm can find better solutions in terms of the objective function value, convergence speed and the number of solutions with lower objective function than the mean value, compared with other versions of ABC algorithm and other heuristic algorithms. Three test cases with different size and characteristics are used to evaluate the efficiency of the propose algorithm. The obtained results using the proposed IABC algorithm are compared with the most recent proposed algorithms and it is observed that the IABC converges to a feasible solution with the lower total cost in a reasonable time in comparison with the previously reported algorithms. The numerical results substantiate that:

•
The obtained results by the proposed IABC algorithm has small diversity and in most cases the algorithm converges to optimal or near optimal solutions. In other words, the variance of the obtained solutions is small.

•
The algorithm converges in relatively small number of iterations. This means that the algorithm has a good converge speed which enables it to be used in large systems. • In test system I, the obtained value for the objective function is less than the average value in 66% of the trial runs. This is 54% for test system II and 56% for test system III, which means that the proposed algorithm is able to attain solutions lower than the mean value, in more than half of the trials.

•
The obtained results are also feasible which indicates that the algorithm has the capability of attaining solutions which are both optimal and feasible.
The better solution results, especially in large test systems, confirms the applicability of the proposed algorithm for dealing with the real world systems. From the application perspective, the proposed method results in an hourly saving of $205.14 per hour which means $1,797,033 saving in each year for small scale 7-unit test system. The hourly saving of $694.4 is obtained for a 48-unit case, which equals more than $6 million saving annually. As a future work, the proposed method could be extended to solve unit commitment problem considering CHP units.