Linear programming control of a group of heat pumps

Abstract

Background

For a new district in the Dutch city Meppel, a hybrid energy concept is developed based on bio-gas co-generation. 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 large-scale 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 large-scale 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.

Background

Introduction

To reduce fossil fuel usage and increase the share of renewable energy, the Netherlands plans for 14 % renewable energy of the total energy consumption by the year 2020 and 100 % in the year 2050 [1]. However, integration of renewable energy in existing distribution grids may cause problems with power stability in some regions at times when there is a large amount of energy offered by solar PV installations or wind turbines and low demand [2]. Smart grids are seen as a possible solution and therefore the Dutch government started the Topteam Energy which defined the innovation contract smart grids [3] with the following goals:

• 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 so-called 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 bi-directional 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.

For load balancing purposes, Stadler et al. [7] introduce three possible device control strategies for refrigerators which are also applicable for heat pumps.

• 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. Price-based 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 time-varying 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 price-based 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 price-based control. The only minor disadvantage of price-based 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 large-scale 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 NP-complete. 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 NP-complete 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 large-scale 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

In this section, we present a mathematical description of the studied model and possible applications of this model. The model is set up to be generally applicable for any type of heating appliance. Typical appliances for heating water are electrical and gas heating systems, heat pumps, and combined heat and power units (microCHP). The heated water is stored in buffers to be prepared for the demand of inhabitants. In our model, a house consists of two local heating systems, one for space heating and the other for tap. A schematic overview of the model is presented in Fig. 1. It consists of:

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

  Fig. 1

  

  Schematic picture of a house with two separated heating systems for space heating and tap

• 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.

The operational variables of the converters x _c,t and the states of charge of buffers s _c,t are restricted by the following constraints.

$$\begin{array}{@{}rcl@{}} \!\!\!\!\!\!\!s_{\text{\textit{c,t}}+1} = s_{\text{\textit{c,t}}} + H_{c} x_{\text{\textit{c,t}}} - D_{\text{\textit{c,t}}} &\! \text{for} & \!c \in {\mathcal{C}},\, t \in {\mathcal{T}} \end{array} $$

((1))

$$\begin{array}{@{}rcl@{}} 0 \!\le s_{\text{\textit{c,t}}} \!\le U_{c} & \!\text{for} & \!c \in {\mathcal{C}},\, t \in \!\{1,\dots,T\!+1\} \end{array} $$

((2))

$$\begin{array}{@{}rcl@{}} x_{\text{\textit{c,t}}} \in \{0,1\} & \text{for} & \!c \in {\mathcal{C}},\, t \in {\mathcal{T}} \end{array} $$

((3))

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}$).

In this paper, we consider the objective function of minimizing the peak:

$$ \begin{array}{rl} \text{minimize} &m \\ \text{where} &m \ge \sum_{c \in \mathcal{C}} E_{c} x_{\text{\textit{c,t}}} \ \ \text{for} \ \ t \in \mathcal{T} \end{array} $$

((4))

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.

Following this, for the even more distant future, we only need a rough planning. In order to explain the idea of rough planning, let us consider the state of charge equation, e.g., for time intervals t=8,9 and 10.

$$\begin{array}{@{}rcl@{}} s_{c,9} &=& s_{c,8} + H_{c} x_{c,8} - D_{c,8} \\ s_{c,10} &=& s_{c,9} + H_{c} x_{c,9} - D_{c,9} \\ s_{c,11} &=& s_{c,10} + H_{c} x_{c,10} - D_{c,10} \end{array} $$

((5))

We sum these equations and after simplification, we obtain

$$ s_{c,11} = s_{c,8} + H_{c} (x_{c,8} + x_{c,9} + x_{c,10}) - (D_{c,8} + D_{c,9} + D_{c,10}) $$

((6))

The rough plan for converter c for time intervals t, t+1,…,t ^′ is now defined by $x_{\text {\textit {c,t}}..t^{\prime }} = \sum _{i=t}^{t^{\prime }} x_{c,i}$, that is we replace time intervals t, t+1,…,t ^′ by one block of time intervals t..t ^′. During this block, converter c consumes $E_{c} x_{\text {\textit {c,t}}..t^{\prime }}$ electricity and produces $H_{c} x_{\text {\textit {c,t}}..t^{\prime }}$ heat. Using this notation, the state of charge equation for a block 8..10 of time intervals t=8,9 and 10 is

$$ s_{c,11} = s_{c,8} + H_{c} x_{c,8..10} - D_{c,8..10} $$

((7))

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.

In our example, we consider rough planning variables x _c,6..7, x _c,8..10, and x _11..15. In summary, all state of charge equations are

$$\begin{array}{@{}rcl@{}} s_{\text{\textit{c,t}}+1} &=& s_{\text{\textit{c,t}}} + H_{c} x_{\text{\textit{c,t}}} - D_{\text{\textit{c,t}}} \ \ \text{for} \ \ t = 1, \dots, 5 \\ s_{c,8} &=& s_{c,6} + H_{c} x_{c,6..7} - D_{c,6..7} \\ s_{c,11} &=& s_{c,8} + H_{c} x_{c,8..10} - D_{c,8..10} \\ s_{c,16} &=& s_{c,11} + H_{c} x_{c,11..15} - D_{c,11..15} \end{array} $$

((8))

for every $c \in \mathcal {C}$.

The capacity constraints of buffers remain the same, so

$$ L_{\text{\textit{c,t}}} \le s_{\text{\textit{c,t}}} \le U_{\text{\textit{c,t}}} \ \ \text{for} \ \ t \in \{1,2,3,4,5,6,8,11,16\} $$

((9))

The operational constraints of converters now are

$$\begin{array}{rcll} &x_{\text{\textit{c,t}}} &\in \{0,1\} &\text{for}~t \in \{1,2,3\} \\ 0 \le &x_{\text{\textit{c,t}}} &\le 1 &\text{for}~t \in \{4,5\} \\ 0 \le &x_{c,6..7} &\le 2 \\ 0 \le &x_{c,8..10} &\le 3 \\ 0 \le &x_{c,11..15} &\le 5 \end{array} $$

and the objective is

$$\begin{array}{crclr} \text{minimize} \ m \\ \text{where} & m &\ge& \sum_{c \in \mathcal{C}} E_{c} x_{\text{\textit{c,t}}} &\!\! \text{for}~t\! \in \{1,2,3,4,5\} \\ & 2m &\ge& \sum_{c \in \mathcal{C}} E_{c} x_{c,6..7} \\ & 3m &\ge& \sum_{c \in \mathcal{C}} E_{c} x_{c,8..10} \\ & 5m &\ge& \sum_{c \in \mathcal{C}} E_{c} x_{c,11..15} \end{array} $$

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

Both methods of direct control (global and time scale MILP control) are applied to a specific case involving 135 houses of the Meppel project. Each house is equipped with two heat pumps, one for domestic hot water and one for space heating. The applied energy system in Meppel is explained in detail in [18]. For the present paper, we investigate the quality of load balancing of heat pump electricity demand. We compare both types of direct control with results of reference PI control that is determined by each heat pump by itself, based on signals from house and storage thermostats. To generate heat demand profiles, the following approach is followed:

• 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

Suitable methods to determine space heating demand are listed in [19] and include modelling of thermal network, radiant time series, and transfer function methods. In the following, we develop a so-called grey-box model, i.e., a simplified model of a house based on the thermal network approach. An advantage of this approach is that dynamic heat demand of a house is described by only a few differential equations which are easy to integrate into smart grid control algorithms. However, model parameters may not correlate very well with true physical characteristics of the concerned building. For the purpose of the present paper, absolute accuracy of the space heating demand model is of less importance, hence we determined model parameters by physical estimation and not, e.g., by applying identification techniques using simulated or measured data. Validation of accuracy of the applied space heating demand model is part of our future work. The thermal network representation of the simplified model is shown in Fig. 2. For the indoor part, the model contains lumped thermal masses for the interior zone and floor structure. We assume an air-based heating system which gives heating input directly into the zone. Towards the exterior, there are two lumped thermal masses, one for the interior walls and one for the exterior walls of the house envelope. The main heat loss contributions to the ambient are through the envelope walls, windows, and roof.

Fig. 2

Applied thermal network model

The model equations are formulated as follows:

$${} {\fontsize{8}{12}{\begin{aligned} C_{\text{wo}} \cdot \frac{dT_{\text{wo}}}{\text{dt}} &= \frac{(T_{a}-T_{\text{wo}})\cdot A_{\text{wall}}}{R_{wo-a}} + \frac{(T_{\text{wi}}-T_{\text{wo}})\cdot A_{\text{wall}}}{R_{wi-wo}} \\ C_{\text{wi}} \cdot \frac{dT_{\text{wi}}}{\text{dt}} &= \frac{(T_{\text{wo}}-T_{\text{wi}})\cdot A_{\text{wall}}}{R_{\text{wi-wo}}} + \frac{(T_{z}-T_{\text{wi}}) \cdot A_{\text{wall}}}{R_{z-\text{wi}}} \\ C_{f} \cdot \frac{dT_{f}}{\text{dt}} &= \frac{(T_{z}-T_{f})\cdot A_{f}}{R_{z-f}} + \frac{(T_{cs}-T_{f})\cdot A_{f}}{R_{f- \text{cs}}} \\ C_{z} \cdot \frac{dT_{z}}{\text{dt}} &= D_{\text{\textit{c,t}}} + q_{\text{vent}} + q_{\text{inf}} + q_{\text{gain}} + \frac{(T_{f}-T_{z})\cdot A_{f}}{R_{z-f}} \\ & \quad+ \frac{(T_{\text{wi}}\,-\,T_{z})\cdot A_{\text{wall}}}{R_{\text{wall}}} + \frac{(T_{a}-T_{z})\cdot A_{\text{windows}}}{R_{\text{windows}}}\\ &\quad+ \frac{(T_{a}-T_{z})\cdot A_{\text{roof}}}{R_{\text{roof}}} \end{aligned}}} $$

((10))

The terms in these equations are explained in Table 1. Besides heat loss to the ambient, infiltration- and ventilation-associated heat losses (q _inf and q _vent) are determined by the following equations:

$$\begin{array}{@{}rcl@{}} q_{\text{inf}}=\phi_{\text{air,inf}}\cdot \rho_{\text{air}}\cdot c_{p, \text{air}}\cdot (T_{a}-T_{z}) \\ q_{\text{vent}}=\phi_{\text{air,vent}}\cdot \rho_{\text{air}}\cdot c_{p, \text{air}}\cdot (T_{a}-T_{z}) \cdot (1-\gamma_{\text{hr}}) \end{array} $$

((11))

Table 1 Nomenclature energy system characterization

The used terms are also explained in Table 1. Ventilation air flow values ϕ _air,vent in m ³/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.

The heat demand D _c,t is calculated by solving the model equations including a simple zone temperature control. For this, setpoint schedules for the zone temperature are defined in section “House and household case information”. The demand D _c,t for every house c and time interval t is determined by the following rules:

• 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 non-negative 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

We consider a total of 135 households and we define three types of households (see Table 2) living in semi-detached and detached houses. Houses will be built in three construction phases and later phases will have a better thermal insulation due to tightening regulations. In this study, we assume that in each phase, 30 semi-detached and 15 detached houses will be built. Rc values of semi-detached houses in the respective construction phases are 3.5, 5.0, and 7.5 m ² K/W and the Rc values of detached houses are 5.0, 7.5, and 10.0 m ² K/W.

Table 2 Types of households, the number of household types in both types of houses, schedules of higher temperature setting, and hot water demands

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 15-min time step which is required as a minimum time step for controlling the thermal storage and to calculate heat pump operation times. To obtain 15-min weather data, we applied linear interpolation of the 1-h 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

In Fig. 3, we show the total electric energy consumption of the three methods of control (reference control, global MILP control, and time scale MILP control). The reference control, which resembles PI 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.

Fig. 3

Case results

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.

Another important quality aspect of the heat pump operation schedules obtained from the optimized control is the experienced thermal comfort within the houses. This is influenced by the upper bound U _c of the state of charge s _c,t (Eq. (3)). In Fig. 4, we show the total electric energy consumption of the heat pumps sorted from highest to lowest consumption. The MILP control results in a more horizontal distribution than the reference control. In fact, depending on the value of the upper bound U _c , a whole continuum of slopes between the slope of the reference control and MILP control is possible. For the reference control, the value of U _c is approximately zero, i.e., the interior temperature equals the defined setpoints in time, except for the warming up and cooling down times when the setpoint is increased or decreased from a previous state. During those times, the interior temperature needs some time to change due to thermal mass effects within the house. Part of our future research is to investigate the quality of the interior temperature in relation to the value of the upper bound U _c .

Fig. 4

Sorted case results

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.