Two-phase optimizing approach to design assessments of long distance heat transportation for CHP systems

Applied Energy 182 (2016) 164–176

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Two-phase optimizing approach to design assessments of long distance heat transportation for CHP systems Piotr Hirsch ⇑, Kazimierz Duzinkiewicz, Michał Grochowski, Robert Piotrowski ´ sk University of Technology, Gdan ´ sk, Poland Faculty of Electrical and Control Engineering, Gdan

h i g h l i g h t s
New method for long distance heat transportation system effectivity evaluation.
Decision model formulation which reflects time and spatial structure of the problem.
Multi-criteria and complex approach to solving the decision-making problem.
Solver based on simulation-optimization approach with two-phase optimization method.
Sensitivity analysis of the optimization procedure elements.

a r t i c l e

i n f o

Article history: Received 1 March 2016 Received in revised form 1 August 2016 Accepted 17 August 2016

Keywords: Pipeline optimization Heat transportation Decision-making problem Life cycle cost modeling

a b s t r a c t Cogeneration or Combined Heat and Power (CHP) for power plants is a method of putting to use waste heat which would be otherwise released to the environment. This allows the increase in thermodynamic efficiency of the plant and can be a source of environmental friendly heat for District Heating (DH). In the paper CHP for Nuclear Power Plant (NPP) is analyzed with the focus on heat transportation. A method for effectivity and feasibility evaluation of the long distance, high power Heat Transportation System (HTS) between the NPP and the DH network is proposed. As a part of the method the multi-criteria decisionmaking problem, having the structure of the mathematical programming problem, for optimized selection of design and operating parameters of the HTS is formulated. The constraints for this problem include a static model of HTS, that allows considerations of system lifetime, time variability and spatial topology. Thereby variation of annual heat demand within the DH area, variability of ground temperature, insulation and pipe aging and/or terrain elevation profile can be taken into account in the decision-making process. The HTS construction costs, pumping power, and heat losses are considered as objective functions. In general, the analyzed optimization problem is multi-criteria, hybrid and nonlinear. The two-phase optimization based on optimization-simulation framework is proposed to solve the decision-making problem. The solver introduces a number of assumptions concerning the optimization process. Methods for problem decomposition, scalarization and relaxation are proposed and optimization procedures for the decomposed problem are discussed. The methodology is tested on a sample case study of the NPP planned to be built in Northern Poland. The sensitivity analysis of the problem is also provided. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction The European energy policy, concluded in the Strategic Energy Technology (SET) [1] Plan, aims at the reduction of greenhouse gas emissions (low-carbon society) and enhancing energy supply security. One of low carbon energy technologies supported in the

⇑ Corresponding author. E-mail addresses: [email protected] (P. Hirsch), kazimierz.duzinkiewicz@ pg.edu.pl (K. Duzinkiewicz), [email protected] (M. Grochowski), robert. [email protected] (R. Piotrowski). http://dx.doi.org/10.1016/j.apenergy.2016.08.107 0306-2619/Ó 2016 Elsevier Ltd. All rights reserved.

SET Plan is nuclear energy. Nowadays it provides nearly 30% of the electricity consumed in the EU [2]. Nuclear energy is also used in small scale for district and process heating purposes, mostly in cogeneration applications. Benefits coming from nuclear cogeneration (simultaneous heat and power production in one plant) include: increased thermal efficiency of the Combined Heat and Power (CHP), Nuclear Power Plant (NPP), and delivery of carbonfree heat for residential and industrial heating. Nuclear reactors are basically heat generators and are technically viable to be employed as heat sources. The main obstacle in nuclear cogeneration development is effective heat transportation.

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

Today, only about 33% of fission energy produced in a standard NPP is transformed into electricity [3]. The remaining part of energy is dumped into the environment as waste heat. One of reasons why nuclear cogeneration is not used on a larger scale, despite the advantages it provides, is, in general, long distances between the NPP site and crowded urban areas, which make providing sufficient market for heat delivery impossible, or extremely difficult. Heat transportation over long distances, exceeding typical district heating cases, is often considered as ineffective because of high thermal losses. For example, in the recent study on roles of the innovative technologies in the zero emission scenario, nuclear energy was only considered for electricity production [4]. However, nowadays conditions for such projects have significantly changed. Modern insulation materials have thermal conductivity at such a low value that it is possible to transport the heat over very long distances e.g. 150 km, with only a few percent of the transported energy being lost [3]. Given the technological outbreak, it is reasonable to reassess the prospect of supplying nuclear and conventional heat to District Heating (DH) networks or industry over long distances. Due to high complexity of the discussed matter, functional decomposition of NPP heat extraction and transportation is proposed. The issue related to heat extraction includes construction and operation of a subsystem intended to extract the heat from the NPP. Heat may be extracted either by recovering the waste heat in a condenser system, or by utilizing the steam bleed from the turbine system. To provide enough power to the condenser systems, the turbine systems have to be redesigned. Therefore the latter option seems to be more feasible, as it does not require any changes in the existing NPP turbine systems. The other part of the decomposed problem relates to the construction and operation of a long distance Heat Transportation System (HTS), which is the main topic of this paper. Heat transportation is an important issue not only for NPP CHP, but also for conventional CHP or similar applications related to district or technological heating. The approach presented in the paper applies to an arbitrary HTS, regardless of the heat source and destination.

1.1. Survey of related works Recent study [5] shows, that conversion of conventional thermal power plants into CHP plants is highly attractive from economic and environmental perspectives. However, it was assumed that power plants were located at ‘manageable distance’ to nearest town. Regarding the nuclear cogeneration, preliminary economic _ studies have shown that the CHP NPP located by Lake Zarnowiec, Northern Poland, would have potential to become a primary source of heat in the surrounding area [6–8]. More detailed, but still rough economic estimations were performed for the French Nogent-surSeine NPP located 110 km east of Paris [3]. The payback time was expected to be less than 10 years, with the gain equal to 360 M€/ year. The optimal piping diameter for a 1500 MW heat power transportation system has been found equal to 2000 mm, as the result of thermal energy cost vs. pipeline diameter minimization. In [9] a 77 km long 1000 MW HTS has been analyzed for possible CHP production at the Loviisa 3 NPP project in Finland. The simulations presented by the authors focused on analyzing the behavior and safety issues of the discussed large-scale HTS. Methods for optimal design of long distance fuel (oil/gas) pipeline systems are already well known. A survey of these methods is presented in [10]. Three typical optimization problems have been identified, which are: the design of the distribution pipeline network, and its operational optimization and batch scheduling. The first two problems are to some extent similar to the optimal HTS design task. It has been stated that modern population-based opti-

165

mization techniques, such as evolutionary and genetic algorithms, reveal promising performance in dealing with problems of this type. Sample utilization of such techniques is shown in [11], where the evolutionary algorithm is used to solve the problem of minimizing the consumption of electricity for pumping, the number of pumping stations, and the pipeline mass required for an oil pipeline project. The decision variables specified for this optimization problem were: the outer diameter and wall thickness of the pipeline, and the suction and discharge pressures of the pump stations. Another method for optimal sizing of a pipeline for transporting multisized solid-liquid mixtures is presented in [12]. Time dependences of various costs over the life cycle, with the pipeline diameter treated as a decision variable and objective function in the form of the summed construction, and the pumping energy costs have been studied. The solution of the given optimization problem was obtained by equating the objective function first derivative to zero. In [13] an overview of natural gas transportation via pipeline systems is presented. The gas transmission problem involves designing and operation of the system. The former issue refers to decisions on pipeline length and diameter, and locations of compressor stations, while the latter applies to the compressors operating costs. Both, the steady-state and transient state approaches are reviewed. In [14] a design of the CO2 transportation pipeline network is discussed. Authors find the optimal pipeline diameter by comparing capital and operational costs of several design options. Analyzing methods for DH construction and operation may be helpful for determining the HTS design method specification. In [15] the importance of DH optimization is stressed. The oversized and badly insulated pipeline system in Narva city (Estonia) is a source of significant heat losses, which could be decreased up to three times if the DH network was optimized. The analysis of four types of pre-insulated pipeline systems for small-scale DH retrofit is presented in [16]. The comparison of the scenarios took into account pressure losses, heat losses, and investment costs. Steady-state heat losses in buried pipes have been investigated in [17]. The influence of soil temperature throughout the year and the temperature-dependence of the conductivity coefficient were taken into account to achieve high accuracy of the results obtained using the finite element method. In [18] the minimization of short-term operation costs via supply temperature optimization, done using the method based on dynamic programming, has been presented. The authors considered some major simplifications such as: fixed flow distribution, heat loss independence from the ground temperature, and fixed temperature drop across consumers. In [19] the annual cost of the DH piping network was established in the function of pressure drop per unit length and minimized. Design parameters of the network such as temperature regimes were treated as parameters and used to test different case studies. 1.2. Problem description and main contributions The overview of the related works shows, that CHP NPPs have the potential to become important heat suppliers for DH. Adequate methods and tools are needed to evaluate technical and economic feasibility of such systems. Design and operation of a long distance, high power HTS is a part of this issue. The first approach to the problem was presented in the conference paper [7], with the main focus on proposition of the mathematical model for the decisionmaking process. Upon the overview of methods developed for fuel transportation and DH networks, a desired characteristic of the HTS design assessment method can be formulated. First of all, the proposed method should take advantage of optimization techniques to minimize HTS construction costs, including: costs of materials and

166

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

labor for both the pipeline and the pumping stations, together with the electric power consumed by pumping stations and the heat power lost on transportation. Inner pipe diameter, supply and return temperatures, insulation thickness and number of pumping stations are identified as the important design variables. The estimation of HTS cost effectiveness has to take into account the presumed system lifetime and the variability of the accompanying parameters, such as: pipeline roughness, thermal conductivity of insulation, and heat and electricity prices, all considered over the period. Other important aspects that need to be considered are: annual variability of DH heat power demand and ground temperature, as well as the impact of the terrain elevation profile on the number and location of pumping stations. The proposed methodology requires precise specification of a number of inputs. Firstly, the localizations of NPP and DH are assumed to be known. This allows to define the terrain elevation profile and the predicted annual ground temperature distribution. Secondly, the annual DH heat demand trajectory has to be specified in the form of structured graph. Further, the range of possible temperatures, pressures and flow velocities for the analyzed case is to be identified. The main objective of the paper was to propose a method for effectivity and feasibility evaluation of the long distance HTS used in the NPP CHP over its life cycle. The proposed method includes novel decision model formulation, which takes into account time and spatial structure of the project. This requires a multi-level problem discretization. The decision problem is solved by using the proposed two-phase optimization method based on simulation-optimization framework. The presented approach allows optimal selection of parameters of the long distance HTS, as well as economic and technical evaluation of the project. Sensitivity analysis of the elements of the optimization procedure against variability of parameters, e.g. electricity and heat prices or pipe roughness completes the case study.

2. Design assessments as optimizing decision problem In the paper the decision model for optimal HTS design is proposed. It consists of the mathematical model of objectives and constraints, and optimization problem formulation. The mathematical model allows evaluation of different HTS designs, with respect to objectives and constraints, for a given case. The designs are chosen from the decision space by the optimizer, with the aim of minimizing the specified objective functions. The mathematical model of the constraints consists of two parts: hydraulic and thermal. Transient states influence on the evaluation can be neglected because of the long evaluation period. Therefore, only steady states for hydraulic and thermal models are considered. When constructing

(a)

Q& [MWh ]

1

the mathematical model, the following assumptions were made: the minor head losses are neglected, the density of water and the friction factor are temperature independent, the heating medium is incompressible, and the pumping impact on the heating medium temperature is neglected. These assumptions allow to reduce the computational burden without significant loss of accuracy. The HTS life cycle expectancy is several decades. During that time some characteristics of the system are changing slowly, due to processes such as insulation and pipe aging, or energy market shifts. It is convenient to approximate the impact of these processes as fixed over some intervals. This approach leads to Long Horizon Time Discretization (LHTD) of the HTS life cycle LT, with DsT time step, into kt intervals indexed by t. As a result, for selected parameters, such as pipeline roughness, thermal conductivity, and the price of heat and electricity, it is possible to acquire different constant values for each life cycle interval. Depending on the decision-maker needs and knowledge, the LHTD time step can either be fixed or variable. Two important inputs to the HTS, namely the annual heat demand and the annual ground temperature, vary greatly throughout the year. The maximum heat demand is usually associated with the lowest ground temperature, which results in extreme values of flows, temperatures and losses, and translates into most demanding technological requirements for the given HTS design. It is noteworthy, however, that the maximum value of the heat demand generally occurs only during the short part of the year, while being much lower during the remaining period. Therefore, on the one hand, the optimally designed HTS needs to be capable of transporting extreme heat load, but, on the other hand, it should ensure that the system is not oversized. Short Horizon Time Discretization (SHTD) of the decision model according to the discretized, structured graph of annual heat demand (see Fig. 1) is proposed as the solution to the above problem. In this way, M sub-models with variable time step DsY indexed by y, are established for each of kt intervals of the LHTD life cycle. Each submodel corresponds to a different heat demand level (different part of the year), and is characterized by different flows and temperatures. During the computation, each sub-model operates on the same vector of decision variables. The results of computations are merged by weighing in accordance to the part of the year to which the given sub-model applies. The total duration of the single set of the M sub-models is one year. Although the duration of each LHTD interval may be decades, its parameters are assumed to be constant over that period and, consequently, it may be covered by only one M sub-model, the results of which are multiplied according to the LHTD time step. In the long distance HTS system considered in this study the terrain elevation profile (see Fig. 2) can have significant effect on

Q& [MWh ]

(b)

2

3

4

5

1

t [h]

2

3

4

5

t [h]

Fig. 1. Sample structured graphs of annual heat demand of district heating area, (a) long period of high and more balanced heat demand and (b) significant step changes in heat demand, as a result of heating season changes, y stands for sub-model number.

167

Terrain elevation profile [m]

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

Δl

Fig. 4. Pipeline section scheme with associated variables/parameters.

Distance [km]

Fig. 2. Sample HTS terrain profile with scheme of spatial discretization, dotted red line indicates approximation result. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the nodal pressure of the system and on the locations of pumping stations. It implies the need for spatial model discretization in order to efficiently deal with system topology. The considered HTS is regarded as a network composed of nodes and links between them. Its structure is shown in Fig. 3. The links are pipeline sections, pump stations, or heat exchangers. As the length of the pump station link and the heat exchanger link is assumed to be negligibly small, this kind of links is complemented by the pipeline sections. The total length L (km) of the HTS pipeline in one direction, is discretized with the step Dl (m) into k nodes, indexed by i. Two directions are denoted with u.

(W), and the heat power Qk,2 at the last node [k, 2] is equal to 0. Heat powers and temperatures at other nodes are calculated based on heat losses in the appropriate pipeline section at given y, t. A scheme of the pipeline section having the length Dl is shown in Fig. 4, where: U is the inner pipe diameter (m), s represents pipe walls thickness (m), and d is insulation thickness (m). Tg,k is the ground temperature (K). If the pipeline section length Dl is small enough, the heating medium temperature for this section can be approximated by the constant value equal to Ti1 (K), which is the node temperature at section inlet. The heat flux transferred from the pipeline to the environment (heat loss) DQi (W) over the i pipeline section is described by relation (1) [6,7]:

DQ i Q i1  Q i ¼ ¼ Dl Dl

2.1. Mathematical models The primary notation scheme used in the mathematical model description has the following form: i, u, y, t are the subscripts, where i = 1, 2, . . ., k stands for the node number, and u = 1, 2 stands for the supply or return line respectively. y = 1, 2, . . ., M is the submodel number and t = 1, 2, . . ., kt specifies the interval number. Each variable dependent on i is also dependent on u, y and t, e.g.: the pressure at each node for each sub-model in each interval Pi,u,y,t. For this reason the reduced notation is introduced that omits u, y, t for these variables, unless necessary, e.g. giving Pi. 2.1.1. Thermal model The thermal model used in the paper is characterized by heating medium mass flow, and heat power and temperature at each node. Based on the initial assumptions, at node [1,1] the temperature T1,1 is equal to the supply temperature Ts (K), while at node [k, 2] the temperature Tk,2 is equal to the return temperature TR (K). Consequently, the temperature drop obtained in the NPP heat exchanger is DTNPP = T1,1  Tk,2 (K). Similarly, the heat power Q1,1 at the first node is equal to the heat power absorbed from NPP Q NPPy;t

T i1  T gy 1 pUhcw

þ pU

1 in hcg

þ lnð2Upok=pUÞ þ lnð2Upink=Uo Þ

ð1Þ

ist

where Uin = U + 2s + 2d is the diameter of the pipe with insulation layer (m), Uo = U + 2s is the outer pipe diameter (m), hcw and hcg are convective heat transfer coefficients for water and ground (W/ (m2K)), respectively. kis and kp are thermal conductivity coefficients (W/(mK)) of insulation and pipe walls material, respectively. The heating medium temperature at node i can be calculated from the heat loss over the pipeline section i:

DT i ¼ T i1  T i ¼

DQ i _ y;t C p m

ð2Þ

In Eq. (2) Cp is specific water heat capacity (J/(kgK)) at 100 °C, _ y;t is the mass flow of the heating medium (kg/s). Assuming and m the absence of the heating medium phase change, the mass flow _ y;t can be calculated, once the heat power and temperature drop m on either DH or NPP heat exchanger are known:

_ y;t ¼ m

Q DHy

ge DT DHy;t C p

Fig. 3. Scheme of HTS structure.

¼

Q NPPy;t

ge DT NPP C p

ð3Þ

168

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

In Eq. (3) Q DHy is the heat delivered to the DH network (W). It is equal to the heat demand. DT DHy;t is the temperature drop on the DH exchanger (K) and ge is the heat exchanger efficiency. The amount of heat absorbed from NPP Q NPPy;t can be calculated as the sum of the demanded heat power Q DHy and the total heat power lost in transportation for the y sub-model DQ Ly;t (W):

Q NPPy;t ¼ Q DHy þ DQ Ly;t

ð4Þ

DQ Ly;t ¼ DQ LSy;t þ DQ LRy;t ¼ Q 1;1  Q k;1 þ Q 1;2  Q k;2

ð5Þ

In the above equations DQ LSy;t and DQ LRy;t are heat power losses (W) in the supply and return line, respectively. 2.1.2. Hydraulic model The heating medium hydraulics is characterized by two variables: flow velocity v (m/s) and pressure P (Pa). The analyzed HTS pipeline has two important features: fixed diameter and the structure that does not include branching nodes (Fig. 1). As a consequence, the flow velocity is assumed constant for all nodes at a given time step. The flow velocity can be obtained from the mass flow:

v y;t ¼

_ y;t _ y;t 4m m ¼ qA qp  U2

ð6Þ

where q is the density of water (kg/m3) and A is the cross-sectional area (m2). The pressure P1,1 at the supply line inlet was assumed to be equal to the nominal pressure pn (Pa). Because of the neglected minor pressure losses, P1,2 is equal to Pk,1. Nodal pressures are calculated forward according to the following formula [20]:

Pi1 Pi þ zi1 þ hg ¼ þ z þ hly;t qg qg i

ð7Þ

where qPg is the pressure head (m), z is the elevation head (m) and hg, hl represent either head gain or loss at corresponding nodes (m). The hydraulic scheme of the pipeline section is shown in Fig 5. The elevation head can be determined from the terrain elevation profile. The head losses over the pipeline section i are described by the Darcy-Weisbach equation:

ð8Þ 2

where g is the acceleration due to gravity (m/s ) and f is the DarcyWeisbach friction coefficient, which for fully turbulent flow (Reynolds number > 4000) can be calculated from Eq. (9):

ð9Þ

In Eq. (9) n is inner surface roughness of the pipe (mm). If the pressure at a given node is expected to drop below the minimal values pmin,S or pmin,R (MPa), which are equal to the saturation pressures at the corresponding lines, increased by a safe margin, a pump station is added at that location. The head gain vs. flow rate pump characteristic is assumed to be linear, thus the head added by the pump station is constant regardless of the flow rate. Finally, the pressure at node i is calculated from:

Pi ¼ P i1  qg  ðzi  zi1 þ hg  hly;t Þ

ð10Þ

2.2. Optimization problem formulation Each sub-model introduced by SHTD is characterized by the mass flow, head loss, heat absorbed from NPP, number of pumping stations, main pump pressure gain, heat power at each node, heating medium temperature at each node, and pressure at each node. All of these variables, together with certain parameters, such as the inner pipe diameter and insulation thickness, constitute the decision space xd. However, not all of these variables are of equal importance for the decision making process. Based on HTS characterization, three categories of variables can be identified. The first category includes variables crucial to HTS design, in terms of both: decision maker preferences and technological. They are collected into a so-called design vector, denoted xd1 . The second category, composing the simulation vector xd2 , encloses variables of operational nature. These variables are related to the design vector via the mathematical model. The last category is formed to include auxiliary variables, implemented to deal with logical conditions and the mixed-integer nature of the system. This category is called auxiliary and denoted xd3 . According to the above categorization, the decision space is represented as:

 xd ¼ xd1

xd2

xd3

T

ð11Þ

Due to technological limitations, the following bounds on xd can be considered:

xdlb 6 xd 6 xdub

Dlv y;t 2g  U 2

hly;t ¼ f

  1 n pffiffiffi ¼ 2log10 3:7U f

ð12Þ

where xdlb and xdub are the lower and upper bounds, respectively. The decision maker can also specify additional inequality constraints associated with the operation of the HTS. For instance, the maximum heating medium flow velocity can be specified to prevent pipe damage due to vibrations and erosion. The maximum allowed pressure drop and heat losses over pipe section length Dl are the other examples. A common feature of this type of constraints is that they apply to the sub-model related with the peak heat demand of the DH. The remaining sub-models are characterized by lower flow velocities, and lower heat and pressure losses (maximum heat demand is associated with minimal temperatures), therefore in those cases the inequality constraints are consequently met. The introduced inequality constraints have the following vector form:

8 g 1 ðxd Þ > > > > < g 2 ðxd Þ gðxd Þ ¼ .. > > . > > : g Nic ðxd Þ

ð13Þ

where g ic ðxd Þjic¼1;2;...;Nic are nonlinear inequality constraints, with Fig. 5. Scheme of pipeline section hydraulics [20].

their number equal to Nic. The constraints on maximum heating medium flow velocity vmax (m/s), maximum pressure drop DP/Dlmax

169

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

(Pa/km) and maximum heat losses DQ/Dlmax (W/km), in relation to the mathematical model Eqs. (1)–(10), take the following form:

gðxd Þ ¼

8 > > > > > > < > > > > > > :

4Q DH

6 v max

1;kt

pqC p DT DH1;k U2 t

f qv

2 1;kt

2U T S T g1

6 DP=Dlmax

lnðU =UÞ

ð14Þ

lnðUin =Uo Þ 6 DQ =Dlmax

o 1 1 pUhw þpUin hg þ 2pkp þ 2pkis

The second objective of HTS design is pumping power reduction. Like for the heat loss reduction, the total consumed electric power WTL (MW h) is estimated using the annual pumping power

 h , which is the outcome of the summation of results from W At MW year the decomposed models (Eq. (19)). WTL is calculated using Eq. (20).

W At ¼

y¼1

kt

Relations between the simulation and design variables are described by the mathematical model and can be represented by the set of nonlinear equality constraints:

 hðxd Þ ¼ hj ðxd Þj¼1;...;kt M

ð15Þ

where hj(xd) is the set of constraints imposed by the single submodel y, t (the right side indices y, t are omitted):

8 _  c4  Q NPP ¼ 0 m > > > > > _2¼0 P l  c2  m > > > > > Q  Q DH þ ðQ 1;1  Q k;1 Þ þ ðQ 1;2  Q k;2 Þ ¼ 0 > > > NPP > > Q 1;1  Q NPP ¼ 0 > > > > > Q  Q > 2:k;1 1:k1;1  c 1 ðT 1:k;1  T g Þ ¼ 0 > > > > > Q 1:k1;2  Q 2:k;2 þ c1 ðT 1:k;2  T g Þ ¼ 0 > > > > > Q k;2 ¼ 0 > > <  TS ¼ 0 T 1;1 hj ðxd Þ ¼ > ðQ Q 2:k;1 Þ > > T 2:k;1  T 1:k1;1  1:k1;1 ¼0 > _ C p m > > > > ðQ 1:k1;1 Q 2:k;1 Þ > T 1:k1;2  T 2:k;2 þ ¼0 > _ > C p m > > > > T  T ¼ 0 > R k;2 > > > > P1;1  Pn ¼ 0 > > > > > > P2:k;1  P1:k1;1  c3 ðz2:k;1  z1:k1;1 Þ  Pl ¼ 0 > > > > > P 1;2  P k;1 ¼ 0 > > : P2:k;2  P1:k1;2  c3 ðz2:k;2  z1:k1;2 Þ  Pl ¼ 0

" M X kw

W TL ¼

gp

#

 DsY y ðny;t  pg þ pmgy;t Þ  V_ y;t

kt X ½DsT t  W At 

ð19Þ

ð20Þ

t¼1

where kw = 24 is the unit conversion coefficient, gp is a pump efficiency, ny,t is the number of pumping stations set for the given sub-model, pg is the pressure added by the pump station (MPa), p is the pressure added by the main pump station (Pa) and V_ y;t mg y;t

is the volumetric flow rate (m3/s). In addition to the aforementioned heat losses and electric power consumption, the HTS construction cost are evaluated:

C C ¼ 2L  f d ðUÞ þ cps  n1;kt

ð16Þ

ð21Þ

where Cc is the overall HTS construction cost (m.u.), fd(U) is the pipeline material and labor cost, as a function of pipeline diameter, and cps is the cost of pumping station construction (m.u.). With the objective functions Eqs. (17)–(21) the objective vector is given as:

Jðxd Þ ¼ ½ C C ðxd Þ DQ TL ðxd Þ W TL ðxd Þ T

ð22Þ

Hence, the general optimization problem can be formulated as follows (23)–(26):

min Jðxd Þ

ð23Þ

s:t: hðxd Þ ¼ 0

ð24Þ

where expressions c1,2,3,4 are introduced for conciseness and are defined as follows:

gðxd Þ 6 0

ð25Þ

2pkis Dl logð1 þ 2d UÞ    1 n 8Dl 2 c2 ¼ 4log10 3:7U qp2 U5

xdlb 6 xd 6 xdub

ð26Þ

c3 ¼ qg

Proposition of the solver to the decision problem described in previous chapter was preceded by problem analysis. The presented optimization problem Eqs. (23)–(26) are nonlinear, constrained, hybrid and multi-criteria, what makes it generally troublesome to solve. However, taking into account the nature of the considered decision model, certain features of the optimization problem can be utilized for problem simplification. In order to reduce the size of the design vector, and to effectively solve the decision problem, a simulation-optimization framework is proposed (Fig. 6). It is convenient to decompose the decision problem into several cooperating structures, which can be called optimizer, simulator and pumping stations module. Additionally, the starting points generation module for simulator can be extracted. Each structure operates on different part of the decision space, while the other part is parameterized. The simulator solves d1 the set of nonlinear equality constraints (Eq. (24)) for given x

xd

c1 ¼

3. Simulation-optimization framework methodology

1 c4 ¼ ge C p ðT S  T R Þ Heat loss reduction is one of the objectives of the proposed methodology for optimal HTS design. The total amount of the heat lost throughout the HTS lifetime, DQTL (GJ), is estimated using the total heat power lost on transportation for y sub-model DQ Ly;t .

 GJ is calculated as the sum of Firstly, the annual heat loss DQ At year losses generated by each decomposed model (Eq. (17)). Then the DQTL is established from Eq. (18).

DQ At ¼

M h X

kQ  DsY y  DQ Ly;t

i

ð17Þ

y¼1

DQ TL ¼

kt X 

Ds T t  DQ A t



ð18Þ

t¼1

where: kQ = 8.64  105 is the unit conversion coefficient, DsT t is the time step of interval t (years) and DsY y is the time step of sub-model y (days).

(dash means that the variable is treated as a parameter), returning xd2 as the result. The optimizer finds xd1 which satisfies the inequality constraints and bounds (Eqs. (25) and (26)), in order d2 and to minimize the objective vector Eq. (23), having given x d3 . The pumping stations module helps in dealing with the x mixed-integer property of the system. Thereby, the simulator is

170

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

STARTING POINTS GENERATION MODULE

STARTING POINTS

OPTIMIZER min

PUMPING STATIONS MODULE Eq. (45)

SIMULATOR )

)≤0

Fig. 6. Framework scheme of the solver.

the system of nonlinear equality equations, while the optimizer operates only on xd1 . The inner pipe diameter, the return and supply temperatures, the insulation thickness and the number of pumping stations were specified as the decision variables. However, the supply and return temperatures belong to the set of real numbers, while the remaining quantities are integers. The mixed-integer property of the problem can be treated using additional assumptions and approximations which will leave only real decision variables, starting with the inner pipe diameter, which in general belongs to the set of integers but may be treated as a real number due to the relaxation method. The real number solution is then rounded to the nearest value in the set of available pipes diameters. The insulation thickness represents another case, as there are only few insulation options available (most often two standardized thicknesses). Hence, it is effective to treat the insulation thickness as a parameter (fix its value) and to perform the optimization for each case separately, comparing the results at the final step. Regarding the number of pumping stations, it can be effectively solved with introduction of the pumping stations module to the optimizersimulator structure (Fig. 6). This enables transformation of the number of pumping stations and the main pump pressure gain into the vector xd3 of auxiliary variables. These considerations allow to simplify the hybrid nature of the optimization problem, as the remaining three decision variables: inner pipe diameter, return temperature, and supply temperature, belong now to the set of real numbers. Having the other variables included in the simulation vector, the design space takes the following form:

xd1 ¼ ½ U T S

T R T

_ hl xd2 ¼ ½ m

Q NPP

 T xd3 ¼ n pmg

ð27Þ Q

T

P

T

ð28Þ ð29Þ

In order to treat the multi-criteria nature of the optimization problem, the weighted-sum scalarization method is proposed. It fits the described problem, as electric power consumption and heat losses are naturally transformed into operation costs by multiplication by electricity/heat weight factors, allowing summation of lifetime operation costs and construction costs. Hence, the optimization problem is transformed into the single-objective problem, the solution of which is one point of the Pareto optimal set. The values of the weighting factors are specified by the decision maker prior to the optimization. As a result of application of the weighted-sum method, the scalar objective function is obtained:

Jðxd Þ ¼ wh  DQ TL ðxd Þ þ we  W TL ðxd Þ þ wc  C C ðxd Þ

ð30Þ

where wh, we and wc are weighting factors representing decisionmaker’s preferences, while corresponding to heat and electricity price coefficients (m.u./GJ), (m.u./MW h). Changes of prices may be taken into account by incorporating wh, we and wc to Eqs. (17)–(21), hence different price coefficients for each LHTD interval may be specified as Eqs. (31)–(33). Finally, the objective function is given by Eq. (34).

DQ 0TL ¼

kt X M h X

wht  kQ  DsT t  DsY y  DQ Ly;t

i

ð31Þ

t¼1 y¼1

W 0TL

" # kt  X kw

M _ ¼ sumy¼1 wet  ny;t  pg þ pmgy;t  DsT t  DsY y  V y;t t¼1

gp

ð32Þ kt X 



 wct 2L  f d ðUÞ þ cps  n1;kt

ð33Þ

Jðxd Þ ¼ DQ 0TL ðxd Þ þ W 0TL ðxd Þ þ C 0C ðxd Þ

ð34Þ

C 0C ¼

t¼1

The size of the decision problem is described by three parameters: the number of LHTD intervals kt, the number of SHTD submodels M, and the number of spatial discretization nodes k, multiplied by two (for the supply and return line). The total number of variables is equal to the sum of the decision, operational and auxiliary variables. The length of vector xd1 is equal to three. Each of M _ sub-models for each of kt intervals is described by mass flow m, head loss hl, heat absorbed from NPP QNPP, number of pumping stations n and the main pump pressure gain pmg, as well as by heat power Q, temperature T and pressure P for each node. In total, there are n1 = 3  M  kt + 3  M  kt  k  2 = 3M  kt(1 + 2k) operational variables and 2  M  kt auxiliary variables. The analysis of the equality constraints vector h(xd) shows, that each hj(xd) aggregates n2 = 3 + 6k constraints, giving total of M  kt(3 + 6k) = 3M  kt(1 + 2k) = n1 equality constraints. Hence, the total number of equality constraints is equal to the total number of simulation variables, what means that the system of equations in the simulator is square. Flowchart of the optimization procedure is shown in Fig. 7. The individual elements and corresponding equations of framework scheme and flowchart are introduced and discussed in details in Sections 3.4. 3.1. Optimizer The optimization is performed using the interior point method, where the problem with inequality constraints Eqs. (23), (25) and

171

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

each iteration, the Newton step is attempted by solving the Karush-Kuhn-Tucker equations of (35)–(37). If Hessian of the Lagrangian, approximated using the Broyden–Fletcher–Goldfarb–S hanno (BFGS) method, is not positively defined, the step cannot be taken, and the conjugate gradient step is attempted. At each iteration of the algorithm, the productivity of the step is checked against the merit function:

Start

Set starting points:

 d2 ; x  d3 Þ  l /v ðxd1 ; sÞ ¼ Jðxd1 ; x

Nic X

  d2 Þ þ s ð38Þ logðsic Þ þ v gðxd1 ; x

ic¼1

with the parameter v > 0 increasing with the iteration number. If the attempted step does not decrease the merit function, it is rejected, and the new step is attempted. For the Newton step, backtracking is used to determine the new, shorter step (at the same direction), while the conjugate gradient is computed for the new, smaller trust region. The stopping criteria have the form of lower relative bounds on the step size and the objective function value change.

Find starting points for simulator , using the sequential method (see section 3.2.1)

Find

as a solution to

3.2. Simulator

using Eqs. (41-44)

Find

Due to problem decomposition, the system of equations tackled in the simulator acquires the diagonally block structure. Each interval and sub-model is independent, what results from considering only the steady-states and xd1 parametrization. Therefore, the objective functions are the only components that link the sub-models within the structure. This provides opportunity for decomposition and simultaneous computation, as a consequence. Two approaches to simulator designing are proposed in the paper. One of them takes advantage of the knowledge of physical nature of the system to resolve it sequentially, while the other utilizes systemic approach, leaving room for numerical decomposition.

from Eq. (45)

Find new using interiorpoint method, Eqs. (35-38)

3.2.1. Sequential method The system of nonlinear equality constraints hðxd1 ; xd2 Þ consists  of hj ðxd1 ; xd2 Þj¼1;...;kt M component systems. Due to parametrization  of the decision vector xd , individual hj ðxd ; xd Þ systems 1

No

Stopping criteria met? Yes Stop

Fig. 7. Optimization procedure flowchart.

(26) is transformed into Eqs. (35)–(37) by introducing the logarithmic-barrier function which reduces it to a sequence of nonlinear equality constrained problems: N ic X d2 ; x  d3 Þ  l min Jðxd1 ; x logðsic Þ xd ;l 1

xd1

min

ð35Þ

ic¼1

 d2 Þ þ s ¼ 0 gðxd1 ; x 6 xd1 6 xd1max

1

2

d xd ¼ x 1

1

are independent, thus can be solved in parallel. Utilizing the  d1 ; xd2 Þ component knowledge of the physical nature of hj ðx j¼1;...;kt M

ð36Þ ð37Þ

in which the sequence of barrier parameters l converges to zero, and the slack variables s are restricted to positive real numbers. The algorithm combines line search (Newton step) with trustregion steps (using conjugate gradient), as presented in [21]. At

equations, it is possible to formulate certain assumptions, that allow to break the chain of mutual dependences. Thereafter, further equations can be solved in a sequence. The thermal and hydraulic models are linked by the mass flow, as it is determined in the former part and utilized in the latter. However, to determine the mass flow for each sub-model, the heat power and the temperature drop on NPP or DH heat exchanger must be known (Eq. (3)). In the considered case, only the demanded heat power at DH side and the temperature drop at NPP side are known a priori. The former quantity is one of inputs to the decision model, while the latter is found by the optimizer. In order to calculate either the nodal temperature or the heat power (Eqs. (1) and (2)) in the system, needed to obtain temperature drop at DH side or heat power at NPP side, the mass flow it to be known. One way to break the chain of mutual dependences is to approximate the heat losses (Eq. (1)) over the whole length of the HTS for both, the supply and return lines, by putting Dl = L and assuming the heating medium temperature for the resulting section equal to TS for both lines. The resultant approximation of the e L is far from accurate, as the assumption of small heat losses D Q y;t

e NPP from Dl is neglected, but it allows further approximation of Q y;t _ ~ y;t can be determined using Eq. (4). At this point, the mass flow m Eq. (3). Based on the inaccurate mass flow approximation, the e i;u;y;t e i;u;y;t from Eq. (2) and the heat power Q nodal temperature T

172

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

from Eq. (1) for HTS can be calculated. The temperature at the last e k;2 which is obtained in the above way difnode of the return line T fers significantly from the return temperature TR. However, at this point it is possible to evaluate the heat losses again, this time as the sum of differences between the heat power at the first and last nodes of the supply and return line:

e 1;1  Q e k;1 eS ¼ Q DQ y;t e 0 ¼ DQ e S þ DQ eR DQ y;t y;t Ly;t

ð39Þ

e L is higher than The first approximation of heat losses D Q y;t e 0 , as the former value is obtained at the assumption of conDQ Ly;t stant temperature for the whole line, while the latter takes into account temperature drops over small sections. Consequently, e0 e 0 and Q ~_ 0y;t approximated based on D Q the mass flow m Ly;t NPP y;t is e0 ~_ y;t . The T lower than m i;u;y;t calculation shows that the temperature at the last node converges to TR. The sequential recalculation of heat losses, mass flow, nodal temperature, and heat power, is repeated until the inaccuracy, measured as Eq. (40), drops to the desired level.

e ¼ absj Te k;2  T R j

ð40Þ

3.2.2. Two-phase method This method presents a more general approach to the solution  of the system hðxd1 ; xd2 Þx ¼x ¼ 0. Unlike the rigid structural d1

d1

decomposition which treats each sub-model separately, this approach leaves room for decomposition based on structural properties of the system as a whole. This method is also more flexible than the sequential approach, as it is capable of solving an extended problem in which individual sub-models are no longer independent. Such a case arises, for example, when heat storage is considered. Such a case arises for example, when heat storage is considered. d1 ; xd2 Þ is nonlinear and The problem of finding the roots of hðx large scale. Therefore, it was decided to utilize the trust-regiondogleg algorithm, where the original problem is approximated by the quadratic model, that is trusted to adequately represent the original problem within the trust region. The model (merit) function used to decide whether the xd2;kþ1 is better or worse than xd2;k , is defined as follows:

mk ðdÞ ¼

1 T hðxd2;k ÞT hðxd2;k Þ þ d J h ðxd2;k ÞT hðxd2;k Þ 2 1 T þ d J h ðxd2;k ÞT J h ðxd2;k Þd 2

ð41Þ

d1 was dropped where k is the iteration subscript. The notation of x d1 ; xd2 Þ: for simplification and J h ðxd2;k Þ is the n1  n1 Jacobian of hðx

where Dk > 0 is the trust region radius and D is the diagonal scaling matrix. The Powell dogleg procedure is used to compute the step d [22,23]. d1 ; xd2 Þ ¼ 0 depends on the startThe effectiveness of solving hðx ing point selection. In general, the closer xd20 is to the actual roots of the system xd2 , the faster the algorithm will converge. As the system of equations in the simulator needs to be solved a number of times for each optimizer iteration step, it is crucial for the efficiency of the whole procedure to provide a starting point in the vicinity of the root. Therefore starting point generation becomes d1 ; xd2 Þ ¼ 0, the first phase of the procedure solving the system hðx while finding of the roots themselves is its second phase. 3.3. Starting points generation module There are several possible approaches to the problem of starting point selection. The least effective way is to pick a starting point at random. Selecting a fixed starting point, either arbitrary or based upon expected values of operational variables, is not satisfactory  d1 ; xd2 Þ either. This is due to the fact, that each hj ðx describes j¼1;...;kt M a sub-model working with different conditions. Efficient starting point generation requires utilization of the information contained within the design vector xd1 . This requirement is met by the sequential method described in Section 3.2.1. It can provide the ~d to hðx d1 ; xd2 Þ ¼ 0, which is depenapproximate solution xd20 ¼ x 2 dent on the actual design vector values and adjusted to each d1 ; xd2 Þ. The exact xd is obtained in few steps in the second hj ðx 2 phase of the method. 3.4. Pumping stations module The number of pumping stations is an integer and cannot be handled by the relaxation or parameterization approach. However, it can be effectively related to the existing operational variables by applying logical conditions in the form of Eq. (45). The presented conditions prevent the pressure at each node, in either the supply or return line, from dropping below chosen minimal values. However, in that case binary variables need to be introduced to represent the logical conditions as constraints. To avoid that difficulty, the problem of finding the number and locations of pumping stations is separated from the simulator and introduced as a new module (Fig. 6). This can be applied because the operational variables xd2 are independent of the number of pumping stations. As a result, the simulator deals only with real numbers while comput d2 . ing xd2 , and xd3 is determined upon the basis of parametrized x The relations representing the pumping stations module for the single sub-model are as follows:

if P2:k;1 þ nsy;t  pg < pmin;S then nsy;t ¼ nsy;t þ 1 if P2:k;2 þ nsy;t  pg þ nry;t  pg < pmin;R

3

then nry;t ¼ nry;t þ 1

7 6 6 rh ðx ÞT 7 2 d2;k 7 6 7 6 J h ðxd2;k Þ ¼ 6 7 . .. 7 6 5 4 T rhMkT ðxd2;k Þ

ny;t ¼ nsy;t þ nry;t þ 1

2

rh1 ðxd2;k Þ

T

ð42Þ

To obtain the next step, the solution of the sub-problem is sought:

ð45Þ

pmgy;t ¼ P1;1  Pk;2 þ ðnsy;t þ nry;t Þ  pg where pmin,S and pmin,R are the limits on the minimal pressure in the supply and return line, respectively, while nsy;t and nry;t are the numbers of pumping stations in, respectively, the supply and return line. 4. Case study

min mk ðdÞ d

s:t: kD  dk 6 Dk

ð43Þ ð44Þ

The presented methodology has been verified and analyzed based on the case study of the NPP planned to be built by Lake _ Zarnowiec, Northern Poland. The main heat customer of the poten-

173

P. Hirsch et al. / Applied Energy 182 (2016) 164–176 Table 1 Case study parameter values.

120

Terrain elevation profile [m]

100

80

60

40

20

0 0

5

10

15

20

25

30

35

40

Distance [km]

Parameter

Symbol

Value

HTS length Spatial discretization step Number of nodes Pipe roughness Insulation thermal conductivity Insulation thickness Nominal pressure Minimal supply line pressure Minimal return line pressure Pump pressure gain Density of water Specific water heat capacity System lifetime Lifetime discretization step Number of intervals Sub-model duration

L (km) Dl (m) k n (mm) kis (W/(mK)) d (m) pn (MPa) pmin S (MPa) pmin R (MPa) pg (MPa) q (kg/m3) Cp (J/(kgK)) LT (years) DsT (years) kt DsY (days)

Number of sub models Ground temperature

m Tg (K)

Heat price Electricity price Construction cost weight factor Pumping station cost Pump efficiency Heat exchanger efficiency Maximum flow velocity Maximum pressure drop per 1 km Maximum heat loss per 1 km

ph (€/GJ) pe (€/MW h) pc cps (M€)

40 500 81 [0.1, 0.15, 0.2] [0.02, 0.03, 0.05] [0.2 0.3] 1.5 0.2 0.02 1 958 4220 60 [20 20 20] 3 [25, 30, 35, 90, 45, 18, 122] 7 [274, 275, 277, 281, 283, 287, 291] [12, 18, 24] [60, 64, 68] [1, 1, 1] 3.5 0.75 0.8 5 1 0.2

_ Fig. 8. Gdynia-Zarnowiec terrain elevation profile.

tial CHP NPP would be Gdynia city and its surroundings (about 300’000 residents), which is approximately 40 km distant from the intended NPP site. Knowledge of the location of NPP and DH allows specification of the decision model inputs. The Gdynia-Z_ arnowiec terrain elevation profile is shown in Fig. 8. The annual ground temperature at 1.5 m below the ground surface and the DH heat power demand in the form of structured graphs are shown in Fig. 9. Further, bounds on the design vector xd1 are specified, with regard to integration with the existing DH network:

2

3

0:1

2

U

3

2

2

ð46Þ

where the supply and return temperatures are specified in (K). The pipeline construction cost as a function of inner pipe diameter is estimated by Eq. (47). The initial point is specified in Eq. (48).

f d ðUÞ ¼ 1:5  107  U2 þ 2  106  U þ 5  105

ð47Þ

3 3 2 0:7 U 7 6 7 6 ¼ 4 T S 5 ¼ 4 418:15 5 328:15 TR 2

ð48Þ

The parameter values for the case study are collated in Table 1. Three fixed LHTD intervals, of 20 years each, are considered. The increasing values of pipe roughness and insulation thermal con-

ductivity for subsequent intervals represent aging of the HTS. Also, the variability of heat and electricity prices is taken into account, by ph and pe weighting factors, while all three components of the objective function are treated equally (corresponding weighting factors are equal to prices). Seven variable time step sub-models are considered, according to the structured graph shown in Fig. 8. As the result, the decision problem size is given by n1 = 10269 and n2 = 489, representing the number of operational variables and equality constraints of single sub-model, respectively. Two possible insulation thickness cases are considered. It is assumed that thicker insulation results in 20% pipeline construction cost increase in Eq. (47). For the purposes of this case study, the impact of heat convection in water and ground, as well as that of pipe wall thickness and pipe wall conductivity in Eq. (1) are omitted as negligibly small compared to the insulation impact. The presented methodology was implemented and solved in MATLAB R2013b environment, with the use of Optimization Toolbox. 4.1. Results discussion

18 250

16 14

200

12 10

150

8 100

6 4

50

Ground temperature Tg [°C]

DH heat power demand QDH [MW]

DPmax (bar/km) DQmax (MW/km)

20

300

2 0

vmax (m/s)

3

7 6 7 6 7 6 4 383:15 5 6 4 T S 5 6 4 423:15 5 363:15 323:15 TR

xd10

gp ge

0

50

100

150

200

250

300

350

0

Days Fig. 9. Structured graph of annual ground temperature and heat demand.

The results of the optimization procedure for the analyzed case are presented in Table 2. It can be observed, that the difference between the supply and return temperatures is maximized, resulting in TS and TR values equal to the upper and lower bound, respectively. The construction costs are the main expenditure being responsible for about 50% of the total costs. The results shown in the third part of Table 2 refer to the y = 1, t = 3 sub-model which is characterized by the highest heat power demand (y = 1) and the worst insulation and pipe quality due to aging (t = kt). Therefore, the corresponding results for other sub-models are significantly lower. For example for y = 7, t = 3 sub-model, only the main pumping station is in operation (n = 0) with the flow velocity v7,3 = 0.86 m/s. The results for the both examined insulation thicknesses are similar, as the increased construction costs are balanced by the decreased heat power loss costs. However, the second case

174

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

Table 2 Optimization results. Variable Design variables Inner pipe diameter Supply temperature Return temperature

Symbol

Results for d = 0.2

Results for d = 0.3

U (m) TS (K) TR (K)

0.6207 423.15 323.15

0.6322 423.15 323.15

C0 C DQ0 TL W0 TL J

157.72 112.68 45.67 316.07

190.30 84.06 40.89 315.25

5 2.08 0.79 0.094 3.75

5 2 0.71 0.07 2.8

1.24

0.93

1.96

1.47

700

300 0.55

400 0.6

410

0.65

Φ [m]

0.7

420

0.75

TS [K]

Fig. 11. Objective function for TR = 323.15 K.

(d = 0.3) has slightly better performance. The flow velocity, the pressure drop per kilometer and the heat loss per kilometer shown in Table 2 are within a similar range to reported in other related projects [3,9]. Both, the pressure drop per kilometer and the heat losses are small enough to make the heat transportation for even longer distances technically possible. Given that the transported heat is worth 4507 M€, using the assumed heat price, the optimally designed HTS for the considered case is highly cost effective. The pressure profile and the pumping station locations for the case d = 0.2 are shown in Fig. 10 and are marked by blue (supply), green (return) and black (main pumping station) dots. Significant impact of the terrain topology on the overall pressure profile and pumping station locations can be observed. In particular, different pumping station locations are needed for the supply line and the return line. The objective function surface graph, flattened to three dimensions, is shown in Fig. 11. Many parallel flat valleys, hindering the optimization process, can be observed. Each valley corresponds to the specific number of pumping stations. The optimization problem is non-convex, as shown in the zoom on the objective functions represented in two dimensions in Fig. 12. For this reason, the presented optimization algorithm is sensitive

345

200

150 335

100 J C'c

Δ Q'TL

325

50

W'TL

315

Objective function components [M€]

DQ R1;kt (MW) (%)

Return line heat loss

J [M€]

390

Highest load sub-model parameters Number of pumping stations n1,kt Flow velocity v1,kt (m/s) Pressure drop per 1 km DP1,kt (bar/km) Heat loss per 1 km DQ1,kt (MW/km) Supply line heat loss DQ S1;kt (MW) Heat loss rate

500

400

Objective function J [M€]

Costs (M€) Construction costs Heat power loss costs Pumping power costs Total costs

600

0

0.55

0.6

0.65

0.7

0.75

Inner pipe diameter Φ [m] Fig. 12. Objective functions for TR = 323.15 K and TS = 423.15 K.

to the choice of starting points and can converge in different valleys for different starting points. 4.2. Sensitivity analysis

x 106 2 Supply line Return line

1.8

Pressure profile [Pa]

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

30

35

Distance [km] Fig. 10. Pressure profile with marked pumping station locations.

40

The presented optimization procedure is analyzed for sensitivity to parameters in reference to the considered case study with insulation thickness equal to 0.2 m. The results are shown in Table 3. Five parameters are considered: electricity and heat prices, insulation thermal conductivity, annual heat demand, and pipe roughness. For each parameter value of the objective function, the inner pipe diameter and the pumping station number are presented. In the effect of variations of electricity and heat prices, corresponding objective function is scaled. In each case, the construction costs are the same; hence, the optimal pipe diameters for the same pumping station number are equal (the optimal solution is in the same valley). It can be observed that the increase of insulation thermal conductivity has virtually the same effect as the heat price increases, and results in both cases in pipe diameter reduction and the increased number of pumping stations. The parameter with strongest influence on the pipe diameter, and consequently on the total costs, is the annual heat demand. Therefore accurate estimation of this parameter is crucial for proper HTS evaluation. The ratio of the heat transported (demanded) to the

175

P. Hirsch et al. / Applied Energy 182 (2016) 164–176 Table 3 Optimization procedure sensitivity to parameter values.

x 106 2

40%

20%

0%

Electricity price pe J (M€) 296.7 U (m) 0.604 n 6

306.8 0.621 5

316.1 0.621 5

324.2 0.643 4

331.9 0.643 4

Heat price ph J (M€) 270.2 U (m) 0.643 n 4

293.4 0.643 4

316.1 0.621 5

338.6 0.621 5

361.1 0.621 5

Insulation thermal conductivity kis J (M€) 269.1 292.8 U (m) 0.641 0.642 n 4 4

316.1 0.621 5

339.2 0.622 5

362.4 0.623 5

20%

Supply line Return line

40% 1.8

Pressure profile [Pa]

1.6 1.4 1.2 1 0.8 0.6 0.4

Annual heat demand QDH J (M€) 244.7 U (m) 0.512 n 5 Pipe roughness n J (M€) 309.6 U (m) 0.608 n 5

281.7 0.571 5 313.2 0.615 5

316.07 0.621 5

348.6 0.665 5

316.07 0.621 5

379.9 0.706 5

318.5 0.626 5

0.2 0

0

5

15

20

25

30

35

40

Distance [km]

320.6 0.63 5

total cost is very significant – it shows that 40% increase of the transported heat results in only 20% increase of costs. Figs. 13 and 14 present nodal pressure and pumping station locations for extreme values of pipe roughness. It can be observed that the pumping station locations are virtually the same, regardless of pipe roughness. This is due to the optimization which balances higher pressure losses caused by increased pipe roughness by increasing the pipe diameter accordingly.

10

Fig. 14. Pump station locations for +40% n.

Table 4 Optimization time comparison. Starting points selection method

Optimization time (s)

Two-phase Fixed, based on experience Random, (0, 1000) interval

109.1 267.9 381.7

5. Conclusions 4.3. Comparison of starting points selection methods Starting point generation is the first phase of simulator problem solution. Table 4 compares the optimization times corresponding to particular starting points selection methods. The computations were performed on a computer with Intel Core I7-2630QM processor. Utilization of the sequential method for starting points generation decreases the optimization times by more than twice, as compared to the fixed starting point selected on the basis of experience. In comparison with random starting point selection from (0, 1000) interval, the two-phase method is 3.5 times less time consuming. x 10 6 2 Supply line Return line

1.8

Pressure profile [Pa]

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

5

10

15

20

25

30

Distance [km] Fig. 13. Pump station locations for 40% n.

35

40

Utilizing the heat produced at NPP for DH may lead to the increase of thermodynamic efficiency of the NPP. However, the distance of the NPP from large urban centers which can have sufficiently high heat demand often exceeds 100 km. Heat transportation via pipeline systems on such long distances is not practiced. Therefore, there is a need for technical feasibility and economic effectivity evaluation of such CHP NPP projects. In this paper the evaluation method is proposed in the form of the mathematical programming problem for optimized selection of design and operating parameters of the long distance HTS. The constraints to the problem are introduced with the aid of the static HTS model which takes advantage of long and short horizon time discretization to represent the system life cycle and annual time variability. In this way the HTS analysis can comprehend phenomena such as pipe and insulation aging, variability of heat and electricity prices, as well as annual heat demand and ground temperatures. In addition, the proposed method takes into account spatial topology of the HTS to evaluate its impact on the nodal pressure and localization of pumping stations. Three objective functions are considered, which are HTS construction cost, total pumping power, and total heat loss. The proposed solver to the optimization problem is based on certain problem characteristics. In particular, the decision space is categorized to specify design, operational and auxiliary vectors, each enclosing a different set of variables. Based on this categorization, the framework of the solver is proposed. It includes problem decomposition and introduction of modular structure. The optimizer, the starting points generation module, the solver, and the pumping stations module are introduced. Each module operates on different part of the decision space, while the remaining part is treated as a parameter. As a result, the individual modules are simplified. This approach, together with relaxation of some variables, enables to overcome hybrid nature of the problem. To deal with multi-criteria, scalarization with a weighted sum method is

176

P. Hirsch et al. / Applied Energy 182 (2016) 164–176

proposed. Optimization procedures for the optimizer and solver are discussed and described. The proposed approach is verified by the HTS case study. Finally, the sensitivity analysis regarding the uncertainty of parameters and the discussion of starting point selection methods, are presented. Future work in the field of effectivity and feasibility analysis of long distance HTSs may extend the decision-making problem with the possibility of energy storage. A possibility of pipeline system division to several smaller lines can also be evaluated. Acknowledgments This work was supported by the National Centre for Research and Development – Poland, under Strategic Research Project No. SP/J/10/176450/12. The authors wish to express their thanks for the support. References [1] European Commission. A European strategic energy technology plan (SETPlan)—towards a low carbon future. COM 2007:723. Final, 2007. [2] Fütterer MA, Carlsson J, de Groot S, Deffrennes M, Bredimas A. European energy policy and the potential impact of HTR and nuclear cogeneration. Nucl Eng Des 2014;271:73–8. http://dx.doi.org/10.1016/j.nucengdes.2013.11.013. [3] Safa H. Heat recovery from nuclear power plants. Int J Electr Power Energy Syst 2012;42:553–9. http://dx.doi.org/10.1016/j.ijepes.2012.04.052. [4] Tokimatsu K, Konishi S, Ishihara K, Tezuka T, Yasuoka R, Nishio M. Role of innovative technologies under the global zero emissions scenarios. Appl Energy 2016;162:1483–93. http://dx.doi.org/10.1016/j.apenergy.2015.02.051. [5] Colmenar-Santos A, Rosales-Asensio E, Borge-Diez D, Collado-Fernández E. Evaluation of the cost of using power plant reject heat in low-temperature district heating and cooling networks. Appl Energy 2016;162:892–907. http:// dx.doi.org/10.1016/j.apenergy.2015.10.161. [6] Ren´ski A, Duzinkiewicz K. Development and implementation of an exemplary method of system analysis of the nuclear power plant with light water reactor in partial cogeneration work. Gdan´sk University of Technology: The National Centre for Research and Development research project for years 2012–2014, final report, project leader: A. Ren´ski; 2014 [in Polish]. [7] Hirsch P, Grochowski M, Duzinkiewicz K. Pipeline system for heat transportation from nuclear power plant - an optimizing approach. IEEE 2015:1044–9. http://dx.doi.org/10.1109/MMAR.2015.7284023. [8] Minkiewicz T, Ren´ski A. The possibility to use a nuclear power plant as a source of electrical energy and heat. Acta Energy 2014:114–8. http://dx.doi.org/ 10.12736/issn.2300-3022.2014310.

[9] Paananen M, Henttonen T. Investigations of a long-distance 1000 MW heat transport system with apros simulation software. Finland: Espoo; 2009. [10] Wang Y, Tian C-H, Yan J, Huang J. A survey on oil/gas pipeline optimization: problems, methods and challenges. Suzhou: IEEE; 2012. p. 150–5. http://dx. doi.org/10.1109/SOLI.2012.6273521. [11] Fettaka S, Thibault J. Pipeline optimization using a novel hybrid algorithm combining front projection and the non-dominated sorting genetic algorithmII (FP-NSGA-II). Cancun, Mexico: IEEE; 2013. p. 697–704. http://dx.doi.org/ 10.1109/CEC.2013.6557636. [12] Asim T, Mishra R, Kollar LE, Pradhan SR. Optimal sizing and life-cycle cost modelling of pipelines transporting multi-sized solid–liquid mixtures. Int J Press Vessels Pip 2014;113:40–8. http://dx.doi.org/10.1016/j. ijpvp.2013.11.003. [13] Ríos-Mercado RZ, Borraz-Sánchez C. Optimization problems in natural gas transportation systems: a state-of-the-art review. Appl Energy 2015;147:536–55. http://dx.doi.org/10.1016/j.apenergy.2015.03.017. [14] Luo X, Wang M, Oko E, Okezue C. Simulation-based techno-economic evaluation for optimal design of CO2 transport pipeline network. Appl Energy 2014;132:610–20. http://dx.doi.org/10.1016/j.apenergy.2014.07.063. [15] Hlebnikov A, Dementjeva N, Siirde A. Optimization of Narva District heating network and analysis of competitiveness of oil shale CHP building in Narva. Oil Shale 2009;26:269–82. http://dx.doi.org/10.3176/oil.2009.3S.09. [16] Razvan - Corneliu L, Daniela P, Gabriel IP. Establishing the optimal solution for retrofit of district heating networks. Case study, Iasi. Romania: IEEE; 2014. p. 1083–8. http://dx.doi.org/10.1109/ICEPE.2014.6970075. [17] Dalla Rosa A, Li H, Svendsen S. Method for optimal design of pipes for lowenergy district heating, with focus on heat losses. Energy 2011;36:2407–18. http://dx.doi.org/10.1016/j.energy.2011.01.024. [18] Ikonen E, Selek I, Kovacs J, Neuvonen M, Szabo Z, Bene J, et al. Short term optimization of district heating network supply temperatures. Cavtat, Croatia: IEEE; 2014. p. 996–1003. http://dx.doi.org/10.1109/ ENERGYCON.2014.6850547. [19] Jie P, Kong X, Rong X, Xie S. Selecting the optimum pressure drop per unit length of district heating piping network based on operating strategies. Appl Energy 2016;177:341–53. http://dx.doi.org/10.1016/j.apenergy.2016.05.095. [20] Duzinkiewicz K. Integrated control of drinking water supply and distribution systems. Kraków: AGH University of Science and Technology Press; 2005 [in Polish]. [21] Waltz RA, Morales JL, Nocedal J, Orban D. An interior algorithm for nonlinear optimization that combines line search and trust region steps. Math Program 2006;107:391–408. http://dx.doi.org/10.1007/s10107-004-0560-5. [22] Powell MJD. A FORTRAN subroutine for solving systems of nonlinear algebraic equations. In: Rabinowitz P, editor. Numerical methods non-linear Algebr Equ. Gordon and Breach; 1970. p. 115–61 [n.d.]. [23] Nocedal J, Wright SJ. Numerical optimization. 2nd ed. New York: Springer; 2006.