 ORIGINAL ARTICLE
 Open Access
Linear programming control of a group of heat pumps
 Jiří Fink^{2}Email author,
 Richard P. van Leeuwen^{2, 1},
 Johann L. Hurink^{2} and
 Gerard J.M. Smit^{2}
https://doi.org/10.1186/s1370501500619
© Fink et al. 2015
 Received: 23 February 2015
 Accepted: 14 October 2015
 Published: 17 November 2015
Abstract
Background
For a new district in the Dutch city Meppel, a hybrid energy concept is developed based on biogas cogeneration. The generated electricity is used to power domestic heat pumps which supply thermal energy for domestic hot water and space heating demand of households. In this paper, we investigate direct control of the heat pumps by the utility and how the largescale optimization problem that is created can be reduced significantly.
Methods
Two different linear programming control methods (global MILP and time scale MILP) are presented. The latter solves largescale optimization problems in considerably less computational time. For simulation purposes, data of household thermal demand is obtained from prediction models developed for this research. The control methods are compared with a reference control method resembling PI on/off control of each heat pump.
Results
The reference control results in a dynamic electricity consumption with many peak loads on the network, which indicates a high level of simultaneous running heat pumps at those times. Both methods of mix integer linear programming (MILP) control of the heat pumps lead to a much improved, almost flat electricity consumption profile.
Conclusions
Both optimization control methods are equally able to minimize the maximum peak consumption of electric power by the heat pumps, but the time scale MILP method requires much less computational effort. Future work is dedicated on further development of optimized control of the heat pumps and the central CHP.
Keywords
 Smart grid
 Linear programming control
 Heat pumps
 Thermal storage
 Load balancing
 Demand side management
Background
Introduction

reduce costs for integration of renewable energy into distribution grids

increase consumer awareness and energy savings

contribute to a competitive energy market and consumer choice resulting in reductions in energy prices
It is widely recognized that smart grids are required to balance energy production with loads and energy storage. This is made possible by sophisticated control strategies and communication networks.
To increase the share of renewable energy in today’s energy system, it is not always necessary to integrate the renewable energy into existing grids. With socalled smart microgrids, it is possible to connect or disconnect a region to or from the main distribution grid. In Meppel, a small city in the Netherlands, there is a plan to build a district (Nieuwveenselanden) with an energy system called MeppelEnergie [4–6] that can be defined as a smart microgrid and which has the goal to become almost 100 % renewable. The heart of the system is a combined heat and power (CHP) unit which is supplied with biogas from a nearby wastewater treatment facility. Part of the houses in the district are connected to the district heating system, while other houses have a heat pump. Both the heat for the district heating system and power for the heat pumps are supplied by the CHP. This is a typical case where a smart grid is required for load balancing purposes, otherwise the district would need frequent power exchanges with the existing power grid. Without a smart grid, this would probably involve large peak loads and strengthening of the existing grid would be required, involving high costs for new cabling and bidirectional transformers. For different cases in Germany involving regional solar PV generation, these effects are explained by Nykamp [2]. Therefore, the heat pumps should be scheduled in such a way that they only consume, if possible, the electricity produced by the CHP. If this is not possible, the remaining energy has to be bought on the electricity market at minimal cost.
The planning of a group of heating systems may have many objectives in practice. In the Meppel project, biogas is transported through a dedicated pipeline and electric energy is distributed through a private cable and converted by heating supply systems. In general, generators and transport equipment have to be dimensioned for the maximal peak consumption. Thus, the main objective is minimizing the maximal consumption which may decrease investments in the system and will also lead to maximum self consumption of the generated electricity by the biogas CHP.
This paper develops a smart grid control method with the goal to minimize peak consumption. We further demonstrate the importance of this method for the Meppel case. After a review of related work in section “Related work”, a more detailed problem formulation which leads to an algorithm called global mix integer linear programming (MILP) control is given in section “Problem statement and global MILP control”. As this algorithm requires a lot of computational power, we develop an algorithm called time scale MILP control in section “Time scale MILP control”. The Meppel case is explained in section “Case application”. The simulation results are presented in section “Results and discussion”. Finally, we draw conclusions in section “Conclusions”.
Related work
The problem of minimizing peak consumption of a group of heat pumps is a typical smart grid control problem. As demonstrated by Nykamp [2], integration of renewable energy (e.g., by wind turbines, solar PV, and biomass conversion) requires implementation of control strategies for load balancing and load shifting in electrical grids. This is either due to the stochastic nature of the energy production or due to supply limitations of the biomass. In the Meppel case, the supply of biogas to the CHP which generates the electricity for the heat pumps is limited. Besides that, the biogas supply is approximately constant in time and it is not possible to divert the biogas supply to other consumers. These are totally different circumstances than existing natural gas CHP systems. The natural gas supply system is widespread within the Netherlands, and for CHP’s on natural gas, there are usually no limitations to size or consumption patterns.

1. Autonomous control is performed by the device itself, e.g., by sensing the grid frequency, interpreting the difference with a nominal frequency, and making decisions to start or stop the heat pump. In our case, we expect poor results from this type of control because there are also constraints on thermal comfort which are likely to create conflicts in time with reaching grid stability.

2. Pricebased control is the dominating control strategy for scheduling devices within a smart grid. The heat pump determines an operational schedule based on dynamic price information which is sent to the heat pumps by the utility. However, to make a sensible decision about future schedules, it is required that a heat pump controller is able to predict future heat demand of the building. For this purpose, model predictive control is implemented within the heat pump controller. For climate control of a single building with a heat pump and timevarying electricity prices, model predictive control is developed and demonstrated by Halvgaard et al. [8]. For groups of devices, a smart grid control strategy based on the principles of pricebased control is developed within the TRIANA smart grid control methodology by Molderink [9], Bakker [10], and Bosman [11]. Bosman demonstrates the use of this type of control for a group of microCHP’s in [11], while Bakker also demonstrates results for a group of heat pumps in [10]. These examples demonstrate effectiveness of pricebased control. The only minor disadvantage of pricebased control is the need for an intelligent heat pump controller for each house.

3. Direct control involves making control decisions for each heat pump on the utility level. This has the advantage that the heat pump controllers can be simple. On the other hand, the utility has to solve largescale optimization problems (in our case, minimizing peak electricity consumption) if there are a large number of heat pumps to control. Besides that, the optimization has to include heat demand predictions for each house at every time step of the calculation. In the Meppel case, the number of houses with a heat pump connected to a single biogas CHP is limited to 135. For mathematical optimization, we estimate that this is a feasible challenge and that is why we chose this method for the Meppel case.
The mathematical background of the minimizing peak problem is presented in [12] which proves that the problem of minimizing peak is NPcomplete. A somewhat similar problem was considered by Bosman et al. [13, 14] who studied a microCHP planning problem and also proved that minimizing peak is NPcomplete in their model [15]. Bosman et al. [11] also present a dynamic programming algorithm for the microCHP planning problem whose time complexity is O(T ^{3C+1}) where T is the number of time intervals and C is the number of microCHPs.
The main contribution of this paper is that we demonstrate by development of the time scale control method and the application within the Meppel case that the number of variables and equations which are part of the largescale optimization problem on the utility level can be reduced considerably while the control objective of minimizing the peak consumption is still reached.
Methods
Problem statement and global MILP control

a supply which represents some source of energy (e.g., electricity, gas),

a converter which converts the energy into heat (hot water),

a buffer which stores heat for later usage, and

a demand which represents the consumption profile of heat.
In principle, the presented model can consider arbitrary types of energy but in this paper, we use electricity and heat to distinguish consumed and produced energy. This simple model of a local heating system cannot only be applied for heating water but has many other applications in smart grids (e.g., control of fridges and freezers) and inventory management. We come back on this at the end of this section.
We consider a discrete time model for the considered problem, meaning that we split the planning period into T time intervals of the same length. We consider a set \({\mathcal {C}} = \{1,\dots,C\}\) of C heating systems and a set \({\mathcal {T}} = \{1,\dots,T\}\) of T time intervals. Note that the heating of a house is split into two independent heating systems (see Fig. 1). In this paper, the letter c is always an index of a heating system (either space heating or tap) and t is an index of a time interval. For mathematical purposes, we separate a heating system into a converter, a buffer, and demand; see Fig. 1. We say “a converter c,” “a buffer c,” or “a demand c” to refer to the devices of the heating system c∈.
We consider a simple converter which has only two states: In every time interval, the converter is either turned on or turned off. The amount of consumed electricity is E _{ c } and the amount of produced heat (or any other form of energy) is H _{ c } during one time interval in which the converter \(c \in \mathcal {C}\) is turned on. If the converter is turned off, then it consumes and produces no energy. Let x _{ c,t }∈{0,1} be the variable indicating whether the converter \(c \in \mathcal {C}\) is running in time interval \(t \in \mathcal {T}\).
The state of charge of a buffer \(c \in \mathcal {C}\) at the beginning of time interval \(t \in \mathcal {T}\) is denoted by s _{ c,t } which represents the amount of heat in the buffer. Note that s _{ c,T+1} is the state of charge at the end of planning period. The state of charge s _{ c,t } is limited by an upper bound U _{ c }.
The amount of consumed heat by the inhabitants of the house from heating system \(c \in \mathcal {C}\) during time interval \(t \in \mathcal {T}\) is denoted by D _{ c,t }. This amount is assumed to be given and is called the demand of heating system c. In this paper, we study offline problems, so we assume that demands D _{ c,t } are given for the whole planning period.
Equation 1 is the charging equation of the buffer. During time interval \(t \in \mathcal {T}\), the state of charge s _{ c,t } of a buffer \(c \in \mathcal {C}\) is increased by the production of the converter which is H _{ c } x _{ c,t } and it is decreased by the demand D _{ c,t }. Equations (2) and (3) ensure that the domains of variables s _{ c,t } and x _{ c,t }, respectively, are taken into account. Note that the initial state of charge s _{ c,1} can be fixed (e.g., by setting \(s_{c,1} = \frac {U_{c}}{2}\)).
Since E _{ c } x _{ c,t } is the amount of consumed electricity by a converter c in time interval t, the sum \(\sum _{c \in \mathcal {C}} E_{c} x_{\text {\textit {c,t}}}\) is the amount of electricity consumed by all converters in time interval t. Furthermore, the inequality and the objective function (4) guarantees that the value of the variable m is the maximal consumption of electricity during one time period within the whole planning period.
In the following, we give some other possible applications of this model.
Fridges and freezers: A fridge essentially works in the opposite way than heating, so it may be modelled similarly. However, we have to be careful with the correct interpretation of all parameters. The state of charge of the buffer again represents the temperature inside the fridge, but a higher state of charge means a lower temperature. The converter does not produce heat to the fridge but it decreases the temperature inside the fridge, so the converter increases the state of charge of the buffer (fridge). The demand decreases the state of charge of the fridge due to thermal loss and usage of the fridge by humans.
Inventory: The considered heating problem is also related to inventory control problems [16]. A buffer may represent an inventory and a converter may represent orders. However, this leads to a situation, where only a limited capacity of inventory is given and it is only possible to order a fix amount of goods which is not a typical situation in inventory management.
Note that the objective function and all constrains are linear and operational state variables are binary, so constraints and the objective (1)–(4) form an instance of MILP. This instance can be solved by any MILP solver (see, e.g., [17]) and we call this approach global MILP control. However, as the number of binary variables may get too large for planning many houses over a long planning horizon, this method may get computationally expensive. Therefore, in the following sections, an algorithm which significantly reduces the number of variables is given.
Time scale MILP control
The method presented in the previous section creates one large instance of MILP and solves it by a MILP solver. This method gives us an optimal solution for the whole planning period but it may not be suitable for practical purposes. First, finding an optimal solution requires a lot of computational power. Next, the prediction of demand for the distant future may be very inaccurate.
Therefore, we consider an online control in which the decision which converters will be running is made only for the coming time interval. On the other hand, we cannot ignore the future completely. Indeed, we should take more care about the near future time intervals than the distant ones because the current decision has stronger impact on the near future and the prediction is in general more accurate for the near future.
As the general formal notation of time scale control may be hard to understand, we use an example to present the approach. Assume that the upcoming time interval has index 1. The decision which converters \(c \in \mathcal {C}\) will be running during the coming time interval 1 needs to be made, meaning that the values of variables x _{ c,1} need to be decided. Since also the influence of the very near future needs to be detailed, we also consider binary variables, e.g., for the next two time intervals, i.e., variables x _{ c,2} and x _{ c,3}. For the further future, the plan does not need to be so precise, so,e.g., another two time intervals, we relax the integral constraints, meaning that we require 0≤x _{ c,4},x _{ c,5}≤1. The reason for relaxing these variables is to decrease the number of integral variables which has the dominating effect on the computational time required to solve an MILP instance. From the practical point of view, these relaxed variables (e.g., x _{ c,4}) can signify the probability that a converter c will run in time interval 4, and so \(\sum _{c \in \mathcal {C}} H_{c} x_{c,4}\) is the expected demand of electricity.
where D _{ c,8..10}=D _{ c,8}+D _{ c,9}+D _{ c,10} is the cumulative demand for time intervals 8,9 and 10. In this example, we replace the three variables x _{ c,8}, x _{ c,9}, and x _{ c,10} by one aggregated variable x _{ c,8..10} which is constrained by bounds 0≤x _{ c,8..10}≤3. In this way, we can cover a longer planning horizon without requiring too many variables in an MILP instance. Furthermore, note that only the sum of demands D _{ c,8..10} for distant future time intervals is important. The practical consequence is that the time where a significant demand occurs does not have to be predicted precisely, e.g., in the morning, it is sufficient to predict the amount of hot water demand for evening showers but the exact time when inhabitants will take a shower can be approximated.
for every \(c \in \mathcal {C}\).
This instance of MILP problem can be solved by any MILP solver. The values of variables x _{ c,1} of an optimal solution are used to determine which converters \(c \in \mathcal {C}\) should run in the coming time interval 1. For the next time interval 2, a similar instance of an MILP problem is created by shifting indices of time intervals by one and refining the predicted demands D _{ c,t }.
There is no general rule how time intervals should be split into blocks since it is strongly influenced by the particular case studies. In this study, we use one specific choice to study the potential of this approach.
Case application

develop a thermal model to determine house space heating demand

apply measured weather data (we investigate 1 week with a high space heating demand)

define various typical house and household profiles to generate a variety of space heating and domestic hot water demand profiles

simulate space heating demand of the various households

determine reference PI control results

input domestic hot water and space heating demands into the global and time scale MILP control algorithms and determine results
In the next sections, the individual steps of the approach are explained in more details.
Space heating thermal model
Nomenclature energy system characterization
Term  Signification 

T _{ a }  Ambient temperature 
T _{ z }, C _{ z }  Zone (room) temperature and thermal capacity 
T _{wo}, C _{wo}  Outside wall temperature and thermal capacity 
T _{wi}, C _{wi}  Inside wall temperature and thermal capacity 
T _{ f }, C _{ f }  Zone floor temperature and thermal capacity 
T _{cs}  Cellar or creeping space temperature 
R _{roof}  Thermal resistance of roof 
R _{window}  Thermal resistance of windows 
R _{wo−a }  Thermal resistance between outside wall and ambient 
R _{wiwo}  Thermal resistance between outside and inside wall 
R _{ z−wi}  Thermal resistance between zone and inside wall 
R _{ z−f }  Thermal resistance between zone and floor 
R _{ f−cs}  Thermal resistance between floor and concrete structure 
q _{vent}  Ventilation heat flow 
q _{inf}  Infiltration heat flow 
q _{gain}  Internal gain heat flow 
D _{ c,t }  Heating load flow (heating Demand) 
ϕ _{air,inf}  Infiltration air flow [ m ^{3}/h] 
ϕ _{air,vent}  Ventilation air flow [ m ^{3}/h] 
ρ _{air}  Air density 
c _{ p,air}  Specific heat capacity of air 
The used terms are also explained in Table 1. Ventilation air flow values ϕ _{air,vent} in m ^{3}/h are defined by schedules given in section “House and household case information”. Infiltration air flows ϕ _{air,inf} are related to leakages of the building and are assumed as constant values. We defined some variations including heat recovery ventilation with a certain heat recovery efficiency γ _{hr}.
Besides heat losses, there are scheduled heat gains (q _{gain}) due to resident and electric appliance dissipations. The applied schedule is given in section “House and household case information”. For the present paper, solar energy absorbed by the interior is neglected because we study a relatively cold and cloudy week.

If the zone temperature T _{ z } equals the setpoint, then the demand D _{ c,t } is the amount of energy which the heating system has to generate to keep the zone temperature constant.

If the setpoint is increased, then the demand D _{ c,t } is increased to raise the zone temperature in a given warmup speed \(\frac {dT_{z}}{\text {dt}}\).

If the zone temperature T _{ z } is above the setpoint (e.g., due to decreasing the setpoint or natural heating by internal gains), the demand D _{ c,t } is the minimal nonnegative amount of energy which keeps the zone temperature above the setpoint.
Application of weather data
As we explained in the previous section, we limited our investigation to heat loss due to temperature differences between the zone temperature and ambient temperatures. Effects of solar gains and wind speeds on the heat demand are excluded. Weather data containing hourly average ambient temperatures are obtained from the website of the Dutch national weather institute ([20]). We choose data of weather station Hoogeveen which is close to Meppel. We investigate the coldest week of 2012 as this week results in a relatively high heat demand for space heating.
House and household case information
Types of households, the number of household types in both types of houses, schedules of higher temperature setting, and hot water demands
Type of household  Young couple  Young family  Elderly people  

Number of persons in a household  2  4  2  
Number of houses  semidetached  27  54  9 
detached  12  27  6  
Higher setpoint  weekdays  17–22  8–22  10–23 
weekend  9–23  9–23  10–23  
Hot water on weekdays  morning  15 MJ  8 MJ  4 MJ 
afternoon  0 MJ  4 MJ  4 MJ  
evening  20 MJ  24 MJ  20 MJ  
Hot water on weekend  morning  8 MJ  4 MJ  4 MJ 
afternoon  4 MJ  4 MJ  4 MJ  
evening  24 MJ  32 MJ  24 MJ 
We define a lower setpoint for the zone temperature (18 °C) during working hours and night and a higher setpoint temperature (20 °C) otherwise. We also define domestic hot water demand for morning, afternoon, and evening peaks. See Table 2 for schedules of the temperature setpoints and energy demands for hot water.
Simulation results of heating demand
The model equations given in section “Space heating thermal model” are solved using a 15min time step which is required as a minimum time step for controlling the thermal storage and to calculate heat pump operation times. To obtain 15min weather data, we applied linear interpolation of the 1h data values.
We assume a heat pump coefficient of performance (COP) of 4.5 for space heating and 2.5 for domestic hot water generation during the simulation. With that, we obtain as average electricity input 116.4 kW for space heating and 19.3 kW for domestic hot water, leading to 135.7 kW total average electricity demand. If the MILP control performs well on minimizing peaks, we expect control schedules for the heat pumps to give results close to this average electricity demand.
Results and discussion
As seen on Fig. 3, both methods of MILP control of the heat pumps lead to a much improved, almost flat electricity consumption profile. It is also demonstrated that the result of time scale MILP control is very close to global MILP control. The average for both is 135.6 kW. The peak reduction compared to reference control is 98 % for global MILP control and 96 % for time scale MILP control. If we consider this almost equal performance and take into account the much reduced computational effort of time scale MILP control, we prefer this method for future algorithm development.
The energy system described in section “Introduction” contains a central CHP. The obtained flat electricity demand profile of the heat pumps may be generated entirely by the central CHP but this depends on the possibility to simultaneously store or use the produced heat by the CHP in a connected district heating system. Part of future work is to study possible CHP control schemes more closely, focusing on stable CHP operation, optimal sustainability, and optimal profitability in relation to varying external grid prices.
Conclusions
In this paper, we investigate how well a group of heat pumps for 135 different households can be controlled to minimize peaks within the electricity distribution network by two methods of direct control. The significance of this research is that this type of control enables integration of renewable energy within a microelectricity grid. Without this control, power consumption peaks would occur frequently and this results into higher investments needed for network strengthening.
From simulations with a simplified thermal network model of the houses, we obtained heat demand data for space heating. For the simulations, we defined occupancy schedules and different insulation properties of the houses. We also defined domestic hot water demand profiles for three types of households. We compare the total electricity consumption of all the heat pumps for the following control methods: reference (PI) control where each heat pump determines itself when to switch on or off, and two forms of direct control: global MILP control and time scale MILP control. The latter has the advantage of a much reduced scale for solving the optimization problem.
In the defined case of 135 houses, MILP control decreases electricity peaks by 96 % compared to reference control. The difference between global and time scale MILP control is small. As time scale MILP control is computationally much more efficient, we propose to use time scale MILP control. The influence of the chosen time scaling on computational time and the quality of the achieved result is part of future work. Although the results look promising, due to the high level of peak reduction, future work will be dedicated to investigate the resulting thermal comfort as a result of the obtained heat pump planning. Also part of future work is to investigate methods of reaching social fairness in heat pump planning and integrating this work into the TRIANA smart grid control method [21] to investigate possible CHP control schemes with a focus on stability, sustainability, and profitability in relation to varying electricity grid prices.
Declarations
Acknowledgements
The authors would like to thank the Dutch national program TKISwitch2SmartGrids for supporting the project Meppelenergy and the STW organization for supporting the project iCare 11854.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Authors’ Affiliations
References
 Meer Duurzame Energie in de Toekomst. visited October 2014. http://www.rijksoverheid.nl/onderwerpen/duurzameenergie/meerduurzameenergieindetoekomst. Accessed 26 Oct 2015.
 Nykamp S (2012) Integrating renewables in distribution grids—storage, regulation and the interaction of different stakeholders in future grids. PhD thesis, University of Twente (2013). doi:10.3990/1.9789036500579.
 Innovation Contract Smart Grids. visited October 2014. http://topsectorenergie.nl/wpcontent/uploads/2013/10/InnovatieContractSmartGrids2012.pdf. Accessed 26 Oct 2015.
 Nieuwveense Landen. visited March 2014. http://www.meppelwoont.nl/nieuwveenselanden/. Accessed 26 Oct 2015.
 RVO Meppel heats new housing development with biogas. visited March 2014 (2012) http://www.rvo.nl/sites/default/files/Meppel%20heats%20new%20housing%20development%20with%20biogas.pdf. Accessed 26 Oct 2015.
 Meppel Energie. visited March 2014. http://www.meppelenergie.nl/nieuwveenselanden. Accessed 26 Oct 2015.
 Stadler M, Krause W, Sonnenschein M, Vogel U (2009) Modelling and evaluation of control schemes for enhancing load shift of electricity demand for cooling devices. Environ Modelling & Softw 24(2): 285–295.View ArticleGoogle Scholar
 Halvgaard R, Poulsen NK, Madsen H, Jorgensen JBEconomic model predictive control for building climate control in a smart grid. Proceedings of 2012 IEEE PES innovative Smart Grid Technologies (ISGT).Google Scholar
 Molderink A (2011) On the three step methodology for smart grids. PhD thesis, University of Twente. 10.3990/1.9789036531702.
 Bakker V (2011) Triana: A control strategy for smart grids: forecasting, planning and realtime control. PhD thesis, University of Twente. 10.3990/1.9789036533140.
 Bosman M (2012) Planning in smart grids. PhD thesis, University of Twente. 10.3990/1.9789036533867.
 Fink J, Hurink JL (2015) Minimizing costs is easier than minimizing peaks when supplying the heat demand of a group of houses. Eur J Oper Res 242: 644–650.View ArticleGoogle Scholar
 Bosman MGC, Bakker V, Molderink A, Hurink JL, Smit GJM (2009) The microCHP scheduling problem In: Proceedings of the Second Global Conference on Power Control and Optimization, PCO 2009, Bali, Indonesia. AIP Conference Proceedings, 268–275.Google Scholar
 Bosman MGC, Bakker V, Molderink A, Hurink JL, Smit GJM (2011) Controlling a group of microCHPs: planning and realization In: First International Conference on Smart Grids, Green Communications and IT Energyaware Technologies, ENERGY 2011, Venice, Italy, 179–184.Google Scholar
 Bosman MGC, Bakker V, Molderink A, Hurink JL, Smit GJM (2010) On the microCHP scheduling problem In: Proceedings of the 3rd Global Conference on Power Control and Optimization PCO, 2010, Gold Coast, Australia. AIP Conference Proceedings, 367–374.Google Scholar
 Axsäter S (2006) Inventory control. International series in operations research and management Science, Vol. 90. Springer, New York.Google Scholar
 Cook WJ, Cunningham WH, Pulleyblank WR, Schrijver A (2011) Combinatorial optimization, Vol. 33. Wiley. com, Rhode Island.Google Scholar
 van Leeuwen RP, Fink J, de Wit JB, Smit GJ (2015) Upscaling a district heating system based on biogas cogeneration and heat pumps. Energy, Sustainability Soc 5(1): 16. 10.1186/s137050150044x.View ArticleGoogle Scholar
 Mitchel JW, Braun JE (2013) Principles of heating, ventilation, and air conditioning in buildings. Wiley, New York.Google Scholar
 KNMIweather records station Hoogeveen. http://www.knmi.nl. Accessed 26 Oct 2015.
 Toersche HA, Bakker V, Molderink A, Nykamp S, Hurink JL, Smit GJM (2012) Controlling the heating mode of heat pumps with the TRIANA three step methodology In: Innovative Smart Grid Technologies (ISGT), 2012 IEEE PES, 1–7. http://ieeexplore.ieee.org/iel5/6170475/6175527/06175662.pdf.