A lagrange relaxation-based alternating iterative algorithm for non-convex combined heat and power dispatch problem

Electric Power Systems Research 177 (2019) 105982

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

A lagrange relaxation-based alternating iterative algorithm for non-convex combined heat and power dispatch problem

T

Jianhua Chena,⁎, Yao Zhangb a b

State Grid Jibei Electric Power Company Limited, Xicheng District, Beijing, China China Nuclear Power Engineering Co., Ltd., Beijing, China

ARTICLE INFO

ABSTRACT

Keywords: Economic load dispatch Alternating iterative Combined heat and power Non-convex operating area

This paper proposes a Lagrange relaxation-based alternating iterative (AI) algorithm for economic dispatch problem with non-convex operating characteristic of combined heat and power (CHP) units. In this algorithm, the non-convex operating area of each CHP unit is divided into multiple convex piecewise sub-areas. Then, a bigM based method is adopted to transform the piecewise constraints to continuous linear constraints. Then, an alternating iterative algorithm is proposed to decompose the non-concave bilinear term in the objective to linear term. And the original combinatorial optimization problem is converted into the conventional quadratic programming problem for each unit, which can be easily solved by comparing the symmetry axis and the upper and lower bounds. The proposed algorithm can reach the global optimum point with a high efficiency, which is suitable for online application. Numerical results proved the validity and efficiency of the proposed method.

1. Introduction Combined heat and power (CHP) units, which produce both electricity and heat to consumers with high efficiency, are now playing an increasingly important role in power system. In North China Grid, almost all the large thermal units have CHP function. However, economic load dispatch (ED) problem with CHP units introduces technical challenge because of the non-linear and non-smooth operating characteristics in both the objective functions and the constraints. Few work has been done in the non-convex CHP economic dispatch area, and they can be classified into two classes: the mathematical programming methods [1–5] and the meta-heuristic optimization methods [6–17], which are briefly listed in Table 1. In the aspect of mathematical programming methods, [1] proposes a fast Newton-based iteration algorithm to solve the static economic dispatch problem, in which the separability of the cost function and the constraints is explored. In [2], a clear physical interpretation of the heat-power feasible region constraints of co-generation units in the dispatch process has been given. Based on this, a Lagrange multiplier based equal incremental rate method is adopted to solve the static optimization problem. However, the method is difficult to take the transmission capacity constraints into consideration. In [3], a surrogate sub-gradient multiplier updates method is proposed to improve the convergence speed of Lagrange relaxation method. In [4], a multi-

⁎

objective CHP dispatch problem containing cost and emission minimization with the considerations of demand response is proposed. And the ε-constraint method is used to solve it. In [5], an optimal robust short-term scheduling method considering the integrated heat and power microgrids is proposed and the robust optimization is applied to provide the optimal solutions avoiding the price deviation. Considering the complexity of conventional optimization process, the meta-heuristic optimization methods have been studied by some researchers. The genetic algorithm is applied in [6–9] to solve the cogeneration economic dispatch problem. In [6], a self-adaptive realcoded genetic algorithm is implemented, and the self-adaptation is achieved by tournament selection along with simulated binary crossover. In [7], a steady state real-coded genetic algorithm is proposed, in which the worst performer is determined and thrown out and replaced by a new member after every function evaluation. In [8], a hybrid real coded genetic algorithm is formed to improve the local search capability of genetic algorithm by hybridizing real coded genetic algorithm with ‘uniform random’ local search. In [9], an improved genetic algorithm using novel crossover and mutation to solve the combined heat and power economic dispatch problems. An oppositional teaching learning-based optimization method for combined heat and power dispatch is presented in [10] to solve combined heat and power dispatch problem with bounded feasible operating region. In [11], a harmony search algorithm is proposed to solve combined heat and power

Corresponding author. E-mail address: [email protected] (J. Chen).

https://doi.org/10.1016/j.epsr.2019.105982 Received 9 November 2018; Received in revised form 24 May 2019; Accepted 2 August 2019 Available online 13 August 2019 0378-7796/ © 2019 Elsevier B.V. All rights reserved.

2009

2012

2014

2019

2014

2011

2016

2016

2016

2015

2015

2019 2019

SARGA [6]

SSGA [7]

RCGAu [8]

IGA-NCM [9]

OTLBO [10]

HAS [11]

GSO [12]

Gwo [13]

COA [14]

MPSO [15]

CSO [16]

Woa [17] Proposed AI method in this paper

Year

A self-adaptive real-coded genetic algorithm is implemented to increase the probability towards the global optimum and prevent premature. A steady state real-coded genetic algorithm is proposed, which has a higher precision for complex problems than binary GA A hybrid real coded genetic algorithm is formed to improve the local search capability of genetic algorithm An improved genetic algorithm using novel crossover and mutation to enhance the ability of global and local searches in solution space An oppositional teaching learning-based optimization method is firstly used to solve combined heat and power dispatch problem A harmony search algorithm is firstly proposed to solve combined heat and power economic dispatch problem A group search optimization is firstly proposed to solve the complex nonconvex combined heat and power economic dispatch problem A grey wolf optimization for combined heat and power dispatch with cogeneration systems is firstly proposed A cuckoo optimization algorithm is firstly implemented to solve combined heat and power dispatch problem. A modified particle swarm optimization is proposed, where Gaussian random variables are used in the velocity term to improve the search efficiency and obtain the global optimum A novel crisscross optimization algorithm is implemented to solve the large scale CHPED problem Whale optimization algorithm is employed to solve CHPED problem A Lagrange relaxation-based alternating iterative algorithm for non-convex CHP dispatch problem is proposed to solve the non-convex and non-smooth characteristics of CHP units in both the objective functions and the constraints

Main contribution

Table 1 Comparison of recent research works in CHP economic dispatch area.

(1) Could reach the global optimal point fast enough for online application (2) Only needs to run once, and the results are repeatable. (3) The results are insensitive to the parameter settings.

Have rarely no restrictions on the shape of the cost function and constraints, hence they can solve some special problems where the mathematical programming methods cannot do.

Pros

(1)Time-consuming and the results are not repeatable, i.e. different results for each trial. (2)The results are very sensitive to the parameter settings.

Cons

J. Chen and Y. Zhang

Electric Power Systems Research 177 (2019) 105982

2

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

for CHP units, the model is NP-hard to solve. An alternating iterative method is proposed to realize decoupling optimization for multivariable model. As the decoupled models are all single-unit for heat and power, they can be easily solved by comparing the symmetry axis and the upper and lower bounds for each unit instead of large-scale optimization, which could improve the efficiency again compared with the conventional methods. (3) When the valve point effect for conventional units is considered, a big-M based piecewise linearization method is proposed to convert the non-linear and non-smooth cost function into the conventional mixed-integer quadratic function. Compared with the conventional piecewise linearization method [22], the proposed big-M based method could reduce the scale of the optimization problem, which is important for improving the efficiency of the ED problems. The proposed method possesses the advantages of the mathematical optimization technique, which could reach the global or very near global optimum point and only needs to run once. Besides, the algorithm runs fast enough for online application requirement. Meantime, the results are insensitive to the parameter settings. The paper is organized as follows. Section 2 introduces the characteristics of the CHP units and formulates the dispatch model with CHP units. Section 3 presents the proposed optimization algorithm. Numerical simulations are described in Section 4. Section 5 concludes the paper.

Fig. 1. A typical 250 MW CHP unit operating characteristic of North China Grid.

economic dispatch problem. In [12], a group search optimization is proposed to solve the complex non-smooth non-convex combined heat and power economic dispatch problem, where the valve-point effects and prohibited operating zones for conventional units are considered. In [13], a grey wolf optimization for combined heat and power dispatch with cogeneration systems is proposed. In [14], cuckoo optimization algorithm is implemented to solve energy production cost minimization in a combined heat and power (CHP) generation system 15. resents a modified particle swarm optimization to solve the non-smooth nonconvex combined heat and power economic dispatch problem, where Gaussian random variables are used in the velocity term to improve the search efficiency and obtain the global optimum. In [16], a novel crisscross optimization algorithm is implemented to solve the large scale CHPED problem. Whale optimization algorithm is employed in [17] to solve CHPED problem, and the valve-point loading impact is studied. The meta-heuristic programming methods have the merits of rarely no restrictions on the shape of the cost function and constraints, hence they can solve some special problems where the mathematical programming methods cannot do. However, the meta-heuristic programming methods are usually time-consuming and the results are not repeatable, i.e. different results for each trial. Besides this, the results are very sensitive to the parameter settings. For online application, the ED algorithm with CHP unit consideration is required to reach the global optimum points as fast as possible, and what’s more, the results must be repeatable. However, the algorithm with both the two characteristics simultaneously is still lack. This paper presents a new method to solve the ED problem with non-linear and non-smooth operating characteristic of CHP units, which has a high efficiency for online application. In this algorithm, the nonconvex operating area of each CHP unit is divided into multiple convex piecewise sub-areas. Then, a big-M based method is adopted to transform the piecewise constraints to continuous linear constraints. Then, a Lagrange relaxation-based alternating iterative algorithm is proposed to decompose the non-concave bilinear term in the objective to linear term. Based on these, the CHP dispatch problem is decomposed into two sub-problems, the heat dispatch sub-problem and power dispatch sub-problem. And the original combinatorial optimization problem is converted into the conventional quadratic programming problem for each unit, which can be easily solved by comparing the symmetry axis and the upper and lower bounds. The contributions of the paper can be summarized as follows: (1) A big-M based method is proposed to equivalently converting the non-linear and non-smooth operating regions constraints of CHP units into continuous convex linear constraints. Compared with the conventional method, the proposed big-M based method could reduce the scale of the optimization problem greatly, which is important for improving the efficiency of the ED problems. (2) As the objective function contains a non-concave bilinear term

2. Problem formulations and modeling 2.1. CHP unit output characteristic Fig. 1 illustrates a typical output characteristic of a 250 MW CHP unit, which is non-convex and non-linear. And it can be depicted by (1) and (2).

pjmin (hj )

pj

pjmax (hj )

(1)

hjmin (pj )

hj

hjmax (pj )

(2)

where j denotes the index of each CHP unit; pj and hj are respectively the power and heat supply of unit j; pjmin , pjmax are respectively the

active power output lower and upper bound of unit j; hjmin , hjmax are respectively the heat output lower and upper bound of unit j. It can be seen that the power output limits of CHP units are functions of the heat output, and the heat output limits are functions of the power output. The non-convex and non-linear feasible region constraints make the economic dispatch problem complicated. Hence, Eqs. (1) and (2) should first be linearized. Firstly, the dashed lines are used here to divide the operating area into multiple convex sub-areas, the heat-power feasible region constraint in sub-area l can be depicted by (3):

_k j, l hj + _ j, l

k¯j, l hj + ¯j, l , hjmin ,l

pj

hj

hjmax ,l

(3)

where l represents the index of each sub-area of CHP unit j; k 、 k and 、 are the mutual dependency coefficients of power and heat conmax straint. hjmin are respectively the heat output lower and upper , l , h j, l bound in sub-area l of unit j. The heat-power feasible region constraint in all the sub-areas can be depicted by (4):

_k j,1 hj + _ j,1 or , ..., or ,

pj

k¯j,1 hj + ¯j,1, hjmin ,1

_k j, l hj + _ j, l or , ..., or ,

pj

k¯j, l hj + ¯j, l , hjmin ,l

_k j, nj hj + _ j, n j

pj

k¯j, nj hj + ¯j, nj , hjmin , nj

hjmax ,1

hj hj

hjmax ,l hj

hjmax , nj

(4)

Where nj is the total number of sub-areas of CHP unit j, l = 1,2,…,nj. “or” implies that each CHP unit can operate in only one of the sub3

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

areas.

3. Solution procedure

2.2. Economic dispatch model with non-convex CHP units

In this paper, we propose a new big-M and Lagrange relaxationbased alternating iterative (AI) algorithm. Firstly, a big-M based coding scheme [18] is applied here to transform the piecewise constraint (4) to continuous linear constraint. Then, the Lagrange relaxation method [19] is used to decompose the problem into master problem and singleunit sub problem by relaxing the coupling constraints. With the proposed alternating iterative method, the sub problem is further decomposed into power dispatch sub problem and heat dispatch sub problem, both of which are the conventional quadratic programming problem without coupling constraints. And they can be easily solved by comparing the symmetry axis and the upper and lower bounds for each unit. Firstly, based on the big-M method [18], the piecewise constraint (4) can be equivalently converted to mixed integer linear constraint (12).

The ED problem with CHP units is formulated as follows: 1) Objective function The cost function is comprised of two parts: the cost of conventional thermal units and the cost of CHP units. The cost of conventional thermal units can be generally depicted by a quadratic function: n

f (pi ) = i=1

(ai pi2 + bi pi + ci )

(5)

where i denotes the index of conventional thermal unit; n is the total number of units, i = 1,2,…,n; ai , bi, ci are the production cost coefficients of unit i. When the valve point effect is considered, an additional cost item should be added to f (pi ) as follows, which makes the cost function nonsmooth.

hjmax ,l

k j, l h j +

n

f (pi ) = i=1

(ai pi2 + bi pi + ci + ei | sin fi (pi

pi )|)

k j, l h j +

(6)

m j =1

(aj pj2 + bj pj + cj + dj h2j + ej h j + f j + gj pj h j )

(7)

where m is the total number of CHP units, j = 1,2,…,m; dj , ej , f j , gj are the heat cost coefficients of unit j. 1) Constraints 2) The power-load balance constraint is: n

m

pi + i=1

pj = D

(8)

j =1

Where D is the load demand of the system. 3) The heat balance constraint is: m

hj = S

(9)

j=1

Where S is the heat demand of the system. 4) unit capacity constraint:

pi

pi

p¯i

(10)

Where _pi 、 p¯i are respectively the lower and upper generation output of unit i. 5) transmission capacity constraint: n

TLl

m

(kli pi ) + i=1

(klj pj ) j =1

TLl , l = 1, …, L

j, l j, l

M (xj, l-1

+ M (xj, l-1

x j, l + 1) xj, l + 1)

x j, l )

pj (12)

where M is a big real number; x j, l is a binary variable. From the big-M method, it is known that when M is big enough, the result is insensitive to the value of M. However, in practical applications, M can be conservatively taking the value of max (pjmax , hjmax ) , which is a sufficient but not necessary condition. As can be seen from (12), the right-hand side and the left-hand side of the first constraint are complementary, i.e. if x j, l = 1, the right-hand side constraint will become effective; while if x j, l = 0, the left-hand side constraint will become effective, which means that the constraints on both sides cannot take effect at the same time. The first and the second constraints of (12) implies that the generator can operate in only one of the sub-areas each time, i.e. the interval with x j, l =0 and x j, l + 1 = 1. And the last constraint makes sure that only in this sub-area the heat-power feasible region constraint be effective. The proposed big-M based model (12) is very useful for reducing the scale of the optimization problem. For example, if CHP unit j has 5 piecewise feasible sub-areas, i.e. nj = 5, then four binary variables x j,1, xj,2 , xj,3 , x j,4 are needed in (12). As each of them could take 0 or 1 value, there are totally 24 = 16 combinations. However, subject to the third constraint of (12), only five combinations are effective, i.e. 0000, 0001, 0011, 0111, 1111. And the remaining 11 combinations are not allowed. Hence the scale of the problem is greatly reduced. Likewise, when nj = 10, only 11 combinations are effective out of totally 210 = 1024 combinations. And when nj becomes larger, the scale of the optimization problem will be reduced more. This is important for improving the efficiency of the proposed alternating iterative method and will be mentioned later. Substituting (12) for (4), the ED model in section II.B becomes a mixed integer programming problem. However, as the objective function (7) contains a bilinear term pj h j , the model is still NP-hard to solve. A Lagrange relaxation-based alternating iterative method is proposed for it. The details of the algorithm are described as follows. 1) As there are only three system-wide constraints: (8), (9) and (11) and all the other constraints refer to only one unit, the Lagrange relaxation method is used here to decompose the problem into many single-unit optimization problems [1]. The Lagrange dual expression of the dispatch problem in this paper can be expressed as (13).

The cost of CHP units is comprised of three parts: the power production cost, the heat production cost and the mutual dependency cost between heat and power:

f (pj , h j ) =

hj hjmax + M (1 ,l xj, l-1 x j, l xj,0 = 0,xj, nj = 1

Mxj, l

(11)

Where TLl 、 TLl are respectively the lower and upper limit of transmission interface l with bus loads excluded. L is the total number of transmission interfaces, l = 1,2,…,L. kli is the quasisteady-state generation distribution shift factor of unit i to transmission interface l. 6) CHP unit output characteristic constraints are given by Eq. (4). Eqs. (4–11) constitute a nonlinear piecewise optimization model, which cannot be directly solved efficiently. Hence, a new Lagrange relaxation-based alternating iterative (AI) algorithm is proposed in Section 3.

4

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

max: (wl , wl, v, ), wl

0, wl

n

0

m

(w, v, ) = inf

f (pi ) +

i=1 m

j=1

n

f (pj ,

h j)

v( i=1

mi = ceil(K

m

pi +

j=1

pj

D) i, k

hj

( n

m

wl (TLl l=1

i=1

(kli pi )

j=1

L

(klj pj ))

i, k

n

wl (

j=1

ei | sin fi (p¯i, k

l=1

i=1

(klj pj ) TLl )|p

p¯i, k =

pi )|

ei | sin fi (p¯i, k

p¯i, k

p¯i, k

0, wl

0

n

L

(f (pi )

vpi +

(f (pj , h j )

vpj

i=1

(wl l=1

pi )|

k = 1, mi Kfi p¯i k = mi

wl ) kli pi )

m

+ j=1

(wl

wl ) klj pj ) + C |p

(14)

(wl

hjmax ,l

s. t . Mx j, l xj, l-1 xj, l xj,0 = 0,x j, nj = 1 k j, l h j + k j , l hj +

As f (pi ) is a quadratic function, the sub-problem is a conventional quadratic programming problem without coupling constraints. Hence, the optimal power output can be reached by just comparing the symmetry axis of the quadratic objective function expression and the upper and lower bounds of each unit, i.e.

v+

(wl l=1

wl ) kli)/2ai

(16)

p¯i, k + M (1

k j, l h j +

+

i, k

M (yi, k

1

xj, l + 1)

x j, l )

pj

x j, l + 1)

(19)

j, l

M (xj, l-1

v + gj h j +

xj, l + 1)

(wl

wl ) klj ) pj + cj )

l=1

pj

k j, l h j + xj, l + 1)

j, l

(20)

4) At the kth iteration, the power output of unit j is set to be the results of k-1th iteration pjk 1 and we could traverse every effective combination of the binary variables x j, l . The optimization model (19) then becomes a conventional quadratic programming model. And the optimal heat output of unit j can be easily reached by comparing the symmetry axis and the upper and lower bounds for each unit. 5) The heat output and x j, l of each unit is set to be the results of kth iteration, and the optimization model (20) then becomes a conventional quadratic programming model. The optimal power output can also be reached by comparing the symmetry axis and the upper and lower bounds for each unit as like in step 4). 6) This procedure is repeated and repeated until the objective

yi, k )

yi, k + 1)

+ M (1

+ M (xj, l-1

yi, k 1 yi, k yi,0 = 0,yi, mi = 1 i, k pi

+ M (xj, l-1

+ gj pj ) h j + f j )

L

v ) pi + ci + pi

j, l

M (xj, l-1

f2 (pj , h j ) = min(aj pj2 + (bj

When the valve point effect is considered, a big-M based piecewise linearization method is proposed to convert the non-linear and nonsmooth cost function (6) into a conventional mixed-integer quadratic function as follows.

Myi, k

j, l

hj

hjmax ,l

As it was mentioned above, only a few combinations of the binary variables are effective in this algorithm, we could traverse every combination of them. And in each combination, as the binary variables are given, the optimization model becomes a conventional quadratic optimization model. The power dispatch sub-problem is as follows:

L

(bi

(18)

h j,

f1 (pj , h j ) = min(dj hj2 + (ej

(15)

s. t . constraint (9)

s. t . p¯i, k

wl ) klj pj

As f (pj , contains a bilinear term pj the model is NP-hard to solve. Hence, an alternating iterative method is proposed as follows to solve it. Firstly, we decompose Eq. (18) into two single-unit sub-problems: the heat dispatch and the power dispatch. And the heat dispatch subproblem is as follows:

wl ) kli pi )

l=1

f (pi ) = ai pi2 + (bi

(wl

h j)

L

pi = max pi , min p¯i ,

hj +

vpj

s. t . constraint (11)

De}

Where, C is a constant. As can be seen, Eq. (14) is the sum of many subexpression, each with a unit, i.e. Eq. (14) is separable. Thus, we could decompose it into several parallel single-unit sub-problems. And each sub-problem only includes the variables and constraints with respect to one unit. Iterating between the master and sub problems, the optimal solution can be gradually approached. 2) From step 1), we know that the sub-problem is a single-unit problem. And the sub-problem for conventional thermal unit is:

vpi +

1

l=1

l=1

min(f (pi )

¯i, k i, k p

L

min f (pj , h j )

L

hj +

pi )|

1

1

where, K is the number of segments between two consecutive valve points. The principle of the conversion is basically same with (12), and will not be discussed again at here. Subject to the third constraint of (17), only a few combinations of the binary variables are effective in this algorithm, hence we could traverse every combination of them. And in each combination, as the binary variables are given, the optimization model becomes a conventional quadratic optimization model and (16) can be easily applied to solve it. 3) The sub-problem for CHP unit is:

Where, inf{} is the infimum of the expression; De represents the singleunit constraints, including constraints (10) and (12). wl , wl , v and represent the Lagrange relaxation coefficients corresponding to the coupling equality and inequality constraints (8), (9) and (11). Eq. (13) can be transformed to the following form.

(wl , wl, v , ) = inf

pi + k

De} (13)

max: (wl , wl, v, ), wl

)

= f (pi ) = ei | sin fi (p¯i, k

(kli pi )

m

+

pi )

S)

j=1 L

=

fi (p¯i

(17)

Where i, k , i, k are the piecewise cost coefficients; yi, k is a binary variable representing whether unit i operates in segment k; _pi, k , p¯i, k are respectively the lower and upper bounds of segment k; mi represents the number of segments for unit i. Parameters i, k , i, k , mi and p¯i, k can be calculated as follows [22]. 5

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

Table 2 Comparison of performance with AI and other methods.

SARGA [6] SSGA [7] RCGAu [8] IGA-NCM [9] AI

CPU time(s)

Costmin

Costmax

Costmean

Number of iteration

2.81 3.15 5.68 2.66 0.30

9257.083 9257.8153 9257.0869 9257.075 –

9380.302 9426.561 9409.1808 9257.9014 –

9276.9606 9308.2587 9300.9727 9257.1553 9257.075

288 334 732 255 19

Table 3 Optimization results for each unit.

P1(MW) P2(MW) P3(MW) h2(MWth) h3(MWth) h4(MWth)

function is no longer reduced. The flow chart of the algorithm is as follows (Fig. 2). From the iteration procedure, it can be seen that:

f2 (hjk , pik )

f1 (hjk 1, pik 1 )

(21)

f2 (hjk , pik 1 )

(22)

Adding (18) to (19), the following expression is reached:

f1 (hjk , pik 1 ) + f2 (hjk , pik )

f1 (hjk 1, pik 1 ) + f2 (hjk , pik 1 )

And when the common term

gj pjk 1 hjk

RCGA-IMM [20]

The marginal cost

0 160 40 40 75 0

0 160 40 40 75 0

50 26.78 40.305 11.56 5.09 23.4

heuristic optimization methods needs to run several times to avoid the impact of a few extreme optimization results. Meantime, the number of iterations for different methods is also different, which are respectively 288, 334, 732, 255 and 19 with SARGA, SSGA, RCGAu, IGA-NCM and AI. As the number of runs for AI method is the least, the calculation time is also the least. Hence, the computation efficiency of AI method is the best, which takes only 0.30 s, while the other four algorithms take respectively 2.81 s, 3.15 s, 5.68 s and 2.66 s. The convergence characteristic of the proposed AI method is shown in Fig. 3. It can be seen from Fig. 3 that the global optimal solution can be approached by just about 19 iterations, which is fast and suitable for online application. From Table 3, as the marginal cost of conventional thermal generator #1 and the heat-only unit #4 are higher than CHP units #2 and #3 at the optimal point, hence, p1 and h4 take the value of 0. Meantime, as the marginal cost of power output for CHP unit #2 is lower than CHP unit #3, p3 takes the value of 40, which is the minimum value in the heat-power feasible region. Fig. 4 depicts the iterative process of p2 and h2 with the AI method. It can be seen from Fig. 4 that even when the start point is far from the optimal point, the algorithm still has a good convergence. And the first few iterative steps are relatively large, which gradually become smaller when it approaches the optimal point. Fig. 5 and Table 4 show the fuel cost of the system with different transmission interface flow constraints for unit #1 and #2. It can be seen from Fig. 5 that the fuel cost is influenced greatly by

Fig. 2. Flow chart of the proposed algorithm.

f1 (hjk , pik 1 )

AI

(23)

at both sides is eliminated, we

get f (pik, t , hik, t ) f (pik, t 1 , hik, t 1) , which means that the objective function is gradually decreased in every iteration and the optimal point can be reached with finite iteration. 4. Numerical tests

The performance of the proposed Lagrange relaxation-based alternating iterative method was respectively tested with a 4-unitsystem [2] and a 48-unit large scale system [21]. The simulation tool used is Matlab on a 2.0 GHz personal computer with 2 G RAM. 4.1. Simulation results on a four-unit test system The test system is comprised of a conventional thermal generator, two CHP units and a heat-only unit. The power demand of the system is 200 MW, and the heat demand is 115MWth. M takes the value of 247. Table 2 lists the comparison of the performance with different algorithms. Table 3 gives the optimization results for each unit. From Table 2, we can see that the proposed AI method gives the minimum fuel cost compared with the mean value of SARGA, SSGA, RCGAu and IGA-NCM, which is also the global optimal solution as reported in [9,10]. What’s more, the proposed AI method is repeatable, while the other meta-heuristic optimization methods are not. Thus, the proposed AI method only needs to run once, and the other meta-

Fig. 3. Convergence characteristics of the proposed algorithm. 6

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

Fig. 4. Iterative process of the CHP unit #2.

Fig. 6. Fuel cost of the system with different load demands and efficiencies of electric heating equipment.

Fig. 5. Fuel cost of the system with different transmission interface capacity. Table 4 Comparison of results with different transmission interface capacities. Transmission interface capacities(MW)

100

120

150

Obj. value($) P1(MW) P2(MW) P3(MW) h2(MWth) h3(MWth) h4(MWth)

10,280.89 0 100 100 5.35 109.65 0

9825.39 0 120 80 5.47 109.53 0

9352.11 0 150 50 31.37 83.63 0

Fig. 7. Fuel cost of the system with different proportion of wind power for heating.

Sometimes, part of the electricity, which is used for heating, comes from wind power instead of thermal units. Fig. 7 gives the fuel cost of the system with different proportion of wind power for heating, where the efficiency of the electric heating equipment is supposed to be 50%. It can be seen from Fig. 7 that, when the proportion of wind power for heating is low, the fuel cost increases along with load demand. However, when the proportion of wind power for heating is high, the fuel cost decreases when the load demand increases. This is because most of the customers’ heat demand is come from wind power, and although the load demand is still increasing, the heating cost and the total fuel cost is decreasing.

the transmission interface flow constraints. And when the transmission interface capacity is higher than 160 MW, the fuel cost is $9257.1, which is the same with Table 2 where the transmission interface is not considered. Meanwhile, when the transmission interface capacity is lower than 160 MW, the fuel cost will start to increase, which shows the importance of the transmission interface capacity for reducing the cost. Table 4 shows that when the transmission interface capacity increases, the power output and heat output of unit #2 are gradually increasing, which is just the same value as the transmission interface capacity; while the power output and heat output of unit #3 are decreasing. This is because the marginal cost of unit #2 are lower than units #3, hence, the power output of unit #2 should be higher than unit #3 if not being limited by the transmission interface power flow constraints. From the perspective of environmental protection, Chinese government wants customers use more electricity for heating, instead of just heating by coal. This will result in the growth of load demand and decline in heat demand of the system. Fig. 6 gives the fuel cost of the system with different load demands and efficiencies of electric heating equipment. From Fig. 6, the fuel cost increases along with load demand, and which is opposite with the efficiency of electric heating equipment. This implies that a high efficiency is necessary not only for customers but also for electric companies.

4.2. Simulation results on a 48-unit test system [21] The test system is comprised of 26 conventional power units, 12 CHP units and 10 heat-only units. The power demand of the system is 4700 MW, and the heat demand is 2500MWth. M takes the value of 2700. The optimization results comparison is listed in Table 5. Table 6 Table 5 Comparison of performance with AI and other methods.

COA [14] OTLBO [10] MPSO [15] GSO [12] CSO [16] AI

7

CPU time(s)

Costmin

Costmax

Costmean

81.3027 36.43 9.2516 9.51269 95.21 3.29

116,789.91535 116,579.239 116,465.5395 116,457.9578 115,967.7205 –

117068.2693 116,649.4473 116,482.4404 116,473.2183 116,047.2154 –

116,835.5491 116,613.6505 116,471.3609 116,463.6522 115,995.8842 115,851.0504

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang

Table 6 Optimization results for each unit. Output

COA

OTLBOa

MPSO

GSO

CSO

AI

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 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

628.0052 223.367 224.2239 160.4662 109.41 158.5215 108.0529 109.6997 159.5992 40.0587 40.0024 55.0956 92.0298 447.8708 224.4969 73.9883 159.5506 109.6294 159.3689 159.3796 158.4264 159.3175 41.4182 40.1996 90.8237 92.4161 84.0293 97.5254 81.3338 46.9334 10.2119 38.5893 82.9338 40.487 81.7879 44.8094 11.8972 54.0504 106.2509 124.6536 104.8896 80.7382 40.0726 21.4735 105.8367 75.4084 105.2044 79.1119 40.4439 28.4089 444.7157 59.5783 60 119.8104 119.665 427.0172 59.1377 59.2648 119.4493 118.8751

628.3199 225.3313 223.9653 159.8516 109.915 159.7795 109.8946 109.9321 159.9569 40.897 41.3115 55.1748 92.4003 448.8359 225.7871 75.46 160.1192 110.3532 159.819 159.7765 159.737 160.1751 40.114 40.3042 92.4149 92.5012 85.9857 98.5005 81.7197 48.9055 10.0832 39.311 82.0236 40.1104 81.3039 45.67 13.8709 30.3881 107.5951 125.4997 105.1942 82.6853 40.0346 21.9568 105.3622 75.0938 104.9667 79.8936 41.6554 17.9018 445.0937 59.9967 59.9974 119.8834 119.5231 428.7605 59.9957 59.9638 119.5025 119.444

179.7773 299.9621 290.3748 179.5655 160.0408 119.7704 162.9561 159.7361 119.8672 114.8494 114.9528 90.8763 92.7706 92.7706 226.5089 299.2174 159.9777 109.8667 159.7483 170.6698 160.6908 109.9663 114.837 116.0113 92.493 92.4025 82.4557 44.3773 97.1621 44.6356 10 35.7983 82.5995 44.2751 91.0334 45.0766 10.0002 35.1411 105.6169 78.8068 113.8706 79.0298 40.0004 20.3628 105.6984 8.7186 110.4315 79.4109 40.0004 20.0641 522.1121 60 59.9999 120 119.9996 385.8784 59.9997 60 119.9995 119.9996

448.9126 150.5151 80.766 160.0923 109.9592 159.852 160.1104 159.8453 160.0219 114.9957 115.1906 92.6482 55.042 269.4783 299.4636 299.7375 159.9635 159.8998 159.7568 60.0218 160.0075 159.9142 114.9146 40.116 93.87 93.6315 89.9223 49.6516 81.2954 46.6966 10.0212 37.7288 92.038 50.4524 95.2834 52.3657 10.0683 45.7741 109.8046 83.3599 104.961 80.8014 39.9976 21.2295 110.9901 84.0301 112.7913 85.6985 40.0211 24.8763 458.7095 59.9975 60 119.9632 119.999 422.7929 59.9792 59.9974 120 120

538.5512 224.4094 224.4313 159.7336 60.0128 60.0051 109.8705 109.8672 159.7445 114.81 77.4097 55.0031 55 628.3206 149.6535 360 109.8693 60.0014 109.8408 159.7348 159.7384 60 77.4041 114.8152 92.404 92.4144 82.9843 40 81.4314 40 10 35.0003 81.4693 40 81.0691 40 10.0007 35 105.9126 75 105.0423 75 40 20 105.0634 75 104.8372 75 40 20 464.9032 59.9999 60 119.9994 119.9999 474.2431 59.9997 59.9998 119.9997 119.9998

619.8595 299.1993 299.1993 109.1237 109.1237 109.1237 109.1237 109.1237 109.1237 40 40 55 55 619.8595 299.1993 299.1993 109.1237 109.1237 109.1237 109.1237 109.1237 109.1237 40 40 55 55 81 40 81 40 10 35 81 40 81 40 10 35 104.8 75 104.8 75 40 20 104.8 75 104.8 75 40 20 470.4 60 60 120 120 470.4 60 60 120 120

Fig. 8. Convergence characteristics of the proposed algorithm on 48-unit test system.

system in section A. 5. Conclusions This paper proposes a Lagrange relaxation-based alternating iterative method to formulate the economic dispatch problem with nonconvex operating characteristic of combined heat and power units. In this algorithm, the non-convex operating area of each CHP unit is divided into multiple convex piecewise sub-areas. And a big-M based method is adopted to transform the piecewise constraints to continuous linear constraints. Then, with the alternating iterative method, the nonconcave bilinear term in the objective is decomposed to linear term, and the original optimization problem is converted into the conventional quadratic programming problem for each unit, which can be easily solved by comparing the symmetry axis and the upper and lower bounds. Numerical tests proved the validity and efficiency of the proposed method. Conflict of interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled. References [1] F.J. Rooijers, R.A.M. Amerongen, Static economic dispatch for Co-generation systems, IEEE Trans. Power Syst. 9 (August (3)) (1994) 1392–1398. [2] Tao Guo, M.I. Henwood, M. Ooijen, An algorithm for combined heat and power economic dispatch, IEEE Trans. Power Syst. 11 (November (4)) (1996) 1778–1784. [3] A. Sashirekha, J. Pasupuleti, N.H. Moin, C.S. Tan, Combined heat and power (CHP) economic dispatch solved using Lagrangian relaxation with surrogate subgradient multiplier updates, Electr. Power Energy Syst. 44 (2013) 421–430. [4] M. Nazari-Heris, S. Abapour, B. Mohammadi-Ivatloo, Optimal economic dispatch of FC- CHP based heat and power micro-grids, Appl. Therm. Eng. (2016), https://doi. org/10.1016/j.applthermaleng.2016.12.016. [5] M. Nazari-Heris, B. Mohammadi-Ivatloo, G.B. Gharehpetian, et al., Robust shortterm scheduling of integrated heat and power microgrids, IEEE Syst. J. (June) (2018) 1–9. [6] P. Subbaraj, R. Rengaraj, S. Salivahanan, Enhancement of combined heat and power economic dispatch using self adaptive real-coded genetic algorithm, Appl. Energy 86 (6) (2009) 915–921. [7] J.D. Dyer, R.J. Hartfield, G.V. Dozier, J.E. Burkhalter, Aerospace design optimization using a steady state real-coded genetic algorithm, Appl. Math. Comput. 218 (9) (2012) 4710–4730. [8] B.A. Sawyerr, A.O. Adewumi, M.M. Ali, Real-coded genetic algorithm with uniform random local search, Appl. Math. Comput. 228 (2014) 589–597. [9] D. Zou, S. Li, X. Kong, et al., Solving the combined heat and power economic dispatch problems by an improved genetic algorithm and a new constraint handling strategy, Appl. Energy 237 (2019) 646–670.

gives the optimization results for each unit. It can be seen from Table 5 that optimization cost with the proposed AI method is 115,851.0504, which is the minimum fuel cost compared with the mean cost of COA, OTLBO, MPSO, GSO and CSO. Meantime, the running time is only 3.29 s with the proposed method, while the time are respectively 81.30 s, 36.43 s, 9.25 s, 9.51 s and 95.21 s with the other five methods. Hence, the computation efficiency of AI method is the best. The convergence characteristic of the proposed AI method is shown in Fig. 8. It can be seen from Fig. 8 that the optimal points can be approached by just about 20 iterations, which is similar with the 4-unit 8

Electric Power Systems Research 177 (2019) 105982

J. Chen and Y. Zhang [10] P.K. Roy, C. Paul, S. Sultana, Oppositional teaching learning based optimization approach for combined heat and power dispatch, Int. J. Electr. Power Energy Syst. 57 (May) (2014) 392–403. [11] E. Khorram, M. Jaberipour, Harmony search algorithm for solving combined heat and power economic dispatch problems, Energy Convers. Manage. 52 (2011) 1550–1554. [12] M. Basu, Group search optimization for combined heat and power economic dispatch, Int J Electr Power Energy Syst 78 (2016) 138–147. [13] N. Jayakumar, S. Subramanian, S. Ganesan, E.B. Elanchezhian, Grey wolf optimization for combined heat and power dispatch with cogeneration systems, Int. J. Electr. Power Energy Syst. 74 (2016) 252–264. [14] M. Mehdinejad, B. Mohammadi-Ivatloo, R. Dadashzadeh-Bonab, Energy production cost minimization in a combined heat and power generation systems using cuckoo optimization algorithm, Energy Effic. (2016) 1–16. [15] M. Basu, Modified particle swarm optimization for non-smooth non-convex combined heat and power economic dispatch, Electr Power Compon Syst 43 (2015) 2146–2155. [16] A. Meng, P. Mei, H. Yin, X. Peng, Z. Guo, Crisscross optimization algorithm for solving combined heat and power economic dispatch problem, Energy Convers. Manage. 105 (2015) 1303–1317. [17] M. Nazari-Heris, M. Mehdinejad, B. Mohammadi-Ivatloo, et al., Combined heat and power economic dispatch problem solution by implementation of whale optimization method, Neural Comput. Appl. 31 (2) (2019) 421–436. [18] Tao Ding, Rui Bo, Wei Gu, Big-M based MIQP method for economic dispatch with disjoint prohibited zones, IEEE Trans. Power Syst. 29 (March (2)) (2014) 976–977. [19] F. Antonio, G. Claudio, Solving nonlinear single-unit commitment problems with ramping constraints, Oper. Res. 54 (4) (2006) 767–775.

[20] A. Haghrah, M. Nazari-Heris, B. Mohammadi-Ivatloo, Solving combined heat and power economic dispatch problem using real coded genetic algorithm with improved Mühlenbein mutation, Appl. Therm. Eng. 99 (2016) 465–475. [21] M. Nazari-Heris, B. Mohammadi-Ivatloo, G.B. Gharehpetian, A comprehensive review of heuristic optimization algorithms for optimal combined heat and power dispatch from economic and environmental perspectives, Renewable Sustainable Energy Rev. 81 (2) (2018) 2128–2143. [22] M.Q. Wang, H.B. Gooi, S.X. Chen, S. Lu, A mixed integer quadratic programming for dynamic economic dispatch with valve point effect, IEEE Trans. Power Syst. 29 (5) (2014). Jianhua Chen received the B.S. degree in electrical engineering and automation from Xi’an Jiaotong University of Xi’an, Xi’an, China, in 2008 and the Ph.D. degree in electrical engineering from Tsinghua University of Beijing, Beijing, China, in 2013. Currently, he is a senior Electrical Engineer at State Grid Jibei Electric Power Company, Beijing, China. His research interests include large-scale wind power dispatch and control theory, system security, and robust optimization.

9