A cooling coffee cup model

Here we show an example (found in examples/coffee_cup) where we examine the changes in temperature of a cooling coffee cup that follows Newton’s law of cooling:

\[\frac{dT(t)}{dt} = -\kappa(T(t) - T_{env})\]

This equation tells how the temperature \(T\) of the coffee cup changes with time \(t\), when it is in an environment with temperature \(T_{env}\). \(\kappa\) is a proportionality constant that is characteristic of the system and regulates how fast the coffee cup radiates heat to the environment. For simplicity we set the initial temperature to a fixed value, \(T_0 = 95 ^\circ\text{C}\), and let \(\kappa\) and \(T_{env}\) be uncertain parameters. We give the uncertain parameters in the following distributions:

\[ \begin{align}\begin{aligned}\kappa &= \mathrm{Uniform}(0.025, 0.075),\\T_{env} &= \mathrm{Uniform}(15, 25).\end{aligned}\end{align} \]

Using a function

There are two approaches to creating the model, using a function or a class. The function method is easiest so we start with that approach. The complete for this example can be found in examples/coffee_cup/uq_coffee_function.py. We start by importing the packages we use:

import uncertainpy as un
import chaospy as cp                       # To create distributions
import numpy as np                         # For the time array
from scipy.integrate import odeint         # To integrate our equation

To create the model we define a Python function coffee_cup that takes the uncertain parameters kappa and T_env as input arguments. Inside this function we solve our equation by integrating it using scipy.integrate.odeint, before we return the results. The implementation of the model function is:

# Create the coffee cup model function
def coffee_cup(kappa, T_env):
    # Initial temperature and time array
    time = np.linspace(0, 200, 150)            # Minutes
    T_0 = 95                                   # Celsius

    # The equation describing the model
    def f(T, time, kappa, T_env):
        return -kappa*(T - T_env)

    # Solving the equation by integration
    temperature = odeint(f, T_0, time, args=(kappa, T_env))[:, 0]

    # Return time and model output
    return time, temperature

We could use this function directly in UncertaintyQuantification, but we would like to have labels on the axes when plotting. So we create a Model with the above run function and labels:

# Create a model from the coffee_cup function and add labels
model = un.Model(run=coffee_cup, labels=["Time (min)", "Temperature (C)"])

The next step is to define the uncertain parameters. We use Chaospy to create the distributions, and create a parameter dictionary:

# Create the distributions
kappa_dist = cp.Uniform(0.025, 0.075)
T_env_dist = cp.Uniform(15, 25)

# Define the parameter dictionary
parameters = {"kappa": kappa_dist, "T_env": T_env_dist}

We can now calculate the uncertainty and sensitivity using polynomial chaos expansions with point collocation, which is the default option of quantify. We set the seed to easier be able to reproduce the result.

# Set up the uncertainty quantification
UQ = un.UncertaintyQuantification(model=model, parameters=parameters)

# Perform the uncertainty quantification using
# polynomial chaos with point collocation (by default)
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(seed=10)

The complete code becomes:

import uncertainpy as un
import chaospy as cp                       # To create distributions
import numpy as np                         # For the time array
from scipy.integrate import odeint         # To integrate our equation


# Create the coffee cup model function
def coffee_cup(kappa, T_env):
    # Initial temperature and time array
    time = np.linspace(0, 200, 150)            # Minutes
    T_0 = 95                                   # Celsius

    # The equation describing the model
    def f(T, time, kappa, T_env):
        return -kappa*(T - T_env)

    # Solving the equation by integration
    temperature = odeint(f, T_0, time, args=(kappa, T_env))[:, 0]

    # Return time and model output
    return time, temperature


# Create a model from the coffee_cup function and add labels
model = un.Model(run=coffee_cup, labels=["Time (min)", "Temperature (C)"])

# Create the distributions
kappa_dist = cp.Uniform(0.025, 0.075)
T_env_dist = cp.Uniform(15, 25)

# Define the parameter dictionary
parameters = {"kappa": kappa_dist, "T_env": T_env_dist}

# Set up the uncertainty quantification
UQ = un.UncertaintyQuantification(model=model, parameters=parameters)

# Perform the uncertainty quantification using
# polynomial chaos with point collocation (by default)
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(seed=10)

Using a class

The model can also be created as a class instead of using a function. Most of the code is unchanged. The complete for this example is in examples/coffee_cup/uq_coffee_class.py. We create a class that inherits from Model. To add the labels we call on the constructor of the parent class and give it the labels.

# Create the coffee cup model
class CoffeeCup(un.Model):
    # Add labels to the model by calling the constructor of the parent un.Model
    def __init__(self):
        super(CoffeeCup, self).__init__(labels=["Time (s)", "Temperature (C)"])

We can then implement the run method:

    # Define the run method
    def run(self, kappa, T_env):
        # Initial temperature and time array
        time = np.linspace(0, 200, 150)            # Minutes
        T_0 = 95                                   # Celsius

        # The equation describing the model
        def f(T, time, kappa, T_env):
            return -kappa*(T - T_env)

        # Solving the equation by integration
        temperature = odeint(f, T_0, time, args=(kappa, T_env))[:, 0]

        # Return time and model output
        return time, temperature

Now, instead of creating a model from a model function, we initialize our CoffeeCup model:

# Initialize the model
model = CoffeeCup()

While the rest is unchanged:

# Create the distributions
kappa_dist = cp.Uniform(0.025, 0.075)
T_env_dist = cp.Uniform(15, 25)

# Define the parameters dictionary
parameters = {"kappa": kappa_dist, "T_env": T_env_dist}

# Set up the uncertainty quantification
UQ = un.UncertaintyQuantification(model=model, parameters=parameters)

# Perform the uncertainty quantification using
# polynomial chaos with point collocation (by default)
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(seed=10)