# Global sensitivity analysis of Acausal circulation model - Sobol Method

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
# 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 :)``````

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

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

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?