Supervisory control of a combined heat and power plant by economic optimization

(1)

Hamburg University of Applied Sciences Faculty of Life Sciences

SUPERVISORY CONTROL OF A COMBINED HEAT AND

POWER PLANT BY ECONOMIC OPTIMIZATION

MUSTAFA GÖKSEL DELİKAYA

A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF

THE REQUIREMENTS FOR THE DEGREE OF

MASTER OF ENGINEERING (M.ENG) IN

RENEWABLE ENERGY SYSTEMS BY

Submitted on the 31st of March 2015

This master thesis has been completed at the Fraunhofer Institute for Solar Energy Systems (ISE) in Freiburg, Germany.

Supervisors

Prof. Dr.-Ing. Gerwald Lichtenberg (HAW Hamburg) Mehmet Elci, M.Sc. (Fraunhofer ISE)

(2)

(3)

iii

ABSTRACT

Using combined heat and power (CHP) units within district heating systems (DHS) has been an effective way of meeting residential energy demand in Germany. Generally speaking, electricity fed in to the grid by the CHP is usually sold at a fixed price in today’s electricity market. Assuming that the share of renewable energies will be higher in the near future, it can be anticipated that the electricity prices will highly fluctuate due to the uncertainties within the renewable energy sources, such as wind speed and solar irradiance. Therefore, control mechanisms for heat and power producing plants are expected to switch their operation strategy from heat-driven to power-driven operation. A power-driven operation makes sure that the CHPs are shut down when the electricity market is not competitive enough to produce electricity. In this master’s thesis, a power-driven operation is achieved through an economic optimization. The optimization problem, which is formulated as a discrete optimization problem, is to find out the best ON/OFF operation trajectory of the units involved in a DHS; namely a CHP, a boiler and a storage tank. A simplified model capturing the power-based dynamics of a physical DHS model is implemented at simulation and modeling tool Dymola (Dynamic Modeling Laboratory). Optimization tool GenOpt (Generic Optimization Program) with particle swarm optimization (PSO) algorithm is used to solve the discrete optimization problem. The implementation of the model is verified by several test cases. Finally, a future scenario of the year 2023 is approximated in order to compare the financial gains and grid interactivity of the power-driven and the heat-driven operation. In addition, the effect of varying the storage size on plant gains and grid interactivity is investigated and discussed.

(4)

(5)

v

ACKNOWLEDGEMENTS

I am deeply grateful to my supervisor Prof. Dr.-Ing. Gerwald Lichtenberg for taking me his student and providing me all the supervision and encouragement I needed. His guidance was mentoring me on my first steps at my academic career. I especially thank him in making me realize the keywords and the underlying problem of this master’s thesis, which were not clear in the beginning. It always amazed me how he approached to the problems and proposed creative solutions. It has always been enlightening to have a conversation with him. I would like to thank all his advices on academic writing and corrections for my preliminary report as well.

In addition, I would like to thank to my second supervisor PhD student Mr. Mehmet Elci for letting me involved in such an interesting project. I would like to thank him for spending his valuable time with me anytime I needed his supervision. It was very inspiring to have discussions and trying different approaches to overcome the problems that are faced during this study. I am also thankful for all his corrections, which helped me a lot to give the final form of this thesis.

I would like to thank to PhD student Mr. Kai Kruppa for his advices and corrections, which helped me to structure the modeling chapter of this thesis in a more understandable way.

Lastly, I would like to thank to my colleagues, Ms. Sunah Park, Mr. Tom Cordes, Ms. Friederike Rautenberg, Mr. Kai Iking, Mr. Marc Eisenbarth and Mr. Max Walch for the great working environment, the friendship and their emotional support throughout my work.

(6)

(7)

vii

List of Figures

Fig. 1.1 State trajectory of a DES ... 1

Fig. 1.2 Ramade Wonham (RW) Framework ... 1

Fig. 2.1 Schematic representation of a DHS ... 7

Fig. 2.2 Site plan of District “Gutleutmatten” ... 8

Fig. 3.1 Single node representation of the storage tank ... 12

Fig. 3.2 Energy flow within the system boundary ... 17

Fig. 3.3 Input-outputs and disturbances of the system ... 20

Fig. 4.1 Implementation of the system at “Dymola” ... 21

Fig. 4.2 Implementation of the minimum operation time constraint ... 22

Fig. 4.3 Representation of the storage tank with a single node ... 23

Fig. 4.4 Operation range of the storage tank ... 24

Fig. 4.5 Hysteresis band for the storage tank ... 24

Fig. 5.1 Coupling a simulation program with GenOpt ... 27

Fig. 5.2 PSO algorithm ... 30

Fig. 6.1 Load and prices for Test Case 1 ... 32

Fig. 6.2 Optimization results of Test Case 1 ... 32

Fig. 6.3 Model with 8 input variables ... 33

Fig. 6.4 Model with 2 input variables ... 34

Fig. 6.5 Results of enumeration and optimization ... 35

Fig. 6.10 Optimization results of 1st half-day (1) ... 41

Fig. 6.11 Optimization results of 1st half-day (2) ... 41

Fig. 6.12 Illustration of the constraint violation ... 42

Fig. 6.13 Optimization results of whole day (1) ... 44

Fig. 6.14 Optimization results of whole day (2) ... 44

Fig. 6.16 Optimization results of 1st case ... 46

Fig. 6.17 Optimization results of 2nd case ... 47

Fig. 7.1 Decision Tree for selection of the operating states at heat-driven operation ... 50

Fig. 7.2 Operation state distribution in the tank with respect to temperature levels ... 51

Fig. 7.3 Electricity prices and residual load of year 2011 ... 52

(10)

x

Fig. 8.1 Results of the future scenario (2) ... 54 Fig. 8.2 Comparison of LGMC_abs values ... 57 Fig. 8.3 Optimization results of the same summer scenario with different storage sizes ... 58

(11)

xi

List of Tables

Table 3.1 Operating States ... 11

Table 3.2 Parameters of the system ... 13

Table 3.3 Disturbances of the system ... 14

Table 3.4 Outputs of the system ... 14

Table 3.5 Active equations depending on the operating state ... 16

Table 6.1 Optimization parameters for Test Case 1 ... 31

Table 6.3 Optimization Parameters for Test Case 3 ... 38

Table 6.6 Comparison of overall gains ... 47

Table 7.1 Power Capacities ... 51

Table 8.1 Optimization parameters for comparison scenario ... 53

Table 8.2 Parameters values for initialization ... 53

Table 8.3 Result of the future scenario (1) ... 54

(12)

(13)

xiii

List of Abbreviations

Abbreviation

Description

BB Branch and Bound

CHP Combined Heat and Power

DE Differential Evolution

DES Discrete Event Systems

DHS District Heating System

Dymola Dynamic Modeling Laboratory

GA Genetic Algorithm

GenOpt Generic Optimization Algorithm

GLM Gutleutmatten (name of a project)

HVAC Heating Ventilation and Air Conditioning

IP Integer Programming

ISE Institute for Solar Energy Systems

KP Knapsack Problem

MIP Mixed Integer Programming

PSO Particle Swarm Optimization

SA Simulated Annealing

TRNSYS Transient Systems Simulation Program

(14)

(15)

xv

List of Symbols

Symbol

Unit

Description

𝐴 𝑚2 _{Heat Transfer Area}

𝑏CHP Boolean signal showing the CHP’s state (on/off)

𝑏Boiler Boolean signal showing the Boiler’s state (on/off)

𝑏St Boolean signal showing the Storage Tank’s state (on/off)

𝑐c Cognitive velocity constant

𝑐el 𝑐𝑐. 𝑘𝑘ℎ⁄ Electricity price

𝑐fuel 𝑐𝑐. 𝑘𝑘ℎ⁄ Fuel price

𝑐heat 𝑐𝑐. 𝑘𝑘ℎ⁄ Heat price

𝑐p 𝐽 (𝑘𝑘 ∙ 𝐾⁄ ) Specific heat capacity

𝑐s Social velocity constant

𝑑𝑑 𝐽 Change of internal energy

𝑑𝑑k 𝐽 Change of kinetic energy

𝑑𝑑e 𝐽 Change of potential energy

𝐺 𝑘 Residual Load

𝐺̅ 𝑘 Mean residual load over an evaluation period

𝐺k(𝑡) the 𝑘th component of the best position vector ever reached

𝐻p Prediction horizon

𝐽CHP 𝑐𝑐. Gain of the CHP

𝐽Boiler 𝑐𝑐. Gain of the Boiler

𝐽St 𝑐𝑐. Gain of the Storage Tank

𝐽T 𝑐𝑐. Total System Gains

𝐿𝐺𝐿𝐿abs Absolute load grid matching coefficient

𝑛i Number of input variables

𝑝 Number of time steps

𝑑el 𝑘 Electrical power production

𝑑ik(𝑡) the 𝑘th component of the best position vector of particle 𝑖

𝑑load 𝑘 National heat load

𝑑PV 𝑘 Power production by photovoltaic systems

𝑑wind 𝑘 Power production by wind turbines

(16)

xvi

𝑑2 Penalty value for the temperature constraint

𝑄 𝐽 Heat exchange within a control mass

𝑄ch 𝑘 Heat charged to the storage tank

𝑄CHP,th 𝑘 CHP nominal thermal power

𝑄CHP,el 𝑘 CHP nominal electrical power

𝑄disch 𝐽 Heat discharged from the storage tank

𝑄fuel,CHP 𝑘 Amount of fuel consumption at CHP

𝑄fuel,Boiler 𝑘 Amount of fuel consumption at Boiler

𝑄loss 𝐽 Amount of heat loss

𝑄loss,Boiler 𝑘 Amount of heat loss at Boiler

𝑄loss,CHP 𝑘 Amount of heat loss at CHP

𝑄loss,St 𝑘 Amount of heat loss at Storage Tank

𝑄̇loss 𝑘 Amount of heat loss per unit time

𝑄net 𝐽 Net heat transfer

𝑄̇net 𝑘 Net heat transfer per unit time

𝑄St,disch 𝑘 Nominal discharge power of the storage tank

𝑟1(𝑡), 𝑟2(𝑡) Random coefficients between 0 and 1

𝑠CHP Boolean signal for the minimum operation time constraint

𝑠St Boolean signal for the temperature constraint

𝑇amb ℃ Ambient temperature

𝑇high ℃ High temperature level

𝑇low ℃ Low temperature level

𝑇s ℃ Temperature of the storage tank

𝑐min ℎ𝑜𝑜𝑟𝑠 Minimum operation time

𝑑 𝑘 (𝑚_⁄ 2_{∙ 𝐾)} _{Overall heat transfer coefficient}

𝑜b 𝑘 Heat supply by the boiler

𝑜CHP 𝑘 Heat supply by the CHP

𝑜CHP,St 𝑘 Heat supplied to the storage tank

𝑜load 𝑘 Heat load

𝑜st 𝑘 Heat supply by the storage tank

𝑉 𝑚3 _{Volume of the storage tank}

𝑣ik(𝑡) 𝑘th velocity vector component of the particle 𝑖

𝑘 𝐽 Work done on a control mass

𝑘el 𝑘𝑘ℎ Total electricity production for an evaluation period

𝑥ik(𝑡) 𝑘th position vector component of the particle 𝑖

(17)

xvii

∆𝑇 ℃ Temperature difference

𝜂CHP,th Thermal efficiency of the CHP

𝜂CHP,el Electrical efficiency of the CHP

(18)

(19)

1

1. INTRODUCTION

Supervisory control theory has been initiated in 1989, which deals with comprehensive ways of solving control problems with discrete event systems (DES) (Ramadge and Wonham 1989, p. 81). These systems can be encountered at logistic systems, traffic systems, manufacturing workcells and many other industrial processes (Thistle 1996, p. 25). Possible event sequence of such systems is illustrated in Fig. 1.1.

Fig. 1.1 State trajectory of a DES (Ramadge and Wonham 1989, p. 82)

As shown in Fig. 1.1, such systems include state transitions at discrete time periods. These states indicate certain events to be performed for a certain amount of time. At many industrial processes including the DESs, optimizing the order of these events make enormous changes in terms of timing and efficiency of the process. According to Ramadge Wonham (RW) framework as represented in Fig 1.2, a supervisor controls the events generated by the plant and terminate them when needed.

Fig. 1.2 Ramade Wonham (RW) Framework (Morgenstern and Schneider n.d., p. 3)

In this thesis a supervisory control of a combined heat and power (CHP) plant is achieved by an economic optimization. Thus, an optimization tool acts as a supervisor and optimizes the best trajectory of possible discrete events rather than controlling them.

(20)

2

1.1. Motivation

Concurrent production of electrical and thermal energy by the CHPs has been accepted as a good way of increasing the overall efficiency of district heating systems (DHS). German parliament has approved a new CHP law in 2008, which aims to increase the total share of electricity production by the CHP units up to 25 % by 2020 (BMWi 2012, p. 2). Considering this information one can claim that the CHPs will have a significant role in electricity production in the near future. On the other hand, Germany is one of the European leaders in renewable energies and many supporting schemes and incentives for renewable energy based plants are in force in Germany. Therefore, it can be anticipated that wind and solar power production will also take a higher share in electricity production, thus playing a significant role in variation of electricity prices in the future. Thus, producing electricity at peak hours is desired for financial reasons (Streckiene et al. 2009, p. 2308). Therefore, an optimal operation of a DHS will be significant in future’s highly fluctuating electricity market in Germany. By an optimal operation of a DHS it is meant that the system responds to the changes in electricity prices and distributes the load among the heat production units within the system, so that the highest financial gain is obtained and the heat demand is fully covered. The DHS should be also able to cover the heat load while producing electricity. This is why combining a heat storage with the CHP units is a widely preferred solution for a flexible operation of a DHS, since the heat storage is able to store the excess heat produced by the CHP unit especially at peak hours, when producing electricity is more favorable than other time periods of a day (Zhao et al. 1998). Therefore, an economic optimization has been carried out in this thesis in order to achieve an optimal operation of a DHS.

1.2. Literature Survey

In order to analyze all subcomponents of the DHS, it is common to model such systems within a simulation environment for carrying out an optimization process. Simulation-based optimization is a promising way of analyzing the performance of a system with complex natures. There are many simulation tools, such as “EnergyPlus”, “TRNSYS” and “Dymola”, which can be coupled with an optimization tool to optimize various parameters in heating, ventilation and air conditioning (HVAC) systems of buildings and district heating systems in general. Parametric simulation method is an alternative technique to optimization, at which each single variable is manually parameterized to see the effect on the performance of the systems while all other variables are kept constant (Nguyen et al. 2013, p. 1044). Nevertheless without doing an optimization, achieving a good performance of such systems might be a cumbersome task considering the effort and the time spent. A simulation-based optimization can be summarized in following steps:

(21)

3

• mathematical modeling and formulation of the optimization problem, setting constraints and design variables

• determining a tool for implementation and simulation of the model

• determining an optimization tool and an algorithm to be coupled with the simulation tool and executing the optimization process

• controlling whether the global optima is found

Arguably, while planning a DHS, it must be analyzed what kind of heat production units should be utilized and how they should share the heat demand. Furthermore, temperature levels, at which heat should be transferred, must be well specified in order to operate the DHS in a cost-effective way. Until today, many parameters have been investigated and several approaches have been taken to minimize the operational costs of the DHS. Those can be categorized into three groups (Benonysson et al. 1995, p. 299):

• optimal distribution of heat load between different heat producing units and a heat storage unit, if applicable

• minimization of supply temperatures of each heat production unit, so that the cost savings through less fuel consumption can be achieved

• optimizing the DHS dynamics, such as start-stop times of pumps, the temperature levels at each heat production unit and as well as load distribution among the units concurrently

Modeling and expressing all physical behaviors of the components in a DHS brings also complexities, i.e. extra constraints to be handled, together. That’s why actual models should be replaced with simplified models for optimization purposes. A simplified model is a reduced order model, which only captures main behaviors of the actual model. It can be then used for special optimization purposes. It is asserted that optimizing physically expressed actual models require undesirably large execution times with respect to simplified models (Chandan et al. 2012, p. 3070). Furthermore, it is pointed out that simplified models perform similar results to those obtained by optimizing actual models (Nguyen et al. 2013, p. 1045). However, it is also necessary to note that simplified models lead to an uncertainty within the model. Therefore, the assumptions being made to simplify the actual models should be realistic to get accurate and reliable simplified models. To meet this end, in this master’s thesis a reduced order model, which represents the main features of a detailed DHS, is modeled and the aforementioned first optimization approach, which focuses on sharing the heat demand optimally among the units involved in the system, is taken to minimize the overall system costs.

As mentioned, it is aimed to share out the heat demand optimally among the heat production units while considering the fluctuations in electricity prices. This optimization

(22)

4

problem has been modeled as a discrete optimization problem, where the discrete variables define the modes of operation of the DHS. These modes or operating states are introduced in Section 3.1. With this in mind, in literature there are two well-known problems that are linked to the discrete optimization problems:

First one is the Travelling Salesman Problem (TSP) that is related to discrete optimization problems. There are many articles written and solution methods are proposed for the TSP. Consider a region of n cities. A salesman is supposed to stop by each of n cities, aiming to take the shortest route and return to the starting point and each random route is a possible solution of the TSP (Laporte 1992, pp. 231-232). Similarly the Knapsack Problem (KP) is the second typical example of discrete optimization problems. The optimization problem is finding the best combination of items to be packed in a knapsack, in which the total weight must be lower than or equal to a certain limit and usually each item has a monetary value disproportional to its weight (Suzuki 1978, p. 162). Both of these problems can be subject to some other constraints.

A literature study has shown that aforementioned problems can be generally called as Integer Programming (IP) problems or 0-1 IP problems more specifically. An unconstrained linear 0-1 IP problem can be mathematically represented as

𝑚𝑚𝑥 � 𝑐j∙ 𝑥j , 𝑛 𝑗=1 (𝟏. 𝟏𝟏) subject to � 𝑠j∙ 𝑥j 𝑛 𝑗=1 ≤ d , (𝟏. 𝟏𝟏) 𝑥j 𝜖 {0,1}𝑛 . (𝟏. 𝟏𝟏)

Considering the KP, 𝑐_j would be the monetary value of each item; 𝑠_j could be interpreted as the weight of the item; and 𝑑 would be the maximum allowed weight in the knapsack; 𝑥_j would stand for whether the item is selected or not and 𝑛 would be the number of items. There are many different formulations IP problems depending on the constraints of the problem. When 𝑥_j can take both integer and non-integer values, then this particular problem is called mixed integer programing (MIP).

There are many methods to solve IP and MIP problems. These can be grouped under two main headlines: stochastic and deterministic ones. Stochastic methods provide a solution of coping with generally nonlinear, high dimensional complex systems, which are not suitable for other classical deterministic optimization methods (Spall 2004, p. 170). Many of them are nature-inspired, and they make use of iterative trial and error processes to converge an optimal solution. Thus, stochastic methods do not assure finding the global

(23)

5

optima. Deterministic methods, on the contrary, make use of the analytical properties of the problem to converge to the global solution (Lin et al. 2012). There are no stochastic elements and assumptions within the algorithms, they rely on linear algebra and are mainly based on computation of gradient of the variables for convergence (Cavazzuti 2013, p. 77). One advantage of deterministic methods is that they guarantee finding the global optima eventually. However, in comparison to stochastic methods, they are more compute-intensive and therefore require higher execution times (Collet and Rennard 2006).

Some examples of both these methods are explained below.

1.2.1. Deterministic Methods

One of the well-known algorithms to solve the IP problems is Branch and Bound (BB) method. It is an efficient algorithm for systematically enumerating of discrete problems. The algorithm mainly includes two steps. First step is branching, which is a relaxation of the feasible region by dividing it into several sub-regions. Over these sub-regions upper and lower bounds are determined, which is known as the bounding step of the algorithm. The idea is eliminating the branches by comparing upper and lower bound values of each sub-region, thus reducing the search space.

Cutting Planes method aims simply to tighten the feasible solution region by adding linear

inequality constraints, which cut the search plane into two groups and omit some of the non-integer values. By doing this recursively, an non-integer solution is found. Although some problems are solved using cutting planes method, it is asserted that it behaves differently on similar problems and shows weak performance in reducing the size of search space. Nevertheless when it’s combined with BB, the method is known as “branch and cut”, it demonstrates a good performance on solving IP problems (Sherali and Driscoll 2000).

Benders Decomposition is another solution method for solving especially integer

programming problems. It is a class of multistage optimization algorithms. This method, unlike the other traditional approaches, divides variables in the optimization problem into two sets: master problem variables and sub-problem variables. First a first-stage master problem with an arbitrary number of variables is solved. Then some of the set of solutions are eliminated by the information interpreted from the results of sub-problem. This elimination procedure, which is also known as Bender’s cut, is done iteratively, thus leading to the optimality (Taskin 2010).

1.2.2. Stochastic Methods

Simulated Annealing (SA) method is a local search meta-heuristic method, which is used

(24)

6

slow cooling, reaches its best lattice energy state by avoiding the crystal defects. The name itself stems from annealing in metallurgy. The SA method inspires from this thermodynamic behavior and establishes a connection in the algorithm (Henderson et al. 2003, p. 288).

Genetic Algorithm (GA) finds solution to optimization problems using principles inspired

by genetics and natural selection, such as crossover, selection, and mutation. Each possible solution or chromosome is assigned to a fitness value. Then, relatively “fit” chromosomes are selected for reproduction of new individuals. A crossover operator is used for exchanging the genetic information between each chromosome, and mutation operator is used to keep enough diversity among the chromosomes to avoid premature convergence (Deep et al. 2009, p. 506). In literature many applications of the GAs are available for both integer and mixed integer programming problems. A similar algorithm to the GA is Differential Evolution (DE) algorithm. The main difference in finding better solution is that, the GA makes use of selection operator, whereas the DE relies on mutation operator (Karaboga and Ökdem 2004, pp. 53-54).

Among all of these methods, it is not intended to show pros or cons of each algorithm, but to give a general overview on solution methods for IP problems to the reader. There is no clear consensus that any of the algorithms mentioned above outperforms the other one though. However, a greater emphasis will be placed on explaining Particle Swarm

Optimization (PSO) algorithm, which is a type of stochastic algorithms, in Section 5.2.

Because it is the algorithm being used for solving discrete optimization problem described within this thesis.

1.3. Thesis Objectives

Main objective of this thesis is to carry out an economic optimization of a district heating system. It concerns following tasks:

• mathematical modeling of the economic optimization problem

• implementing constraints and improving the already existing DHS model, which included the main units of the DHS

• constructing the cost function of constrained optimization problem (The coupling of the model with the optimization tool had already been carried out.)

• plausibility analysis of the model via several test cases

• comparing the power-driven operation with the heat-driven operation with respect to overall financial gains and grid interactivity

(25)

7

2. SYSTEM

2.1. Structure of a District Heating System

DHSs are highly utilized systems for space heating and water heating in industrial and residential buildings. The DHS consists of heat production units, which could be either heat-only plants or CHP plants coupled with a boiler and a storage tank. The CHPs are power generation plants, which make use of the thermal energy in exhaust gases, cooling systems or other thermal energy waste streams to enhance the overall efficiency of the system. CHPs can be powered by fossil fuels, such as gas, oil or coal. Wood or biogas in the case of biomass boilers or heat from earth of in the case of geothermal energy can be used to heat up the water in the boiler. The water is then channeled to the various places of consumption, which are industrial, public, or private buildings. A district heating system can be designed to serve only a few residences, as well as a small village or a part of a city. Using CHP plants within the DHS is a prevalent way of supplying local energy demand for public authorities or industries. A schematic representation of a DHS is given in Fig. 2.1.

Fig. 2.1 Schematic representation of a DHS

As seen in Fig. 2.1, the CHP generates two types of energy from a single authority: electricity and heat for local use. The principle is as follows: the boiler produces a superheated steam. This steam then powers a turbine, which is connected to a generator to produce electricity. Similarly a CHP can be powered by a combustion engine to generate energy. Electrical energy produced is sent to the power grid to supply local users. The heat given off by the turbine, which is usually lost, is recovered to provide heat for a local

(26)

8

community or an industrial facility. Thus, the CHP offers a higher level of efficiency than that of solely heat or power producing plants.

2.2. District “Gutleutmatten”

In this thesis, “Gutleutmatten” (GLM), which is a part of the city Freiburg, is the district in question. The project will be completed by 2016 with the joint partnership of Fraunhofer Institute for Solar Energy Systems (ISE), Badenova Wärmeplus GmbH, and Solites. To meet residential heating demand, decentralized solar thermal systems and a district heating system will be utilized. On a total of 40,000 𝑚2 heated floor area, 32 buildings including around 500 households will be built (Oliva et al. 2014, para. 1). All buildings will have a decentralized storage unit, with a total capacity of 200 𝑚3 (Elci et al. 2014, para. 7). Based on the energy standard “KfW-Effizienzhaus 55” as cited in (Elci et al. 2014, para. 6), specific heating demand of 35 kWh/(m2a) is assumed for each building. Main power supply units will be a biogas-fired CHP with an electrical nominal power of 600 𝑘𝑘_el and thermal capacity of 654 𝑘𝑘_th and a natural gas-fired boilers with total capacity of 4000 𝑘𝑘_𝑠 (Elci et al. 2014, para. 6). The district heating station for District GLM already in operation and it is located close to the district, and it is also capable of supplying heat to other buildings outside the district GLM. A site plan of the district GLM is shown in Fig. 2.2.

Fig. 2.2 Site plan of District “Gutleutmatten”(Elci et al. 2014)

2.3. System Boundaries

Before describing the model, system boundaries are clearly defined and they stay consistent for any of the test cases, which will be discussed in the next sections. The system consists of a boiler, a CHP and a storage tank unit. First and foremost fuel burned at the boiler and at the CHP is regarded as consumptions, while the electricity fed into the grid and the supplied heat to consumer are regarded as gains in financial terms.

(27)

9

It is also necessary to address following assumptions that define the system boundaries: • The electricity produced by the CHP unit is sold at the energy exchange. However,

the electricity consumption of pumps and other auxiliary units in the DHS is neglected.

• Gas, electricity and heat prices as well as heat demand and technical specifications of the CHP and the boiler are assumed to be known. They are the disturbances of the system, which are neither controlled nor predicted.

• Investment and maintenance costs are ignored and they are not playing any role for the economic optimization. Both the boiler and the CHP are expected to operate at any periods of a year. Decommission of the heat production units due to maintenance is ignored.

• There is a continuous fuel supply available to run the CHP and the boiler. • The DHS must be able to cover the heat demand under any circumstances.

(28)

(29)

11

3. MODELING

3.1. Operating States

As mentioned before, discrete variables of the optimization problem are the modes of operation or the operating states of the DHS.

Maintenance services for the CHPs are regularly performed depending on the yearly operation time of the units. Maintenance costs per unit operation time is generally fixed, thus the CHPs are generally operated at nominal power in order to minimize the proportion of the maintenance costs with respect to the gain of the CHP per unit operation time. That’s why; it is assumed that there are two operating possibilities for the CHP, either running at nominal power or “shutdown state”. On the other hand, the boiler and the storage tank are flexible to supply heat at any level. Considering there are three main units, one can argue that there are eight possible combinations that these three units can share the heat demand. Table 3.1 shows all possible states, at which a DHS can operate. “ON” stands for that the unit supplies heat or electricity, whereas “OFF” stands for that the unit is shut down.

Table 3.1 Operating States

States CHP Boiler Storage Tank

1 ON OFF OFF 2 ON ON OFF 3 OFF OFF ON 4 ON OFF ON 5 ON ON ON 6 OFF ON OFF 7 OFF ON ON

8 OFF OFF OFF

3.2. Storage Tank Model

Storage tanks, which are operated properly, are generally stratified with respect to the temperature levels of the water in the tank. The water, which is fed into the tank, is placed to the corresponding part of the storage tank depending on its temperature. Generally, the upper part of the tank stores relatively hotter water, which is then supplied to provide the needed heat for domestic use. Since only input-output relations are considered and all hydraulic parts of the DHS are neglected for the simplicity, the storage tank is simplified and represented as a single mass rather than a stratified model. The storage tank is basically charged by the CHP. This charging power depends on the heat load. As long as the heat

(30)

12

load is less than the CHP nominal power, then excess heat produced by the CHP is stored in the storage tank. Fig. 3.1 shows a single node representation of the storage tank. Additional to the charging and the discharging affect, a heat loss through its outer wall is considered. It’s assumed that the temperature of the storage tank, ambient temperature, and properties of the water are known to calculate the temperature in the next time step.

Fig. 3.1 Single node representation of the storage tank

Following formulations show how the charging and the discharging processes affect the temperature of the storage tank. Amount of net heat transfer 𝑄_net into the body can be calculated as

𝑄net= 𝑄ch− (𝑄disch+ 𝑄loss) , (𝟑. 𝟏)

where 𝑄_ch and 𝑄_disch stand for the charging and the discharging powers respectively and 𝑄loss stands for the amount of the heat lost to the environment.

Assume that the temperature of the single mass represented in Fig. 3.1 at time 𝑐 is denoted by 𝑇(𝑡). When this mass is subject to a heat transfer that amounts to 𝑄_net until time (𝑐 + 1), then the temperature at that time 𝑇(𝑡+1)_{can be calculated by using 1}st_{law of}

thermodynamics, which is

𝑄 − 𝑘 = 𝑑𝑑 + 𝑑𝑑k+ 𝑑𝑑e , (𝟑. 𝟐)

where 𝑄 is the heat given to the system; 𝑘 is the work done on the system; 𝑑𝑑 is the change of internal energy; 𝑑𝑑_k is the change of kinetic energy and 𝑑𝑑_e is the change of potential energy.

The system boundary in this case is the mass including all heat transfers. In thermodynamics, everything within a closed system is called a control mass. In this case, there is no work done on the system, and 𝑑𝑑_k and 𝑑𝑑_e are generally negligible for a stationary control mass. Internal energy change of the control mass at constant pressure can be represented as

(31)

13

where 𝜌 is the density [𝑘𝑘 𝑚⁄ 3] of the control mass; 𝑉 is the volume [𝑚3]; of the control mass 𝑐_p is the specific heat capacity [𝐽 (𝑘𝑘 ∙ 𝐾)⁄ ] of the control mass and ∆𝑇 is the temperature difference [𝐾] due to heat transfer.

It is now possible to apply the 1st law of thermodynamics to the control mass. The net heat transfer would be then

𝑄net= 𝜌 ∙ 𝑉 ∙ 𝐿p∙ �𝑇(𝑡+1)− 𝑇(𝑡)� . (𝟑. 𝟒)

When the unknown term is get alone on the left side, the final form of the equation will be 𝑇(𝑡+1)₌ 𝑄net

𝜌 ∙ 𝑉 ∙ 𝑐p+ 𝑇

(𝑡)_{. (𝟑. 𝟓)}

Let us now derive the term 𝑄̇_loss. According to the Newton’s law of cooling it can be written as

𝑄̇ = 𝑑 ∙ 𝐴 ∙ ∆𝑇 , (𝟑. 𝟔) where 𝑄̇ is the heat removed or gained [𝑘] by convection per unit time; 𝑑 is the overall heat transfer coefficient [𝑘 (𝑚⁄ 2∙ 𝐾)]; 𝐴 is the heat transfer surface area [𝑚2] and ∆𝑇 is the temperature difference [𝐾] between ambient and the control mass.

When the corresponding temperatures are placed in Eq. (3.6), following formulation can be derived as

𝑄̇loss= 𝑑 ∙ 𝐴 ∙ �𝑇(𝑡)− 𝑇amb� , (𝟑. 𝟕)

where 𝑇_amb is the ambient temperature.

3.3. Mathematical Model of the System

For simplicity, Table 3.2, Table 3.3, and Table 3.4 are given to describe the parameters, the outputs, and the disturbances of the system, which will then be used in the equations describing the mathematical model. Additional to the parameters (𝜌, 𝑑, 𝐴, 𝑉 and 𝑇_amb) which are explained in the previous section, other parameters are also listed in Table 3.2.

Table 3.2 Parameters of the system

Parameters Description

𝑸𝐂𝐂𝐂,𝐭𝐭 Nominal thermal power of CHP

𝑸𝐂𝐂𝐂,𝐞𝐞 Nominal electrical power of CHP

(32)

14

Table 3.3 Disturbances of the system

Disturbances Description

𝒖𝐞𝐥𝟏𝐝 Heat Load

𝒄𝐞𝐞 Electricity price

𝒄𝐟𝐟𝐞𝐞 Fuel price

𝒄𝐭𝐞𝟏𝐭 Heat Price

Table 3.4 Outputs of the system

Outputs Description

𝒖𝐂𝐂𝐂 Heat supply by CHP (only to consumer)

𝒖𝐂𝐂𝐂,𝐒𝐭 Heat supply to storage tank

𝒖𝐝𝐭 Heat supply by storage tank

𝒖𝟏 Heat supply by boiler

𝑻𝐝 Temperature of storage tank

𝑷𝐞𝐞 Electricity production

The outputs of the system depend on the operating states of the system. In order to formulate the whole system mathematically, the outputs of the system are separately defined for each operating state and unit.

Let us first consider the CHP. As known, it can charge the storage tank, supply heat to the consumer, do both actions simultaneously, or it can be shut down.

i. When CHP charges storage tank and supplies to the consumer ;

𝑜CHP = 𝑚𝑖𝑛(𝑄CHP,th , 𝑜load) 1 , (𝟑. 𝟖)

𝑜CHP,St= 𝑚𝑚𝑥((𝑄CHP,th− 𝑜load) , 0) 2 , (𝟑. 𝟗)

𝑑el = 𝑄CHP,el . (𝟑. 𝟏𝟏) ii. When CHP is shut down;

𝑜CHP = 0 , (𝟑. 𝟏𝟏)

𝑜CHP,St= 0 , (𝟑. 𝟏𝟐)

𝑑el = 0 . (𝟑. 𝟏𝟑) iii. When CHP supplies heat only to the consume;

1

min(a,b) outputs the minimum of a and b.

2

(33)

15

𝑜CHP = 𝑄CHP,th (𝟑. 𝟏𝟒)

and Eq. (3.12) and Eq. (3.13) also hold true in this case.

The boiler can either supply heat or not. The supplied heat is the remaining heat demand after the CHP and the storage tank supplied to the consumer. The two equations for those actions are

𝑜b= 𝑜load− (𝑜CHP+ 𝑜st) , (𝟑. 𝟏𝟓)

𝑜b= 0 . (𝟑. 𝟏𝟔)

In order to derive the equations for the storage tank, three cases are separately considered. Temperatures are calculated based on the formulations derived in Section 3.2.

i. Storage tank is charged, but no discharge is allowed at the same time.

The net heat transfer is then the difference between supplied heat by the CHP and heat loss to the environment. The corresponding equations are

𝑜st= 0 , (𝟑. 𝟏𝟕)

𝑄̇net= 𝑜CHP,St− 𝑄̇loss = 𝑜CHP,St− 𝑑 ∙ 𝐴 ∙ �𝑇s(𝑡)− 𝑇amb�. (𝟑. 𝟏𝟖)

When 𝑄_net in Eq. (3.18) is substituted into Eq. (3.5), then following equation is obtained:

𝑇_s(𝑡+1)=�−𝑑 ∙ 𝐴 + 𝜌 ∙ 𝑉 ∙ 𝑐p� ∙ 𝑇s

(𝑡)_{+ 𝑜}

CHP,St+ 𝑑 ∙ 𝐴 ∙ 𝑇amb

𝜌 ∙ 𝑉 ∙ 𝑐p . (𝟑. 𝟏𝟗) ii. Storage tank is neither charged nor discharged.

Eq. (3.6) holds true in this case as well. The net heat transfer is basically only the heat loss, which is

𝑄̇net= −𝑄̇loss= −𝑑 ∙ 𝐴 ∙ �𝑇s(𝑡)− 𝑇amb�. (𝟑. 𝟐𝟏)

When 𝑄_𝑛𝑛𝑡 in Eq. (3.20) is substituted into Eq. (3.5), then following equation is obtained:

𝑇_s(𝑡+1) =�−𝑑 ∙ 𝐴 + 𝜌 ∙ 𝑉 ∙ 𝑐p� ∙ 𝑇s

(𝑡)_{+ 𝑑 ∙ 𝐴 ∙ 𝑇} amb

𝜌 ∙ 𝑉 ∙ 𝑐p . (𝟑. 𝟐𝟏) iii. Storage tank is not charged, but discharged.

In this case the net heat transfer is sum of the discharge power and the heat loss, which can be written as

(34)

16

The discharge power of the storage tank is the remaining heat demand after the CHP supplied to the consumer provided that it is less than the maximum discharge power of the storage tank. Heat supply by the CHP (𝑜_CHP) can also be “0”. The discharge power can be written as follows:

𝑜st= min�(𝑜load− 𝑜CHP) , 𝑄St,disch�. (𝟑. 𝟐𝟑)

When 𝑄_net in Eq. (3.23) is substituted into Eq. (3.5), then following equation is obtained:

𝑇s(𝑡+1) =�−𝑑 ∙ 𝐴 + 𝜌 ∙ 𝑉 ∙ 𝑐p� ∙ 𝑇s (𝑡)_{− 𝑜}

𝑆𝑡+ 𝑑 ∙ 𝐴 ∙ 𝑇amb

𝜌 ∙ 𝑉 ∙ 𝑐p . (𝟑. 𝟐𝟒)

As can be anticipated, these equations are state dependent. Therefore, depending on the operating state some of the equations might not be valid. It is shown in Table 3.5, when these equations hold true. Each column has 6 green boxes, indicating that the corresponding equation on the left column is used for calculations. It is worth noting that these 6 equations, which are active at each state, are used to calculate the system outputs given in Table 3.4.

Table 3.5 Active equations depending on the operating state

States 1 2 3 4 5 6 7 8 Eq. (3.8) Eq. (3.9) Eq. (3.10) Eq. (3.11) Eq. (3.12) Eq. (3.13) Eq. (3.14) Eq. (3.15) Eq. (3.16) Eq. (3.17) Eq. (3.19) Eq. (3.21) Eq. (3.22) Eq. (3.24)

(35)

17

3.4. Formulation of the Cost Function

Before formulating the cost function, it is useful to show some fundamental mathematical expressions of individual units. Fig 3.2 shows all the inputs and outputs (power-based) of each single unit. 𝑄_fuel,CHP and 𝑄_fuel,Boiler are fuel consumptions by the CHP and the boiler

respectively. 𝑄_loss,CHP, 𝑄_loss,Boiler and 𝑄_loss,St are heat losses of individual units as shown. Other variables are already described in Table 3.3.

Fig. 3.2 Energy flow within the system boundary

General energy balance equations of each single unit can be listed as

𝑄CHP,th= 𝑜CHP+ 𝑜CHP,St , (𝟑. 𝟐𝟓)

𝑄fuel,CHP= 𝑄CHP,th+ 𝑄loss,CHP , (𝟑. 𝟐𝟔)

𝑄fuel,Boiler= 𝑜b+ 𝑄loss,Boiler , (𝟑. 𝟐𝟕)

𝑜CHP,St= 𝑜st+ 𝑄loss,St . (𝟑. 𝟐𝟖)

As already shown in Section 3.2, heat loss within the storage tank has been mathematically formulated. However for the boiler and the storage tank efficiencies are taken into account to consider losses. The amount of energy that needed to be extracted can be found by the equations

𝑄fuel,CHP= 𝑄_𝜂CHP,th CHP,th =

𝑑el

(36)

18 and

𝑄fuel,Boiler=𝑜_𝜂b

b , (𝟑. 𝟑𝟏)

where 𝜂_CHP,th is the thermal efficiency of the CHP; 𝜂_CHP,el is the electrical efficiency of the CHP and 𝜂_b is the efficiency of the boiler.

Based on these expressions, the gain of each unit can be formulated as follows:

𝐽CHP= � [𝑑el∙ 𝑐el+ 𝑜CHP∙ 𝑐heat] − 𝑄fuel,CHP∙ 𝑐fuel� ∗ 𝜏 , (𝟑. 𝟑𝟏)

𝐽Boiler= �𝑜b∙ 𝑐heat− 𝑄fuel,Boiler∙ 𝑐fuel� ∗ 𝜏 , (𝟑. 𝟑𝟐)

𝐽St= [𝑜st∙ 𝑐heat] ∙ 𝜏 , (𝟑. 𝟑𝟑)

where 𝜏 is the operation time in hours [h]; 𝑐_el , 𝑐_heat and 𝑐_fuel are the previously described prices in cents per kilowatt-hour [ct./kWh]; 𝑑_el , 𝑜_CHP, 𝑄_fuel,CHP, 𝑜_b , 𝑄_fuel,Boiler and 𝑜_st are in kilowatts [kW]; 𝐽_CHP , 𝐽_Boiler and 𝐽_St are the gains of the CHP, the boiler and the storage tank in cents [ct.].

As can be understood, the gains are basically calculated on the basis of extracting consumption by production. It is necessary to note that heat supply to the consumer (𝑜_CHP,St) is not counted as a gain of the CHP, since it is not directly sold to the consumer. That’s why; the actual heat sold to the consumer 𝑜_st is counted as the gain of the storage tank.

There are mainly two constraints within the model that should be considered before deriving the cost function.

• First one is the temperature constraint for the storage tank. The temperature of the storage tank must not be charged any more when it is full.

• Secondly, a minimum operation time for the CHP should be taken into consideration. This is to make sure that the CHP is warmed up and, is running at optimal temperature for a certain minimum time before it is shut down again. Similarly when CHP is shut down, it should wait at least for minimum operation time before operating again. Moreover, with minimum operation time is desired since the life cycle of the CHP can be increased with a proper operation schedule.

Then the final form of the cost function is written as follows:

𝐽T= 𝑚𝑚𝑥 �( 𝐽CHP(𝑡) ∙ 𝑏CHP(𝑡) + 𝐽St(𝑡)∙ 𝑏St(𝑡) 𝐻𝑝

𝑡=1

+ 𝐽_Boiler(𝑡) ∙ 𝑏_Boiler(𝑡) ), (𝟑. 𝟑𝟒𝟏)

which is subject to the temperature constraint

(37)

19

and the minimum operation time constraint, which is only active

𝑤ℎ𝑒𝑛 𝑏_CHP(𝑡) − 𝑏_CHP(𝑡−1)≠ 0, 𝑐ℎ𝑒𝑛 � �𝑏CHP(𝑘) − 𝑏CHP(𝑘−1)� = (𝑡+𝑡min )

𝑘=𝑡

0 , (𝟑. 𝟑𝟒𝟏)

where 𝐽_T is the overall gain over a prediction horizon; 𝑇_high is the allowed maximum temperature level for storage tank; 𝑇_s is the temperature of the storage tank; 𝑐_min is the minimum operation time period, at which the CHP should maintain its mode of operation; 𝑏CHP, 𝑏St, 𝑏Boiler∈ {0,1}𝑝 are binary variables, 𝑝 is the number of time steps in the prediction

horizon and 𝐻_p is the prediction horizon.

As can be deduced, the optimization problem is maximizing the overall system gains over a prediction horizon. It is worth mentioning that the combination of binary variables3𝑏_CHP, 𝑏_St and 𝑏_Boiler are indicating whether the unit is operating, thus forming operating states of the system. Consider i.e. {𝑏_CHP, 𝑏_St, 𝑏_Boiler} = {1,0,0}, which would be the “Operating State 1” according to the Table 3.1. Therefore, optimization problem can be also regarded as finding the best sequence of operating states over a prediction horizon so that overall system gains are maximized.

Every time when these binary variables change, then the operating state and the related equations change as well. Although the gains of each single unit are second-based, the binary variables are not expected to change at any second. Thus, a number of time steps should be defined to fix the maximum number of input changes over a prediction horizon.

As can be understood from Eq. (3.34c), the minimum operation time constraint is not always active. The constraint is only active when 𝑏_CHP ∈ {0,1} changes its value. That would mean that the CHP switches its operation mode from “ON” to “OFF” or vice versa. When it is the case, then the constraint is active and Eq. (3.34c) makes sure that, the CHP maintains its operation for a certain minimum time period. After that, the constraint is inactive again until 𝑏CHP variable takes another value. There might be cases, where one constraint must be

violated in order not to violate the other constraint. In those cases, the temperature constraint will be regarded as a soft constraint, which can be violated in some rare occasions.

3

(38)

20

3.5. Model Identification

Figure 3.3 shows the input-output relation of the model.

Fig. 3.3 Input-outputs and disturbances of the system

The operating state, which is the function of binary variables, is the input variable. It is assumed that the operating states of the system are discrete variables and they are assumed to be constant for a certain time step. Likewise, electricity price and heat load information can only be obtained in discrete-time by authorities. That’s why all the outputs of the system are in discrete-time as well.

It is also necessary to mention that it is a time-dependent linear model. The gains of the system 𝐽_CHP, 𝐽_St and 𝐽_Boiler as derived previously are in the following form:

𝐽(𝑡)_{= 𝑚}

1(𝑡)∙ 𝑜1(𝑡)+ 𝑚2(𝑡)∙ 𝑜2(𝑡)+ ⋯ , (𝟑. 𝟑𝟓)

where 𝑚₁ and 𝑚₂ can are time-variant constants (prices); and 𝑜₁ and 𝑜₂ are the outputs (power production and heat consumption or production) which form the gains. As can be deduced, the gains of the system are linear combinations of these variables.

(39)

21

4. IMPLEMENTATION

For the implementation of the mathematical model, Modelica based tool Dymola (Dynamic Modeling Laboratory) has been utilized. Modelica is a widely used programming language for modeling of large and complex systems, for instance, process oriented systems, air conditioning systems, and many other applications involving electrical, hydraulic, thermal and control subsystems (Otter 2009, p. 7). In Modelica, systems can be mathematically expressed by algebraic, differential or discrete equations. For graphical editing and visualizing simulation results, commercial Modelica environment Dymola is used.

Dymola has two main windows: modeling and simulation. One can very quickly make

changes in the model, and see the results by plotting or animating them at simulation window.

As explained before the model is a discrete-time model with discrete input variables. Implementation of the model is achieved by defining each operation state and related mathematical formulations. A screenshot from the main window of the implemented model is given in Fig. 4.1.

Fig. 4.1 Implementation of the system at “Dymola”

As seen, three main units, the CHP, the boiler, and the storage tank are implemented. The operating states of the system are stored in the table “Operating_States”. At each time step, information about the operating state is conveyed to the main units as integer variables.

(40)

22

All parameters, variables and formulations are written as code within these units. “Heat_Load” block enables all the units to communicate each other. The information of how much power each single unit supplies is gathered at this block and then distributed to the single units. By this arrangement, an energy balance is set up between the heat demand and the total heat supply. Electricity and the heat demand are given to the model as text (txt.) files. Electricity prices are stored in the table “Electricity_Price”, whereas the heat demand is stored in “Heat_Load” block. “Overal_Gain” operator within the blue box on the top sums up the overall system gains continuously.

Constraints are implemented within the model since the optimization tool, which is going to be explained in Section 5, does not support constraint handling within the tool. Thus, constraints are handled with a penalty function approach within the model. For each constraint, a corresponding Boolean signal is defined. By this way, one can see in the simulation window whether a constraint is violated based on the signal value becoming “True” or “False”. In the next sections, how the temperature and minimum operation time constraint has been implemented is explained.

4.1. Implementing minimum operation time constraint

Two timers are used to implement the constraint. Fig. 4.2 illustrates how the minimum operation time is implemented.

Fig. 4.2 Implementation of the minimum operation time constraint

4.2. Implementing storage tank model

As seen in Fig. 4.2, the operating state is given as input to “Integer_to_Boolean” block. At this block integer input is converted to a Boolean variable, where “True” stands for that the CHP is running and “False” stands for that the CHP is shut down. Then these signals are conveyed to the timers. The timer “ON” starts counting, when the CHP operates. On the contrary, the timer “OFF” starts counting, when the CHP is shut down. Then the information

(41)

23

from the timers is conveyed to the next block to check whether the minimum operation time is reached. Boolean signals “On_Timer” and “Off_Timer” show initially “False” as default. When the threshold value, which is the minimum operation time, is reached, then the corresponding Boolean output signal switches to “True”. By this, one can easily figure out whether minimum operation time is reached. For instance, if the Boolean signal “Off_Timer” shows “False” while the CHP is running, then that would mean the minimum operation time is not reached, and the constraint is violated. A proper operation would be running the CHP after “Off_Timer” switches to “True”.

Fig. 4.3 Representation of the storage tank with a single node

It is necessary to define high and low temperature levels to indicate whether there is heat available in the storage tank. The output and the input of the storage tank model are power-based signals. The temperature of the storage should be interpreted as follows:

It is assumed that the discharge temperature of the storage tank is constant. However, the temperature which is defined at this heat capacitor is just an indicator, showing the availability of the storage tank. For instance, when the temperature of the heat capacitor is equal to a low limit, then it is assumed that the storage tank is empty. On the other hand, when the temperature of the heat capacitor reaches to a high limit, then this would mean that the storage tank is full and should not be charged any more.

Therefore, as seen in Fig. 4.4, a low temperature level 𝑇_low is defined to cut off the discharge power. Likewise, a high temperature level 𝑇_high is defined to avoid charging the storage tank.

(42)

24

Fig. 4.4 Operation range of the storage tank

As seen in Fig. 4.4, the storage tank is assumed to be empty, when its temperature is less than low temperature level. Whenever the temperature drops below this level, the storage tank does not supply heat. However that’s not regarded as a constraint but a usual occasion that can happen in real plants, as the ambient temperature mainly determines the temperature of the storage tank in winter. The interval between the high and the low temperature levels is the normal operation region for the storage tank. When the temperature exceeds the high temperature level, then the Boolean signal shows “True” as output for a certain time. Fig. 4.5 illustrates how it is defined.

Fig. 4.5 Hysteresis band for the storage tank

When the temperature of the heat capacitor exceeds the level "T_high", the Boolean output becomes “True”. As soon as the temperature drops to the level "T_high− ∆𝑇", then the Boolean operator shows output as “False” again. The reason defining such a hysteresis band is to make sure that the temperature of the storage tank drops to a certain level (𝑇_high− ∆𝑇).

(43)

25

This ∆𝑇 value is lower than 0.1℃ and it is implemented due to the following reason. When temperature of the storage tank reaches to the level 𝑇_high, then the temperature constraint is violated and the system should automatically drop the temperature of the storage tank based on the control input given by the optimizer. However; it is not the case in simulation environment. The solver at Dymola is not able to capture the changes happening in a very short time span (milliseconds level). That’s why the storage tank is repetitively charged and discharged (CHP operates intermittently), when the temperature of the storage tank reaches to the level 𝑇_high. That’s why this very small hysteresis band is implemented to make sure that the CHP does not change its mode of operation so quickly that simulator do not recognize. Therefore, by this arrangement, the CHP is decommissioned for a short time (until the storage tank temperature drops by ∆𝑇) then the simulation programs already captures that the CHP is shut down and the minimum operation time constraint makes sure that the CHP does not charge right after the Boolean signal shows “False” again.

Computation of the cost function is done within the model. To do that, the constraints are combined with the cost function in such a way:

𝐽 = 𝐽_CHP(𝑡) ∙ 𝑏_CHP(𝑡) + 𝐽_Storage(𝑡) ∙ 𝑏_St(𝑡)+ 𝐽_Boiler(𝑡) ∙ 𝑏_Boiler(𝑡) − 𝑑1∙ 𝑠CHP(𝑡) − 𝑑2∙ 𝑠St(𝑡) , (𝟒. 𝟏)

where 𝑐 is the subscript denoting the time instant; 𝑑₁ and 𝑑₂ are penalty functions for the CHP and the storage tank respectively; 𝑠_CHPis the Boolean signal for minimum operation time constraint; 𝑏_CHP, 𝑏_St, and 𝑏_Boiler are signals indicating whether the CHP unit operates and 𝑠_St is the the Boolean signal for the temperature constraint.

It can be seen in Eq. (4.1) that the first three terms are already defined instantaneous gains from each individual unit. As can be understood, the constraints 𝑑₁ and 𝑑₂ are not always active. It depends on the previous actions (states) of the system whether a constraint is going to be active or not. Therefore when a Boolean signal becomes “True”, then the corresponding constraint in the Eq. (4.1) is active. Without doubt, the penalties 𝑑₁ and 𝑑₂ decrease the total gains of the system. This is why violating a constraint would mean that the system gain drops; therefore a violation of one of the constraints leads to an unfeasible result. How to select proper values for 𝑑₁ and 𝑑₂ is discussed in Section 6.5.

(44)

(45)

27

5. OPTIMIZATION

5.1. Optimization Tool

There are many projects at Fraunhofer ISE, which are conducted by different subdivisions. All ongoing works should be coordinated and organized in a way that everyone can get involved in them and be able to contribute to them when needed. It is common that most of the student works are part of a large-scale project. Thus, at Fraunhofer ISE each member of the project groups is supposed to work more or less with the same simulation, optimization and programming tools. As an optimization tool GenOpt, generic optimization program, is selected, as it can be easily coupled with Dymola and it is generally adopted tool at Fraunhofer ISE.

GenOpt is a tool that is used for minimization or maximization of cost functions calculated by an external simulation program. Being a freely downloadable java program, GenOpt is mainly an interface between any text-based simulation program and an optimization algorithm (Coffey et al. 2010, p. 1086). Any simulation program can be coupled with GenOpt, provided that the simulation program is able to read text (txt.) files and calculated cost function must be also written in a text (txt.) file by the simulation program (Wetter 2009, p. 9). GenOpt is not suitable for optimization problems, in which cost function is quadratic, or has a gradient (Wetter 2015). The standard library of GenOpt offers many optimization algorithms, and it’s also possible to extend library by the user with new optimization algorithms (Wetter 2009, p. 70). Fig 5.1 shows how GenOpt can be coupled with a simulation program.

Fig. 5.1 Coupling a simulation program with GenOpt (Wetter 2015)

As seen in Fig. 5.1, GenOpt consists of two parts: simulation and optimization. The simulation part reads the input files, stores the simulation results, and writes necessary

(46)

28

output files, while the optimization part contains the algorithms, controls the optimization process and stops it when the convergence criteria is fulfilled. There are several steps to do an optimization with GenOpt (Wetter 2009, pp. 69-70):

Firstly, to carry out an optimization it should be defined, where the input and the output files should be stored and overwritten. That’s the initialization part. Then user should select the optimization algorithm and related parameters in a “command” file. User can specify in a “configuration” file, how the simulation should end, what the error indicators should be etc. Lastly user should define an input template, which in turn will be overwritten by GenOpt. The input template will be the real input file after the optimization stops. GenOpt overwrites the keywords assigned to the parameters in the template file by replacing them with numerical values obtained from optimization results.

It is worth mentioning that the optimization constraints can only be implemented as penalty or barrier functions at GenOpt (Wetter 2009, p. 15).

5.2. Optimization Algorithm

The debate over how to select an optimization algorithm for district heating systems has drawn a large amount of attention. It has been asserted one should make a detailed analysis of the optimization problem before selecting the algorithm. Type of input variables (discrete, continuous or both) and constraints (equality or inequality) are the main decision factors. Moreover, nature of the cost function, whether it is linear, nonlinear, continuous, discontinuous, or has a local minima etc. should be analyzed (Nguyen et al. 2013, p. 1050). As mentioned earlier, there are many ways of solving discrete optimization problems. However at GenOpt, only Particle Swarm Optimization (PSO) algorithm can be used for solving discrete optimization problems. It is one of the stochastic optimization methods, which is based on simulation of a social behavior of fish schools, bird flocks or other biological groups, which, in general, take action by information exchange among members in a group (Nema et al. 2008, p. 1412). As other stochastic algorithms the PSO does not require gradient of the cost function.

Scientist interested in investigating how a bird swarm can flock synchronously, how they change direction, scatter and get together suddenly etc. Research in this area showed that, the behavior of a fish school or a bird swarm is not a coincidently but a deliberately taken action. As cited in (Kennedy and Eberhart 1995, p. 1943), a sociobiologist E.O Wilson has stated, “In theory at least, individual members of the school can profit from the discoveries and previous experience of all other members of the school during the search for food items, whenever the resource is unpredictably distributed in patches.” It is further claimed that birds are trying to keep the optimum distance between them, so that they could communicate and

(47)

29

orient themselves quickly to a location where they can find food or where they can escape from a predator (Kennedy and Eberhart 1995, p.1943).

In 1995 Kennedy and Eberhart have been inspired by these findings and introduced the PSO, which has roots in swarming intelligence, particularly bird flocking and fish schooling (Pal et al. 2011, p. 663). It is a population based optimization algorithm. In the PSO, population is called swarm, and the individuals are called particles. A swarm is regarded as a disorganized collection of individuals, which tend to cluster together.

In the PSO, each particle has a velocity and a position vector. Position vectors can be regarded as solution vectors of the optimization. The algorithm aims to find the best position vector that is reached throughout all iterations. The positions vectors are function of both velocity vectors and the previous position vector of the particles.

Consider e.g. n-dimensional optimization problem and let the set of solutions be 𝑥 ∈ ℝ𝑛 (particles’ position vectors), with 𝑥_𝑖𝑘 denoting 𝑘th position vector component of the particle 𝑖. Similarly 𝑣_𝑖𝑘 represents the velocity vector component of the particle with the same notation. At each 𝑐th iteration, the particles update their velocity and position vectors according to following formulations found in (Nema et al. 2008) and (Pal et al. 2011):

𝑣_ik(𝑡) = 𝑣_ik(𝑡−1)+ 𝑐c∙ 𝑟1(𝑡−1)∙ �𝑑ik(𝑡−1)− 𝑥ik(𝑡−1)� + 𝑐s∙ 𝑟2(𝑡−1)∙ �𝐺ik(𝑡−1)− 𝑥ik(𝑡−1)� , (𝟓. 𝟏) 𝑥_ik(𝑡)= 𝑥_ik(𝑡−1)+ 𝑣_ik(𝑡−1) , (𝟓. 𝟐) 𝑑_ik(𝑡)= 𝑚𝑟𝑘 𝑚𝑖𝑛 𝑥_j(𝑙):0≤𝑙≤𝑡−1�𝑓 �𝑥j (𝑙)_��₄_{, (𝟓. 𝟑)} 𝐺_k(𝑡)= 𝑚𝑟𝑘 𝑚𝑖𝑛 𝑃_ik(𝑡):1≤𝑖≤𝑁�𝑓 �𝑑ik (𝑡)_{�� , (𝟓. 𝟒)}

where 𝑐_c 𝑚𝑛𝑑 𝑐_s are acceleration constants; 𝑟₁ 𝑚𝑛𝑑 𝑟₂ are random numbers distributed in [0,1] that attain the stochastic swarm behavior; 𝑓 is the cost function;𝑑_ik is the 𝑘th component of the best position vector that particle 𝑖 reaches by its own experience; 𝐺_k is the 𝑘 th component of the best position vector reached in the entire swarm of N particles. Fig 5.2 indicates how the algorithm works based on the equations above.

4

(48)

30

Fig. 5.2 PSO algorithm (own figure inspired by (Nema et al. 2008, p. 1412))

As can be seen optimization is terminated when maximum number of iterations is reached. Then the best position vector 𝐺(𝑡) is specified as the output of the PSO algorithm.

Supervisory control of a combined heat and power plant by economic optimization