Calibration of a Dynamical System using Collocation Methods

Collocation Methods are techniques where the loss function minimizes the derivatives instead of the states or observed measurements in the data. This technique can be quite useful when the model is known to be unstable or takes a lot of time to simulate. This is because this suite of techniques does not involve simulating the model at all as derivatives are estimated from the data and derivatives for those time points are obtained from the model using the right hand side function. This also has a limitation - if we have partial observability, i.e., if we have data to only a few states or observed variables, there is no guarantee whether we can obtain the derivatives from the model. In this case, estimation of rest of the states from the data is internally handled and if it is not possible, it errors out gracefully.

In this tutorial, we will show how to use collocation methods for calibration of a dynamical system namely Continuous Stirred Tank Reactor. This dynamical system is a common model representing a chemical reactor.

Julia Environment

For this tutorial we will need the following packages:

using JuliaSimModelOptimizer
using ModelingToolkit
import ModelingToolkit: D_nounits as D, t_nounits as t
using ModelingToolkitStandardLibrary.Blocks: RealInput, TimeVaryingFunction
using OrdinaryDiffEq
using CSV, DataFrames
using DataSets
using DataInterpolations: AkimaInterpolation, CubicSpline
using Plots

Data Setup

We recorded some simulation data from a MPC Controller for this model using JuliaSimControl as illustrated here. Let us use this data for calibration.

cstr_dataset = dataset("juliasimtutorials/cstr_mpc_data")

data = open(IO, cstr_dataset) do io
    CSV.read(io, DataFrame)
end

first(data, 5)

Model Setup

f = CubicSpline(data[!, "model.F"], data[!, "timestamp"])
q = CubicSpline(data[!, "model.Q̇"], data[!, "timestamp"])

@variables Cₐ(t)=0.8 Cᵦ(t)=0.5 Tᵣ(t)=134.14 Tₖ(t)=130.0
@named finput = RealInput()
@named qinput = RealInput()
F = finput.u
Q̇ = qinput.u

ps = @parameters K0_ab=1.287e12 K0_bc=1.287e12 K0_ad=9.043e9 R_gas=8.31446e-3 E_A_ab=9758.3 E_A_bc=9758.3 E_A_ad=8560.0 Hᵣ_ab=4.2 Hᵣ_bc=-11.0 Hᵣ_ad=-41.85 Rou=0.9342 Cp=3.01 Cpₖ=2.0 Aᵣ=0.215 Vᵣ=1.0 m_k=2.5 T_in=130.0 K_w=4032.0 C_A0=(5.7+4.5)/2.0*1.0

eqs = [
    D(Cₐ) ~ F*(C_A0 - Cₐ)-(K0_ab * exp((-E_A_ab)/((Tᵣ+273.15))))*Cₐ - (K0_ad * exp((-E_A_ad)/((Tᵣ+273.15))))*abs2(Cₐ)
    D(Cᵦ) ~ -F*Cᵦ + (K0_ab * exp((-E_A_ab)/((Tᵣ+273.15))))*Cₐ - (K0_bc * exp((-E_A_bc)/((Tᵣ+273.15))))*Cᵦ
    D(Tᵣ) ~ (((K0_ab * exp((-E_A_ab)/((Tᵣ+273.15))))*Cₐ*Hᵣ_ab + (K0_bc * exp((-E_A_bc)/((Tᵣ+273.15))))*Cᵦ*Hᵣ_bc +
                (K0_ad * exp((-E_A_ad)/((Tᵣ+273.15))))*abs2(Cₐ)*Hᵣ_ad)/(-Rou*Cp)) + F*(T_in-Tᵣ) + (((K_w*Aᵣ)*(-(Tᵣ-Tₖ)))/(Rou*Cp*Vᵣ))
    D(Tₖ) ~ (Q̇ + K_w*Aᵣ*(Tᵣ-Tₖ))/(m_k*Cpₖ)
]

@named model = ODESystem(eqs, t, [Cₐ, Cᵦ, Tᵣ, Tₖ], ps, systems=[finput, qinput])

fsrc = TimeVaryingFunction(f; name = :F)
qsrc = TimeVaryingFunction(q; name = :Q̇)

eqs = [
    ModelingToolkit.connect(model.finput, fsrc.output),
    ModelingToolkit.connect(model.qinput, qsrc.output),
]

@named cstr = ODESystem(eqs, t; systems = [model, fsrc, qsrc])
sys = structural_simplify(complete(cstr))

Defining Experiment and InverseProblem

In order to create an Experiment, we will use the default initial values of the states and parameters of our model. These are our initial guesses which will be used to optimize the inverse problem in order to fit the given data.

experiment = Experiment(data, sys, loss_func = meansquaredl2loss)

Experiment for cstr 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, 0.245)

Once we have created the experiment, the next step is to create an InverseProblem. This inverse problem, requires us to provide the search space as a vector of pairs corresponding to the parameters that we want to recover and the assumption that we have for their respective bounds. Here we are trying to calibrate the values of m_k which is the mass of the coolant and Vᵣ which is the volume of the reactor.

prob = InverseProblem(experiment, [sys.model.m_k => (0.01, 10.0), sys.model.Vᵣ => (0.01, 20.0)])

InverseProblem with one experiment with 2 elements in the search space.

Calibration

SingleShooting

Let us first try to solve this problem using SingleShooting. To do this, we first define an algorithm alg and then call calibrate with the prob and alg.

alg = SingleShooting(maxiters = 1000)
r = calibrate(prob, alg)

Calibration result computed in 25 seconds and 723 milliseconds. Final objective value: 0.000586378.
Optimization ended with Success Return Code and returned.

┌───────────┬──────────┐
│ model₊m_k │ model₊Vᵣ │
├───────────┼──────────┤
│    4.9996 │   10.137 │
└───────────┴──────────┘

The true values of the parameters are - m_k has a value of 5.0 and Vᵣ has a value of 10.0. We can see that the calibrated parameters do not match the true values.

Collocation

Now, let us try to do calibration process using a collocation method. For this tutorial, we will use SplineCollocation with Akima Interpolation. To do this, we have to define the algorithm very similar to how we did with SingleShooting.

alg = SplineCollocation(maxtime = 100, interp = AkimaInterpolation, cutoff = (0.05, 0.0))

SplineCollocation method, optimizing with OptimizerWithAttributes.
maxiters = 1000 and maxtime = 100.0

Before we do calibration, we can visualize whether the estimated states match the data and estimated derivatives match from the model. This is an important step to check before calibration to visually verify whether the fit is accurate or not. To do this we use CollocationData, and call it with the algorithm, experiment and inverse problem we defined and also pass in the true value of parameters.

collocated_data = CollocationData(alg, experiment, prob, [5.0, 10.0])

CollocationData{Matrix{Float64}, Matrix{Float64}, SubArray{Float64, 2, Matrix{Float64}, Tuple{Vector{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, Matrix{Float64}, Vector{Float64}, Vector{Symbol}}([0.8 0.8052892673432875 … 0.7048304177472274 0.7048313103871084; 0.5 0.5426393075394137 … 0.600011095818671 0.6000098454780807; 134.14 133.95722321035586 … 140.3875192805388 140.41321977282766; 130.0 131.4078100362802 … 139.5705447235497 139.5891515963195], [0.8 0.8052892673432875 … 0.7048304177472274 0.7048313103871084; 0.5 0.5426393075394137 … 0.600011095818671 0.6000098454780807; 134.14 133.95722321035586 … 140.3875192805388 140.41321977282766; 130.0 131.4078100362802 … 139.5705447235497 139.5891515963195], [0.937234043923505 -7.691087460299609 … 0.04468692019304399 0.04048900961303037; 9.865361724911505 8.48138795566868 … -0.009133820515586422 -0.008292030580633991; -68.93781267896443 -2.133990863123465 … 5.312040707167938 4.810111715345272; 357.97996119058547 219.09492441530134 … 3.3781908729521772 2.881009440522007], [4.861675340721061 -4.718286695138286 … 0.00017277072479236581 0.00018473453717859903; 9.50355086651151 7.720606305991691 … -0.0002569768771312814 -0.00024305344258346274; -61.83744485913109 -2.720539531995406 … 5.395580494179455 4.877559717277788; 332.6378917951758 214.60324990095154 … 3.97747880826698 3.457772655238498], [0.0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045  …  0.2, 0.205, 0.21, 0.215, 0.22, 0.225, 0.23, 0.235, 0.24, 0.245], [Symbol("model.Cₐ"), Symbol("model.Cᵦ"), Symbol("model.Tᵣ"), Symbol("model.Tₖ")])

To plot it, we do:

For states,

plot(collocated_data, vars = "states", layout = (4, 1), size = (1000, 1000))

For derivatives,

plot(collocated_data, vars = "derivatives", layout = (4, 1), size = (1000, 1000))

We can see the derivatives match well except around the beginning. This is expected with collocation methods as they won't fit the data perfectly around the edges. To mitigate this, we use cutoff argument in the alg which tells us how much data to cut in the beginning and end for computing loss.

We can now proceed with calibration.

r = calibrate(prob, alg)

Calibration result computed in 2 seconds and 986 milliseconds. Final objective value: 1.51174.
Optimization ended with Success Return Code and returned.

┌───────────┬──────────┐
│ model₊m_k │ model₊Vᵣ │
├───────────┼──────────┤
│   5.00969 │  10.2065 │
└───────────┴──────────┘

We can see that the calibrated parameters are quite close to the true values!