Solution interpolation gives NaN while no NaN in solution

[baggepinnen]

Here's an example that builds a hybrid cont./disc. system and successfully produces a solution. The sol.u contains only finite values, but sol(0) is NaN

dt = 0.5 # Sample interval
@variables r(t)
clock = Clock(t, dt)
k = ShiftIndex(clock)

function plant(; name)
    @variables x(t)=1 u(t) y(t)
    D = Differential(t)
    eqs = [D(x) ~ -x + u
           y ~ x]
    ODESystem(eqs, t; name = name)
end

function filt(; name) # Reference filter
    @variables x(t) u(t) y(t)
    a = 1 / exp(dt)
    # eqs = [x(k) ~ a * x(k - 1) + (1 - a) * u(k-1)
    #        y ~ x(k-1)]
    eqs = [x(k) ~ a * x(k - 1) + (1 - a) * u(k - 1)
           y ~ x(k)]
    ODESystem(eqs, t, name = name)
end

function controller(kp; name)
    @variables y(t) r(t) ud(t) yd(t)
    @parameters kp = kp
    eqs = [yd ~ Sample(clock)(y)
           ud ~ kp * (r - yd)]
    ODESystem(eqs, t; name = name)
end

Tf = 15
kp = 1
@named f = filt()
@named c = controller(kp)
@named p = plant()

connections = [r ~ (t >= 1) + (t >= 6)          # reference signal
               f.u ~ r             # reference to filter input
               f.y ~ c.r           # filtered reference to controller reference
               Hold(c.ud) ~ p.u    # controller output to plant input (zero-order-hold)
               p.y ~ c.y]          # plant output to controller feedback

@named cl = ODESystem(connections, t, systems = [f, c, p])
model = complete(cl)
ssys = structural_simplify(model)
prob = ODEProblem(ssys, [], (0, Tf),
    [model.f.x(k - 1) => 1, model.f.u(k - 1) => 1]);
sol_mtk = solve(
    prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8, kwargshandle = KeywordArgSilent);

julia> sol_mtk(0)
1-element Vector{Float64}:
 NaN

julia> sol_mtk.u
146-element Vector{Vector{Float64}}:
 [1.0]
 [1.0]
 [0.9851011821935]
 [0.9429484364997274]

[ChrisRackauckas]

Is this MTK related at all? Can this not be simplified down to just OrdinaryDiffEq.jl?