From 207adaaa269329ae6c3935ca5793c28d6263101c Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 10 Jul 2023 15:39:03 -0400 Subject: [PATCH 1/7] Add DDEs support in `structural_simplify` --- src/systems/systemstructure.jl | 8 ++++++- test/dde.jl | 44 ++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) create mode 100644 test/dde.jl diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 4a9730ac80..7a719cf674 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -253,6 +253,7 @@ end function TearingState(sys; quick_cancel = false, check = true) sys = flatten(sys) ivs = independent_variables(sys) + iv = only(ivs) eqs = copy(equations(sys)) neqs = length(eqs) dervaridxs = OrderedSet{Int}() @@ -287,10 +288,15 @@ function TearingState(sys; quick_cancel = false, check = true) isalgeq = true statevars = [] for var in vars + if istree(var) && !ModelingToolkit.isoperator(var, Symbolics.Operator) + args = arguments(var) + length(args) == 1 || continue + isequal(args[1], iv) || continue + end set_incidence = true @label ANOTHER_VAR _var, _ = var_from_nested_derivative(var) - any(isequal(_var), ivs) && continue + isequal(_var, iv) && continue if isparameter(_var) || (istree(_var) && isparameter(operation(_var)) || isconstant(_var)) continue diff --git a/test/dde.jl b/test/dde.jl new file mode 100644 index 0000000000..ca4ab7861e --- /dev/null +++ b/test/dde.jl @@ -0,0 +1,44 @@ +using ModelingToolkit, DelayDiffEq +p0 = 0.2; +q0 = 0.3; +v0 = 1; +d0 = 5; +p1 = 0.2; +q1 = 0.3; +v1 = 1; +d1 = 1; +d2 = 1; +beta0 = 1; +beta1 = 1; +tau = 1; +function bc_model(du, u, h, p, t) + du[1] = (v0 / (1 + beta0 * (h(p, t - tau)[3]^2))) * (p0 - q0) * u[1] - d0 * u[1] + du[2] = (v0 / (1 + beta0 * (h(p, t - tau)[3]^2))) * (1 - p0 + q0) * u[1] + + (v1 / (1 + beta1 * (h(p, t - tau)[3]^2))) * (p1 - q1) * u[2] - d1 * u[2] + du[3] = (v1 / (1 + beta1 * (h(p, t - tau)[3]^2))) * (1 - p1 + q1) * u[2] - d2 * u[3] +end +lags = [tau] +h(p, t) = ones(3) +tspan = (0.0, 10.0) +u0 = [1.0, 1.0, 1.0] +prob = DDEProblem(bc_model, u0, h, tspan, constant_lags = lags) +alg = MethodOfSteps(Tsit5()) +sol = solve(prob, alg) + +@parameters p0=0.2 p1=0.2 q0=0.3 q1=0.3 v0=1 v1=1 d0=5 d1=1 d2=1 beta0=1 beta1=1 tau=1 +@variables t x₀(t) x₁(t) x₂(..) +D = Differential(t) +eqs = [D(x₀) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (p0 - q0) * x₀ - d0 * x₀ + D(x₁) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (1 - p0 + q0) * x₀ + + (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (p1 - q1) * x₁ - d1 * x₁ + D(x₂(t)) ~ (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (1 - p1 + q1) * x₁ - d2 * x₂(t)] +@named sys = System(eqs) +h(p, t) = ones(3) +tspan = (0.0, 10.0) +prob = DDEProblem(sys, + [x₀ => 1.0, x₁ => 1.0, x₂(t) => 1.0], + h, + tspan, + constant_lags = [tau]) +alg = MethodOfSteps(Tsit5()) +sol = solve(prob, alg) From 023d3709b0194412161276aa5d860934fcab0409 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 10 Jul 2023 16:59:20 -0400 Subject: [PATCH 2/7] Add DDE lowering --- src/systems/diffeqs/abstractodesystem.jl | 139 +++++++++++++++++++++-- src/systems/diffeqs/odesystem.jl | 5 + src/systems/systemstructure.jl | 6 +- test/dde.jl | 8 +- 4 files changed, 139 insertions(+), 19 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index af23769cef..d4a9a2112e 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -120,9 +120,14 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param implicit_dae = false, ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) : nothing, + isdde = false, has_difference = false, kwargs...) - eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)] + if isdde + eqs = delay_to_function(sys) + else + eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)] + end if !implicit_dae check_operator_variables(eqs, Differential) check_lhs(eqs, Differential, Set(dvs)) @@ -136,15 +141,58 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param p = map(x -> time_varying_as_func(value(x), sys), ps) t = get_iv(sys) - pre, sol_states = get_substitutions_and_solved_states(sys, - no_postprocess = has_difference) + if isdde + build_function(rhss, u, DDE_HISTORY_FUN, p, t; kwargs...) + else + pre, sol_states = get_substitutions_and_solved_states(sys, + no_postprocess = has_difference) - if implicit_dae - build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre, states = sol_states, - kwargs...) + if implicit_dae + build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre, + states = sol_states, + kwargs...) + else + build_function(rhss, u, p, t; postprocess_fbody = pre, states = sol_states, + kwargs...) + end + end +end + +function isdelay(var, iv) + isvariable(var) || return false + if istree(var) && !ModelingToolkit.isoperator(var, Symbolics.Operator) + args = arguments(var) + length(args) == 1 || return false + isequal(args[1], iv) || return true + end + return false +end +const DDE_HISTORY_FUN = Sym{Symbolics.FnType{Tuple{Any, <:Real}, Vector{Real}}}(:___history___) +function delay_to_function(sys::AbstractODESystem) + delay_to_function(full_equations(sys), + get_iv(sys), + Dict{Any, Int}(operation(s) => i for (i, s) in enumerate(states(sys))), + parameters(sys), + DDE_HISTORY_FUN) +end +function delay_to_function(eqs::Vector{<:Equation}, iv, sts, ps, h) + delay_to_function.(eqs, (iv,), (sts,), (ps,), (h,)) +end +function delay_to_function(eq::Equation, iv, sts, ps, h) + delay_to_function(eq.lhs, iv, sts, ps, h) ~ delay_to_function(eq.rhs, iv, sts, ps, h) +end +function delay_to_function(expr, iv, sts, ps, h) + if isdelay(expr, iv) + v = operation(expr) + time = arguments(expr)[1] + idx = sts[v] + return term(getindex, h(Sym{Any}(:ˍ₋arg3), time), idx, type = Real) # BIG BIG HACK + elseif istree(expr) + return similarterm(expr, + operation(expr), + map(x -> delay_to_function(x, iv, sts, ps, h), arguments(expr))) else - build_function(rhss, u, p, t; postprocess_fbody = pre, states = sol_states, - kwargs...) + return expr end end @@ -485,6 +533,30 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys), observed = observedfun) end +function DiffEqBase.DDEFunction(sys::AbstractODESystem, args...; kwargs...) + DDEFunction{true}(sys, args...; kwargs...) +end + +function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = states(sys), + ps = parameters(sys), u0 = nothing; + eval_module = @__MODULE__, + checkbounds = false, + kwargs...) where {iip} + f_gen = generate_function(sys, dvs, ps; isdde = true, + expression = Val{true}, + expression_module = eval_module, checkbounds = checkbounds, + kwargs...) + f_oop, f_iip = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) + f(u, p, h, t) = f_oop(u, p, h, t) + f(du, u, p, h, t) = f_iip(du, u, p, h, t) + + DDEFunction{iip}(f, + sys = sys, + syms = Symbol.(dvs), + indepsym = Symbol(get_iv(sys)), + paramsyms = Symbol.(ps)) +end + """ ```julia ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys), @@ -577,7 +649,7 @@ end """ u0, p, defs = get_u0_p(sys, u0map, parammap; use_union=false, tofloat=!use_union) -Take dictionaries with initial conditions and parameters and convert them to numeric arrays `u0` and `p`. Also return the merged dictionary `defs` containing the entire operating point. +Take dictionaries with initial conditions and parameters and convert them to numeric arrays `u0` and `p`. Also return the merged dictionary `defs` containing the entire operating point. """ function get_u0_p(sys, u0map, parammap; use_union = false, tofloat = !use_union) eqs = equations(sys) @@ -802,6 +874,55 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan end end +function DiffEqBase.DDEProblem(sys::AbstractODESystem, args...; kwargs...) + DDEProblem{true}(sys, args...; kwargs...) +end +function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [], + h = (u, p) -> zeros(length(states(sts))), + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + callback = nothing, + check_length = true, + kwargs...) where {iip} + has_difference = any(isdifferenceeq, equations(sys)) + f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, + has_difference = has_difference, + check_length, kwargs...) + cbs = process_events(sys; callback, has_difference, kwargs...) + if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing + affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...) + discrete_cbs = map(affects, clocks, svs) do affect, clock, sv + if clock isa Clock + PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt) + else + error("$clock is not a supported clock type.") + end + end + if cbs === nothing + if length(discrete_cbs) == 1 + cbs = only(discrete_cbs) + else + cbs = CallbackSet(discrete_cbs...) + end + else + cbs = CallbackSet(cbs, discrete_cbs) + end + else + svs = nothing + end + kwargs = filter_kwargs(kwargs) + + kwargs1 = (;) + if cbs !== nothing + kwargs1 = merge(kwargs1, (callback = cbs,)) + end + if svs !== nothing + kwargs1 = merge(kwargs1, (disc_saved_values = svs,)) + end + DDEProblem{iip}(f, u0, h, tspan, p; kwargs1..., kwargs...) +end + """ ```julia ODEProblemExpr{iip}(sys::AbstractODESystem, u0map, tspan, diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index cc60ac28fb..4195d6fce4 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -181,6 +181,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; checks = true, metadata = nothing, gui_metadata = nothing) + dvs = filter(x -> !isdelay(x, iv), dvs) name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) deqs = scalarize(deqs) @@ -258,6 +259,10 @@ function ODESystem(eqs, iv = nothing; kwargs...) push!(algeeq, eq) end end + for v in allstates + isdelay(v, iv) || continue + collect_vars!(allstates, ps, arguments(v)[1], iv) + end algevars = setdiff(allstates, diffvars) # the orders here are very important! return ODESystem(Equation[diffeq; algeeq; compressed_eqs], iv, diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 7a719cf674..0e1f9d7d10 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -288,11 +288,7 @@ function TearingState(sys; quick_cancel = false, check = true) isalgeq = true statevars = [] for var in vars - if istree(var) && !ModelingToolkit.isoperator(var, Symbolics.Operator) - args = arguments(var) - length(args) == 1 || continue - isequal(args[1], iv) || continue - end + ModelingToolkit.isdelay(var, iv) && continue set_incidence = true @label ANOTHER_VAR _var, _ = var_from_nested_derivative(var) diff --git a/test/dde.jl b/test/dde.jl index ca4ab7861e..1b943235cb 100644 --- a/test/dde.jl +++ b/test/dde.jl @@ -1,4 +1,4 @@ -using ModelingToolkit, DelayDiffEq +using ModelingToolkit, DelayDiffEq, Test p0 = 0.2; q0 = 0.3; v0 = 1; @@ -33,12 +33,10 @@ eqs = [D(x₀) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (p0 - q0) * x₀ - d0 (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (p1 - q1) * x₁ - d1 * x₁ D(x₂(t)) ~ (v1 / (1 + beta1 * (x₂(t - tau)^2))) * (1 - p1 + q1) * x₁ - d2 * x₂(t)] @named sys = System(eqs) -h(p, t) = ones(3) -tspan = (0.0, 10.0) prob = DDEProblem(sys, [x₀ => 1.0, x₁ => 1.0, x₂(t) => 1.0], h, tspan, constant_lags = [tau]) -alg = MethodOfSteps(Tsit5()) -sol = solve(prob, alg) +sol_mtk = solve(prob, alg) +@test Array(sol_mtk) ≈ Array(sol) From 8e422ee897e3285412f7a56f1f1ff98166bf423a Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 11 Jul 2023 11:24:07 -0400 Subject: [PATCH 3/7] Fix TearingState construction --- src/systems/systemstructure.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 0e1f9d7d10..4e330619d5 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -253,7 +253,7 @@ end function TearingState(sys; quick_cancel = false, check = true) sys = flatten(sys) ivs = independent_variables(sys) - iv = only(ivs) + iv = length(ivs) == 1 ? ivs[1] : nothing eqs = copy(equations(sys)) neqs = length(eqs) dervaridxs = OrderedSet{Int}() @@ -292,7 +292,7 @@ function TearingState(sys; quick_cancel = false, check = true) set_incidence = true @label ANOTHER_VAR _var, _ = var_from_nested_derivative(var) - isequal(_var, iv) && continue + any(isequal(_var), ivs) && continue if isparameter(_var) || (istree(_var) && isparameter(operation(_var)) || isconstant(_var)) continue From 7d594a1a72dfc619172834128926604bac69d869 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 11 Jul 2023 16:34:35 -0400 Subject: [PATCH 4/7] Symbolic support for DDE histories --- src/systems/diffeqs/abstractodesystem.jl | 25 ++++++++++++++++++++---- test/dde.jl | 21 ++++++++++++++------ 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index d4a9a2112e..4f8bb2bd59 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -651,7 +651,12 @@ end Take dictionaries with initial conditions and parameters and convert them to numeric arrays `u0` and `p`. Also return the merged dictionary `defs` containing the entire operating point. """ -function get_u0_p(sys, u0map, parammap; use_union = false, tofloat = !use_union) +function get_u0_p(sys, + u0map, + parammap; + use_union = false, + tofloat = !use_union, + symbolic_u0 = false) eqs = equations(sys) dvs = states(sys) ps = parameters(sys) @@ -660,7 +665,11 @@ function get_u0_p(sys, u0map, parammap; use_union = false, tofloat = !use_union) defs = mergedefaults(defs, parammap, ps) defs = mergedefaults(defs, u0map, dvs) - u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true) + if symbolic_u0 + u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = false, use_union = false) + else + u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true) + end p = varmap_to_vars(parammap, ps; defaults = defs, tofloat, use_union) p = p === nothing ? SciMLBase.NullParameters() : p u0, p, defs @@ -676,13 +685,14 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; eval_expression = true, use_union = false, tofloat = !use_union, + symbolic_u0 = false, kwargs...) eqs = equations(sys) dvs = states(sys) ps = parameters(sys) iv = get_iv(sys) - u0, p, defs = get_u0_p(sys, u0map, parammap; tofloat, use_union) + u0, p, defs = get_u0_p(sys, u0map, parammap; tofloat, use_union, symbolic_u0) if implicit_dae && du0map !== nothing ddvs = map(Differential(iv), dvs) @@ -874,11 +884,14 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan end end +function generate_history(sys::AbstractODESystem, u0; kwargs...) + build_function(u0, parameters(sys), get_iv(sys); expression = Val{false}, kwargs...) +end + function DiffEqBase.DDEProblem(sys::AbstractODESystem, args...; kwargs...) DDEProblem{true}(sys, args...; kwargs...) end function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [], - h = (u, p) -> zeros(length(states(sts))), tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); callback = nothing, @@ -888,7 +901,11 @@ function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [], f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, has_difference = has_difference, + symbolic_u0 = true, check_length, kwargs...) + h_oop, h_iip = generate_history(sys, u0) + h = h_oop + u0 = h(p, tspan[1]) cbs = process_events(sys; callback, has_difference, kwargs...) if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...) diff --git a/test/dde.jl b/test/dde.jl index 1b943235cb..ef546d4ce8 100644 --- a/test/dde.jl +++ b/test/dde.jl @@ -19,14 +19,18 @@ function bc_model(du, u, h, p, t) end lags = [tau] h(p, t) = ones(3) +h2(p, t) = ones(3) .- t * q1 tspan = (0.0, 10.0) u0 = [1.0, 1.0, 1.0] prob = DDEProblem(bc_model, u0, h, tspan, constant_lags = lags) -alg = MethodOfSteps(Tsit5()) -sol = solve(prob, alg) +alg = MethodOfSteps(Vern9()) +sol = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) +prob2 = DDEProblem(bc_model, u0, h2, tspan, constant_lags = lags) +sol2 = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) -@parameters p0=0.2 p1=0.2 q0=0.3 q1=0.3 v0=1 v1=1 d0=5 d1=1 d2=1 beta0=1 beta1=1 tau=1 +@parameters p0=0.2 p1=0.2 q0=0.3 q1=0.3 v0=1 v1=1 d0=5 d1=1 d2=1 beta0=1 beta1=1 @variables t x₀(t) x₁(t) x₂(..) +tau = 1 D = Differential(t) eqs = [D(x₀) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (p0 - q0) * x₀ - d0 * x₀ D(x₁) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (1 - p0 + q0) * x₀ + @@ -35,8 +39,13 @@ eqs = [D(x₀) ~ (v0 / (1 + beta0 * (x₂(t - tau)^2))) * (p0 - q0) * x₀ - d0 @named sys = System(eqs) prob = DDEProblem(sys, [x₀ => 1.0, x₁ => 1.0, x₂(t) => 1.0], - h, tspan, constant_lags = [tau]) -sol_mtk = solve(prob, alg) -@test Array(sol_mtk) ≈ Array(sol) +sol_mtk = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) +@test sol_mtk.u[end] ≈ sol.u[end] +prob2 = DDEProblem(sys, + [x₀ => 1.0 - t * q1, x₁ => 1.0 - t * q1, x₂(t) => 1.0 - t * q1], + tspan, + constant_lags = [tau]) +sol2_mtk = solve(prob2, alg, reltol = 1e-7, abstol = 1e-10) +@test sol2_mtk.u[end]≈sol2.u[end] atol=1e-5 From 5b678b290a763b2f3249a02a611d7e19eaceb8f0 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 11 Jul 2023 16:38:47 -0400 Subject: [PATCH 5/7] Fix test --- src/systems/diffeqs/odesystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 4195d6fce4..ab66628a54 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -181,16 +181,16 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; checks = true, metadata = nothing, gui_metadata = nothing) - dvs = filter(x -> !isdelay(x, iv), dvs) name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) deqs = scalarize(deqs) @assert all(control -> any(isequal.(control, ps)), controls) "All controls must also be parameters." iv′ = value(scalarize(iv)) - dvs′ = value.(scalarize(dvs)) ps′ = value.(scalarize(ps)) ctrl′ = value.(scalarize(controls)) + dvs′ = value.(scalarize(dvs)) + dvs′ = filter(x -> !isdelay(x, iv), dvs′) if !(isempty(default_u0) && isempty(default_p)) Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", From 045782f2d30ac229f7913ee5452a9ff3e91816bb Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 11 Jul 2023 16:42:42 -0400 Subject: [PATCH 6/7] Fix tests --- test/dde.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/dde.jl b/test/dde.jl index ef546d4ce8..aad39fe27f 100644 --- a/test/dde.jl +++ b/test/dde.jl @@ -19,14 +19,14 @@ function bc_model(du, u, h, p, t) end lags = [tau] h(p, t) = ones(3) -h2(p, t) = ones(3) .- t * q1 +h2(p, t) = ones(3) .- t * q1 * 10 tspan = (0.0, 10.0) u0 = [1.0, 1.0, 1.0] prob = DDEProblem(bc_model, u0, h, tspan, constant_lags = lags) alg = MethodOfSteps(Vern9()) sol = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) prob2 = DDEProblem(bc_model, u0, h2, tspan, constant_lags = lags) -sol2 = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) +sol2 = solve(prob2, alg, reltol = 1e-7, abstol = 1e-10) @parameters p0=0.2 p1=0.2 q0=0.3 q1=0.3 v0=1 v1=1 d0=5 d1=1 d2=1 beta0=1 beta1=1 @variables t x₀(t) x₁(t) x₂(..) @@ -44,8 +44,8 @@ prob = DDEProblem(sys, sol_mtk = solve(prob, alg, reltol = 1e-7, abstol = 1e-10) @test sol_mtk.u[end] ≈ sol.u[end] prob2 = DDEProblem(sys, - [x₀ => 1.0 - t * q1, x₁ => 1.0 - t * q1, x₂(t) => 1.0 - t * q1], + [x₀ => 1.0 - t * q1 * 10, x₁ => 1.0 - t * q1 * 10, x₂(t) => 1.0 - t * q1 * 10], tspan, constant_lags = [tau]) sol2_mtk = solve(prob2, alg, reltol = 1e-7, abstol = 1e-10) -@test sol2_mtk.u[end]≈sol2.u[end] atol=1e-5 +@test sol2_mtk.u[end] ≈ sol2.u[end] From d76d73fe11a81a94439618ea0a0cc9fa8165757c Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 11 Jul 2023 19:18:58 -0400 Subject: [PATCH 7/7] Fix `isdelay` for NonlinearSystem --- src/systems/diffeqs/abstractodesystem.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 4f8bb2bd59..392aa0e1c7 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -159,6 +159,7 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param end function isdelay(var, iv) + iv === nothing && return false isvariable(var) || return false if istree(var) && !ModelingToolkit.isoperator(var, Symbolics.Operator) args = arguments(var)