This is building on a different discourse post I made previously:
And on Github:
https://github.com/SciML/Catalyst.jl/issues/398
I am trying to solve a chemical system of equations that combines both equilibrium and kinetic species. Below is the sample code I am using:
using OrdinaryDiffEq
using Catalyst
using DifferentialEquations
using Plots
using Latexify
using ModelingToolkit
@parameters t
@variables H(t), OH(t), HOCl(t), NH3(t), NH2Cl(t), NHCl2(t), BrCl(t), OCl(t), OH(t), NH4(t), Br(t)
I’ve defined the equilibrium species:
pH2 = [
log10(H) + log10(OH) + 14 ~ 0,
log10((NH3*H)/NH4) + 9.25 ~ 0,
log10((OCl*H)/HOCl) + 7.54 ~ 0,
];
@named pHsystem = ODESystem(pH2,t,[H,OH,NH3,NH4,OCl,HOCl],[], );
I’ve also defined the kinetic species using Catalyst.jl and then coverted it to an ODE system
@parameters k1, k2, k3, k4, k5
disinfectionmodel = [
Reaction(k1, [HOCl, NH3], [NH2Cl]),
Reaction(k2, [NH2Cl], [HOCl, NH3]),
Reaction(k3, [HOCl, NH2Cl],[NHCl2]),
Reaction(k4, [NHCl2], [HOCl, NH2Cl]),
Reaction(k5, [HOCl, H, Br], [BrCl])
]
var= [H, OCl, HOCl, NH3, NH2Cl, NHCl2, BrCl, Br]
rates =[k1, k2, k3, k4, k5]
@named dis = ReactionSystem(disinfectionmodel, t, var, rates)
dis2= convert(ODESystem, dis)
I then combined the systems, simplified, and tried to define the ODEProblem
@named combined = extend(dis2, pHsystem)
flattenmodel = ModelingToolkit.flatten(combined)
sys=ModelingToolkit.alias_elimination(flattenmodel)
### Parameters [k1 k2 k3 k4 k5];
p = [4.2*10^6, 2.1*10^-5, 278, 6.4*10^-7, 1.32*10^6];
tspan = (0.0, 500.0);
##Initialization in order of appearance in reaction network
u0=[Hu, OClu, HOClu, NH3u, 0.0, 0.0, 0.0, Bru, OHu, NH4u];
alg = Tsit5();
op=ODEProblem(sys, u0, tspan, p)
I receive the following stacktrace:
ArgumentError: Equations (11), states (10), and initial conditions (10) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.
Stacktrace:
[1] check_eqs_u0(eqs::Vector{Equation}, dvs::Vector{Any}, u0::Vector{Float64}; check_length::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\abstractsystem.jl:854
[2] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Float64}, parammap::Vector{Float64}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:477
[3] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Float64}, tspan::Tuple{Float64, Float64}, parammap::Vector{Float64}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:557
If I try add the kwarg check_length=false, the stacktrace changes to this:
DimensionMismatch("new dimensions (10, 10) must be consistent with array size 121")
Stacktrace:
[1] (::Base.var"#throw_dmrsa#234")(dims::Tuple{Int64, Int64}, len::Int64)
@ Base .\reshapedarray.jl:41
[2] reshape
@ .\reshapedarray.jl:45 [inlined]
[3] reshape
@ .\reshapedarray.jl:116 [inlined]
[4] restructure(x::Matrix{Float64}, y::Matrix{Float64})
@ ArrayInterface ~\.julia\packages\ArrayInterface\61qJ7\src\ArrayInterface.jl:403
[5] (ODEFunction{true, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV} where {F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV})(sys::ODESystem, dvs::Vector{Any}, ps::Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}, u0::Vector{Float64}; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, eval_module::Module, steady_state::Bool, checkbounds::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{5, Symbol}, NamedTuple{(:ddvs, :linenumbers, :parallel, :has_difference, :check_length), Tuple{Nothing, Bool, Symbolics.SerialForm, Bool, Bool}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:260
[6] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Float64}, parammap::Vector{Float64}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:has_difference, :check_length), Tuple{Bool, Bool}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:479
[7] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Float64}, tspan::Tuple{Float64, Float64}, parammap::Vector{Float64}; kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:check_length,), Tuple{Bool}}})
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:557
[8] #ODEProblem#221
@ ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\diffeqs\abstractodesystem.jl:535 [inlined]
[9] top-level scope
@ In[18]:2
[10] eval
@ .\boot.jl:360 [inlined]
[11] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1094
I understand the current issue is that the system is overdefined with more equations than states. I’m not sure if I’m making a mistake setting up the problem or omitting something to reduce the system of equations equal to the unknowns (structural simplify didn’t work, I received the following stacktrace):
ExtraEquationsSystemException: The system is unbalanced. There are 10 highest order derivative variables and 11 equations.
More equations than variables, here are the potential extra equation(s):
0 ~ -7.54 - log10((H(t)*OCl(t)) / HOCl(t))
Stacktrace:
[1] error_reporting(sys::ODESystem, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
@ ModelingToolkit.StructuralTransformations ~\.julia\packages\ModelingToolkit\eSpDL\src\structural_transformation\utils.jl:64
[2] check_consistency(sys::ODESystem)
@ ModelingToolkit.StructuralTransformations ~\.julia\packages\ModelingToolkit\eSpDL\src\structural_transformation\utils.jl:103
[3] structural_simplify(sys::ODESystem; simplify::Bool)
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\abstractsystem.jl:815
[4] structural_simplify(sys::ODESystem)
@ ModelingToolkit ~\.julia\packages\ModelingToolkit\eSpDL\src\systems\abstractsystem.jl:814
[5] top-level scope
@ In[19]:4
[6] eval
@ .\boot.jl:360 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1094
Any help is greatly appreciated!