Optimised Heat Pump Management for Increasing Photovoltaic Penetration into the Electricity Grid

Abstract

Advanced control of heat pumps with thermal storage and photovoltaics has recently been promoted as a promising solution to help decarbonise the residential sector. Heat pumps and thermal storage offer a valuable flexibilisation mean to integrate stochastic renewable energy sources into the electricity grid. Heat pump energy conversion is nonlinear, leading to a challenging nonlinear optimisation problem. However, issues like global optimum uncertainty and the time-consuming methods of current nonlinear programming solvers draw researchers to linearise heat pump models that are then implemented in faster and globally convergent linear programming solvers. Nevertheless, these linearisations generate some inaccuracies, especially in the calculation of the heat pump’s coefficient of performance ($\mathit{C}\mathit{O}\mathit{P}$). In order to solve all of these issues, this paper presents a heuristic control algorithm (HCA) to provide a fast, accurate and near-optimal solution to the original nonlinear optimisation problem for a single-family house with a photovoltaic system, using real consumption data from a typical Swiss house. Results highlight that the HCA solves this optimisation problem up to 1000 times faster, yielding an operation that is up to 49% cheaper and self-consumption rates that are 5% greater than other nonlinear solvers. Comparing the performance of the HCA and the linear solver intlinprog, it is shown that the HCA provides more accurate heat pump control with an increase of up to 9% in system Operating Expense OPEX and a decrease of 8% in self-consumption values.

1. Introduction

The implementation of heat pumps (HP) with thermal energy storage (TES) appears to be a promising solution to meet the demand for hot water while promoting an efficient and low-carbon energy network [1]. Additionally, smart technologies offer the possibility of optimising heat pump scheduling coupled with TES, allowing for increased integration of decentralised renewable electricity generation into the grid. This would also reduce the need for expensive grid reinforcements and would reduce residential energy expenditures [2,3,4]. Thus, research regarding optimal control of heat pumps applied to residential energy systems is becoming a relevant topic in recent years.

There are currently many different approaches in the literature related to the subject of optimised energy management in residential energy systems. Among these, heuristic control strategies offer the simplest method to manage energy systems. Heuristics-based algorithms rely on simple rules based on experimental knowledge which allow solving detailed discrete or continuous problems using low computational resources. They cannot guarantee optimality, but they can generate solutions reasonably close to the global optimum. Riesen et al. [5] presented a simple rule-based control algorithm, valid only with a feed-in limit, that operated a residential photovoltaic (PV) and battery storage system. Employing a simulated PV clear-sky production profile as an input, this algorithm achieved nearly the same electricity cost savings as a linear MATLAB solver with exact forecasts, outperforming the latter when real forecasts were used. Thygesen and Karlsson [4] developed a rule-based control algorithm with a PV production forecast seven hours ahead, varying the set point of a domestic hot water (DHW) tank coupled with a ground source heat pump, attaining a gain of 7% in the annual self-consumption rate. Alimohammadisagvand et al. [6] developed three different rule-based control algorithms based on previous, current and future hourly electricity prices to minimise the life-cycle cost of the same kind of system using different TES sizes and set points. Bee et al. [7] highlighted the decreased in energy exchanges with the grid using a rule-based algorithm aiming at maximising PV self-consumption and considering both building and domestic hot water heating. Rule-based algorithms have also been used in [8,9] showing a limited optimisation capability compared to other kinds of optimisation algorithms such as Dynamic Programming. Furthermore, a different heuristic approach has been developed in [10] to perform load shifting strategies applied to domestic appliances of a household, concluding that these flexible loads had hardly any impact on the annual self-consumption. Salpakari and Lund [8] showed the same results after utilising a more complex Dynamic Programming algorithm. Furthermore, other kinds of heuristic solvers such as nature-inspired particle swarm optimisation and evolutionary algorithms were successfully tested either alone [11] or as a part of a two-step hybrid optimisation [12] to solve optimisation problems determining the optimal control and design of residential renewable energy systems.

More advanced control methods have also been developed using mathematical programming. The nature of the optimal control problem (OCP) for the system of study is defined by many authors as nonlinear, due to Coefficient of Performance ($\mathit{C}\mathit{O}\mathit{P}$) dependency on the water supply temperature and HP part-load efficiency [8,13,14,15]. Nonlinear OCPs are often addressed using different nonlinear programming (NLP) solvers minimising targets such as energy costs [16,17], HP power consumption [13] and user discomfort [16]. However, studies conducted in [13,16,17] presented results only over one-day periods, restricted by the time-consuming mathematical background of these solvers. This gives a motivation to simplify the original OCP into a linear model since linear programming (LP) solvers can guarantee global convergence and require low computational efforts. Investigations conducted in [12,18] performed LP and Mixed Integer Linear Programming (MILP) techniques using linear $\mathit{C}\mathit{O}\mathit{P}$ expressions with correction factors based on the theoretical Carnot cycle efficiency. Performance maps given by manufacturers were also used in [13,15,19,20]. In [20,21,22,23], a linear $\mathit{C}\mathit{O}\mathit{P}$ was also assumed by considering the water supply temperature constant to find the optimal operation of HPs coupled with underfloor heating systems. Likewise, Beck et al. [2] employed an MILP strategy with a constant $\mathit{C}\mathit{O}\mathit{P}$, assuming that heat source and sink temperatures are constant. These considerations can be supported by long HP operation cycles and smooth water supply temperature changes induced by underfloor heating or building large thermal inertia. However, the relatively low thermal inertia of residential DHW tanks and the discrete DHW consumption profile imply a high variability of the tank temperature through time. Wanjiru et al. [24] considered this premise but chose a constant $\mathit{C}\mathit{O}\mathit{P}$ value for an air-to-water heat pump (AWHP) water heater delivering hot water service to a DHW tank between 45 °C and 50 °C in a farmhouse in South Africa. It is expected that this hypothesis will barely affect the optimal scheduling decision for heat pumps, but it may not provide an accurate measurement of its thermal energy provision [20]. Moreover, linearisation of HP models that neglect its part-load dependency may also lead to underestimation of the optimal heat pump sizing [14].

The limitations previously observed in NLP and LP solvers have motivated the design and development of a heuristic algorithm using accurate HP modelling that could provide a near-optimal solution using low computational resources. The novelty of this study lies in the introduction of a unique heuristic control algorithm (HCA) which aims to minimise system operating expenses, including costs due to electricity exchange with the grid and HP. None of the reviewed heuristic and optimisation algorithms present neither a similar algorithm design nor such objective function. The present algorithm focuses on finding the right balance between optimality and simplicity. Using a purely heuristic approach, and by pre-selecting the possible operation points of the HP, the proposed algorithm is able to compute control trajectories close to the optimal solution in a much shorter time than standard NLP solvers.

In order to test the performance of our HCA, this paper conducts an assessment and comparison of this algorithm with other optimisation solvers. Assessment criteria rely mainly on the optimisation of heat pump management for enhancing PV self-consumption while minimising system operating expenditures. Additional indicators, such as the algorithm’s running time, reliability and the system’s energy efficiency, are considered in this evaluation. Thus, the HCA’s performance is compared with outputs given by MATLAB NLP solvers GlobalSearch and fmincon and the MATLAB LP solver intlinprog. Simulations are carried out over short and long periods (48 h and one month respectively) in the summer and winter using real consumption data from a typical Swiss house as well as calculated PV production. Details on the model, data used and simulation framework are given in the next section.

2. Methodology and Simulation Framework

This section describes the considered system, the modelling, capital costs and lifespan of its components and all the data input—e.g., weather conditions, household electricity and domestic hot water demand profiles along with electricity prices—required for the simulation. It also presents the optimisation of the HP operation using numerical solvers and the proposed HCA.

2.1. System Components and Modelling

The system under study is shown in Figure 1. It consists of a grid-connected PV system coupled with a variable-speed AWHP and a domestic hot water (DHW) storage unit for a single-family house. The grid-connected PV system generates electricity that can be either consumed by the household’s electricity demand or injected into the grid. PV electricity can also be used by the AWHP which, absorbing additional energy from the outside air, is able to provide enough thermal energy to the hot water tank to cover the household’s DHW demand when required. The electricity grid supplies the household’s demand in case of insufficient renewable energy generation and can absorb PV excess respecting the imposed feed-in regulation.

Modelling of the grid-connected PV system includes the model of a grid-connected inverter and uses the PVLIB toolbox [25] to simulate the DC PV production. Regarding the combination of the HP and the TES unit, the thermal capacities of the water supply and return circuit are neglected. Therefore, the water supply temperature at the outlet of the HP equals the inner tank temperature. All the elements of this system are assumed to comply with their model in all circumstances. System faults are not considered and coordination between all system agents is assumed to be ideal. If not, resilient control mechanisms could be in principle implemented [26,27].

2.1.1. Inverter

After comparing several inverter models in the literature [9,28,29], we selected the model for high-efficiency inverters proposed by [29]. In Notton’s paper, this model performed accurately when compared to a broad database of 21 inverters from the market. The other models didn’t show testing procedures as rigorous as this one (to the best of the author’s knowledge) but could also be valid for this study. The power at the output of the inverter is described by:

$$P_{\mathrm{PV},\mathrm{AC}}={\eta }_{\mathrm{inv}}·P_{\mathrm{PV},\mathrm{DC}}\mathrm{where}{\eta }_{\mathrm{inv}}=\frac{p}{p+p_0+k{p}^2},$$

(1)

where p is the ratio between the $P_{\mathrm{PV},\mathrm{AC}}$ and the nominal inverter power, $p_0$ and k are numerical constants described in Reference [29].

2.1.2. Heat Pump

The heat generated ${\dot{Q}}_{\mathrm{hp}}$ by a HP depends on the $\mathit{C}\mathit{O}\mathit{P}$ and the input electrical power $P_{\mathrm{hp}}$ as in Equation (2):

$${\dot{Q}}_{\mathrm{hp}}\left(t\right)=\mathit{C}\mathit{O}\mathit{P}\left(t\right)·P_{\mathrm{hp}}\left(t\right).$$

(2)

Nonlinear and linear models of a monovalent variable-speed AWHP are here considered. The corresponding $\mathit{C}\mathit{O}\mathit{P}$s will be then used in the OCP formulations defined in Section 2.2.3, Section 2.2.4 and Section 2.2.5. Various models of the $\mathit{C}\mathit{O}\mathit{P}$ value evolution with both ambient and water source temperature are used by the community. A quadratic formulation is proposed in [20]. In [13,30], the authors used a linear model while, in [2], the $\mathit{C}\mathit{O}\mathit{P}$ is assumed constant. There is a trade-off between modelling simplicity and model accuracy. However, in [13], the authors claim that increasing the order of the linear fit does not bring any useful accuracy with respect to the control problem. For this reason, the linear approximation of Equation (3) is used:

$$\mathrm{linear}:\mathit{C}\mathit{O}\mathit{P}\left(T_{\mathrm{amb}}\left(t\right),T_{\mathrm{tank}}\left(t\right)\right)=d_0+d_1·T_{\mathrm{amb}}\left(t\right)-d_2·T_{\mathrm{tank}}\left(t\right),$$

(3)

where $T_{\mathrm{amb}}\left(t\right)$ and $T_{\mathrm{tank}}\left(t\right)$ are the ambient temperature and the tank temperature at simulation time t, respectively, and $d_i$ are coefficients listed in

Table A1. Definition of MATLAB main classes, parameters and variables used in the simulation.

Pos	Property	Unit	Description	Value
DHW Tank (dtank)
Inverter (inv)
Environment (en)
1	$P_{\mathrm{PV},\mathrm{DC}}$	W	Photovoltaic DC-production	$M_{525601x1}$
2	$P_{\mathrm{load}}$	W	Household electric consumption profile	$M_{525601x1}$
Temps (ti)
3	$t_0$	min	Starting time of simulation	$t_0ϵ$ [1, 525,601]
4	$t_{\mathrm{end}}$	min	End of simulation time	$t_{end}ϵ$ [1, 525,601]
5	$t_{\mathrm{step}}$	min	Simulation time step	1
Heat Pump (hp)
6	$d_0$	Unitless	Independent $\mathit{C}\mathit{O}\mathit{P}$ coefficient	5.5930
7	$d_1$	1/°C	$T_{\mathrm{amb}}$-dependent $\mathit{C}\mathit{O}\mathit{P}$ coefficient	0.0569
8	$d_2$	1/°C	$T_{\mathrm{tank}}$-dependent $\mathit{C}\mathit{O}\mathit{P}$ coefficient	−0.0661
9	$a_1$	Unitless	1st order part-load efficiency coefficient	8.3350
10	$a_2$	Unitless	2nd order part-load efficiency coefficient	−38.0747
11	$a_3$	Unitless	3rd order part-load efficiency coefficient	104.6758
12	$a_4$	Unitless	4th order part-load efficiency coefficient	−159.6927
13	$a_5$	Unitless	5th order part-load efficiency coefficient	121.4477
14	$a_6$	Unitless	6th order part-load efficiency coefficient	−35.9697
Instantaneous Control Algorithm (ins)
15	$\mathit{C}\mathit{O}\mathit{P}$	Unitless	Heat pump coefficient of performance	$M_{525601x1}$
16	$P_{\exp }$	W	Exporting grid flux	$M_{525601x1}$
17	$P_{\mathrm{hp}}$	W	Heat pump electric power	$M_{525601x1}$
18	$P_{\mathrm{imp}}$	W	Importing grid flux	$M_{525601x1}$
19	$P_{\mathrm{load}}$	W	Household electricity consumption	$M_{525601x1}$
20	$P_{\mathrm{loss}}$	W	Power curtailment due to $\mathit{F}\mathit{i}\mathit{L}$	$M_{525601x1}$
21	$P_{\mathrm{PV},\mathrm{AC}}$	W	Photovoltaic-AC production	$M_{525601x1}$
22	$T_{\max }$	°C	Maximum heat pump supply temperature	65
23	$T_{\min }$	°C	Minimum tank temperature	55
Optimisation (opt)
24	$\mathit{C}\mathit{O}{\mathit{P}}^{opt}$	Unitless	Optimised heat pump coefficient of performance	$M_{96x1}$
25	$\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$	Unitless	$\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ = $\mathit{C}\mathit{O}{\mathit{P}}^{opt}·f_c^{opt}$	$M_{96x1}$
26	$costvec$	CHF	Matrix of cost values for all tested heating states	$M_{qxj}$
27	$f_c^{opt}$	Unitless	Optimised HP part-load correction factor	$M_{96x1}$
28	$G_f^{opt}$	W	Optimised net grid power flux $P_{\mathrm{PV},\mathrm{AC}}^{opt}-P_{\mathrm{load}}^{opt}-P_{\mathrm{hp}}^{opt}$	$M_{96x1}$
29	${\widehat{P}}_{hp}^{opt}$	Unitless	Optimised heat pump power profile normalized to $P_{nom}$	$M_{96x1}$
30	${\widehat{P}}_{{\mathrm{hp}}_{\mathrm{old}}}^{opt}$	Unitless	${\widehat{P}}_{\mathrm{hp}}^{opt}$ vector taken from the last decision as reference	$M_{96x1}$
31	$P_{hp,pv}^{opt}$	W	Amount of PV electricity consumed in the optimised HP operation	$M_{96x1}$
32	$P_{\mathrm{load}}^{opt}$	W	Load electric consumption forecast	$M_{96x1}$
33	$P_{\mathrm{PV},\mathrm{AC}}^{opt}$	W	Photovoltaic-AC production Forecast	$M_{96x1}$
34	$P_{\mathrm{PV},\mathrm{AC}}^{opt}$	W	$t_{\mathrm{step}}^{opt}$-min PV production forecast	$M_{96x1}$
35	${\dot{Q}}_{\mathrm{hp}}^{opt}$	$W_{th}$	Optimised heat pump thermal power	$M_{96x1}$
36	$T_{\max }^{opt}$	°C	Maximum temperature value of $T_{\mathrm{tank}}^{opt}$ vector	Var

in the Appendix A. An improvement to the $\mathit{C}\mathit{O}\mathit{P}$ formulation is to include $\mathit{C}\mathit{O}\mathit{P}$ part-load dependency, which can be added to Equation (3) as a correction factor $f_c$. It is based on a 6th order polynomial fitting curve of the normalized part-load efficiency data for alternating current (AC) inverter compressors reported by [31]:

$$\mathrm{nonlinear}:\mathit{C}\mathit{O}{\mathit{P}}_m\left(T_{\mathrm{amb}}\left(t\right),T_{\mathrm{tank}}\left(t\right),P_{\mathrm{hp}}\left(t\right)\right)=\mathit{C}\mathit{O}\mathit{P}\left(t\right)·f_c\left(P_{\mathrm{hp}}\left(t\right)\right),$$

(4)

where ${\displaystyle f_c\left(P_{\mathrm{hp}}\left(t\right)\right)=\sum _{n=1}^{n=6}{a}_n·{\left(\mathit{P}\mathit{L}\mathit{R}\left(t\right)\right)}^n}$ in which the variable $\mathit{P}\mathit{L}\mathit{R}$ designates the HP part-load ratio, defined as the quotient of the HP’s electric power $P_{\mathrm{hp}}$ and its nominal electric power $P_{\mathrm{nom}}$. Values of coefficients $a_n$ can be also found in Table A1 in the Appendix A as well as an illustration of the fitting curve in Figure A1.

2.1.3. Domestic Hot Water Tank

A physics-based model of a DHW tank is created considering the volume $\left(V\right)$, diameter $\left(D\right)$ and standard heat losses $\left({\dot{Q}}_{l,st}\right)$ data from [32]. The tank’s height $\left(H\right)$ and surface area $\left(A\right)$ are calculated based on geometric equations $H=\frac{4·V}{1000·\pi ·D^2}$ and $A=\pi ·D·H+2·\pi ·{\left(\frac{D}{2}\right)}^2$, respectively. Furthermore, ${\dot{Q}}_{l,st}$ is defined using a linear fit ${\dot{Q}}_{l,st}=0.075·V+109.2$ upon heat loss data for volumes ranging between 300 and 1000 litres. H and D are expressed in meters, V in litres and ${\dot{Q}}_{l,st}$ in watts. All these variables are used in Equation (5) to calculate the tank’s conductive heat transfer coefficient $\left(U\right)$, in $\frac{\mathrm{W}}{{\mathrm{m}}^2·\mathrm{K}}$, assuming a design temperature difference $(\Delta T_{\mathrm{design}})$ of 45 °C, in accordance with the set of procedures accepted in [33] to determine standing heat losses for factory-insulated storage water heaters:

$$U=\frac{{\dot{Q}}_{l,st}}{A·\Delta T_{\mathrm{design}}}.$$

(5)

The tank’s energy dynamics are modelled applying the energy conservation principle, considering the tank as a homogeneous control volume. A basic formulation is presented in Equation (6):

$${\dot{E}}_{\mathrm{tank}}\left(t\right)={\dot{Q}}_{\mathrm{hp}}\left(t\right)-{\dot{Q}}_l\left(t\right)-{\dot{Q}}_d\left(t\right),$$

(6)

where ${\dot{E}}_{\mathrm{tank}}$ is the energy variation of the TES unit, ${\dot{Q}}_{\mathrm{hp}}$ designates the HP’s thermal power transferred into the tank and ${\dot{Q}}_l$ and ${\dot{Q}}_d$ are heat losses due to conductive heat transfer with the environment and domestic hot water demand, respectively. All are expressed in watts. Due to the water incompressibility assumption, the stored energy variation in the tank can be simplified into equation ${\dot{E}}_{\mathrm{tank}}\left(t\right)={\rho }_w·c_w·V·{\dot{T}}_{\mathrm{tank}}\left(t\right)$; $c_w$ and ${\rho }_w$ the specific heat capacity and density of water, with values of 4180 $\mathrm{J}/\mathrm{kg}·\mathrm{K}$ and 1 kg/L, respectively. Thermal energy transferred from the heat pump into the tank was calculated using the equation ${\dot{Q}}_{\mathrm{hp}}\left(t\right)=P_{\mathrm{hp}}\left(t\right)·\mathit{C}\mathit{O}{\mathit{P}}_m\left(t\right)$. The heat released through the reservoir walls into the surroundings is expressed as ${\dot{Q}}_l\left(t\right)=U·A·\left(T_{\mathrm{tank}}\left(t\right)-T_{\mathrm{room}}\right)$, assuming a constant tank room temperature ($T_{\mathrm{room}}$) of 20 °C. Energy transfer associated with the supply of domestic hot water is formulated as ${\dot{Q}}_d\left(t\right)={\dot{m}}_d\left(t\right)·c_w·\left(T_{ws}-T_{\mathrm{cold}}\right)$. In our case, the hot water consumption temperature ($T_{ws}$) has a predefined value of 55 °C while the cold water ($T_{\mathrm{cold}}$) enters the tank at 15 °C. The term ${\dot{m}}_d\left(t\right)$ represents the domestic hot water mass flow requested by consumers in kg/s.

2.2. HP Operation Optimisation

Optimal operation of the HP relies on two distinct procedures: an algorithm defining the (cost-based) optimal operation scheduling for the next 48 h and an instantaneous control algorithm (ICA) that effectively controls the HP and keeps its safe operation as close as possible to the scheduled plan. This approach is referred as hierarchical control in [34,35]. Optimal operation scheduling is based on the load, consumption and weather forecasts and is determined by optimisation solvers used in this work or by the novel HCA developed in this study. Note that this optimisation model is based on the input of perfect forecast data. Figure 2 shows all the elements involved in the simulation framework object of study.

2.2.1. Data Input

The electricity consumption profile consists of real electricity consumption measurements in 15-min time steps from a Swiss single-family house. The annual electricity consumption accounts for a total of 3.85 MWh. Hot water mass flow usage data is modelled from real temperature measurements in one-minute time steps registered in a DHW tank of another Swiss family house over most of the year 2016. A water consumption profile was obtained by multiplying the temperature drops higher than the heat loss (equal to 0.4 °C/h) by 20 L/°C.

The PV production is simulated using the PV-LIB toolbox [25] for a PV system comprising SunPower SPR-315 modules installed facing south with a tilt of 27°. Input data comprises real environmental data collected from the Mühleberg (CH) weather station in 10-min time steps over the year 2015. The retail electricity price is chosen at a fixed value of 20 cts/kWh, based on price values surveyed by [36] in 2015. Monetary values are expressed in Swiss francs (CHF). The feed-in tariff (FiT) price is set at a constant revenue of 6 cts/kWh following prices experienced in countries with large renewable energy penetration such as Germany [37]. The same retail and feed-in tariff electricity prices are considered as forecast in the optimisation. A feed-in limit (FiL) on the active power of 70% of the nominal PV array capacity is chosen based on the German regulation for PV electricity storage, §9-2.b of [38]. The capital expenditure (CAPEX) for the HP is determined in CHF by using the linear function $5680+1.24·P_{nom}$ presented in [12], $P_{nom}$ being the nominal electric power of the HP in W. Lifespan values of the HP account for 100,000 h and 50,000 cycles, based on data provided by HP manufacturers [39].

2.2.2. Parameters of Optimal Control Problem

For the arrangement shown in Figure 1, the main decision variable is $P_{\mathrm{hp}}^{opt}$. However, this decision variable implicitly controls other variables such as the exporting and importing power fluxes with the grid ($P_{\exp }^{opt}$ and $P_{\mathrm{imp}}^{opt}$ respectively, expressed in W), the power curtailed due to the feed-in limit ($P_{\mathrm{loss}}^{opt}$, in W), and the temperature reached in the tank ($T_{\mathrm{tank}}^{opt}$, in  °C). Therefore, for each discrete state ($t^{opt}$), the model tries to optimise these five variables. The control time step ($t_{\mathrm{step}}^{opt}$) accounts for 30 min, throughout a 48 h time horizon. Consequently, the optimisation problem is established with 96 different discrete states and 480 decision variables.

2.2.3. Nonlinear Optimisation Problem

The objective function of the NLP optimisation is to minimise the operative expenses due to the electricity exchanges with the grid, which is formulated in Equation (7):

$$f_{\mathit{N}\mathit{L}\mathit{P}}=\sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(P_{\mathrm{imp}}^{opt}\left(t^{opt}\right)·C_b-P_{\exp }^{opt}\left(t^{opt}\right)·C_s\right)·t_{\mathrm{step}}^{opt},$$

(7)

where $C_b$ and $C_s$ are the buying and selling price of electricity, in $\frac{\mathrm{CHF}}{\mathrm{W}·\min }$. This objective function is subjected to Equation (6) using the nonlinear $\mathit{C}\mathit{O}{\mathit{P}}_m$ expression given by Equation (4). A set of nonlinear constraints is completed by adding Equation (8), which forbids the system to export and import electricity simultaneously:

$$P_{\mathrm{imp}}^{opt}\left(t^{opt}\right)·P_{\exp }^{opt}\left(t^{opt}\right)=0.$$

(8)

Equation (9) formulates the power flow balance in the switchboard node shown in Figure 1, which gives a linear constraint for each time step $t^{opt}$. The left part of the equation is the incoming active power which equals the right part that expresses the outgoing active power:

$$P_{\mathrm{imp}}^{opt}\left(t^{opt}\right)+P_{\mathrm{PV}},{\mathrm{AC}}^{opt}\left(t^{opt}\right)=P_{\mathrm{load}}^{opt}\left(t^{opt}\right)+P_{\mathrm{hp}}^{opt}\left(t^{opt}\right)+P_{\exp }^{opt}\left(t^{opt}\right).+P_{\mathrm{loss}}^{opt}\left(t^{opt}\right)$$

(9)

Furthermore, the physical limitations of the system were captured by the bounds $l{b}_{\mathit{N}\mathit{L}\mathit{P}}\le \left[P_{\mathrm{hp}}^{opt},T_{\mathrm{tank}}^{opt},P_{\mathrm{loss}}^{opt},P_{\mathrm{imp}}^{opt},P_{\exp }^{opt}\right]\le u{b}_{\mathit{N}\mathit{L}\mathit{P}}$, with $u{b}_{\mathit{N}\mathit{L}\mathit{P}}$ and $l{b}_{\mathit{N}\mathit{L}\mathit{P}}$ being the upper and lower limits with values of $\left[P_{\mathrm{nom}},65,\infty ,\infty ,\mathit{F}\mathit{i}\mathit{L}\right]$ and $\left[0,55,0,0,0\right],$ respectively.

2.2.4. Mixed-Integer Linear Optimisation Problem

The $\mathit{M}\mathit{I}\mathit{L}\mathit{P}$ objective function comprises three parts. The first one (Equation (10a)) aims to minimise the grid exchange cost (difference between the purchased electricity from the grid and the excess PV energy sold to the grid). The second one (Equation (10b)) minimises the operating cost of the HP, composed of the running cost (cost associated with the time degradation of the HP components) and the switching cost (cost associated with each switch on of the HP). The last one (Equation (10c)) allows the tank temperature to drop below 55 °C with an additional cost:

$$\begin{array}{ccc}\hfill f_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}=& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(P_{\mathrm{imp}}^{opt}\left(t^{opt}\right)·C_b-P_{\exp }^{opt}\left(t^{opt}\right)·C_s\right)·t_{\mathrm{step}}^{opt}\hfill & \hfill \mathrm{grid}\mathrm{exchange},\end{array}$$

(10a)

$$\begin{array}{ccc}\hfill +& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(R^{opt}\left(t^{opt}\right)·C_{\mathrm{run}}\right)·t_{\mathrm{step}}^{opt}+S^{opt}\left(t^{opt}\right)·C_{sw}\hfill & \hfill \mathrm{HP}\mathrm{operation},\end{array}$$

(10b)

$$\begin{array}{ccc}\hfill +& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(e^{opt}\left(t^{opt}\right)+\mathit{O}{\mathit{B}}^{opt}\left(t^{opt}\right)\right)·C_{ob}\hfill & \hfill \mathrm{temperature}\mathrm{boundary},\end{array}$$

(10c)

where $e^{opt}$ is a slack continuous variable used for constraint relaxation and $S^{opt}$ and $R^{opt}$ are binary variables that take values equal to 1 during HP switching and running periods respectively and 0 when the HP is off. $\mathit{O}{\mathit{B}}^{opt}$ also refers to a binary variable that is equal to 1 when the tank temperature drops below 55 °C (out of desirable boundaries) and 0 when it stays above this set point. The high cost of 10 CHF was chosen for $C_{ob}$ to allow the system to experience temperatures below 55 °C only if there is no possibility to keep the tank temperature within the suitable range 55–65 °C. Constants $C_{sw}$ and $C_{\mathrm{run}}$ set a switching and running cost on the HP, based on the HP’s CAPEX and lifespan values from Section 2.2.1. $C_{sw}$ was computed as the ratio between the HP’s CAPEX and HP’s nominal number of cycles while $C_{\mathrm{run}}$ has as divider the number of operating minutes.

This objective function is also constrained by the same power flow balance as the nonlinear case (Equation (9)). Equation (6) is also chosen to represent HP-DHW storage energy dynamics. In this case, the variable $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ is defined as the linear $\mathit{C}\mathit{O}\mathit{P}$ formulation given by Equation (3) with a constant tank temperature of 60 °C. Additional linear inequalities (Equations (11)–(15)) are designed to provide a more detailed optimisation. Equations (11)–(13) express the relation between the HP running and switching variables. Equation (14) forces the HP to operate in power ranges higher than or equal to 20% of its nominal power. Finally, Equation (15) activates variable $\mathit{O}{\mathit{B}}^{opt}\left(t^{opt}\right)$ every time the temperature values go below the tank temperature’s set point:

$$\begin{array}{cc}\hfill -S^{opt}\left(t^{opt}\right)-R^{opt}(t^{opt}-1)+R^{opt}\left(t^{opt}\right)& \le 0,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill 2·S^{opt}\left(t^{opt}\right)+R^{opt}(t^{opt}-1)-R^{opt}\left(t^{opt}\right)& \le 1,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill P_{\mathrm{hp}}^{opt}\left(t^{opt}\right)-M·R^{opt}\left(t^{opt}\right)& \le 0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill -P_{\mathrm{hp}}^{opt}\left(t^{opt}\right)+0.2·P_{\mathrm{nom}}·R^{opt}\left(t^{opt}\right)-e^{opt}\left(t^{opt}\right)& \le 0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill -T_{\mathrm{tank}}^{opt}\left(t^{opt}\right)-M·\mathit{O}{\mathit{B}}^{opt}\left(t^{opt}\right)& \le -55,\hfill \end{array}$$

(15)

where M is a constant with a high value (M = ${10}^4$) used to modify the restriction degree of constraints. The upper and lower system boundaries ($u{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$ and $l{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$) are defined in

Table 1. Lower ($l{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$) and upper ($u{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$) system boundaries at each time step for the $\mathit{M}\mathit{I}\mathit{L}\mathit{P}$ optimisation problem.

	${\mathit{P}}_{\mathbf{hp}}^{\mathit{opt}}$	${\mathit{T}}_{\mathbf{tank}}^{\mathit{opt}}$	${\mathit{P}}_{\mathbf{loss}}^{\mathit{opt}}$	${\mathit{P}}_{\mathbf{imp}}^{\mathit{opt}}$	${\mathit{P}}_{\mathbf{exp}}^{\mathit{opt}}$	${\mathit{S}}^{\mathit{opt}}$	${\mathit{R}}^{\mathit{opt}}$	${\mathit{e}}^{\mathit{opt}}$	$\mathit{O}{\mathit{B}}^{\mathit{opt}}$
$l{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$	0	10	0	0	0	0	0	0	0
$u{b}_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$	$P_{\mathrm{nom}}$	65	∞	∞	$\mathit{F}\mathit{i}\mathit{L}$	1	1	∞	1

.

2.2.5. Heuristic Optimal Control Problem

This version of the OCP is able to take advantage of the different benefits presented in the NLP and MILP methods. On the one hand, it can adopt the same set of constraints used by NLP in Section 2.2.3. On the other hand, its objective function $f_{\mathit{H}\mathit{C}\mathit{A}}$ (presented in Equation (16)) admits binary variables, thereby minimising electricity and HP operating costs such as $f_{\mathit{M}\mathit{I}\mathit{L}\mathit{P}}$ introduced in Equation (10). The HCA objective function also comprises three parts. The two first (Equation (16a,b)) are similar to those in Equation (10). The last one (Equation (16c)) is a heuristic factor described in Equation (17):

$$\begin{array}{ccc}\hfill f_{\mathit{H}\mathit{C}\mathit{A}}=& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(P_{\mathrm{imp}}^{opt}\left(t^{opt}\right)·C_b-P_{\exp }^{opt}\left(t^{opt}\right)·C_s\right)·t_{\mathrm{step}}^{opt}\hfill & \hfill \mathrm{grid}\mathrm{exchange},\end{array}$$

(16a)

$$\begin{array}{ccc}\hfill +& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}\left(R^{opt}\left(t^{opt}\right)·C_{\mathrm{run}}\right)·t_{\mathrm{step}}^{opt}+S^{opt}\left(t^{opt}\right)·C_{sw}\hfill & \hfill \mathrm{HP}\mathrm{operation},\end{array}$$

(16b)

$$\begin{array}{ccc}\hfill +& \sum _{t^{opt}=1}^{t_{\mathrm{end}}^{opt}=96}h\left(t^{opt}\right)\hfill & \hfill \mathrm{heuristic}\mathrm{factor}.\end{array}$$

(16c)

The heuristic factor h assists the solver to make decisions considering the HP’s part-load efficiency. The constant $\alpha$ of Equation (17) has a value of $3.33\times {10}^{-4}{\mathrm{W}}^{-1}{\min }^{-1}$:

$$h\left(t^{opt}\right)=\left\{\begin{array}{cc}\frac{1}{\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}\left(t^{opt}\right)·\alpha }·C_b\hfill & ,\mathrm{if}\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}\left(t^{opt}\right)>0,\hfill \\ 0\hfill & ,\mathrm{if}\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}\left(t^{opt}\right)\le 0.\hfill \end{array}\right.$$

(17)

2.2.6. Optimisation Solvers

To solve the nonlinear OCP, we used the potentiality of MATLAB gradient-based nonlinear solvers fmincon and GlobalSearch. However, their large computational burden and their dependency on the starting point ($x_0$) raise the interest to linearize the original OCP into a simpler version shown in Section 2.2.4, for which the fast MILP MATLAB solver intlinprog was selected. Moreover, the heuristic OCP was fed into the novel HCA developed in this work. A flowchart of this HCA algorithm can be found in Figure 3. Other variables included in this flowchart are listed in Table A1 in the Appendix A.

Following Figure 3, the HCA starts initializing all variables of the optimisation and takes the current tank temperature ($T_{\tan \mathrm{k}}\left(t\right)$) from the ICA, described in Section 2.2.7, as an input. It then determines the iteration (i) when the tank temperature would drop below the lower limit ($T_{\min }$ = 55 °C). When this happens, the HCA computes all the possible states from the beginning of the optimisation $j=1$ until $j=i$, increasing the HP’s power from 0 to 100% in 20% power steps (discrete power levels are designated by the variable $q=0\cdots 5$). With this, the cost of each alternative is evaluated according to the objective function expressed in Equation (16). At the end, a certain power value ($q_{\min }·P_{\mathrm{step}}$) at a particular discrete state ($j_{\min }$) is chosen, which yields the lowest value of the objective function. Then, the tank temperature continues to decrease according to the DHW consumption profile and tank conductive heat losses. The next time the tank temperature falls below its threshold, the same iterative process is repeated until the tank temperature profile meets the temperature constraint on the whole time horizon.

One peculiarity addressed by the solver is that, at times of large DHW consumption, the system could possibly not provide sufficient TES and HP capacity to be able to maintain the temperature above the set point by running the HP in all possible previous states. At this point, the maximum operation variable $\mathit{MO}\left(i\right)$ is set to 1 and the solver has no other choice than to proceed to the next iteration to provide newly available states to use the HP in order to recover the tank temperature as soon as possible. To sum up, the overall benefits and limitations of each solver are presented in

Table 2. Benefits, limitations and optimal control problems (OCP) fed into the different optimisation solvers. $x_0$ refers to the starting point provided to solvers fmincon and GlobalSearch.

Solver	OCP	Benefits	Limitations
fmincon	NLP	Detailed model	Only continuous Equations
		Global Optimisation	Largest computational time
			Optimality not guaranteed
			High $x_0$ dependency
GlobalSearch	NLP	Detailed model	Only continuous Equations
		Global Optimisation	Largest computational time
		Low $x_0$ dependency	Optimality not guaranteed
intlinprog	MILP	Simplest model	Accuracy
		Lowest computational time
		Continuous/Discrete Equations
		Optimality guaranteed
		Global Optimisation
HCA	NLP	Detailed model	Optimality not guaranteed
	+ Heuristics	Low computational time	Sequential optimisation
		Continuous/Discrete Equations

.

2.2.7. Instantaneous Control Algorithm (ICA)

ICA refers to the rule-based control routine that governs the present system state based on the control signal sent by the optimisation solver every 12 h (see Figure 2). The fact of having different optimisation and control time steps (30 min and 1 min respectively) may create a mismatch between the DHW demand profiles used in the optimisation and experienced by the ICA. This mismatch can lead to situations in which the algorithm, running with the optimised HP schedule, has to face unpredicted DHW demand. Therefore, to be able to address these operational disturbances and preserve the HP lifespan, three protection measures were implemented in the ICA.

The first measure prevents the HP from starting unless if the tank temperature is at least 2 °C lower than its upper-temperature limit. The second measure keeps the HP in the off-state for at least 10 min after the tank reaches its maximum temperature (65 °C). The third measure overcomes unexpected real-time tank temperature drops (below 55 °C) by calculating whether there would be sufficient HP power scheduled for the next 15 min to keep temperatures within the feasible range and, if not, it uses the available HP power to raise the tank temperature up to 56 °C (value chosen to avoid short-cycling) in a period of at least 15 min.

2.3. Evaluation Metrics

The OPEX and self-consumption variables are used as the main indicators to assess the performance of the different optimisation solvers (case i = opt), as presented in Section 3.1, and also to evaluate the operation of the ICA under the control of the HCA and intlinprog solvers in particular (case i = no superscript), as shown in Section 3.2. Equation (18) shows the OPEX function, which includes the cost of buying and selling electricity as well as the HP running and switching costs:

$$\begin{array}{c}\hfill {\mathrm{OPEX}}^i=\sum _{t_0^i}^{t_{\mathrm{end}}^i}[\left(P_{\mathrm{imp}}^i\left(t^i\right)·C_b-P_{\exp }^i\left(t^i\right)·C_s+R^i\left(t^i\right)·C_{\mathrm{run}}\right)·t_{\mathrm{step}}^i+S^i\left(t^i\right)·C_{sw}].\end{array}$$

(18)

To quantify the PV electricity integration capability of these solvers, a self-consumption variable is used with its definition established in Equation (19):

$$S{C}^i(\%)=\frac{{\sum }_{t_0^i}^{t_{\mathrm{end}}^i}[\min (P_{\mathrm{load}}^i\left(t^i\right)+P_{\mathrm{hp}}^i\left(t^i\right),P_{\mathrm{PV}}^i,\mathrm{AC}\left(t^i\right)]}{{\sum }_{t_0^i}^{t_{\mathrm{end}}^i}[P_{\mathrm{PV}}^i,\mathrm{AC}\left(t^i\right)]}·100.$$

(19)

To measure the gain of self-consumption accomplished by the adoption of optimal HP management strategies, a reference case is determined, consisting of a situation in which the same HP energy profile calculated by the solvers would have to be imported from the grid. For this case, the self-consumption rate is calculated neglecting the $P_{\mathrm{hp}}^i$ term placed in the numerator of Equation (19).

3. Results and Discussion

This section assesses the HCA’s performance with regard to OPEX minimisation and self-consumption enhancement. Indicators such as system energy efficiency, solver accuracy and running time are also considered. The MATLAB solvers intlinprog, fmincon and GlobalSearch are used as benchmarks. Firstly, the HCA’s effectiveness is evaluated in 48 h horizon short-term simulations in winter and summer for a system size of 1 kW HP, 600 l TES and 3 kWp PV. This size is chosen to reproduce a typical household configuration with a large DHW tank to enhance PV self-consumption. The second section restricts the analysis of long-term system operation to cases controlled by the HCA and intlinprog using the same system sizing applied to short-term simulations in winter and summer. fmincon and GlobalSearch are excluded from this case since they could not always find a reasonable local optimum or did not converge at all. Those issues arise as a consequence of the relatively flat search path experienced by NLP solvers when computing the values of the objective function defined in the nonlinear OCP (Section 2.2.3) throughout the feasible domain determined by the applicable set of constraints. The CPU used for simulations comprises an Intel Core i5-3230 (2.60 GHz) processor with 8 GB RAM on a 64-bits operating system.

3.1. Short-Term Simulations

The evaluation is carried out for a 48 h optimisation period including PV production profiles corresponding with one clear-sky day and one cloudy day. The initial point $x_0$ of the optimisation using solvers fmincon and GlobalSearch was given by the result of intlinprog. This procedure avoids using a non-optimal initial point $x_0$ and reduces these solvers’ running times. For the sake of clarity, this section is focused solely on comparing the performance of the different optimisation solvers. The performance of the ICA based on optimal HP energy management signals will be analysed in Section 3.2.

Figure 4 displays the main system variables computed by the different optimisation solvers for summers’ simulations. Regarding temperature values, it can be seen that, in all cases, the tank temperature profiles end at the set point of 55 °C.

Analysing the HP power scheduling, we note that fmincon decides to operate the HP at mid-power levels intermittently over most of the PV production period to enhance the part-load efficiency values, keeping supply temperatures low. This results in the highest $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values of all the solvers while the PV self-consumption is increased. Nevertheless, this power profile implies continuous HP switching, which would reduce its lifespan rapidly. For GlobalSearch, the solver yields a HP behaviour similar to that of fmincon, except that it runs the HP during some periods at very low power values with $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values lower than 1. The reason for both deficiencies relies on the inability of derivative-based NLP solvers to utilize discrete variables. On the one hand, this means that fmincon and GlobalSearch can accept $P_{hp}^{opt}\left(t\right)$ only as a continuous decision variable. Hence, although there is a penalty in efficiency, these routines switch the HP on at very low power values as low power ranges are allowed and this inefficient operation incurs no cost in their objective function (Equation (7)). On the other hand, this limitation makes it difficult to assign a cost to the HP switching and running time, leading to detrimental intermittent HP power profiles.

For the HCA, the solver is able to concentrate all the HP power working at a higher power level compared to the NLP solvers, shaving PV injection peaks more effectively. Despite the fact that working at higher power levels results in lower HP $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values, the HCA control strategies provide the best match between HP consumption and PV generation, accomplishing the highest share of PV self-consumption (44%) during summer simulations.

Table 3. OPEX ($OPE{X}^{opt}$), PV self-consumption ($S{C}^{opt}$) and mean run time (${\overline{t}}_{run}$) results of the heuristic control algorithm (HCA), intlinprog, fmincon and GlobalSearch solvers for a single 48-h optimisation for winter and summer. The self-consumption rate compared to a reference case ($S{C}_{ref}^{opt}$) is shown, considering that the entire HP load is met only by importing grid electricity. Upper/lower uncertainty limits ($\Delta t_{run,max}^{+}$ and $\Delta t_{run,max}^{-}$, respectively) denote, with respect to its mean value, the maximum run time variability observed after repeating each optimisation twenty times, except for GlobalSearch, which was run only five times.

	Winter				Summer
	HCA	int.	fmin.	GS	HCA	int.	fmin.	GS
$S{C}_{ref}^{opt}$ (%)	66				21
$S{C}^{opt}$ (%)	92	96	99	99	44	39	39	39
$OPE{X}^{opt}$	7.10	6.28	7.87	7.54	1.45	0.97	2.84	2.45
${\overline{t}}_{run}$ (s)	80	<1	1047	76,891	57	2	920	62,448
$\Delta t_{run,max}^{+}$ (s)	4	<1	18	13,676	2	<1	9	8215
$\Delta t_{run,max}^{-}$ (s)	3	<1	8	15,425	3	<1	7	11,360

shows that self-consumption rates for the other solvers stay 5% below the HCA’s result.

The HP scheduling given by intlinprog also follows a similar profile as the HCA’s. However, the definition of a linear $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ with a constant supply temperature has as consequences that the solver considers neither the reduction of $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values when increasing the tank temperature nor the decrease of $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ when the HP is working at full load due to part-load ratio. All these factors lead to an overestimated $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ profile, which contributes to reducing the HP running periods with respect to those generated by the HCA.

Looking at the ${\mathrm{OPEX}}^{opt}$ values in Table 3, it can be seen that the HCA offers the most cost-competitive solution during the summer (1.45 CHF) among MATLAB nonlinear solvers (fmincon and GlobalSearch lead to 96% and 69% more expensive optimal control strategies, respectively). These results are only improved by the 33% cheaper intlinprog operation. In winter, NLP solvers present the same behaviour as for the summer period, achieving the highest $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values due to a combination of low supply temperatures and HP operation at optimal part-load efficiency ranges. For this case, they achieve local consumption of almost 100% of the PV generation. Moreover, the HCA and intlinprog attain the lowest OPEX values since they can run the HP more smoothly than the NLP solvers while promoting the self-consumption of locally generated PV electricity. Note that, despite the fact that intlinprog provides the cheapest operation strategy for summer and winter short-term simulations, those results are based on overestimated $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ values. Section 3.2 will present some of the impacts that this $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ inaccuracy produces in the real-time control of the system.

Studying solver calculation time, values in Table 3 suggest that the implementation of current NLP algorithms cannot provide suitable real-time optimal control due to their large computational burden. This statement depends very much on the hardware and on the efficiency of the implementation of such nonlinear solver. Thus, better implementation of such solvers may draw a different conclusion. However, considering the MATLAB implementation, only the HCA and intlinprog can effectively deliver control trajectories in a computation time smaller than the simulation time-step, i.e., enabling real-time control. Note that all results gathered in Table 3 are totally deterministic, except running time values, which depend partly on external factors related to computer processing hierarchy, and GlobalSearch output, which may vary according to the set of trial starting points chosen randomly by a scatter search algorithm integrated in this solver [40]. As a consequence of the first exception, short-term winter and summer simulations were repeated twenty times for the HCA, intlinprog and fmincon while they were repeated five times for GlobalSearch in order to obtain a mean run time value and its uncertainty range shown in Table 3. Due to the second exception, the GlobalSearch ${\mathrm{OPEX}}^{opt}$ and $S{C}^{opt}$ values produced in the fastest run out of the five executed for winter and summer periods are shown. The GlobalSearch main variable profiles depicted in Figure 4 correspond to this case.

3.2. Long-Term Simulations

Long-term simulations are performed to optimise the energy balance of the system during a winter (1 to 31 January) and a summer month (1 to 30 June) with the aim of minimising the system’s OPEX values and enhancing self-consumption by using the optimal HP scheduling calculated every 12 h by the HCA and intlinprog. In the case of DHW, low and high monthly consumption profiles are provided to simulations in summer and winter conditions, respectively. These long-term simulations also allow us to check the robustness of the algorithm with respect to the time dependent environmental conditions and demand profiles.

Running simulations for the use case defined in this study, it is observed that intlinprog records execution times significantly lower than the HCA. As in Section 3.1, summer and winter long-term simulations are repeated twenty times for both solvers, obtaining for each of them, run time data samples of 1200 and 1240 values that correspond to simulations run during the summer and winter months respectively. Distribution of these data samples is represented through several boxplots shown in Figure 5. Note that each of these run time values is the time the HCA or intlinprog spends to calculate an optimal HP management strategy for a 48-h period, which occurs twice per day during the whole simulation month (the ICA running time is not included in these values). The larger PV generation experienced in summer made intlinprog evaluate an increased number of possibilities to allocate HP power scheduling for self-consumption enhancement. This leads to longer running time values compared to results given for the winter month. Even so, intlinprog can converge in 75% of all cases below three seconds and one second in summer and winter seasons, respectively.

For the HCA, the third quartile of the running time during the winter month (January) reached a value of 133 s, whereas 54 s was recorded for the summer month (June). This time difference depends on the HCA’s number of decisions. In winter, factors such as low ambient temperatures and high DHW consumption decrease the HP’s $\mathit{C}\mathit{O}\mathit{P}$ and increase the heat demand, thereby raising the number of times the algorithm needs to operate the HP in order to maintain the tank temperature above the set point (55 °C). Additionally, it is observed that allocating the HP load within the PV generation period also helped the HP to increase its efficiency, reducing the HCA’s number of decisions. This can be explained by the correlation of solar energy availability with higher ambient temperatures. Furthermore, 

Table 4. The system OPEX and self-consumption rate (SC) achieved when using the HCA and intlinprog solvers in winter and summer simulations. ($SC\mathrm{ref}$) denotes the self-consumption rate considering a reference case in which the entire HP load is met only by importing grid electricity.

	Winter		Summer
	HCA	Intlinprog	HCA	Intlinprog
SCref (%)	64		28
SC (%)	88	96	46	47
OPEX (CHF)	166.66	153.18	37.78	36.27

highlights that the HCA is able to provide an optimal control strategy that can adjust the HP power to self-consume as much PV electricity as intlinprog. The self-consumption rate experienced in the winter simulation controlled by the HCA is just 8% lower than is the case when intlinprog is used. During the summer, although the PV resource increases significantly, this difference reduces to less than 1%.

Thus, these small differences in self-consumption contribute to obtain a tight comparison regarding OPEX values: the HCA could control real-time system performance during the winter at a cost that is 9% higher than intlinprog, whereas for summer this difference is reduced to 4%.

Regarding solver accuracy, Figure 6 shows that considering an LP approach to solve this optimisation problem presents several drawbacks. When tank temperatures are higher/lower than 60 °C, intlinprog’s $\mathit{C}\mathit{O}\mathit{P}$ is overestimated/underestimated (yellow area). Even when temperatures remain close to 60 °C, instantaneous $\mathit{C}\mathit{O}{\mathit{P}}_m$ values are frequently lower than those of intlinprog due to part-load efficiency reasons (blue area). Consequently, according to 3th protection measure (see Section 2.2.7), the ICA has to react using auxiliary heat pump capacity to provide additional heating when the tank temperature drops unexpectedly below the set point (here on 25 January around 10:30 h (purple area)). Moreover, Figure 7 quantifies the impact of the previous shortcomings over system real-time operation by showing values of the supplementary heat pump energy used by the ICA (additional to the power predicted in the optimisation) as well as the DHW deficit durations (periods when tank temperature is below the set point) for winter and summer months.

It is highlighted that cases managed by the HCA undergo certain periods of DHW deficit (80 h and 0.5 h for winter and summer months, respectively). However, this deficit is a consequence of two external factors: the heavy DHW peaks experienced, especially during the winter, which lead to temperature drops larger than the 10 K desirable temperature band, and the time step difference between the simulation and the optimisation, which causes tank temperatures to sometimes fall below the set point as the ICA experiences unexpected DHW demand peaks during a few minutes that could not be anticipated by the optimisation, which uses forecast input data with a 30-min time resolution. Nevertheless, these unexpected real-time temperature drops make real $\mathit{C}\mathit{O}\mathit{P}$ values higher than $\mathit{C}\mathit{O}{\mathit{P}}^{opt}$ values at that time, enabling the HP to deliver a higher thermal energy output than predicted during the next minutes, which helps to restore or sometimes even surpass the optimised temperature value calculated by the HCA at the end of that time step. Thus, although the DHW deficit periods are inevitable, the HCA scheduling remains valid even for the aforementioned time step difference as it avoids the use of any supplementary HP power ($t_{aux}$, $E_{hp,aux}$ and $Q_{hp,aux}$ are 0).

For the case of intlinprog, in addition to the external factors mentioned in the case of the HCA, $\mathit{C}\mathit{O}{\mathit{P}}_m$ overestimation leads to a deficit of HP thermal energy provision, which contributes to an increase of DHW deficit duration with respect to the HCA case (108 h and 7 h for winter and summer months, respectively). This forces the ICA to use the supplementary HP power to compensate for this deficit, providing instantaneously around 2 kWh of thermal energy for more than one hour for the winter month, while 4 kWh and two hours are registered for the summer month. This higher activity of the ICA in summer compared to winter can be explained from the fact that the HP operates for less time in summer than in winter, which makes it less probable that the ICA can find any scheduled HP power in the next 15 min when trying to restore unexpected temperature drops, being forced to use more supplementary HP power.

It is worthwhile to note that the optimised HP scheduling used by the ICA is refreshed every 12 h. If a longer refreshing period were chosen, intlinprog inaccuracies would have a greater misleading effect on real-time system operation.

4. Conclusions

In this paper, a novel heuristic control algorithm (HCA) is developed to optimise HP management with very little computational efforts by enhancing PV self-consumption while minimising the operating costs of a typical Swiss single-family house equipped with a 1-kW variable-speed AWHP, a 600-L DHW storage unit and a 3-kWp PV system. Inputs to this optimisation problem are exact weather forecasts, household consumption patterns, electricity prices and a feed-in limit. Furthermore, the HCA output is used by an instantaneous control algorithm (ICA) to govern real-time system operation. Within this framework, this research presents a comparison of system performance under HCA control to same-use cases controlled by the nonlinear and linear optimisation solvers fmincon, GlobalSearch and intlinprog over different scenarios, considering measures such as the system OPEX, PV self-consumption and energy efficiency as well as solver running time and reliability.

In short-term winter and summer simulations, the HCA converges up to 1000 times faster than fmincon and GlobalSearch. Furthermore, the HCA allows for the definition of a nonlinear and discrete AWHP model, enabling the algorithm to run the HP more smoothly than the fragmented and low-efficiency HP operation produced by NLP solvers. All these benefits make the HCA capable of providing optimised HP management up to 49% cheaper than the nonlinear solvers. Additionally, the HCA achieves the highest PV self-consumption share of all solvers during the summer (44%), and reached a 92% rate over the winter period (7% below the best case).

However, long-term winter and summer simulations confirm the reliability and robustness of the HCA to converge in 100% of the cases. In contrast, NLP solvers are not able to find a reasonable local optimum due to their limitation to solve optimisation problems with a nearly flat search space. Thus, long-term simulations using HCA and intlinprog show that the HCA’s strategy self-consumes as much PV electricity as intlinprog, registering an increase of up to 9% in system operating costs with respect to this linear solver. Moreover, HCA control with a more precise HP representation leads to fewer ICA interventions. In contrast, the intlinprog control requires more ICA interventions to ensure that the tank temperature stays within the desired range, as the linearisation of the HP model frequently induces $\mathit{C}\mathit{O}{\mathit{P}}_m^{opt}$ miscalculations, overestimating in some cases the HP’s thermal energy provision. A further advantage of the HCA is that it can be easily implemented on lean hardware.

Future work would focus mainly on the development of the HCA code architecture towards a more efficient and faster version. One solution would be to reduce the temporal search space of the HCA for faster optimisation. Furthermore, the HCA’s performance will be tested under new conditions such as implementing a broader database of DHW consumption profiles, a new building model to optimise HP operation for DHW and heating purposes, and even a stratified DHW tank model and real forecasts to evaluate its effects on system operation.