Development of a dynamic model for start-up optimization of coal-red power plants

(1)

Development of a dynamic model for

start-up optimization of coal-red power

plants

EVA ANDERSSON

Master's Degree Project

Stockholm, Sweden June 2013

(2)

Abstract

The expansion of renewable energies and the deregulations of the energy market are in-creasing the demand of regulating power. In absence of hydropower, thermal power is often used for this purpose instead. This thesis focuses on the Vattenfall owned power plant Jän-schwalde located in Germany. The goal is to optimize the start-up procedure, so that the start-up time is reduced without causing too much thermal stress on the important thick-walled components in the boiler. By reducing the start-up time, the plant can become more protable, exible and better suited to regulate the electricity market.

A start-up model was built in Dymola and validated against measurement data and a sim-ulation model of Jänschwalde, which is too complex to use for optimization purposes. The JModelica.org platform was used for the optimization part of the project. It was possible to nd optimal solutions for the start-up process of the Jänschwalde power plant, but the convergence of the optimization algorithm was very dependent on the optimization options used and the scaling of the plant model. Further work includes development of the compo-nents used in the start-up model, rening the discretization and scaling for the optimization problem.

(3)

List of Figures

1 Illustration of a coal-red power plant. Coal is burnt in the furnace to heat up water in the heat exchangers. The hot steam is used to drive a steam turbine, which generates electricity to the power grid [Ignou,2013]. . . 7 2 The Carnot Cycle presented in an temperature-entropy (T-S) diagram. The heat

QIN is added to the system and QLOSS is lost to a cooling medium. The heat

engine operates between TH and TL (with TH > TL) and the total eciency is

given by (1). . . 8 3 The Rankine Cycle presented in an temperature-entropy (T-S) diagram. Heat

QIN is added to the system in several steps of the steam cycle and QLOSS is lost

when the steam is condensed. The steam is reheated to improve the eciency of the plant. . . 9 4 Measurement data of a hot start-up. The two upper subplots shows the

tempera-ture and pressure for live steam and reheated steam. In the middle, the generator and ring power is presented next to the live steam mass ow rate. The two lower subplots shows revolutions per minute (RPM) for the turbine and the opening level of the high-pressure and low-pressure bypass valves. . . 11 5 Model of a heat exchanger, where heat is transferred from ue gas to steam

through a wall segment. . . 12 6 A cross-sectional illustration of a header. . . 14 7 Life-time in terms of load cycles N of components related to stress amplitude ∆σ,

at T = 100◦_C_{and T = 550}◦_C_{. . . 15}

8 An element ne in the interval [ti, ti+1],containing three Radau points q1, q2 and q3. 17

9 An overview of the complex model of Jänschwalde, owned by Vattenfall. The complete steam cycle is included in the model. The water is heated, used to drive the turbine, reheated and nally condensed and cooled. The total number of equations and unknowns is 16,927. . . 23 10 The superheater section is one subcomponent from the model shown in Figure

9. The gure shows four superheaters (SH) and three reheaters (RH) from the complex model of the Jänschwalde power plant. The circles represent connectors for ue gas (green), water (blue) and coal (brown), and are used to link the component to other submodels or boundary conditions. . . 24 11 The furnace section in Figure 9, which contains the combustion process and the

evaporator. . . 25 12 The ECO & SH1 section from Figure 9, which shows the economizer and one

superheater. . . 25 13 The gure shows how well the process models P1, P1D, P0 and P0D, generated

in the System Identication Toolbox in Matlab, ts the reference data. . . 28 14 Model used for start-up optimization. The heat exchangers modelled are the

superheaters (SH), reheaters (RH) and the evaporator (Evap). Other components include wall and volume components to calculate stress levels in headers and separators and components to represent boundary conditions. The model has one input du connected to the integrator component which describes how the load is changed du

dt. There are three outputs to the right for live steam temperature,

header stress and separator stress. . . 31 15 Temperatures of in- and outowing water, for the complex superheater and the

(5)

16 Temperatures of in- and outowing ue gas, for the complex superheater and the start-up superheater. The temperature curves for inowing ue gas are equal for

both superheaters. . . 34

17 Comparison between reference data, complex simulation model and start-up model at steady-state with full load. . . 35

18 Comparison of the complex simulation model of Jänschwalde (dotted lines) and measurement data (solid lines) from the actual plant, when the system was ramped from 100 % load to 60 % and back again to 100 % load. . . 36

19 Simulation of the start-up model when the system was ramped from 100 % load to 60 % and back again to 100 % load. The start-up model is not aected by disturbances, compared to the systems shown in Figure 18. . . 37

20 Plot of coal ow u and the control signal du dt, with blocking factors [1, 1, 2, 6]. . . 39

21 Initial optimization with blocking factors [1,1,2,6]. The weight parameters in (20) are chosen as α = 0.01, β = 0.01. The steady-state value Tref is reached for the live steam exiting the last superheater SH4. . . 41

22 Initial optimization with blocking factors [1,1,2,6]. The weight parameters in (20) are chosen as α = 0.0004, β = 0.01. The steady-state value for Tref is not reached during 2000 seconds. . . 41

23 Initial optimization with blocking factors [1,1,2,6]. The weight parameters in (20) are chosen as α = 0.0025, β = 0.01. The nal value of the steam temperature in SH4 is 1040 K. . . 42

24 Optimization with uniformly distributed blocking factors [2, 2, 2, 2, 2, 2, 2, 2] and ne= 16. The stress amplitude is approximately 135 MPa. . . 43

25 Optimization with a high blocking factor at the end of the interval [2, 2, 2, 2, 2, 2, 6] and ne= 16. The stress amplitude is approximately 210 MPa. . . 44

26 Optimization with blocking factors [1, 1, 1, 1, 1, 5] and ne= 10. The stress ampli-tude is approximately 200 MPa. . . 44

27 Optimization with blocking factors [5, 1, 1, 1, 1, 1] and ne= 10. The stress ampli-tude is approximately 125 MPa. . . 45

List of Tables

1 Start-ups at Jänschwalde during 2012 . . . 10

2 Optimization options for JMUs . . . 21

3 Optimization options for the IPOPT solver . . . 21

4 Process models used in the System Identication Toolbox . . . 28

5 Steam temperatures out of each heat exchanger component at a steady-state simulation with full load. . . 35

6 Specication of values used for Tref and constraints. . . 39

7 Summary of parameter values used for dierent problem set-ups. . . 40

(6)

1 Introduction

Recent development in the electricity market has started and will continue to change how power plants are operated. The increase in renewable energies is increasing the demand of regulating power needed, to balance production and consumption of electricity, so that the frequency of the electricity grid is kept constant. In addition to this, the deregulations of the electricity markets makes it possible for the consumers to also produce electricity. The Swedish market has access to hydropower, which is a very good regulating power since it can be switched on and o almost instantly. In other countries it is common the use thermal power plants instead. Besides having a large environmental impact, thermal power plants require a longer start-up time.

This thesis will focus on the start-up procedure of the German power plant Jänschwalde, owned by Vattenfall. The goal is to nd a model that describes the behaviour of the power plant suciently well and can be used for optimization. By formulating the start-up procedure as an optimization problem, one can also take into account economical and environmental factors, for example by reducing the fuel consumption during the start-up. During start-ups there are certain components, that can be exposed to heavy thermal stress. Too much thermal stress, will shorten the components life-time and therefore constraints on stress levels are needed.

Several studies on start-up optimization have been done in the past, but most of them consider Combined Cycle Power Plants, for example [Lind & Sällberg, 2012], [Casella et al.,2011] and [Casella, Donida & Åkesson, 2011]. The Combined Cycle Power Plants consists of both a gas turbine and a steam turbine, and are run with gas. The model derived in this thesis will be dierent from previous work, mainly because Jänschwalde is run with brown coal, which intro-duces delays in the overall system. Another dierence is that Jänschwalde only has turbines that are driven by steam. A very detailed model of Jänschwalde is under development at Vattenfall. However, the high complexity of it makes it unsuitable for optimization purposes. There exists some studies on coal-red boilers, including for example [Franke, Rode & Krüger, 2003]. The models and optimization methods from those studies are unfortunately not public.

1.1 Aim of Thesis

The aim of the thesis is mainly to develop a dynamic model of Jänschwalde, that can be used for optimization purposes. Another goal is to try out the optimization techniques within the JModelica.org platform [JModelica.org,2013], to see if they can be applied to a dynamic model of a coal-red power plant. The nal goal is to use dierent optimization options and try to nd optimal solutions for the start-up problem.

Since the actual power plant is very complex, there are many simplications that must be done in order to model it. The model developed within this thesis only contains the most important components relevant for the start-up, which are located in the boiler.

1.2 Structure of thesis

In Chapter 2 a brief introduction to coal-red power plants is given. The mathematical equations that are used for the modelling are explained and derived. It is also shown how the dynamic optimization problem can be transferred to a non-linear programming problem, which can be solved by the chosen optimization algorithm.

(7)

The relevant programs and how they are used are briey explained in the Chapter 3. The chap-ter also includes examples on how the JModelica.org platform is used for optimization.

In the beginning of Chapter 4 the complex simulation model of Jänschwalde, owned by Vat-tenfall, is described. The components used in the start-up model are explained and the most important code is included. The start-up model is compared to the complex model and mea-surement data from the real plant in Chapter 5.

Chapter 6 includes the optimization set-ups used and the outcomes are presented in Chapter 7. Finally, the model and the optimization outcomes are evaluated in Chapter 8. Further improvements and sources of error are also discussed.

(8)

2 Background

This chapter is divided into two parts. First an introduction to coal-red power plant is given, where the steam cycle and the start-up procedure are presented. The second part describes the mathematical theories, that are needed to understand the dynamic model and the optimization techniques used.

2.1 Power Plant

The subject of study is the power plant Jänschwalde, which is owned by Vattenfall[Vattenfall,2012]. It is the largest lignite power plant in Germany, consisting of 6 blocks with 500 MW installed capacity each [Saarinen, Boman & Funkquist, 2011]. The power plant started operating in 1981 and has an electrical eciency of 36 %.

2.1.1 Overview

A basic illustration of a simple coal-red power plant is shown in Figure 1. The coal is burnt in the furnace to produce ue gas and ash. The heat from the ue gas is used to warm the feed-water and convert it to steam in the boiler. The feedfeed-water goes through dierent types of heat exchangers, including the economizer, the evaporator and several superheaters. The superheated steam is then expanded in the steam turbine. The turbine is connected to a generator, which is used to generate electricity to the power grid. Finally the steam is condensed and cooled down, before it re-enters the boiler [Ignou,2013].

Since the power plant at Jänschwalde is much more complicated than this example, a more detailed description of the steam cycle will be given in Section 2.1.2. In this thesis the main focus will be on the boiler, because it contains the components that are exposed to the highest levels of stress.

Figure 1: Illustration of a coal-red power plant. Coal is burnt in the furnace to heat up water in the heat exchangers. The hot steam is used to drive a steam turbine, which generates electricity to the power grid [Ignou,2013].

(9)

2.1.2 The Steam Cycle

The eciency of the power plant can be evaluated by looking at the Rankine cycle, which is a modication of the Carnot cycle. The Carnot cycle is an ideal process, where all processes are reversible and no heat is lost to the surrounding environment [Alvarez, 2006]. It is a theoretical model, which maximizes the eciency of a heat engine. The Carnot cycle consists of four thermodynamic steps, which are shown in a temperature-entropy diagram in Figure 2. In this example the heat engine is working in the two-phase region, meaning that the water exists in both liquid and vaporized states at the same time.

Figure 2: The Carnot Cycle presented in an temperature-entropy (T-S) diagram. The heat QIN is added

to the system and QLOSS is lost to a cooling medium. The heat engine operates between TH and TL

(with TH> TL) and the total eciency is given by (1).

Carnot cycle

1 → 2: First the medium is expanded at a constant temperature TH, by absorbing heat QIN.

2 → 3: The next step consists of an adiabatic expansion of the gas, in which no heat is lost to the surrounding environment. The temperature decreases to TL and the entropy is constant

since no heat is added to or subtracted from the system.

3 → 4: The medium is compressed at the constant temperature TL and the heat QLOSS is

transferred to a cooling medium.

4 → 1: Finally the medium is adiabatically compressed, until the temperature TH is achieved

and the cycle is complete.

The Carnot eciency η of the heat engine is given by the work done by the system W , divided by the heat added to the system QIN:

η = W QIN = QIN− QLOSS QIN = 1 −QLOSS QIN = 1 − TL TH . (1)

Equation (1) shows that it is desirable to have a large dierence between TL and TH, to get a

high Carnot eciency.

A coal-red power plant is better described by the Rankine cycle with reheat [Ignou,2013]. A Rankine cycle that operates in the one- and two-phase region, is shown in Figure 3.

(10)

Figure 3: The Rankine Cycle presented in an temperature-entropy (T-S) diagram. Heat QIN is added

to the system in several steps of the steam cycle and QLOSS is lost when the steam is condensed. The

steam is reheated to improve the eciency of the plant.

Rankine cycle with reheat

1 → 2First the water is compressed by a pump and the temperature is increased.

2 → 3The temperature of the water is further increased in the pre-heaters and the economizer, by addition of heat QIN.

3 → 4 Next, the water is evaporated into steam. The evaporator operates in the two phase region, containing both water and steam.

4 → 5The steam temperature is further increased in the superheaters.

5 → 6In the next step, steam is expanded in the high pressure turbine to generate electricity. 6 → 7After the high pressure turbine, the steam is reheated and led to the low pressure turbine. 7 → 8The steam is expanded again and electricity is generated in the low pressure turbine. 8 → 1Finally the steam is condensed completely.

The eciency of the actual process can never be higher than the Carnot eciency [Alvarez, 2006]. The best ways to achieve a high eciency is to have small energy losses throughout the cycle, and a very high dierence between the maximal and minimal water temperature in the cycle. A dierence between Figure 2 and 3, is that the water is superheated to only exist in a vaporized state, so that a very high maximal temperature can be reached. The life-time of the turbine can be improved by using superheated steam which is very dry, to generate electricity. Another dierence is the reheating step, which is another way to increase the eciency of the cycle. The steam is only partially expanded in the high pressure turbine, reheated and then completely expanded in the low-pressure turbine.

2.1.3 Start-up Procedure

The start-up procedure of Jänschwalde can be divided into three dierent phases.

(i) First the ignition oil burners are started up, to quickly start the combustion process in the furnace. Then the coal mills are started up and the combustion medium is gradually

(11)

Table 1: Start-ups at Jänschwalde during 2012

Start-up type Stand-still time Occurrences 2012

Hot t ≤ 12h 4

Warm 12h < t < 48 h 6

Cold 48h ≤ t 15

shifted from oil to coal during the start-up. Meanwhile the temperature and pressure of the live steam are ramped up. During this stage the power plant is not producing any electricity, since the live steam bypasses the steam turbine.

(ii) The next phase is to accelerate the steam turbine, which is started once you have reached stable steam conditions. The turbine is accelerated until it reaches 3000 revolutions per minute.

(iii) In the last phase, the bypass valve is gradually closed so that the live steam enters the steam turbine. The production of electricity begins and the load is increased until the nal output capacity is reached.

The length of the start-up is highly dependent of the stand-still time of the power plant, since it takes a long time for the metal parts to cool down. Start-ups are therefore categorized into hot,warm and cold start-ups according to Table 1. In Figure 4 measurement data are presented for a hot start. The ignition oil burners are initiated in the interval t ∈ [15, 30] min and the coal mills are started when t ∈ [40, 70] min, which is represented by the change in ring power. Stable steam conditions for the live steam mass ow and pressure are reached for t ≈ 100 min, but one can observe that the turbine is started at t ≈ 75 min which indicates an overlap between phase 1 and 2. In the lower right corner, the opening of the high- and low-pressure valves are shown. When the valves are open, the steam is bypassed the turbine. The valves are closed at t ≈150 min, meaning that the steam goes through the turbine, which generates electricity. The nal phase of the start-up includes increasing the load (ring power) and live steam ow rate, which is done in the interval t ∈ [150, 210] min.

2.2 Theory

A thermal power plant consists of several heat exchangers, which warms up the water gradually during the steam cycle. Therefore the theory behind a heat exchanger component will be given in this section. Since the start-up process is limited by stress in thick-walled components, the stress calculation will be presented as well. Finally the dynamic optimization problem is set up and transferred to a non-linear programming problem.

2.2.1 Balance Equations

In Figure 5 a model of a heat exchanger is shown, which can be used in thermal power plants. The component consists of two sides, containing ue gas and steam which are divided by a wall. The ue gas is obtained from the combustion process and is assumed to have a higher temperature than the wall.

(12)

Figure 4: Measurement data of a hot start-up. The two upper subplots shows the temperature and pressure for live steam and reheated steam. In the middle, the generator and ring power is presented next to the live steam mass ow rate. The two lower subplots shows revolutions per minute (RPM) for the turbine and the opening level of the high-pressure and low-pressure bypass valves.

(13)

Hot ue gas with enthalpy hgas,in, mass ow rate qgas,in and pressure pgas,in enters the gas side

of the heat exchanger. Since the wall temperature is lower than the gas temperature, heat Qgas

will be transferred to the wall according to the second law of thermodynamics. The ue gas then exits the heat exchanger with enthalpy hgas,out, mass ow rate qgas,outand pressure pgas,out.

The process of the steam side works in the corresponding way, except that heat Qwater is added

to the system.

Figure 5: Model of a heat exchanger, where heat is transferred from ue gas to steam through a wall segment.

The dynamics of each side of the heat exchanger are given by the balance equations, dMm

dt = qm,in− qm,out, (2a)

dUm

dt = qm,inhm,in− qm,outhm,out+ Qf low,m, (2b) pm,in− pm,out= fm

q2_m,in ρm

, (2c)

with m ∈ {gas, steam}, mass of the matter Mm, internal energy Um, density of the medium

ρm and a constant fm [Krüger, Franke & Rode, 2004]. fm describes the pressure drop and can

be found experimentally. For the balance equations heat ow rate Qf low,m = dQm_dt , is used to

describe at which rate the heat is transferred.

For this application it is more convenient to use pressure p and enthalpy h as state variables instead of mass M and internal energy U. The pressure p and temperature T for dierent stages of the process are measured in the real plant. It is however easier to use enthalpy h instead of temperature T , because the right-hand side of the energy balance equation is expressed in terms of enthalpy.

The mass of the matter Mm can be expressed in terms of volume Vm and density of the matter

(14)

way: dMm dt = Vm dρm dt (3a) = Vm ∂ρm ∂hm p dhm dt + ∂ρm ∂pm h dpm dt = qm,in− qm,out, (3b) dUm dt = d(Vmρmhm− pmVm) dt = Vm d(ρmhm) dt − dpm dt (3c) = Vm hm dρm dt + ρm dhm dt − dpm dt (3d) = Vm hm ∂ρm ∂hm p dhm dt + hm ∂ρm ∂pm h dpm dt + ρm dhm dt − dpm dt (3e) = Vm hm ∂ρm ∂hm p + ρm dhm dt + hm ∂ρm ∂pm h − 1 dpm dt (3f) = qm,inhm,in− qm,outhm,out+ Qf low,m. (3g)

In section 4.2.4 the mass balance Vm ∂ρm ∂hm p dhm dt + ∂ρm ∂pm h dpm dt = qm,in− qm,out, (4)

and the energy balance Vm hm ∂ρm ∂hm p + ρm dhm dt + hm ∂ρm ∂pm h − 1 dpm dt

= qm,inhm,in− qm,outhm,out+ Qf low,m, (5)

is used to model the dynamics of gas and water in the heat exchangers.

To calculate how much heat is transferred between the gas and water side, equations for the wall segment must be included. The change in temperature of the wall, depends on the heat capacity cp,wall, the mass Mwall and how much heat is added or subtracted.

Mwallcp,wall d(Twall) dt = X m∈gas,steam Qf low,m, (6) where Qm is given by

Qf low,m= αwall,mAwall,m(Tm− Twall), (7)

with the heat transfer coecient αwall,m and the area of exposure Awall,m between the wall and

the medium m [Alvarez, 2006]. 2.2.2 Stress Calculation

During the start-up, the thick-walled components are aected by thermal stress σT and

mechan-ical stress σm, due to changes in temperature and pressure [Meinke, 2012]. The thick-walled

components that are considered most critical during start-up are the separator and the headers. A header consists of a large hollow pipe and many smaller pipes, which can be seen in Figure 6. The steam enters or exits the header through the smaller pipes, which introduces extra stress at the exit points. In the stress calculations the worst case scenarios are considered. The total

(15)

stress σ is given by the mechanical stress σm and the stress induced by changes in temperature

for the material σT

σ = σm+ σT. (8)

The thermal stress σT is given by

σT = αT

βLTET

1 − ν (Tm− Ti). (9)

The elasticity module ET, the thermal expansion coecient βLT and the Poisson's ratio ν are

given by the material properties. Ti corresponds to the inner temperature of the wall, Tm to the

mean temperature and αT is a constant.

The mechanical stress σm is given by:

σm= (αm0+ fuαb)p

dm

sb

, (10)

where dm corresponds to the inner diameter of the header, sb is the thickness of the wall and

p is the pressure on the wall. αm0 and αb are constants. To take into account the extra stress

introduced on the exit point, a form factor fu is introduced:

fu = 1.5 dm sb 1 +1−ν₂ 2_ETp (dm_sb)3 U 100. (11)

The constants αT, U, αb, αm0 are dependent on the material and the geometry. In the specic

header problem the parameter values αT = 2.0, U = 10−5, αb = 2.9 and αm0 = 3.2 have

been used. The same equations and parameters have been used to calculate stress levels in the separator, since the separators geometry is similar to the headers geometry. The main dierence is that the separator has a larger internal radius and consists of fewer small pipes.

Figure 6: A cross-sectional illustration of a header.

The life-time of the components can be linked to stress levels by using a Wöhler digram. The life-time is dened as the number of cycles N the material can tolerate, which depends on the stress amplitude ∆σ, given by

(16)

∆σ = σmax− σmin

2 , (12)

were σmax and σmin are the maximal and minimal values of the total stress during the start-up.

In Figure 7 a Wöhler curve is shown, which has been developed to analyse life-times of compo-nents in thermal power plants [Meinke, 2012]. The total number of start-ups for 2012 was 25, which can be seen in Table 1. If one assumes that a power plant is operated for 60 years and is started up 25 times/year, the total number of start-ups becomes 1500. Figure 7 shows that the stress amplitude ∆σ must be below 700 MPa, to not wear out a header that operates around 550◦C. For the separator, a slightly higher stress amplitude is tolerated, since it operates at ≈ 350◦C.

Figure 7: Life-time in terms of load cycles N of components related to stress amplitude ∆σ, at T = 100◦_C

and T = 550◦_C_.

2.2.3 Optimization Problem

A dynamic model often consists of dierential and algebraic equations, and together they form a Dierential Algebraic (DAE) system [Ljung & Glad, 2004]. A dierential equation describes the relationship between a function and its derivative [Edsberg, 2008]. Two typi-cal dierential equations are the balance equations for mass and energy, seen in (4) and (5). The algebraic equations ensures the physical and thermodynamical relations in the model [Biegler, Cervantes & Wächter, 2002], for example (7), which describes how heat is transferred between media. A general dierential equation

(17)

can be rewritten by introducing x(t) = (z(t), y(t))T _{to separate the dierential variables z, whose}

derivatives appear in the DAE system, from the algebraic variables y. Equation (13) is then transformed to:

F ( ˙z(t), z(t), y(t), u(t), t) = 0, ∀t ∈ [t₀, tf]. (14)

By dierentiating (14) with respect to time, an Ordinary Dierential Equation (ODE) system can be obtained. The number of times n the system must be dierentiated, before an ODE sys-tem is achieved corresponds to the dierential index of the DAE [Ljung & Glad, 2004]. Therefore an ODE system corresponds to a DAE system of index 0.

A general formulation of a DAE optimization problem is given by [Wächter, 2002]: min

z(t),y(t),u(t),tf,p ϕ(z(tf), y(tf), tf, p) (15a)

s.t dz(t)

dt = F (z(t), y(t), u(t), t, p) (15b)

0 = G(z(t), y(t), u(t), t, p) (15c)

z(0) = z0 (15d)

0 = Hs(z(ts), y(ts), u(ts), ts, p) for s ∈ [1, ..., ns] (15e)

The variable bounds are given by:

zL ≤ z(t) ≤ zU

yL ≤ y(t) ≤ yU

uL ≤ u(t) ≤ uU

pL ≤ p ≤ pU tL_f ≤ tf ≤ tUf

ϕ is a scalar objective function,

F are the right hand sides of dierential equation constraints, G are algebraic state prole vectors of index one, Hs are additional point conditions at xed times ts,

z are dierential state prole vectors,

z0 are the initial values of z,

y are algebraic state prole vectors,

u are control prole vectors,

p is a time-independent parameter vector,

tf is the nal time.

In the general formulation the algebraic equations have been collected from (14) to form (15c). Point conditions have been dened at certain xed times and all variables used are bounded by upper and lower values. To solve the problem numerically, the dynamic optimization problem is transferred to an non-linear programming (NLP) problem. This is done by introducing a discretization so that the state and control proles are approximated by a family of polynomials on nite elements (t0 < t1 < . . . < tne < tf) [Biegler, Cervantes & Wächter, 2002]. The

dis-cretization is further rened by placing Radau collocation points ncp within each element, see

(18)

t

_i

q

1

q

2

q

3

t

_i+1

Figure 8: An element nein the interval [ti, ti+1],containing three Radau points q1, q2 and q3.

Interior Point Method

After the discretization, the NLP problem can be formulated in the following way: min

x∈Rn f (x) (16a)

s.t c(x) = 0 (16b)

xL≤ x ≤ xU. (16c)

where x = (dz

dt i,q, zi, yi,q, ui,q, t, p)T and q = 1, . . . , ncpand i = 1, . . . , ne[Wächter & Biegler, 2006].

The objective function f(x) : Rn_{→ R and equality constraints c(x) : R}n_{→ R}m _{should be twice}

continuously dierentiable and m < n. The variables x have an upper xU _{∈ R}n_{and lower bound}

xL∈ Rn_{and x}U _{≥ x}

L [Wächter & Biegler, 2006].

To demonstrate the interior point method, the following cost function will be used min

x∈Rn f (x) (17a)

s.t c(x) = 0 (17b)

x(i)≥ 0 for i ∈ I. (17c)

with I ∈ 1, ..., n is the set of indices of bounded variables [Wächter, 2002]. In (16c) the variables had upper and lower bounds, but it is easier to demonstrate the algorithm with the new prob-lem formulation only consisting of variables with an lower bound of zero [Wächter, 2002]. The solution and algorithm for the problem (17) can be applied to the problem (16), for example if the rst problem is redened with slack variables [Wächter, 2002].

(19)

The interior point method is also referred to as a barrier method, and is dened accordingly min x∈Rn ϕµ(x) := f (x) − µ n X i=1 ln(x(i)) (18) s.t c(x) = 0. (19)

A barrier function ϕµ(x)is introduced where an logarithmic barrier term is added to the problem

[Wächter, 2002]. In that way the variable bounds from (17c) are included in the objective function. If the variables approach the variable bounds in (17c), the barrier function becomes arbitrarily large and therefore the local solution x?µ is in the interior of the set, meaning that it

fulls x(i)

?µ > 0 for all i ∈ I. When the barrier function µ → 0 the solution x_?µ converges to x_?

(20)

3 Programs and Tools

The dynamic start-up model has been developed in the Modelica environment Dymola. The JModelica.org platform was used for the optimization part, with the Modelica language extension Optimica. The programs and how they are used are briey described in this chapter.

3.1 Modelica

Modelica is an open programming language, which is used to model and simulate complex physi-cal systems. The language is object-oriented, allowing you to inherit properties from submodels. Modelica is equation-based which makes it easier to combine dierent kinds of physical systems. A Modelica model can be a combination of a variety of mechanical, electrical, thermodynamical or other types of components. Another benet with using the Modelica language is that it can handle DAEs. If the modelling is instead done in Simulink [MathWorks,2013], you need to manually convert the DAEs to ODEs, to be able to simulate the model [Modelica.org,2013].

3.2 Dymola

Dymola is a commercial Modelica environment, which for example can be used for modelling and simulation. In Dymola you can switch between working in graphical and code-based envi-ronment. The changes you make in the graphical environment are automatically converted into Modelica code. In the modelling environment it is also possible to browse the open-source Mod-elica standard library and use those components to build models. The simulation environment is mainly used to simulate models and plot results. Other features include 3D-animations of results and the possibility to export models to formats that can be handled by other simulation software, such as Simulink [Modelica.org,2013].

3.3 Optimica

Optimica is an extension of the Modelica language, which makes it possible for the user to formulate dynamic optimization problems [Åkesson, 2008]. The Optimica extension includes the class optimization with class attributes such as objective, startTime and nalTime. Optimica can be used to express dierent types of optimization problems, for example minimize a cost function or minimize the nal time. It is also possible to specify constraints for the optimization problems, that must be fullled for the optimal solution.

3.4 JModelica

JModelica is an open source platform, that can be used to optimize and simulate Modelica models. The project started at the Department of Automatic Control at Lund University and now the development is done by Modelon AB and in academia [JModelica.org,2013]. JModelica supports the extension Optimica, which makes it possible to formulate and solve many types of dynamic optimization problems.

Three les are generally needed to formulate and solve a dynamic optimization problem, the Modelica model (.mo), the Optimica specications (.mop) and a Python script le (.py). The Modelica and Optimica code is compiled in the Optimica compiler to generate C code. The C code is then compiled and connected to numerical solvers, to solve the optimization problem [Åkesson, Gäfvert & Tummescheit 2009].

(21)

Simulation and Optimization

JModelica has three dierent types of model objects FMUModel, JMUModel and CasadiModel. These objects can be used for simulation and optimization. In this project only the rst two object types are used and will therefore be described briey.

The main dierences between a FMUModel and an JMUModel, is that only the JMUModel can be used for optimization, since it supports the Optimica attributes. A JMUModel can only be created in the JModelica.org environment and used together with Modelica and Optimica code[Modelon AB,2012]. FMUModels, on the other hand, can be created from other modelling and simulation environments since it follows the FMI(Functional Mock-up Interface) standard. The FMUModel is used to simulate the model and the result is later used as an initial guess for the optimization problem, which is handled by an JMUModel. To begin with the functions compile_fmu and load_fmu are imported in the Python script le.

from pymodelica import compile_fmu

from pyfmi import load_fmu

init_fmu = compile_fmu ('simModel ',[' Heatexch_opt .mo ',' Water_poly ']) init_model = load_fmu ( init_fmu )

init_res = init_model . simulate ( start_time =0, final_time =4000)

The function compile_fmu needs the model name ('simModel')and the Modelica le ('Hea-texch_opt') as input. Additional packages that contain functions or models needed for the simulation, can be added as arguments. In the examples shown in this section, the media model for water ('Water_poly') is added in this way. Finally the model can be simulated, by calling the function simulate with start and nal time. Additional options for compilations and simulations of FMUs can be found in the JModelica.org Userguide [Modelon AB,2012].

The next step is to nd an optimal solution to the problem. As previously, the functions needed are imported and a compilation is done. The function compile_jmu needs an additional argument, which is the Optimica le ('heatopt.mop').

from pymodelica import compile_jmu

from pyjmi import JMUModel

jmu_name = compile_jmu (" heatopt . Opt1 ",[' Heatexch_opt .mo ','heatopt . mop ',' Water_poly '])

opt_result = JMUModel ( jmu_name ) res = opt_result . optimize ()

For the optimization problem in this report it was necessary to dene optimization options, to improve the convergence of the optimization algorithm. The options are dened by using the function optimization_options and an example in which the number of elements n_e are changed is shown below.

opts = opt_result . optimize_options () opts ['n_e ']=20

res = opt_result . optimize ( options = opts )

The options used for this thesis project can be found in Table 2 and all available options can be found in the JModelica User Guide [Modelon AB,2012].

IPOPT

The JModelica platform uses a simultaneous optimization algorithm with based orthogonal collocation using Lagrange polynomials on Radau points [JModelica.org,2013]. The process

(22)

Table 2: Optimization options for JMUs

Options Description

n_e Number of elements that the optimization interval is divided in. The default number of elements is 50.

n_cp Number of collocation points for each element in n_e. The default number of collocations points is 3.

blocking_factors Blocking factors are used to keep the control prole con-stant for nearby elements. Blocking factors are not used in the default case, but can improve the convergence of the op-timization algorithm. The blocking factors [2,1,6,1] makes u0 = u1, u3 = u4 = u5 = u6 = u7 = u8, while u2 and u9 is

not dependent on nearby elements. The sum of the elements must equal the number of elements n_e.

result_le_name Writes results to a text le.

init_traj An initial trajectory is used as a rst guess for the optimiza-tion. problem

Table 3: Optimization options for the IPOPT solver

Options Description

tol Decides to convergence tolerance. If the (scaled) NLP error is smaller than the value of tol, the optimization algorithm will be successfully terminated. For this project 5 × 10−4 _is

used and the default values is 10−8_.

max_iter Maximum iterations allowed for the optimization problem. In this project 100-300 is used and the default value is 3000.

makes it possible to transfer dierential algebraic equations (DAE) used in the dynamic model, into a large-scale NLP problem [Åkesson, Gäfvert & Tummescheit 2009]. The NLP problem is solved with the IPOPT(Interior Point OPTimizer) solver [IPOPT,2013], which is an open-source solver appropriate for large-scale problems. There are several options for the IPOPT solver that can be changed, to suit ones optimization problem. The options are changed in the Python script le, in the following way:

opts [' IPOPT_options ']['tol ']=5e -4

In Table 3 descriptions of the IPOPT options used in this thesis project is presented. In this project the IPOPT solver MA27 have been used, instead of the standard solver MUMPS. MA27 is better at handling large-scale problems and is also more robust.

3.5 Media Models

In Modelica it is convenient to use media models, to describe the behaviour of the dierent media used in the model. Media models are very useful when working with uid models, because they can describe the non-linear behaviour of the uid in a good way. For thermodynamic systems a typical problematic process to describe is when water is evaporated to steam, because the water

(23)

exists in two phases.

In the Modelica standard library you can nd media models for air, water, ideal gases, incom-pressible media and comincom-pressible liquids. If the medium only consists of one substance you usually need two states, to derive the thermodynamical properties. If your medium consists of several substances, you might need the composition of the medium to use media models. Some of the more common gas mixtures, for example air, has a default composition that is used if the user doesn't dene a composition.

The media models combine equations, tables and substance specic constants to dene the medium properties. It is therefore important to verify that your model operates in a region that is supported by the media model.

In this project it has been crucial to have a good media model for water. The standard model in Modelica WaterIF97, is unfortunately not supported by the optimization algorithms in JMod-elica. Instead the media model Water_poly developed by Modelon is used, where polynomial expansions around the phase limit is used to calculate the thermodynamical properties. Wa-ter_poly is not as extensive as WaterIF97, but sucient to describe the variables and constants needed for a simple model of a heat exchanger.

(24)

4 Modelling

The complex model of Jänschwalde is presented in this chapter and the important components are pointed out. Further on the components developed within this thesis project are described and the complete start-up model is shown at the end of the chapter.

4.1 Complex model of Jänschwalde

Figure 9: An overview of the complex model of Jänschwalde, owned by Vattenfall. The complete steam cycle is included in the model. The water is heated, used to drive the turbine, reheated and nally condensed and cooled. The total number of equations and unknowns is 16,927.

In Figure 9 a complex model of a block in the Jänschwalde power plant is shown, which was initially developed by Modelon in a previous Research and Development project. The model has been further developed by Vattenfall Research and Development in collaboration with the University of Rostock. The model consists of 16,927 unknowns and therefore requires a great amount of simulation time. The model has been validated against measurement data from the real plant and managed to describe the process well in steady-state, especially when 100 % load was used as input [Saarinen, Boman & Funkquist, 2011]. The model can not be used for start-up optimization because of it's complexity, but has been used to validate components and simulation models developed within this thesis project.

Since the goal of this thesis project was to nd a model that could be used for start-up op-timization, the main focus was to develop a model with less complexity. The strategy was to:

• Only include the components that are most important for the start-up. • Reduce the complexity of the components to achieve a short simulation time. • Use simpler models to describe media.

(25)

The furnace section in Figure 11 includes the combustion process and the evaporator. The model includes two combustor components, since the combustion is incomplete in the rst component. In the evaporator water is evaporated to steam, which continues on to the superheater section. Most components that occur in the start-up model are from the superheater section, shown in Figure 10. The purpose of the superheater section is to gradually heat up the steam so it can drive the turbines. The steam from the evaporator enters the component SHS, but then it continues to the section ECO and SH1, shown in Figure 12. After SH1 the steam goes back to the superheater section. The economizer is not included in the start-up model, because it is not considered crucial for the start-up scenario.

Figure 10: The superheater section is one subcomponent from the model shown in Figure 9. The gure shows four superheaters (SH) and three reheaters (RH) from the complex model of the Jänschwalde power plant. The circles represent connectors for ue gas (green), water (blue) and coal (brown), and are used to link the component to other submodels or boundary conditions.

(26)

Figure 11: The furnace section in Figure 9, which contains the combustion process and the evaporator.

(27)

4.2 Components

The modelling has been done in Dymola, with the language Modelica 3.2. The majority of the components used have been developed within this thesis project. The simulation models include some components from Blocks in the Modelica Standard Library as well. A summary of the components can be found Appendix A.

4.2.1 Connector

A very basic component that is needed for modelling is a connector. Connectors are used to link dierent components together. In this project two types of connectors have been created called owport and owport_heat. An example of owport, which is used to transport uids and gases is shown below.

connector flowport "1-D transfer of fluid / gas "

flow Water_poly.Units.MassFlowRate m_flow " Mass flow rate ";

stream Water_poly.Units.SpecificEnthalpy h_outflow " Specific Enthalpy ";

Water_poly.Units.Pressure p " Pressure ";

end flowport ;

The connector contains the three variables mass ow rate m_flow, specic enthalpy h_outflow and pressure p. The package Units in Water_poly is used to describe the variables and the type MassFlowRate is shown as an example below.

type MassFlowRate =

Modelica.SIunits.MassFlowRate( min =-1000, max =1000, nominal =100) ;

It is advantageous to use Units because it helps the user to detect errors in equations. A warn-ing will be given durwarn-ing the compilation phase if there are any dimensional errors in the model, which could imply incorrect equations. The nominal value is used by the optimization algorithm for scaling, and can be modied to suit ones optimization problem. For dimensionless variables the declaration Real can be used.

Another interesting thing to observe is the usage of the type stream. In older versions of Modelica only two types ow and potential have been used to describe interaction between components. The old approach can still be used for many types of systems, for example electrical and mechan-ical systems, but is insucient when dealing with thermodynamic systems [Franke et al., 2009]. The main problem is that the ow can be bi-directional in uid models and that the enthalpy is dependent on the ow direction. With stream variables it is possible to formulate the balance equations, so that the result is dependent on the ow direction. Another approach could be to introduce Boolean unknowns to represent the ow direction, but then the system would be non-linear and singular around zero ow [Franke et al., 2009]. An example of how stream variables are used is presented in section 4.2.4.

4.2.2 Boundary conditions

To be able to simulate the start-up model one needs to dene boundary conditions for the water and gas side, also called sinks and sources when working in Modelica environments. Two com-ponents are used, called MassFlowBoundary and PressureBoundary.

For MassFlowBoundary, four parameters must be specied, mass ow rate, specic enthalpy, pressure and hydraulic conductance. The component also includes Boolean parameters, which

(28)

gives the user a possibility to choose a nominal value or an external signal as input. The code associated with the mass ow rate is demonstrated below as an example.

parameter Water_poly.Units.HydraulicConductance G=0 " HydraulicConductance ";

parameter Water_poly.Units.MassFlowRate m_flow0 =0 " Nominal mass flowrate "

annotation ( Dialog ( enable =not use_mdot_in ));

parameter Boolean use_mdot_in =false

"If true, mass flow rate is set from input connector mdot_in else from m_flow0 "

annotation ( Evaluate = true,HideResult = true,choices ( __Dymola_checkBox =true), Dialog ( group =" Options "));

Modelica.Blocks.Interfaces.RealInput mdot_in if use_mdot_in ;

protected

Modelica.Blocks.Interfaces.RealInput m_flow ;

equation

port.m_flow = - m_flow + ( port.p - p0)*G;

if not use_mdot_in then

m_flow = m_flow0 " Flow rate set by parameter ";

end if;

connect( m_flow, mdot_in );

In the graphical interface the Boolean variable use_mdot_in, can be changed to "True" by using a check box. In that case the connector RealInput mdot_in will appear, which can be connected to a another component that describes the signal. Some useful components that can be used to evaluate the dynamic response of the system is Ramp and Step from Sources in the Modelica Standard Library. In the equation section you can see that m_flow is equal to mdot_in, if the Boolean variable is True. Otherwise m_flow is given by the parameter m_flow0, which can be changed in the code or in the graphical interface.

The example includes a parameter G, that represents the hydraulic conductance. This parame-ter was introduced to relax the boundary conditions and makes larger models easier to initialize. The downside is that mass ow rate that exits the component (port.m_flow) can deviate from the value dened in the component. However, if a small value of G is used (≈ 0.01) the deviation is often acceptable.

The structure of the component PressureBoundary is similar to MassFlowBoundary, but it only has three variables, pressure p, specic enthalpy h_outflow and hydraulic resistance R. In this case the hydraulic resistance is used to relax the boundary conditions.

4.2.3 Combustor

The combustion process in a power plant is very complex and depends on many dierent factors, of which only some can be measured. One complication is that coal, air and ue gas consists of several dierent substances. In the complicated model this problem was solved by introducing transformation matrices between fuel and ue gas. The transformation matrices insure that the chemical reactions are fullled when fuel and air are burnt. A very simple example is the chemical reaction C + O2− > CO2, which requires one mole coal C and one mole oxygen O2

to produce one mole carbon dioxide CO2. If there is a lack of oxygen, only a part of the coal

(29)

(a) The gas temperature is used as input and gas

enthalpy is used as output. (b) The mass ow rate of coal is used as input andmass ow rate of ue gas is used as output. Figure 13: The gure shows how well the process models P1, P1D, P0 and P0D, generated in the System Identication Toolbox in Matlab, ts the reference data.

combustion process, a very accurate description of the reality is obtained. The downside is that the matrices introduce large systems of equation, because the combustion consists of several chemical reactions. In the complex model seven chemical reactions are used, to give an approxi-mate description of the real process. This is still a bit too complex for the start-up model, since it requires much computational time.

Instead a simple model of a combustor was developed by using the System Identication Toolbox in Matlab with data from the complex model of the Jänschwalde power plant. Two sets of data were used to identify the relationship between ue gas temperature and enthalpy, see Figure 13a, and mass ow rate of coal and ue gas, see Figure 13b. The data sets were divided in estimation and validation data and then dierent process models in Table 4 were used for identication.

Table 4: Process models used in the System Identication Toolbox

Model Description P0 K P0D K exp(−T d · s) P1 K 1+T p1·s P1D K exp(−T d·s) 1+T p1·s

When the four dierent process models were tried out, it turned out that a constant relation between input and output gave the best performance for both cases. The two dierent constants obtained by the system identication were used in the combustor component together with typical equations for a boundary component described in previous section.

(30)

4.2.4 Heat exchanger

The component heatexchanger was built to represent the superheaters (SH), reheaters (RH) and the evaporator in the plant. The equations used are described by (4) to (7), in Section 2.2.1. The component consists of four connectors water_in, water_out, gas_in and gas_out, and each component has three variables representing mass ow rate, pressure and enthalpy.

The dynamics of the component are given by the balance equations for mass and energy. The relevant code for the waterside and the wall component is given below.

Water_poly.Units.SpecificEnthalpy h " Enthalpy of water ";

Water_poly.Units.Pressure p " Pressure of water ";

Water_poly.Units.Density rho = Water_poly.Functions.density_ph( p,h ) " Water Density ";

Real ddhp = Water_poly.Functions.ddhp_ph( p,h )

" Water density derivative w.r.t h at constant p";

Real ddph = Water_poly.Functions.ddph_ph( p,h )

" Water density derivative w.r.t p at constant h";

equation

// water side

V*( ddhp *der(h)+ ddph *der(p))= water_in.m_flow + water_out.m_flow ; // mass balance

V*(der(h)*( rho +h* ddhp ) + der(p)*(h* ddph-1 ))= water_in.m_flow *

actualStream( water_in.h_outflow )+ water_out.m_flow *actualStream( water_out.h_outflow )- Q_L ; // energy balance

water_in.p = water_out.p ;

// wall between gas and water side M* c_p_wall *der( T_wall )= Q_L + Q_G ;

Q_L = A_wall * alpha_wall_water *( T-T_wall );

initial equation

p=p0;

T_wall = T_wall_0 ; h=h0;

This component includes initial conditions for the states pressure p, specic enthalpy h and wall temperature T _wall. It also includes the constant parameters volume of water tank V , mass of metal wall M, specic heat capacity of metal wall c_p_wall, area of metal wall in contact with water A_wall and coecient of heat transfer alpha_wall_water. Functions from W ater_poly were used to evaluate the density ρ, and the partial derivatives of ρ. For the ue gas side a static mass and energy balance was used, because there were no good media models simple enough for the optimization phase available.

In the energy balance the usage of stream variables is demonstrated. The operator actualStream means

actualStream( water_in.h_outflow ) ==

if water_in.m_flow > 0 then inStream( water_in.h_outflow )

else water_in.h_outflow

according to [Franke et al., 2009]. If the inStream operator is dened in the following way

inStream( water_in.h_outflow )= h_water_in ; water_in.h_outflow = h_water_out ;

(31)

you get an energy balance which is correct also if the ow shifts direction. The ow entering the component is by convention dened as positive, while the ow leaving a component is dened as negative.

In the heatexchanger model one can also observe that the pressure is equal for the uid en-tering and exiting the component. The pressure drop is instead handled in the component V alve_linear. Except for inducing a pressure drop, the component is needed when connecting several heat exchangers together, so that all ow variables entering a heat exchanger have the same sign.

4.2.5 Header and Separator

To calculate stress levels in the header and separator two components are used. The rst com-ponent is called Volume and is used to describe the dynamics of the water. The comcom-ponent is similar to the water side of the heatexchanger component, but has a connector called ow-port_heat that is connected to the Wall_component.

The Wall_component is a cylindrical component, discretized in the radial direction. The tem-perature T and heat ow rate Q_flow, from the connector are used as boundary conditions for the inlet of the cylinder. The components is discretized in two elements and Fourier's equation is used to describe the temperature distribution. The stress calculations from Section 2.2.2, is used to calculate stress that occur because of pressure and thermal changes.

The model library developed within this project also contains a more detailed version of a separator, inspired by the component BensonSeparator developed for a library called Ther-malPower(version 1.5). Unfortunately it could not be used for optimization purposes, since it made it hard for the problem to converge.

4.2.6 Turbine

In the turbine steam is expanded, which causes a drop in pressure for the medium. The process can be described by the Stodola equation [Stodola, 1945], which relates the mass ow rate of the steam m_flow_in, the steam temperature T _in, and pressure for the entering and exiting steam p_in and p_out.

m_flow_in =Kt*Heatexch_opt.Functions.regRoot(( p_in ^2- p_out ^2) / T_in );

Kt is a turbine specic constant, that can be estimated from experimental data. The function regRoot(x)from the Modelica standard library is used instead of sqrt(x), to calculate the square root. regRoot(x) is an approximation of sqrt(x), that makes sure that the derivative of √x is nite and smooth in x = 0.

4.3 Model used for optimization

Figure 14 shows the model that is used for start-up optimization. The components used in the start-up model include one evaporator, ve superheaters, three reheaters, two volume compo-nents, two wall compocompo-nents, one coal combustor, one sink for ue gas, two sources and sinks for water. A description of these components is given in Appendix A. The biggest dierences between the start-up model and the actual plant are:

• The separator component has been simplied. The real plant has a separator with one inlet and two outlets for water and steam. The purpose of the separator is to split steam and

(32)

Figure 14: Model used for start-up optimization. The heat exchangers modelled are the superheaters (SH), reheaters (RH) and the evaporator (Evap). Other components include wall and volume components to calculate stress levels in headers and separators and components to represent boundary conditions. The model has one input du connected to the integrator component which describes how the load is changed du

dt. There are three outputs to the right for live steam temperature, header stress and separator

(33)

water, so that the steam continues to the rst superheater and water is recirculated back to the evaporator. In the start-up model the separator only has one outlet for steam/water, because it made the optimization problem easier to solve.

• Thermal radiation is not modelled. In reality the high temperatures cause radiation, that aects the lower components in the boiler, mainly the evaporator and SHS.

• Spray water is not modelled, which is used to control the steam temperatures between superheaters.

• Heat losses to surroundings are not included in the model.

• System Identication techniques have been used to simplify complex components, like the combustor.

• Time delays do not exist in the start-up model

• There are no material models for metals, so constants are used to describe the material properties in walls.

(34)

5 Validation

5.1 Superheater

A model of a heat exchanger from Section 4.2.4 has been compared to a superheater from the complex model. The complex superheater consists of 1,082 equations and unknowns and the simplied start-up superheater has 29 equations and unknowns. In Figure 15 the steam tem-peratures are presented and in Figure 16 gas temtem-peratures are shown, for both components. To compare the dynamics of the components two step tests were done, rst the ue gas temper-ature was increased by 200◦_C _{at T =500 s and secondly the mass ow rate of the ue gas was}

increased by 20 % at T =1000 s. The same boundary conditions were applied for pressure and en-thalpy for the two components. The rst dierence one can observe is a temperature dierence of 2◦_C _{for the water entering the components. This eect probably occurs because dierent}

media models are used for water, which gives slightly dierent temperatures even though same enthalpy and pressure are used for the incoming water. A dierence in temperature can also be seen on the gas side, because of dierences in the gas medium models.

Figure 15: Temperatures of in- and outowing water, for the complex superheater and the start-up superheater.

When the temperature is changed at T =500 s, the systems behave similar and the temperature of the out-owing water is increased by approximately 20◦_{C. The complex system has a shorter}

rise and settling time. At T =1000 s, when the mass ow rate of the ue gas is increased, a much larger response is registered for the complex component. At rst glance you would think that the start-up model is insucient for describing these kind of changes. In Section 5.2, it is shown that when several of the simpler heat exchangers are connected together, they can describe the

(35)

Figure 16: Temperatures of in- and outowing ue gas, for the complex superheater and the start-up superheater. The temperature curves for inowing ue gas are equal for both superheaters.

complex model owned by Vattenfall quite well. In the complex model the steam temperature is reduced by using spray water. It shows that the start-up model, without spray water, is a sucient description of the large model, including spray water.

It was possible to change parameters in the start-up model, so that a higher response was re-ceived when the mass ow rate changed, but then the start-up model could not be used to describe the real plant. Some explanations for the deviations between models include dierent media and material models and that the wall in the complicated model is discretized.

5.2 Start-up model

To validate the complete start-up model, a steady-state simulation for the full load case was done and compared with reference data from the plant and the result is shown in Table 5. Data from the actual plant have been used to represent the boundary conditions for pressure and mass ow rate. Boundary conditions for enthalpy have been picked from a steady-state simulation of the detailed model.

The overall performance of the model is considered good, since the relative errors are small, which can be seen in Table 5. The large dierence for the component SHS, is probably because thermal radiation is neglected.

(36)

Table 5: Steam temperatures out of each heat exchanger component at a steady-state simulation with full load. Heat exchangers SHS SH1 SH2 SH3 SH4 RH1 RH2 Reference ◦_C ₃₇₃ ₄₀₆ ₄₅₁ ₄₉₂ ₅₃₄ ₄₈₁ ₅₄₀ Start-up model ◦_C ₃₅₉ ₄₁₃ ₄₅₂ ₄₉₆ ₅₃₆ ₄₈₃ ₅₄₄ Absolute error -14 7 1 4 2 2 4 Relative error -3.8 % 1.7 % 0.22 % 0.81 % 0.37 % 0.41 % 0.74 %

Figure 17: Comparison between reference data, complex simulation model and start-up model at steady-state with full load.

temperatures from a full load simulation of the complex model. The conclusion is that the boiler of the power plant can be fairly well represented by the start-up model. There exists some devia-tions in steam temperatures for both the complex model and the start-up model, but both models can be used to simulate the behaviour of the real process. The start-up model consists of 457 unknowns and equations and is therefore possible to use for optimization purposes in JModelica. In [Saarinen, Boman & Funkquist, 2011] ramp tests were done to compare the complex simula-tion model of Jänschwalde with measurement data from the plant, see Figure 18. The systems were ramped from 100 % load to 60 % load during 19 minutes. The system was given time to stabilize at 60 % load and then ramped up again to full load during 20 min. In Figure 19 the same ramps were applied to the start-up model, to ramp up and down the coal ow. The feed-water ow rate was also ramped to replicate the real experiment.

When looking at the system responses for the complex model and the start-up model, it shows that both systems respond immediately to load changes, when live steam temperatures are

(37)

eval-Figure 18: Comparison of the complex simulation model of Jänschwalde (dotted lines) and measurement data (solid lines) from the actual plant, when the system was ramped from 100 % load to 60 % and back again to 100 % load.

uated. The measurement data shows that the real system has a time delay of approximately 20 min when the system is ramped down. However the response is not as slow for the real sys-tem, when the load is ramped up again. If the delay is neglected all three systems have about the same rise time. The settling time can not be evaluated for the real system, since the steam temperatures uctuate because of disturbances in the real plant.

For all systems the steam temperatures in SH1, SH2, SH3 and SH4 decreases when the load is ramped down. For SHS the steam temperature increases for the start-up model, while it is approximately constant in Figure 18. When the load is ramped up, the complex model has an overshoot that is larger compared to the start-up model and the measurement data. The superheater component also had a very large step response compared to the start-up model, when the mass ow rate was changed in Section 5.1. It is therefore possible that the response is too high for the superheater components, which causes some of the overshoot seen in Figure 17.

(38)

Figure 19: Simulation of the start-up model when the system was ramped from 100 % load to 60 % and back again to 100 % load. The start-up model is not aected by disturbances, compared to the systems shown in Figure 18.

(39)

6 Optimization

6.1 Problem specications

The optimization problem is dened by an objective function, constraints and an optimization interval. The objective function ϕ is dened by a quadratic cost function, which is minimized over an optimization interval. Since ϕ should be minimized, deviations between the live steam temperature TSH4 and the reference value Tref are penalized. Changes in the control variable

du

dt are also penalized, since it is desirable to reach steady-state ( du

dt = 0) quickly. The objective

function is dened as ϕ = Z tf t0 α(TSH4− Tref)2+ β du dt 2 dt, (20)

with α ≥ 0, β ≥ 0. The constraints are given by the stress levels in the header and the separator σmin,header ≤ σheader ≤ σmax,header, (21)

σmin,separator ≤ σseparator ≤ σmax,separator. (22)

Calculations of σheader and σseparator are done in the Dymola model shown in Figure 14. If one

has dened minimum and maximum attributes for the variables in the Dymola model, these attributes will also be used as constraints in the optimization in JModelica.

The objective function

ϕ = Z t_f

t0

dt, (23)

was also considered for the optimization problem, with the constraints dened by (21) and (22). The objective function (23) actually minimizes the nal time tf, while (20) minimizes deviations

from the steady-state values and changes in the control signal over the interval. However, it was not possible to get the IPOPT solver to converge with the objective function in (23).

6.2 Optimization options

The optimized solution is highly dependent on the discretization, which is decided by the num-ber of nite elements neand collocation points ncp used. A more rened discretization generally

improves the accuracy of the solution, but also lengthens the computation time. As mentioned in section 3.4, blocking factors can be used to help the IPOPT solver converge to a solution. In this project blocking factors are used to keep the control signal du

dt constant for several elements.

An example is shown in Figure 20, with the blocking factors [1,1,2,6]. The strategy for this project was to use ne = 10 and lock the control signal for several elements for the rst rounds

of the optimization. The number of elements ne was increased gradually and the control signal

was allowed to change more freely as the optimization continued.

The optimization result also depends on the parameters α and β in the cost function described by (20). In Chapter 7 some dierent set-ups of the parameters are tried out. The values used for the reference value Tref and the constraints in (21) and (22) are presented in Table 6. With

the constraint boundaries used, the materials should be able to handle at least 105 _start-ups

(40)

Figure 20: Plot of coal ow u and the control signal du

dt, with blocking factors [1, 1, 2, 6].

Table 6: Specication of values used for Tref and constraints.

Parameter Value Tref 806 K σmin,header -200 MPa σmax,header 400 MPa σmin,separator 800 MPa σmax,separator 1400 MPa 6.3 Initialization

It was very hard for the IPOPT solver to converge to a solution, without using an initial trajec-tory as a start-guess. To solve this problem, a simulation of the model was done with constant load. It is important that the initial guess is within the constraints dened by the optimization problem, so that is why a steady-state simulation was used. The nal values of the state-variables were used as initial values for the optimization problem. In this way, some of the large transients that occur during the initialization phase can be avoided.

The input trajectory from the simulation problem was used as a starting guess. It is easier for the optimization algorithm to converge to a solution, if the starting guess is close to the optimal solution. That is why numerical optimization techniques often consist of an iterative process. Generally it is easier to start with a very rough discretization. If an optimal solution is found, it can be used to initiate the same optimization problem, with a rened discretization.

(41)

7 Results

Variations of weight parameters in cost function

The rst step was to investigate how dierent values of α and β in the objective function (20) aected the outcome of the optimization. The blocking factors were set to [1,1,2,6] and then the parameters α and β were varied for the dierent problem set-ups. For the 6 experiments, the same optimization and IPOPT options were used: ne= 10, ncp = 3, max_iter = 300 and

tol = 5e − 4.

Table 7: Summary of parameter values used for dierent problem set-ups.

Experiment α β Final value of TSH4≈ Tref Number of iterations

1 10−2 10−1 N/A N/A 2 10−2 4 × 10−2 Yes 39 3 10−2 10−2 Yes 41 4 1.6 × 10−3 10−2 No 40 5 2.5 × 10−3 10−2 No 34 6 4 × 10−4 10−2 No 40 7 10−4 10−2 No 36

When the parameter settings for Experiment 1 was used from Table 7, the maximum number of iterations was reached by the IPOPT solver, so the optimization algorithm was terminated unsuccessfully. The IPOPT solver managed to converge to a solution for the other experiment set-ups and the number of iterations needed was between 36 and 41. The total computation time for those cases was between 84 and 98 seconds.

The result for Experiment 3 is shown in Figure 21. The steady-state value for the steam tem-perature in SH4 (Tref) was reached at T ≈ 800 s. The stress amplitude for the header was ≈

150 MPa, while the stress amplitude was very small for the separator. The maximum coal ow was very high, 43 % over the steady-state value at T =200 s.

In Figure 22 the optimization result is shown for Experiment 6 in Table 7. The steady-state value Tref for the temperature is not reached during the optimization interval of 2000 s. The

stress amplitudes are smaller, since the steam temperature increases at a smaller pace.

An interesting result was obtained for Experiment 5. The IPOPT solver converged to a solu-tion, with very high temperature levels, which are shown in Figure 23. The nal temperature of the steam exiting superheater SH4, is 234 K over the reference temperature Tref. The stress

Development of a dynamic model for start-up optimization of coal-red power plants