# Step response error for delayed system

**URL:** <https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425>\
**Category:** New to Julia\
**Created:** [June 8, 2022, 11:03am UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425 "2022-06-08T11:03:00Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![lukasdos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lukasdos/32/37010_2.png) [@lukasdos](https://discourse.julialang.org/u/lukasdos)\
**Post date:** [June 8, 2022, 11:03am UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/1 "2022-06-08T11:03:00Z")

</div>

Hi, I am trying to understand to basic functions from ControlSystems package. Especially step() is mystery for me. In some cases this function predicts the simulation time to plot the whole response well, in others it simulates only few seconds (if not specified final time). For stable process, of course.

However, my biggest issue is about simulation of delayed systems. Here I share two similar cases - one works fine, the other one does not.

```julia
using Plots;
using ControlSystems;

WORKS = 1; # switch between 0,1: 1-works fine
            # 0-doesn't work

s = tf("s");

if WORKS == 1
    P = 1 / (0.85*s + 1)*exp(-0.15*s);
elseif WORKS == 0
    P = 1 / (0.85*s + 1)*exp(-0.14*s);
end

display(nyquistplot(P)); #always works

y, t = step(P);
display(plot(t, y'));

```

The only difference is in dead time of the process, but for the numbers closer to 0 an error occures after calling step() function.

```julia
#=
ERROR in the case of WORKS=0

┌ Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=1.4000000000000001. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:...\.julia\packages\SciMLBase\chsnh\src\integrator_interface.jl:366
ERROR: BoundsError: attempt to access 21-element Vector{Tuple{Vector{Float64}, Vector{Float64}}} at index [22]
Stacktrace:
 [1] getindex
   @ .\array.jl:861 [inlined]
 [2] _lsim(sys::DelayLtiSystem{Float64, Float64}, u!::Any, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, x0::Vector{Float64}, alg::DelayDiffEq.MethodOfSteps{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DiffEqBase.NLFunctional{Rational{Int64}, Rational{Int64}}, false}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:abstol, :reltol), Tuple{Float64, Float64}}})
   @ ControlSystems C:...\.julia\packages\ControlSystems\J1JEd\src\delay_systems.jl:204
 [3] lsim(sys::DelayLtiSystem{Float64, Float64}, u::ControlSystems.var"#327#328", t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}; x0::Vector{Float64}, alg::DelayDiffEq.MethodOfSteps{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DiffEqBase.NLFunctional{Rational{Int64}, Rational{Int64}}, false}, abstol::Float64, reltol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ControlSystems C:...\.julia\packages\ControlSystems\J1JEd\src\delay_systems.jl:98
 [4] step(sys::DelayLtiSystem{Float64, Float64}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ControlSystems C:...\.julia\packages\ControlSystems\J1JEd\src\delay_systems.jl:247
 [5] step
   @ C:...\.julia\packages\ControlSystems\J1JEd\src\delay_systems.jl:240 [inlined]
 [6] #step#185
   @ C:...\.julia\packages\ControlSystems\J1JEd\src\timeresp.jl:66 [inlined]
 [7] step(sys::DelayLtiSystem{Float64, Float64})
   @ ControlSystems C:...\.julia\packages\ControlSystems\J1JEd\src\timeresp.jl:66
 [8] top-level scope
   @ c:\...\why_step_doesnt_work.jl:17
=#

```

---

<div class="post-metadata">

**Author:** ![albheim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albheim/32/34660_2.png) [@albheim](https://discourse.julialang.org/u/albheim)\
**Post date:** [June 8, 2022, 11:39am UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/2 "2022-06-08T11:39:58Z")

</div>

Regarding the choice of what time to pick, I think the heuristic is just very simple, and could certainly be improved.

> <https://github.com/JuliaControl/ControlSystems.jl/blob/dae20a11ff53fa2eb4d438f6cebdc01cc9392bdc/src/timeresp.jl#L358-L378>

You can simply add a desired final time to the step function if you want a different simulation time, e.g. `y, t = step(P, 5)`

Regarding the error, it says in the warning that the simulation was aborted. This leads to invalid indexing when it tries to access the data from the failed simulation to generate the output. It can similarly be solved by changing the variable `dtmin` that the warning mentions as `step(P, dtmin=1e-18)`. I’m not sure why this becomes that hard, maybe it is how ControlSystems internally represents and transforms the delay systems?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 8, 2022, 12:13pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/3 "2022-06-08T12:13:44Z")

</div>

Whenever I get problems with the solver hitting `dtmin`, I try `force_dtmin=true`

```julia
res = step(P, 5, force_dtmin=true)
display(plot(res));

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/e/ceee88436f92087d46abee92ad123f0fa2254971.png)

I don’t see why this particular problem would be hard on the solver, but the `force_dtmin` hack appears to have worked.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 8, 2022, 12:18pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/4 "2022-06-08T12:18:27Z")

</div>

Is the discontinuity declared?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 8, 2022, 12:20pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/5 "2022-06-08T12:20:57Z")

</div>

It uses DelayDiffEq.jl with a delay of 0.14s in the troublesome example. The discountinuity is due to the delay, not due to the dynamics

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 8, 2022, 12:22pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/6 "2022-06-08T12:22:37Z")

</div>

Is the discontinuity declared in the DDEProblem definition?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 8, 2022, 12:24pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/7 "2022-06-08T12:24:48Z")

</div>

How do one declare the discontinuity? The definition looks like this

```julia
prob = DDEProblem{true}(dde_param, u0, h!,
                (T(t[1]), T(t[end])),
                p,
                constant_lags=sort(Tau),# Not sure if sort needed
                neutral=true, # We have derivatives on RHS (d(t)) 
                callback=cb)
# Important to stop at t since we can not access derivatives in SavingCallback
sol = DelayDiffEq.solve(prob, alg; tstops=t, saveat=t, kwargs...)

```

where `Tau` contains the delays, 0.14 in this case.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 8, 2022, 12:31pm UTC](https://discourse.julialang.org/t/step-response-error-for-delayed-system/82425/8 "2022-06-08T12:31:26Z")

</div>

That looks right. It would be good to get an MWE of this for DelayDiffEq.jl