# How to implement a Bessel filter with ModelingToolkit?

How can I implement the following Bessel filter with ModelingToolkit?

\frac{w_1²}{0.618 s²~+~1.3617 w_1 s~+~w_1²}

I know there is this filter component in the standard library: Basic Blocks · ModelingToolkitStandardLibrary.jl

But I don’t think I can use it directly…

There are a lot of blocks in the standard library, but I could not find a generic transfer function block…

UPDATE:
I found this tutorial how to convert a transfer function into a differential equation: Single Diff Eq → Transfer Function

I will try this approach…

1 Like

I’m sure you’re aware that it’s straightforward to convert a transfer function to a linear statespace system. A statespace system is a system of linear differential equations. You should be able to take it from there.

Well, I will manage, but it is no so very user friendly for people who try to convert Simulink models to ModelingToolkit…

Sure you can, just scale both numerator and denominator so that the denominator is monic (this does not change the input-output relationship) and then identify the parameters with those of the SecondOrder type

demonstrates one particular way of obtaining an MTK model for a transfer function.

Thanks!

Using the link in the first post I found the following solution:

Step one: convert the transfer function into a differential equation:
0.618 ~ \ddot y~+~1.3617 ~w_1~ \dot y ~+~ w_1^2~y~=~w_1^2~u

Step two: implement the differential equation with ModelingToolkit:

# Bessel filter with step signal as input
function model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
println("Defining the model...")

@variables t u(t)=0 y(t)=y0 yd(t)=yd0
@parameters w=w1

D = Differential(t)
eqs = [D(y)  ~ yd,
D(yd) ~ 1/0.618 * (w*w * u - 1.3617w * yd - w*w*y),
u     ~ if_else(t > step_time, step_size, 0.0)]

@named sys = ODESystem(eqs, t)
println("Model defined...")

sys = structural_simplify(sys)
println("Model simplified!")
sys, u, y
end


This works fine and has the following step response:

In contrast to a Butterworth filter a Bessel filter has no overshoot. Simulating high-order transfer functions this way is prone to numerical difficulties, the method I the post I linked mitigates this by performing numerical balancing of the statespace system.

Well, second order is not high order, or would you see that differently?

No, second order is not really all that high. You may still experience numerical difficulties due to poor scaling though, even for this simple system. Consider what happens if you increase w1, for w1=1e6, the solver exits with

┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=1.9999999999999998, and step error estimate = 8.930759781475992. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/s9wrq/src/integrator_interface.jl:599
retcode: DtLessThanMin


Already for w1=1000, you can see part of the problem in the solution:

model, u, y = model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
plot(sol)


The derivative state becomes enormous, much larger that the position state.

Compare this to a balanced representation:

using ControlSystemsMTK, ControlSystemsBase

function balanced_model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
println("Defining the model...")

@variables t u(t)=0 y(t)=y0 yd(t)=yd0
@parameters w=w1

H = tf(w1^2, [0.618, 1.3617w1, w1^2])
Hss = ss(H) # This performs balancing
@named filter = ODESystem(Hss)

eqs = [filter.input.u ~ u
u     ~ ifelse(t > step_time, step_size, 0.0)]

@named sys = ODESystem(eqs, t, systems=[filter])
println("Model defined...")

sys = structural_simplify(sys)
println("Model simplified!")
sys, u
end

model, u = balanced_model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
@show length(sol.t)
plot(sol)


Here, the two state variables are of comparable magnitude, making the job for the solver much easier.

To plot the correctly scaled filter output in this case, you’d have to plot plot(sol, idxs=model.filter₊output₊u), since the output is scaled w.r.t. the state.

1 Like

Thank you very much for explaining the balanced solution, but I am hesitant to add more packages to my project unless I really need them. I am using already way too many packages… And wind turbines don’t operate in the MHz range… is that a challenge?

1 Like

I am sorry for the digression from the topic of this thread but reading the code

plot(sol, idxs=model.filter₊output₊u)


I have to share my confusion or perhaps even frustration regarding the separation symbols in MTK. So, dots or sub_pluses? . or ₊? In the above code they are combined and yet the REPL output (as well as the generated plot legend) reads:

julia> model.filter₊output₊u
sys₊filter₊output₊u(t)


But clearly the sub_plus cannot be used everywhere:

julia> model₊filter₊output₊u
ERROR: UndefVarError: model₊filter₊output₊u not defined


If I remember well, you were also involved in the discussion of the choice of the namespace separation char. The conclusion seemed to go back to visualising the separation char as ., but now it does not even work when writing the code.

julia> model.filter.output.u
ERROR: ArgumentError: System sys: variable filter does not exist
Stacktrace:
 getvar(sys::ODESystem, name::Symbol; namespace::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:358
 getvar
@ ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:311 [inlined]
 getproperty(sys::ODESystem, name::Symbol; namespace::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:309
 getproperty(sys::ODESystem, name::Symbol)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/abstractsystem.jl:308
 top-level scope
@ REPL:1


Yes, I don’t really like this either. In this case, the plus was used since the function written by Uwe returned the simplified system rather than the original system. The intended usage is to access names in the original, unsimplified system, in which case you never need the subplus.

Unfortunately, you still occasionally need to call complete on the unsimplified system for this to work properly 