Pattern-Based Set Partitioning Algorithm for the Integrated Sustainable Operation of a District Heating Network

Abstract

District heating is a system of distributing heated water from centralized facilities to local homes and buildings. In this paper, we model the distribution planning problem as a supply chain planning problem and propose an explicit column generation algorithm to handle large-scale data and nonlinear constraints. The algorithm is successfully applied to a Korean district heating company and computational experiments show that the integrated operation of the district heating network increases the total profit compared to previous isolated networks.

1. Introduction

Due to global warming, unusual extreme weather conditions have dramatically increased the demand for cooling and heating during summer and winter, respectively, and energy efficiency is becoming a key issue in heating and cooling systems which heavily rely on fossil fuels such as natural gas and coal. Existing heating and cooling systems have been based on individual systems in homes and buildings, but their efficiency is not sufficiently high due to the use of small heat sources.

In order to solve this problem, district heating systems which efficiently use central heat sources and supply hot water to areas of demand through networks are becoming more important in terms of the effective utilization of energy, and many countries, including Austria and Sweden, operate such networks [1,2,3]. In South Korea, the Korea District Heating Corporation (KDHC) (https://www.kdhc.co.kr/), GS Power (http://www.gspower.co.kr/) and other energy companies operate district heating networks. KDHC mainly supplies hot water for the Seoul metropolitan area, resulting in sales of KRW 1.7 trillion (USD 160 million) as of 2017.

The operation of a district heating network can be viewed as supply chain planning, producing hot water for heating from various heat sources with different characteristics and then supplying heat to the required areas. Generally, the heat sources used for district heating include combined heat and power (CHP) plants, boilers and incinerators, and recently, solar heating, geothermal heat pumps, and fuel cells have been tested and implemented as alternative sources to fossil fuel-based heat generation. According to the 2017 annual report by the Korean district heating industry, CHP is the single largest heat supply source for the local heat demand, supplying 43% of the annual heat demand. Incinerators and the excess heat from electricity generation and manufacturing processes provide 42% of the heat demand, while fossil fuel-based heat-only boilers contribute about 14%. Alternative heat sources, including fuel cells, heat pumps and solar heating only provide 1.45% of the annual heat demand. Considering the supply mix of heat sources, it is important to optimize and balance the operation of different heat sources in the current Korean district heating network.

From a supply chain planning perspective, the optimization problem involving the determination of the number and quantity of products manufactured at different plants with different production costs and operating characteristics is equivalent to the concept of optimizing the economic operation of a district heat network. For example, if there is a discrepancy between the heating demand and the hot water supply, the over-supply of hot water can be stored in a storage tank in a similar way to holding inventory in supply chain planning.

District heating companies in South Korea supply hot steam water (105 $°$C~115 $°$C) to customer facilities. Heat exchangers in each customer facility absorb heat from the supplied hot steam water and supply both hot water (45 $°$C~60 $°$C) and space heating (45 $°$C~55 $°$C) using the heat absorbed by the heat exchangers. To supply space heating services to customers located in the Seoul metropolitan area in South Korea, KDHC initially established a number of heat production facilities near demand sites, such as Gangnam and Suwon, so that it became possible to supply the necessary heat to those demand sites. In the beginning, these facilities relied heavily on heat-only boilers (HOBs), which produce heat by burning fossil fuels, and they did not have any centralized planning system. Each production site supplied heat to local areas independently, and the operation of the HOBs was rather simple. Heat transfer between production facilities was conducted in a restricted way. However, combined heat and power plants dramatically changed the operation and economics of the district heating system. KDHC decided to convert its independent operating mode (see Figure 1a) into the integrated operating network (Figure 1b).

The existing independent network does not supply hot water for heating between the heat production facilities but supplies it locally in the corresponding area to a heat production site. If the heat production costs and heating demand do not vary much between the areas, then there is no major problem associated with operating independently. However, as new types of heat production facilities, such as combined heat and power plants, have been introduced, several issues have arisen with the existing stand-alone operations. Combined heat and power plants produce hot steam by burning natural gas in a gas turbine which is then supplied to a steam turbine and a heat generator to produce heat and electricity altogether. District heating companies with cogeneration plants need to effectively balance the production of heat and electricity. With the advent of large-scale combined heat and power plants, there has been a significant change in the heat production cost and the sales of electric power. Electric power sales are highly dependent on the market price of electricity, and there is a large variation in heat demand. This implies that the heat production from CHPs may far exceed the heat demand from an associated demand area. In this case, the excess heat may be either stored in a thermal storage or supplied to other connected areas. To fully utilize large-size CHPs, it is necessary to have both mechanisms—heat storage and network connectivity. If multiple CHPs and thermal storage sites exist, it will be complicated determining the optimal levels of heat production, transfer, and storage plan in the network. The planning system should consider system-wide profits by comparing different options with CHP and storage. As a result, the independent operation of each production facility has been converted into the integrated operation of a single network in Korea to better manage the entire district heating network.

However, the integrated approach has shifted the paradigm of regional heating network operations with the problem of a rapid increase in size. The existing district heating network is characterized by an independent operation which produces large quantities of hot water for heating with less concern about the cost of heat production, stores it in heat storage tanks, and responds to heating demand with the stored hot water. On the other hand, an integrated operation requires hot water heated in low-price areas to be moved to high-price areas, and the sales of electricity and heat needs to be optimized.

In fact, the cost-minimization operating philosophy has shifted to one of profit maximization. In the existing independent system with heat-only boilers, the operation of the facilities does not affect revenue, but rather the production costs. Since the heat demand is exogenous and needs to be satisfied at any time, the heat revenue is linked to the given heat demand which implies that the minimization of heat production costs is sufficient to maximize profit. The isolated production sites are independent of each other. In the integrated system with multiple CHP plants, the integration affects both revenue and cost structures concurrently. CHPs can adjust the ratio of heat to power production. If the electricity market is strong, they will generate more electricity which, in turn, will increase electricity revenue, but the heat supply from the CHPs will decrease. Other heat sources are required to make up for the decrease in the heat supply. If multiple CHPs interconnected in a single network exist, any system-wide planning system should optimize network-level profits by comparing the efficiency of CHP operation, heat transfer, and storage. By balancing the production of heat and electricity and by transferring heat between distant facilities, the company can maximize network-wide profits. In real-world data, there is a significant temporal difference between electricity prices and the demand for heat (Figure 2). Since combined heat and power plants need to produce more electricity during peak hours, if the electricity and heat demands mismatch temporally, any excess heat from the combined heat and power plants should be transferred to other demand areas or stored in accumulators for future use. Therefore, the integrated operation of a district heating network needs to make decisions about power generation, heat production, and heat flow between production facilities, considering the production costs of all heat sources, electricity sales, and heat sales.

Most district heating operation methods are sufficient from the viewpoint of minimizing operating costs via the utilization of boilers or incinerators, and the majority of the district heating networks are small in scale. On the other hand, for large-area district heating networks, like the Korea District Heating Corporation, which consist of many cogeneration plants, boilers, incinerators, and other heat sources, there is a strong need for integrated optimization, unlike regional, independent small-scale operations. The optimization of a district heating network that considers the integrated benefits discussed in this study is similar to that of Unit Commitment (specifically, Profit-based Unit Commitment or PBUC) [4]. PBUC concerns the optimization of a generator’s operation plan to maximize electricity sales; thus, the operation of cogeneration power plants in a district heating network can be modeled after PBUC. However, unlike PBUC, a district heating network is characterized by the fact that many cogeneration plants simultaneously produce heat and electricity, making it more difficult to obtain solutions. Because existing PBUC-related optimization algorithms only consider characteristics in terms of power generation, they should be extended and applied to local heating network optimization [5]. In this paper, we aim to optimize the operating plan so that it maximizes the overall profit in a large district heating network.

The remainder of the paper is as follows. Section 2 provides a literature review and discusses our contributions to the existing literature. Section 3 presents the mathematical programming model for the integrated sustainable operations of a district heating network. In Section 4, we propose a pattern-based set partitioning algorithm to solve the aforementioned problem in Section 3. Section 5 illustrates the real-world application of the proposed approach to address the optimization of operations in the districted heating network in South Korea. Finally, we conclude with a discussion in Section 6.

2. Literature Review

A district heating system consists of heat production and storage facilities connected to demand nodes through long-distance networks. As district heating systems have become popular in urban areas and new technologies have been introduced, researchers have paid attention to the benefits of combining various technologies in the networks, including smart and sustainable heat production resources.

Lund et al. [6] and Lund et al. [7] proposed fourth generation district heating models that included smart energy and smart thermal grids. They argued that future district heating systems will be connected as smart networks with distributed renewable energy sources. Jiang et al. [8] proposed an integrated district heating system with wind, solar, gas, and electricity resources. Based on multiple operation strategies, they compared their proposed integrated system with traditional natural gas-based district heating systems.

Li et al. [9] developed deterministic and robust models to evaluate an integrated system of wind power and combined heat and power units. They tested whether the integration was beneficial for efficient operation. Dou et al. [10] evaluated a district heating system design with multiple land-use scenarios and demonstrated that a district heating system with waste heat can achieve economic and environmental benefits. Schweiger et al. [11] studied the integrated operation of electricity generation and a district heating system. As renewable electricity production technologies are prevalent, the system suffers from a lack of flexibility due to uncertain operations. By integrating electricity and district heating systems, it is expected that the system would be able to stabilize energy production and supply.

To analyze and evaluate district heating systems, it is important to apply simulation and optimization models. In addition, because district heating systems are interconnected through networks, it is difficult to design and manage multiple facilities in an integrated way. Combining heat and power facilities is quite complicated due to the multiple parameters that need to be considered, and thus, advanced optimization models are necessary. To analyze a system in a logical way, optimization models have been proposed and tested with genuine operational data obtained from local systems.

Sartor et al. [12] analyzed a biomass-based combined heat and power system and developed simulation models for the accurate estimation of combined heat and power (CHP) connected to a district heating network. Buoro et al. [13] identified the optimal operation strategies for a district heating system with a combined heat and power plant based on mixed integer programming. They showed that the integrated operation of various energy sources can result in improved economic and environmental performance. Wang et al. [14] proposed a mathematical programming-based optimal algorithm for the planning and operation of a CHP system with an energy storage system and renewable energy sources. Vesterlund and Dahl [15] studied a district heating network containing loops. They derived a process integration model to analyze the impacts of the loops and the behaviors of the network. Li et al. [9] proposed an iterative algorithm to solve the combined heat and power dispatch problem considering both the electricity and district heating systems.

Ameri and Besharati [16] proposed a mixed integer linear programming model for the design and management of combined cooling, heating, and power networks. The proposed model was shown to minimize costs and greenhouse gas emissions. Carpaneto et al. [17] studied an algorithm for coordinating energy sources in a district heating network that included solar thermal energy. They tested the algorithm when designing and planning a district heating system. Olsthoorn et al. [18] provided a detailed review of modeling and optimization in district heating system operations. District heating systems were classified based on the optimization objectives, including operational efficiency, costs, and environmental effects. Morvaj et al. [19] developed mixed integer programming models to design a distributed energy system with district heating networks under multiple design objectives. The models were shown to be useful for evaluating the impacts of parameter changes, including CHP constraints, network heat losses, and improved modeling of thermal storage. Facci et al. [20] studied the optimal operation of CHPs and developed dynamic programming-based heuristics. Rong and Lahdelma [21] proposed a linear programming-based heuristic for a transmission-constrained multi-site CHP system. Bordin et al. [22], Mitra et al. [23] and Fang and Lahdelma [24] developed mixed integer programming formulations for CHPs and district heating network operations. Salgado and Pedro [25], Olsthoorn et al. [18] and Kumbartzky et al. [26] reviewed optimization models for district heating systems with multi-sources, including CHPs.

Most previous research has applied mixed integer programming models to optimize the design and operation of a district heating system. However, as the size of a network becomes bigger, it is important to develop scalable algorithms, even for dozens of network nodes and facilities. CHPs and network connectivity are subject to complicated non-linear logical constraints, meaning a standard mixed integer programming model could have millions of variables and constraints, even for a small number of nodes and for short planning periods (i.e., of a few days). Therefore, optimization research on district heating systems has focused on detailed short-term planning, referring to reference [25]. This paper proposes a pattern-based set partitioning model to deal with a large-scale district heating system interconnected in a network.

3. Model

In this section, we present a mathematical programming model for the integrated sustainable operation of a district heating network. Specifically, we present the models related to the heat production facilities themselves as well as the structure of the district heating network and the objective function.

3.1. Modeling of Heat Production Facilities

Based on the heat production of a heat production facility (i) at a specific time (t) being $h_{i,t}$, the heat production facilities considered in this study were classified into three types. (type-A, type-B and type-C).

First, for type-A facilities, such as incinerators, optimization is difficult during the planning stage, and waste heat recovery and incineration plans are used to produce heat; thus, in practice, an external plan should be implemented, rather than optimization. Thus, for type-A heat production facility $i\in PA$ as an external heat source, the constraint on the external heat source simply takes the external heat production plan ($S_{i,t}$) into account as follows:

$$h_{i,t}=S_{i,t},\hspace{1em}\forall i\in PA,t\in T,$$

(1)

where $PA$ is the set of type-A facilities (external heat sources) and T is the time index set whose element has a value from 0 to 8640 h (i.e., approximately 360 days).

A type-B heat production facility is a heat source dedicated to heat production, such as a boiler. Here, we let $PB$ be the set of type-B facilities. A heat production facility in $PB$ can control the production quantity according to the heating demand. Basically, a heat source for only heat production should abide by the minimum and maximum heat production capacity per hour—$C_i^{MIN}$ and $C_i^{MAX}$, respectively. We let the decision variable $y_{i,t}$ indicate the operating state of the heat source, where $y_{i,t}$ is 1 if the heat source $i\in PB$ is activated at time $t\in T$ and 0 otherwise. The heat production constraints can then be expressed as

$$\begin{array}{cc}\hfill h_{i,t}\ge C_i^{MIN}{y}_{i,t},& \forall i\in PB,t\in T\end{array}$$

(2)

$$\begin{array}{cc}\hfill h_{i,t}\le C_i^{MAX}{y}_{i,t},& \forall i\in PB,t\in T.\end{array}$$

(3)

In addition to the heat production capacity, we also considered the start-up and shut-down times. Depending on the equipment, it may take time to start and stop the heat source. The former is called the minimum operation time limit, $T{S}_i^B$, and the latter is called the minimum stop time, $T{I}_i^B$. For example, when the heat source is activated, the heat source must remain in operation for $T{S}_i^B$ hours from its start time (i.e., the condition must be maintained for at least the duration of $T{S}_i^B$ hours). Let us introduce $z_{i,t}$, which is 1 if the heat source is started at time t and 0 otherwise, and $w_{i,t}$, which is 1 if the heat source is stopped and 0 otherwise. To express the conditions discussed above, the corresponding constraints were modeled as follows.

$$\begin{array}{cc}\hfill y_{i,t+dt}\ge z_{i,t},& \forall i\in PB,t\in T,dt\in \{0,1,\dots ,T{S}_i^B\}\end{array}$$

(4)

$$\begin{array}{cc}\hfill y_{i,t+dt}\le w_{i,t},& \forall i\in PB,t\in T,dt\in \{0,1,\dots ,T{I}_i^B\}\end{array}$$

(5)

$$\begin{array}{cc}\hfill y_{i,t-1}+z_{i,t}=y_{i,t}+w_{i,t},& \forall i\in PB,t\in T.\end{array}$$

(6)

Constraint (4) indicates that the heat source should be maintained for $T{S}_i^B$ hours when the heat source is activated at a specific time (t), and constraint (5) means that if the heat source is stopped, it should remain in the stopped state for $T{I}_i^B$ hours. Constraint (6) is a logical constraint that connects variables between the active state and the suspended state.

Type-C heat production facilities ($PC$) are facilities that produce electricity and heat at the same time. Combined heat and power production (CHP) and community energy systems (CES) are examples of type-C heat production. When it comes to cogeneration power plants in $PC$, constraints (2) to (6) can be modeled similarly to $PB$ for heat production. However, unlike heat-dedicated facilities, cogeneration facility operations also include electric power production. Therefore, constraints (2) and (3) should be amended as follows:

$$\begin{array}{cc}\hfill x_{i,t}\ge P_i^{MIN}{y}_{i,t},& \forall i\in PC,t\in T\end{array}$$

(7)

$$\begin{array}{cc}\hfill x_{i,t}\le P_i^{MAX}{y}_{i,t},& \forall i\in PC,t\in T,\end{array}$$

(8)

where $P_i^{MIN}$ and $P_i^{MAX}$ represent the minimum and maximum electricity production capacities respectively. In the above constraints, $x_{i,t}$ indicates the basic power output at time t of the cogeneration facility $i\in PC$. A typical cogeneration plant generates additional heat and power from a steam turbine connected to the gas turbine when the output is determined by a gas turbine. In steam turbines, a combination of equipment for producing electricity and equipment for producing heat is used. The ratio of heat and power production is determined by the proportion of high-temperature steam from the operation of the gas turbine into each facility in the steam turbine. Let $p_{i,t}$ and $h_{i,t}$ be the power and heat production quantity of heat source $i\in PC$ at time t, respectively. In practice, they are determined by the steam distribution ratio ($r_{i,t}$) in the steam turbine and the gas turbine power production $x_{i,t}$, and hence, they are indeed functions of $x_{i,t}$ and $r_{i,t}$ (i.e, $p_{i,t}=F(x_{i,t},r_{i,t})$ and $h_{i,t}=G(x_{i,t},1-r_{i,t})$). However, it is impossible to change the steam distribution ratio dynamically over time. Indeed, the cogeneration plant operates according to a predetermined operating mode. Typical operating modes are as follows:

• Operation Mode I (maximum heat operation): The maximum steam is distributed from the steam turbines to the heat production facility ($r_{i,t}=0$). Thus, electricity production is minimized, and this mode normally applies to seasons when demand for heating is high, such as during winter.
• Operation Mode III (maximum power operation): As much steam is distributed as possible from the steam turbines to the power production facility ($r_{i,t}=1$). It is mainly used during summer when heating demand is low and the sales unit price for electricity is high.
• Operation Mode V-P (power production priority mode): The steam turbines provide a balanced supply of steam to the heat and power production facilities, with a higher priority for the power production facility ($r_{i,t}=0.75$).
• Operation Mode V-H (heat production priority mode): The steam turbines supply steam to the heat and power production facilities equally, but the heat production facility is supplied first ($r_{i,t}=0.25$). The operating mode V is mainly applied to seasonal switches, such as spring and fall.

Given the operating mode (m), the heat and power outputs can be determined by the outputs from the gas turbines as follows:

$$\begin{array}{cc}\hfill p_{i,t}={\alpha }_{i,m}^p{x}_{i,m,t}^2+{\beta }_{i,m}^p{x}_{i,m,t}+{\gamma }_{i,m}^p{y}_{i,m,t},& \forall i\in PC,m\in M,t\in T\end{array}$$

(9)

$$\begin{array}{cc}\hfill h_{i,t}={\alpha }_{i,m}^H{x}_{i,m,t}^2+{\beta }_{i,m}^H{x}_{i,m,t}+{\gamma }_{i,m}^H{y}_{i,m,t},& \forall i\in PC,m\in M,t\in T\end{array}$$

(10)

$$\begin{array}{cc}\hfill x_{i,t}=\sum _{m\in M}{x}_{i,m,t},& \forall i\in PC,t\in T\end{array}$$

(11)

$$\begin{array}{cc}\hfill x_{i,m,t}\le P_i^{MAX}{y}_{i,m,t},& \forall i\in PC,m\in M,t\in T,\end{array}$$

(12)

where M is the set of operation modes (Modes I, II, V-P and V-H); $x_{i,m,t}$ is the electricity production at time t for cogeneration facility i with operation mode m; $y_{i,m,t}\in \{0,1\}$ is the operating state of facility i with operation mode m at time t; and $\alpha$,$\beta$ and $\gamma$ are weights (parameters) for the terms $x_{i,m,t}^2$, $x_{i,m,t}$ and $y_{i,m,t}$, respectively. Constraints (9) and (10) refer to power and heat production according to the gas turbine output, and efficiency factors are given for each operation mode. In general, the efficiency coefficient is modeled as a nonlinear function, which was assumed to be a quadratic function in this study. Constraints (11) and (12) are logical constraints that relate the gas turbine output at heat source i with mode m to the gas turbine output at heat source i.

There is also a minimum start-up time constraint in the operation mode. The constraints can be modeled in the same way as constraints (4)–(6) which implies that, if the specific operation mode m is determined, then the corresponding operation mode m should be used during the minimum operation time $T{S}_{i,m}^C$ for facility $i\in PC$.

$$\begin{array}{cc}\hfill y_{i,m,t+dt}\ge z_{i,m,t},& \forall i\in PC,m\in M,t\in T,dt\in \{0,1,\dots ,T{S}_{i,m}^C\}\end{array}$$

(13)

$$\begin{array}{cc}\hfill y_{i,m,t-1}+z_{i,m,t}=y_{i,m,t},& \forall i\in PC,m\in M,t\in T.\end{array}$$

(14)

In the constraints above, the decision variable $z_{i,m,t}$ is 1 if source i with mode m is started at time t and 0 otherwise. Constraint (13) implies that the corresponding operation mode should be active during $T{S}_{i,m}^C$ when the specific operation mode m is activated at time t. Constraint (14) is a logical constraint that associates variables $y_{i,m,t}$ and $z_{i,m,t}$.

3.2. Modeling of a District Heating Network

The district heating network consists of heat production nodes, heat demand nodes, and lines for connecting heat production nodes with demand nodes. We let N and E be the set of the nodes and edges, respectively. The network was modeled as flow balance constraints which are frequently used in network modeling. Let heat demand at time t of heat demand node k be $D_{k,t}$ and let the heat flow from node j to node k be $f_{j,k,t}$. In this case, the flow balance constraints can be expressed as follows:

$$\sum _{\{j\in N|(j,k)\in E\}}{f}_{j,k,t}+\sum _{i\in PS\left(k\right)}{h}_{i,t}+a_{k,t}=\sum _{\{j\in N|(k,j)\in E\}}{f}_{k,j,t}+D_{k,t}+b_{k,t},\forall k\in N,t\in T,$$

(15)

where $PS\left(k\right)$ is the set of production facilities in node k; $a_{k,t}$ indicates the heat that is dissipated in the storage tank k at time t; and $b_{k,t}$ is the heat stored in storage tank k at time t. In constraint (15), the left-hand side indicates hot water inflow and heat production at node k, and the right-hand side includes hot water outflow and heat demand at node k. Furthermore, $a_{k,t}$ and $b_{k,t}$ express the thermal balance according to the heat storage tank. Let the decision variable $r_{k,t}$ be the quantity of heat stored in the heat storage tank k at time t. Then, the following constraints can be derived:

$$\begin{array}{cc}\hfill r_{k,t-1}+b_{k,t}=r_{k,t}+a_{k,t},& \forall k\in N,t\in T\end{array}$$

(16)

$$\begin{array}{cc}\hfill a_{k,t}\le C{A}_k^{MAX},& \forall k\in N,t\in T\end{array}$$

(17)

$$\begin{array}{cc}\hfill b_{k,t}\le C{B}_k^{MAX},& \forall k\in N,t\in T\end{array}$$

(18)

$$\begin{array}{cc}\hfill r_{k,t}\le IN{V}_k^{MAX},& \forall k\in N,t\in T\end{array}$$

(19)

$$\begin{array}{cc}\hfill r_{k,t}\ge IN{V}_k^{MIN},& \forall k\in N,t\in T\end{array}$$

(20)

$$\begin{array}{cc}\hfill r_{k,T}\ge IN{V}_k^{TARGET},& \forall k\in N\end{array}$$

(21)

$$\begin{array}{cc}\hfill r_{k,0}=IN{V}_k^{INIT},& \forall k\in N.\end{array}$$

(22)

Constraint (16) is known as the inventory balance constraint. This is a relational expression for the heat stock stored at time t in the heat storage tank. Constraints (17) and (18) represent the maximum capacities for the amount of heat that is stored and the heat dissipated per hour, given that $C{A}_k^{MAX}$ and $C{B}_k^{MAX}$ represent the maximum allowable amounts of heat that are stored and dissipated per hour, respectively. Constraints (19) and (20) present the minimum and maximum capacity ($IN{V}_k^{MIN}$ and $IN{V}_k^{MAX}$, respectively) of the heat stock stored in the heat storage tank. Constraint (21) says that the inventory level in the storage tank at the end of the planning period must be above the target inventory level, $IN{V}_k^{TARGET}$, and constraint (22) specifies the initial inventory level of the heat stored in the storage tank at the beginning of the planning period, $IN{V}_k^{INIT}$.

3.3. Modeling the Objective Function

The objective function of this problem maximizes total profit over the planning horizon; the profit is calculated by subtracting the total costs for heat source operation from the revenue, which include power and heat sales. Unlike power sales and operating costs, heat sales are determined by multiplying the costs of the goods sold per calorie unit once the heat demand is determined. The operating costs differ depending on the heat source, but in the case of external heat sources (type-A facilities) and heat-exclusive equipments (type-B facilities), the operating costs are determined according to the unit production cost per calorie. The production costs of cogeneration plants (type-C facilities) are determined by the output of the gas turbine, regardless of the operation mode. Power sales were assumed to be determined by the market transaction price ($SM{P}_t$) in this study, although a complex method is applied depending on the trading and settlement characteristics in the electricity market. Therefore, the objective function can be expressed as follows:

$${\sum }_{k\in N}{\sum }_{t\in T}R·D_{k,t}+{\sum }_{i\in PC}{\sum }_{t\in T}SM{P}_t·p_{i,t}-\left\{{\sum }_{i\in PA\cup PB}{\sum }_{t\in T}{H}_i·h_{i,t}+{\sum }_{i\in PC}{\sum }_{t\in T}\left({\mu }_i{x}_{i,t}^3+{\eta }_i{x}_{i,t}^2+{\theta }_i{x}_{i,t}+{\kappa }_i{y}_{i,t}\right)\right\}.$$

(23)

The first term of the objective function indicates the heat sales according to the heating demand based on the heat sale unit price (R), and the second term presents the electricity sales calculated according to the unit price of the electricity sales ($SM{P}_t$) and the quantity of electric power production. The third term indicates the total production costs from the external heat sources and the heat-exclusive heat sources, where $H_t$ is the unit production cost for source $i\in PA\cup PB$. Finally, the last term represents the production costs according to the output from the cogeneration facilities. Generally, the cost of unit production per heat source is determined from the fuel consumed by the generators. The relationship between fuel input and gas turbine output in a cogeneration plant is generally a nonlinear function, and in this study, we specifically incorporated the cubic function ${\mu }_i{x}_{i,t}^3+{\eta }_i{x}_{i,t}^2+{\theta }_i{x}_{i,t}+{\kappa }_i{y}_{i,t}$ for the production costs of cogeneration facilities.

4. Proposed Algorithm

In this section, we propose a set partitioning-based algorithm to address the problem discussed in Section 3. The optimization problem for the integrated district heating operations discussed above is known to be very difficult to solve, because the size of the problem tends to be very large in practice and the model contains nonlinear functions. We note that many heuristic algorithms, including neighborhood search methods and the Lagrangian relaxation method, have been applied to the unit commitment problem for generators within a district heating network [5]. However, these approaches require a lot of time to derive a solution, and when a new constraint is added, the solution itself needs to be recalculated. Therefore, when they are applied to tackle real-world problems, a large amount of maintenance resources are required. One of the ways to address the difficulty is to relax the nonlinear function with a linear function; however, this method also has limitations in finding the solution when the problem size becomes large. To overcome the difficulties discussed above, we extended the explicit column generation algorithm, which has been successfully applied to the existing PBUC problems, to solve the concomitant problem.

The issue that increases the complexity of the problem is that we have a nonlinear function in the objective function and constraint functions for the cogeneration plants in $PC$. Therefore, we defined the heat and power generation patterns for this issue and solved the problem by selecting the optimal patterns after adding the operation patterns that reflect the local heating network operator’s know-how as input parameters for the optimization model. Operation pattern q of the cogeneration plant includes the operating mode and gas turbine output information for 24 h in the daytime operation pattern, and in this study, the network operation specialist, in practice, provides patterns satisfying the constraints (7)–(14), which are associated with type-C heat production facilities. For this, we introduced binary decision variable ${\varphi }_{i,q,d}$, which has a value of 1 when operation pattern q is selected on day d for heat source $i\in PC$, and it is 0 otherwise. Constraints (7)–(14) can then be simplified as follows:

$$\begin{array}{cc}\hfill p_{i,t}=\sum _{q\in Q}{P}_{i,q,Z\left(t\right)}{\varphi }_{i,q,D\left(t\right)},& \forall i\in PC,t\in T\end{array}$$

(24)

$$\begin{array}{cc}\hfill h_{i,t}=\sum _{q\in Q}{H}_{i,q,Z\left(t\right)}{\varphi }_{i,q,D\left(t\right)},& \forall i\in PC,t\in T,\end{array}$$

(25)

where Q is the set of operation patterns provided by the network operators; $Z\left(t\right)$ is a function for converting time to a specific time zone in $\{1,2,\dots ,24\}$ in 24 h; $D\left(t\right)$ is a function for converting time $t\in T$ to day $d\in D$, where D is the set of days, and $P_{i,q,tz}$ and $H_{i,q,tz}$ are the power and heat production of heat source $i\in PC$ at time zone $tz$ from operation pattern q, respectively. Because the operation pattern is a daily pattern, it is necessary to convert the time into days. Constraint (24) presents the power production at time t according to the selected operation patterns, and constraint (25) represents the heat production according to the operation pattern. Furthermore, because the selected operation pattern at day $d\in D$ for each heat source $i\in PC$ should be, at most, one, we added constraint (26), as shown below.

$$\sum _{q\in Q}{\varphi }_{i,q,d}\le 1,\hspace{1em}\forall i\in PC,d\in D.$$

(26)

Finally, we replaced the nonlinear function in the objective function (23) with a linear function as follows:

$$\sum _{k\in N}\sum _{t\in T}R·D_{k,t}+\sum _{i\in PC}\sum _{q\in Q}\sum _{d\in D}R{P}_{i,q,d}·p_{i,t}-\left(\sum _{i\in PA\cup PB}\sum _{t\in T}{H}_i·h_{i,t}+\sum _{i\in PC}\sum _{q\in Q}\sum _{d\in D}C{P}_{i,q,d}·{\varphi }_{i,q,d}\right).$$

(27)

In the updated objective function (27), $R{P}_{i,q,d}$ denotes power sales from source i on day d when operation pattern q is selected, and $C{P}_{i,q,d}$ represents the fuel costs for source i on day d when operation pattern q is selected. Given the operation pattern, the output and the operation mode are known for each time zone. Therefore, even if the nonlinear function is used, it is possible to pre-compute the operation pattern data in the pre-processing stage during the operation of the algorithm. Therefore, the problem of interest is transformed into the Set Partitioning Problem based on the operation pattern. The final formulation based on the set partitioning model is summarized as follows:

Model 

SetPartitioning

$$\begin{array}{cc}Minimize\hfill & \sum _{k\in N}\sum _{t\in T}R·D_{k,t}+\sum _{i\in PC}\sum _{q\in Q}\sum _{d\in D}R{P}_{i,q,d}·p_{i,t}-\left(\sum _{i\in PA\cup PB}\sum _{t\in T}{H}_i·h_{i,t}+\sum _{i\in PC}\sum _{q\in Q}\sum _{d\in D}C{P}_{i,q,d}·{\varphi }_{i,q,d}\right)\hfill \\ SubjectTo\hfill & h_{i,t}=S_{i,t},\hspace{1em}\forall i\in PA,t\in T\hfill \\ & h_{i,t}\ge C_i^{MIN}{y}_{i,t},\hspace{1em}\forall i\in PB,t\in T\hfill \\ & h_{i,t}\le C_i^{MAX}{y}_{i,t},\hspace{1em}\forall i\in PB,t\in T\hfill \\ & y_{i,t+dt}\ge z_{i,t},\hspace{1em}\forall i\in PB,t\in T,dt\in \{0,1,\dots ,T{S}_i^B\}\hfill \\ & y_{i,t+dt}\le w_{i,t},\hspace{1em}\forall i\in PB,t\in T,dt\in \{0,1,\dots ,T{I}_i^B\}\hfill \\ & y_{i,t-1}+z_{i,t}=y_{i,t}+w_{i,t},\hspace{1em}\forall i\in PB,t\in T\hfill \\ & p_{i,t}=\sum _{q\in Q}{P}_{i,q,Z\left(t\right)}{\varphi }_{i,q,D\left(t\right)},\hspace{1em}\forall i\in PC,t\in T\hfill \\ & h_{i,t}=\sum _{q\in Q}{H}_{i,q,Z\left(t\right)}{\varphi }_{i,q,D\left(t\right)},\hspace{1em}\forall i\in PC,t\in T\hfill \\ & \sum _{q\in Q}{\varphi }_{i,q,d}\le 1,\hspace{1em}\forall i\in PC,d\in D\hfill \\ & \sum _{\{j\in N|(j,k)\in E\}}{f}_{j,k,t}+\sum _{i\in PS\left(k\right)}{h}_{i,t}+a_{k,t}=\sum _{\{j\in N|(k,j)\in E\}}{f}_{k,j,t}+D_{k,t}+b_{k,t},\hspace{1em}\forall k\in N,t\in T\hfill \\ & r_{i,t-1}+b_{i,t}=r_{i,t}+a_{i,t},\hspace{1em}\forall k\in N,t\in T\hfill \\ & a_{i,t}\le C{A}_i^{MAX},\hspace{1em}\forall k\in N,t\in T\hfill \\ & b_{i,t}\le C{B}_i^{MAX},\hspace{1em}\forall k\in N,t\in T\hfill \\ & r_{i,t}\le IN{V}_i^{MAX},\hspace{1em}\forall k\in N,t\in T\hfill \\ & r_{i,t}\ge IN{V}_i^{MIN},\hspace{1em}\forall k\in N,t\in T\hfill \\ & r_{i,T}\ge IN{V}_i^{TARGET},\hspace{1em}\forall k\in N\hfill \\ & r_{i,0}=IN{V}_i^{INIT},\hspace{1em}\forall k\in N\hfill \\ & {\varphi }_{i,q,d},y_{i,t},z_{i,t},w_{i,t}\in \{0,1\}\hfill \\ & p_{i,t},h_{i,t},f_{j,k,t},r_{i,t},a_{i,t},b_{i,t}\in \mathbb{R}\hfill \end{array}$$

The SetPartitioning problem above can be solved by applying commercial optimization software, including CPLEX and Gurobi, even for large-scale problems in practice which is not possible for the earlier formulation presented in Section 3.

5. Case Study: Integrated Regional Heating Network Operation

In this section, the effectiveness of the proposed solution is verified based on a real-world case applied to an actual district heating network in South Korea, and the effect of integrated sustainable operation is analyzed. The data used in the tests is based on annual data and includes data from all of the regional offices of the Korea District Heating Corporation.

The constraints for some of the heat source facilities (constraints 4 and 5) were relaxed because the data needed to include time-specific optimization for 8760 h. FICO’s XpressMP 7.4 version was applied to find the optimal solution. For the case study, an optimization algorithm was implemented on a PC equipped with an Intel Sandbridge Quad Core i7 CPU (3.4 GHZ) and 4 GB of memory.

Table 1. Sample size for the case study.

Feature	Number
Heat demand nodes	24
Heat production nodes	38
Links (connecting two nodes)	3844
Heat production plants ($PA+PB$)	85
Heat cogeneration plants	9
Service coverage	1.5 million houses, 2.3 thousand buildings
Annual heat supply	13,013,000 Gcal
Annual electricity production	8,027,000 MWh

shows the size of the problem instance with 330,690 constraints, 361,936 variables, and 1,020,260 coefficients.

In order to express the benefits of integrated operation, we limited the total amount of interconnections between nodes and confirmed the change in the integrated profit, as presented in Figure 3 and

Table 2. Changes in integrated profit.

Changes in Capacity for Total Heat Flow between Nodes	Total Heat Flow between Nodes	Total Profits
0%	100	100
10%	81	90
20%	72	79
30%	63	68
40%	54	54

.

We set the algorithm to be terminated after derivation of the optimal solution within 3.5% within 150 s. Therefore, the optimal basic case was terminated by deriving a solution with a 1.95% difference from optimality in 150 s. The biggest difference between integrated operation and independent operation is how active heat exchange between the heat source nodes occurs. Therefore, we evaluated the change in integrated profit by restricting the maximally connectable calories using the following inequality, with the integrated profit and set at total calorific value for the solution derived from the proposed formulation in this study set at 100.

$$\sum _{(j,k)\in E}{f}_{j,k,t}\le TARGET\times CUT\_RATIO,\hspace{1em}\forall t\in T.$$

(28)

In constraint (28), $TARGET$ is the sum of the combined calories when there is no restriction on the combined calories, and $CUT\_RATIO$ is limited to the total calories in the combined case, which is limited to 10% to 40% in this case study. From the results, we can see that the integrated profit increased significantly with every 10% increase in the connected calories. In particular, when we switched from a 40% limit to a 30% limit, the combined profit increased by 14%, and when we subjected the combined calories from the 10% limit to the 0% limit, the combined profit increased by 10%. This means that the integration profit increases sharply when converted from independent operation to integrated operation. On an annual basis, it can be seen that the efficiency of the entire district heating network is increased by an increase in the integrated profit due to the integrated operation, in contrast to the existing independent operation. Since the Korea District Heating Corporation’s 2017 operating profit was KRW 161 billion, a 10% increase in profits is worth KRW 16 billion (or equivalently, USD 14.8 million).

To verify the effectiveness of the set partitioning approach, an optimal solution from the approach was compared to a restricted solution. When planning annual operations, the traditional approach of the company was to decrease the size of the problem by restricting the CHP operation modes. For example, CHP mode III and CHP mode I were disabled in planning for the winter and summer seasons, respectively, considering the heat demand in each season. Since heat demand during winter season is enormous, it is unnecessary to run CHPs in Mode III. By restricting the number of CHP modes and patterns in each month, planners can quickly consider alternative scenarios. They need to sacrifice solution quality a bit due to the reduced solution space arising from the restricted modes and patterns. In contrast, the proposed set partitioning approach quickly considers a vast set of operational patterns so as to maximize profit.

The solution comparison in Figure 4 shows that the overall solution gap between the optimal and restricted solution was 2.5%, which is about USD 18 million. The gap between the solutions was large during the summer season and smaller during the winter season. This is because, in winter, most efficient CHPs need to run at full workload, which implies that the problem is rather simple—just run the plants at full capacity. In the spring, summer, and fall seasons, the heat demand becomes lower, and the number of feasible planning options increases drastically. Therefore, the proposed approach is an improvement over the restricted solution during these low to mid-level demand seasons.

6. Discussion and Conclusions

In this study, we proposed a mathematical programming model and its solution in practice for a district heating network with combined heat and power plants. Unlike typical district heating systems, the Korea District Heating Corporation (KDHC)’s network is a large-scale single district heating network connecting dozens of sources of heating demand and heating supply nodes dispersed along major cities in Korea. This case study is the first in which optimization has been attempted by switching from the existing independent operating model to an integrated operating model. The model and solution algorithm proposed in this study have, in fact, been applied to the KDHC which has installed an integrated operation center in Pankyo to integrate and operate the entire district heating network of the company. By incorporating the proposed algorithm in this paper, the operation efficiency of the KDHC district heating network greatly increased, and it was possible to construct a flexible system that works in practice.

The optimal operation of district heating systems with combined heat and power plants has been extensively studied in the past. Most recent approaches have applied mixed integer programming to model the complicated operation of a diverse range of heat sources and the interconnectivity of production and demand nodes. To deal with nonlinear objective functions and constraints, piecewise linear approximation techniques have been used, but the size of the formulations used in previous research increased drastically, which means they are impractical for real-world applications that involve longer planning periods (e.g., annually). In addition, the detailed modeling of mixed integer programming approaches tends to generate black box solutions for planners because the formulations consist of a huge set of logical constraints which is sensitive to outside parameters. To manage the complexity of the problem and to generate practically acceptable solutions, we derived operation patterns from CHPs, which encapsulated the detailed operational characteristics of CHPs into a single operation pattern. This encapsulation successfully removed some of the constraints from the formulation which contributed to the enhanced efficiency of the proposed set partitioning algorithm.

The proposed model and subsequent computational analysis have the following implications.

• First, by applying the proposed solution algorithm to the annual planning of the case study company, complicated problems such as monthly or annual planning can be solved efficiently. Previous mixed-integer programming models have required millions of variables and constraints even for shorter weekly planning scenarios, but the proposed pattern-based algorithm encapsulates many complex constraints into the patterns, which contributes to efficient optimization. The proposed solution algorithm can quickly solve daily or weekly planning problems.
• Second, when we tightly interconnected the distant production and demand nodes, the overall profit showed the potential to increase a lot. When we increased the link capacity of the network, the system-wide profits increased accordingly. When there are multiple competing CHPs and thermal storages connected in a single district heating network, the proposed algorithm can effectively balance heat production, transfer, and storage decisions. Integratative operations of district heating networks can be managed using the proposed solution algorithm.
• Third, the annual planning application can be used to identify any bottleneck links which may restrict flows between nodes. If we ensure greater capacity in those bottleneck links, it is expected that the productivity of the network can be improved. Network-level profits may be negatively impacted by just a few bottleneck links. We need to focus on this small number of critical links which constrain the network connectivity.
• Fourth, the case study demonstrated that, when we consider more patterns, the overall profit can increase accordingly. More patterns allow more realistic scenarios. The comparison in the case study showed that the optimization algorithm is particularly useful in lower-demand seasons because, in higher-demand seasons, the only available option is to just run the CHPs at full capacity. However, during seasons like spring and fall, we need to carefully compare the available options and scenarios.

In this paper, we focused on introducing a practical adaptive solution to nonlinear functions and profit maximization for large-scale network problems. It is expected that the algorithm can be further improved in order to reduce the time required to derive the solution.