**Need help with data science or mathematical modeling?**I do consulting work in Norway. Read about my previous work experience and reach out to me for more information.

# Heating a hot tub

- 23. January 2022 (modified 26. January 2022)
- #optimization

Norwegian electricity prices have soared this winter. My friend Marius wondered if one could use estimated future electricity prices to heat up a hot tub in an optimal way. In this article we’ll formulate this problem as an optimization problem and propose two approximation algorithms.

## Problem statement

Assume that we want to heat up an outdoor hot tub to a comfortable 40 degrees Celsius 48 hours into the future.

Let us start by reviewing the thermodynamics of heating and cooling water:

**Heating.**Heating water is a simple process: if you add twice as much energy you get twice as much increase in temperature. In other words, the specific heat capacity is practically constant in the range we’re considering. The same applies to time—if we supply the same amount of energy per time unit (watts), then doubling the time will double the temperature increase.**Cooling.**Cooling water follows Newton’s law of cooling: the rate of heat loss of a body is proportional to the difference in temperature between the body and its environment.

Let \(T(t)\) denote the temperature at a given time \(t\). If there is no heating, then the temperature is governed by the differential equation

where \(\lambda\) is a constant (which must be estimated in practice). This equation has a closed solution

Let \(x(t)\) be a function with codomain \(\{0, 1\}\), indicating whether we supply energy to the system. If we do, then the increase in temperature per time unit is denoted by \(h\). Then the system is governed by

and we can fire up a Runge-Kutta integrator to simulate the system.

The goal is to find an optimal function \(x^{\star}(t)\) so that \(T(t=48) \geq 40\) in the cheapest way possible. There’s no loss in saying “at least 40 degrees” instead of “exactly 40 degrees”, since when we minimize the total price we won’t overheat the hot tub because it would cost more.

Let \(p(t)\) be the price at time \(t\), and recall that \(x(t)\) denotes whether we’re heating the hot tub or not. We introduce a slack variable \(\alpha\), which allows us to find a feasible solution even if it’s impossible to heat the hot tub to 40 degrees. In that case we want to heat it as much as possible. The optimization problem becomes

where \(\{0,1\}^{\mathbb{R}}\) denotes the set of all functions \(x: \mathbb{R} \to \{0,1\}\). We’re looking for an optimal function \(x(t)\) which switches the heating on and off to produce the cheapest possible way to heat the hot tub.

Here’s one way to heat up a hot tub to 40 degrees Celsius, given future electricity prices \(p(t)\). The green areas show when the heating is turned on. We overheat twice when the price is low, and end up with 40 degrees. This particular strategy is not optimal. Later we’ll see that we can get the cost down to \(268\), which is much lower than the value of \(685\) shown in the figure below.

## Solving a discretized version of the problem

One way to obtain a practical solution to the problem is to discretize it and use Mixed Integer Programming (MIP). Let’s assume that \(T_{\text{env}} = 0\), so without heat added to the system it undergoes a simple exponential decay.

Let’s discretize the dynamics of the problem as

where \(0 < k < 1\) is a decay constant measuring the heat loss and \(H\) is the heat (measured in degrees Celsius) added to the system in one discrete time step. For instance, when \(k=0.99\) the heat loss is slow, but when \(k=0.9\) the heat loss is rather large.

**The final temperature is linear in the inputs.**
Let \(f\) be a function mapping from a set of discrete time periods to a final temperature in the last discrete time step.
For instance \(f(\{x_3, x_8\})\) is the final temperature achieved if we heat in discrete time periods \(x_3\) and \(x_8\).
It’s easy to see that the final temperature is linear:

If heat \(H\) is added at time \(t\), and \(t=1, 2, \ldots, n\), then the contribution to the final temperature by the heat added at time \(t\) is \(H_t = H k^{n-t}\). For instance, if \(x_{n-2} = 1\), so heat is added in the third-to-last time period, then the temperature remaining in period \(t=n\) would be \(H_{n-2} = H k^{2}\).

**A Mixed Integer Program.**
Defining \(H_t := H k^{n-t}\), one possible MIP formulation is:

The parameter \(\alpha\) is a slack variable, introduced to ensure infeasibility even if it’s impossible to heat the hot tub to 40 degrees. As always, \(\epsilon\) is a small positive number, ensuring that we prioritize heating the hot tub before we prioritize lowering the total price \(\sum_t p_t x_t\).

**Examining solutions.**
Let’s generate a random price \(p(t)\), discretize it into hours, and examine solutions as we vary the heat loss parameter \(k\).

When \(k = 0.9\), the heat loss is so large that it makes sense to heat up just before using the hot tub. If we spend money on cheap electricity earlier on, the heat would just disappear anyway.

Conversely, if \(k = 0.99\), then there is almost no heat loss. A good strategy is essentially to supply heat to the hot tub when it’s cheap to do so.

A middle ground is found when \(k=0.95\). We supply heat in both the cheap periods, but more in the last one.

## Solving a continuous version of the problem

The MIP above solves a knapsack-like problem. In knapsack problems, the greedy strategy (ordering by heat per price, \(\text{ratio}_t := H_t / p_t\), and picking the largest value until \(\sum_t H_t x_t \geq 40\)) is not necessarily optimal. However, if we discretize further and further so each element becomes arbitrarily small, the difference between the optimal solution and the greedy solution becomes smaller and smaller. In the limit, the greedy strategy is optimal under the discretization that we used.

Let’s examine the continuous problem once again. Recall that we set \(T_{\text{env}} = 0\), so the system is governed by

If we supply heat at time \(t\), the temperature we supply will have decreased by a factor of \(\exp \left( \lambda (t - 48) \right)\) at the end of the time period. The continuous version of the ratio \(H_t / p_t\) is

The term containing \(T(t)\) is troublesome, because we \(T(t)\) depends on \(x(t)\), which we haven’t found yet. However, if we are willing to assume that \(T'(t) \gg 0\) when we heat the system, we can approximate the term \(\left[ h - \lambda T(t) \right]\) as \(h\) and get rid of the term \(\lambda T(t)\). Let’s make that assumption and see where it takes us. We define

Now the problem is to choose a function \(x(t)\) that separates out the intervals of \(t\) with the highest \(\text{ratio}(t)\) from the rest.

We achieve this using root-finding, given an unknown \(0 < \beta < \infty\), we define:

If \(\beta = 0\), we do not heat at all. As \(\beta \to \infty\), we heat all the time. It should be clear that the amount of heat we supply is a non-decreasing function of \(\beta\), and this means we can use root-finding.

We need to find a value of \(\beta\) so that we end up with 40 degrees in the end. To do this, we use a root finding algorithm (e.g. bisection). For any candidate \(\beta\) we run a simulation and return the final temperature of the system.

**Answer for the initial example.**

With my initial hand-made strategy (first figure in this article), the original example had a cost of \(685\). The strategy has a cost of \(268\), and is shown below:

**More examples.**

Here are three more examples with a different price function \(p(t)\) and a different value for \(h\). Each figure shows a good strategy, with different values for \(\lambda\).

## Notes

- Instead of approximating \(\left[ h - \lambda T(t) \right]\) as \(h\), we could replace \(T(t)\) with a constant \(T_{\text{const}}\). For instance, we could choose \(T_{\text{const}} = 20\), representing a value between the initial and final temperature.
- Instead of approximating \(\left[ h - \lambda T(t) \right]\) as \(h\), we could use an iterative approach. First we set \(T_1(x) = 0\), then we use root finding to find the optimal \(\beta\) and a resulting \(x_1(t)\). Then we use that \(x_1(t)\) to compute \(T_2(t)\) by simulation. This function \(T_2(t)\) is then used to find \(x_2(t)\), and so forth.