1 Introduction

The integration of renewable energy sources (RES) in the power system and the economic and environmental benefits they offer are not coming without any problems. The continuous variations of these sources cause stability and balancing problems in the power system. Power generation units should be scheduled in advance to provide reserves in the system for the case of mismatch between the scheduled production and consumption, increasing in that way the total generation cost. As the integration of RES will be increasing, covering higher percentage of the demand, the need for balancing power will also be increasing. This is the reason why a lot of research is being conducted in order to find the optimal allocation of reserves in power systems [1].

Units that provide balancing services are characterized by quick start up times and high ramp rates. Hydro and gas turbines are two typical types of units that can offer such services. Combined heat and power (CHP) units can also be possibly used for providing reserves in the system. These units usually operate in district heating (DH) systems providing heat to the network and generating power that can be sold to the electricity market or used to cover some internal load of the system. There are also heat accumulators, heat boilers and in some cases heat pumps which make the system quite flexible to power and heat load changes. Furthermore, the efficiency of the CHP units is very high and this is the reason why the use of DH is promoted in Europe and other countries. According to [2], by the year 2050, DH should cover 50% of the total heat demand in contrast to 10% as per today.

In order to estimate the capability of a DH system to provide reserves and follow the variations of the load and the RES generation, a model that simulates the hour by hour operation of such system is needed. Such short term operation scheduling models involve the decisions about which units are going to operate each hour and cover the heat and power demand at minimum cost. In the case of DH systems, the coupling of power and heat production in the CHP units increases the complexity of this scheduling problem. The first studies in the literature used non linear deterministic models of the CHP operation [3, 4] without considering any RES generation. In order to tackle with the non linearities, various algorithms have been proposed to solve the problem, among them particles warm optimization [5], artificial immune systems [6] and differential evolution [7].

The introduction of competitive electricity markets over the last decades, has brought back the interest on CHP operation modeling. The new objective has been the maximization of profits in the electricity markets. In [8], for example, an algorithm is suggested to solve the optimal hourly CHP operation as a function of the electricity price and in [9] a deterministic model is suggested to minimize the CHP operation costs when selling the power to the day ahead market. In [10] both day ahead and balancing markets are considered for a portfolio of CHP and wind power while in [11] a framework for trading CHP generation in day ahead, intraday and reserves markets is proposed. Furthermore, new modeling techniques that help make optimal decisions under uncertainty have found application in the operation scheduling of power systems in general, and in DH systems in particular. For example, dynamic stochastic programming is used in [12] to determine the optimal bidding strategy of a CHP operator in various electricity markets and robust optimization in [13] to optimally schedule the operation of CHP units in the day ahead market. In [14] stochastic programming is used to estimate the value of the electric boilers and heat pumps in a DH system as it is done also in [15], while in [16] a multi-objective chance constrained model is proposed to minimize the CHP operation costs, the total amount of primary energy consumption and the \(\hbox {CO}_2\) emissions. Stochastic programming is also used in [17] where a multi-stage medium-term model for CHP operation is suggested considering some risk measure. [18, 19] provide a deep review of the models and techniques that have been used so far for CHP and DH operation scheduling.

Regarding the integration of RES into the power system there is some recent research that investigates the capability of CHP/DH systems to provide load following and reserves services. [20] uses a commercial software to model a CHP system which balances the variable wind generation while in [21] a model is proposed for optimizing the operation of a DH system with a solar thermal plant. In [22], a deterministic model is used to determine the provision of control power by a DH system. A unit commitment and power/heat dispatch model for CHP operation is used (but not presented) in [23] for wind power integration studies. [24] suggest that RES integration should be seen as part of a broader smart energy system. Electricity-to-heat conversion systems for wind power integration is considered in [25]. A linear dispatch model is used in [26] to estimate the flexibility that CHP units, heat storage and electric boilers are offering for wind power integration in China. Finally, in [27] the potential of power-to-heat technology in Germany is discussed.

In this paper we model the short-term operation of a DH system that provides both power and heat to the demand while it tries to maximize the wind power penetration in the power system by avoiding the wind power curtailment. The DH system is fully modeled with CHP units, heat only boilers, electric boilers, heat pumps and heat storage. The proposed model is mixed integer linear and the solution provides the optimal unit commitment and power/heat dispatch with hourly resolution. In brief, the purpose of this study is:

  1. 1)

    To provide a short term unit commitment and energy (both power and heat) dispatch model for a DH system.

  2. 2)

    To verify the load following capabilities in hourly and sub-hourly resolution of a DH system for various levels of wind power penetration in the power system.

  3. 3)

    To estimate the value of the heat pumps/electric boilers and the heat storage.

The main contribution of this paper includes the proposal of an optimization model that incorporates all major units of a modern DH system. The majority of the works presented above consider only the CHP units with or without heat storage. Furthermore, the objective varies, and the combination of DH systems with RES finds application in recent only works. Regarding the modeling part, it should be added that the back pressure steam turbine is modeled to operate in two modes which has not been presented before. Another contribution regards the case study where the unit parameter data and the heat load are coming from a real DH system. The same applies for the power load and the wind generation, allowing us to obtain realistic results and make better conclusions regarding the integration of wind power with the help of DH systems. Finally, one more contribution is the simulation of the DH operation with 15-min period data in the case study and the comparison with the hourly operation. Such study has not been conducted in the previous works.

The rest of this paper is organized as follows. Section 2 presents the operational model of a district heating system as explained above. In Section 3, a case study is conducted to show the application of the model in a modified real system. Finally, Section 4 concludes the paper and suggests further modifications to the suggested model.

2 Model description

Before mathematically formulating the model in Section 2.3, a brief description of the DH main components and their operation is given in Section 2.1. This is done for better understanding of the model equations. Furthermore, in Section 2.2, the assumptions made for the model are stated.

2.1 District heating system brief description

A DH system which consists of CHP units, heat only boilers, electric boilers, heat pumps and heat accumulators is shown in Fig. 1. The CHP units are the main heat and power generation units. There are three main types of turbines used with different characteristics: back-pressure steam turbines, extraction condensing steam turbines and gas turbines. The back-pressure steam turbines operate with fixed power-to-heat ratio which results in less flexibility. The feasible operation zone is represented by an inclined line shown in Fig. 2b. In some installations, there is the possibility for the steam to bypass the turbine and operate in heat boiler mode. The gas turbines also operate with fixed power-to-heat ratio. In this case they can also operate as only power units. The extraction condensing steam turbines are characterized for the high flexibility they offer. They can operate in any point of the feasible operation zone shown in Fig. 2a. The heat only boiler burn some type of fuel to produce heat. They are usually used as back up of the CHP units. The electric boilers or immersion heaters and the heat pumps are consuming power to produce heat. The difference is that the heat pumps need a heat source like e.g. a lake and they are much more efficient compared to the electric boilers. Finally the heat accumulators can store hot water for many hours and with very low losses. This gives the possibility to the DH operator, to decouple the power and the heat production.

Fig. 1
figure 1

Representation of district heating system

Fig. 2
figure 2

Power-heat load charts of back pressure and extraction condensing steam turbine

2.2 Assumptions

The following assumption are made:

  1. 1)

    The power flow in the grid is not modeled. We assume that there are not any power flow limits or line losses. This assumption is not strong if a local system is considered which is typical for DH systems. However, if a wider system should be modeled, it is easy to incorporate the network constraints into the the model. In that case, if the size of the system is quite big, the problem may be difficult to solve needing some special decomposition techniques to be applied [28].

  2. 2)

    The DH network is not also modeled. It is just represented with a simple heat balance equation. This assumption is also not strong as the heat losses are indirectly considered as part of the heat load.

  3. 3)

    The power and heat loads and the wind power generation are considered to be perfectly forecasted.

2.3 Formulation

2.3.1 Objective function

The objective function (1) stands for the minimization of the aggregated operational costs over all time periods which are denoted by t and over all units which are denoted by g. These include the cost of fuel consumption for power and heat generation and the start up and shut down costs of the units. The objective function also includes the cost of wind power curtailment and the costs of power and heat load shedding. These costs are incorporated into the objective function in order to avoid as much as possible these conditions.

$$\begin{aligned}&\underset{t=1}{\overset{N_{T}}{\sum }}\Biggl [\underset{g=1}{\overset{N_{G}}{\sum }}\left( \lambda _{g}^{fuel}P_{g,t}^{fuel}+c_{g}^{stop}Z_{g,t}\right. \nonumber \\&\quad +\,c_{g}^{start,h}Y_{g,t}^{h}+c_{g}^{start,w}Y_{g,t}^{w}+c_{g}^{start,c}Y_{g,t}^{c}\Bigr )\nonumber \\&\quad +\,c^{w}W_{t}^{spil}+c^{l}P_{t}^{shed}+c^{h}Q_{t}^{shed}\Biggr ] \end{aligned}$$
(1)

where \(P_{g,t}^{fuel}\) is the fuel consumption; \(Z_{g,t}\) is a binary variable which is equal to 1 when a unit is shutting down; \(Y_{g,t}^{h}\), \(Y_{g,t}^{w}\) , \(Y_{g,t}^{c}\) are binary variables which are equal to 1 if there is a hot, warm or cold start up of the unit, respectively; \(W_{t}^{spil}\) is the wind power that is curtailed; \(P_{t}^{shed}\) is the power load not served and \(Q_{t}^{shed}\) is the heat load not served. Parameters \(\lambda _{g}^{fuel}\), \(c_{g}^{stop}\), \(c_{g}^{start,h}\), \(c_{g}^{start,w}\), \(c_{g}^{start,c}\), \(c^{w}\), \(c^{l}\) , \(c^{h}\) are referring to the costs of fuel, shut down, start up (hot, warm and cold), wind power curtailment, power load shedding and heat load shedding respectively.

2.3.2 Unit technical limit constraints and fuel consumption functions

  1. 1)

    Extraction condensing steam turbine

    $$\begin{aligned}&\beta _{g}^{el}P_{g,t}+\beta _{g}^{th}Q_{g,t}\leqslant \beta _{g}^{el}p_{g}^{\mathrm {max}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (2)
    $$\begin{aligned}&\beta _{g}^{el}P_{g,t}+\beta _{g}^{th}Q_{g,t}\geqslant \beta _{g}^{el}p_{g}^{\mathrm {min}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (3)
    $$\begin{aligned}&P_{g,t}\geqslant \alpha _{g}Q_{g,t},\quad \forall g,t \end{aligned}$$
    (4)
    $$\begin{aligned}&Q_{g,t}\leqslant q_{g}^{\mathrm {max}},\quad \forall g,t \end{aligned}$$
    (5)
    $$\begin{aligned}&P_{g,t}^{fuel}=\frac{1}{\eta _{g}}\left( \beta _{g}^{el}P_{g,t}+\beta _{g}^{th}Q_{g,t}\right) ,\quad \forall g,t \end{aligned}$$
    (6)

    The feasible operational zone of an extraction condensing unit as shown in Fig. 2a is defined by the set of constraints (2)–(5), where \(P_{g,t}\) and \(Q_{g,t}\) are the power and heat generation; \(U_{g,t}\) is a binary variable that defines when the unit is operating. The marginal fuel consumption for power \(\beta _{g}^{el}\) and heat \(\beta _{g}^{th}\) production, the power generation upper \(p_{g}^{\mathrm {max}}\) and lower \(p_{g}^{\mathrm {min}}\) limits, the heat generation upper limit \(q_{g}^{\mathrm {max}}\) and the power-to-heat ratio \(\alpha _{g}\) are technical parameters of the units. The fuel consumption is given by (6) where {\(\eta _{g}\)} is the efficiency of the unit.

  2. 2)

    Back pressure steam turbine

    $$\begin{aligned}&P_{g,t}=\alpha _{g}Q_{g,t}^{m1},\quad \forall g,t \end{aligned}$$
    (7)
    $$\begin{aligned}&p_{g}^{\mathrm {min}}M_{g,t}^{1}\leqslant P_{g,t}\leqslant p_{g}^{\mathrm {max}}M_{g,t}^{1},\quad \forall g,t \end{aligned}$$
    (8)
    $$\begin{aligned}&q_{g}^{\mathrm {min}}M_{g,t}^{2}\leqslant Q_{g,t}^{m2}\leqslant q_{g}^{\mathrm {max}}M_{g,t}^{2},\quad \forall g,t \end{aligned}$$
    (9)
    $$\begin{aligned}&Q_{g,t}=Q_{g,t}^{m1}+Q_{g,t}^{m2},\quad \forall g,t \end{aligned}$$
    (10)
    $$\begin{aligned}&\underset{\tau =t}{\overset{t+t_{g}^{dn}-1}{\sum }}\left( M_{g,\tau }^{1}+M_{g,t}^{2}-1\right) \leqslant 0,\quad \forall g,t \end{aligned}$$
    (11)
    $$\begin{aligned}&U_{g,t}=M_{g,t}^{1}+M_{g,t}^{2},\quad \forall g,t \end{aligned}$$
    (12)
    $$\begin{aligned}&P_{g,t}^{fuel}=\frac{1}{\eta _{g}}\left( P_{g,t}+Q_{g,t}\right),\quad \forall g,t \end{aligned}$$
    (13)

    A back pressure unit can operate in two modes, as CHP and as heat only boiler. When operating as CHP, the power generation is connected to the heat production through a fixed power-to-heat ratio as defined by (7) where \(Q_{g,t}^{m1}\) is the heat generation of the unit in CHP mode. The generation limits are given by (8) where \(M_{g,t}^{1}\) is a binary variable that defines when the unit is operating in CHP mode. Because of the fixed power-to-heat ratio there is no need to define limits for the heat generation. When operating in heat boiler mode, the heat production is limited by its lower \(q_{g}^{\mathrm {min}}\) and upper \(q_{g}^{\mathrm {max}}\) limits (9). \(M_{g,t}^{2}\) is a binary variable that defines when the unit is operating in heat boiler mode and \(Q_{g,t}^{m2}\) is the heat generation of the unit in heat boiler mode. The total heat production is equal to the sum of productions in the two modes (10), however (11) restricts the operation of the unit in only one mode each period. Furthermore, when the unit is operating in the heat boiler mode, it needs some time until it can switch back to the CHP mode. This is not applied for the opposite direction, from CHP to heat boiler mode. This time can be considered equal to the minimum time \(t_{g}^{dn}\) that the unit should be kept off-line after a shut down although in this case the boiler is still operating and the plant is not totally switched off. Therefore the switch back to the CHP mode is a hot start up as it is explained below. This condition is also satisfied by (11). Constraint (12) is used to assign the operation status (on/off) of the unit. When the unit operates either in CHP or heat boiler mode is considered to be on. The fuel consumption is given by (13).

  3. 3)

    Gas turbine

    $$\begin{aligned}&P_{g,t}\geqslant \alpha _{g}Q_{g,t},\quad \forall g,t \end{aligned}$$
    (14)
    $$\begin{aligned}&p_{g}^{\mathrm {min}}U_{g,t}\leqslant P_{g,t}\leqslant p_{g}^{\mathrm {max}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (15)
    $$\begin{aligned}&P_{g,t}^{fuel}=\frac{1}{\eta _{g}}\cdot \frac{\alpha _{g}+1}{\alpha _{g}}P_{g,t},\quad \forall g,t \end{aligned}$$
    (16)

    The technical limits for a gas turbine are similar to the back pressure unit (14)–(15). Here an inequality can be used to express the relation between power generation and heat production (14). With this, it is possible to model the case when the heat is released to the environment through a cooling system instead of the DH network. By fixing the heat production variable to zero, the gas turbine is generating only power. The fuel consumption is given by (16).

  4. 4)

    Heat only boiler

    $$\begin{aligned}&q_{g}^{\mathrm {min}}U_{g,t}\leqslant Q_{g,t}\leqslant q_{g}^{\mathrm {max}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (17)
    $$\begin{aligned}&P_{g,t}^{fuel}=\frac{Q_{g,t}}{\eta _{g}},\quad \forall g,t \end{aligned}$$
    (18)

    The output of a heat boiler is restricted by its lower and upper limits (17) . The fuel consumption is equal to the heat production divided by the efficiency of the boiler (18).

  5. 5)

    Electric boiler

    $$\begin{aligned}&Q_{g,t}\leqslant q_{g}^{\mathrm {max}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (19)
    $$\begin{aligned}&Q_{g,t}=\eta _{g}P_{g,t}^{con},\quad \forall g,t \end{aligned}$$
    (20)

    The output of an electric boiler is only restricted by its upper limit (19) while the heat production is equal to the power consumption \(P_{g,t}^{con}\) multiplied by the efficiency of the unit (20).

  6. 6)

    Heat pump

    $$\begin{aligned}&q_{g}^{\mathrm {min}}U_{g,t}\leqslant Q_{g,t}\leqslant q_{g}^{\mathrm {max}}U_{g,t},\quad \forall g,t \end{aligned}$$
    (21)
    $$\begin{aligned}&Q_{g,t}=cop_{g}P_{g,t}^{con},\quad \forall g,t \end{aligned}$$
    (22)

    Similarly, the output of a heat pump is restricted by its lower and upper limits (21) while the heat production is equal to the power consumption multiplied by the coefficient of performance (COP) \(cop_{g}\) of the unit (22).

  7. 7)

    Heat accumulator

    $$\begin{aligned}&v^{\mathrm {min}}\leqslant V_{t}\leqslant v^{\mathrm {max}},\quad \forall t \end{aligned}$$
    (23)
    $$\begin{aligned}&-q^{flow}\leqslant V_{t}-V_{t-1}\leqslant q^{flow},\quad \forall t \end{aligned}$$
    (24)

    The heat content \(V_{t}\) of the accumulator is restricted up by the capacity \(v^{\mathrm {max}}\) of the accumulator and down by a limit \(v^{\mathrm {min}}\) that is applied for safety reasons (23). Furthermore there is an upper limit \(q^{flow}\) of the heat that flows into and out of the accumulator (24).

2.3.3 Energy balance constraints

$$\begin{aligned}&V_{t}=\left( 1-f\right) V_{t-1}+\underset{g=1}{\overset{N_{G}}{\sum }}Q_{g,t}-q_{t}^{load}+Q_{t}^{shed},\quad \forall t \end{aligned}$$
(25)
$$\begin{aligned}&P_{t}^{wind}=w_{t}-W_{t}^{spil},\quad \forall t \end{aligned}$$
(26)
$$\begin{aligned}&\underset{g=1}{\overset{N_{G}}{\sum }}P_{g,t}+P_{t}^{wind}=\underset{g=1}{\overset{N_{G}}{\sum }}P_{g,t}^{con}+P_{t}^{self}+p_{t}^{load}-P_{t}^{shed},\quad \forall t \end{aligned}$$
(27)

Constraint (25) is used to define the heat balance of the system. According to this, the accumulator’s heat content at the end of a period t is equal to the heat content at the end of the previous period \(t-1\) multiplied by a loss factor \(1-f\), plus the difference between the total heat production, the heat load shedding and the heat consumption \(q_{t}^{load}\). When there is no storage, this constraint simply sets the total heat production plus any heat shedding equal to the heat load. Heat losses of the storage, which are defined with parameter f, are typically very low. The heat load shedding variable is used to avoid any conditions that will render the model infeasible. Constraint (26) is just fixing the wind power generation variable \(P_{t}^{wind}\) to the forecast \(w_{t}\) minus any possible wind power curtailment. Finally, constraint (27) defines the power balance of the system. According to this, the total power generation of the units plus the wind power generation should be equal to total consumption of the units (electric boilers and heat pumps), plus the internal consumption of the system \(P_{t}^{self}\) and the power load \(p_{t}^{load}\), minus any possible load shedding. Any power consumption that is not part of the load can be assigned to variable \(P_{t}^{self}\), like e.g. the pumping power in the DH network.

2.3.4 Unit commitment constraints

$$\begin{aligned}&{\left\{ \begin{array}{ll} Y_{g,t}\leqslant U_{g,t}\\ Y_{g,t}\leqslant 1-U_{g,t-1}\\ Y_{g,t}\geqslant U_{g,t}-U_{g,t-1} \end{array}\right. }\quad \forall g,t \end{aligned}$$
(28)
$$\begin{aligned}&{\left\{ \begin{array}{ll} Z_{g,t}\leqslant U_{g,t-1}\\ Z_{g,t}\leqslant 1-U_{g,t}\\ Z_{g,t}\geqslant U_{g,t-1}-U_{g,t} \end{array}\right. } \quad \forall g,t \end{aligned}$$
(29)

The set of constraints (28) and (29) assign the correct values to the binary variables that define start up \(Y_{g,t}\) and shut down \(Z_{g,t}\) status of the units. When a unit is starting up variable \(Y_{g,t}\) takes a value of 1. Similarly, when a unit is shutting down, variable \(Z_{g,t}\) becomes 1. This form of the constraints allows for the relaxation of variables \(Y_{g,t}\) and \(Z_{g,t}\) which can be set as positive variables instead of binary, as they can only take the values of 0 or 1.

2.3.5 Ramping constraints

$$\begin{aligned}&-t^{per}r_{g}^{p\downarrow }\leqslant \begin{array}{c} P_{g,t}-P_{g,t-1}\leqslant t^{per}r_{g}^{p\uparrow }\end{array},\quad \forall g,t \end{aligned}$$
(30)
$$\begin{aligned}&-r_{g}^{q\downarrow }\leqslant \begin{array}{c} Q_{g,t}-Q_{g,t-1}\leqslant r_{g}^{q\uparrow }\end{array} \end{aligned}$$
(31)

Constraints (30)–(31) define the ramp up and ramp down limits for both the power and heat generation. \(t^{per}\) is the length of the period in minutes. \(r_{g}^{p\uparrow }\)and \(r_{g}^{p\downarrow }\) are the upward and downward ramp rates for power generation given in MW/min. \(r_{g}^{q\uparrow }\) and \(r_{g}^{q\downarrow }\) are the upward and downward ramp rates for heat generation given in \(\hbox {MW}_{\mathrm{th}}\)/hour.

2.3.6 Start up type constraints

$$\begin{aligned}&Y_{g,t}^{h}\leqslant \underset{\tau =t-T_{g}^{w}+1}{\overset{t}{\sum }}Z_{g,\tau },\quad \forall g,t \end{aligned}$$
(32)
$$\begin{aligned}&Y_{g,t}^{w}\leqslant \underset{\tau =t-T_{g}^{c}+1}{\overset{t-T_{g}^{w}}{\sum }}Z_{g,\tau },\quad \forall g,t \end{aligned}$$
(33)
$$\begin{aligned}&Y_{g,t}=Y_{g,t}^{h}+Y_{g,t}^{w}+Y_{g,t}^{c},\quad \forall g,t \end{aligned}$$
(34)

The set of constraints (32) and (33) are used to select the correct type of the unit start up which depends on the prior reservation time of the unit. If a unit starts up less than \(T_{g}^{w}\) hours after it has been shut down, then this is considered a hot start up and binary variable \(Y_{g,t}^{h}\) takes a value of 1. If it starts up after \(T_{g}^{w}\) hours but less than \(T_{g}^{c}\), then this is considered a warm start up and binary variable \(Y_{g,t}^{w}\) takes a value of 1. Finally, if the start up takes place more than \(T_{g}^{c}\) hours after the unit has been shut down there is a cold start up. This last case does not need to be modeled separately as constraint (34) satisfies that if there is a start up it will be one of hot, warm or cold type. With every type of start up a different cost is associated, being less for a hot start up and more for a cold start up. This is briefly summarized in Table 1 where \(T^{off}\) is the prior reservation time in hours.

Table 1 Unit start up type modeling

2.3.7 Minimum up and down time constraints

$$\begin{aligned}&\underset{t=1}{\overset{l_{g}}{\sum }}\left( 1-U_{g,t}\right) =0,\quad \forall g \end{aligned}$$
(35)
$$\begin{aligned}&\underset{\tau =t}{\overset{t+t_{g}^{up}-1}{\sum }}\left( U_{g,\tau }-Y_{g,t}\right) \geqslant 0,\nonumber \\&\quad \forall g,\,t=l_{g}+1,...,N_{T}-t_{g}^{up}+1 \end{aligned}$$
(36)
$$\begin{aligned}&\underset{\tau =t}{\overset{N_{T}}{\sum }}\left( U_{g,\tau }-Y_{g,t}\right) \geqslant 0,\nonumber \\&\quad \forall g,\,t=N_{T}-t_{g}^{up}+2,...,N_{T} \end{aligned}$$
(37)
$$\begin{aligned}&\underset{t=1}{\overset{f_{g}}{\sum }}U_{g,t}=0,\quad \forall g \end{aligned}$$
(38)
$$\begin{aligned}&\underset{\tau =t}{\overset{t+t_{g}^{dn}-1}{\sum }}\left( 1-U_{g,\tau }-Z_{g,t}\right) \geqslant 0,\nonumber \\& \quad \forall g,\,t=f_{g}+1,...,N_{T}-t_{g}^{dn}+1 \end{aligned}$$
(39)
$$\begin{aligned}&\underset{\tau =t}{\overset{N_{T}}{\sum }}\left( 1-U_{g,\tau }-Z_{g,t}\right) \geqslant 0,\nonumber \\&\quad \forall g,\,t=N_{T}-t_{g}^{dn}+2,...,N_{T} \end{aligned}$$
(40)

Constraints (35)–(40) define minimum up and down times of the units. \(t_{g}^{up}\) is the minimum time in hours that a unit should operate after its starting up and similarly \(t_{g}^{dn}\) is the minimum time in hours that a unit should be kept off-line after it has been shut down. In (36), \(l_{g}=\mathrm {min}\left\{ N_{T},\left( t_{g}^{up}-t_{g}^{up,0}\right) u_{g}^{0}\right\}\) are the hours in the beginning of the planning horizon that the unit is restricted to operate while in (39), \(f_{g}=\mathrm {min}\bigl \{ N_{T},\left( t_{g}^{dn}-t_{g}^{dn,0}\right) \left( 1-u_{g}^{0}\right) \bigr \}\) the hours that the unit is restricted to be off-line due to the initial conditions. \(N_{T}\) is the cardinality of the period set T, \(u_{g}^{0}\) the initial operating status of unit g, and \(t_{g}^{up,0}/t_{g}^{dn,0}\) the hours that unit g was on-line/off-line before the planning period.

Table 2 Case study system parameters
Fig. 3
figure 3

Annual load and heat demand profile

3 Case study

The purpose of the case study is first to verify the applicability of the proposed model and second to derive some conclusions regarding the integration of wind power with the help of DH systems. For this reason a realistic DH system is modeled and various scenarios are tested. The system configuration is described in Section 3.1 while the results are presented and analyzed in Section 3.2.

3.1 System configuration

A DH system is considered which consists of 32 units. The type of each unit and their parameters are given in Table 2. In this table, BP stands for back pressure, EC for extraction condensing, GT for gas turbine, HB for heat boiler, EB for electric boiler and HP for heat pump. The symbol “-” means that the specific parameter is not applicable to that type of unit. The \(33^{\mathrm{rd}}\) unit is a hypothetical heat boiler used in some scenarios as it is explained below. There is also a heat storage with a capacity ranging from 0 to 1500 \(\hbox {MWh}_{\mathrm{th}}\) according to the scenario.

Data of load and heat demand with hourly resolution are used. These data are coming from real power and DH networks but are scaled down to 70% of the system’s power capacity and 60% of the system’s (of 32 units) heat capacity respectively. The annual profile of these two parameters is depicted with different scales in Fig. 3. As it can be seen, heat demand is ranging from approximately 120 \(\hbox {MWh}_{\mathrm{th}}\) in summer up to 1200 \(\hbox {MWh}_{\mathrm{th}}\) in winter. The load profile on the other hand is not presenting so big fluctuations among seasons. It ranges between approximately 280 and 420 MWh.

Wind power generation data with 15-min resolution are also used from an operating wind park. Initially, the wind power series are averaged in order to get a series with hourly resolution and then they are normalized between 0 and 1. Then four new wind power series are created by up-scaling the normalized series. The available annual wind generation of these new series is equal to 25%, 50%, 75% or 100% of the yearly load. This is done in order to test various scenarios of wind power penetration. The cases of 100% and 50% of the annual load are depicted in Fig. 4. As it can be seen, the annual wind power generation is equal to the annual load when the installed wind power capacity is 4-5 times the maximum load. For the 50% case, this should be about 2 times the maximum load.

Fig. 4
figure 4

Annual load and wind power profile

3.2 Results

Various scenarios are tested. We start with a annual simulation of the DH operation under various wind power penetration scenarios and then we compare the results with that of a DH system without electric boilers or heat pumps. We also compare the operation costs for DH systems with different heat storage capacities. Finally, wind power generation data with 15-min resolution are used in a new simulation and the results are compared with the results of the initial simulation.

3.2.1 Simulation of DH daily operation

The daily operation of the DH system consisted of units 1–32 for a period of one year is simulated. This is done for the four wind power generation series which were analyzed before plus for the case of no wind power generation. The system has no heat storage capacity. The cost of wind power curtailment is set to 1000 €/MWh, while the costs of power load and heat load shedding are set to 2000 and 1000 €/MWh respectively. The internal power consumption of the DH system is mostly driven by the need for water pumping and is set to be 0.5% of the heat generation (\(P_{t}^{self}=0.005Q_{g,t}\)) according to [29]. Table 3 presents the main results of the annual simulation for the different wind power generation scenarios. The annual operation cost is decreasing with the increase of the wind power capacity till a level of 30%–40%. After that level, although the wind power generation keeps increasing, the cost of wind power curtailment is getting very high affecting the total operation cost. The level where the cost starts increasing depends on the choice of the cost of wind power curtailment. A lower cost will increase that level but also will increase the wind power curtailment. Power generation from the CHP units and the gas turbine is decreasing with the increase of the wind power level. The power consumption from the electric boilers and the heat pumps follows the opposite direction while the heat generation is fixed. Regarding the power and heat load shedding, they are negligible in all cases.

The weakly operation of the DH system for the first week of the year is shown in Fig. 5 for the power generation and in Fig. 6 for the heat generation and for levels 0%, 50% and 100%. Regarding the load, this is mostly covered by the CHP units and the wind power. When wind power is available, the heat pumps and the electric boilers start operating to produce heat. This can be easily seen in hours 1–40 for the 50% or 100% wind power level. Wind power generation is high, exceeding the load. This extra power is used by the heat pumps and the electric boilers to produce heat. This results in the shut down of the CHP units for some hours during which the extra heat needed is produced by boilers.

Table 3 Aggregated results of DH system annual operation
Fig. 5
figure 5

Power generation and load for first week of year

Fig. 6
figure 6

Heat generation and heat demand for first week of year

A question regarding back pressure CHP units is if it profitable to operate them in boiler mode instead for example just shutting them down and let some only heat boiler to produce that heat. Figure 7 shows the annual heat production from the back pressure CHP units (both in CHP and boiler mode) and the heat only boilers. The CHP operates most of the time in CHP mode as it is expected. The production of heat in boiler mode is much smaller but in the same order of magnitude with the heat only boilers. What is also interesting is the fact that with the increase of wind power, the production of heat in CHP mode is reducing while that in boiler mode is increasing.

Fig. 7
figure 7

Annual heat production from back pressure CHP units and boilers

3.2.2 Value of heat storage

In this case a heat storage is added to the DH system. The technical parameters of the heat storage are given as: \(f=0.02\); \(v^{\mathrm {min}}=200\); \(v^{\mathrm {max}}=500\); \(q^{flow}=400\). The annual operation costs compared to the costs in the previous case are shown in Fig. 8. The reduction of the operation costs is negligible. To test if the storage capacity affects this, the weekly operation of a DH system with various levels of heat storage is simulated. Capacities are ranging from 0 MWh to 1500 \(\hbox {MWh}_{\mathrm{th}}\) . Results are presented in Fig. 9. It seams that the heat storage is not improving the operation of the system. Only for very high levels of wind power there is a slight reduction of the operation cost. However it should be noted that as the simulation was done in a daily basis, the use of the heat storage content is not optimal due to the short planning period. A longer planning period (weekly) should be used for optimal use of the heat storage content.

Fig. 8
figure 8

Comparison of annual operation cost for DH system with and without heat storage

Fig. 9
figure 9

Weakly cost of energy production for various heat storage capacities

3.2.3 Value of electric boilers and heat pumps

In this case the electric boilers and the heat pumps in the initial DH system are replaced with a heat only boiler. Its technical parameters are given in Table 2 (unit #33). Its capacity is equal to the total capacity of the heat pumps and the electric boilers. A yearly simulation is run as before. Figure 10 presents the operation costs, the wind power generation and the wind power curtailment for both systems and for all wind power generation levels. Till the level of 25%, the results are similar. But as the percentage of wind power is increasing, the DH system which uses electric boilers and heat pumps is much more flexible, having reduced operation cost and allowing higher penetration of wind power.

Fig. 10
figure 10

Comparison of annual operation cost, wind power generation and wind power curtailment of DH system with heat pumps/electric heaters and DH system with heat boilers

3.2.4 Fifteen minute resolution simulation

As wind power generation may present high variations in sub-hourly scale, it is important to check if the results from the hourly simulation are still valid under such wind power variations. For this reason a new simulation is run with 15-min periods. The planning horizon of the day consists now of 96 periods. As it has been referred, wind power generation data have 15-min resolution. Figure 11 presents the available wind power scaled to cover up to 100% of the load for the first week of the year. Both the 15-min and 1-hour resolution series of wind power generation are given.

Load and heat demand data which have hourly resolution are expanded in a way that each previous hourly period consists of four new periods with the same value. This is a simplification but as both load and heat demand do not present very frequent variations, this simplification is not affecting the results so much. Special care should be given to the system parameters like costs and minimum times because the first are measured in €/MWh while the second in hours. Therefore all costs which are in per MWh basis should be divided by four, while minimum times should be multiplied by four. Similarly, upward and downward ramp rates for heat generation should be divided by four while the period length should be set to 15 minutes. Annual simulation results are given in Table 4. A comparison between these and the results in Table 3 shows that the differences are very small. This can also be seen in Fig. 12 which presents the weekly power generation with 15-min resolution. A comparison with Fig. 5 shows that the differences in power generation are small. The reason for this is the fact that CHP units have high ramp rates, typically 10–30 MW/min. The same applies for the electric boilers and the heat pumps. However, as seen in Fig. 11, wind power generation variations are usually below 100 MW/quarter in this case study. Therefore, a model with hourly period may be used for analyzing wind power integration with the help of DH systems.

Fig. 11
figure 11

Available wind power for period of one week with 1-hour and 15-min resolution

Table 4 Annual results for simulation with 15-min resolution
Fig. 12
figure 12

Power generation schedule with 15-min period

4 Conclusion

This paper presents a short term operation model for a DH system which can also be used for wind power (and other RES) integration studies. It has been shown with a case study that simulates the daily operation of a realistic system for a period of one year, that units of a DH system have the capability to follow the load and wind power generation variations even for very high wind power penetration although the cost is very high in this case. The use of electric boilers, heat pumps and heat storage can reduce the generation costs and increase the wind power generation. Furthermore, it has been shown that a model with hourly resolution provides good results and there is not need for higher resolution modeling.

The presented model is deterministic and assumes perfect forecast of the load, the heat demand and the wind power generation. In order to take uncertainty into account, this model can be extended into a stochastic one where the uncertain parameters will be described by a set of scenarios with their respective probabilities.