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:

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)

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 464 milliseconds.
┌─────────┬─────────┐
│ mass1₊C │ mass2₊C │
├─────────┼─────────┤
│ 14.0739 │ 14.0755 │
│ 13.9747 │ 13.9968 │
│ 14.0436 │ 14.0237 │
│ 13.9796 │ 13.9773 │
│ 13.9921 │ 13.9861 │
│ 13.9698 │ 13.9698 │
│ 13.9831 │ 13.9874 │
│ 14.1027 │ 14.1016 │
│ 13.9991 │ 13.9881 │
│ 14.0407 │ 14.0341 │
│ 14.0128 │ 13.9945 │
│ 13.9941 │ 13.9869 │
│ 14.0339 │ 14.0357 │
│  13.998 │ 14.0034 │
│ 14.0071 │ 14.0125 │
│ 13.9858 │ 13.9895 │
│    ⋮    │    ⋮    │
└─────────┴─────────┘
      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))

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