Optimal planning of cogeneration production with provision of ancillary services

Electric Power Systems Research 95 (2013) 47–55

Contents lists available at SciVerse ScienceDirect

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

Optimal planning of cogeneration production with provision of ancillary services ˇ Petr Havel a,∗ , Tomáˇs Simoviˇ cb a b

Department of Control Engineering, Faculty of Electrical Engineering, Czech Technical University, Technická 2, Praha, Czech Republic Department of Optimization, Simulation, Risk and Statistics (OSIRIS), EDF R&D, Paris, France

a r t i c l e

i n f o

Article history: Received 2 November 2011 Received in revised form 10 June 2012 Accepted 5 July 2012 Available online 29 September 2012 Keywords: Power generation planning Cogeneration Ancillary services Combined heat and power Optimal scheduling

a b s t r a c t Optimal planning of a cogeneration plant operation for a day-ahead horizon with emphasis on provision of ancillary services is presented in this paper. A first-principle mathematical model is used to formulate an optimization problem in the form of a mixed integer and linear program (MILP). The resulting decision support tool recommends the plant’s operator the best strategy how to trade electricity, provide contracted ancillary services and operate the plant in each hour in order to maximize profit. Thanks to the cooperation with the company Taures, a.s., this tool has been in operation at two large cogeneration plants for more than two years. © 2012 Elsevier B.V. All rights reserved.

1. Introduction Large combined heat and power (CHP) plants (cogeneration plants) are complex systems capable of providing heat for a district system or industrial needs, electrical energy and ancillary services (AS) for the transmission system. Efficient operation of a CHP plant in a liberalized market environment requires effective decision support tools capable of determining a cost effective operation of a CHP plant to meet time varying demands of heat, electrical energy and ancillary services (regulation reserves). This paper describes an industrial application of such a tool developed for and currently used at two CHP plants in the Central European region for a day-ahead production planning and electricity trading. Five different demands are considered (heat, electrical energy and three distinct types of AS: primary, secondary and tertiary reserve [1]) along with the possibility to trade additional electrical energy on the day-ahead market. With electric market deregulation and the increasing introduction of renewable sources of energy, the provision of AS can be a significant source of revenues for operators of larger CHP plants and therefore needs to be considered in production planning [2]. Hence one of the main contributions of this paper is the inclusion of the provision of AS in a CHP production planning model. The paper is organized as follows: after a brief survey of relevant literature, the proposed model of a CHP plant is described in Section 3. Section 4 presents the inclusion of ancillary services into

∗ Corresponding author. Tel.: +420 224 357 343; fax: +420 224 916 648. E-mail address: [email protected] (P. Havel). 0378-7796/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.epsr.2012.07.020

the model. Finally, Section 5 describes an optimization problem formulation and in Section 6 some remarks on implementation and sample optimization results are presented. 2. A short survey on CHP unit commitment The problem addressed in this paper is similar to the one known in literature as a unit commitment1 (UC) problem where in our case UC of boilers, turbines and other parts of the technology in a single plant is considered. A good general survey of UC problems is given in [3] while [4] is specifically focused on short term cogeneration planning. 2.1. Modeling approaches A formulation of the UC problem for a CHP plant requires a model that describes the operation of the system. Two distinct approaches can be used: a white box (first principles) approach using balance equations and black box approach using a feasible operating region. A feasible operating region (Fig. 1) can be represented by a polytope defined by coordinates (p, q, c) where (p, q) is a working point of the turbine yielding p MW of electricity and q MW of heat, and c is the cost function representing either directly fuel costs or an amount of steam necessary for operation at the

1 In our survey, in line with the terminology used in [3], if a problem considered in a given publication does not consider the on/off states of equipment and the related features (start up constraints, minimum up and down times, etc.) we label it as “economic dispatch” (ED) problem. We consider an economic dispatch problem as a sub problem of the UC problem.

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

48

Nomenclature Variable Description [units] cjstart costs for a single start-up of the unit j [EUR] stop cj

costs for a single shut-down of the unit j [EUR] fuel price per MWh of combustion heat (including price of CO2 allowances) [EUR/MWh] Cfuel (t) total fuel costs at a time sample t [EUR] cshort (t) expected day-ahead price of electrical energy for the hour t [EUR/MWh] Csusd (t) total cost for start-ups and shut-downs of all units at time t [EUR] hIN boiler feedwater enthalpy [kJ/kg] boiler hOUT boiler output steam enthalpy [kJ/kg] boiler enthalpy of condensate outflow from a heat hOUT cond exchanger [kJ/kg] hcool enthalpy of cooling water injected to a reduction station [kJ/kg] hi steam enthalpy at the extraction point i of a turbine [kJ/kg] junct hi enthalpy of mass flow to the junction i [kJ/kg] hIN enthalpy of steam inflow to a reduction station RS [kJ/kg] hOUT enthalpy of steam outflow from a reduction station RS [kJ/kg] enthalpy of steam inflow to a heat exchanger [kJ/kg] hIN steam hIN enthalpy of water inflow to a heat exchanger [kJ/kg] water hOUT enthalpy of water outflow from a heat exchanger water [kJ/kg] boiler feedwater flow [t/h] Mboiler OUT Mcond condensate outflow from a heat exchanger [t/h] flow of cooling water injected to a reduction station Mcool [t/h] junct Mi mass flow to the junction [t/h] Mj steam flow through the jth extraction point of a turbine [t/h] IN steam inflow to a heat exchanger [t/h] Msteam IN MRS steam mass inflow to a reduction station [t/h] OUT MRS steam mass outflow from a reduction station [t/h] IN water inflow to a heat exchanger [t/h] Mwater OUT Mwater water outflow from a heat exchanger [t/h] Paux auxiliary power consumption [MWe] IN Pboiler fuel consumption of a boiler [MWt] IN Pboiler,op boiler fuel consumption at a base operating point [MWt] IN Pboiler,Si boiler fuel consumption when the ancillary service i is deployed [MWt] OUT Pboiler thermal output of a boiler [MWt] OUT Pboiler,min boiler power output at a minimal plant electric output [MWt] OUT Pboiler,max boiler power output at a maximal plant electric output [MWt] OUT Pboiler,op boiler power output at a base operating point [MWt] cfuel

OUT Pboiler,Si

Pcor pj (t) Pjmax

boiler power output when the ancillary service i is deployed [MWt] correction of a turbine electric output [MWe] output of the unit j in the discrete time step t [MW] technologically maximal power output of the unit j [MW]

PkIN

Rjdown

boiler fuel consumption at a characteristic point k [MWt] boiler power output at a characteristic point k [MWt] volume of energy contracted on a longer horizon (prior to the day-ahead optimization) [MWe] total maximal possible electricity supply [MWe] auxiliary power consumption at the plant maximal electric output [MWe] maximal possible electricity supply from turbine TG [MWe] total minimal possible electricity supply [MWe] auxiliary power consumption at the plant minimal electric output [MWe] minimal possible electricity supply from turbine TG [MWe] auxiliary power consumption at the base operating point [MWe] plant base operating point when no ancillary services is activated [MWe] turbine electric output at the base operating point [MWe] plant operating point when ancillary service i is fully deployed [MWe] auxiliary power consumption when the ancillary service i is deployed [MWe] total contracted volume of the ancillary service i [MWe] volume of the ancillary service i provided at a turbine TG [MWe] recommended optimal volume of electrical energy to be traded on the day-ahead market [MWe] turbine electric output [MWe] demanded thermal power for district heating [MWt] ramp-down limit for a unit j [MW/time sample]

Rj

ramp-up limit for a unit j [MW/time sample]

down rboiler up rboiler Sjdown

ramp-down limit of a boiler [MW/min] ramp-up limit of a boiler [MW/min]

PkOUT c Plong

Pmax aux Pmax TG Pmax

Pmin aux Pmin TG Pmin aux Pop c Pop TG Pop

PSi aux PSi c PSi TG PSi c Pshort

PTG Qdemand

up

up

Sj T

Tjdown up

Tj

TSi ub uTG

vj (t) wj (t)

power limit for a transition from on-state to offstate [MW] power limit for a transition from off-state to onstate [MW] length of the optimization horizon [number of time samples] minimum down time of a unit j [number of time samples] minimum up time of a unit j [number of time samples] required deployment time of the ancillary service i [min] unit commitment of a boiler (binary decision variable) [–] unit commitment of a turbine (binary decision variable) [–] binary variable indicating start-up of the unit j in time t [–] binary variable indicating shut-down of the unit j in time t [–]

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

Fig. 1. Feasible convex operating region of a CHP plant [4].

working point (p, q) [5–8]. For more complex plants this approach cannot guarantee the practicability of transitions between adjacent points of the feasible operating region (e.g. when AS are activated), as these points can be potentially achieved with radically different equipment configurations. This representation cannot be therefore used for all types of CHP plants. Besides, it is applicable only if sufficient measured historical data across the entire feasible region is available. Mass and energy balance equations can be used to model more complicated systems but result in a more complex steady state mathematical model. Non linear representation has been used in a number of papers [8–12] to model the input/output enthalpy function of boilers and turbines. On the other hand, a linear [13] or piece wise linear [14,15] dependence seemed to be sufficient to model equipment in other publications. A simple dynamic model has been proposed in [16] based on the framework of mixed logic dynamical systems (MLD). However, this framework and dynamic models in general, seem to better suit to online optimization and control applications which is not the scope of this paper. 2.2. Solution methods A range of methods for finding the optimal UC solutions can be found in the literature, the main ones together with their features are listed below. An advantage of dynamic programming (DP) in solving CHP problems is that it is an exact method that allows for non-linear characteristics of technological equipment. DP was used in [11] to solve an economic dispatch (ED) problem for a cogeneration plant including heat storage and ramping constraints, with sub problems being solved by a commercial non linear solver, based on sequential quadratic programming (SQP). An SQP method customized for the CHP ED problem is used in [8]. Starting from an initial solution, inequality constraints in the problem are relaxed and iteratively transformed to equality constraints to solve an easier problem. In [12] a Newton method is used to solve the ED problem of a simple CHP model. In recent years, an increasing number of publications use a MILP formulation solved by Branch & Bound method. In [14] a customized B&B method is used to solve a profit maximization UC problem considering time dependent start-up costs with significant improvement reported over a standard B&B method. In [5] a customized B&B method to solve a MILP subproblem used for a long term CHP production planning problem is considered. Commercial B&B software is sufficient to solve the profit maximization problem presented in [13] and [15]. Lagrangian relaxation (LR) methods have been applied to problems that are too large to be effectively solved by a B&B approach. In [6] the LR method was used to

49

decompose the CHP planning problem into a heat dispatch and a power dispatch that are performed iteratively. In [7] a LR approach led to better results for longer planning horizons and more complicated heat diagrams compared to a commercial B&B solver. Heuristic or metaheuristic methods seem to be adapted in cases where equipment has to be modeled using 3rd or 4th order polynomials ([9] and [10], for example). Lately, most publications on the CHP planning problem include also electricity market in various forms. Ancillary services, however, have not yet been dealt with in a CHP production planning problem in comprehensive manner to our knowledge. In [7] secondary reserve is considered but is committed using a priority list and not with a global optimization algorithm. In [15] ancillary services are given as a percentage of turbine power output, no other constraints seem to be taken into account. The authors of [14] consider an aggregated reserve and note that AS provision may change unit commitment. However, the primary and secondary reserve is not set apart and no other constraints besides reserve demand are considered. In [2] it is noted that provision of ancillary services can change the operational configuration of a CHP plant in a significant manner compared to operation without ancillary services. Nevertheless, the focus of [2] is neither the unit commitment problem nor an industrial decision-aid tool. 3. Model of CHP plant The choice of a solution method was guided by the necessity to deliver an application to be routinely used in the industry. The complexity of the studied CHP plant neither allows the use of the explicit feasible operating P–Q region approach nor simple blackbox models mentioned in the survey. Therefore, a first-principle modeling approach based on mass and energy balance equations in the closed steam cycle was chosen to model the CHP technology and operational costs. In this approach the P–Q region is expressed implicitly by the set of equations and technical limits of individual technological units in the CHP plant, mainly available power ranges of boilers. Such an approach makes it possible to model a condensing extraction turbine with multiple extractions. Another advantage is the possibility to model operational states that have not been used by the plant operators before and are thus not contained in the historical measurements, although they can be currently economically beneficial due to changing market opportunities or fuel price variations. Consequently, a MILP formulation of the problem, combined with a general purpose high-performance commercial solver was selected for two main reasons: • Compared to the past, very efficient MILP solvers are available today. • The optimality gap indicator of the Branch & Bound method indicates the “quality” of the solution obtained and a convergence to the global optimum is guaranteed. This is important if the tool is used to make financially sensitive decisions. A nonlinear physical behavior of some of the plant components is incorporated into the MILP model by two means: • The nonlinear characteristic is directly approximated by a linear characteristic with an acceptable error. It is assumed that the steam/water parameters, expressed by the enthalpy value, at inputs/outputs of the CHP plant technological components are constant or at least do not change significantly.

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

50

• A piece-wise linear approximation is used with as many characteristic points as necessary to reach an acceptable error. However, this approach increases the optimization problem size and computational time. In the following text, models of CHP plant components will be described together with additional equations that represent other constraints of the resulting optimization problem (ramping rates, limits on fuel consumption, ranges of regulation reserves for AS provisioning etc.). The equations in the Sections 3.1–3.3 represent constraints for each modeled time step (e.g. an hour, a quarter of hour) while the equations in Sections 3.4 and 3.5 represent coupling between the time steps. 3.1. Models of technological components Models of basic CHP plant components are described by their input/output terminals and associated values of steam/water enthalpies. 3.1.1. Boiler Each steam boiler is characterized by a user-defined piecewise linear curve representing the boiler efficiency, i.e. the boiler outOUT (power in MW effectively delivered to the feedwater to put Pboiler IN generate steam) vs. boiler input Pboiler (fuel consumption in MW). In the optimization, the boiler is modeled with a special order set (SOS) of type 2 constraints: OUT = Pboiler



IN Pboiler =



k PkOUT k PkIN

(1a)

k

k − ub = 0

k

k ≥ 0, k ∈ SOS 2 where k form SOS of type 2, which means that not more than two adjacent k may take non-zero values. The variables PkOUT and PkIN represent measured or known characteristic points of the boiler operation (the first one P1OUT represents minimal boiler output, the last one PNOUT maximal boiler output). The variable ub is a binary decision variable (ub = 1: the boiler is on, ub = 0: the boiler is off). The relation between the boiler power output and the steam generation is OUT = Mboiler (hOUT − hIN ) Pboiler boiler boiler

(1b)

where Mboiler is the feedwater inflow (equal to the steam outflow is the output steam enthalpy and if dynamics is neglected), hOUT boiler hIN is the feedwater enthalpy. boiler 3.1.2. Turbine As a general case, an extraction condensing turbine is assumed. It has a potential of flexible power generation not so tightly linked with the heat supply which enables provision of AS. The turbine may have extraction points at different pressures for different purposes, e.g. a supply of steam for district heating or industrial use or steam for boiler water reheating circuit. The electrical output of an extraction condensing turbogenerator (Fig. 2) can be approximated as P TG =

i=0

TG

steam boiler

h1

M1

Mn+1,h n+1 condenser

h2 … hn Mn M2

Fig. 2. A schematic representation of a condensing extraction turbine with multiple extractions.

the extraction point i and n is a number of extraction points. The enthalpy hn+1 is the enthalpy of the last steam flow Mn+1 going to the condenser. The term Pcor represents a correction in the form Pcor = acor M0 + bcor uTG which accounts for mechanical losses and other non-modeled phenomena. The coefficients acor , bcor can be either estimated or tuned on real historical measurements utilizing, e.g., a least-square fit function, the variable uTG is a binary decision variable determining the turbine state (on/off). Besides the energy balance (2a) the mass balance must hold true M0 =

n+1 

Mj .

(2b)

j=1

k 

n 

M0,h 0

⎛

⎝M0 −

i 

⎞

Mj ⎠ (hi − hi+1 ) − Pcor ,

(2a)

j=1

where M0 is an input steam mass flow to the turbine, Mj the steam flow through the jth extraction point, hi steam enthalpy in

The chosen modeling approach (2a) enables to linearly model almost any turbine with arbitrary number of extraction points under the assumption that the enthalpies hi are constant or at least do not depend significantly on the steam flow. If the enthalpy dependency on the steam flow cannot be neglected (thermodynamic efficiency is variable) the turbine or each of its sections has to be modeled by a piece-wise linear model similarly to the boilers model (1a). This, of course, increases optimization complexity and solve time. 3.1.3. Steam pressure reduction and cooling station A reduction and cooling station is typically used to bypass a turbine when there is a high heat demand and simultaneously a low requirement for electricity generation, e.g. in case of negative AS activation (downward regulation). The reduction and cooling station is modeled by the following energy and mass balance constraints IN IN OUT OUT hRS + Mcool hcool = MRS hRS MRS

(3a)

OUT IN = MRS + Mcool MRS

(3b)

IN where MRS

is steam mass inflow, hIN RS

its enthalpy, Mcool is flow of OUT steam outflow and injected cooling water, hcool its enthalpy, MRS OUT hRS output steam enthalpy. 3.1.4. Heat exchanger Heat exchangers (Fig. 3) are typically used for two reasons: to preheat feedwater for boilers or to heat water for the district heating. A mass and energy balance model of the heat exchanger is IN IN IN OUT OUT OUT OUT Msteam hIN steam + Mwater hwater = Mcond hcond + Mwater hwater

(4a)

IN OUT Msteam = Mcond

(4b)

IN OUT Mwater = Mwater

(4c)

where the upper index IN denotes inflow, OUT outflow and the lower index cond denotes the steam condensate. A heat exchanger may have more than one input of water, e.g. when condensates

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

where ux are selected binary variables representing on/off states of devices that have significant influence on the auxiliary power consumption and kx are weight coefficients determined from least-squares fitting on real data of auxiliary power consumption. Similarly, the variables vy are selected continuous variables (e.g. flows to the cooling unit, flows through the steam boilers, flow of hot water to the district) and ky the respective coefficients.

IN IN M steam , hsteam

steam IN IN M water , hwa ter

OUT OUT M water , hwa ter

cold water

hot water

condensate

51

OUT OUT M con d , hcond

3.4. Ramping constraints

Fig. 3. Schematic diagram of a heat exchanger.

from more devices are collected, which can be modeled by inserting more terms to left-hand side of (4a) and (4c). The exchanger used for district heat supply in form of hot water may be modeled by modification of (4a) IN OUT OUT hIN Msteam steam − Mcond hcond = Qdemand

Ramping constraints limit changes on production of units between the time steps. Such a unit can be, for instance, a steam boiler which is a typical bottleneck for flexible electricity generation. The ramping constrains are defined according to [17]: pj (t) − pj (t − 1)

(5)

up

+Sj (uj (t) − uj (t − 1))

where Qdemand is the heat demand in MW. pj (t − 1) − pj (t)

The basic components of the plant (boilers, turbines, heat exchangers etc.) are connected to each other by junctions to form a steam cycle of a CHP plant (Fig. 4). Each junction is described by a balance equation junct junct hi Mi

=0

(6)

i junct

where Mi

are mass flows to the junction, with positive sign junct

for inflows and negative sign for outflows, and hi corresponding to the flow

(8)

+Pjmax (1 − uj (t))

3.2. Steam cycle model



up

≤ Rj uj (t − 1)

≤ Rjdown uj (t) +Sjdown (uj (t − 1) − uj (t)) +Pjmax (1 − uj (t

(9)

− 1))

where pj (t) is an output of the unit j in the discrete time step t, up Pjmax is the maximum power output of the unit, Rj and Rjdown are allowed maximal outputs changes between two time steps (rampup up and ramp-down limits) and Sj , Sjdown are limits for a transition from off-state (zero production) to on-state and vice versa.

is the enthalpy

junct Mi .

3.5. Minimum up and down times and start-up and shut-down costs

3.3. Auxiliary power consumption The auxiliary power consumption includes electrical power supplied to the water pumps, generator circuits and other pieces of equipment. Its value Paux is modeled as a linear model

In order to avoid undesired fast switching of boilers and turbines in time, constraints modeling minimum up and down times and start-up and shut-down costs are employed. Their formulation for each unit j was implemented based on recommendations in [18]

P aux =

vj (t) − wj (t) = uj (t)



kx ux +

x



ky vy

(7)

y

t 

vj (q) ≤ uj (t),

(10)

∀t ∈





up

Tj , . . . , T

(11)

up +1 j

q=T −T

steam PK boiler

t 

wj (q) ≤ 1 − uj (t),

turbine TG1



∀t ∈ Tjdown , . . . , T



(12)

q=T −T down +1 j

reducon and cooling staon

R

condenser heat exchangers

district heating system

Csusd (t) =



stop

(cjstart vj (t) + cj

wj (t))

(13)

j

where vj (t) and wj (t) are binary variables indicating start-ups and up shut-downs of the unit j in time t, Tj and Tjdown are minimum up and minimum down times, T is a length of the optimization stop horizon (CHP operation planning horizon), cjstart and cj are costs for a single start-up and shut-down of the unit j. Csusd (t) are total cost for start-ups and shut-downs of all units at time t. 4. Provision of ancillary services

feedwater tank Fig. 4. Example of a cogeneration steam cycle.

Although it differs from country to country and market from market [19], in our approach it is assumed that a major volume of AS

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

52

is traded on longer horizons2 (typically in yearly or monthly “baseload” contracts) and only an energy market is interesting and has high liquidity in a day-ahead horizon. Hence, the optimizer should find such a solution that trades additional energy on the day-ahead market and suggests a technically feasible generation schedule of the energy and contracted AS while maximizing profit. 4.1. Available regulation range

condensing regime is applicable in the plant. Minimal possible generation is often limited by minimal outputs of boilers or capacity of reduction stations that provide turbine steam bypass. Of course, both cases are also limited by the heat demand. From the dynamic point of view, provision of AS and available regulation range is typically limited by the rate of change of steam production. This is included in the optimization problem as follows (for each of the boilers): up

Provision of contracted AS means that the CHP plant is operated in such a way, that in each time period (e.g. an hour) it has enough generation capacity and speed to change the power output on request according to specifications of different contracted AS types such as primary, secondary and tertiary control reserve. Let’s assume the CHP plant provides N types of ancillary services {S1, S2,· · ·, Si,· · ·, SN} that have required times of their deployment {TS1 , TS2 ,· · ·, TSN } and their contracted volumes are {PS1 , PS2 ,· · ·, PSN }. The deployment time is measured between the activation command sent for the given AS and real increase/decrease of the generation to the contracted volume of a given AS. For instance, S1 can be secondary upward regulation with the deployment time TS1 = 10 min c = 5 MW, S2 secondary downward regand contracted volume PS1 c = −12 MW, S3 tertiary upward regulation ulation with volume PS2 etc. In order to model technological feasibility of AS deployment, constraints (1)–(7) will be formulated for different operating points of the plant: c is sup• a base operating point when the contracted electricity Pop plied to the grid (no AS activated), • operating points when a service Si is fully deployed, with the c + Pc , corresponding required electricity supply PSi = Pop Si • points of minimal Pmin and maximal Pmax possible electricity supply under a given commitment of boilers and turbines and given heat demand.

Hence, in the following text, the subscripts op, Si, min and max denote variables belonging to the above listed operating points. For purposes of energy and AS contracts the plant is evaluated as a virtual block so, in addition to modeling mass and energy balance in each of the operating points, the supply of energy to the grid must satisfy following contractual constraints



TG aux c Pop − Pop = Pop

(14a)

∀TG



TG aux c c PSi − PSi = Pop + PSi = PSi ,

∀TG



TG aux c Pmax − Pmax ≥ Pop +

∀TG

 ∀TG



c PSi ,

∀i

(14b)

∀i : PSic > 0

(14c)

∀i : PSic < 0

(14d)

i

TG aux c Pmin − Pmin ≤ Pop +



c PSi ,

i

The constraint (14a) models the base operating point and (14b) provision of individual AS. The constraint (14c) ensures that there is enough generation capacity for the case when all positive AS (upward regulation) are activated simultaneously. Similarly, the constraint (14d) handles the negative case. From the balance point of view, the maximal generation is typically limited by capacity of boilers or the cooling unit in case the

2 In the AS market considered in this paper, AS suppliers receive payments for capacity made available to the TSO and for energy actually provided (or received).

OUT OUT Pboiler,Si − Pboiler,op ≤ TSi · rboiler ,

∀i : PSic > 0

(15a)

OUT OUT down Pboiler,op − Pboiler,Si ≤ TSi · rboiler ,

∀i : PSic < 0

(15b)

up

OUT OUT Pboiler,max − Pboiler,op ≤ max(TSi ) · rboiler ,

∀i : PSic > 0

(15c)

OUT OUT down − Pboiler,max ≤ max(TSi ) · rboiler , Pboiler,op

∀i : PSic < 0

(15d)

up

down are the maximal ramp-up rate The parameters rboiler and rboiler and ramp-down rate of the boiler, respectively, e.g. in MW/min. Contrary to the ramping constraints (8) and (9), which represent the energy supply ramping constraints between the time steps c (t) and P c (t + 1)), the constraints (15) treat power (between Pop op output changes within one time step.

4.2. Associated costs Since the revenues from the AS are predetermined prior the dayahead optimization, only costs associated with AS provision are modeled. Typically, but not necessarily, activation of upward regulation causes increased fuel costs and downward regulation leads to fuel savings compared to the base operating point. The fuel costs Cfuel (t) in each time step are approximated as

⎡

zop



IN Pboiler,op (t)+

⎢ boilers ⎢    ⎣+ z P IN

Cfuel (t) = cfuel · ⎢

Si

zop +



i

boiler,Si

⎤

⎥ ⎥ ⎥ ⎦ (t)

(16)

boilers

zSi = 1

i

where cfuel is a fuel price per MWh of its combustion heat (including price of CO2 allowances) and zSi are user-defined weights from range 0,1 corresponding to how often the respective AS is expected to be activated (probability of activation). Depending on the control area characteristic, dispatch rules or customs in the control area or AS typical activation frequency the weights can be set by experience to make a best possible estimate of costs (16). If the costs associated with the AS activation can be neglected (for instance, if the revenues from AS are high enough) the weights zSi are set to zero to make the optimization problem simpler and (16) then represents only fuel costs for generation at the base operating point (zop = 1). 5. Optimization problem formulation As it was described before, purpose of the proposed optimization task is to optimize a CHP plant operation for a 24-h horizon on a day-ahead basis. The main parameters and inputs to the optimization problem are • technical parameters of the plant (production characteristics of turbines, boilers, enthalpies in the steam cycle, ramping rates etc.), • planned outages of individual components of the plant, • predicted heat demand Qdemand (t),

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

53

• volumes of energy P c (t) and AS volumes P c (t) already conlong Si tracted on a longer horizon (prior the day-ahead optimization), • expected day-ahead energy market opportunities (expected price of energy cshort (t)), • fuel price. The objective of the optimization is to maximize the profit in terms of revenues vs. fuel costs. The revenues consist of revenues from heat supply and traded energy and AS. The revenues from heat are given by the heat demand and usually a regulated tariff is applied, so they have no impact on optimality. Also, revenues from already traded energy and AS are fixed. So the optimizer can increase the profit by revenues from recommended optimal volume c Pshort to be traded on the day-ahead market and by minimizing the fuel costs by optimal allocation of the contracted energy and AS to the plant’s technology in each time step. This can be formalized as

max

T  



c cshort (t)Pshort (t) − Cfuel (t) − Csusd (t)

(17)

t=0

Fig. 5. Cost characteristic of the studied CHP plant.

subject to c c = Plong (t) + Pshort (t),

(18)

constraints (1)–(7) and (14)–(16) applied to each time step t in the optimization horizon T and constraints (8)–(13) applied to the optimization horizon as a whole. In the objective (17) the term c (t) represents revenues from additional sales of electricshort (t)Pshort cal energy on the day-ahead market, Cfuel (t) are fuel costs associated with the plant operation and Csusd (t) are total start-up and shutdown costs if any of the units are recommended to be started or shut down during the optimization horizon. The main outputs of the optimization task are • recommended volume P c (t) to be traded on the day-ahead short market, • commitment of boilers, turbines and other pieces of the equipment in each time step, • expected fuel costs, • optimal technological realization of the base operating point c (t) and recommended dispatch of regulation reserves P c (t) Pop Si (including allocation of energy and AS to individual turbines). 6. Model implementation and results The optimization problem for the real plant was programmed in Matlab with the help of Yalmip toolbox [20]. CPLEX 11.1 was used as an external MILP solver. The size of the problem for a 24 h horizon is up to 70 000 constraints and 20 000 variables with 2000 variables being binary (AS provision makes the scheduling problem larger since the constraints (1)–(7) have to be replicated for each AS type). The model parameters were determined from turbine and boilers certificates, enthalpy measurements and/or tuned on the historical operation data and adjusted according to the operators’ experience. The technological model (the set of balancing equations and parameters) was verified against historical measurements and was tested by the plant operators before the optimization module was built upon it. In practice, a typical computational time to get an optimal solution for a one-day horizon is in order of minutes, in most cases less than 5 min for a large CHP plant (for an optimality gap equal or less than 1%). Since the real optimization results and values cannot be presented due to the confidentiality issues the capability of the optimization tool will be illustrated on an example of a similar CHP plant consisting of two steam boilers and one turbine where the

boilers work to a common steam header. A condenser, a reduction station and a steam reheat cycle are also included. Fig. 5 shows a cost characteristic of the studied plant that was computed using the described modeling approach. Two scenarios for the day ahead planning will be considered: one without ancillary services and the second with ancillary services. What both scenarios have in common is the heat demand in Fig. 6 and c the baseload energy diagram Plong (t) = 20 MW already traded on longer-horizon markets. The optimizer is asked to recommend the c optimal energy hourly diagram Pshort (t) to be additionally traded on the day-ahead market based on expected prices that were set to 75 D /MWh in the peak hours (8:00–20:00) and 50 D /MWh in the off-peak hours. For the first scenario, the recommended optimal generation c aux (t) of the plant is depicted in TG (t) = P c level Pop (t) + Pshort (t) + Pop long Fig. 7a. Such a generation is realized by two boilers, a gas-fired one (Fig. 7b) and a coal-fired one (Fig. 7c). The much more expensive gas-fired boiler is started only for the morning peak of the heat demand and is operated at its minimum. Since the prices of electricity are relatively low the turbine is optimally operated in the backpressure mode (the electricity trading and generation is correlated with the heat demand) except for a few afternoon hours where the condensing unit is partially utilized (Fig. 7d). The second scenario differs from the first one in that sense, that also AS have been traded on longer-horizons. In our example, the previously contracted volumes are: positive tertiary reserve c = 17 MW in the peak period, P c = 9 MW off-peak, and negaPS3 S3 c = −24 MW peak period and P c = −21 MW tive tertiary reserve PS4 S4 off-peak. These volumes appear in the constraint (14) and the optic (t) to be traded on mizer is asked again to recommend energy Pshort the day-ahead market while satisfying this constraint. The result is

80

Heat demand [%]

c (t) Pop

Heat demand 60 40 20 0

0

4

8

12 hours

16

Fig. 6. Heat demand for the two studied scenarios.

20

24

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

Turbogenerator output [MW]

(a)

60 40 20 0

0

4

8

12

16

20

24

Boilers output

(b) 100 80 60 40 20 0

0

4

8

12

16

20

24

100 80 60 40

0

base point

0

4

8

12

16

20

24

Condeser output [%]

100 80 60 40 20 0

0

4

8

12

16

20

24

40 20 0

0

4

8

12

16

20

24

12

16

20

24

12

16

20

24

Boilers output 100 80 60 40 20 0

0

4

8

100 80 60 40

base point AS- activation AS+ activation

20 0

0

4

8

Condenser and reduction station output 100 80 60 40 20 0 0

4

8

0

4

8

12

16

20

24

12

16

20

24

(d) Reduction station [%]

(e) Reduction station [%]

60

(c)

Condenser and reduction station output

(d) Condeser output [%]

Boiler #2 output [%]

Boiler #2 output [%]

(c)

20

Turbogenerator electric output

(b) Boiler #1 output [%]

Boiler #1 output [%]

(a)

Turbogenerator electric output

Turbogenerator output [MW]

54

100 80 60 40 20 0

0

4

8

12

16

20

24

100 80 60 40 20 0

hours

hours

Fig. 7. Production plan for scenario #1.

Fig. 8. Production plan for scenario #2.

depicted in Fig. 8a where besides the optimal base operating point TG (t) of the turbine also the generation level P (t) for the positive Pop S3 reserve deployment (blue line) and PS4 (t) for the negative reserve deployment (green line) is shown. Because of the AS constraints the resulting optimal base generation diagram is lower than the one in scenario #1 to enable larger regulation range and is not correlated with the heat demand anymore. Also, operation of the plant’s technology is different: the gas-fired boiler #1 is turned on for a longer period to leave enough steam generation flexibility for positive AS provision (Fig. 8b). The condensing regime of the turbine (Fig. 8d) is used to satisfy activation of the positive AS and, quite surprisingly, also for the negative AS in the afternoon hours. This is

caused by the fact that the steam boilers run at their minimums and cannot go lower so the generated heat surplus has to be cooled in the condenser. The reduction station is used in hours of higher heat demand and for activation of negative AS to bypass the steam from boilers directly to the heating system exchangers so that the turbine electricity output remains at the desired (low) level (Fig. 8e). 7. Conclusions An industrial application of mathematical modeling, identification and advanced mathematical optimization of a plant was presented in this paper with emphasizes on AS provision. The

P. Havel, T. Sˇ imoviˇc / Electric Power Systems Research 95 (2013) 47–55

devised AS optimization model is quite complex and increases the optimization problem size several times; however, a simpler model would not able to describe the technical capabilities of the plant with enough fidelity so that all possible AS market opportunities can be exploited and reliable provision of AS can be guaranteed. Without such a model the sale and provision of AS would have to be more conservative in order to avoid high penalties for failing to provide the entire range of contracted AS later in real time when the AS are activated by the transmission system operator. The resulting decision-support software tool helps the CHP plant’s operator to run the technology in the most efficient way while the market opportunities are fully exploited. The tool has been in operation at two larger CHP plants consisting of several boilers and turbines for more than 2 years. It has turned out that the linear or piecewise linear model of the plant and its components has an acceptable error compared to other uncertainties in operation and in input data. Also, there have been no complaints on the model fidelity from the plant’s operators so far. Compared to the previously used approach of the plant production planning (MS Excel & operator’s experience) the presented optimization approach based on first-principle mathematical model automatically considers all possible operation scenarios to find the optimal one, even though such a scenario has never been applied before in the past operation of the plant. Besides other benefits, it has allowed the plant operator to run the plant with fewer boilers during the spring/autumn period of changing heat demand, yielding substantial economic and environmental benefits. Hence it helps the CHP plant to remain competitive in nowadays dynamically changing energy markets. Acknowledgments This work was supported by the Ministry of Education of the Czech Republic under the project No. 1M0567 and the company Taures, plc. References [1] UCTE Operation Handbook, Policy 1: Load-Frequency Control and Performance (Online) https://www.entsoe.eu/resources/publications/systemoperations/operation-handbook/ (accessed 10.06.12).

55

[2] M. Variny, O. Mierka, Improvement of part load efficiency of a combined cycle power plant provisioning ancillary services, Applied Energy 86 (June (6)) (2009) 888–894. [3] G. Sheble, G. Fahd, Unit commitment literature synopsis, IEEE Transactions on Power System 9 (1994) 128–135. [4] F. Salgado, P. Pedrero, Short-term operation planning on cogeneration systems: a survey, Electrical Power System Research 78 (2008) 835–848. [5] A. Rong, R. Lahdelma, An efficient envelope-based branch and bound algorithm for non-convex combined heat and power production planning, European Journal of Operational Research 183 (2007) 412–431. [6] G. Tao, M.I. Henwood, M. van Ooijen, An algorithm for combined heat and power economic dispatch, IEEE Transactions on Power System 11 (1996) 1778– 1784. [7] E. Thorin, H. Brand, C. Weber, Long-term optimization of cogeneration systems in a competitive market environment, Applied Energy 81 (2005) 152– 169. [8] M.A. Gonzalez Chapa, J.R. Vega Galaz, An economic dispatch algorithm for cogeneration systems, in: IEEE 2004 Power Engineering Society General Meeting, 6–10 June 2004, 2004, pp. 989–994. [9] S. Chen, M. Tsay, H. Gow, Scheduling of cogeneration plants considering electricity wheeling using enhanced immune algorithm, International Journal of Electrical Power and Energy System 27 (2005) 31–38. [10] M. Tsay, Applying multi-objective approach for operation strategy of cogeneration systems under environmental constraints, Electrical Power and Energy System 25 (2003) 219–226. [11] E. Dotzauer, Algorithms for short-term production planning of cogeneration plants, Ph.D. Thesis, Linkoping University, Linkoping, Sweden, 1997. [12] S. Ashok, R. Banerjee, Optimal operation of industrial cogeneration for load management, IEEE Transactions on Power Systems 18 (2) (2003). [13] E. Gomez-Villalva, A. Ramos, Optimal energy management of an industrial consumer in liberalized markets, IEEE Transactions on Power System 18 (2003) 716–723. [14] S.W. Ilerhaus, J.F. Verstege, Optimal operation of industrial CHP-based power systems in liberalized energy markets, in: IEEE Power Tech. Conference, vol. BPT99-352-12, August 29–September 2, 1999. [15] J.M. Yusta, P.M. De Oliveira-De Jesus, H.M. Khodr, Optimal energy exchange of an industrial cogeneration in a day ahead market, Electric Power System Research 78 (2008) 1764–1772. [16] G. Ferrari-Trecate, et al., Modeling and control of co-generation power plants: a hybrid system approach, IEEE Transactions on Control, System and Technology 12 (2004) 694–705. [17] M. Carrion, J. Arroyo, A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem, IEEE Transactions on Power Systems 21 (2006) 1371–1378. [18] K. Hedman, R. O’Neill, S. Oren, Analyzing valid inequalities of the generation unit commitment problem, in: Power Systems Conference and Exposition, PSCE’09, 2009, pp. 1–6. [19] R. Raineri, S. Ríos, D. Schiele, Technical and economic aspects of ancillary services markets in the electric power industry: an international comparison, Energy Policy 34 (2006) 1540–1555. [20] J. Löfberg, YALMIP: a toolbox for modeling and optimization in MATLAB, in: Proceedings of the CACSD Conference, Taipei, Taiwa, 2004.