Step response error for delayed system

New to Julia
1

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.

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.

#=
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
=#
2

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

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?

3

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

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

image
image779×577 15.2 KB

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

1 Like
4

Is the discontinuity declared?

5

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

6

Is the discontinuity declared in the DDEProblem definition?

7

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

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.

8

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