# Modeling and Simulation in Python

``````

In :

# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

``````

## Inferring thermal resistance and thermal mass of a wall

This case study is based on Gori, Marincioni, Biddulph, Elwell, "Inferring the thermal resistance and effective thermal mass distribution of a wall from in situ measurements to characterise heat transfer at both the interior and exterior surfaces", Energy and Buildings, Volume 135, 15 January 2017, Pages 398-409, which I downloaded here.

The authors put their paper under a Creative Commons license, and make their data available here. I thank them for their commitment to open, reproducible science, which made this case study possible.

The goal of their paper is to model the thermal behavior of a wall as a step toward understanding the "performance gap between the expected energy use of buildings and their measured energy use". The wall they study is identified as the exterior wall of an office building in central London, not unlike this one.

The following figure shows the scenario and their model:

On the interior and exterior surfaces of the wall, they measure temperature and heat flux over a period of three days. They model the wall using two thermal masses connected to the surfaces, and to each other, by thermal resistors.

The primary methodology of the paper is a Bayesian method for inferring the parameters of the system (two thermal masses and three thermal resistances).

The primary result is a comparison of two models: the one shown here with two thermal masses, and a simpler model with only one thermal mass. They find that the two-mass model is able to reproduce the measured fluxes substantially better.

Tempting as it is, I will not replicate their method for estimating the parameters. Rather, I will

1. Implement their model and run it with their estimated parameters, to replicate the results, and

2. Use SciPy's `leastsq` to see if I can find parameters that yield lower errors (root mean square).

`leastsq` is a wrapper for some venerable FORTRAN code that runs "a modification of the Levenberg-Marquardt algorithm", which is one of my favorites (really).

Implementing their model in the ModSimPy framework turns out to be straightforward. The simulations run fast enough even when we carry units through the computation. And the results are visually similar to the ones in the original paper.

I find that `leastsq` is not able to find parameters that yield substantially better results, which suggest that the estimates in the paper are at least locally optimal.

First I'll load the units we need from Pint.

``````

In :

m = UNITS.meter
K = UNITS.kelvin
W = UNITS.watt
J = UNITS.joule
degC = UNITS.celsius
s = UNITS.second

``````
``````

Out:

second

``````

``````

In :

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

Q_in
Q_out
T_int
T_ext

2014-10-05 16:30:00
10.994
6.840
16.92
14.68

2014-10-05 16:35:00
10.952
6.012
16.92
14.69

2014-10-05 16:40:00
10.882
7.040
16.93
14.66

2014-10-05 16:45:00
10.798
8.880
16.93
14.59

2014-10-05 16:50:00
10.756
10.491
16.94
14.50

``````

The index contains Pandas `Timestamp` objects, which is good for dealing with real-world dates and times, but not as good for running the simulations, so I'm going to convert to seconds.

``````

In :

timestamp_0 = get_first_label(data)

``````
``````

Out:

Timestamp('2014-10-05 16:30:00')

``````

Subtracting the first `Timestamp` yields `Timedelta` objects:

``````

In :

time_deltas = data.index - timestamp_0

``````
``````

Out:

TimedeltaIndex(['0 days 00:00:00', '0 days 00:05:00', '0 days 00:10:00',
'0 days 00:15:00', '0 days 00:20:00', '0 days 00:25:00',
'0 days 00:30:00', '0 days 00:35:00', '0 days 00:40:00',
'0 days 00:45:00',
...
'2 days 23:10:00', '2 days 23:15:00', '2 days 23:20:00',
'2 days 23:25:00', '2 days 23:30:00', '2 days 23:35:00',
'2 days 23:40:00', '2 days 23:45:00', '2 days 23:50:00',
'2 days 23:55:00'],
dtype='timedelta64[ns]', length=864, freq=None)

``````

Then we can convert to seconds and replace the index.

``````

In :

data.index = time_deltas.days * 86400 + time_deltas.seconds

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

Q_in
Q_out
T_int
T_ext

0
10.994
6.840
16.92
14.68

300
10.952
6.012
16.92
14.69

600
10.882
7.040
16.93
14.66

900
10.798
8.880
16.93
14.59

1200
10.756
10.491
16.94
14.50

``````

The timesteps are all 5 minutes:

``````

In :

np.all(np.diff(data.index) == 300)

``````
``````

Out:

True

``````

Mark the columns of the `Dataframe` with units.

``````

In :

data.Q_in.units = W / m**2
data.Q_out.units = W / m**2
data.T_int.units = degC
data.T_ext.units = degC

``````

Plot the measured fluxes.

``````

In :

plot(data.Q_in, color='C2')
plot(data.Q_out, color='C0')
decorate(xlabel='Time (s)',
ylabel='Heat flux (W/\$m^2\$)')

``````
``````

``````

Plot the measured temperatures.

``````

In :

plot(data.T_int, color='C2')
plot(data.T_ext, color='C0')
decorate(xlabel='Time (s)',
ylabel='Temperature (degC)')

``````
``````

``````

### Params and System objects

Here's a `Params` object with the estimated parameters from the paper.

``````

In :

params = Params(
R1 = 0.076 * m**2 * K / W,
R2 = 0.272 * m**2 * K / W,
R3 = 0.078 * m**2 * K / W,
C1 = 212900 * J / m**2 / K,
C2 = 113100 * J / m**2 / K)

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

R1
0.076 kelvin * meter ** 2 / watt

R2
0.272 kelvin * meter ** 2 / watt

R3
0.078 kelvin * meter ** 2 / watt

C1
212900.0 joule / kelvin / meter ** 2

C2
113100.0 joule / kelvin / meter ** 2

``````

I'll pass the `Params` object `make_system`, which computes `init`, packs the parameters into `Series` objects, and computes the interpolation functions.

``````

In :

def make_system(params, data):
"""Makes a System object for the given conditions.

params: Params object

returns: System object
"""
R1, R2, R3, C1, C2 = params

init = State(T_C1 = Quantity(16.11, degC),
T_C2 = Quantity(15.27, degC))

ts = data.index
t_end = ts[-1] * s

return System(init=init,
R=Series([R1, R2, R3]),
C=Series([C1, C2]),
T_int_func=interpolate(data.T_int),
T_ext_func=interpolate(data.T_ext),
unit_temp=degC,
t_end=t_end, ts=ts)

``````

Make a `System` object

``````

In :

system = make_system(params, data)

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

init
T_C1    16.11 degC
T_C2    15.27 degC
dtype: o...

R
0    0.076 kelvin * meter ** 2 / watt
1    0.2...

C
0    212900.0 joule / kelvin / meter ** 2
1   ...

T_int_func
<function interpolate.<locals>.wrapper at 0x7f...

T_ext_func
<function interpolate.<locals>.wrapper at 0x7f...

unit_temp
degC

t_end
258900 second

ts
Int64Index([     0,    300,    600,    900,   ...

``````

Test the interpolation function:

``````

In :

system.T_ext_func(0), system.T_ext_func(150), system.T_ext_func(300)

``````
``````

Out:

(14.68, 14.684999999999999, 14.69)

``````

### Implementing the model

Next we need a slope function that takes instantaneous values of the two internal temperatures and computes their time rates of change.

The slope function gets called two ways.

• When we call it directly, `state` is a `State` object and the values it contains have units.

• When `run_ode_solver` calls it, `state` is an array and the values it contains don't have units.

In the second case, we have to apply the units before attempting the computation. `require_units` applies units if necessary:

The following function computes the fluxes between the four zones.

``````

In :

def compute_flux(state, t, system):
"""Compute the fluxes between the walls surfaces and the internal masses.

state: State with T_C1 and T_C2
t: time in seconds
system: System with interpolated measurements and the R Series

returns: Series of fluxes
"""
# unpack the temperatures
T_C1, T_C2 = state

# compute a series of temperatures from inside out
T_int = system.T_int_func(t)
T_ext = system.T_ext_func(t)

T = [require_units(T_int, system.unit_temp),
require_units(T_C1, system.unit_temp),
require_units(T_C2, system.unit_temp),
require_units(T_ext, system.unit_temp)]

# compute differences of adjacent temperatures
T_diff = np.diff(T)

# compute fluxes between adjacent compartments
Q = T_diff / system.R
return Q

``````

We can test it like this.

``````

In :

compute_flux(system.init, 0, system)

``````
``````

Out:

0    -10.657894736842135 delta_degC * watt / kelvin...
1    -3.0882352941176463 delta_degC * watt / kelvin...
2    -7.564102564102562 delta_degC * watt / kelvin ...
dtype: object

``````

Here's a slope function that computes derivatives of `T_C1` and `T_C2`

``````

In :

def slope_func(state, t, system):
"""Compute derivatives of the state.

state: position, velocity
t: time
system: System object

returns: derivatives of y and v
"""
Q = compute_flux(state, t, system)

# compute the net flux in each node
Q_diff = np.diff(Q)

# compute the rate of change of temperature
dQdt = Q_diff / system.C
return dQdt

``````

Test the slope function with the initial conditions.

``````

In :

slopes = slope_func(system.init, system.ts, system)

``````
``````

Out:

0     3.555499973097458e-05 delta_degC * watt / joule
1    -3.844086774341106e-05 delta_degC * watt / joule
dtype: object

``````
``````

In :

for y, slope in zip(system.init, slopes):
print(y, slope*s)

``````
``````

16.11 degC 3.555499973097458e-05 delta_degC * second * watt / joule
15.27 degC -3.844086774341106e-05 delta_degC * second * watt / joule

``````

Now let's run the simulation, generating estimates for the time steps in the data.

``````

In :

results, details = run_ode_solver(system, slope_func)
details

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

success
True

message
The solver successfully reached the end of the...

``````

Here's what the results look like.

``````

In :

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

T_C1
T_C2

0
16.11 degC
15.27 degC

2589
16.19153931824483 degC
15.119152775835897 degC

5178
16.2528855986394 degC
14.860249399364765 degC

7767
16.294553283278788 degC
14.58496288165028 degC

10356
16.317098239731408 degC
14.348725465924332 degC

``````
``````

In :

def plot_results(results, data):
plot(data.T_int, color='C2')
plot(results.T_C1, color='C3')
plot(results.T_C2, color='C1')
plot(data.T_ext, color='C0')
decorate(xlabel='Time (s)',
ylabel='Temperature (degC)')

plot_results(results, data)

``````
``````

``````

These results are similar to what's in the paper:

.

To get the estimated fluxes, we have to go through the results and basically do the flux calculation again.

``````

In :

def recompute_fluxes(results, system):
"""Compute fluxes between wall surfaces and internal masses.

results: Timeframe with T_C1 and T_C2
system: System object

returns: Timeframe with Q_in and Q_out
"""
Q_frame = TimeFrame(index=results.index, columns=['Q_in', 'Q_out'])

for t, row in results.iterrows():
Q = compute_flux(row, t, system)

Q_frame.row[t] = (-Q.magnitude,
-Q.magnitude)

return Q_frame

Q_frame = recompute_fluxes(results, system)

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

Q_in
Q_out

0
10.6579
7.5641

2589
9.93106
12.6725

5178
9.30414
17.6391

7767
8.90193
17.1162

10356
8.79081
16.4555

``````

Let's see how the estimates compare to the data.

``````

In :

def plot_Q_in(frame, data):
plot(frame.Q_in, color='gray')
plot(data.Q_in, color='C2')
decorate(xlabel='Time (s)',
ylabel='Heat flux (W/\$m^2\$)')

plot_Q_in(Q_frame, data)

``````
``````

``````
``````

In :

def plot_Q_out(frame, data):
plot(frame.Q_out, color='gray')
plot(data.Q_out, color='C0')
decorate(xlabel='Time (s)',
ylabel='Heat flux (W/\$m^2\$)')

plot_Q_out(Q_frame, data)

``````
``````

``````

These results are also similar to what's in the paper (the bottom row):

``````

In :

def compute_error(frame, data):
model_Q_in = interpolate(Q_frame.Q_in)(data.index)
error_Q_in = model_Q_in - data.Q_in

model_Q_out = interpolate(Q_frame.Q_out)(data.index)
error_Q_out = model_Q_out - data.Q_out

return np.hstack([error_Q_in, error_Q_out])

``````
``````

In :

errors = compute_error(Q_frame, data)

``````
``````

Out:

array([-0.33610526, -0.37832695, -0.39254863, ...,  2.26392988,
-0.7752054 ,  1.99865932])

``````

Here's the root mean squared error.

``````

In :

print(np.sqrt(np.mean(errors**2)))

``````
``````

2.17995317290043

``````

### Estimating parameters

Now let's see if we can do any better than the parameters in the paper.

In order to work with `leastsq`, we need a version of `params` with no units.

``````

In :

params = [0.076, 0.272, 0.078, 212900, 113100]

``````
``````

Out:

[0.076, 0.272, 0.078, 212900, 113100]

``````

Here's an error function that takes a hypothetical set of parameters, runs the simulation, and returns an array of errors.

``````

In :

def error_func(params, data):
"""Run a simulation and return an array of errors.

params: Params object or array
data: DataFrame

returns: array of float
"""
print(params)
system = make_system(params, data)
system = remove_units(system)

results, details = run_ode_solver(system, slope_func)
Q_frame = recompute_fluxes(results, system)
errors = compute_error(Q_frame, data)
print('RMSE', np.sqrt(np.mean(errors**2)))
return errors

``````

Testing `error_func`.

``````

In :

errors = error_func(params, data)

``````
``````

[0.076, 0.272, 0.078, 212900, 113100]
RMSE 2.17995317290043

Out:

array([-0.33610526, -0.37832695, -0.39254863, ...,  2.26392988,
-0.7752054 ,  1.99865932])

``````

Pass `error_func` to `leastsq` to see if it can do any better.

``````

In :

best_params, details = leastsq(error_func, params, data)

``````
``````

[7.600e-02 2.720e-01 7.800e-02 2.129e+05 1.131e+05]
RMSE 2.17995317290043
[7.600e-02 2.720e-01 7.800e-02 2.129e+05 1.131e+05]
RMSE 2.17995317290043
[7.600e-02 2.720e-01 7.800e-02 2.129e+05 1.131e+05]
RMSE 2.17995317290043
[7.60000011e-02 2.72000000e-01 7.80000000e-02 2.12900000e+05
1.13100000e+05]
RMSE 2.17995317290043
[7.60000000e-02 2.72000004e-01 7.80000000e-02 2.12900000e+05
1.13100000e+05]
RMSE 2.17995317290043
[7.60000000e-02 2.72000000e-01 7.80000012e-02 2.12900000e+05
1.13100000e+05]
RMSE 2.17995317290043
[7.60000000e-02 2.72000000e-01 7.80000000e-02 2.12900003e+05
1.13100000e+05]
RMSE 2.17995317290043
[7.60000000e-02 2.72000000e-01 7.80000000e-02 2.12900000e+05
1.13100002e+05]
RMSE 2.17995317290043

``````
``````

In :

details

``````
``````

Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

fvec
[-0.33610526315786693, -0.378326947375756, -0....

nfev
6

fjac
[[-0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0...

ipvt
[1, 2, 3, 4, 5]

qtf
[-0.33610526315786693, -0.378326947375756, -0....

cov_x
None

mesg
The cosine of the angle between func(x) and an...

ier
4

``````
``````

In :

details.mesg

``````
``````

Out:

'The cosine of the angle between func(x) and any column of the\n  Jacobian is at most 0.000000 in absolute value'

``````

The best params are only slightly different from the starting place.

``````

In :

best_params

``````
``````

Out:

array([7.600e-02, 2.720e-01, 7.800e-02, 2.129e+05, 1.131e+05])

``````

Here's what the results look like with the `best_params`.

``````

In :

system = make_system(best_params, data)
system = remove_units(system)

results, details = run_ode_solver(system, slope_func, t_eval=system.ts)
Q_frame = recompute_fluxes(results, system)
errors = compute_error(Q_frame, data)
print(np.sqrt(np.mean(errors**2)))

``````
``````

2.17995317290043

``````

The RMS error is only slightly smaller.

And the results are visually similar.

``````

In :

plot_Q_in(Q_frame, data)

``````
``````

``````
``````

In :

plot_Q_out(Q_frame, data)

``````
``````

``````

Exercise: Try starting the model with a different set of parameters and see if it moves toward the parameters in the paper.

I found that no matter where I start, `leastsq` doesn't move far, which suggests that it is not able to optimize the parameters effectively.

Notes on working with degC.

Usually I construct a `Quantity` object by multiplying a number and a unit. With degC, that doesn't work; you get

``OffsetUnitCalculusError: Ambiguous operation with offset unit (degC).``
``````

In :

#16.11 * C

``````

The problem is that it doesn't know whether you want a temperature measurement or a temperature difference.

You can create a temperature measurement like this.

``````

In :

T = Quantity(16.11, degC)

``````
``````

Out:

16.11 degC

``````

If you convert to Kelvin, it does the right thing.

``````

In :

T.to(K)

``````
``````

Out:

289.26 kelvin

``````

When you subtract temperatures, the results is a temperature difference, indicated by the units.

``````

In :

diff = T - 273.15 * K

``````
``````

Out:

16.11 delta_degC

``````

If you convert a temperature difference to Kelvin, it does the right thing.

``````

In :

diff.to(K)

``````
``````

Out:

16.11 kelvin

``````
``````

In [ ]:

``````