Parameter substitutions in ModelingToolkit.jl

Specific Domains Modelling & Simulations
1

Hi,

I’m getting to grips with ModelingToolkit.jl. Let’s take the ODE building example in the tutorial:

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

First question: can I substitute the parameters in eqs with fixed constants? I couldn’t figure out how to do this.

The reason I’m asking this quesiton:
Suppose I have some fixed, nominal values for the tree parameters. I want to build a 7 state ODEProblem whose vector field looks like:

[eqs with fixed parameter values;
eqs with free parameter values;
l2 difference of the two eqs]

When I solve such a system, the final state vector at time t will give me the integral (over time) of the L2 difference between the Lorenz system with fixed parameter values, and with free parameter values. In particular, I’m hoping that the fixed parameter values enter the ODEProblem just as constants, so that there are only three free parameters in the ODEProblem

The question is, how do I substitute fixed parameter values into my array eqs? Or is there another way to generate an ODEFunction of the form I want?

I want to do this so I can efficiently take automatic derivatives of the L2 difference between the fixed=parameter and variable=parameter solutions, as a function of the variable parameters. This is equivalent to solving the adjoint sensitivity problem, where the cost functional depends on both the fixed-parameter and variable-parameter solutions.

Thanks in advance!

1 Like
2

More simply: given an equation with some variables, how do I substitute a different set of variables/constants in the equation?

3

We don’t have a function for that yet… @shashi have you thought of a way to implement it?

4

I would do it like this.

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

function  eqs(σ, ρ, β)
    [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
end

eqs(1, ρ, β)  # σ is known, etc.
3 Likes
5

I just want to preface this with the fact that I have almost no programming experience, and know way less than either of those who have replied, so definitely listen to them first and foremost. I just want to add what I did, which is more like substituting a value for your parameters instead of redefining the equation to redefine a parameter-value. I am continuing with your example of the lorentz attractor:

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
de    = ODESystem(eqs, t, [x, y, z], [σ, ρ, β])
ode_f = ODEFunction(de)
u₀    = [x₀, y₀, z₀]
tspan = (your_custom_timespan)
p     = [σ_val, ρ_val, β_val]
prob  = ODEProblem(ode_f, u₀, tspan, p)
sol   = solve(prob)

So when you define p with parameter values, you are substituting in your fixed constants. Another advantage, which is what I found useful about this way of defining things is that you can redefine parameters without redefining your original equations. So if you want to double the sigma value and run again, you can do

p    = [2σ_val, ρ_val, β_val]

or, even more compact,

p[1] = 2σ_val

and then run

prob  = ODEProblem(ode_f, u₀, tspan, p)
sol   = solve(prob)

again to get the new solution, all while using what you have already defined.

\n
\n

So that brings me to why I was googling and stumbled upon this post, and this is more directed at the two gentlemen who already replied, @shashi and @ChrisRackauckas (Thank you so much for your work, you are amazing):

The problem is that I feel like I have defined what the parameters are 3 times:
The “@parameters” macro,
In the equations,
In the “de” definition

and their values a seperate time, when defining p

I do realize that this is to get a nice structure and not necessary, but I don’t want to abandon the nice structure at all. I am imagining something like the following:

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
p     = Dict(σ=>σ_val, ρ=>ρ_val, β=>β_val) # constructing a dictionary linking parameters and their values
                                           # (maybe Pair is a better datatype with more natural syntax, I don't know)

de    = ODESystem(eqs, t, [x, y, z], p)    # this would read the dictionary keys
ode_f = ODEFunction(de)
u₀    = [x₀, y₀, z₀]
tspan = (your_custom_timespan)
prob  = ODEProblem(ode_f, u₀, tspan, p)    # This would index the dict keys to get values
sol   = solve(prob)

From here a redefinition is just as easy, and so the value of the nice structure is conserved:

p[σ]  = 2σ_val
prob  = ODEProblem(ode_f, u₀, tspan, p)
sol   = solve(prob)

Now I don’t that my suggestion is a perfect solution at all. Something maybe close to ideal would be to link parameters with values when they are declared in the macro (p=>p_val Pairs perhaps?), with the possibility of easy redefinition, and never having to pass them as arguments in “prob” and “de” definitions. I know you are doing a lot of amazing work and don’t want to waste your time on syntax-polishing if you have better things to do, so feel free to ignore this suggestion if it is hard to implement, unfeasible or plain stupid ^-^

6

And after some reflection I fialize the the @ode_def macro does pretty much everything I want. I have not seriously concidered it as it is part of ParameterizedFuntions, which will be depricated (right?). Is there a chance that that functionality will be moved to ModellingToolkit or DifferentialEquations?

7

No. In fact, I recently updated so its backend is ModelingToolkit and re-added it to DifferentialEquations . It’s not going to go anywhere because it works for what it does, and it does what it does well. I just don’t plan to add much to it, but feel free to continue using it if you like it. I probably won’t showcase it as much in future documentation though, but it is nice.

4 Likes
8

Awesome. Thanks

9

Thanks again for the helpful comments! I implemented a more general way of doing parameter transformations by hijacking the source code of modelingtoolkitize(). I’m putting it here for completeness.

Basically: defines a struct for any given transformation, with functions for both the transformation and the inverse transformation. I’ve provided examples for the log(abs(…)) transformation and the fix_parameters transformation. Given such a struct, and an ODEProblem, the function transform_problem(ODEProblem, TransformationStruct)…does what it says on the tin.

Note that I need the (inverse) transformation functions to operate both on Arrays of numbers, and arrays of operations. So e.g. the inverse transform for fixing parameters is coded to be agnostic to the type T of the array (T,1) input.

struct TransformationStructure{T<:Function,U<:Function}
    name::Union{String,Nothing}
    p_transform::T
    inv_p_transform::U
end

function logabs_transform(p0)
    """
    flip the signs of negative parameters, and then do the transform p -> log(p)
    """
    name = "logabs"
    is_positive = convert.(Int64, p0 .>= 0)
    is_positive[is_positive .== 0] .= -1
    pos(p) = p.* is_positive 
    p_transform(p) = log.(p.*is_positive)
    inv_p_transform(logp) = exp.(logp).*is_positive
    return TransformationStructure(name, p_transform, inv_p_transform)
end

function fix_params(p0, indices)
    name = nothing
    indices |> unique! |> sort!
    not_indices = setdiff(collect(1:length(p0)), indices)
    p0 = deepcopy(p0)
    function p_transform(p)
        return deleteat!(deepcopy(p), indices)
    end

    function inv_p_transform(p)
        out = [p[1] for el in p0]
        out[not_indices] .= p 
        out[indices] .= p0[indices]
        return out
    end
    return TransformationStructure(name, p_transform, inv_p_transform)
end



function transform_problem(prob::ODEProblem, tr::TransformationStructure; unames = nothing, pnames=nothing)
   

    @parameters t
    @derivatives D'~t

    newp0 = tr.p_transform(prob.p)

    unames === nothing && (unames = [Variable(:x,i)(t) for i in eachindex(prob.u0)])
    pnames === nothing && (pnames = [Variable(:p,i)() for i in eachindex(newp0)])
    
    pnames = transform_names!(pnames, tr)

    
    vars = reshape([Variable(Symbol(unames[i]))(t)   for i in eachindex(prob.u0)], size(prob.u0))
    
    params = reshape([Variable(Symbol(pnames[i]))()  for i in eachindex(newp0)], size(newp0))

    lhs = [D(var) for var in vars]
    if DiffEqBase.isinplace(prob)
        rhs = similar(vars, Any)
        prob.f(rhs, vars, tr.inv_p_transform(params), t)
    else
        rhs = prob.f(vars, tr.inv_p_transform(params), t)
    end

    eqs = vcat([lhs[i] ~ rhs[i] for i in eachindex(prob.u0)]...)
    de = ODESystem(eqs,t,vec(vars),vec(params))
    nprob = ODEProblem(de, vars .=> prob.u0, prob.tspan, params .=> newp0)
    return de, nprob
end


function transform_names(nv, tr::TransformationStructure)
    nv = tr.p_transform(nv)
    names = repr.(nv)
    nnv = [Variable(Symbol(names[i]))() for i in eachindex(names)]
    return nnv
end
2 Likes