Global Model Calibration of Thermal Conduction

Thermal conduction model

Here we present an example of integrating the ModelingToolkitStandardLibrary with JuliaSimModelOptimizer by using the thermal conduction example. This example demonstrates the thermal response of two masses connected by a conducting element. The two masses have the same unknown heat capacity but different initial temperatures (T1=100 [°C], T2=0 [°C]). The mass with the higher temperature will cool off while the mass with the lower temperature heats up. They will each asymptotically approach the calculated temperature T_final that results from dividing the total initial energy in the system by the sum of the heat capacities of each element. The goal of this example is to find the unknown heat capacitance and also quantify its uncertainty.

Julia environment

For this example, we will need the following packages:

ModuleDescription
JuliaSimModelOptimizerThe high-level library used to formulate our problem and perform automated model discovery
ModelingToolkitThe symbolic modeling environment
ModelingToolkitStandardLibraryLibrary for using standard modeling components
OrdinaryDiffEqThe numerical differential equation solvers
DataFramesFor converting simulation into a Dataframe
PlotsThe plotting and visualization library
using JuliaSimModelOptimizer
using ModelingToolkit
using ModelingToolkitStandardLibrary.Thermal
using OrdinaryDiffEq
using DataFrames
using Plots

Model setup

@parameters t

C = 15.0
@named mass1 = HeatCapacitor(C = C, T = 373.15)
@named mass2 = HeatCapacitor(C = C, T = 273.15)
@named conduction = ThermalConductor(G = 10)
@named Tsensor1 = TemperatureSensor()
@named Tsensor2 = TemperatureSensor()

connections = [
    connect(mass1.port, conduction.port_a),
    connect(conduction.port_b, mass2.port),
    connect(mass1.port, Tsensor1.port),
    connect(mass2.port, Tsensor2.port),
]

@named model = ODESystem(connections, t, systems = [mass1, mass2, conduction, Tsensor1, Tsensor2])
sys = structural_simplify(model)

\[ \begin{align} \frac{\mathrm{d} mass1_{+}T\left( t \right)}{\mathrm{d}t} =& mass1_{+}der_{T}\left( t \right) \\ \frac{\mathrm{d} mass2_{+}T\left( t \right)}{\mathrm{d}t} =& mass2_{+}der_{T}\left( t \right) \end{align} \]

Data Setup

Now, lets simulate to generate some data.

prob = ODEProblem(sys, [mass1.C => 14.0, mass2.C => 14.0], (0, 100.0))
sol = solve(prob, Tsit5())
data = DataFrame(sol)
first(data, 5)
5×3 DataFrame
Rowtimestampmass1₊Tmass2₊T
Float64Float64Float64
10.0373.15273.15
20.125024364.972281.328
30.385166351.991294.309
40.715191341.149305.151
51.13432333.041313.259

Defining Experiment and InverseProblem

We next use the data to setup an Experiment.

experiment = Experiment(data, sys)
Experiment for model with no fixed parameters or initial conditions.
The simulation of this experiment is given by:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)

Next, we specify an InverseProblem of interest to us where we specify the parameters we want to recover and their bounds.

invprob = InverseProblem(experiment, sys, [mass1.C => (10, 20), mass2.C => (10, 20)])
InverseProblem with one experiment optimizing 2 parameters.

Parametric Uncertainty Quantification

We now run the parameteric uncertainty quantification using the parametric_uq function. We specify the InverseProblem and the StochGlobalOpt algorithm. We also specify the sample size for the analysis.

ps = parametric_uq(invprob, StochGlobalOpt(maxiters = 1000), sample_size = 50)
Parametric uncertainty ensemble of length 50 computed in 13 seconds and 445 milliseconds.
┌─────────┬─────────┐
│ mass1₊C │ mass2₊C │
├─────────┼─────────┤
│ 13.9659 │ 13.9828 │
│ 13.9642 │ 13.9788 │
│ 13.9349 │ 13.9437 │
│ 14.0141 │ 14.0082 │
│ 13.9338 │ 13.9213 │
│ 14.0005 │ 14.0065 │
│ 13.9477 │  13.963 │
│ 14.0673 │ 14.0708 │
│ 13.9735 │ 13.9778 │
│  14.012 │ 14.0053 │
│ 13.9227 │ 13.9212 │
│  13.976 │  13.985 │
│ 13.9825 │ 13.9836 │
│ 14.0206 │ 14.0204 │
│ 14.0147 │ 14.0199 │
│ 13.9924 │ 13.9971 │
│    ⋮    │    ⋮    │
└─────────┴─────────┘
      34 rows omitted

Visualization

We can see the results of the analysis by plotting its histogram.

plot(ps, layout = (2, 1), bins = 10, max_freq = [20, 20], legend = :best, size = (1000, 1000))
Example block output

We can see the distribution for the parameters is concentrated near their true values.