Introducing environmental constraints in generation expansion problems

NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS Numer. Linear Algebra Appl. 2000; 00:1–6 Prepared using nlaauth.cls [Version: 2002/09/18 v1.02]

Introducing environmental constraints in generation expansion problems Rachel Marques Marcato1 and Claudia Sagastiz´abal1,2,∗ 1

IMPA Instituto de Matem´ atica Pura e Aplicada. Estrada Dona Castorina 110, Jardim Botˆ anico, Rio de Janeiro 22460-320, RJ, Brazil. 2 On leave from INRIA-Rocquencourt, France. E-mail [email protected], [email protected]

SUMMARY Generation expansion planning problems should take into account not only satisfaction of future demand but also environmental constraints, such as government imposed limits on pollution produced by coal based power plants. For a vertically integrated system of power plants, we model those features as a two-stage stochastic programming problem with recourse that can be solved by decomposition methods. We present results for different demand scenarios, corresponding to low, medium, and high economic future growth, as well as different environmental cases, with stringent and non-stringent c 2000 John Wiley & emission limits for the thermal plants composing the electric mix. Copyright Sons, Ltd. key words: Generation expansion problems, environmental constraints, two-stage stochastic programming problems

1. Introduction. Motivation The generation expansion planning (GEP) problem consist in determining an optimal investment plan for the construction of new power plants which ensures an economic and reliable supply to the future electricity demand. Generally, discounted expected costs of investment and operation are minimized subject to constraints depending on uncertain data, such as future growth of electricity demand, environmental restrictions, and primary fuel costs. Historically, GEP stems from centralized systems, with assets mostly belonging to the government. In spite of the worldwide trend of liberalization of energy markets, GEP tools did not loose relevance. On the contrary, they are now of interest also for private and institutional agents, acting on vertically integrated systems. For example, a minimum-cost GEP solution can help a private investor evaluating if a given plant can bring the desired return; see [6]. After Kyoto’s protocol was signed in 1997, limiting pollution became a crucial issue for the energy optimization area. For generation scheduling problems, pollution controls took the

∗ Correspondence

to: 1 IMPA Instituto de Matem´ atica Pura e Aplicada. Estrada Dona Castorina 110, Jardim Botˆ anico, Rio de Janeiro 22460-320, RJ, Brazil.

c 2000 John Wiley & Sons, Ltd. Copyright

2

´ R. M. MARCATO AND C. SAGASTIZABAL

form of green certificates, to be traded in a stock market. For expansion planning problems, pollution constraints of coal plants can be imposed on the operation of the coal mix by setting a limit to the maximum allowed for the total CO2/sulphur emission (which can be expressed as a function of the amount of coal burnt by each plant). The corresponding optimization problem is a multi-stage mixed-integer programming problem, solved by decomposition techniques; see [14], [15], [5, Ch. 10]. Regarding the Brazilian hydro-thermal power system, research in the subject started in the 80’s, evolving from the deterministic continuous model in [16] to the stochastic 0-1 model currently being used by ELETROBRAS [9]. The intermediate model [8], including hydrological uncertainties in full detail but too heavy from the computational point of view, should also be mentioned; see also [11]. In the international context, related works can be tracked back at least to [2], [4], [3]; and more recently, [17], [1], [12]; see also references therein. We present in this paper a model along the lines of [3] which takes into account different scenarios for pollution limits and for the future demand. As mentioned, we do not place the problem in a competitive economic environment, but rather in a vertically integrated system. We consider a portfolio of investments composed by candidate projects of two types: projects adding capacity to existent plants (with a concave piecewise linear cost curve for thermal plants) and projects for constructing new plants. The new coal plants can be standard, or using the so-called flue-gas desulphuration (fgd) technique to reduce pollution. Coal of both types (high and low sulphur), can be bought in limited quantities on the spot market in each period. In addition, long term contracts are available for ensuring the supply of low sulphur coal in large quantities along the planning horizon time of 10 years. However, these (takeor-pay) contracts involve the commitment of a minimum purchase in each period. Generation costs include operating and fuel costs. Load constraints are given by blocks (base, intermediate, peak), with a reserve margin in the peak block. The resulting model is a two-stage stochastic programming problem with recourse that we solve using a cutting-planes method. We assess our model by applying the methodology to a prototype system composed by nuclear, gas, oil, and coal plants. We present results for different demands scenarios (low, medium and high economic growth) and environmental cases (non-stringent and stringent). Our paper is organized as follows. In Sec. 2 we give the mathematical description of the optimization problem. Section 3 contains the solution method. Finally, we give our numerical experience in Sec. 4 and end with some concluding remarks in Sec. 5.

2. Formulation of the problem Our generation expansion problem minimizes the discounted sum of investments and operating costs of generating power in order to satisfy (uncertain) future electricity demand and (uncertain) environmental constraints, as well as (deterministic) operational constraints, such as maximum limit of generation for a given plant. In our model, decisions are taken every two-years, but there is annual information over a planning horizon of 2N years. Uncertainties are dealt with in two stages: the first stage, deterministic, corresponds to the first N1 periods, i.e., to the first 2N1 years, while the second stage, or recourse stage, consists of the remaining N − N1 periods, covering years 2N1 + 1 to 2N . c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

3

ENVIRONMENTAL CONSTRAINTS IN GENERATION EXPANSION PROBLEMS

Uncertainties regarding future demand and pollution limits are modelled in a scenario tree as in Figure 1, which shows that all uncertainties are revealed at the end of the N1th -period. From this period on, the optimization process selects optimal decisions for each one of the Q scenarios in the last N − N1 periods of the planning horizon. PSfrag replacements t=0 ... N1 N1 + 1 . . . N q=1 q=2 q=Q 1

Figure 1. Scenario tree representing uncertainties

Each branch q in the scenario tree has a probability denoted by pq , corresponding to the joint probability of the scenario uncertainties. For example, if scenario q = 1 corresponds to the case of low economic growth and stringent pollution control, each one with probability 0.3, then p1 = 0.09. 2.1. Objective function In our generation expansion problem, in addition to generation, operation, and fuel costs, there are the cost of expanding the mix and the cost of buying low sulphur coal with two long term contracts. Both low and high sulphur coal can be bought in the spot market, at a variable price. Contracts offer the possibility of a large supply of nonpolluting coal, at a fixed price. The expansion portfolio has projects indexed in a set M , with the coal plants ranging in the subset I. Projects include plants that are already built but can have their capacity expanded, as well as candidate plants, not yet built. The variables cap, cont, spot, and gen below denote, respectively, added capacity, low sulphur coal bought with a contract (contract 1 or contract 2), coal bought in the spot market (low or high sulphur coal, denoted by 1 and 2, respectively), and generation (at ν = base, peak, and intermediate demand) of a given project at a given time period. We also use modelling binary variables, x and y. For each variable, a subindex (i) specifies the project and a superindex (t) the time period. For second stage variables, a second subindex (q) indicates the corresponding scenario. Finally, in order to give an abstract formulation of the problem, we gather first stage variables in the vector p, and second stage variables in the vector p˜, as follows:
pti = x1 ti , cap1 ti , x2 ti , cap2 ti contt1 , contt2 , spot1 ti , spot2 ti , y1 , y2 , {genν ti }3ν=1 for i ∈ M and t ≤ N1 , and p˜ti,q = x1,q ti , cap1,q ti , x2,q ti , cap2,q ti contt1,q , contt2,q , spot1,q ti , spot2,q ti , {genν,q ti }3ν=1




for i ∈ M , q = 1, . . . , Q, and t ≥ N1 + 1. c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

´ R. M. MARCATO AND C. SAGASTIZABAL

4

2.1.1. Expansion cost Typically, thermal expansion costs depend on the size of the considered project, as well as on the technology applied (for instance, combined cycle and open cycle plants have different costs). There are many different ways of modelling expansion costs. The range goes from simple piecewise concave costs (as considered in this paper) to sophisticated cost assignment models based on game theory (as in [7]). The simple model we present here is suitable for a centralized power system, like Brazil’s, with a high predominance of hydro-plants in the mix. In our model, projects of the expansion portfolio that correspond to coal based plants have a piecewise concave cost for capacity additions, as shown in Figure 2. As shown in this figure, the values fi , bi , ci , and di vary with each project i. Such variation allows us to consider different sizes and technologies for different projects.

PSfrag replacements GW 106 $

0 fi bi ci di Figure 2. Capacity expansion cost for ith thermal project

We start with first stage expansion costs. The formulation of the cost concavity requires the introduction of auxiliary binary variables. More precisely, we set x1 ti ∈ {0, 1} to 1 if the ith project is built at time period t, while x1 ti = 0 otherwise. Similarly, x2 ti ∈ {0, 1} is 1 if the addition of capacity is bigger than ci , and is 0 otherwise. With this notation, if the respective added capacities are cap1 ti ∈ [bi , ci ] and cap2 ti ∈ [ci , di ], the cost of adding a total capacity cap1 ti + cap2 ti for the ith project at time period t is equal to Eit (pti ) = (fi + s1i bi )x1 ti + s1i (cap1 ti − bi x1 ti ) + s2i (cap2 ti − bi x2 ti ),

(1)

where s1i ≥ s2i correspond to the slopes in Figure 2 and the variables must satisfy the following constraints x2 ti ≤ x1 ti ,

bi x1 ti ≤ cap1 ti ≤ ci x1 ti ,

and ci x2 ti ≤ cap2 ti ≤ di x2 ti

(2)

for all i ∈ I and t ≤ N1 . Second stage expansion costs are formulated similarly, but using an extra subindex to identify the scenario. More precisely, for all i ∈ I, q ≤ Q and t ≥ N1 + 1, Eq ti (˜ pti,q ) = (fi + s1i bi )x1,q ti + s1i (cap1,q ti − bi x1,q ti ) + s2i (cap2,q ti − bi x2,q ti ),

(3)

with variables x1,q ti , x2,q ti , cap1,q ti , and cap1,q ti satisfying a relation like in (2). c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

ENVIRONMENTAL CONSTRAINTS IN GENERATION EXPANSION PROBLEMS

5

2.1.2. Coal costs For the two long term contracts supplying low sulphur coal in large quantities for the whole thermal mix, the variable contt1 (contt2 ) represents the consumption of low sulphur coal at time period t that was bought via the first (second) contract, expressed in millions of BTU. In addition, each thermal plant can buy coal in the spot market, in the amounts spott1,i and spott2,i corresponding, respectively, to high and low sulphur coal bought by plant i at time period t. The total coal cost of the ith thermal plant at time period t ≤ N1 is given by Cit (pti ) = α1 contt1 + α2 contt2 + α3 spott2,i + α4 spott2,i ,

(4)

where α corresponds to unit costs, possibly including conversion and scaling factors. For t ≥ N1 + 1, the cost depends on each scenario q ≤ Q: Cq ti (˜ pti,q ) = α1 contt1,q + α2 contt2,q + α3 spott2,i,q + α4 spott2,i,q .

(5)

The two contracts are modelled with binary variables: y1 ∈ {0, 1} is 1 if the first contract is purchased, and is 0 otherwise. Likewise for the second contract, with variable y2 . The decision of contracting or not is made at the first time period (this is the reason why the variables y1 , y2 do not depend on t), but contracts impose minimum and maximum purchased amounts at each period: min1 y1 ≤ contt1 ≤ max1 y1 and likewise for second stage variables

and min2 y2 ≤ contt2 ≤ max2 y2 ,

contt1,q

and

(6)

contt2,q .

2.2. Operation costs Other costs are those corresponding to operation and fuel costs, both for plants already existing in the mix, and for plants that are built along the planning horizon. Consider the ith project and denote by opi its operating cost (expressed in 104 $); by f ueli its fuel cost (in 106 $), estimated to increase a 2% per year; by genν ti , its generation at time period t, during the ν th demand component, corresponding to ν = base, peak, and intermediate demand (in GWh). The resulting cost is 3 X Oit (pti ) = opi + 1.022(t−1) f ueli genν ti (7) ν=1

for i ∈ M and t ≤ N1 , and Oq ti (˜ pti,q ) = opi + 1.022(t−1) f ueli

3 X

genν,q ti

(8)

ν=1

for i ∈ M , t ≥ N1 + 1 and q ≤ Q. The total cost is given by the addition of all the first stage costs and the expected value of second stage costs above, discounted by an annual rate τ . In an abstract manner, our objective function is given by fo (p) + f˜o (˜ p) where we defined N1 X

 X 1 t t t t t t E (p ) + C (p ) + O (p ) , i i i i (1 + τ )2(t−1) i∈M i i t=1 Q N  X X X 1 t t t t t t f˜o (˜ p) = p E (p ) + C (p ) + O (p ) , q q q q i i,q i i,q i i,q (1 + τ )2(t−1) q=1 i∈M t=N +1 fo (p) =

(9)

1

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

´ R. M. MARCATO AND C. SAGASTIZABAL

6

and expansion, coal, and operation costs are as defined in (1)-(3), (4)-(5), and (7)-(8), respectively. Coal costs (4)-(5) are defined only for i ∈ I, we defined Cit = Cq ti = 0 for i ∈ M \I and all t to simplify the summation in (9). 2.3. Problem constraints Except for the binary variables, which range in the interval [0, 1], and the purchased coal in the spot market, which range in [0, spottmax ], variables vary in [0, +∞). We express such constraints in the condensed form p ∈ Pbox and p˜ ∈ P˜box , (10) ˜ where Pbox and Pbox are polyhedral sets of appropriate dimensions. In addition to the above box constraints, there are general equality and inequality constraints. Among these constraints, some of them correspond to relations that couple first and second stages. We start with equality constraints. First, each demand component ν must be satisfied at every time period t: X X gerν ti = demtν and gerν,q ti = demtν,q for each scenario q. (11) i∈M

i∈M

For thermal plants, at every time period burnt coal equals generation, via a conversion factor βi . Burnt coal can be either high sulphur (spott1,i ), or low sulphur (spott2,i ) bought in the spot market, or low sulphur supplied by one of the contracts. Introducing the auxiliary variable lowit for the amount of low sulphur coal supplied to plant i at time period t by both contracts yields the relations: βi

3 X

genν ti = lowit + spott1,i

X

and

ν=1

lowit + spott2,i + contt1 + contt2 = 0.

(12)

i∈M

Similar equality constraints appear for second stage variables, mutatis mutandis. Since all these constraints are affine and do not couple stages, we gather (11) and (12) (with its corresponding variant for second stage variables) in the polyhedral sets p ∈ Paff= and p˜ ∈ P˜aff= . (13) Suppose some project ic is compulsory, i.e., it has to be built during the planning horizon. For each scenario q ≤ Q the corresponding modelling binary variables must satisfy the relation N1 X

N X

x1 tic +

t=1

x1,q tic = 1,

(14)

t=N1 +1

coupling first and second stage variables. Even though these are just linear relations, for more generality we write coupling equality constraints like (14) in the form cE (p) + cE (˜ p) = 0.

(15)

We now consider inequality constraints, starting with those that do not couple stages. All capacity binary variables satisfy (2). For t ≤ N1 , contract variables satisfy (6). Also, for t ≤ N1 , generation cannot be bigger than the current capacity of the plant: 3 X

γν i gerν ti ≤ cap0i +

ν=1

c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

t−1 X

(x1 ti + cap1 ti + cap2 ti ) ,

(16)

j=1

Numer. Linear Algebra Appl. 2000; 00:1–6

ENVIRONMENTAL CONSTRAINTS IN GENERATION EXPANSION PROBLEMS

7

where γ is a conversion factor and cap0i is the initial capacity of the plant. Pollution constraints are expressed in the form of maximum limits of the total sulphur emission allowed for the coal mix for each environmental scenario, we denote such limits by polmax q . Clearly, the emission produced when using low sulphur coal (denoted by `) is less than when using high sulphur coal (~). Moreover, for thermal plants with fgd, both levels of emission are lower, we denote them by `fgd and ~fgd , respectively. As a result, for all t ≥ N1 + 1 and q ≤ Q, environmental constraints have the form X X X X spot1,q ti + `fgd ` lowq ti + ~ lowq ti + ~fgd spot1,q ti ≤ polmax q . (17) i∈Inotfgd

i∈Inotfgd

i∈Ifgd

i∈Ifgd

Constraints (16) and (17) are affine and not coupling, we again express them using polyhedral sets as follows: p ∈ Paff≤ and p˜ ∈ P˜aff≤ . (18) For not compulsory projects, capacity binary variables must satisfy a relation like (14) with inequality, coupling first and second stage variables. The same holds for contract inequalities (6) when t ≥ N1 + 1. Similarly, for t ≥ N1 + 1, the maximum generation constraint (16) has the form P3 PN1 t 0 + j=1 (x t + cap1 ti + cap2 ti ) ν=1 γν,q i gerν,q i ≤ capi Pt−1 1 i (19) + j=N1 +1 (x1,q ti + cap1,q ti + cap2,q ti ) . We write constraints (14) with inequality, (6) for second stage variables, and (19) in the general form cI (p) + cI (˜ p) ≤ 0. (20) Putting together equations (9), (18), (15), and (20), we obtain the following optimization abstract problem:  min fo (p) + f˜o (˜ p)      p ∈ P := lin(Pbox ) ∩ Paff= ∩ Paff≤ (21) p˜ ∈ P˜ := lin(P˜box ) ∩ P˜aff= ∩ P˜aff≤    c (p) + c˜E (˜ p) = 0   E cI (p) + c˜I (˜ p) ≤ 0 . In the feasible set P × P˜ above, we use the sets in (13) and (18); as for the box constraints, instead of the sets in (10), we use their linear relaxation, obtained when replacing each 0-1 variable by a continuous variable ranging in the interval [0, 1]. We make this approximation of the original mixed-integer linear programming problem in order to ensure convergence of the nonsmooth optimization method applied for solving (21).

3. Solution method Problem (21) has a structure suitable for the application of a decomposition method. Specifically, we identify in (21) a partition of variables; namely, easy ones, corresponding to first stage variables p, and complicating variables p˜, for second stage decisions. A similar partition can be identified in the objective function and the constraints of (21). Furthermore, constraints involving second stage variables (constraints (14), (6) for t ≥ N1 + 1, and (19)) do not couple c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

´ R. M. MARCATO AND C. SAGASTIZABAL

8

different scenarios. We observe, in addition, that constraints present a strong coupling along time steps. For these reasons, rather than a method based on decomposing the problem at each time stage, we apply a Benders’-like decomposition, as described below. Consider the following bi-level problem, equivalent to (21):  min fo (p) + v opt (p) (22) p∈P where we defined the optimality function  min f˜o (˜ p)    p˜ ∈ P˜ opt v (p) =  c˜ (˜ p) = −cE (p)   E c˜I (˜ p) ≤ −cI (p) .

(23)

For our relaxed problem (21), where all the data is convex, v opt is a convex nondifferentiable function with separable structure (by scenarios). The reformulation of (21) using two levels yields a decomposition method, composed by a master program (22), and a slave problem (23), decomposable in scenarios. At each iteration, the master program sends a value pk to the Q slave problems. The solution of the slave problems gives the value of v opt (pk ) that returns to the master where a new point pk+1 is chosen. Since the master program is a constrained nonsmooth optimization problem, pk+1 will be chosen by using a nondifferentiable black-box method ([5, Ch. 9]), i.e., using the objective value fo (pk ) + v opt (pk ), and a subgradient. To compute the subgradient information for the optimality function, we use the primal-dual information obtained when solving (23), as shown in the following lemma (for completeness, we include a proof of this result in Appendix II). Lemma 3.1. Consider (23) written with p = pk , and suppose the constraints satisfy Slater ˜ k , the following holds: condition. Letting p˜k be a solution with associated Lagrange multiplier λ v opt (pk ) = f˜o (˜ pk )

and

˜ k ∈ ∂v opt (pk ), s(pk ) = Jc(pk )T λ

(24)

where we used the transposed Jacobian matrix Jc(pk ) = (JcE (pk ), JcI (pk )). Therefore, at the coordination step, i.e., when the master program chooses pk+1 , both v (pk ) and the subgradient s(pk ) are available after solving the second stage problem (23). The decomposition technique in two levels that uses a cutting-planes method for solving (22) is called Benders decomposition when the original problem is linear, see [Ch. 10][5]. For our problem (21), fo (p) = ha, pi for some vector a and P is polyhedral. At each iteration, the objective function in (22) is then approximated by opt

ha, pi + vˇkopt (p), where vˇkopt is a cutting-planes model for v opt (·), built with the black-box information accumulated so far:  vˇkopt (p) = maxj≤k nv opt (pj ) +D hs(pj ), p − pj i Eo ˜ j , p − pj = maxj≤k f˜o (˜ pj ) + Jc(pj )T λ . c 2000 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls

Numer. Linear Algebra Appl. 2000; 00:1–6

ENVIRONMENTAL CONSTRAINTS IN GENERATION EXPANSION PROBLEMS

9

A straightforward application of the cutting-planes method would give for the master program a linear programming problem (LP) of the form:    min ha, pi + r min ha, pi + vˇkopt (p) p∈P,r∈< D E ⇔ ˜ j , p − pj for j ≤ k. p∈P  r ≥ fo (˜ pj ) + Jc(pj )T λ However, at this point, it is important to notice that, depending on pk , the feasible set in (23), corresponding to second stage decisions, i.e. n o p˜ ∈ P˜ : c˜E (˜ p) = −cE (pk ) and c˜I (˜ p) ≤ −cI (pk ) , may be empty. In this case, the value of v opt (pk ) and a subgradient cannot be computed. Instead, a new pk+1 , for which problem (23) written with p = pk+1 is feasible, should be generated by the master program. The new iterate is computed by introducing in the master LP a feasibility cut. The idea is, before trying to compute v opt (pk ), to define the feasibility function   min kzk1   ˜ + − z+ ≥ 0   I  p˜ ∈ P , zE +, zE , − c˜E (˜ p) + zE − zE = −cE (pk ) (25) v f eas (pk ) =  c˜I (˜ p) + zI+ ≤ −cI (pk )    + − +  − zE , zI ) ∈