Incorrect solution for DiscreteSystem

[ValentinKaisermayer]

Consider the following model:

using ModelingToolkit, DifferentialEquations, Plots

@variables t
function ZOH(dt; name)
    @variables u(t)=0.0 y(t)=0.0
    U = DiscreteUpdate(t; dt=dt)
    eqs = [
        U(y) ~ u
    ]
    DiscreteSystem(eqs, t, name=name)
end
function Sin(;name, omega=1, amplitude=1)
    @variables y(t)=0.0
    @parameters omega=omega amplitude=amplitude
    eqs = [
        y ~ amplitude * sin(omega * t)
    ]
    DiscreteSystem(eqs, t, name=name)
end

dt = 1
@named zoh_block = ZOH(dt)
@named sin_block = Sin()
@named model = DiscreteSystem([zoh_block.u ~ sin_block.y], systems=[zoh_block, sin_block])
prob = DiscreteProblem(model, [], (0, 10.0), [])
sol = solve(prob, FunctionMap(), adaptive=false; dt)

plot()
plot!(sol)
plot!(0:0.01:10, sin, label="sin(t)")

The solution is incorrect:

julia> sol
retcode: Success
Interpolation: left-endpoint piecewise constant
t: 11-element Vector{Float64}:
  0.0
  1.0
  2.0
  3.0
  ⋮
  7.0
  8.0
  9.0
 10.0
u: 11-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.8414709848078965]
 [0.0, 0.0, 0.9092974268256817]
 [0.0, 0.0, 0.1411200080598672]
 ⋮
 [0.0, 0.0, 0.6569865987187891]
 [0.0, 0.0, 0.9893582466233818]
 [0.0, 0.0, 0.4121184852417566]
 [0.0, 0.0, -0.5440211108893698]

Although zoh_block₊y(t) is correct, the other two (zoh_block₊u(t) and sin_block₊y(t)) are wrong (constant to zero).

 zoh_block₊u(t) ~ sin_block₊y(t)
 Difference(t; dt=1, update=true)(zoh_block₊y(t)) ~ zoh_block₊u(t)
 sin_block₊y(t) ~ sin_block₊amplitude*sin(sin_block₊omega*t)

[baggepinnen]

A few comments. I believe you intended to implement a zero-order hold as opposed to first-order hold? In any case, I think the implementation U(y) ~ u would model a unit delay? I.e.,

y(t+dt) = u(t)

what you'd really like is to communicate (pseudo code)

y(t to (t+dt)) = u(t)

but there is currently no way to do this. The semantics for hybrid systems (both discrete and continuous) is still on the discussion stage https://github.com/SciML/ModelingToolkit.jl/issues/894

The way Modelica does this is with the when block in combination with the previous/pre operator (similar to DiscreteUpdate) like so

    when {sampleTrigger, initial()} then
      ySample = u;
    end when;
    y = pre(ySample);

[baggepinnen]

As for the rest of the issue (states are 0). There are currently a lot of work to do on discrete systems

• https://github.com/SciML/ModelingToolkit.jl/issues/1330

indicates that algebraic equations are currently not handled and simplification does not work on discrete systems.

[ValentinKaisermayer]

A few comments. I believe you intended to implement a zero-order hold as opposed to first-order hold?

You are right, I wanted a zero-orderhold. I fixed that in the OP.

So DiscreteUpdate does nothing and the behavior is simply because of the time discrete evaluation of sin(t)?

[ValentinKaisermayer]

However, I think for ODESystems, DiscreteUpdate is really a ZOH (updates variables at certain discrete timestamps): https://github.com/SciML/ModelingToolkit.jl/blob/9cc6f21c1ffb1299d1a2d120c848f036d5f232c9/src/systems/diffeqs/abstractodesystem.jl#L129

For the DiscreteSystem, DiscreteUpdate seems to be ignored completely.

[ValentinKaisermayer]

It seems that the system:

julia> equations(model)
3-element Vector{Equation}:
 foh_block₊u(t) ~ sin_block₊y(t)
 Difference(t; dt=1, update=true)(foh_block₊y(t)) ~ foh_block₊u(t)
 sin_block₊y(t) ~ sin_block₊amplitude*sin(sin_block₊omega*t)  

julia> states(model)
3-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 foh_block₊u(t)
 sin_block₊y(t)
 foh_block₊y(t)

is incorrectly converted into a DiscreteProblem.

Using the method: https://github.com/SciML/ModelingToolkit.jl/blob/9cc6f21c1ffb1299d1a2d120c848f036d5f232c9/src/systems/discrete_system/discrete_system.jl#L253 the following function is generated:

:(function (ˍ₋out, ˍ₋arg1, ˍ₋arg2, t)
      let (var"foh_block₊u(t)" =  @inbounds(ˍ₋arg1[1]), 
      var"sin_block₊y(t)" =  @inbounds(ˍ₋arg1[2]), 
      var"foh_block₊y(t)" =  @inbounds(ˍ₋arg1[3]), 
      sin_block₊omega = @inbounds(ˍ₋arg2[1]), 
      sin_block₊amplitude =  @inbounds(ˍ₋arg2[2]))

         @inbounds begin
                  ˍ₋out[1] = var"sin_block₊y(t)"
                  ˍ₋out[2] = var"foh_block₊u(t)"
                  ˍ₋out[3] = (*)(sin_block₊amplitude, (sin)((*)(sin_block₊omega, t)))
                  nothing
              end
      end
  end)

Which seems to be not in the correct order (if the states are ordered correctly in states(model). That explains why foh_block₊y(t) has a non-zero value in the solution but the others not.

[baggepinnen]

Interesting, I'll take a look at this issue today.

[baggepinnen]

I think the problem is once again that there is no simplification step running when discrete problems are created, similar to how it's run for ODEProblem. This leads to the equations being in the wrong order. https://github.com/SciML/ModelingToolkit.jl/issues/1330#issuecomment-978903559 summarizes what needs to be done.

I'll dig a bit deeper into how discrete and continuous equations coexist when an ODEProblem is created, which I believe is a separate issue.

[baggepinnen]

If the system is treated as an ODE system with discrete states, alias_alimination rewrites the difference equation such that the difference term ends up on the rhs. That leads to isdifferenceeq not identifying the equation as a difference eq. and treats it as an algebraic equation. The difference callbacks are never generated.

Here's an exmaple

using ModelingToolkit, OrdinaryDiffEq

@variables t
function ZOH(dt; name)
    @variables u(t)=0.0 y(t)=0.0
    z = Shift(t; dt=dt)
    eqs = [
        z(y) ~ u
    ]
    ODESystem(eqs, t, name=name)
end
function Sin(;name, omega=1, amplitude=1)
    @variables y(t)=0.0
    @parameters omega=omega amplitude=amplitude
    eqs = [
        y ~ amplitude * sin(omega * t)
    ]
    ODESystem(eqs, t, name=name)
end
function Plant(;name)
    @variables y(t)=0.0 x(t)=0.0 u(t)=0.0
    eqs = [
        Differential(t)(x) ~ -x + u
        y ~ x
    ]
    ODESystem(eqs, t, name=name)
end

dt = 1
@named zoh_block = ZOH(dt)
@named sin_block = Sin()
@named plant = Plant()
@named model = ODESystem([zoh_block.u ~ sin_block.y, plant.u ~ zoh_block.y], systems=[zoh_block, sin_block, plant])

model = alias_elimination(model) # wrong after this step

# model = structural_simplify(model)
prob = ODEProblem(model, [], (0, 10.0), [])
sol = solve(prob, Rosenbrock23())

[ValentinKaisermayer]

Can you print the equations? Before and after

[baggepinnen]

julia> equations(model)
6-element Vector{Equation}:
 zoh_block₊u(t) ~ sin_block₊y(t)
 plant₊u(t) ~ zoh_block₊y(t)
 Difference(t; dt=1, shift=true)(zoh_block₊y(t)) ~ zoh_block₊u(t)
 sin_block₊y(t) ~ sin_block₊amplitude*sin(sin_block₊omega*t)
 Differential(t)(plant₊x(t)) ~ plant₊u(t) - plant₊x(t)
 plant₊y(t) ~ plant₊x(t)

julia> equations(alias_elimination(model))
3-element Vector{Equation}:
 0 ~ sin_block₊y(t) - Difference(t; dt=1, shift=true)(zoh_block₊y(t))
 0 ~ sin_block₊amplitude*sin(sin_block₊omega*t) - sin_block₊y(t)
 Differential(t)(plant₊x(t)) ~ zoh_block₊y(t) - plant₊x(t)

[ChrisRackauckas]

Old discrete interface was removed.