Symbolic support for DDE histories

Author: YingboMa

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...)

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