Change of variables in ModelingToolkit.jl

General Usage
,
1

Hi,

I am trying to estimate the parameters of a simple RC network in ModelingToolkit.jl.

image

The voltage over C can be defined by the following ODE:

\frac{dT_i}{dt} = \Phi_h(t) + \frac{T_a(t) - T_i(t)}{C_i R_{ia}}

Estimating C_i and R_{ia} turns out to be quite difficult due to the fact that R_{ia} and C_i can be small (~0.1) and the singularity present at 0.

I have found that thus replacing C_i and R_{ia} by their log and replacing \frac{1}{C_i R_{ia}} by \exp{(-\log(C) - \log(R))} significantly improves stability:

\frac{dT_i}{dt} = \Phi_h(t) + \exp{(-\log(C) - \log(R))}\left[T_a(t) - T_i(t)\right]

The question is, how to insert this into the model defined by modelingtoolkit?

NB: I am using ModelingToolkitStandardLibrary.jl, so essentially this model: RC Circuit Model

2

Hi!

How are you doing the parameter estimation exactly (I mean which package do you use)? If it is some ā€œexternalā€ tool (with respect to the ODE solution) then you could just wrap the ODE solution in another function that accepts the log-parameters and exponentiates them before solving.

Something like

function wrapped_model(logC, logR)
    C, R = exp(C), exp(R)
    return model(C, R)
end

where model is the something that solves the system you create with ModelingToolkit.jl

2 Likes
3

Thatā€™s indeed a good option, but it wouldnā€™t allow me to use the algebraic simplification. Iā€™m not sure how much that affects numerical stability, letā€™s try it

Ps. I am using Turing.jl

1 Like
4

Just add those equations to the system and then use the observation functionality so your interaction is in terms of those chosen variables. Can you share what you tried?

3 Likes
5

Yeah, thatā€™s the big downside. If the instability is related to the changed structure of the equation itself my suggestion wonā€™t help much. At least in my experience with slightly different ODEs it worked well to just wrap the model function (but I had no actual singularity in the model, just a constraint that some parameters need to stay positive ā€“ hence using the log variables helped the optimization algorithm to quickly hone in on the right orders of magnitude while automatically satisfying the constraint).

But I was also curious how to do this within the ModelingToolkit/ODE ecosystem itself. Thanks for the pointer @ChrisRackauckas !

1 Like
6

Thanks for your answer Chris. Iā€™m not entirely sure what you mean by observation functionality? I may need some pointers to the right docs.

Sure, here is the code:

Code

Components

@named Ci = Capacitor(C=5)
@named Ria = Resistor(R=1)
@named Qh = Current()
@named Ta = Voltage()
@named gnd = Ground()

Control signals

value_vector = 300 * randn(Int(t_end))
Ī¦_h(t) = max(t >= t_end ? value_vector[end] : value_vector[Int(floor(t)) + 1], 0)
@register_symbolic Ī¦_h(t)
@named f_h = TimeVaryingFunction(Ī¦_h)

T_bias = T_zeroC + 20 # [K]
T_a(t) = T_bias # + 20*sin((Ļ€)*t)
@register_symbolic T_a(t)
@named f_a = TimeVaryingFunction(T_a)

Connections

eqs = [
    # component connections
    ModelingToolkit.connect(Ci.p, Ria.p),
    ModelingToolkit.connect(Qh.n, Ci.p),
    ModelingToolkit.connect(Ta.p, Ria.n),
    ModelingToolkit.connect(gnd.g, Ci.n),
    ModelingToolkit.connect(gnd.g, Qh.p),
    ModelingToolkit.connect(gnd.g, Ta.n),

    # input to source connections
    ModelingToolkit.connect(f_h.output, Qh.I), 
    ModelingToolkit.connect(f_a.output, Ta.V),
];

mtk_components = [Ci, Ria, Qh, Ta, f_h, f_a, gnd];

Problem setup

@named model = ODESystem(eqs, t, systems = mtk_components)
sys = structural_simplify(model)

# initial conditions
u0 = [T_bias]

prob = ODEProblem(sys, u0, (0, t_end))

Generate test data

function generate_data(nr_samples::Int64, t_end::Float64, noise::Float64, f_input::Function, prob::ODEProblem)
    # sample points
    t_int = t_end/nr_samples
    t_sample = 0:t_int:t_end

    # sample input signal
    Ī¦_h_sample = f_input.(t_sample)

    # sample output signal
    d_sol = Vector(solve(prob, Tsit5(); saveat=t_int))
    d_noise = noise * randn(size(d_sol))
    T_i_sample = d_sol + d_noise

    return t_sample, Ī¦_h_sample, T_i_sample, T_a.(t_sample)
end

n_samples = Int64(floor((100/40) * t_end))
t_int = t_end/n_samples

# test data samples
t_sample, Ī¦_h_sample, T_i_sample, T_a_sample = generate_data(n_samples, t_end, 2.0, Ī¦_h, prob)

Turing.jl model

using Turing
using DifferentialEquations


@model function LGDS(x::AbstractVector, prob)
    # Prior distributions.
    Ļƒ ~ InverseGamma(2, 3)
    C ~ Gamma(1, 4)
    R ~ Gamma(1, 1)

    # ODE solution
    p = [R, C]
    T_ode = solve(prob, Tsit5(); p=p, saveat=t_int, save_idxs=1)
    if length(T_ode) != length(x)
        return
    end

    # observations
    x ~ MvNormal(T_ode.u, Ļƒ^2 * I)

    return nothing
end

model = LGDS(T_i_sample, prob)

# sample 3 independent chains
n_chains = 3
chain = sample(model, Turing.NUTS(), MCMCThreads(), 1000, n_chains; progress=true)

If you guys have any other suggestions Iā€™m all ears ofcourse.

The main problem with the Turing.jl model right now is that the chains converge to wildly different values. If I create 10 chains I will have 10 non-overlapping distributions.

7

Note that you donā€™t want to just directly use the .u because thatā€™s dependent on the chosen states in the model. If you instead use T_ode[[observable1, observable2]] then youā€™ll get the solution in terms of those observables. Observables can be defined as equations of the model variables and just added to the model, for example:

@variable observable1 observable2
obs = [observable1 ~ Ci.V^2 + Ci.i^2, observable2 ~ Ci.C - Ria.V]
@named model = ODESystem([eqs;obs], t, systems = mtk_components)
1 Like
8

Thanks, appreciate the help. Could you (or someone reading this) please give a hint as to what the observables should be for my problem? Iā€™m a bit lost.

I think I also need to add two parameters, logC and logR. Then I can sample them from a (standard) normal distribution and C and R will be log-normally distributed.

9

You can also try JuliaSimModelOptimizer (JuliaSimModelOptimizer: Automated Solution of Complex Inverse Problems ƂĀ· JuliaSimModelOptimizer). With the above example code, you can try something like

using JuliaSimModelOptimizer

data = DataFrame(:timestamp => t_sample, Symbol("Ta.i") => T_i_sample)

experiment = Experiment(data, sys)
invprob = InverseProblem(experiment, sys, [sys.Ria.R => (1e-3, 1e3, :log10), sys.Ci.C => (1e-3, 1e3, :log10)])

ps = parametric_uq(invprob, StochGlobalOpt(), sample_size=10)

As a quick explanation, the experiment above is used to encode that the generated data corresponds to the system and the column names for the data match the corresponding variables (which can be states or observed) in the model. In the invprob, the third argument describes what parameters to estimate based, their lower and upper bounds and to do the search in log space. The last line actually builds and solves the optimization problem that corresponds to the above sample_size times with initial conditions distributed via Latin hypercube sampling. I hope that this can be helpful :smiley:

P.S. I am involved in the development of JuliaSimModelOptimizer

10

Thanks, Iā€™ll try it out. I am also interested in estimating the uncertainty of my variable estimates since downstream optimizations will require it, would that be possible with this package?

11

Yes, see parametric_uq.