Step response using ModelingToolkit

Specific Domains Modelling & Simulations
,
1

I have a non-linear system defined using ModelingToolkit. Now I want to plot the step response. I use the following code:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations, Plots

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

# helper for step function
@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
    return ifelse(condition, trueval, falseval)
end
ModelingToolkit.@register if_else(x, trueval, falseval)
@eval ModelingToolkit using ..Main: if_else

# initial state
U0  = 9.0         # m/s
Pg_0 = 2.302e6    # W
step_time = 200.0 # s
step_size = 23000 #  ΔPg in W, about 1% of initial value

@variables t ω(t) Pr(t) λ(t) Pg(t)
@parameters J R ρ A U Pg0
D = Differential(t)

eqs = [J * D(ω) * ω ~ Pr - Pg,
       λ            ~ ω * R/U,
       Pg           ~ Pg0 + if_else(t > step_time, step_size, 0.0),
       Pr           ~ 0.5 * ρ * A * U^3 * cp(λ)]

@named sys = ODESystem(eqs, t)
sys = structural_simplify(sys)

u0  = [ω    => 1.302, # initial value
       D(ω) => 0.0]   # initial estimate

# system parameters
p = [ J => 4.0470e+07,
      R => 63,
      ρ => 1.225,
      A => 1.2469e+04,
      U => U0,
      Pg0 => Pg_0]

# simulation parameters
tspan = (0.0, 1000.0)
ts    = 0:0.01:1000.0

# run the simulation
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rosenbrock23(), dt=0.01, saveat=ts)

# plot the result
X = sol.t
OM = sol(X, idxs=ω)
plot(X, OM)

I get the following plot as result:
Figure_1

Which looks nice, but it is very wrong. The step happens at 200s, but already at 184 s the output changes. Expected result: A transition from the upper level to the lower level similar to a first order system.

What am I doing wrong?

Remark: Without the step function the output is correct.

2

There seem to be a step generator in mtk standard library, check first example here

Not an answer to why your way didn’t work, but maybe checking if this works and how it’s implemented could help.

1 Like
3

Thanks! I looked at your example, and found they use the Rodas5() solver. After switching the solver I already get a much better result:
Figure_1

But it still has some distortion before t=200s, which is not present if I use Simulink to simulate this example:
Figure_1

Any idea why, and how to avoid this distortion?

UPDATE: Rodas4() gives an even better result, but there is still an overshoot before t=200s.

4

These solver settings give a very accurate result without overshoot:

sol = solve(prob, Rodas4(), dt=0.001, reltol = 1e-9, saveat=ts)

But how can I be sure the choice of the solver and the settings are good?

5

Even better:

solve(prob, Rodas4(), dt=0.001, abstol=1e-9, reltol = 1e-9, saveat=ts)

Figure_1

6

More or less the same issue was reported here: MTK -- weird solution with step inputs · Issue #2163 · SciML/ModelingToolkit.jl · GitHub

7

The standard Hermite interpolation used for ODE methods in OrdinaryDiffEq.jl is not applicable to the algebraic variables. Thus, for the following mass-matrix methods, use the interpolation (thus saveat ) with caution if the default Hermite interpolation is used. All methods which mention a specialized interpolation (and implicit ODE methods) are safe.

I would double check the interpolation of Rodas4 with algebraic variables. Rodas4’s may not satisfy the algebraic variable conditions and we should mark that in the docs. Rodas4P and Rodas5P should have one that does though IIRC?

1 Like
8

I do not like this kind of “info” in the docs, because it raises more questions than it answers: How shall I know if a solver uses a specialized interpolation or not? And why should I need interpolation in the first place in I specify dt and save the results only at multiples of dt?

9

You only specified the initial dt and used saveat instead of tstops so it’s using the interpolation. This is something that’s mentioned in most of the tutorials and follows the same convention seen in most libraries.

It’s all documented, and it has an overload to the show method that says whether or not.

10

Thanks for the quick reply.

Nevertheless I find the documentation of Assimulo so much clearer and easier to understand than the documentation of DifferentialEquations: RodasODE — Assimulo 3.0 documentation

For example it has a separate page for each solver which clearly states which features it supports and which not,

11

Those pages just aren’t made yet.

https://docs.sciml.ai/OrdinaryDiffEq/stable/

You can help finish it.