Global sensitivity analysis of Acausal circulation model - Sobol Method

New to Julia
1

Hi there,

I am having trouble with the following code every time it appears my solver is not converging when I remake the ODE problem. Even when I change the value that I have set to infinity for something else I still get the following errors. I will provide both the code and stacktrace

using ModelingToolkit, DifferentialEquations, Plots, LaTeXStrings, Revise, DiffEqSensitivity, Statistics, OrdinaryDiffEq, LSODA

# CirculationModels is actually a module, but not a public one
# Set up your local modules and then I'll give you access to
# the github repo that I develop this on.

# Switch the following two lines to change from the included
# functions to the module, once you have it set up.
# using CirculationModels
include("CirculationModels.jl") 


include("parameters.jl");

nstep = 100000
t = LinRange(0, τ, nstep)
##
plotly()

## Start Modelling
@parameters t

## Named Elements

@named LV = shiChamber(V₀=10.0, τ=τ, Eₘᵢₙ=Eₘᵢₙ, Eₘₐₓ=Eₘₐₓ, τₑₛ=τₑₛ, τₑₚ=τₑₚ, Eshift=0.0, Ev=Inf)    
@named AV = ResistorDiode(R=1e-2)
@named MV = ResistorDiode(R=1e-2)

@named Rp = Resistor(R=R_p)
@named Rd = Resistor(R=R_d)

@named L1 = Inductance(L=L)

@named C1 = Capacitor(C=C)

@named ground = Ground(P=0.0)

## Connect the system 
circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, Rp.in, L1.in)
    connect(Rp.out, L1.out, C1.in, Rd.in)
    connect(Rd.out, MV.in)
    connect(MV.out, LV.in)
    connect(C1.out, ground.g)
    # Capacitor is connected to cavity
    # pressure which is same as 
    #C.out.p ~ 0.0
]

## Compose the whole ODAE system
@named _circ_model = ODESystem(circ_eqs, t)

##
@named circ_model = compose(_circ_model,
                          [LV, AV, MV, Rp, Rd, L1, C1, ground])

## And simplify it
circ_sys = structural_simplify(circ_model)



## Setup ODE 
prob = ODAEProblem(circ_sys, [70.0, 0.0, 0.0], (0.0, 20.0))


@time sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-12)

plot(sol, vars=[LV.p, Rp.in.p, Rd.out.p, C1.Δp, LV.in.q, LV.out.q], tspan=(18 * τ, 19 * τ))
##
plot(sol, vars=[AV.Δp, MV.Δp], tspan=(18 * τ, 19 * τ))
##

plot(sol, tspan=(18 * τ, 19 * τ) )

## Reparameterisation 

## Function wrapper for Reparameterisation

using ModelingToolkit:varmap_to_vars    


"""
reparameterize the ODEs defined in ModelingToolkit system [sys]
and [prob] with the new [parameter_map] and new initial
conditions [u0]
"""
function reparameterize(sys, prob, parameter_map, u0)
    # transform parameter map to array:
    params = varmap_to_vars(parameter_map, ModelingToolkit.parameters(sys))
    # redefine ODE problem with new parameters
    remake(prob; p=params, u0=u0)
end
##

params = varmap_to_vars([LV.V₀ => 10,LV.Eₘᵢₙ => Eₘᵢₙ,LV.Eₘₐₓ => Eₘₐₓ,LV.τ => 60/HR,LV.τₑₛ => τₑₛ,LV.τₑₚ => τₑₚ,LV.Eshift => 0.0,LV.Ev => 1.0e6,AV.R => 1e-3,MV.R => 1e-3,Rp.R => R_p,Rd.R => R_d,L1.L => L,C1.C => C,ground.P => 0.0], ModelingToolkit.parameters(circ_sys))




parameters(circ_sys)

##

equations(circ_sys)

## Parameter Map 
p = [
    LV.V₀ => 10
    LV.Eₘᵢₙ => Eₘᵢₙ
    LV.Eₘₐₓ => Eₘₐₓ
    LV.τ => 60/HR
    LV.τₑₛ => τₑₛ
    LV.τₑₚ => τₑₚ
    LV.Eshift => 0.0
    LV.Ev => Inf
    AV.R => 1e-3
    MV.R => 1e-3
    Rp.R => R_p
    Rd.R => R_d
    L1.L => L
    C1.C => C
    ground.P => 0.0
    ]

## Initial condition
u0 = [
    50
    0
    0
    ]

## And reparameterise
newprob = reparameterize(circ_sys, prob, p, u0)

## Solve it
@time newsol = solve(newprob, Tsit5(), reltol=1e-6, abstol=1e-12)



savetime = LinRange(19, 20, 100)

function circ(params)
    newprob = remake(prob, p=params)
    sol = solve(newprob,Tsit5(); saveat=savetime, reltol=1e-6, abstol=1e-12)
    [mean(sol[LV.p]),maximum(sol[LV.p])]
end

@time sobol_result = gsa(circ,Sobol(),[[8.0, 12.0], [1.0, 1.5], [2.0, 3.0], [0.8, 1.0], [0.3, 0.5], [0.4, 0.6], [0.0, 0.1], [10000000, 10000000.1], [0.0001, 0.001], [0.0001, 0.001], [0.25, 0.35], [0.9, 1.1], [0.003, 0.005], [0.05, 0.15], [0.0, 0.1]], N=1000);
# GSA does not like Inf for bounds 
#This will give two matrices, ST and S1, which have the Total Order Indices (TO) and the First Order Indices (FO), respectively.

sobol_result.ST
sobol_result.S1

And the Error it produces

ERROR: LoadError: The nonlinear solver failed with the return code MAXITERS_EXCEED.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] nlsolve_failure(rc::Symbol)
    @ ModelingToolkit.StructuralTransformations C:\Users\harry\.julia\packages\ModelingToolkit\r5ZaU\src\structural_transformation\utils.jl:259
  [3] numerical_nlsolve
    @ C:\Users\harry\.julia\packages\ModelingToolkit\r5ZaU\src\structural_transformation\utils.jl:265 [inlined]
  [4] macro expansion
    @ C:\Users\harry\.julia\packages\ModelingToolkit\r5ZaU\src\structural_transformation\codegen.jl:172 [inlined]
  [5] macro expansion
    @ C:\Users\harry\.julia\packages\SymbolicUtils\0KTj4\src\code.jl:304 [inlined]
  [6] macro expansion
    @ C:\Users\harry\.julia\packages\RuntimeGeneratedFunctions\KrkGo\src\RuntimeGeneratedFunctions.jl:129 [inlined]
  [7] macro expansion
    @ .\none:0 [inlined]
  [8] generated_callfunc
    @ .\none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)})(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Float64)
    @ RuntimeGeneratedFunctions C:\Users\harry\.julia\packages\RuntimeGeneratedFunctions\KrkGo\src\RuntimeGeneratedFunctions.jl:117
 [10] ODEFunction
    @ C:\Users\harry\.julia\packages\SciMLBase\Luo0F\src\scimlfunctions.jl:334 [inlined]
 [11] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\harry\.julia\packages\OrdinaryDiffEq\WD8cC\src\perform_step\low_order_rk_perform_step.jl:689      
 [12] perform_step!
    @ C:\Users\harry\.julia\packages\OrdinaryDiffEq\WD8cC\src\perform_step\low_order_rk_perform_step.jl:675 [inlined]
 [13] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq C:\Users\harry\.julia\packages\OrdinaryDiffEq\WD8cC\src\solve.jl:478
 [14] #__solve#496
    @ C:\Users\harry\.julia\packages\OrdinaryDiffEq\WD8cC\src\solve.jl:5 [inlined]
 [15] #solve_call#37
    @ C:\Users\harry\.julia\packages\DiffEqBase\8q10H\src\solve.jl:61 [inlined]
 [16] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#381"), Symbol("##arg#382"), Symbol("##arg#383"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x025f90df, 0x00e3b476, 0xc637b9ee, 0x8d0a1d09, 0xa1943bca)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#62"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :reltol, :abstol), Tuple{LinRange{Float64}, Float64, Float64}}}) 
    @ DiffEqBase C:\Users\harry\.julia\packages\DiffEqBase\8q10H\src\solve.jl:87
 [17] #solve#38
    @ C:\Users\harry\.julia\packages\DiffEqBase\8q10H\src\solve.jl:73 [inlined]
 [18] circ(params::Vector{Float64})
    @ Main c:\Users\harry\Desktop\PhD\Year 1\Sem 1\Valve and Ventricle analysis\AcausalModel.jl:145
 [19] (::GlobalSensitivity.var"#27#28"{typeof(circ)})(i::Int64)
    @ GlobalSensitivity .\none:0
 [20] iterate
    @ .\generator.jl:47 [inlined]
 [21] collect_to!(dest::Vector{Vector{Float64}}, itr::Base.Generator{UnitRange{Int64}, GlobalSensitivity.var"#27#28"{typeof(circ)}}, offs::Int64, st::Int64)
    @ Base .\array.jl:728
 [22] collect_to_with_first!(dest::Vector{Vector{Float64}}, v1::Vector{Float64}, itr::Base.Generator{UnitRange{Int64}, GlobalSensitivity.var"#27#28"{typeof(circ)}}, st::Int64)
    @ Base .\array.jl:706
 [23] collect(itr::Base.Generator{UnitRange{Int64}, GlobalSensitivity.var"#27#28"{typeof(circ)}})
    @ Base .\array.jl:687
 [24] gsa(f::typeof(circ), method::Sobol, A::Matrix{Float64}, B::Matrix{Float64}; batch::Bool, Ei_estimator::Symbol, distributed::Val{false}, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ GlobalSensitivity C:\Users\harry\.julia\packages\GlobalSensitivity\U9GZJ\src\sobol_sensitivity.jl:54
 [25] gsa(f::Function, method::Sobol, A::Matrix{Float64}, B::Matrix{Float64})
    @ GlobalSensitivity C:\Users\harry\.julia\packages\GlobalSensitivity\U9GZJ\src\sobol_sensitivity.jl:28
 [26] gsa(f::Function, method::Sobol, p_range::Vector{Vector{Float64}}; N::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ GlobalSensitivity C:\Users\harry\.julia\packages\GlobalSensitivity\U9GZJ\src\sobol_sensitivity.jl:187
 [27] gsa(f::Function, method::Sobol, args::Vector{Vector{Float64}}; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, 
NamedTuple{(:N,), Tuple{Int64}}})
    @ DiffEqSensitivity C:\Users\harry\.julia\packages\DiffEqSensitivity\xMT4T\src\DiffEqSensitivity.jl:75
 [28] top-level scope
    @ .\timing.jl:210 [inlined]
 [29] top-level scope
    @ c:\Users\harry\Desktop\PhD\Year 1\Sem 1\Valve and Ventricle analysis\AcausalModel.jl:0
in expression starting at c:\Users\harry\Desktop\PhD\Year 1\Sem 1\Valve and Ventricle analysis\AcausalModel.jl:149
```

Any help would be greatly appreciated I am new to Julia so just looking to improve :)
2

Did you double check the parameter ordering with the bounds? If you run the solve separate from the GSA, do you have an issue?

3

The bounds in the GSA match the parameter map p above so I don’t think the issue is there when using the remake it does through throw an error however when using the function reparametrise to remake the ODE the solve runs absolutely fine.

Thanks for the help :slight_smile:

4

So it’s impossible to isolate the parameters for which this occurs? If you do @show prob.p right before the solve and then separately solve with those parameters it’s fine?