Miss-compilation when Hold is applied to delayed discrete variable

[baggepinnen]

The following models a discrete-time delay $y(k) = u(k-d)$. At $t = 2$, a step in the input occurs from 0 to 1, it should take $d$ seconds for the step to appear in the output, instead, it takes $d-1$ seconds, and the behavior for $d=0$ and $d=1$ are identical.

The problem is exposed by the following analysis

1. We use d discrete-time state variables to store the delayed signal. Let's assume d = 1 for simplicity
2. At time t, the input to the discrete partition changes (in the example t = 2 and the change is from 0 to 1)
3. At time t, the c2d variable is thus changed by c2d_obs(integrator.u, p, t)
4. At time t, the disc update stores the changed signal value in the discrete-time state variable which is used as memory
5. At time t, the d2c variable is set to the value of the discrete-time state variable. This is the incorrect step, the discrete update correctly represents x(k) = f(x(k-1), c2d(k), but the d2c update in this case has to output the old discrete-time state, not the new, since the Hold operator was applied to the delayed input signal input(k-d).

There are two ways to handle this issue

1. The compiler can realize that an additional state variable is required in order to keep the old value in memory until the d2c pass is run. Right now, the compiler correctly handles cases when the d2c involves state variables at time k, but not at k-d for some positive d.
2. Split the d2c update in two, one that is run before the discrete affect, and one after. This would reduce the number of state variables, but adds complexity instead.

MWE

using ModelingToolkit, Plots, OrdinaryDiffEq
using ModelingToolkit: D_nounits as D
using ModelingToolkit: t_nounits as t
k = ShiftIndex(Clock(t, 1))
@mtkmodel DelayModel begin
    @variables begin
        input(t) = 0
        delay(t) = 0
        x(t) = 0
    end
    @structural_parameters begin
        d
    end
    @equations begin
        input ~ (t >= 2)
        delay(k) ~ input(k-d)
        D(x) ~ (-x + Hold(delay))/1e-3 # Simple fast lowpass filter to be able to show a continuous-time output variable
    end
end

plot()
for d = 0:3
    @mtkbuild m = DelayModel(; d)
    prob = ODEProblem(m, [m.delay(k-3)=>0, m.delay(k-2)=>0, m.delay(k-1)=>0], (0.0, 10.0))
    sol = solve(prob, Tsit5(), kwargshandle = KeywordArgSilent)
    plot!(sol, lab="d = $d")
end
display(current())

the plot shows the result for a few different values of $d$

[YingboMa]

# Codegen issue
# Periodic callback of dt = 1
# args: [input(k-1)] # currect behavior
input(k) = (t >= 2)
delay(k) = input(k-1)
-> # Current wrong implementation does
# args: [input]
input = (t >= 2)
delay = input

# Codegen result
function discrete_update(out, u, p, t)
    out[1] = (>=)(t, 2) # input(t)
end
function d2c(disc_unknowns, ...)
    var"Shift(t, -1)(input(t))" = disc_unknowns[1]
    var"delay(t)" = var"Shift(t, -1)(input(t))"
end
...
discrete_update(disc_unknowns, u, p, t)
d2c(disc_unknowns, ...)
...
# note that both `input(k-1)` and `input(k)` are aliased at `disc_unknowns[1]`.

# Conclusion: unalias memory.