Sensitivity analysis transformation of ODEs

[ChrisRackauckas]

It would be nice to build the transformation equations for these kinds of things https://docs.sciml.ai/latest/extras/sensitivity_math/ directly from the IR since that would make it much faster by avoiding differencing.

[ChrisRackauckas]

Calculating the Jacobian can be done via

https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/systems/diffeqs/diffeqsystem.jl#L61-L76

which we know works now. All it does on the inside is:

https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/operators.jl#L107-L109

So we can create one of those for the gradients. Note that the syntax

https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/systems/diffeqs/diffeqsystem.jl#L70

lets you choose what variables to differentiate by, so we just need to create it generic like this and then the sensitivity analysis equations are this same Jacobian plus the parameter gradient.

[ArnoStrouwen]

Maybe this is useful:

using ModelingToolkit
using DifferentialEquations

@parameters t #time
@derivatives D'~t
@variables Cₛ(t), Cₓ(t), ν(t) # states
x = [Cₛ, Cₓ, ν]
n_x = length(x)
@parameters uᵢₙ # control input
u = [uᵢₙ]
n_u = length(u)
@parameters μₘₐₓ Cₛ_ᵢₙ # uncertain parameters
p = [μₘₐₓ, Cₛ_ᵢₙ]
n_p = length(p)
@parameters  Kₛ, yₓ_ₛ, m #other constants
μ = μₘₐₓ*Cₛ/Kₛ + Cₛ
σ = μ/yₓ_ₛ + m

eqs = [D(Cₓ) ~ -σ*Cₓ + uᵢₙ/ν*Cₛ_ᵢₙ - uᵢₙ/ν*Cₛ,
       D(Cₛ) ~ μ*Cₓ - uᵢₙ/ν*Cₓ,
       D(ν) ~ uᵢₙ]
sys = ODESystem(eqs)
function sensitivity_transform(sys,p)
    iv  = sys.iv()
    f = du_dt = [eq.rhs for eq ∈ sys.eqs]
    u = [state(iv) for state ∈ states(sys)] 
    sen_symbols = [Symbol("d",state.op.name,"/d",par.op.name) for state ∈ u, par in p]
    du_dp = [Variable(sen_symbol)(iv) for sen_symbol ∈ sen_symbols]
    @derivatives D'~iv
    uext= vcat(u,vec(du_dp))
    df_du = ModelingToolkit.jacobian(f,u)
    df_dp = ModelingToolkit.jacobian(f,p)
    ddu_dpdt = df_du*du_dp + df_dp
    duext_dt = fext = vcat(du_dt,vec(ddu_dpdt))
    eqsext = D.(uext) .~ fext
    sysext = ODESystem(eqsext)
end
sysext = sensitivity_transform(sys,p)

A bit dangerous is that u0/dpcan not be handled at this level.

[ChrisRackauckas]

A bit dangerous is that u0/dp can not be handled at this level.

That's the initial condition to the sensitivity part. But indeed, that looks correct.

[YingboMa]

using ModelingToolkit, OrdinaryDiffEq

@parameters t #time
@derivatives D'~t
@variables Cₛ(t), Cₓ(t), ν(t) # states
x = [Cₛ, Cₓ, ν]
n_x = length(x)
@parameters uᵢₙ # control input
u = [uᵢₙ]
n_u = length(u)
@parameters μₘₐₓ Cₛ_ᵢₙ # uncertain parameters
p = [μₘₐₓ, Cₛ_ᵢₙ]
n_p = length(p)
@parameters  Kₛ, yₓ_ₛ, m #other constants
μ = μₘₐₓ*Cₛ/Kₛ + Cₛ
σ = μ/yₓ_ₛ + m

eqs = [D(Cₓ) ~ -σ*Cₓ + uᵢₙ/ν*Cₛ_ᵢₙ - uᵢₙ/ν*Cₛ,
       D(Cₛ) ~ μ*Cₓ - uᵢₙ/ν*Cₓ,
       D(ν) ~ uᵢₙ]
sys = ODESystem(eqs)
function sensitivity_transform(sys,p)
    p = map(ModelingToolkit.value, p)
    iv  = independent_variable(sys)
    f = du_dt = [eq.rhs for eq in equations(sys)]
    u = states(sys)
    du_dp = [ModelingToolkit.rename(state, Symbol("d", ModelingToolkit.getname(state), "/d", ModelingToolkit.getname(par))) for state in u, par in p]
    @derivatives D'~iv
    uext= [u; vec(du_dp)]
    df_du = ModelingToolkit.jacobian(f, u)
    df_dp = ModelingToolkit.jacobian(f, p)
    ddu_dpdt = df_du*du_dp + df_dp
    duext_dt = fext = [du_dt; vec(ddu_dpdt)]
    eqsext = D.(uext) .~ fext
    sysext = ODESystem(eqsext, iv)
end
sysext = sensitivity_transform(sys, p)
prob = ODEProblem(sysext, [1; 2; 3; zeros(3*2)], (0, 1.0), collect(1:6))
solve(prob, Tsit5())

Works

[ChrisRackauckas]

sol[Differential(μₘₐₓ)(Cₓ)] fails because Symbol(Differential(μₘₐₓ)(Cₓ)) == Symbol("Cₓˍμₘₐₓ(t)") not var"dCₓ/dμₘₐₓ". Can we more simply support du_dp = Differential(par)(state)?

[YingboMa]

No, we cannot. That will make structural simplify fail.

[ChrisRackauckas]

Why would structural simplify care about derivatives w.r.t things other than the independent variable?

[YingboMa]

Ah, right, we can check the independent variable of the differential variable. This would make automatic system detection harder, though. We might want to force users to supply the independent variable.