Improved group search optimization method for solving CHPED in large scale power systems

Energy Conversion and Management 80 (2014) 446–456

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Improved group search optimization method for solving CHPED in large scale power systems Mehrdad Tarafdar Hagh, Saeed Teimourzadeh ⇑, Manijeh Alipour, Parinaz Aliasghary Faculty of Electrical and Computer Engineering, University of Tabriz, Tabriz, Iran

a r t i c l e

i n f o

Article history: Received 23 November 2013 Accepted 24 January 2014

Keywords: Combined heat and power (CHP) Economic dispatch Group search optimization Non-convex optimization Prohibited operating zones

a b s t r a c t This paper presents an improved heuristic algorithm for solving combined heat and power economic dispatch (CHPED) problem in large scale power systems, by employing improved group search optimization (IGSO) method. The basic deficiency of the original GSO algorithm is the fact that, it gives a near optimal solution rather than an optimal one in a limited runtime period. In this paper, some modifications have been applied to the original GSO in order to improve its searching ability. In this work, a sinusoidal term is employed to consider the valve-point effects and Kron’s formula is used in order to consider the transmission losses. Moreover, prohibited operating zones of power-only units are considered. The effectiveness and accuracy of the proposed method for solving the non-linear and non-convex problems is validated by carrying out simulation studies on sample benchmark test cases and three large scale power systems. Numerical results indicate that proposed method has the better performance in comparison with original GSO and some recently published papers. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The conversion of primary fossil fuels, even by the most modern combined cycle plants, can only achieve efficiencies between 50% and 60%. However, the combined heat and power (CHP) units can achieve fuel conversion efficiencies of the order of 90% [1]. In order to obtain the optimal utilization of CHP units, economic dispatch must be applied to achieve the optimal combination of power and heat output of the CHP units while satisfying heat and power demand of whole power system. The CHPED problem is a complicated problem due to mutual dependencies of heat and power capacity of CHP units. The problem even may be more complicated if transmission losses and the effect of valve-points are taken into account. Solving such a non-linear and non-convex optimization problem requires powerful techniques. In order to handle the economic dispatch problem, two types of approaches are broadly used, i.e., Lagrange multiplier based methods and stochastic search methods [2,3]. Lagrange methods usually implement the derivative operations during optimization, i.e., Maclaurin series-based Lagrangian method [4]. The basic deficiency of the Maclaurin series-based Lagrangian method is the fact that, in some cases the cost function may have no derivative or it ⇑ Corresponding author. Tel./fax: +98 935 508 4689. E-mail addresses: [email protected] (M.T. Hagh), teimourzadeh91@ms. tabrizu.ac.ir (S. Teimourzadeh), [email protected] (M. Alipour), [email protected] (P. Aliasghary). http://dx.doi.org/10.1016/j.enconman.2014.01.051 0196-8904/Ó 2014 Elsevier Ltd. All rights reserved.

may be difficult to derive the derivative function. Therefore, deterministic approaches are not proper if the information for defining the search direction at each step is unavailable. Recently, stochastic search techniques are becoming more popular than ever as an alternative to deal with complex and difficult engineering optimization problems [5]. Stochastic weight trade-off particle swarm optimization (PSO) algorithm in [6], hybrid CPSOSQP in [7] and time-varying acceleration coefficients IPSO in [8] are some examples of PSO based methods for solving the economic dispatch problem. In addition, in [9], incremental bee colony algorithm considering prohibited operation zones effect is proposed. A modified harmony search based approach is proposed in [10] in order to solve environmental/economic dispatch problem for realworld power systems. Moreover, the CHPED problem is solved using Lagrange multiplier methods and stochastic search methods in literature. In [11], CHPED problem is solved using Lagrangian relaxation with surrogate sub gradient multiplier updates. In [12], an improved genetic algorithm with multiplier updating (IGA-MU) is utilized to solve CHPED problem. Linear cost function for heat-only and power-only units is assumed in [12]. Optimal solution of CHPED problem using harmony search (HS) algorithm is proposed in [1,13]. However, in [1,13], valve-point effects and cubic cost function for power-only units are not considered. In [14], ant colony search algorithm (ACSA) along with modified ACSA are recommended, in order to solve premature convergence problem of CHPED. In [15], multi objective CHPED problem is solved

447

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

considering wind power generation using PSO algorithm. The authors of [15] have used quadratic cost functions for the CHP units while the objective is minimizing cost and emissions. In [16], self-adaptive real-coded genetic algorithm (SARGA) is implemented and penalty constraint handling strategy is used to handle the equality and inequality constraints. In [14–16], the system loss and valve-point effect are not taken into consideration. Differential evolution (DE), bee colony optimization algorithm and time varying acceleration PSO (TVAC- PSO) algorithm are proposed to solve the CHPED problem considering valve-point effect and network losses in [17–19], respectively. A review of research works related to economic dispatch of CHP units is presented in [20]. Moreover, in [21,22], a benders decomposition approach and a direct solution are utilized to handle the CHPED problem, respectively. Recently, a new, easy-to-implement, reasonably fast and robust evolutionary algorithm has been introduced known as group search optimization algorithm (GSO). The GSO is inspired by group-living, a phenomenon typical of the animal kingdom [23,24]. Most of above reviewed algorithms [6–22], are stochastic schemes and have uncertain performance. Therefore, the reported results could be improved. According to the results of [23,25], GSO based algorithms have better performance in comparison with above mentioned heuristics methods for solving economic power dispatch problem. However, original GSO often converges to local optima and also its convergence speed is almost low. In this work, an improved GSO (IGSO) algorithm is proposed and implemented for some benchmark test cases and CHPED problem. The proposed modifications help the producer member of GSO group to access the optimal solution. Therefore, searching ability and convergence speed of the conventional GSO are improved. At each iteration, appropriate candidate positions are generated for the producer member. These points help the producer member to find the optimal solution more prompt than that of the conventional GSO. Moreover, the paper performs more accurate model of the problem by considering valve-point effects, network losses and prohibited operating zones effect. In addition, in this paper, it is assumed that, the heat does not transfer in long distance transmission and it distributes locally. Comprehensive simulation studies are carried out in order to validate the accuracy of the proposed method. The efficiency of the proposed method is evaluated by using sample benchmark test cases and some large scale power systems. Based on the simulation studies, the proposed method illustrates its efficiency in large scale power systems. The proposed method has better performance in comparison with conventional and some recently published papers. Moreover, convergence speed of MGSO is improved in comparison with the original GSO. The paper is organized as follows: In Section 2, the CHPED problem formulation considering valve-point effect and network losses is discussed. Section 3 presents the proposed IGSO algorithm. In Section 4 the step by step implementation of IGSO for CHPED problem is presented. In Section 5, in addition to some benchmark test cases, the proposed method is implemented on three large scale power systems. Finally Section 6 concludes the findings and contributions of the paper.

2. Problem formulation

min

Np Nh Nc X X X C i ðPpi Þ þ C j ðPcj ; Hcj Þ þ C k ðHhk Þ ð$=hÞ i¼1

j¼1

where Ci, Cj and Ck are the production cost of the conventional power, cogeneration and thermal units, respectively. Np, Nc and Nh are the number of the conventional power, cogeneration and thermal units, respectively. The H and P are the heat and power output of units, respectively. The total cost of the conventional thermal power unit can be expressed as: 2

C i ðPpi Þ ¼ ai ðPpi Þ þ bi Ppi þ ci ð$=hÞ

The objective of the CHPED problem is to minimize the total heat and power production cost. The problem considers both power and heat load demands constraints. There are three types of units including power only, combined heat and power and heat only units. Mathematically, the problem is formulated as follows:

ð2Þ

where ai, bi and ci represent the cost coefficients of ith conventional power unit. The cost function of the cogeneration units can be expressed as a convex function of electrical and thermal power outputs as follows: 2

2

C j ðPcj ; Hcj Þ ¼ aj ðPcj Þ þ bj Pcj þ cj þ dj ðHcj Þ þ ej Hcj þ fj Pcj Hcj ð$=hÞ

ð3Þ

where aj, bj, cj, dj, ej and fj represent the cost coefficients of jth cogeneration unit. The cost function of heat-only units is expressed as: 2

C k ðHhk Þ ¼ ak ðHhk Þ þ bk Hhk þ ck ð$=hÞ

ð4Þ

where ak, bk and ck represent the cost coefficients of the kth heatonly unit. Minimization of the cost function is subjected to following constraints: (a) Power balance constraint Np Nc X X Ppi þ P cj ¼ Pd þ P loss i¼1

ð5Þ

j¼1

where Pd is the total power demand of the system. The active transmission losses of system should be considered in the economic dispatch problem. One of the most common methods to express transmission losses as function of generator power is to derive Kron loss formula which is known as B-matrix method. The transmission loss formula using B-matrix coefficients is [26]:

PLoss ¼

Np X Np X

Ppi Bim Ppm þ

i¼1 m¼1

Np X Nc X

Ppi Bij Pcj þ

Nc X Nc X

i¼1 j¼1

Pcj Bin Pcn

ð6Þ

i¼1 n¼1

(b) Heat balance constraint Nh Nc X X Hpi þ Hcj ¼ Hd i¼1

ð7Þ

j¼1

where Hd is the heat demand of the system. (c) Generation limits of units

Ppi min 6 Ppi 6 Ppi max i ¼ 1; . . . ; Np

ð8Þ

Pcj min ðHcj Þ 6 Pcj 6 Pcj max ðHcj Þ j ¼ 1; . . . ; Nc

ð9Þ

Hcj min ðPcj Þ 6 Hcj 6 Hcj max ðPcj Þ j ¼ 1; . . . ; Nc

ð10Þ

Hhk min 6 Hhk 6 Hhk max

ð11Þ

Ppi min

2.1. Original formulation

ð1Þ

i¼1

i ¼ 1; . . . ; N h

Ppi max

where and are the minimum and maximum values of the ith unit in MW, respectively. Pcj min ðHcj Þ and Pcj max ðHcj Þ are power generation limits of jth CHP unit which are functions of generated heat ðHcj Þ . Hcj min ðPcj Þ and Hcj max ðP cj Þ are heat generation limits of CHP unit j which are functions of generated power ðPcj Þ. It should be mentioned that the heat and power outputs of cogeneration units are non-separable and dependent to each other. Hhk min and Hhk max are minimum and maximum values for kth heat-only unit, respectively.

448

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

2.2. Considering valve-point effects

[24]. Afterwards, the modification procedure of the proposed algorithm will be presented.

The valve-point effects not only results ripples in some heatrate curves but also makes the objective function to be non-convex with multiple minima. For accurate modeling of valve point loading effects, a rectified sinusoidal function [26] is added to the cost function. The ith unit cost function considering valve-point effects is given as:

       2   C i Ppi ¼ ai Ppi þ bi Ppi þ ci þ ki sin qi ppi min  Ppi  ð$=hÞ

ð12Þ

where ki and qi are cost coefficients for modeling valve-point effects. Considering valve-point effects make the CHPED problem non-convex. Therefore, conventional gradient based optimization algorithms could not be efficient. Fig. 1, depicts the cost-power characteristic with and without valve-point effects for ai = 0.00028, bi = 8.1, ci = 550, ki ¼ 300, qi = 0.035, Pmin ¼ 0 and i Pmax ¼ 680. It could be observed from Fig. 1, the cost-power characi teristic considering valve-point effects is higher than that of conventional quadratic function. 2.3. Considering prohibited operating zone effects In practical power systems, due to the mechanical constraints, thermal units have prohibited operating zones. Therefore, the cost-power characteristics of the thermal units become discontinuous. The feasible operating zones for thermal units could be formulated as follows [27]:

8 > Pmin 6 Pi 6 Pli;1 > > < i l U Pi 2 Pi;j1 6 Pi 6 P i;j j ¼ 2; 3; . . . ; nPOZi > > > : PU 6 P 6 Pmax i;nPOZ

i

i

ð13Þ

i;1

where nPOZi is number of prohibited operating zones of the ith thermal unit. P li;j1 and P Ui;j are upper and lower limits of the jth prohibited operating zone, respectively. Fig. 2 depicts cost-power characteristic of thermal units considering prohibited operating zones for ai = 0.00028, bi = 8.1, ci = 550, P min ¼ 0, P max ¼ 680 and i i prohibited operating zones of [90–100] and [400–415]. 3. Optimization algorithm This section, firstly presents a brief overview of group search optimization algorithm. More details of GSO could be found in

3.1. Basics of group search optimization algorithm The GSO algorithm is a novel optimization algorithm which is based on animal searching behavior and their group-living theory. The algorithm is mainly based on the producer–scrounger (PS) model [25,28]. In the proposed framework, the group members search for either ‘‘finding’’ or ‘‘joining’’ opportunities. In other words, the animal scanning mechanisms are employed metaphorically for designing an optimum searching strategy in order to solve continuous optimization problems [24]. The population of the GSO algorithm is called a group and each individual in the population is called a member. At each iteration, a member is defined by its position and head angle. In n-dimensional search space, the ith member of the GSO at the kth searching iteration has a current position X ki 2 Rn and a head angle uki ¼ ðuki1 ; . . . ; ukiðn1Þ Þ 2 Rn1 . Where, R is the set of real numbers and uin is polar angle of ith member relative to the nth dimension. In addition, each member has search direction where, the search direction of the ith member, Dki ðuki Þ could be calculated from its head angle vector, uki , via a polar to Cartesian coordinate transformation as follows [29]: k

k

Dki ðuki Þ ¼ ðdi1 ; . . . ; din Þ 2 Rn 8 n1 Y > > > cosðukiq Þ j¼1 > > > > < q¼1 k n1 Y dij ¼ > sinðukiðj1Þ Þ cosðukiq Þ j ¼ 2; . . . ; n  1 > > > > q¼j > > : sinðuk j¼n iðn1Þ Þ

ð14Þ

The GSO group includes three types of members: producer, scroungers and rangers. The producer finds the opportunities, scroungers join to the opportunities and rangers perform random walk motions. At each iteration, a group member which has the best fitness value, is chosen as the producer. Producer by its vision ability scans the search space for the better states. Vision ability is the ability of testing some points around the producer current position. The producer scans three points around its position in certain distances and head angles. At the kth iteration, the producer behaves as follows:

Fig. 1. Cost-power characteristic with and without valve-point effects.

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

449

Fig. 2. Cost-power characteristic considering prohibited operating zones effects.

(1) The producer scans at zero degree and tests three points toward its position using Eqs. (15)–(17).

X z ¼ X kp þ r1 lmax Dkp ðuk Þ

ð15Þ

X r ¼ X kp þ r1 lmax Dkp ðuk þ r 2 hmax =2Þ

ð16Þ

X l ¼ X kp þ r1 lmax Dkp ðuk  r 2 hmax =2Þ

ð17Þ

where X p is position of the producer, r 1 2 R1 is a normally distributed random number with mean 0 and standard deviation 1 and r 2 2 Rn1 is a uniformly distributed random number in the range of (0, 1), lmax is maximum pursuit distance and hmax is maximum pursuit angle. (2) The producer will then find the best point. If the best point has a better value in comparison with its current position, producer will fly to that point. If not, it will stay in its current position and turn its head using Eq. (18):

ukþ1 ¼ uk þ r2 amax

ð18Þ

1

where amax e R is the maximum turning angle. (3) If the producer cannot find a better area after a iterations, it will turn its head back to zero degree as follows:

ukþa ¼ uk

ð19Þ

where a is a constant. At each iteration, some of group members are selected as scroungers. The scroungers will keep searching for opportunities to join the positions found by the producer member. At the kth iteration, the ith scrounger can be modeled as a randomly walks toward the producer using Eq. (20):

X kþ1 ¼ X ki þ r 3  ðX kp  X ki Þ i

ð20Þ

3.2. Improved group search optimization algorithm According to the results of [23], the conventional GSO algorithm finds a near optimal solution rather than an optimal one in a limited runtime period. In this paper, the problem is solved by applying some modifications. At each iteration, the producer member has the best fitness value. Hence, candidate points in the vicinity of the producer position perhaps are better approach to the global optimal point. Candidate points are new positions for the producer. Therefore, scroungers could follow the candidate points and searching process would be done utilizing multi-producers. However, the candidate points are just generated when, after certain iterations, the fitness value does not change. In other words, at some iterations, the multi-producer searching process is derived either to help the main producer to depart the local minima, or prompt it to accelerate the searching procedure. In this paper, if the fitness value does not change after 10 iterations, multi-producer searching process should be done. Based on the comprehensive simulation studies, it seems that 10 is the best value. If it sets to greater number, intentional delay is added to the searching process. On the other hand, if it sets to a lower number, searching process of ranger members is disrupted. The candidate points are modeled as follows:

   iter10  If X iter 6e p  Xp

ð23Þ

X rptest ¼ X iter p þ r4

   X max  X p iter max iter iter X max

ð24Þ

X lptest ¼ X iter p  r4

   X max  X p iter max iter iter X max

ð25Þ

where X ki is position of ith scrounger at kth iteration and r3 2 Rn is a uniform random sequence in the range of (0, 1). Operator ‘‘°’’ is the Hadamard product, which calculates the entry wise product of the two vectors. At each iteration, some of group members are selected as rangers. Rangers are dispersed from their positions and randomly walk at search space. At the kth iteration, a ranger generates a random head angle ui using Eq. (17), and then it chooses a random distance using Eq. (21) and move to the new point using Eq. (22).

where r4 e R1 is a uniformly distributed random number in the range of (0, 1). Xmax is the maximum value of that variable, iter max is maximum number of iterations and e is a threshold value. Fig. 3, depicts the candidate points and scroungers path. By using Eqs. (23)–(25), several candidate points are generated in vicinity of the main producer. Afterwards, the candidate points are tested by the main producer using Eqs. (15)–(17). In each step, X z ; X r and Xl are appointed as producers rank one to three, respectively. The producers are modeled as follows:

li ¼ ar1 lmax

ð21Þ

X p1 ¼ X z

ð26Þ

X kþ1 ¼ X ki þ li Dki ðukþ1 Þ i

ð22Þ

X p2 ¼ X r

ð27Þ

450

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

 Step 8. Perform scrounging related to the ranked producers: In this step, the scrounger members are divided into four groups and ranked producers are randomly appointed to each group of scroungers. Afterwards, scroungers perform scrounging using Eqs. (29)–(32). Flowchart of the proposed IGSO is depicted in Fig. 4. 4. Implementation of IGSO for CHPED problem The steps of solving the non-convex CHPED problem using IGSO algorithm are as follows:

Fig. 3. Candidate points and scroungers path.

X p3 ¼ X l

ð28Þ

where Xp1, Xp2 and Xp3 are producers rank one to three, respectively. Afterwards, the scrounger members are divided into four groups and producers are randomly appointed to each group of scroungers. Members of each group follow their corresponding producer. Therefore, scrounging procedure at kth iteration is modified as follows:

  X kþ1 ¼ X ki þ r 5  X kp  X ki i

 Step 1. Generate and evaluate initial members: Real power and heat outputs are the decision variables in the CHPED problem. Therefore, arrays of output power and heat of all units considered as members of MGSO group. Each member has sub members. At kth iteration, the ith member position presents as follow:

h i X ki ¼ P pjk ; P Cmk ; HCmk ; Hhnk

ð33Þ

Ppjk ; PCmk ; HCmk and Hhnk are sub members of member X ki . Sub members are initialized as follows:

ð29Þ

ð34Þ

X kþ1 i0

  ¼ X ki0 þ r6  X kp1  X ki0

  Ppjk ¼ Ppjkmin þ r 1 P pjkmax  P pjkmin

ð30Þ

ð35Þ

X kþ1 i00

  ¼ X ki00 þ r7  X kp2  X ki00

  Pcjk ¼ Pcjkmin þ r 2 P cjkmax  P cjkmin

ð31Þ

ð36Þ

X kþ1 i000

  ¼ X ki000 þ r 8  X kp3  X ki000

  Hcjk ¼ Hcjkmin þ r 3 Hcjkmax  Hcjkmin

ð32Þ

kþ1 X kþ1 ; X kþ1 i i0 ; X i00

X kþ1 i000

where and are members of various groups of scroungers, respectively. The r 5 ; r 6 ; r 7 and r8 are random numbers in the range of (0, 1). 3.3. Summary of proposed algorithm Based on the discussions of Sections 3.1 and 3.2, important steps of proposed IGSO are as follows:  Step 1. Generating and evaluating initial members: Similar to other heuristic algorithms, in the proposed method, the initial population should be generated regarding to the upper and lower limits of variables. Afterwards, generated population should be evaluated by objective function.  Step 2. Choosing main producer and perform producing: In this step, the main producer should be chosen. Afterwards, the main producer performs producing using Eqs. (15)–(17).  Step 3. Choosing Scrounger members and perform scrounging: Some of group members are chosen as scroungers and they perform scrounging using Eq. (20).  Step 4. Choosing Ranger members and perform ranging: The rest of group members are chosen as rangers and they perform ranging using Eqs. (21), (22).  Step 5. Decision point: In this step, based on Eq. (23), the future path is identified. It is identified either to continue single producer path or switch to multi-producer path.  Step 6. Generating Candidate points: In this step, candidate points should be generated using Eqs. (24), (25).  Step 7. Set Producers rank one to three: After generating candidate points, the candidate points are tested by the main producer using Eqs. (15)–(17) and producers rank one to three are appointed using Eqs. (26)–(28).

Fig. 4. The flowchart of the IGSO.

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

Hhjk ¼ Hhjkmin

  þ r 4 Hhjkmax  Hhjkmin

451

ð37Þ

where j is the number of unit and r1, r2, r3, r4 are uniformly distributed random numbers in the range of (0, 1). For CHP units, the generated heat is a function of its generated power and the upper and lower bands of generated heat varies by its generated power. Therefore, at each iteration, the upper and lower bands of generated heat of CHP unit should be determined corresponding to its generated power. Fig. 5 depicts power–heat feasible region for a sample CHP unit.  Step 2. Fitness evaluation: Eq. (38) is used to calculate the fitness function, which is the total cost of system. The fitness function should be minimized while satisfying all constraints. At each iteration, inequality constraints should be checked before calculating the fitness value and if they are not in feasible band, they have to be fixed on their limits.

Fitness function ¼

Np Nh Nc X X X C i ðPpi Þ þ C j ðPcj ; Hcj Þ þ C k ðHhi Þ i¼1

j¼1

ð38Þ

k¼1

 Step 3. Producing: Fitness function should be calculated for all members of IGSO group. A group member, which has the best fitness value, would be the producer. Producer performs producing using Eqs. (15)–(17).  Step 4. Scrounging: 40% of IGSO group members are chosen as scroungers and they perform scrounging using Eq. (20).  Step 5. Ranging: 60% of MGSO group are chosen as rangers and they perform ranging using Eq. (22).  Step 6. Improvement procedure: As mentioned, for improving the search ability of this algorithm after certain iterations, if the producer cannot find a better state in comparison with its current position, discussed modifications in Section 3.2 are implemented.  Step7. Check stopping criterion: The terminating criterion is selected to be the maximum number of iterations. The algorithm will be terminated if the maximum iteration number reaches; otherwise it continues from Step3. Based on the above discussion, the flow chart of the proposed method for CHPED problem is depicted in Fig. 6.

Fig. 6. Flow chart of IGSO for CHPED.

5. Case studies

Fig. 5. Power-heat feasible region for a sample CHP unit.

In this section, the proposed method is tested on some test benchmark functions and two large scale test power systems that are selected from other literatures. The proposed method solves CHPED problem considering capacity limits and heat–power mutual dependency constraints. In addition, valve-point effect

452

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

and system losses are taken into account by employing a sinusoidal term and Kron’s loss formula, respectively in large cases. Numerical results indicate that the proposed method has the best performance, in comparison with compared methods. IGSO based methodology is developed by Matlab 7.6 in 2.5 GHz, i5, personal computer. The maximum number of iterations is set to 1000 for all benchmarks and 300 for CHPED problems. In addition, population size for benchmarks is set to 48 [24] and 25 for CHPED problem.

Table 2 Comparison of different algorithms for benchmark functions. Function

GA [24]

PSO [24]

GSO [24]

Proposed IGSO

f1(x) f2(x) f3(x) f4(x) f5(x) f6(x) f7(x)

3388.5516 0.6509 1.0038 9749.9145 3.697 12566.0977 1.0038

37.3582 20.7863 0.2323 0.00011979 0.146 9659.6993 0.2323

49.8359 1.0179 0.0030792 5.7829 0.0016 12569.4882 0.030792

33.8201 0.00345 0.00011 0.05011 0.00358 12569.4431 0.00021

5.1. Benchmarks In this section, some benchmark functions are studied in order to validate the performance of the proposed IGSO. Table 1 presents the definition of benchmark functions [24]. The proposed method is applied to the benchmark functions of Table 1 and the mean results are reported in Table 2. In addition, the results are compared with genetic algorithm (GA), particle swarm optimization algorithm (PSO) and conventional GSO. The default parameter settings are utilized for GA and PSO [24]. From Table 2, it could be observed that, the proposed method has either equal or better performance than that of mentioned methods.

Table 3 Comparison of simulation result for Case I.

5.2. Case I Test system I, suggested by [19], includes 24 units, where units 1–13 are power only units, units 14–19 are cogeneration units and units 20–24 are heat-only units. Data of the conventional thermal generators is based on the 13-unit standard ED test system which has a lot of local minima and is one of the challenging ED test cases [30,31]. The parameters of cost function of cogeneration units are available in [19]. The heat and power demand are 1250 (MWth) and 2350 (MW), respectively. The comparison results of proposed method with CPSO, TVAC-PSO and conventional GSO are reported in Table 3. For this case CPSO [19], TVACPSO [19], conventional GSO and IGSO have computed 59736.26 ($), 58122.74 ($), 58225.74 ($) and 58049.01 ($) as minimum total cost, respectively. The difference between CPSO algorithm and TVACPSO, conventional GSO and IGSO algorithms is mainly due to power outputs of units 2 and 3. In addition, power outputs of unit 12 and unit 13 are involved in the mentioned difference. According to Table 3, IGSO has achieved to the optimal point in 35.54 s, while CPSO and TVACPSO calculated the total cost in 53.36 and 52.25 s, respectively. The 32% difference in computation time shows the better performance of proposed algorithm. Convergence characteristic of conventional and modified GSO algorithms are depicted in Fig. 7. 5.3. Case II This test case includes 26 power-only units, 10 heat-only units and 12 CHP units. Data of this case is obtained by duplicating data

a

Outputa

CPSO [19]

TVACPSO [19]

GSO

Proposed IGSO

P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19

680 0 0 180 180 180 180 180 180 50.5304 50.5304 55 55 117.4854 45.9281 117.4854 45.9281 10.0013 42.1109

538.5587 224.4608 224.4608 109.8666 109.8666 109.8666 109.8666 109.8666 109.8666 77.521 77.521 120 120 88.3514 40.5611 88.3514 40.5611 10.0245 40.4288

627.7455 76.2285 299.5794 159.4386 61.2378 60 157.1503 107.2654 110.1816 113.9894 79.7755 91.1668 115.6511 84.3133 40 81.1796 40 10 35.0970

628.152 299.4778 154.5535 60.846 103.8538 110.0552 159.0773 109.8258 159.992 41.103 77.7055 94.9768 55.7143 83.9536 40 85.7133 40 10 35

H14 H15 H16 H17 H18 H19 H20 H21 H22 H23 H24

125.2754 80.1175 125.2754 80.1174 40.0005 23.2322 415.9815 60 60 120 120

108.9256 75.4844 108.9256 75.484 40.0104 22.4676 458.702 60 60 120 120

106.6588 74.9980 104.9002 74.9980 40 19.7385 469.3368 60 60 119.6511 119.7176

106.4569 74.998 107.4073 74.998 40 20 466.2575 60 60 120 119.8823

Mean ($) Maximum ($) Minimum ($) Time(s)

59853.48 60076.69 59736.26 53.36

58198.311 58359.552 58122.746 52.25

58295.9243 58318.8792 58225.7450 35.54

58156.5192 58219.1413 58049.0197 35.54

Output powers are in (MW) and heats are in (MWth).

of Case I [19]. Where units 1–13 and 14–26 are power-only units and their characteristics are the same as units 1–13 in case I, units 27–32 and 33–38 are CHP units with characteristics the same as units 14–19 in case I, units 39–43 and 44–48 are heat-only units and their characteristics are the same as units 19–24 of case I.

Table 1 Benchmark test functions [24]. Benchmark function f1 ðxÞ ¼

Pn1 i¼1

Pn

ð100ðxiþ1 

2 x2i ÞÞ

þ ðxi  1Þ

2

2

2 i¼1 ðxi

 10 cosð2pxi Þ þ 10Þ   P 0:5   P  f3 ðxÞ ¼ 20 exp 0:2 1n ni¼1 x2i  exp 1n ni¼1 cosð2pxi Þ þ 20 þ e f2 ðxÞ ¼

f4 ðxÞ ¼ f5 ðxÞ ¼

Pn

i¼1 ð

Pn

i¼1 ð½xi

Pn

2

j¼1 xj Þ

2

þ 0:5Þ pffiffiffiffiffiffiffi jxi j

i¼1 xi sin Pn 1 i¼1 ðxi  4000

f6 ðxÞ ¼  f7 ðxÞ ¼

Pi

100Þ2 

Qn

i¼1

xi 100 pffi i

n

Search space

Global minima

30

[30, 30]n

0

30

[5.12, 5.12]n

0

30

[32, 32]n

0

30

[100, 100]n

0

n

0

30

[100, 100]

30

[500, 500]n

12569.5

30

[600, 600]n

0

453

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

Fig. 7. The cost convergence curve of the proposed IGSO and conventional GSO of Case I.

The total power and thermal demands of this case is 4700 (MW) and 2500 (MWth). Table 4 presents the optimal heat and power dispatches using the proposed method. According to Table 4, the total cost obtained by CPSO [19], TVACPSO [19] and conventional GSO algorithms are 119708.9 ($), 117824.89 ($) and 113131.8 ($) respectively, and total calculated cost of IGSO method is 112320.41. This difference in total cost calculation is mostly due to the difference in produced power of units 2, 3 and 14, and generated heat of unit 36 and unit 44. There is 5% difference in total cost and 22% difference in computation time compared with TVACPSO method in this large scale power system. Convergence characteristic of conventional and modified GSO algorithms are depicted in Fig. 8.

5.4. Case III In this case, simulation studies are carried out considering both valve-points and prohibited operating zones effects. System data of this case is the same as case I and extracted from [19]. In addition, in this paper, some prohibited operating zones for power-only units of this system are proposed. Table 5 presents prohibited operating zones for power-only units of case I. According to Table 3 and Table 5, power-only units number 1, 3, 11 and 12 are operating in the prohibited operating zone. By deriving CHEPED problem considering prohibited operating zones, operation problem of mentioned units could be totally solved. Table 6 presents the optimal heat and power dispatches using the

Table 4 Comparison of simulation results for Case II. Output

CPSO [19]

TVAC-PSO[19]

GSO

IGSO

Output

CPSO [19]

TVAC-PSO [19]

GSO

IGSO

P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P20 P21 P22 P23 P24 P25 P26 P27 P28 P29 P30

359.0392 74.5831 74.5831 139.3803 139.3803 139.3803 139.3803 139.3803 139.3803 74.7998 74.7998 74.7998 74.7998 679.881 148.6585 148.6585 139.0809 139.0809 139.0809 139.0809 139.0809 139.0809 74.7998 74.7998 112.1993 112.1993 92.8423 98.7199 92.8423 98.7199

538.5587 75.134 75.134 140.6146 140.6146 140.6146 140.6146 140.6146 140.6146 112.1998 112.1998 74.7999 74.7999 269.2794 299.1993 299.1993 140.3973 140.3973 140.3973 140.3973 140.3973 140.3973 74.7998 74.7998 112.1997 112.1997 86.9119 56.1027 86.9119 56.1027

627.5814 302.5046 225.3696 178.6488 178.2134 159.8844 161.4173 108.776 109.0234 115.1364 114.2308 107.2839 93.0811 0 223.7257 356.9056 109.2667 160.4169 109.6482 160.0005 174.5336 118.6394 40.063 41.2253 55 92.0406 81.3512 40 81.0383 40

629.4952 151.9991 299.2996 159.2254 173.6004 93.4383 160.773 159.351 161.4184 115.2927 112.8994 97.5394 55 0 299.268 225.4102 162.4605 160.9664 164.0177 168.4149 159.5402 110.8099 40.6399 114.3701 92.3275 55 82.1821 40 81.089 40.4281

P31 P32 P33 P34 P35 P36 P37 P38 H27 H28 H29 H30 H31 H32 H33 H34 H35 H36 H37 H38 H39 H40 H41 H42 H43 H44 H45 H46 H47 H48

10.0002 56.7153 109.1877 65.6006 109.1877 65.6006 10.6158 60.5994 111.4458 125.6898 111.4458 125.6898 40.0001 29.8706 120.6188 97.0997 120.6188 97.0997 40.2639 31.6361 357.9456 59.9916 59.9916 120 120 370.6214 59.9999 59.9999 119.9856 119.9856

10.0031 35 95.4799 54.9235 95.4799 54.9235 23.4981 54.0882 108.1177 88.9006 108.1177 88.9006 40.0013 20 112.926 87.8827 112.926 87.8827 45.7849 28.6765 433.9113 60 60 120 120 415.9741 60 60 119.9989 119.9989

10 35.2736 82.878 40 81 40.3336 10.5087 35 104.9965 74.998 104.8209 74.998 40.001 19.2636 105.5564 74.998 104.8032 332.3293 39.5514 20 486.4858 60 60 118.7549 113.2371 212.5981 59.5362 59.9138 113.9272 119.2305

10.6913 35.0696 81 40.1014 81.0922 40.1056 10 35.6838 103.5903 74.998 104.2548 75.3686 40.0999 19.2943 104.8032 75.0858 104.8511 75.086 40 20.3111 428.0157 59.5061 59.9205 114.8048 117.9877 535.65 60 60 107.7179 118.6434

Cost ($) Time (s)

CPSO[19]

TVAC-PSO [19]

GSO

Proposed IGSO

119708.9 93.32

117824.896 89.63

113131.8301 70.65

112320.4159 70.65

454

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

Fig. 8. The cost convergence curve of the proposed IGSO and conventional GSO of Case II.

Table 7, the results of three small test cases based on authors investigations, are presented in rows 2–4 [32]. Considering restrictions in the space of this paper, detailed results of small test cases are not presented in this paper. In addition, simulation results for large scale power systems are presented in rows 4–5. According to Table 7, by increasing the size of test systems, difference between total cost and run time of proposed IGSO, increases in comparison with other compared methods. Therefore, it seems that the proposed IGSO has suitable performance especially in large scale systems.

Table 5 Prohibited operating zones for power-only units of Case I. Generator

Zone 1 (MW)

Zone 2 (MW)

Generator

Zone 1 (MW)

Zone 2 (MW)

1 2 3 4 5 6 7

[100–150] [120–160] [120–160] [70–100] [70–100] [70–100] [70–100]

[610–640] [310–340] [310–340] – – – –

8 9 10 11 12 13

[70–100] [70–100] [60–80] [60–80] [80–100] [80–100]

– – – – – –

6. Conclusion

proposed method, considering prohibited operating zones effects. According to Table 3 and Table 6, considering the prohibited operating zones will increase the total cost. The difference is mostly due to the difference in the produced power of units number 1, 3, 11 and 12. Since units number 1, 3, 11 and 12 were operating in the prohibited operating zone, their corresponding operating points are changed. Hence, it could be concluded that, considering prohibited operating zones constraint increases the total cost. Convergence characteristic of conventional and modified GSO algorithms are depicted in Fig. 9.

In this paper, IGSO algorithm is proposed to solve the CHPED problem. Performance of the proposed IGSO algorithm is compared with conventional GSO, GA and PSO for some benchmark functions. The comparison validates the ability of proposed IGSO to solve such a complicated function. In addition, the proposed approach was successfully implemented to the non-convex CHPED problem. In order to consider practical CHPED model, mentioned problem is solved considering valve-point effect, transmission losses, capacity limits of power units, heat–power mutual dependency of CHP units and prohibited operating zones effect. In the proposed framework, some modifications have been proposed to improve the convergence characteristic of the original GSO algorithm. Three large test cases with lots of local minima, used to demonstrate the applicability of the proposed method. These test systems are large enough to validate the performance of the IGSO. The obtained results for all

5.5. Comparative results This section presents comparative results of implementation the proposed method for small and large test cases. The numerical and comparative results are presented in Table 7. According to

Table 6 Comparison of simulation results for Case III. Output

GSO

IGSO

Output

GSO

IGSO

Output

GSO

IGSO

P1 P2 P3 P4 P5 P6 P7 P8 P9 P10

269.7504 360 77.0947 161.1803 116.4711 160.118 123.1151 162.2646 161.9556 113.852

268.9438 224.1527 294.9102 162.2591 110.2822 159.0552 158.8723 109.653 109.5521 268.9438

P11 P12 P13 P14 P15 P16 P17 P18 P19 H14

116.989 120 114.283 81 40 85.4377 40 10 36.4886 104.803

114.639 114.406 117.4323 118.2459 81.2429 40 81.3534 40 10 35

H15 H16 H17 H18 H19 H20 H21 H22 H23 H24

74.998 107.2898 74.998 40.001 20.6773 467.5825 60 60 119.8791 119.7711

104.94 74.998 104.99 74.998 40.001 20 470.09 60 60 120

Minimum Cost ($) Mean Cost ($) Maximum Cost ($)

GSO

IGSO

58650.283 58706.12 58763.915

58292.0551 58311.8439 58545.4748

455

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

Fig. 9. The cost convergence curve of the proposed IGSO and conventional GSO of Case III.

Table 7 Comparison of simulation results for small and large system. System number

1 2 3 4 5

System type

Small Small Small Large Large

Conventional GSO

CPSO [19]

Cost($)

Proposed IGSO Time (s)

Cost($)

Time (s)

Cost ($)

Time(s)

Cost($)

Time(s)

9257.087 11757.481 10088.517 58049.019 112320.41

1.159 1.25 2.77 35.54 70.65

9257.087 11759.04 10093.58 58225.745 113131.83

1.159 1.25 2.77 35.54 70.65

9257.08 11781.369 10325.3339 59736.26 119708.9

1.42 1.5 3.29 52.25 93.32

9257.07 11758.063 10100.3164 58122.746 117824.896

1.33 1.52 3.25 52.25 89.63

cases using the proposed method are always better than the earlier reported results. Comparison between run time results of recently published methods also indicates that applied modifications have improved convergence characteristic of the proposed method. In addition, based on the simulation studies, it could be concluded that the proposed modification improved the searching ability of the conventional GSO algorithm. On other hand, comparison between conventional and improved GSO, substantiate that the proposed modifications not only improved the convergence characteristic of GSO, but also IGSO is more efficient for large scale power systems. Based on the comparable differences between conventional and modified form of the GSO, it can be concluded that the IGSO appears to be a robust and reliable optimization technique for solving CHPED in large scale power system. References [1] Vasebi A, Fesanghary M, Bathaee S. Combined heat and power economic dispatch by harmony search algorithm. Int J Electr Power Energy Syst 2007;29:713–9. [2] Momoh JA. Electric power system applications of optimization. CRC Press; 2001. vol. 11. [3] Rau NS. Optimization principles: practical applications to the operation and markets of the electric power industry. John Wiley & Sons, Inc.; 2003. [4] Hemamalini S, Simon SP. Dynamic economic dispatch using Maclaurin series based Lagrangian method. Energy Convers Manage 2010;51(11):2212–9. [5] Spall JC. Introduction to stochastic search and optimization: estimation, simulation, and control, vol. 65. Wiley; 2005. [6] Chalermchaiarbha S, Ongsakul W. Stochastic weight trade-off particle swarm optimization for nonconvex economic dispatch. Energy Convers Manage 2013;70(6):66–75. [7] Cai J, Li Q, Li L, Peng H, Yang Y. A hybrid CPSO–SQP method for economic dispatch considering the valve-point effects. Energy Convers Manage 2012;53(1):175–81. [8] Mohammadi-ivatloo B, Rabiee A, Ehsan M. Time-varying acceleration coefficients IPSO for solving dynamic economic dispatch with non-smooth cost function. Energy Convers. Manage. 2012;56:175–83. [9] Özyön S, Aydin D. Incremental artificial bee colony with local search to economic dispatch problem with ramp rate limits and prohibited operating zones. Energy Convers Manage 2013;65(1):397–407.

TVACPSO [19]

[10] Jeddi B, Vahidinasab V. A modified harmony search method for environmental/ economic load dispatch of real-world power systems. Energy Convers Manage 2014;78(2):661–75. [11] Sashirekha A, Pasupuleti J, Moin N, Tan C. Combined heat and power (CHP) economic dispatch solved using Lagrangian relaxation with surrogate subgradient multiplier updates. Int J Electr Power Energy Syst 2013;44:421–30. [12] Su C-T, Chiang C-L. An incorporated algorithm for combined heat and power economic dispatch. Electric Power Syst Res 2004;69:187–95. [13] Khorram E, Jaberipour M. Harmony search algorithm for solving combined heat and power economic dispatch problems. Energy Convers Manage 2011;52:1550–4. [14] Song Y, Chou C, Stonham T. Combined heat and power economic dispatch by improved ant colony search algorithm. Electric Power Syst Res 1999;52:115–21. [15] Piperagkas G, Anastasiadis A, Hatziargyriou N. Stochastic PSO-based heat and power dispatch under environmental constraints incorporating CHP and wind power units. Electric Power Syst Res 2011;81:209–18. [16] Subbaraj P, Rengaraj R, Salivahanan S. Enhancement of combined heat and power economic dispatch using self adaptive real-coded genetic algorithm. Appl Energy 2009;86:915–21. [17] Basu M. Combined heat and power economic dispatch by using differential evolution. Electric Power Comp Syst 2010;38:996–1004. [18] Basu M. Bee colony optimization for combined heat and power economic dispatch. Expert Syst Appl 2011;38:13527–31. [19] Mohammadi-Ivatloo B, Moradi-Dalvand M, Rabiee A. Combined heat and power economic dispatch problem solution using particle swarm optimization with time varying acceleration coefficients. Electric Power Syst Res 2013;95:9–18. [20] Salgado F, Pedrero P. Short-term operation planning on cogeneration systems: a survey. Electric Power Syst Res 2008;78:835–48. [21] Abdolmohammadi HR, Kazemi A. A Benders decomposition approach for a combined heat and power economic dispatch. Energy Convers Manage 2013;71:21–31. [22] Rao PN. Combined heat and power economic dispatch: a direct solution. Electric Power Comp Syst 2006;34:1043–56. [23] Zare K, Haque MT, Davoodi E. Solving non-convex economic dispatch problem with valve point effects using modified group search optimizer method. Electric Power Syst Res 2012;84:83–9. [24] He S, Wu QH, Saunders J. Group search optimizer: an optimization algorithm inspired by animal searching behavior. IEEE Trans Evolut Comput 2009;13:973–90. [25] Moradi-Dalvand M, Mohammadi-Ivatloo B, Najafi A, Rabiee A. Continuous quick group search optimizer for solving non-convex economic dispatch problems. Electric Power Syst Res 2012;93:93–105.

456

M.T. Hagh et al. / Energy Conversion and Management 80 (2014) 446–456

[26] Victoire TAA, Jeyakumar AE. Reserve constrained dynamic dispatch of units with valve-point effects. IEEE Trans Power Syst 2005;20:1273–82. [27] Niknam T, Mojarrad HD, Meymand HZ. Non-smooth economic dispatch computation by fuzzy and self adaptive particle swarm optimization. Appl Soft Comput 2011;11:2805–17. [28] Li L, Liu F. Group search optimization for applications in structural design. Springer; 2011. vol. 9. [29] Mustard D. Numerical integration over the. Math Comput 1964;18:578–89.

[30] Mohammadi-Ivatloo B, Rabiee A, Soroudi A, Ehsan M. Iteration PSO with time varying acceleration coefficients for solving non-convex economic dispatch problems. Int J Electr Power Energy Syst 2012;42:508–16. [31] Meng K, Wang HG, Dong Z, Wong KP. Quantum-inspired particle swarm optimization for valve-point economic load dispatch. IEEE Trans Power Syst 2010;25:215–22. [32] Hagh MT, Teimourzadeh S, Alipour M. Combined heat and power dispatch using modified group search optimization method.