District heating and cooling systems – Framework for Modelica-based simulation and dynamic optimization

Accepted Manuscript District heating systems - Framework for Modelica-based simulation and dynamic optimization Gerald Schweiger, Per-Ola Larsson, Fredrik Magnusson, Patrick Lauenburg, Stephane Velut PII:

S0360-5442(17)30869-1

DOI:

10.1016/j.energy.2017.05.115

Reference:

EGY 10923

To appear in:

Energy

Received Date: 4 November 2016 Revised Date:

24 April 2017

Accepted Date: 17 May 2017

Please cite this article as: Schweiger G, Larsson P-O, Magnusson F, Lauenburg P, Velut S, District heating systems - Framework for Modelica-based simulation and dynamic optimization, Energy (2017), doi: 10.1016/j.energy.2017.05.115. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

RI PT

District heating systems - Framework for Modelica-based simulation and dynamic optimization Gerald Schweigera,∗, Per-Ola Larssonb , Fredrik Magnussonb , Patrick Lauenburgc , Stephane Velutb a AEE

M AN U

SC

- Institute for Sustainable Technologies, 8200 Gleisdorf, Austria b Modelon AB, SE-223 70 Lund, Sweden, c Department of Energy Sciences, Lund University, SE-221 00 Lund, Sweden

Abstract

Future district heating systems (so called 4th generation district heating (4GDH) systems) have to address challenges such as integration of (de)centralized renewable energy sources and storage, low system temperatures and high fluctuation of the supply temperature. This paper presents a novel framework for repre-

D

senting and simplifying on-grid energy systems as well as for dynamic thermohydraulic simulation and optimization of district heating systems. We describe

TE

physically precise and numerically robust models for simulation and continuous optimization. Furthermore, we propose a novel method to decompose a mixedinteger-optimal control problem into two sub-problems, separating the discrete

EP

part from the continuous. Two use cases show the applicability of the framework. An existing district heating system with more than 100 consumers is

AC C

adapted to test the framework based on simulation requirements of 4GDH systems. The second case presents the continuous optimization of a district heating system in a virtual city district. A main advantage of combining equation-based modelling and nonlinear optimization is the possibility of including model coherences based on physical laws into the optimization formulation. Results show that the framework is well-suited for simulating larger scale 4GDH systems and that the solution time of the continuous optimization problem is sufficiently low ∗ Corresponding

author Email address: [email protected] ( Gerald Schweiger)

Preprint submitted to Energy

May 18, 2017

ACCEPTED MANUSCRIPT

RI PT

for real-time applications. Keywords: District Heating, Smart Energy Systems, Modelica, Mixed-Integer-Optimal Control, Dynamic Optimization,

SC

1. Introduction

A major challenge for future energy systems is the design of systems that integrate large shares of fluctuating renewable inputs while improving the overall

M AN U

system efficiency. As fluctuating energy sources such as wind and solar energy expand, other parts of the energy systems must become more flexible to match the available energy from renewable resources with the demand in terms of location, time and quantity. There are a number of options for increasing energy system flexibility, including combining different energy domains, increasing supply and demand flexibility, integrating energy storage technologies and increasing the transmission capacity of the national grid as well as intercon-

D

nections to other countries. Previous research has shown that district heating

TE

infrastructure has the potential to play a key role in sustainable energy systems [1, 2, 3, 4]. In order to enable district heating as an integrated part of an energy system, the present generation of district heating systems will have to be fur-

EP

ther developed into a 4GDH [1, 4]. The simulation and optimization of concrete technical solutions provide the necessary insights and information to support the transformation process towards sustainable energy systems.

AC C

Distinctions between different tools and methods for simulation and opti-

mization can be made on a number of different levels. Tools can be divided into (i) single-domain tools that adress only specific aspects,(ii) multi-domain tools and (iii) urban design and planning tools [5]. Single-domain tools are often restricted and it can be cumbersome to adopt existing and create new models. Multi-domain tools can be further divided into ”general tools” (e.g. Modelica, Matlab, TRNSYS) and ”specific packages” (e.g. SynCity, CitySim) [5]. The advantages of general tools, besides their flexibility, are that most of these tools are structured in libraries. This enables an exchange of ideas and

2

ACCEPTED MANUSCRIPT

RI PT

methods within the scientific community and industrial developments which is a major advantage especially in an interdisciplinary, rapidly evolving research

field. Another advantage is that some of these tools (e.g. Dymola, TRNSYS, Matlab) provide the possibility of co-simulation using the FMI standard (Func-

SC

tional Mock-up Interface)[6].

Another general distinction can be drawn between block diagram modelling (also known as causal modelling) using imperative programming languages and

M AN U

equation-based modelling (also known as acausal modelling). A discussion of differences, applications and pros and cons can be found in [7, 8, 9]. According to [7, 8], equation-based languages are most advantageous for (i) representing the physical structure of systems. It enables the modeller to model the system directly by the means of differential algebraic equations (DAEs). The causality of how to solve the equations is not decided during the modelling stage, which

D

means that the usual need for manual conversion of equations to a block diagram is removed. (ii) Reusability, extensibility, adaptability of models and (iii) simple

TE

to code and to read the code. A disadvantage of these methods is that it is more difficult to go from a mathematical model to the numerical solution algorithm. Highly elaborate tools are required to handle them efficiently. Equation-based

EP

modelling languages/tools are Modelica [10], Neutral Model Format [11], EES [12] or Simscape within the Matlab/Simulink environment.

AC C

According to [13, 9], equation-based languages provide the following advantages for optimization: (i) automatic conversion of simulation models into optimization problems is supported; this makes it more easily applicable for users and it reduces the modelling effort. (ii) They can provide analytic expressions for gradients to be used by Non-Linear Program (NLP) solvers and (iii) the infinite dimensional problem can be automatically converted into a finite dimensional problem using collocation method. Previous research has shown that it is difficult to use tools that are based on imperative languages for optimization [14]. Methods and tools dealing with optimization applied to energy systems can be classified into (i) structural optimization: the main output is the existence of 3

ACCEPTED MANUSCRIPT

RI PT

each component and technology and the spatial layout of the network. Theses kind of optimizations take place as soon as possible in urban planning projects.

(ii) Optimization of design parameters and (iii) Optimization that deals with the

operative control [15, 16]. The first class of methods should be able to analyse existing systems as well as systems that are radically different [17]. Limitations

SC

of the first and second class of methods are that they only consider dynamics to a very limited extent. Sophisticated numerical methods of the third class have often inconvenient application programming interfaces so that such methods are

M AN U

inappropriate for engineering applications [13]. In order to design future energy systems, an interaction of all three methods is required. Structural optimization precedes parameter optimization and operative control.

1.1. Prior research

D

This section provides a literature review that focuses on ”general tools” in the category of multi-domain tools for simulating district heating systems as

TE

well as optimization approaches for district heating and cooling networks of the second and third class.

EP

1.1.1. Modelling

A review of different modelling approaches and tools for district heating systems is given by Olsthoorn et al. [18]. Different simulation tools are compared

AC C

with respect to level of precision as well as computational limitations. Allegrinig et al. [5] provides a review of modelling approaches and tools for the simulation of district-scale energy systems including heat networks, multi-energy systems, renewable energy generation and urban microclimate. Soons et al. [19] made a functional comparison of the simulation environments of Modelica/Dymola, Matlab/Simulink and TRNSYS. The comparison indicates that Modelica performs better in terms of multi-domain modelling, realistic control, modularity and flexibility. Strunik et al. [20] presented a district heating systems model implemented in Matlab using an artificial neural network. An energy and an

4

ACCEPTED MANUSCRIPT

RI PT

exergy-based method were developed to compute the mony flows. They conclude that exergy-based methods give a more realistic picture of the heating cost evaluation. Duquette et al. [21] presented a pipe model implemented in Matlab/Simulink using a steady state heat transfer model combined with a

variable transport delay model. The model is validated using measured data

SC

from a district heating grid and compared to with a transient model regarding computation time. Kruchi et al. [22] developed models in IDA-ICE for simulating district energy-systems including waste-heat, models for heat storage in

M AN U

the ground and pipe models for loops. The pipe is modelled through a perfect mixed water that has the size of the volume within the pipe. Li et al. [23] developed a model of a district heating system in Simscape to analyse heat losses in the distribution network. Approaches that present a quasi-dynamic approach, where the temperature is calculated dynamically while the pressure and flow are calculated in steady state conditions, are found in [24, 25]. The objective of

D

the IEA project Annex 60 is to develop new generation computational tools for building and community energy systems based on the Modelica and FMI stan-

TE

dards [26]. Giraud et al. [27] developed a Modelica based library for district heating simulations. They could reduce the number of equations in the pipe model by a factor of 40 compared to the pipe model of the Modelica Standard

EP

Library. TransiEnt Library [28] is a freely available Modelica based framework to model coupled energy supply grids. However, studies that demonstrate the

AC C

applicability of Modelica as a suitable modelling language to describe and efficiently simulate large-scale district heating systems as well as 4GDH systems, are still lacking. 1.1.2. Optimization Sameti et al. [16] reviewed different optimization approaches and tools for

district heating and cooling networks. A survey of available models and optimization approaches for production planning of district heating systems can be found in [29, 30]. Powell et al. [31] describe a methodd for dynamic optimization to find optimal charging and discharging periods for thermal storage which

5

ACCEPTED MANUSCRIPT

RI PT

is used to shift electrical and cooling loads. The dynamic optimization problem is solved by decomposing the problem into multiple static mixed-integer non-

linear programming (MINLP) problems. A hybrid evolutionar mixed integer linear program (MILP) optmization algorithm is used by Vesterlund et al. [32]

to minimize he total operating costs of district heating networks. Morvay et

SC

al. [33] proposed a MILP model for multi-objective optimisations: structural, design and operational optimizations were performed for distributed energy systems including heating networks. Grosswindhager et al. [34] describe a method

M AN U

for controlling the supply temperature in district heating systems using Fuzzy Direct Matrix Control implemented in Matlab. The model used for control was simplified by aggregating the original network using the Danish method [35]. Sandou et al. [36] present a method for predictive control of complex district heating networks using Sequential Quadratic Programming algorithm implemented in Matlab. It is assumed that the delay is constant for each pipe.

D

Veserlund et al. [37] present a method for optimization of district heating systems with meshed networks. In the first step the district heating system is

TE

simulated in Matlab/Simulink to get the overall thermal energy that has to be supplied to satisfy the consumers demands. This result is an input for the MILP, which calculates an optimal schedule for different plants. Orehounig et

EP

al. [38] extended the energy hub approach ( [39]) with simplified models for district heating systems to analyse decentralized energy systems at a neigh-

AC C

bourhood scale. The problem is formulated as linear Program. A so-called energy hub represents a single entity containing all possible conversions and storage technologies which describes and manages the relation between input and output energy flows. Jiang [40] proposes an optimal operating strategy for district water-heating systems by applying Group Search Optimizer. Guelpa [41] presents an optimization method for minimizing the pumping cost in large district heating networks. A generic algorithm is applied across all different model complexities. A methodology that may be beneficial for larger networks is the Lagrangian relaxation, in which the system-wide constraints are relaxed and the overall 6

ACCEPTED MANUSCRIPT

RI PT

problem is decomposed into a master problem and local unit-specific problems. Due to the separability requirement of the system constraints, it is difficult to

include the distribution network in the formulation. According to [42] and [43] that contain reviews of available softwares for production planning, the most common approach to optimize district heating systems is based on piecewise

SC

linear plant models and a MILP formulation. This results in large-scale integer optimization problems that can efficiently be solved with the available

solvers on the market. The planning methods adopting the MILP formulation

M AN U

contain heavily simplified models, often only modelled as algebraic equations with exceptions for storing dynamics (heat and fuel) and eventually delays in the distribution network. The continuous optimization variables are typically heat flows and the effects of the supply temperature and flow as well as return temperature are not directly considered. This is limiting as these temperatures and flows affect many critical parameters such as amount of energy that can be

D

stored in the net, heat losses in the net and efficiencies for electricity production in steam turbines, see [44] for the large impact of the supply temperature

TE

choice. The impact of critical physical variables cannot be easily and accurately captured without physical models of the plant and distribution network.

EP

1.2. Main contributions

As reported in Allegrini et al. [5], there is much to be done to explore the full

AC C

benefit of innovative district energy systems, including district heating systems. They argue that a shift to fully dynamic models and sophisticated control design would be supportive. Requirements given by the structure of 4GDH systems call for multi-domain tools that must capture all important dynamics within the entire system. We present a novel integrated framework for fully dynamic, thermo-hydraulic

simulation and optimization of 4GDH systems. The framework is based on the high-level equation-based modelling language Modelica and a high-level, largescale continuous dynamic optimization method implemented in JModelica.org and OCT (OptimicaCompiler Toolkit). The discrete optimization method, the 7

ACCEPTED MANUSCRIPT

RI PT

unified network representation and the aggregation algorithms are implemented in Python. The main contributions of this paper are: (i) A unified network rep-

resentation for creating and manipulating of grid-based energy systems. These networks can be automatically translated into executable code for simulation and optimization. (ii) Demonstration of the capabilities of the library for sim-

SC

ulation and continuous optimization of district heating systems implemented

in Modelica. (iii) A method for decomposing a mixed-integer-optimal control problem into two sub-problems: A Mixed Integer Quadratically Constrained

M AN U

Programming (MIQCP) problem and a continuous problem that is solved using a direct, local collocation algorithm. While the method section focuses on the decomposition of the mixed-integer-optimal control problem, the use-cases focus on demonstrating the capabilities of the framework for simulation and continuous optimization.

D

2. Method

TE

2.1. Languages and Tools

The framework presented in this article is based on the modelling language Modelica and the scripting language Python. The advantages of equation-based

EP

modelling together with the fact that Modelica is open source and has a large and fast growing community in different industrial applications and academia [45, 9, 46, 27, 19, 47] makes this language suitable for modelling 4GDH sys-

AC C

tems. Dymola was used to simulate Modelica models. JModelica.org and OCT were used to solve the dynamic optimization problem; these tools support the Modelica language extension Optimica that enables formulation of dynamic optimization problems based on Modelica models [13, 48]. OCT is based on JModelica.org technology, but has several additional features [49]. The compiler of JModelica.org/OCT performs in the first step symbolic transformations. Next, a symbolic representation of the problem is created by CasADI [50], which also provides the needed first- and second-order derivatives by automatic differentiation. This representation is then used to transform the

8

ACCEPTED MANUSCRIPT

RI PT

problem into a NLP using direct, local collocation (details are presented in Section 2.2.2). The direct collocation method transforms the infinite-dimensional

optimal control problem into a finite-dimensional NLP. The NLP is then solved using the interior-point algorithm IPOPT [51] to find a local optimum, since problems in context of this research field are in general nonconvex. Different

SC

NLP solvers are compared in [52]. The discrete problem is implemented using

the open source Python module Pyomo [53] (details are presented in Section 2.2.1). Pyomo is an algebraic modelling language where optimization problems

M AN U

can be formulated similar to mathematical notations. Pyomo has been built to interface with different solvers; in the present work we use the commercial solver Gurobi [54]. 2.2. Optimization method

Production planning of district heating networks can be formulated as an

D

optimization problem, involving both discrete and continuous decision variables, transport delays, nonlinear and dynamical behavior. The resulting optimization

TE

problem is known as a mixed-integer-optimal control problem that can be discretized to become a nonlinear programming problem (MINLP) for which no general, robust and scaleable approach exist [55]. Approaches to solve MINLP

EP

problems do normally not solve simultaneously the discrete and continuous variables. Instead, a related MILP, e.g., Outer Approximation or General Bender

AC C

Decomposition, or a related NLP, eg., Branch and Bound, see [5], is iterated over. There has been recent work on adopting collocation methods for hybrid cases including continuous as well as discrete parts [56, 57]. It is expected that these approaches do not scale-up as well as the proposed decomposition. The authors’ idea is to decompose the MINLP problem into two sub-problems separating the discrete part formulated as MIQCP problem from the continuous part formulated as a NLP (see Figure 1). The sub-problems are normally referred as unit commitment problem (UCP) and Economic Dispatch problem (EDP). The district heating system is simulated with a single unit in order to calculate the overall energy that is required to satisfy consumers’ demand. The MIQCP 9

ACCEPTED MANUSCRIPT

RI PT

method is then used to calculate the optimal schedule for all production units based on a given cost function (similar as in [37]); the status for each producer

is the main result. The optimization horizon within district heating systems can

range from daily to several days. The optimization horizon of the UCP is usually

longer compared to the EDP since start-up and shutdown trajectories have to

SC

be taken into account and the UCP has usually good scale-up properties. The aim of the EDP is to optimize the system based on physical models of the entire

system. Smoothened status trajectories of the units are given by the solution

M AN U

of the UDP; this yields a continuous optimization problem, which is then transformed into a NLP. The major advantages of including model coherences based on physical laws into the optimization formulation are high accuracy and the possibility to impose constraints on physically and operational relevant variables such as temperature, pressure or mass flow, based on physical and operational limitations of the real system. Both problems require a predicted heat load over

D

the entire optimization horizon, which is assumed to be perfect known for this

AC C

EP

TE

work.

Figure 1: Schematic view of the proposed method to decompose the MINLP problem.

10

ACCEPTED MANUSCRIPT

RI PT

2.2.1. Discrete optimization problem The method for solving the discrete problem is based on [58]. The discrete

problem contains integer variables (units on/off), which are in this case binary

variables. The problem can be classified as ”quadratically constrained” since

the problem contains products of variables in the objective function and the

SC

constraints. The problem is to set the amount of heat to be produced at each

time step and for each producer; we get two indexed decision variables: the heat production (qi,k ) and the status of the producer (ui,k ). Two additional

M AN U

binary variables are introduced to formulate start-up and shutdown conditions: start-up variables (yi,k ) and shutdown variables (zi,k ); i is a time step index and k specifies the producer.

The objective function has the form: f (ui,k , qi,k , yi,k , zi,k ) =

I X X i=1

i=1

(1)

D

k∈K

 Xk
I+DD (cfk +cvk )ui,k qi,k +cfk o ·ui,k +Ckup ·yi,k + Ckdown ·zi,k

where cfk is the fuel cost, cvk is a cost proportional to produced heat, cfk o is a

TE

fixed operation cost, Ckup is a start-up cost, Ckdown is a shutdown cost, K is the set of producers and DDk is the shutdown duration.

AC C

EP

The demand constraint is given by: X

ui,k qi,k ≥ qi,D , ∀i ∈ [1, I].

(2)

k∈K

where qi,D is the predicted demand for time step i (based on simulation re-

sults, see Figure 1). The formulation of the start-up and shutdown constraints are based on [59]; the start-up and shutdown trajectories are assumed to be linear, but they could be adjusted to fit the characteristics of a particular producer. 2.2.2. Continuous optimization problem The continuous optimal control problem is solved using JModelica.org’s direct, local collocation algorithm to transform the infinite-dimensional optimal 11

ACCEPTED MANUSCRIPT

RI PT

control problem into a finite-dimensional NLP[48]. See [60] for an overview of different numerical methods for optimal control. This section first describes the

general continuous optimal control problem formulation and then the main steps of the direct collocation used in JModelica.org to transform it into an NLP. The

x,y,u,p

F (t, x(t), ˙ x(t), y(t), u(t), p) = 0,

x(t0 ) = x0 ,

(3b) (3c)

L ≤ (x(t), ˙ x(t), y(t), u(t), p) ≤ U, he (t, x(t), ˙ x(t), y(t), u(t), p) = 0,

He (tf , x(tf ), y(tf ), p) = 0, ∀t ∈ [t0 , tf ].

(3a)

M AN U

s.t.

t0

SC

general continuous optimization problem formulation is Z tf
min. φ(tf , x(tf )) + L x(t), y(t), u(t) dt,

hi (t, x(t), ˙ x(t), y(t), u(t), p) ≤ 0,

Hi (tf , x(tf ), y(tf ), p) ≤ 0,

(3d) (3e)

D

The system dynamics are described by the implicit, index-1 DAE in (3b), where

TE

t ∈ [t0 , tf ] is time, x is the state, y is the algebraic variable, u is the control variable and p is the vector of parameters to be optimized. The objective function is the Bolza functional (15), comprising the Mayer term φ and Lagrange inte-

EP

grand L. In the constraints, the simple bounds (3c) have been separated from the general path constraints (3d), which can be on both equality and inequality

AC C

form. Finally, there are terminal equality and inequality constraints (3e).

The collocation method starts by dividing the time horizon into ne elements.

The time in element i is normalized: tei (τ ) = ti−1 + hi · τ,

∀τ ∈ [0, 1],

(4)

where tei (τ ) is the unnormalized time, τ the normalized time and ti is the right boundary of element i. The time-dependent variables x, y, u are approximated by polynomials in the local time τ using nc collocation points as interpolation

12

ACCEPTED MANUSCRIPT

xi (τ ) =

nc X

xi · lek (τ ),

yi (τ ) =

k=0

nc X

yi · lk (τ ),

RI PT

points: ui (τ ) =

k=1

nc X

ui · lk (τ ),

(5)

k=1

where lek is the Lagrange basis polynomial including the first point and lk without

SC

the first point. The addtional interpolation point for lek has been introduced to enforce continuity of x. The basis polynomials are defined by lek (τ ) =

Y

M AN U

l=0,...,nc , l6=k

τ − τl , τk − τl

lk (τ ) =

Y

l=1,...,nc , l6=k

τ − τl , τk − τl

(6)

(7)

where τk are the collocation points chosen according to the Radau scheme. The polynomials are determined by enforcing all the constraints at the col-

D

location points. In particular, the DAE is transformed into the NLP constraints i = 1, . . . , ne ,

k = 1, . . . , nc .

(8)

TE

F (ti,k , zi,k , p) = 0,

To get the collocation polynomial for x˙ the collocation polynomial xi is

EP

differentiated using (4) and (5) to obtain x˙ i (τ ) =

nC dτ dxi 1 X dlek dxi (τ ) = (τ ) = xi,k · (τ ). hi dτ dtei dtei dτ k=0

(9)

AC C

Based on these discretizations, the continuous optimal control problem is transformed into an NLP, whose solution yields an approximate solution to the optimal control problem. 2.3. Use cases The following section presents two use cases that demonstrate key applica-

tions of the framework: (i) A scale-test for simulating 4GDH systems and (ii) a continuous dynamic optimization of a district heating system in a city district. Both cases are adopted to stress the framework regarding numerics and robustness.

13

ACCEPTED MANUSCRIPT

RI PT

2.3.1. Scale-test for simulating 4GDH systems The scale test is a typical district heating system located in Austria (Figure 2). The system supplies 103 consumers and has a total length of about 14.000

m. The base unit (seen on the left side of the Figure 2) has a nominal thermal capacity of 4 MW. The top unit (seen on the upper right side) has a nominal

SC

thermal capacity of 750 kW. The decentralized solar unit seen on the lower right

side has a nominal power of 55 kWp. The prosumers have solar units with a

M AN U

nominal power in the range of 5 − 12 kWp. Prosumer

Top Unit

C2

Base Unit

C1

System details:

D

Consumer: 103 Total length: ~14.000m Loop base unit: 4MW Top Unit: 750kW Solar Unit: 55kWp Solar Unit Prosumer: 5-12kWp

TE

Solar Unit

EP

• • • • • • •

Figure 2: The scheme of the district heating system (including information about the network and the production untis) with the production untis, prosumers and two consumers (C1, C2)

AC C

are highlighted.

In the base case, the system is simulated with only the base unit to obtain

physical and numerical references. The 4GDH case was designed to test future simulation requirements that are probably given by completely new system settings. Based on [1, 61], the following requirements were defined: decentralized producers, meshed networks, integration of prosumers, temperature waves and high transients. In order to test these requirements we adopted the base case as follows: a decentralized top unit that operates with high gradients during start-up and shutdown for peak load; a temperature step at the base unit was

14

ACCEPTED MANUSCRIPT

RI PT

introduced to simulate a temperature wave. A sophisticated control strategy based on the minimum differential pressure in the entire network was implemented since the critical node could vary during operation.

2.3.2. Continuous dynamic optimization of district heating systems

SC

The dynamic optimization case represents a virtual city district that con-

sists of one producer, 16 consumers and a total length of about 4.200m (Figure 3). We assumed a perfect load prediction over the entire optimization horizon.

M AN U

This case was designed to test (i) the accuracy of the aggregation method, (ii) mass flow dependent delays in dynamic thermo-hydraulic optimizations even for small changes in the mass flow and (iii) the integration of physically relevant constraints (e.g. temperatures and mass flows) based on physical and

AC C

EP

TE

D

operational limitations in the optimization formulation.

Figure 3: The scheme of the district heating system of a virtual city distric. The producer is seen on the left side; the black/blue circles represent the 16 consumers.

3. Overview of proposed framework A schematic view of the framework is presented in Figure 4 and each step for the most complex case, namely the MINLP (steps 1-7 in orange rectangles 15

ACCEPTED MANUSCRIPT

RI PT

in Figure 4) is explained below. • Step 1: The network is created using a unified network representation that includes data of the network, demand and boundary conditions.

• Step 2: The unified network representation is translated into executable

SC

Modelica code. Two simulations are performed: (i) a steady state simulation and (ii) a simulation with the predicted demand for the entire optimization horizon. The first simulation is required for the aggrega-

M AN U

tion method which needs nominal values for mass flows, temperatures and pressure. Results of the second simulation are required as an input for the MIQCP.

• Step 3: The MIQCP is solved. The result is the status and heat production of each unit. The dynamic optimization requires only the discrete variables since the actual heat production of each unit will be a free variable within

D

the dynamic optimization.

TE

• Step 4: The original network is aggregated. Detailed dynamic thermohydraulic optimization of large-scale district heating systems is numerically challenging. Thus, the network topology needs to be simplified to a

EP

size suitable for optimization using so called aggregation methods. Previous works on aggregation methods show that networks can be aggregated up to a very high level even in dynamic operations without losing signifi-

AC C

cant accuracy [35]. The framework includes the so called German method [62]. The aggregation depth is flexible and certain consumers can be excluded from the aggregation.

• Step 5: The aggregated network is simulated to get initial trajectories for the dynamic optimization. Solving large-scale nonconvex optimization problems requires accurate initial guesses of the solution to the problem for. For large scale problems it is troublesome to provide initial guesses manually. The system can be simulated to generate initial trajectories using only initial guesses for the degrees of freedom [48]. Smoothened 16

ACCEPTED MANUSCRIPT

RI PT

status trajectories of the units are given by the solution of the MIQCP. Additional constraints for each unit are lower and upper trajectories dur-

ing start-up, operation and shut-down. Thus, the actual heat production within the dynamic optimization has time varying minimum and maximum constraints also known as path constraints.

SC

• Step 6: The dynamic optimization problem is solved.

• Step 7: The optimal trajectories are applied on the model of the original

EP

TE

D

M AN U

network.

AC C

Figure 4: Schematic view of the framework showing the different steps for simulation, mixed integer problems and continuous problems.

4. Implementation 4.1. Unified network representation The core of the framework is a unified network representation. In general, the representation is applicable for all kinds of grid-based energy systems including district heating systems, gas systems and power systems. Typically such 17

ACCEPTED MANUSCRIPT

RI PT

systems consist of edges and nodes. Edges are transmission or distribution lines, gas pipes or pipes within a district heating system. Nodes are consumers and

producers in any energy domain, storage or hybrid technologies. Such a unified network representation is required for two reasons. Firstly, an automatically

generated simulation or optimization network model based on a unified network

SC

rep. Secondly, several steps require detailed information of the network topology (f.i. see step 2 in Section 4) and other steps change the topology of the network

(f.i. step 4 in Section 4). The unified network representation is implemented in

M AN U

Python and consists of three central modules: network-representation (module NetworkX [63]), aggregation and translation. 4.2. Models for Simulation

The following section presents the main developments in the dynamic thermohydraulic simulation library. Single component tests and validations are pre-

D

sented in [64]. The pipe model for simulation is implemented based on a plugflow approach. The one-dimensional energy balance of a pipe can be described

TE

by the following partial differential equation: ∂T ∂T 1 + v(t) + q(T (x)) = 0 ∂t ∂x Sρcp

EP

(10)

where v is the fluid velocity, S the cross-section area, ρ the fluid density, cp the specific heat capacity of the fluid and q the heat loss to the surroundings

AC C

of the pipe. It depends on the boundary conditions and is often described as a linear function of temperature difference between the boundaries. In that case, Equation 10 is linear and can be integrated explicitly. In the case of time-varying mass flows, it can be solved explicitly by applying the method of characteristics [65]. The temperature at the outlet is given by: − Tτp

T (x = L, t) = Tboundary + (T (x = 0, t − τ ) − Tboundary )e

Tp =

cp ρπR2 kS 18

(11)

(12)

ACCEPTED MANUSCRIPT

RI PT

where L is the pipe length, τ the time-varying delay, Tp is the characteristic decay time, k the heat conductivity of the surrounding and S is shape

factor for the pipe. The pipe model with heat loss is implemented based on a combination of a plug-flow approach and an ideally mixed volume model. The Modelica built-in operator spatialDistribution() [10] provides a robust method

SC

to approximate the solution of one-dimensional partial differential equations. The operator keeps track of the spatial distribution via suitable sampling, interpolation and shifting of the stored distribution [10]. The pipe model uses

M AN U

two spatialDistribution() operators. One for the one-dimensional plug-flow approach and one that uses time as the transported quantity to calculate the time-varying delay τ (Equation 12). The volume model is introduced for two reasons: (i) to avoid a non-linear system of equations when connecting two flowmodels, (ii) the resulting pipe model captures the dynamics of real pipe better, because there is a certain degree of mixing water inside a pipe. Since extremely

D

stiff systems should be avoided, a compromise between accuracy and numerical robustness has to be made.

TE

Different consumer models have been implemented. A consumer model always consist of three major parts: a valve, a heat exchanger and a controller. Two heat exchangers models are available; one that calculates the return tem-

EP

perature and one with a fixed return temperature which is set by the modeller. The valve model is a linear model that provides a pressure drop which is pro-

AC C

portional to the flow rate and to the opening. The opening of the valve is controlled by a PI-Controller or directly calculated based on the given heat demand trajectory. The second method was implemented to improve the simulation performance. A simplified, ideal producer model has been implemented. The overall control strategy for the entire system is based on a minimum differential pressure in the entire network or on a minimum differential pressure at a critical node. It is possible to simulate any district heating system - independent of the network topology - without introducing nonlinear equations or numerical Jacobians.

19

ACCEPTED MANUSCRIPT

RI PT

4.3. Models for continuous Optimization The work is based on previous work of [55, 46]. The formulation of the

optimization problem has to deal with restrictions due to the numerical method that is used: Hybrid constructs are excluded since the right hand side of the dynamics is assumed to be C 2 (twice continuously differentiable). Since this

SC

restriction excludes particular functions, smooth approximations of important functions (e.g. min, max, abs) are integrated in the optimization library.

M AN U

A novel feature is the integration of pressure models into all components. The pipe model contains a combination of a fixed time delay and a finite volume. The goal is to compute the main characteristics of the pipe without having to use a model with a large number of segments, which would increase model complexity. The fixed delay is dependent on the range of the mass flow for each pipe and corresponds to the minimal time delay. The finite volume must capture the flow-dependence of the varying mass flow. The number of segments

D

within the finite volume depends on the geometry of the pipe. The delay τ is the sum of the fixed delay and the delay in the finite volume model. The actual

TE

implementation of the pipe requires that the fixed time delay within the pipe model has to be an integer multiple of one element length. One has to set ne

EP

and nc to get a good compromise between complexity (total number of states) and accuracy (smallest delay). The fixed delay should not be zero for too many pipes. The consumer model is similar to the simulation model. The mass flow

AC C

through the valve is computed using the pressure difference, the valve opening and data from a nominal point. A standard quadratic equation relates mass flow and pressure drop. The opening of the valve is directly calculated based on the given heat demand trajectory. The heat exchanger within the consumer model calculates the heat flow rate using a fixed return temperature which is set by the user. A simplified, ideal producer model has been implemented. 4.4. Continuous Optimization formulation For good numerical properties of the optimization problem (e.g. scaling of the Jacobian), all variables in the optimization problem should be in the same 20

ACCEPTED MANUSCRIPT

RI PT

order of magnitude [48]. In one model, there could be pressures in the range of MPa, temperatures of 250-400 K and pressure ratios of less than 0.1. To avoid excessive differences in magnitude, all variables are scaled to nominal

trajectories. In Modelica, this can be done by using the nominal-attribute. The optimization library provides a unit-package including nominal values. Dynamic

SC

optimization problems are formulated using the Modelica extension Optimica [13]. The decision variables Uj (t) within the dynamic optimization could be f.i. supply temperatures, pressures and mass flows of the producers. In order

M AN U

to define minimum and maximum constraints on the derivatives of the decision variables, the following equation is introduced: Z

Uj (t) =

U˙ j (t)dt

(13)

t

Minimum and maximum constraints can be set in variables attributes min and max. The constraints could also be formulated in the constraint-section

D

which is available through the Optimica extension. However, the first approach

TE

is preferable because it effectively cuts off the search space for the optimization. The Optimica extension provides a class called optimization. Class attributes are objective, which defines the costF unction, startT ime and f inalT ime

EP

which defines the start and end, respectively, of the optimization horizon. A minor cost on the input derivatives is integrated in the cost function for regularity

AC C

reasons. This cost at a certain time instant is formulated as

W (t) =

X

rU˙ j U˙ j (t)

(14)

j=model inputs

where rU˙ j is the weight for derivative U˙ j .

This additional term within the cost function is a trade off and should there-

fore be small compared to the actual cost function. However, a too small ratio can yield a bad numerical performance. Preliminary tests have shown, that the ratio should be about 0.01. The cost function to be minimized has the general form:

21

ACCEPTED MANUSCRIPT

tf

min.

x,y,u,p

RI PT

Z


L x(t), y(t), u(t) dt

t0

(15)

where the optimization interval is [t0 , tf ] and L can be divided into several

different instantaneous costs, referring to physical relevant variables, economics,

SC

ecologic, regularization, etc.

5. Results

M AN U

All cases were run on a laptop with 8 GB RAM and four 2.6 GHz CPUs. 5.1. Scale-test for simulating 4GDH systems

For both cases (base case and 4GDH case), the supply temperature was set to 90 deg.C and the minimum differential pressure to 0.8 bar. The load trajectory is based on building simulations for different types of buildings, scaled with the actual nominal load for each building. The total daily heat load for the entire

EP

2500 2000

AC C

Load [kW]

3000

Total heat load

TE

3500

D

network can be seen in Figure 5.

1500

10000

4

8

12 Time [hours]

16

20

Figure 5: Total heat load

The return temperatures of the consumers are in the range of 60 − 68 deg.C and set to fixed values depending on the building type. Figure 6 shows the

22

24

ACCEPTED MANUSCRIPT

RI PT

pressure at the base unit for the base case as well as for the 4GDH case and the heat production of the top unit in the 4GDH case. From this figure it can be seen that the maximum pressure at the base unit in the base case is about 5.5

bar. The maximum pressure in the 4GDH case is about 1.5 bar lower because

BaseUnit (base case) BaseUnit (4GDH case)

M AN U

Pressure [bar]

5 4 3

10

4

ba

8

D

2

12 Time [hours]

16

TopUnit

1600 1200 800 400 20

TE

Figure 6: Pressure at the base unit for the base case (black, left y-axis) as well as the 4GDH case (black dashed, left y-axis) and the heat production at the TopUnit in the 4GDH case

EP

(orange, right y-axis)

Figure 7 shows the temperature wave propagation. The delay from the producer to consumer C1 (see Figure 2) is about 30 min, to consumer C2 about

AC C

45 min.

Figure 8 shows the change of mass flow direction in a pipe close to the

consumer C2. The pipe implementation handles flow changes as well as the zero flow efficiently and without introducing numerical noise. Table 1 shows the simulation statistics of the two cases. The significant

increase of the physical and engineering systems complexity has little impact on the simulation performance. According to these results we can conclude that the framework is suitable for simulating 4GDH systems.

23

2000

Heat Production [kW]

6

SC

of the top unit which operates during the peak load.

240

ACCEPTED MANUSCRIPT

100

96 94 92

SC

Temperature [deg C]

98

RI PT

T_supply BaseUnit T_supply consumer C1 T_supply consumer C2

90

860

M AN U

88 4

8

12 Time [hours]

16

20

24

Figure 7: Temperature wave propagation. Supply temperature at producer (orange), at con-

D

Mass flow close to consumer C1 Mass flow close to consumer C2

TE

9 8 7 6 5 4 3 2 1 0 1 20

AC C

EP

Mass flow [kg/s]

sumer C1 (black dashed) and at consumer C2 (black with marker).

4

8

12 Time [hours]

16

20

Figure 8: While the flow direction in a pipe close to consumer C1 (see Figure 2) is not changign, the flow direction in a pipe close to the consumer C2 is changing two times during the simulation.

5.2. Continuous dynamic optimization of district heating systems For the sake of simplicity, the objective of the dynamic optimization was to minimize the supply temperature at the producer. The formulation contains two

24

24

ACCEPTED MANUSCRIPT

2764

Variables

49509

Equations

21533

Nontrivial

17035

Nonlinear systems

0

Nummerical jacobians

0

Approx. computation time [min]

8

2875

50816 22093

17511

0

SC

Components

4GDH Case

RI PT

base case

0

M AN U

9

Table 1: Simulation statistics of for the base case and the 4GDH case

degrees of freedom, namely the pressure derivative and the temperature derivative at the producer. Thus, the dynamic optimization problem is formulated as

tf

D

Z

min

TSupply (t) + W (t)dt

(16a)

TE

t0

s.t.

F (t, z(t), p) = 0,

F0 (t0 , z(t0 ), p) = 0

mP rod (t) ≤ mP rodM ax

∀t ∈ [t0 , tf ]]

(16b) (16c) (16d)

dpL Consumer ≤ dpConsumer (t) ∀t ∈ [t0 , tf ]]

(16e)

EP

L TConsumer ≤ TConsumer (t) ∀t ∈ [t0 , tf ]

AC C

where mP rodM ax is the upper limit of the mass flow at the producer (set

to 50 kg/sec); this constraint was introduced to force the supply temperature to increase during the peak loads. This should demonstrate the temperature L propagation within the dynamic optimization. TConsumer is the lower limit

of the supply temperature for all consumers and it was set to 62 deg.C. The lower limit of the differential pressure for all consumers (dpL Consumer ) was set to 0.5 bar. The return temperature at the consumers was set to 43 deg.C. Minimum supply water temperatures and pressure differences for all consumers were introduced based on real network limits to satisfy the consumers’ demand.

25

ACCEPTED MANUSCRIPT

RI PT

Each consumer has a given nominal load. To calculate a load trajectory, a typical residential heat load profile with two daily peaks was assumed (Figure 9).

1.0

Nominal load trajectory

SC

0.8

M AN U

0.6 0.4 0.2 0.00

4

8

12 Time [hours]

16

20

D

Figure 9: Nominal load trajectory

TE

For the dynamic optimization, the original network was aggregated to a depth of 87.5 % which means that 12.5% of the original number of elements remained (see Figure 10). After aggregation, only two consumers were left.

EP

The optimization horizon was 24 hours. The number of elements within the collocation method was set to 100, i.e., the length of each element is 864 seconds.

AC C

As shown in Figure 11, the requested and the obtained load are identical, that means that the consumers receive exactly the amount of load demanded. Figure 12 shows the mass flow dependent delay between the remaining con-

sumers. As expected, the delay increase when the mass flow decrease and vice versa (the mass flow is shown in Figure 13). As expected from Equation 16a, the supply temperature increases only as

the mass flow reaches its constraint of 50 kg/sec (Figure 13). Figure 14 shows a comparison of the results of step 6 and step 7 in 3 (see the blue rectangles ”E” and ”F” for ”Continuous Dynamic Optimization Problem” in Figure 4). The idea is to validate the optimization results and the aggrega26

24

M AN U

(a)

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

Figure 10: Modelica models generated in Python (a) Original network. (b) Network after aggregation.

2500 2000

D

EP

Load [kW]

3000

Consumer 1 load requested Consumer 1 load obtained Consumer 2 load requested Consumer 2 load obtained

TE

3500

1500

AC C

10000

4

8

12 Time [hours]

16

20

Figure 11: Requested and obtained load of the remainging consumers.

tion method by applying the optimal trajectories to the original network with 16 consumers. The pressure trajectory would not be applied in real systems. Thus, the pressure at the production unit is controlled based on the minimum differential pressure, which was set to the same value as dpL Consumer in the optimization

formulation. Only the optimal supply temperature trajectory is applied to the

27

24

ACCEPTED MANUSCRIPT

Delay between remaining two consumers

RI PT

4190 4180

4160

SC

Delay [sec]

4170

4150

41300

4

M AN U

4140 8

12 Time [hours]

16

20

24

Figure 12: Delay between the two remaining consumers.

Producer mass flow

49.6 49.4

D

72 70

TE

49.8

EP

Mass flow [kg/sec]

50.0

74

T_supply Producer

Temperature [deg C]

50.2

68 66

49.2

AC C

49.00

4

64 8

12 Time [hours]

16

20

2462

Figure 13: Mass flow (black dashed, left y-axis) and temperature (orange,right y-axis) at the producer based on the dynamic optimization result.

original network. We refer to this simulation result as the ”validated case”. The orange and grey solid lines show the supply temperatures of the remaining two consumers based on the optimization result (Figure 10b). The two dashed lines show the supply temperatures for the same consumers in the validated case (Figure 10a). The temperature difference at steady-state (i.e between

28

ACCEPTED MANUSCRIPT

69 67 66 65

SC

Temperature [deg C]

68

T_supply Producer T_supply Consumer 1 Sim T_supply Consumer 1 Opt T_supply Consumer 2 Sim T_supply Consumer 2 Opt

RI PT

70

64 63 610

M AN U

62 4

8

12 Time [hours]

16

20

Figure 14: Temperatues at the producer (black), Consumer 1(orange) and Consumer 2 (grey with marker) based on dynamic optimization result. Temperatures in the validated case: Consumer 1 (orange dashed) and Consumer 2 (grey dash-dotted).

D

hour 12 and 14) are due to introduced errors from the aggregation method and different model complexities in the simulation and optimization models. From

TE

the same figure, it can be seen that the heat losses and transport delays are correctly captured within the optimization: the supply temperature at the pe-

EP

riphery (Consumer 2) is lower and has some delay compared to the consumer close to the producer (Consumer 1). Figure 15 shows the difference between the result of the optimization and the validated case. From this figure it can be

AC C

seen, that the maximum difference between both remaining consumers is less than 0.7 K.

The differences of the optimization results and the validated case of the heat

production as well as the mass flow at the producer are presented in Figure 16. The maximum relative error of the heat production is less than 2 %; the maximum relative error of the mass flow at the producer is less than 5 %. According to these results we can conclude that (i) the aggregation method achieves precise results at this level of aggregation, (ii) the dynamic optimization method can handle mass flow dependent delays in a very precise way and (iii) the

29

24

Consumer 1 Consumer 2

4

8

RI PT

0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.80

12 Time [hours]

16

20

24

SC

Absolute error [K]

ACCEPTED MANUSCRIPT

Figure 15: Absolute error between the result of the optimization and the simulation (when

M AN U

3 2 1 0 1 2 3 4 50

Prod m_flow

4

8

D

Relative error in %

Consumer 2 (orange).

12 Time [hours]

16

Prod Q_flow

20

1 2 3 4 24 5

TE

Figure 16: Relative error of the dynamic optimization result compared to the simulation of the original network when applying the results of the optimization: mass flow ad the producer

EP

(black dashed, left y-axis) and heat production (orange, right y-axis).

integration of physically constraints based on operational limitations is possible. The remaining NLP contains about 41.000 variables and 40.000 constraints.

AC C

For comparison, if there are 5 consumers left, the problem has about 90.000 variables and 89.000 constraints. This range of problems can be considered as small- to medium-sized for IPOPT [66]. The calculation time for the whole chain optimizing one day is about 3 minutes.

30

3 2 1 0

Rlative error in %

applying the results of the optmization): Temperature at Consumer 1 (balck dashed) and

ACCEPTED MANUSCRIPT

RI PT

6. Conclusion This paper presents a framework to represent, aggregate, dynamic thermo-

hydraulic simulate and optimize district heating systems based on a high-level equation-based modelling language and a high-level, large-scale dynamic op-

SC

timization method. The high-level description format meets the engineering need. The main advantage of the proposed method to decompose the MINLP

compared to evolutionary algorithms and approaches that solve the MINLP di-

M AN U

rectly is that the solution time is sufficiently low for real system applications. The main drawback of approaches that utilizes reduced models compared to the proposed approach is that they are difficult to use since the models are difficult to interpret physically. The proposed method combines the benefit of the good scale-up properties of the MIQCP method to find the discrete variables and the benefit from NLP to find the best economic dispatch that is highly dependent on physical variables. The Modelica language supports fast proto-

D

typing of physical models in different domains by the use of object-oriented

TE

constructs to enable the reuse of modelling knowledge. A growing community within academia and industry works on Modelica-based simulation of different scales and domains of energy systems. So far, few efforts have been made in the

EP

area of Modelica-based dynamic optimization of energy systems. Compared to most state-of-the-art tools which are based on imperative and domain specific methods, our framework offers a lot of advantages, including flexibility, modular

AC C

expandability, reusability of models and applicability to multiphysics problems. Standard methods often are restricted in their application as they rely on simplified models, static relationships and single-domain approaches and thus are unsuitable to investigate many issues of future 4GDH systems. Two use cases show the applicability of the framework. The first use case

demonstrates the applicability of Modelica as a suitable modelling language to describe and efficiently simulate large-scale district heating systems. The second use case shows the continuous optimization of a district heating system in a virtual city district. The integration of model coherences based on physical laws

31

ACCEPTED MANUSCRIPT

RI PT

into the optimization formulation enables the possibility to impose constraints on physically and technically relevant variables such as temperature, pressure, mass flow or costs, based on physical and operational limitations of the real

system. This improves the reliability of the result significantly. The aggregation method achieves precise results at an aggregation depth of 87.5 % and the

SC

remaining NLP can be considered small to medium sized. These results indicate that the presented method for continuous optimization is suitable for district heating systems with up to 150-200 consumers. The objective formulation can

M AN U

be easily extended to include economical or ecological factors. The solution time is sufficiently low for model predictive control of real-time applications. In our future research we intend to validate the proposed method to decompose the mixed-integer-optimal control problem, perform scale-up studies and integrate models of other domains (e.g. power system) into the framework.

D

7. Acknowledgement

TE

The research hase been supported by Modelon AB and the project FlexEnergySys funded by the Austrian Federal Ministry of Science, Research and Economics.

EP

[1] H. Lund, S. Werner, R. Wiltshire, S. Svendsen, J. E. Thorsen, F. Hvelplund, and B. V. Mathiesen, “4th Generation District Heating (4GDH): Integrat-

AC C

ing smart thermal grids into future sustainable energy systems,” Energy, vol. 68, pp. 1–11, 2014.

[2] G. Schweiger, J. Rantzer, K. Ericsson, and P. Lauenburg, “The potential of power-to-heat in swedish district heating systems,” Energy, 2017.

[3] B. V. Mathiesen, H. Lund, D. Connolly, H. Wenzel, P. A. Østergaard, B. M¨ oller, S. Nielsen, I. Ridjan, P. Karnøe, K. Sperling, and F. K. Hvelplund, “Smart Energy Systems for coherent 100% renewable energy and transport solutions,” Applied Energy, vol. 145, pp. 139–154, 2015.

32

ACCEPTED MANUSCRIPT

RI PT

[4] D. Connolly, H. Lund, B. V. Mathiesen, S. Werner, B. M¨oller, U. Persson, T. Boermans, D. Trier, P. A. Østergaard, and S. Nielsen, “Heat Roadmap Europe: Combining district heating with heat savings to decarbonise the EU energy system,” Energy Policy, vol. 65, pp. 475–489, 2014.

SC

[5] J. Allegrini, K. Orehounig, G. Mavromatidis, F. Ruesch, V. Dorer, and R. Evins, “A review of modelling approaches and tools for the simulation of district-scale energy systems,” Renewable and Sustainable Energy Reviews,

M AN U

vol. 52, pp. 1391–1404, 2015.

[6] T. Blochwitz, M. Otter, M. Arnold, C. Bausch, C. Clauß, H. Elmqvist, A. Junghanns, J. Mauss, M. Monteiro, T. Neidhold, D. Neumerkel, H. Olsson, J. V. Peetz, and S. Wolf, “The Functional Mockup Interface for Tool independent Exchange of Simulation Models,” in 8th International Modelica Conference 2011, pp. 173–184, 2009.

D

[7] J. Kofr´ anek, M. Matej´ ak, P. Privitzer, and M. Tribula, “Causal or acausal

TE

modelling: labour for humans or labour for machines,” in Technical Computing Prague, 2008.

[8] P. Fritzson, Principles of Object-Oriented Modeling and Simulation with

EP

Modelica 3.3. 2nd edition. Piscataway, NJ: Wiley-IEEE Press, 2015.

[9] M. Wetter, M. Bonvini, and T. S. Nouidui, “Equation-based languages A

AC C

new paradigm for building energy modeling , simulation and optimization,” Energy & Buildings, vol. 117, pp. 290–300, 2016.

R - A Unified Object-Oriented Language [10] Modelica Association, “Modelica

for Systems Modeling Language Specification Version 3.3 Revision 1,” 2014.

[11] P. Sahlin and E. F. Sowell, “A Neutral Format for Building Simulation Models Next Generation,” in Second International IBPSA Conference, pp. 147–154, 1989.

33

ACCEPTED MANUSCRIPT

RI PT

[12] S. A. Klein, “Development and integration of an equation-solving program for engineering thermodynamics courses,” Computer Applications in Engineering Education, vol. 1, no. 3, pp. 265–275, 1993.

[13] J. ˚ Akesson, K. ˚ Arz´en, M. G¨afvert, T. Bergdahl, and H. Tummescheit,

“Modeling and optimization with Optimica and JModelica.org Languages

SC

and tools for solving large-scale dynamic optimization problems,” Computers and Chemical Engineering, vol. 34, pp. 1737–1749, 2010.

M AN U

[14] M. Wetter and J. Wright, “A comparison of deterministic and probabilistic optimization algorithms for nonsmooth simulation-based optimization,” Building and Environment, vol. 39, pp. 989–999, 2004. [15] T. Mertz, S. Serra, A. Henon, and J.-M. Reneaume, “A MINLP optimization of the configuration and the design of a district heating network: Academic study cases,” Energy, vol. 117, Part, pp. 450–464, 2016.

D

[16] M. Sameti and F. Haghighat, “Optimization approaches in district heating

2017.

TE

and cooling thermal network,” Energy & Buildings, vol. 140, pp. 121–130,

[17] H. Lund, Renewable Energy Systems: The Choice and Modeling of 100%

EP

Renewable Solutions. Boston: Academic Press, second edition ed., 2014.

[18] D. Olsthoorn, F. Haghighat, and P. A. Mirzaei, “Integration of storage and

AC C

renewable energy into district heating systems : A review of modelling and optimization,” Solar Energy, vol. 136, pp. 49–64, 2016.

[19] F. F. M. Soons, J. I. Torrens, J. L. M. Hensen, and R. A. M. D. Schrevel, “A Modelica based computational model for evaluating a renewable district heating system,” in 9th International Conference on System Simulation in Buildings, pp. 1–16, 2014. [20] D. Strunik and J. Avsec, “Artificial neural networking model of energy and exergy district heating mony flows,” Energy and Buildings, vol. 86, pp. 366 – 375, 2015. 34

ACCEPTED MANUSCRIPT

RI PT

[21] J. Duquette, A. Rowe, and P. Wild, “Thermal performance of a steady state physical pipe model for simulating district heating grids with variable flow,” Applied Energy, vol. 178, pp. 383–393, 2016.

[22] P. Kr¨ auchi and M. Kolb, “Simulation thermischer Arealvernetzungen mit

SC

IDA-ICE,” in Proceedings of Bausim 2012, pp. 205–211, 2012.

[23] Y. Li, Y. Rezgui, and H. Zhu, “Dynamic simulation of heat losses in a district heating system: A case study in Wales,” EEE Smart Energy Grid

M AN U

Engineering (SEGE), pp. 273–277, 2016.

[24] S. Mohammadi, C. Bojsen, and C. Milan, “Identifying the optimal supply temperature in district heating networks - A modelling approach,” in 14th International symposium on district heating and cooling, pp. 1–6, 2014. [25] I. Gabrielaitiene, B. Bøhm, and B. Sunden, “Modelling temperature dy-

D

namics of a district heating system in Naestved , Denmark A case study,” Energy Conversion and Management, vol. 48, pp. 78–86, 2007.

TE

[26] M. Wetter, M. Fuchs, P. Grozman, L. Helsen, F. Jorissen, M. Lauster, M. Dirk, C. Nytsch-geusen, D. Picard, P. Sahlin, and M. Thorade, “IEA

EP

EBC ANNEX 60 Modelica Library - An International Collaboration to Develop a Free Open-Source Model Library for Buildings and Community Energy Systems,” Proceedings of BS2015, pp. 395–402, 2015.

AC C

[27] L. Giraud, R. Bavi`ere, C. Paulus, M. Vall´ee, and J.-f. Robin, “Dynamic Modelling , Experimental Validation and Simulation of a Virtual District Heating Network Dynamic Modelling , Experimental Validation and Simulation of a Virtual District Heating Network,” in ECOS 2015, 2015.

[28] L. Andresen, P. Dubucq, R. Peniche, G. Ackermann, A. Kather, and G. Schmitz, “Status of the TransiEnt Library : Transient simulation of coupled energy networks with high share of renewable energy,” in 11th International Modelica Conference, pp. 695–705, 2015.

35

ACCEPTED MANUSCRIPT

RI PT

[29] E. Dotzauer, Algorithms for Short-Term Production- Planning of Cogeneration Plants. PhD thesis, 1997.

[30] E. Dotzauer, “Produktionsplanering av el och v¨arme- Matematiska modeller och metoder,” tech. rep., 2002.

SC

[31] K. M. Powell, J. Suk, W. J. Cole, K. Kapoor, J. L. Mojica, J. D. Hedengren, and T. F. Edgar, “Thermal energy storage to minimize cost and improve efficiency of a polygeneration district energy system in a real-time electricity

M AN U

market,” Energy, vol. 113, pp. 52–63, 2016.

[32] M. Vesterlund, A. Toffolo, and J. Dahl, “Optimization of multi-source complex district heating network, a case study,” Energy, 2017. [33] B. Morvaj, R. Evins, and J. Carmeliet, “Optimising urban energy systems: Simultaneous system sizing, operation and district heating network layout,” Energy, vol. 116, Part 1, pp. 619 – 636, 2016.

D

[34] S. Grosswindhager, A. Voigt, M. Kozek, and A. V.-c. Models, “Predictive

TE

Control of District Heating Network using Fuzzy DMC,” in International Conference on Modelling, Identification and Control, 2012. [35] H. V. Larsen, B. Bøhm, and M. Wigbels, “A comparison of aggregated

EP

models for simulation and operational optimisation of district heating networks,” Energy Conversion and Management, vol. 45, pp. 1119–1139, 2004.

AC C

[36] G. Sandou, G. Sandou, S. Font, S. Tebbani, A. Hiret, and C. Mondon, “Predictive Control of a Complex District Heating Network,” in IEEE Conference on Decision and Control, pp. 7372–7377, 2006.

[37] M. Vesterlund and J. Dahl, “A method for the simulation and optimization of district heating systems with meshed networks,” Energy Conversion and Management, vol. 89, pp. 555–567, 2015. [38] K. Orehounig, R. Evins, and V. Dorer, “Integration of decentralized energy systems in neighbourhoods using the energy hub approach,” Applied Energy, vol. 154, pp. 277–289, 2015. 36

ACCEPTED MANUSCRIPT

energy carriers,” in 15th PSCC, pp. 22–26, 2005.

RI PT

[39] M. Geidl, “Optimal power dispatch and conversion in systems with multiple

[40] X. S. Jiang, Z. X. Jing, Y. Z. Li, Q. H. Wu, and W. H. Tang, “Modelling

and operation optimization of an integrated energy based direct district

SC

water-heating system,” Energy, vol. 64, 2014.

[41] E. Guelpa, C. Toro, A. Sciacovelli, R. Melli, E. Sciubba, and V. Verda, “Optimal operation of large district heating networks through fast fluid-

M AN U

dynamic simulation,” Energy, vol. 102, pp. 586–595, 2016.

[42] D. H¨ aggst˚ ahl and E. Dotzauer, “Production planning with respect to uncertainties - Simulator based production planning of average sized combined heat and power production plants,” tech. rep., 2004. [43] A. Viana, J. P. D. Sousa, and M. Matos, “A new metaheuristic approach to the unit commitment problem,” in 14th Power Systems Computation

D

Conference, no. June, pp. 24–28, 2002.

TE

[44] K. Steer, A. Wirth, and S. Halgamuge, “Control period selection for improved operating performance in district heating networks,” Energy and

EP

Buildings, vol. 43, no. 23, pp. 605 – 613, 2011. [45] K. Dietl, S. G. Yances, and A. Johnsson, “Industrial application of optimization with Modelica and Optimica using intelligent Python scripting,”

AC C

in 11th International ModelicaConference, pp. 777–786, 2014.

[46] H. Runvik, P.-O. Larsson, S. Velut, J. Funquist, M. Bohlin, A. Nilsson, and S. Modarrez Razavi, “Production Planning for Distributed District Heating Networks with JModelica.org,” in 11th International Modelica Conference,

pp. 217–223, 2015.

[47] M. K¨ ofinger, D. Basciotti, R. R. Schmidt, E. Meissner, C. Doczekal, and A. Giovannini, “Low temperature district heating in Austria : Energetic , ecologic and economic comparison of four case studies,” Energy, vol. 110, pp. 95–104, 2016. 37

ACCEPTED MANUSCRIPT

Processes, vol. 3, pp. 471–496, 2015. [49] Modelon AB, “OPTIMICA Compiler Toolkit.”

RI PT

[48] F. Magnusson and J. ˚ Akesson, “Dynamic Optimization in JModelica.org,”

[50] J. Andersson, A General-Purpose SoftwareFramework for DynamicOpti-

SC

mization. Phd thesis, 2013.

[51] A. W¨ achter and L. T. Biegler, “On the implementation of an interior-

M AN U

point filter line-search algorithm for large-scale nonlinear programming,” Mathematical Programming, vol. 106, no. 1, pp. 25–57, 2006. [52] H. D. Mittelmann, “AMPL-NLP Benchmark,” 2016. [53] Pyomo, “Online Documentation 4.1,” 2016. [54] Gurobi, “Gurobi Documentation,” 2016.

D

[55] S. Velut, P.-O. Larsson, J. Windahl, L. Saarinen, and K. Boman, “Shortterm production planning for district heating networks with JModel-

TE

ica.org,” in 10th International ModelicaConference, pp. 959–968, 2014. [56] K. M. Powell, A. N. Eaton, J. D. Hedengren, and T. F. Edgar, “A contin-

EP

uous formulation for logical decisions in differential algebraic systems using mathematical programs with complementarity constraints,” Processes,

AC C

vol. 4, no. 1, 2016.

[57] M. Fouquet, H. Guguen, D. Faille, and D. Dumur, “Hybrid dynamic optimization of power plants using sum-up rounding and adaptive mesh refinement,” in 2014 IEEE Conference on Control Applications, pp. 316–321, Oct 2014.

[58] J. Rantzer, Robust production planning. Ma thesis, 2015. [59] J. Arroyo and A. Conejo, “Modeling of start-up and shut-down power trajectories of thermal units,” Power Systems, IEEE Transactions on, vol. 19, pp. 1562–1568, Aug 2004. 38

ACCEPTED MANUSCRIPT

RI PT

[60] L. T. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes. Philadelphia, PA: MOS-SIAM, 2010.

[61] L. Brange, J. Englund, and P. Lauenburg, “Prosumers in district heating networks - A Swedish case study,” Applied Energy, vol. 164, pp. 492–500,

SC

2016.

[62] A. Loewen, Entwicklung eines Verfahrens zur Aggregation komplexer

M AN U

Fernw¨ armenetze. PhD thesis, 2001.

[63] A. A. Hagberg, D. A. Schult, and P. J. Swart, “Exploring network structure, dynamics, and function using NetworkX,” Proceedings of the 7th Python in Science Conference, no. SciPy, pp. 11–15, 2008.

[64] R. H¨ agg, Dynamic Simulation of District Heating Networks in Dymola Dynamic Simulation of District Heating Networks in Dymola. Ma thesis, 2016.

D

[65] S. Velut and H. Tummescheit, “Implementation of a transmission line model for fast simulation of fluid flow dynamics,” in 8th International Mod-

TE

elica Conference, p. 8, 2011. [66] P. Larsson, Optimization of low-level controllers and high-level polymer

AC C

EP

grade changes. PhD thesis, 2011.

39