A Novel Method for Economic Dispatch of Combined Heat and Power Generation

The paper proposes a modified Bat algorithm (MBA) for searching optimal solutions of Economic dispatch of combined heat and power generation (CHPGED) with the heat and power generation from three different types of units consisting of pure power generation units, pure heat generation units and cogeneration units. The CHPGED problem becomes complicated and big challenge to optimization tools since it considers both heat and power generation from cogeneration units. Thus, we apply MBA method with the purpose of enhancing high quality solution search ability as well as search speed of conventional Bat algorithm (BA). This proposed approach is established based on three modifications on BA. The first is the adaptive frequency adjustment, the second is the optimal range of updated velocity, and the third is the retained condition of a good solution with objective to ameliorate the search performance of traditional BA. The effectiveness of the proposed approach is evaluated by testing on 7, 24, and 48 units systems and IEEE 14-bus system and comparing results with BA together with other existing methods. As a result, it can conclude that the proposed MBA method is a favorable meta-heuristic algorithm for solving CHPGED problem


Introduction
Economic dispatch of combined heat and power generation (CHPGED) plays an important role in modern power systems. Thus, many researchers have constantly developed its problem formulation and applied different optimization tools for solving such problems. The fuel cost is comprised of heat generation and power generation in which heat from electricity generation process is retained and transmitted to heat loads such as industrial zones and manufacturers as well as for other aims [1]. In combined heat and power systems, pure power generation units and cogeneration units are in charge of providing electricity to electric loads while pure heat generation units and cogeneration units undertake heat generation and heat transmission to heat loads. In general, the operating enhancement of pure power plants is much simpler than that of cogeneration plants; that is because of the complicated output characteristics of a cogeneration plant. In practice, the energy generation characteristics of pure plants are a fuel cost function. For pure plants, this function is either only heat or power output. Meanwhile cogeneration plants comprise both heat and power outputs which shows that the operating range of cogeneration plants is more complicated. Therefore, the CHPGED issue artificial immune system (AIS) method. The objective is to obtain the optimal solution by using of few control parameters and a small number of iterations. However, the disadvantage of this method is the ability to suffer premature convergence when applying the aging operator to eliminate the bad antibodies.
This paper proposes a pertinent method for solving the CHPGED on the basis of the Bat algorithm (BA), considering three modifications that are the adaptive tuning frequency, optimal range of updated velocity, and retained condition of a good solution with the purpose to improve the searching performance of conventional BA. To verify the effectiveness of the proposed method, several simulation cases of the proposed method and conventional BA are analyzed based on four test case with a 7-unit system, 24-unit system, 48-unit system and an IEEE 14-bus system. In addition, transmission power losses, non-differential objective function and all constraints of transmission power networks have been taken into consideration. The main new contributions of this paper are summarized as follows: (i) Conventional BA is first applied for solving the CHPGED and (ii) Pointing out drawbacks of BA and suggesting modifications on conventional BA to develop modified Bat algorithm with better performance.
The remainder of this paper is organized as follows. Section 2 covers the methodologies for the CHPGED problem while Section 3 gives the proposed algorithm for CHPGED. Case study, simulation results, and discussion of the results are handled in Section 4 and Appendix A. Finally, conclusions are reported in Section 5.

Methodologies
Suppose that we have a complex power system that includes the pure power generation units N pp , pure heat generation units N ph , and cogeneration units N c . Here, the main objective of the CHPGED is to minimize total fuel cost of heat generation and electricity generation from available generation units. To reach the objective, the main duty is to determine heat output of cogeneration units and pure heat generation units, and power output of cogeneration units and pure power generation units. In addition, constraints of working limitations from all generation units, power balance and heat balance must be exactly satisfied. Especially, the relationship between heat and power outputs of a cogeneration units shown in Figure 1 is really complex.
Energies 2018, 11,3113 3 of 28 a small number of iterations. However, the disadvantage of this method is the ability to suffer premature convergence when applying the aging operator to eliminate the bad antibodies. In recent years, some most promising methods have proposed for the CHPGED, like the group search optimization [20], oppositional teaching learning based optimization [21], opposition-based group search optimization [22], and eight versions of PSO [23], grey wolf optimization [24], machine learning algorithm [25], and the cuckoo search algorithm [26]. It is observed that the performance of these methods seems to be better than all the methods given in [819] for most cases.
This paper proposes a pertinent method for solving the CHPGED on the basis of the Bat algorithm (BA), considering three modifications that are the adaptive tuning frequency, optimal range of updated velocity, and retained condition of a good solution with the purpose to improve the searching performance of conventional BA. To verify the effectiveness of the proposed method, several simulation cases of the proposed method and conventional BA are analyzed based on four test case with a 7-unit system, 24-unit system, 48-unit system and an IEEE 14-bus system. In addition, transmission power losses, non-differential objective function and all constraints of transmission power networks have been taken into consideration. The main new contributions of this paper are summarized as follows: (i) Conventional BA is first applied for solving the CHPGED and (ii) Pointing out drawbacks of BA and suggesting modifications on conventional BA to develop modified Bat algorithm with better performance.
The remainder of this paper is organized as follows. Section 2 covers the methodologies for the CHPGED problem while Section 3 gives the proposed algorithm for CHPGED. Case study, simulation results, and discussion of the results are handled in Section 4 and Appendix A. Finally, conclusions are reported in Section 5.

Methodologies
Suppose that we have a complex power system that includes the pure power generation units Npp, pure heat generation units Nph, and cogeneration units Nc. Here, the main objective of the CHPGED is to minimize total fuel cost of heat generation and electricity generation from available generation units. To reach the objective, the main duty is to determine heat output of cogeneration units and pure heat generation units, and power output of cogeneration units and pure power generation units. In addition, constraints of working limitations from all generation units, power balance and heat balance must be exactly satisfied. Especially, the relationship between heat and power outputs of a cogeneration units shown in Figure 1 is really complex.   Figure 1 is depicted by using given data of heat outputs and power outputs of a typical cogeneration unit in which the   Figure 1 is depicted by using given data of heat outputs and power outputs of a typical cogeneration unit in which the efficiency of the unit cannot be known as observing the figure but its performance can be evaluated by using a fuel cost function. The observation from Figure 1 in [27] and Figure 5 in [28] can indicate that efficiency of the cogeneration unit is complicated for calculation. In fact, the efficiency is directly influenced by thermal efficiency and heat efficiency, and it is also related to power output and heat output. Simulation results from [29] also indicate the same statement, since points of efficiency are not a mathematical function. In our work, we suppose that feasible operating zone and fuel cost function of cogeneration units are exactly given by taking data from previous studies. Our main task is to determine power and heat outputs so that the values must be inside the feasible operating zone and the fuel cost can be reduced as much as possible.

Cost Objective Function
For the objective of the CHPGED issue, we consider a minimum cost objective function of the heat and electricity power production under the following form: The cost functions of the pure power generator i for the case of neglecting and considering valve point loading effects are respectively mathematically formulated in Equations (2) and (3) [14]: The cost objective function of the kth pure heat generator can be calculated as follows [17]: The cost objective function of the jth cogeneration unit can be calculated as follows [17]: F cj (P cj , H cj ) = a cj + b cj P cj + c cj P 2 cj + k cj H cj + l cj H 2 cj + m cj H cj P cj , where P pi and P cj are the power output of the ith pure power generation unit and jth cogeneration unit, respectively; H cj and H hk the heat output of the j-th cogeneration unit and heat output of the kth pure heat generation unit, respectively; a pi , b pi , c pi , e pi and f pi are the pure electricity generation cost function coefficients of the ith pure power generation unit; a hk , b hk and c hk are the pure heat generation cost function coefficients of the kth pure heat unit; and a cj , b cj , c cj , k cj , l cj and m cj are the heat and electricity generation cost function coefficients of the jth cogeneration unit. Among the definitions, values of power outputs and heat outputs are determined by operators in power plants while coefficients in all cost functions are predetermined by experts in power plants by conducting experiments on the units.

Active Power balance
The active power generated by pure power and cogeneration units must be equal to the active power demand of power systems, and can be considered as follows: where P D is the power demand and P L is the loss power in the transmission line, and can be calculated as follows [14,17]:

Heat balance
The produced total heat due to the pure heat and cogeneration units must satisfy the heat demand and neglect heat loss, and can be considered as follows: In the constraint, heat losses are supposed to be zero because there have not been any studies about heat losses during process of transmitting heat to heat loads. In the work, we have also used that assumption for simplicity. However, if heat losses are a function of heat outputs similar to power loss function or a constant, heat balance constraint will be solved simply and successfully.

Power and heat unbalance
Each generation unit must operate within its upper and lower limit. The power of ith pure power, jth cogeneration, and the kth pure heat and the jth cogeneration units can be expressed, respectively, as follows: in which, the limit functions of power output of the j-th cogeneration unit can be calculated as follows: the limit functions of heat outputs of the j-th cogeneration unit can be calculated as follows:

Slack Units
In this paper, two slack variables including the slack pure power P p1 and the slack pure heat generation unit H h1 are considered to exactly meet the balance constraints of power and heat:

Pure Power Generation Unit
Suppose that the power outputs of (N pp −1) pure power generation units and N c cogeneration units are known. The power output of the first pure power generation unit can be calculated using Equation (6) as follows: where P L can be calculated by using Equation (7).

Pure Heat Generation Unit
The slack pure heat generation unit H h1 is considered using Equation (8) according to the known heat outputs of other pure heat and cogeneration units as:

Conventional Bat Algorithm
The Bat algorithm was developed in 2010 derived from the behavior of bat species during their food search [30]. The main rules of the BA method can be summarized as follows: (i) Each Bat d has a position X d and a velocity V d , and it generates a pulse with the frequency f d , the loudness A d , and the rate r d as well. (ii) The loudness A d generated by a Bat gets the largest value at the beginning of the search process, reduces gradually, and then reaches the lowest one as approaching to its prey.
Each Bat's position is the indices of a solution of optimization problems, which will be gradually enhanced while the number of iterations is increased, and the velocity is used to determine the following positions emblematizing the possibility of changing positions of a Bat. Increasing the number of iterations, normally, the searching process could generate a new position which would be closer to the food source than the previous positions. It means that the quality of a new solution tends to be better than that of the previous solution. To generate new positions, search space around each old solution, called global search and search space around the best solution in the population, called local search are carried out. Basic steps of BA can be summarized in detail as bellow.
Given a set of N B Bats, calling X G d and V G d are the current position and velocity of the dth Bat in the population, d = 1, 2, . . . , N B .
At the beginning of each iteration step, the frequency of each Bat can be randomly initialized as: and, the new velocity and position of each Bat are: where X* is the best position among N B Bats in a population. This solution is determined from the random positions of the initial population for the first iteration and chosen from existing population at the end of each iteration. For a local search, each Bat has a chance to generate a new position which is around the best location, can be described as follows: where ε is a random number selected in range [0, 1], and A G is the average loudness of all Bats. In fact, the loudness and pulse emission rate of each Bat are not fixed values, but they should be rapidly changed as the Bat moves around its prey. It means that once the Bat has discovered a prey, its loudness tends to decreases meanwhile the pulse emission rate increases. As a result, the updated loudness and rate are determined as: in which the coefficients α and γ are arbitrary values in range [0, 1], and kept constant during the search process. In the studies [30][31][32], the two parameters were set to the same value of 0.9 for simplicity, but they should be dependent on the specific applications by experiment. However, analyzing Equation (24), we can see that the pulse emission rate does not increase if the γ is set to 0.9, and the probability for local search is not higher during the search process. The BA's pseudo code was described in detail in [30]. The BA has been more advanced than other algorithms in terms of accuracy and efficiency for many test functions [30] as well as successfully applying for simple cases of optimal operation of power systems [33,34]. However, to deal with more complex problems, some modifications have been suggested in the variants of the BA as mentioned in [35].

Modified Algorithm
In order to improve the searching performance of the BA for solving the CHPGED, in this paper, we propose three new modifications, including adaptive tuning frequency, the optimal range of updated velocity, and new criteria of competitive solutions. The three modifications can be summarized as follows:

Adaptive Frequency Adjustment Corresponding to the Number of Iterations
In the BA, at each iteration, the frequency is randomly produced within its operating frequency range as defined in Equation (19). As a result, an updated velocity generated arbitrarily as defined in Equation (20) could be over-speeded and make the Bat fly away from the optimal point. This problem is more serious at higher iterations where the Bats are closely approached to the prey. In the optimization point of view, at higher iterations when a Bat is nearly reach to the prey, this Bat's flying speed must be slowed down and the movement should be just around the current position (neighboring areas of the prey's location). Thus, by experiment, we suggest that at each iteration, the updated frequency should be adjusted as: In Equation (25), our intention is to reduce the frequency gradually according to the increase of iteration. As the iterations increase, the frequency decreases. On the other hand, the new velocity obtained from Equation (20) and new position calculated from Equation (21) are also narrowed as the iteration increases because the velocity is proportional to the frequency and the position is proportional to the velocity. It is clear that the search space is exploited in large ranges at the beginning of the search process but it contracted as the iterations increase to the maximum value.

Set the Optimal Range of the Updated Velocity
The velocity is one main factor used to update new position for each Bat. If the velocity is high, the Bats' positions, i.e., solutions could move the feasibility area outside whereas the small velocity could limit the discovery areas. For the BA, there is no information discussing the optimal range of velocity, while this term has been constantly improved in PSO [23]. In [23], authors have mentioned that the optimal range of the velocity can be fixed from 10% to 15% of the difference between the and X max d and X min d are the maximum and minimum positions in the searching area, respectively.

Define New Criteria of Competitive Solutions
In the BA, a new solution is retained only if both conditions are satisfied, that is its fitness function is smaller than that of the old solution (of the same Bat) and its assigned random value is equal to or lower than the current loudness [30,31,35]. The selection technique of BA can expressed by the following Algorithm 1.

Algorithm 1. Selection technique of BA
In the algorithm, if the two conditions are true, the loudness A G d becomes much smaller than the previous value because α is fixed at 0.9 and the possibility that better solution is kept to be low. Thus, we suggest that the loudness should be set to 1 at the beginning and afterwards decreased by a very small step size as shown in Equation (28): The equation can increase the value of A G d compared to that of Equation (23), giving an opportunity for the use of a new solution with a better fitness function and replacement of old solutions with worse fitness functions.
The proposed modified bat algorithm (MBA) has been applied for a complex optimal power system operation problem where the combined heat and power outputs have been closely related together. In this problem, the power system was supplied by many kinds of the normal thermal generators and cogeneration sources.

Application
In CHPGED problems, a power system is supplied by N pp pure power generation units, N ph pure heat generation units and N c cogeneration units. Similar to other meta-heuristic algorithms, the proposed MBA needs a set of Bats N B where each Bat has its own velocity and position. The position of the Bat dth, called X d (d = 1 ÷ N B ), is a solution of the problem containing power outputs of (N pp −1) pure power units and heat outputs of (N ph-1 ) pure heat generation units together with heat and power outputs of N c cogeneration units. The maximum iterations are fixed at beginning and the position of Bats (solutions) N B will be updated time by time after each iteration. The search process-based MBA will guide the Bats moving toward the prey (the optimal solution) quickly and stably as shown in the following steps. The flowchart for solving the CHPGED is presented in Figure 2.
value equal to 1, almost all new solutions are nearby the current global best solution and the final optimal solution may possibly fall in a local optimal value. Inversely, few solutions could be newly generated that can be nearby the global best solution if the rd is set to a very low value. Each specific optimization problem should be existed a suitable range of the pulse rates depending on its characteristics. Therefore, in order to discover the suitable, the pulse rates of the MBA for the CHPGED problem will be carried out by varying within range [0, 1] with a step of 0.1. The new solutions need to be checked whether they satisfy their limits, which are mentioned in Equations (9)(12), by using the following models: Step 1: Initialization of the whole population The initialization is implemented as follows: From the initialized values of N B Bats, the fitness function corresponding to the position of each Bat is calculated as follows: where P p1,d nd H h1,d are the power and heat output of the slack pure power and pure heat generation units of the d-th solution, obtained from Equations (17) and (18), respectively; and K p and K h are penalty coefficients for the power and heat outputs of those slack generation units, respectively. The limits of the slack generation units P lim p1 and H lim h1 in Equation (34) are determined as follows: Note that P max p1 and P min p1 are the maximum and minimum power generations of the slack power generation unit, which must be similar to the limits of other pure power generation units. H max h1 and H min h1 are the maximum and minimum heat productions of the slack heat generation unit to be similar to the the limits of other pure heat generation units.
Among all initialized members, the Bat's position with the lowest value of fitness functions is assigned as the global best position (solution), X * .
Step 2: Update New Velocity and Position for Each Bat In this proposed method, each Bat has a position and a velocity where the position represents a solution and the velocity represents an updated step size for a new position. Therefore, the value of the velocity has an impact on the obtained solution while the value of current frequency also plays a very important role to generate an effective velocity. The adaptive frequency f G d is caculated from Equation (25) and the each Bat velocity is updated by applying Equation (20).
Furthermore, the new velocity needs to be checked by the optimal range of the velocity as mentioned in Equations (26) and (27). The constraint for this velocity can be described by: and a new position of each Bat related to the new velocity is determined from Equation (21).
Step 3: Search a New Solution As described in the conventional BA [30,35], there is a possibility that a current X d can be newly produced around the current global best solution by applying Equation (22) if its random value is less than the pulse rate r d , within range [0,1]. This r d is set to a specific value at the beginning of the search process and increased during the following iterations. Obviously, if the r d is set to at a high value equal to 1, almost all new solutions are nearby the current global best solution and the final optimal solution may possibly fall in a local optimal value. Inversely, few solutions could be newly generated that can be nearby the global best solution if the r d is set to a very low value. Each specific optimization problem should be existed a suitable range of the pulse rates depending on its characteristics. Therefore, in order to discover the suitable, the pulse rates of the MBA for the CHPGED problem will be carried out by varying within range [0, 1] with a step of 0.1.
The new solutions need to be checked whether they satisfy their limits, which are mentioned in Equations (9)−(12), by using the following models: The power of the slack power generation unit and the heat of the slack heat generation unit can be calculated by using Equations (17) and (18), and can be checked by applying Equations (35) and (36), respectively.
Step 4: Find New Solution through the New Criteria of Competitive Solutions The current solutions are evaluated through the calculation of the fitness functions, and then the selection of between the current and the initialized solutions) is decided according to the last modification in Section 3.2. We update the new loudness A d according to the current number of iterations G as shown in Equation (28).

Case Study and Discussion
In the section, four study cases with different numbers of units, constraints and objective functions are carried out for investigating the performance of MBA. The description of the four cases in detail is as follows: Case 1: We consider a power system with seven generation units under two conditions of transmission power losses and the non-convex fuel cost function. Case 2: We consider a power system with 24 generation units under two conditions of the neglecting power losses and the non-convex fuel cost function. Case 3: We consider a power system with 48 generation units under two conditions of the neglecting power losses and the non-convex fuel cost function. Case 4: We consider a IEEE 14-bus system with three pure power generation units, two cogeneration units and one pure heat generation unit in which non-convex fuel cost function is considered for slack pure power generation unit.
The power and heat load demands and the information about the number of Bats, number of iterations are listed in Table 1. The control parameters for the BA and the proposed method are chosen at each study case with the purpose of obtaining acceptable results for both the BA and MBA methods when compared to other existing methods. For the CHPGED problem, in this paper some quantitative criteria such as minimum cost ($/h), maximum cost ($/h), average cost ($/h), standard deviation ($/h), and successful rate (number of successful runs/total runs) are considered to compare the performance in which the comparison of costs is used to evaluate the quality of searching process whilst the successful rate reflects the ability of handling the complex problem. The initial pulse rate is considered in range [0.1, 1] with step 0.1.
However, how to verify the effectiveness and robust of the BAs with other methods mentioned in the literature review may not be an easy mission because, in fact, different studies may select the different computer configurations and different values of the population size and maximum number of iterations. In order to adapt to the mentioned differences, the maximum number of fitness evaluations (FE max ) and the scaled computational time (SCT) introduced in [36] and [37] have been used, respectively, re-applied as two relative criteria for fair comparisons when different studies have run their proposed methods in different conditions as different population sizes, different iteration number, and different computer configuration. In fact, for meta-heuristics, FE max is also a comparison criterion to analyze the effectiveness among considered methods [36]. For this study, we consider the FE max = N P × G max [37] to compare the performance of all methods. This is a value that directly concerned with the quality of obtained solution and execution time. On the other hand, the SCT in Equation (42) has been applied to convert the reported CPU time from a different processor to a common value associated with the same processor [37]. A 2.4 GHz processor used in the study. In addition, to show how many times of a proposed method is faster than other ones, SCT should be expressed in per unit (p.u.) as shown in Equation (43) 4.1. Case 1 The power system with seven generation units, including four pure power generation units, two cogeneration units, and one pure heat generation unit, is considered to be the first case. This power system supplies an electricity power demand of 600 MW and heat demand of 150 MWth [14]. The data of this system are given in Tables A1-A3 and in real power loss matrix coefficients in the Appendix A. The obtained results with different pulse rate values through BA and MBA are listed in Tables 2 and 3, respectively. It is pointed out that the best minimum cost from BA is obtained at pulse rate of 0.6 whilst that from MBA is at pulse rate of 0.9; however, the better results from the two methods tend to be in the range from 0.6 to 0.9 while the range from 0.1 to 0.5 results in worse average cost, maximum cost, and standard deviation. When calculating the difference between two values, all the best values of the minimum, average, and maximum cost and standard deviation from the proposed method are less than that from BA by $199.1353, $1080.5437, $431.4505, and $118.6701, respectively. Clearly, the proposed method can reflect higher efficiency than BA in terms of optimal solution and the stable quality of solution. In addition, this case study also investigates the success rate (SR) between the two methods and reports the values at the final column. From the reported values, it can be concluded that the proposed method can handle all constraints more effectively than BA since SR from the proposed method is around 96% but that is only from 70% to 80.6% obtained by BA. The whole search process by using BA and the proposed method at the best run with the lowest fuel cost plotted in Figure 3 from the fourth iteration onwards indicates that the proposed method can find better optimal solution than BA at each iteration. In addition, the same fitness values from the 80th iteration onwards also report that BA has trapped a local optimal solution and its search could not be out of the local zone. On the contrary, the proposed method could tackle the drawback since its fitness value decreased gradually after a low number of iterations. Figure 4 plots the fluctuation of 100 obtained fitness function values by BA at r 0 = 0.6 and the proposed method at r 0 = 0.9. This figure shows that the proposed method has obtained higher peak than BA while the stabilization of 100 runs of the proposed method is much better than those from BA. Most of runs, MBA has yielded better optimal solutions than BA because most points of MBA are below points of BA. On the other hand, there are many solutions having approximate fitness functions with the best solution of MBA while points that are close to the best point of BA are few.  Table 4 illustrates the minimum cost, execution time, and associated information with the comparison of computational time (CPU time) using MBA and other methods. Observing this table, the fuel cost comparison indicates that MBA finds a better optimal solution than most methods excluding two methods, that are GCPSO [23] and GWO [25] with slightly less costs of $33.54 and $66.09 than MBA, respectively. Despite the disadvantage, MBA is still a the strongest method since its CPU time is the fastest among the compared methods, especially 1.01 s for GCPSO [23] and 5.26 s for GWO [25] while that is only 0.1 second for MBA. Moreover, as indicated in the FE max and SCT (s) columns, MBA used 4000 fitness evaluations and spent a runtime of 0.1 second, which are lower than all other methods like BCO [14] and AIS [17], that used 5,000 fitness evaluations and take 6.45 s and  6.62 s, respectively. When compared to GCPSO [23] and GWO [25], FE max and SCT(s) are 20,000 and 0.76 s when using GCPSO, respectively. Whilst there was no population size reported when using GWO, it has taken a quite long time of 5.04 s to search the solution. For the SCTs (PU), MBA is much faster than other methods by 50.41, 64.45 and 80.90 times corresponding to GWO [25], BCO [14] and RCGA [14], respectively. Based on three comparisons, we conclude that the proposed method is a very efficient method for solving the CHPGED problem for a system with seven units. The optimal solution obtained for this study case is given in Table A4 in the Appendix A.   Table 4 illustrates the minimum cost, execution time, and associated information with the comparison of computational time (CPU time) using MBA and other methods. Observing this table, the fuel cost comparison indicates that MBA finds a better optimal solution than most methods excluding two methods, that are GCPSO [23] and GWO [25] with slightly less costs of $33.54 and $66.09 than MBA, respectively. Despite the disadvantage, MBA is still a the strongest method since its CPU time is the fastest among the compared methods, especially 1.01 s for GCPSO [23] and 5.26 s    Table 4 illustrates the minimum cost, execution time, and associated information with the comparison of computational time (CPU time) using MBA and other methods. Observing this table, the fuel cost comparison indicates that MBA finds a better optimal solution than most methods excluding two methods, that are GCPSO [23] and GWO [25] with slightly less costs of $33.54 and $66.09 than MBA, respectively. Despite the disadvantage, MBA is still a the strongest method since its CPU time is the fastest among the compared methods, especially 1.01 s for GCPSO [23] and 5.26 s for GWO [25] while that is only 0.1 second for MBA. Moreover, as indicated in the FEmax and SCT (s)

Case 2
The power system with 24 generation units consisting of 13 pure power generation units, five pure heat generation units, and six cogeneration units, is considered to be the second study case. This power system supplies an electricity power demand of the electrical load demand 2350 MW and a heat load demand of 1250 MWth. The data of this system has been taken from [18] and are given in Tables A5-A7 in the Appendix A. The two control parameters N B and G max are set at 40 and 200 for BA, respectively; and set at 20 and 150 for MBA, respectively. For the study case, SR when applying BA is 74.6% poorer than that applying MBA which is 92.8%. For the minimum fuel costs, when applying BA and MBA, they are $58,198.82 and $57,851.9133, respectively. The fitness function value of the best solution from the 2nd iteration onwards, when using BA and MBA, is plotted in Figure 5. Observing this figure, from the beginning to the 90th iteration, BA and MBA have found approximately good solutions and it is difficult to state which one is better than the other, however, from the 90-th iteration to the end of the search process, BA has failed to improve its search ability because the fitness function values seemed to decrease inconsiderably and it have trapped at a local optimum. Unlikely BA, MBA has still improved its solutions gradually and its fitness function values have been much lower than those of BA. This means that MBA could tackle the local optimum better than BA. Figure 6 plots the 100 fitness function values of BA and MBA. Observing this figure, MBA is superior to BA because most points of MBA are below those of BA. Besides, the differences of fitness function values of BA and MBA are too high, leading to a space above with approximately black color while the space below with approximately blue color. Other information that can be seen is that MBA can search many solutions of similar quality for the best solution. Table 5 lists a comparison of the results yielded by many reference methods. The minimum cost yielded by MBA is $ 57,851.9133, which is less than most methods eliminating GSO [20], OBGSO [22] and GWO [24]. The comparison implies that the proposed method is ranked as the fourth best optimal solution among the compared methods. In addition, to indicate the exact amount of money that MBA can save compared to others, we have calculated the saving cost by determining the difference between the reported cost from others and that from MBA as shown in the last column of Table 5. Observing this figure, MBA obtained a higher cost than GSO [20], OBGSO [22] and GWO [24] [21] and PSO [18], respectively. On the other hand, to apply fair comparisons among the different methods, the values of FE max and SCT (in seconds and PU) have also been calculated and listed in Table 5. It is shown in this table that MBA has employed the second lowest FE max with 3,000 fitness evaluations, which is only higher than that of GWO [24] (1350 fitness evaluations) and significantly less than 1,000,000 of PSO as well as TVAC-PSO [18]. As for the criterion of convergence speed, MBA is outstanding over all other methods, e.g., 75 times faster than GWO [24] and 635.29 times as fast as PSO [18]. In this criterion, BA is also faster than all methods except MBA with the SCT (in PU) of 1.43 (i.e., 1.43 times slower than MBA). Based on the three advantages over other methods including the best optimal solution, the second lowest FE max and the fastest convergence speed, the proposed MBA is absolutely superior to all methods and it is a very efficient method for solving the CHPGED problem in the system with 24 units. The optimal solution obtained by the proposed MBA for case 2 is given in Table A8 in the Appendix A.
that MBA can save compared to others, we have calculated the saving cost by determining the difference between the reported cost from others and that from MBA as shown in the last column of Table 5. Observing this figure, MBA obtained a higher cost than GSO [20], OBGSO [22] and GWO [24] by $5.0733, $8.3933 and $22.6633, respectively. Compared to other methods, MBA can save from $4.3567 to $1884.3467 related to OTLBO [21] and PSO [18], respectively. On the other hand, to apply fair comparisons among the different methods, the values of FEmax and SCT (in seconds and PU) have also been calculated and listed in Table 5. It is shown in this table that MBA has employed the second lowest FEmax with 3,000 fitness evaluations, which is only higher than that of GWO [24] (1350 fitness evaluations) and significantly less than 1,000,000 of PSO as well as TVAC-PSO [18]. As for the criterion of convergence speed, MBA is outstanding over all other methods, e.g., 75 times faster than GWO [24] and 635.29 times as fast as PSO [18]. In this criterion, BA is also faster than all methods except MBA with the SCT (in PU) of 1.43 (i.e., 1.43 times slower than MBA). Based on the three advantages over other methods including the best optimal solution, the second lowest FEmax and the fastest convergence speed, the proposed MBA is absolutely superior to all methods and it is a very efficient method for solving the CHPGED problem in the system with 24 units. The optimal solution obtained by the proposed MBA for case 2 is given in Table A8 in the Appendix.

Case 3
In order to investigate more the capability of the proposed method, in this case study, we consider a large-scale power system with 48 units, consisting of 26 pure power generation units, 10 pure heat generation units, and 12 cogeneration units. This power system supplies a 2500 MWth heat demand and a 4700 MW power demand [18]. The data of the system can be obtained by duplicating the data of System 2 in Case 2. The control parameters N B and G max are fixed at 60 and 200 for BA, respectively; and they are fixed at 30 and 200 for MBA, respectively. Figure 7 plots the fitness function value of the best solution between the 3rd iteration and the 200th iteration obtained when applying BA and MBA, respectively. In general, the best solutions from the two methods could dramatically increase their quality during the first half of the process, afterwards they showed a slight decrease and finally reached the lowest value at the end of the search. However, as shown in Figure 7, the rate of improvement of the fitness values in MBA has been higher than that in BA during the search process, and especially during the last twenty iterations, where MBA's fitness values have been still decreased significantly but BA's fitness values have almost kept constant. Similar to the previous Cases 1 and 2, BA has still become trapped in a local optimum while MBA could tackle the drawback for the case. Figure 8 can confirm the superiority of MBA over BA as in Cases 1 and 2 as most points of MBA are below those from BA. However, there is a big difference in this figure compared to Figures 4 and 6 in that MBA has more stable search than BA for the complicated case with 48 units and non-convex fuel cost function because the points of BA are in the the upper space while the space below belongs to the points of MBA. Table 6 shows the results in terms of minimum cost, computational time, FE max , and SCT (in second and PU) of BA, MBA and other methods. In addition, the saving costs have also been included in the table for further comparisons similar to the two previous cases. As seen from this table, the modifications on MBA are successful for searching for an optimal solution since MBA can obtain better minimum fuel cost, lower FE max , and faster SCT compared to its original BA. In fact, these values from MBA are $115966.0232, 6000 fitness evaluations, and 0.08 s, respectively; which are much lower than $117627.2256, 12,000 fitness evaluations and 0.1 s from BA, respectively. A further investigation of the efficiency of MBA can be executed by focusing on the saving cost and the SCT in PU. Clearly, MBA is much more effective and robust than all methods because it can save huge amounts of money and converges to an optimum significantly faster. For instance, it can save more money and is faster than TVAC-PSO [18]  These compared methods are either well-known original algorithms or improved variants of original algorithms, which have been famous and superior to others in previously published articles. From the three advantages of the MBA over all methods such as the lowest minimum cost, lowest maximum number of fitness evaluations and fastest convergence speed, it can be concluded that MBA is the most efficient method for solving the system up to 48 units considering valve effects on pure power generation units. The optimal solution obtained by the proposed MBA for the case is given in Table A9 in Appendix A.
Energies 2018, 11, 3113 18 of 28 Cases 1 and 2 as most points of MBA are below those from BA. However, there is a big difference in this figure compared to Figure 4 and Figure 6 in that MBA has more stable search than BA for the complicated case with 48 units and non-convex fuel cost function because the points of BA are in the the upper space while the space below belongs to the points of MBA.   Table 6 shows the results in terms of minimum cost, computational time, FEmax, and SCT (in second and PU) of BA, MBA and other methods. In addition, the saving costs have also been included in the table for further comparisons similar to the two previous cases. As seen from this table, the modifications on MBA are successful for searching for an optimal solution since MBA can obtain better minimum fuel cost, lower FEmax, and faster SCT compared to its original BA. In fact,  Cases 1 and 2 as most points of MBA are below those from BA. However, there is a big difference in this figure compared to Figure 4 and Figure 6 in that MBA has more stable search than BA for the complicated case with 48 units and non-convex fuel cost function because the points of BA are in the the upper space while the space below belongs to the points of MBA.   Table 6 shows the results in terms of minimum cost, computational time, FEmax, and SCT (in second and PU) of BA, MBA and other methods. In addition, the saving costs have also been included in the table for further comparisons similar to the two previous cases. As seen from this table, the modifications on MBA are successful for searching for an optimal solution since MBA can obtain better minimum fuel cost, lower FEmax, and faster SCT compared to its original BA. In fact, these values from MBA are $115966.0232, 6000 fitness evaluations, and 0.08 s, respectively; which are

Case 4
In order to evaluate the MBA capability of dealing with the ccomplicated constraints of the CHPGED problem, we have tested both MBA and BA in study Case 4 with the use of an IEEE 14-bus system. The original IEEE 14-bus system shown in Figure 9 covers five generators located at five different buses 1, 2, 3, 6 and 8, 15 load buses, 20 transmission branches, and three transformers and one compensator. For application of the combined heat and power generation units, we have replaced the pure power generation units at bus 2 and bus 3 with two different cogeneration units, which have been taken from [18]. In addition, one pure heat generation unit taken from [18] has also been installed for supplying heat to load in addition to the heat from the two cogeneration units. The power demand of 15 loads is 259 MW while the heat demand of load is supposed to be 400 MWth. The whole data of the IEEE 14-bus system can be obtained by referring to [40] while the entire data of three pure power generation units, two cogeneration units and one pure heat generation unit are given in Tables A10-A12 of Appendix A.   40. Although MBA has used a smaller population size it always get a higher success rate than BA. The highest SR of MBA is approximately 95% at r0 = 0.6, 0.8 and 0.9, but BA only gets 77% at r0 = 0.7. The search process of the best run of BA at r0 = 0.9 and MBA at r0 = 0.8 are depicted in Figure 10. Besides, the best cost of 100 successful runs of BA at r0 = 0.9 and MBA at r0 = 0.8 are shown in Figure 11. Figure 10 indicates that MBA could find better solutions than BA after the 23rd iteration, and the solution of MBA at iteration 120 had a lower cost than that of BA at the final iteration. It is clearly visible that MBA was faster than BA. On the other hand, MBA also shows its better stability since most runs of MBA had a lower cost than those of BA as shown in Figure 11. Thus, we can conclude that MBA can deal with the complicated constraints of CHPGED problems considering a transmission power network and its performance is significantly improved compared to BA. The optimal solution of MBA can be obtained by referring to Table A13 in Appendix A.   For implementing BA and MBA, we have set N B and G max to 20 and 150. As a result, BA could not deal with the constraints successfully whilst MBA could achieve a success rate of lower than 50% but the optimal solutions were of low quality. Thus, we have increased the two parameters to 20 and 200 for MBA but 40 and 200 for BA. Table 7 shows minimum cost, average cost and success rate obtained by BA and MBA for 100 successful runs. As we may observe, the minimum values of BA and MBA can show that MBA could find better optimal solutions than BA corresponding to each value of r 0 . The savings cost also indicates that MBA could improve the results from $81.612 to $353.80956. The best minimum of BA is $12,614.0837 while that of MBA is $12,532.4616. Clearly, MBA could improve results over BA by 0.65%. However, the result can be improved more significantly if we decrease N B of BA to 20 or we increase N B of MBA to 40. Although MBA has used a smaller population size it always get a higher success rate than BA. The highest SR of MBA is approximately 95% at r 0 = 0.6, 0.8 and 0.9, but BA only gets 77% at r 0 = 0.7. The search process of the best run of BA at r 0 = 0.9 and MBA at r 0 = 0.8 are depicted in Figure 10. Besides, the best cost of 100 successful runs of BA at r 0 = 0.9 and MBA at r 0 = 0.8 are shown in Figure 11. Figure 10 indicates that MBA could find better solutions than BA after the 23rd iteration, and the solution of MBA at iteration 120 had a lower cost than that of BA at the final iteration. It is clearly visible that MBA was faster than BA. On the other hand, MBA also shows its better stability since most runs of MBA had a lower cost than those of BA as shown in Figure 11. Thus, we can conclude that MBA can deal with the complicated constraints of CHPGED problems considering a transmission power network and its performance is significantly improved compared to BA. The optimal solution of MBA can be obtained by referring to Table A13 in Appendix A.

Conclusions
This paper has proposed a novel modified Bat algorithm (MBA) developed based on the conventional bat algorithm (BA) with the purpose of finding high quality solutions for CHPGED problems. Challenges to MBA arise from non-convex objective functions and the complicated constraints of CHPGED problems. Complex non-convex objective functions are caused by the consideration of valve point loading effects on pure power generation unit while the complex constraints are from the feasible working zone of cogeneration units and of all components in Figure 11. The best cost values obtained by BA and MBA over 100 runs for Case 4.

Conclusions
This paper has proposed a novel modified Bat algorithm (MBA) developed based on the conventional bat algorithm (BA) with the purpose of finding high quality solutions for CHPGED problems. Challenges to MBA arise from non-convex objective functions and the complicated constraints of CHPGED problems. Complex non-convex objective functions are caused by the consideration of valve point loading effects on pure power generation unit while the complex constraints are from the feasible working zone of cogeneration units and of all components in transmission power network of IEEE 14-bus system. This proposed method has been constructed by performing three modifications. The first is the adaptive frequency adjustment, the second is the optimal range of new velocities, and the third is the retained condition of a good solution. These improvements have been implemented for overcoming the shortcomings of BA. As a result, the proposed method has been more effective than BA for applications to four different CHPGED problem systems in terms of higher quality solution, higher success rate, faster search speed and more stable search ability. The comparisons with other existing methods have shown that the proposed method had the third best solution, the lowest fitness evaluations and the fastest convergence speed for the Case 1 with seven generation units, but the best method regarding its optimal solutions for the Cases 2 and 3, which were systems up to 24 and 48 generation units, respectively. MBA has continued to show its potential search by dealing with IEEE 14 bus system and obtaining outstanding cost as well as success rate as compared to BA. Consequently, we can reach the conclusion that the proposed method is a very effective method for solving combined heat and power systems with different challenges from pure generation units, cogeneration units and transmission power networks. The maximum number of fitness evaluation The cost function of the kth pure heat generation unit F pi (P pi ) The cost function of the ith pure power generation unit G, G max The current and maximum iterations, respectively H cj The heat output of the jth cogeneration unit H D The heat demand H hk The heat output of the kth pure heat generation unit K h The penalty coefficients for the heat outputs of slack generation unit K p The penalty coefficients for the power output of slack generation unit N B Number of Bats N c Number of the cogeneration unit N ph Number of the pure heat generation unit N pp Number of the pure power generation unit P cj The power output of the jth cogeneration unit P D The power demand P L The loss power in the transmission line P pi The power output of the ith pure power generation unit SCT The scaled computational time SR The successful rate ε The random number rand d A random number for solution d Appendix A