-
-
Notifications
You must be signed in to change notification settings - Fork 195
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Sensitivity analysis transformation of ODEs #39
Comments
Calculating the Jacobian can be done via 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 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. |
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 |
That's the initial condition to the sensitivity part. But indeed, that looks correct. |
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 |
|
No, we cannot. That will make structural simplify fail. |
Why would structural simplify care about derivatives w.r.t things other than the independent variable? |
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. |
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.
The text was updated successfully, but these errors were encountered: