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:
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:
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)