Multicriteria optimization approach to design and operation of district heating supply system over its life cycle

District Heating (DH) systems are commonly supplied using local heat sources. Nowadays, modern insulation materials allow for effective and economically viable heat transportation over long distances (over 20 km). In the paper a method for optimized selection of design and operating parameters of long distance Heat Transportation System (HTS) is proposed. The method allows for evaluation of feasibility and effectivity of heat transportation from the considered heat sources. The optimized selection is formulated as multicriteria decision-making problem. The constraints for this problem include a static HTS model, allowing considerations of system life cycle, time variability and spatial topology. Thereby, variation of heat demand and ground temperature within the DH area, insulation and pipe aging and/or terrain elevation profile are taken into account in the decision-making process. The HTS construction costs, pumping power, and heat losses are considered as objective functions. Inner pipe diameter, insulation thickness, temperatures and pumping stations locations are optimized during the decision-making process. Moreover, the variants of pipe-laying e.g. one pipeline with the larger diameter or two with the smaller might be considered during the optimization. The analyzed optimization problem is multicriteria, hybrid and nonlinear. Because of such problem properties, the genetic solver was applied.


Introduction
Space heating is a major component in overall heat consumption in Europe.At the same time heat uses amount to 80% of total energy consumed in residential houses.District Heating (DH) is an attractive solution for space heating, especially when it is supplied by Combined Heat and Power (CHP) plant [1,2].Recovering the waste heat from any kind of thermal power plant allows to greatly increase the plant energy efficiency and is an effective method for cutting carbon emissions.Several studies have shown that conventional CHP plants are highly attractive from economic and environmental perspective [3], however the DH area have to be close to the plant site.Another important Heat Source (HS) for DH systems can be CHP Nuclear Power Plant (NPP).Although the CHP NPP can be a source of a carbon-free heat for residential and industrial use, it is seldom considered as a viable HS for DH.The main obstacle in a nuclear heat utilization is a necessity of a long distance heat transportation, as NPPs are most often located far from dense urban areas.For these reasons, the NPPs are considered solely for electricity production under global zero emission scenarios [4].
Despite the potential investment costs on pipelines and heat losses on transportation, the CHP NPP is still being investigated.In [5] a 77 km long 1000MW HTS has been analyzed for possible CHP production at the Loviisa 3 NPP project in Finland.The simulations presented by the authors showed, that only 11 MW of the heat was lost on transportation.In [1] economic estimations were performed for the French Nogent-sur-Seine NPP located 110 km east of Paris, showing that heat transportation of 1500 MW over long distance can be cost effective, with the payback time less than 10 years.
DH networks are often subject to detailed analyze in order to find optimal pipes sizing [6] or minimize the heat losses [7,8].Similar approaches can be utilized for DH supply system, that would transport heat over long distances, from the HS to the DH area.In [2,9] desired characteristics of the HTS efficiency analysis method are discussed.The method should utilize optimization techniques, to find the design minimizing construction costs, thermal losses and pumping power.The design variables should be: inner pipe diameter, insulation thickness, supply and return temperatures and number of pumping stations.The performance of the HTS changes in time due to the energy market shifts as well as pipe and insulation aging processes.Therefore, it was decided to analyze the HTS over presumed system life cycle.Other important conditions that affect the HTS efficiency are: the heat demand and its annual variability, terrain profile and annual ground temperature distribution.The method presented in the paper is extended with the pipe-lying strategy evaluation, allowing consideration of multiple, smaller pipelines working in parallel.Moreover, the revised mathematical model is presented, resulting in much better computing performance, allowing effective genetic solver application.

Decision model for optimization based efficiency analysis
The decision model for the HTS optimization based efficiency analysis consists of the mathematical model of objectives and constraints, and optimization problem formulation.Upon the decision model formulation, the following assumptions on the HTS operation are considered: the heat power delivered to the DH is equal to the heat demanded; the outlet and inlet temperatures of the HS heat exchanger are the supply and return temperatures, respectively; the pumping stations pressure gain is constant and independent to heating water flow rate.Further assumptions on the decision model computation are: only steady states are considered, the minor pressure losses are neglected.
Evaluation of the HTS design and operation over its life cycle must take into account the HTS operation conditions and parameters variability in time.The time-variant processes affecting the HTS are occurring on significantly different time scales.The pipeline and insulation aging are examples of processes with time constants of several years.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 LC, with ΔτL time step, into ∑   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 be either fixed or variable.The HTS is also subjected to significant annual variability due to annual heat demand and ground temperature changes.To deal with these conditions, the Short Horizon Time Discretization (SHTD) is proposed.The SHTD time base is one year, which is divided with the variable time step ΔτY into ∑   intervals indexed by y, in accordance to discretized, structured graph of annual heat demand (see Fig. 3b.).As a result, ∑   sub-models are created for each of the LHTD intervals.Each sub-model is corresponding with a different heat demand level (different part of the year), and is characterized by different flows and temperatures.Although the duration of each LHTD interval may be decades, its parameters are assumed constant over that period and, consequently, it may be covered by only one set of ∑   sub-models, the results of which are multiplied according to the LHTD time step.
According to the HTS topology scheme shown in Fig. 1, the HTS is a network composed of links and nodes.The links are pipeline sections, pipeline sections with the pumping stations, or heat exchangers.The nodes are allocated according to the spatial discretization, which is introduced in order to deal with the HTS spatiality and the resulting terrain elevation profile (see Fig. 3a.) impact on the nodal pressure and on the pumping stations localization.The terrain elevation profile, given with the base resolution of Δl [m], is chosen as basis of the discretization, with the nodes being allocated at the points of the terrain profile change.Therefore, the HTS of length L [km] is discretized with the variable step Δli [m] that is a multiple of Δl, resulting in N nodes in each direction, indexed by i.The two directions are denoted with the second subscript u.

Thermal and hydraulic mathematical model
The mathematical model formulation utilizes reduced notation for conciseness.The nodal variables (describing a quantity in each node) have the subscripts i,u,y,t, meaning respectively: the node number, the line, the SHTD interval, the LHTD interval.In reduced notation, unless necessary, only the i subscript is listed.
The main purpose of the thermal model is to calculate the heat losses on transportation.The losses must be taken into account to ensure assumed equilibrium of the heat delivered and demanded.Considering the pipeline section of length of Δli [m] with diameter Φ [m], external insulation of thickness δ [m] and thermal conductivity λ W/(m•K)] in which superheated water at temperature Ti [K] is flowing (Fig. 2.), the linear heat losses can be calculated as: where Q is the heat power at given node [W] and Tg is the ground temperature [K].Note that the insulation thermal conductivity and the ground temperature in (1) are discretized according to LHTD and SHTD respectively.Another important feature of the presented linear heat loss model is approximation of superheated water temperature over the section with Ti, which is the section inlet temperature.Different approximation methods, such as inlet-outlet average or logarithmic mean can also be utilized, however if the section length is small enough, the outcomes differences are negligible.The exact value of 'small enough' section length is dependent on the case study and can vary within range from 100 meters up to several kilometers.The loss of the heat translates in direct proportional manner into the temperature loss: where Cp is the specific water heat capacity [J/(kg•K)] at 100°C, and ̇ , is the mass flow of the heating water [kg/s] at a given time step.Assuming the absence of the heating water phase change in the heat exchanger and the heat exchanger efficiency   , the heating water mass flow can be calculated from the amount of the heat exchanged and the temperature drop on the heat exchanger.This applies for both DH and HS side heat exchangers: In Eq. ( 3)   , is the heat delivered to the DH network [W], with the temperature drop of Δ  , [K].Similarly,   , is the heat absorbed from the HS [W] at temperature drop Δ  [K].The heat delivered to the DH network is equal to the heat demand and is known at any time, while the amount of heat needed from HS has to be calculated as a sum of total heat lost on transportation in both directions and the heat delivered.The same applies to temperatures -Δ  is equal to the difference between the supply and return temperatures, which are inputs to the model, while Δ  , has to be calculated using Eqs.(1-3).

Fig. 2. Pipeline section scheme with associated variables and parameters.
The proposed thermal model is characterized by a dense computational nodalization that is required because of the temperature approximation accuracy.At the same time, the spatial discretization based on terrain elevation profile is much sparser.It can be shown, that the thermal model can be calculated only for the selected nodes, without the need of computing all the base nodes in between.
Let us consider the ratio of the heat losses in two successive pipeline sections, as shown in Eq. (1): It can be observed, that the successive heat losses ratio is constant and lower than 1 at any given time step.The same applies to the temperature losses ratio, as the two quantities are proportional.It can be shown, that ∆ +1 ∆  ⁄ also equals to ky,t.Therefore, the heat or the temperature loss at any node can be found from the geometric sequences: ∆  =  −1 ∆ 1 , ∆  =  −1 ∆ 1 .By applying geometric series, the heat power and the temperature at any node can be calculated: Using the extended thermal model (Eqs.(1-6)), the heating water temperature and the transported heat power can be found at any node that is allocated at a multiple of Δl.Proposed thermal model can be solved either as the system of equations, or using sequential procedure as described in [2].
The hydraulic model is a base for calculation of the number of pumping stations required and their location, as well as the electrical power used on pumping.Since considered HTS supply and return lines are buried in a single trench, the terrain elevation has no impact on the number of pumping stations, as its impact on the pressure profile cancels out.However, the terrain has still strong impact on the pumping stations locations.
The pressure drop between consecutive computational nodes Δ  [Pa] is affected by: the linear pressure losses   , [Pa] (minor pressure losses are neglected), the elevation difference pressure    [Pa] and the possible pressure gain Pg [Pa] due to pump station allocation at the inlet node.As a consequence to assumptions of fixed pipeline diameter and closed HTS structure (without branching nodes), the flow velocity is the same for any node at any given time step.Thus, the linear pressure loss over the pipeline section is constant for any given time step and can be described by the Darcy-Weisbach equation: where ρ is the density of water [kg/m 3 ],  , is the heating water flow velocity [m/s] calculated from the mass flow and f is the Darcy-Weisbach friction coefficient, which for fully turbulent flow (Reynolds number>4000) and the inner pipe roughness of ξ [mm] can be calculated as shown in [1].In order to determine the pressure at any node allocated at multiply of Δl, an arithmetic series can be utilized to find the linear pressure losses, while the elevation pressure is equal to the difference between the initial and the nodal elevation pressures, resulting in Eq. ( 8): where z is the HTS elevation head [m], g is the Earth gravity constant [m/s 2 ] and ni is the number of pumping stations allocated between the nodes 1 and i.The pumping stations can be added at the node, when the nodal pressure decreases below the safety value or it can be subject to further analysis.

Optimization problem formulation
As a result of the time discretization, a number of sub-models are established.Each submodel is described by a set of variables.Some of these variables are the same for each submodel (i.e.inner pipe diameter, insulation thickness, supply and return temperatures), while the others, like nodal pressure or temperature, are independent.While all variables together constitute decision space   , the former are visibly more important for the decision making process.Therefore, categorization of the decision space is proposed: the linking variables shall be called design vector, denoted    , while the rest compose simulation vector    .According to this categorization, the decision space is represented as   = [       ]  .
The decision variables are subjected to a set of nonlinear equality constraints ℎ(  )that constitute the mathematical model.Furthermore the upper and lower bounds, as well as inequality constraints (  ), can be considered on the decision space due to technological and operational limitations.The range of considered supply temperatures and maximal heating water flow velocity can be examples of such constraints.
The HTS design and operation is assessed on the aspects of its technical and economic value.Technical value can be represented by terms of thermal and transport efficiency.One way to indicate the HTS thermal efficiency is to estimate its life cycle heat losses on transportation.It can be done by calculating firstly the heat losses of each sub-model ∆  , [W] (Eq.( 9)), followed by annual heat losses ∆   [GJ/year] (Eq.( 10)) and the total heat losses Δ  [GJ] on transportation (Eq. ( 11)).The transport efficiency can be represented by means of energy used on heating water pumping.At first annual pumping energy    [MWh/year] is calculated using Eq. ( 12), then the total pumping energy   [MWh] (Eq.( 13)): where kw=24 is the unit conversion coefficient,   is the pump efficiency,  , is the number of pumping stations set for the given sub-model, pg is the pressure added by the pump station [MPa],   , is the pressure added by the main pump station [Pa] and  ̇, is the volumetric flow rate [m 3 /s].
The economic value of the HTS is estimated using empirical construction cost function: where

Case study and solver
The optimization problem is multicriteria, hybrid and nonlinear.The combination of these properties makes it difficult to solve.In an attempt of simplifying the optimization problem, the weighted-sum scalarization method is proposed.Eq. ( 11) is multiplied by the heat and (13) by the electricity prices, allowing the objective functions summation.Resulting optimization problem, with pipe diameters and insulation thickness being integers and other decision variables being real, is solved using genetic algorithm solver.
The optimization problem formulation can be extended with the possibility of different pipe-lying strategies and variants evaluation.Depending on the maximal heat demand and its annual variation, it may be cost-effective to divide the pipeline into several smaller ones, working alternately or simultaneously.This functionality is introduced by allowing a number of parallel pipelines, each having the same parameters except for inner pipeline diameter, and additional decision variable, that limits the maximal heat provided by the single pipeline.If the heat demand exceed the heat provided by the current working pipelines, additional pipeline is considered, increasing the construction costs accordingly to its diameter.
The presented methodology has been verified and analysed based on the case study of a potential 40 km HTS.The heat source would be intended NPP, located by Lake Żarnowiec, Northern Poland, while the main heat customer would be Gdynia and its surroundings.Knowledge of the NPP and DH location allows specification of the terrain elevation profile (Fig. 3a.) as well as the annual ground temperature and the DH heat demand (Fig. 3b.).The number and length of the SHTD intervals can be read from Fig. 3b  The results of the optimization procedure for the analysed case, together with the decreased and increased heat demand, are presented in Table 1.Each time the maximal heat provided by single pipeline is greater than heat demand, therefore only one pipeline is selected and one pipe diameter shown.It would require greater heat demand for the two smaller pipes case to be viable.It can be observed, that benefits of thicker insulation outweigh the 15% construction costs increase for each of the heat demand cases.

Conclusions
The proposed method for optimized efficiency analysis of the HTS, can be utilized as a part of the decision support system for evaluation of projects associated with the heat transportation.An example of such project is a Combined Heat and Power NPP, where part of the technical and economic feasibility assessment task concerns high power, long distance HTS.With the appliance of the proposed method, the optimal design and operation strategy of the HTS can be found, concerning its construction costs, thermal and transport efficiency, while subjected to a number of operation and technological constraints.The constraints are a result of mathematical model and governing equations, additional can be imposed by decision-maker to ensure specific HTS operation, i.e. heating water flow velocity or heat loss over the pipeline section.Important feature of the proposed method is an ability to evaluate projects in long term or in particular, over the project life cycle.Long Horizon Time Discretization allows consideration of the pipe and insulation aging phenomena, as well as the heat and electricity prices changes.At the same time, conditions of the time scale of the month or even days are still taken into account due to Short Horizon Time Discretization.In addition, the proposed method takes into account the spatial topology of the HTS to evaluate terrain elevation impact on the nodal pressure and localization of the pumping stations.

Fig. 3 .
Fig. 3. a) The case study terrain elevation profile, b) structured graph of annual ground temperature and heat demand.
) = [  (  ) ∆  (  )   (  )]  , leading to general optimization problem formulation: min Cc is the overall HTS construction cost [m.u.], fd(Φ) is the pipeline material and labor costs, as a function of pipeline diameter, and cps is the cost of pumping station construction [m.u.].The objective functions Eqs.(9-14) can be represented as the objective vector