Converting a symbolic transfer function to a TransferFunction

Specific Domains Modelling & Simulations
, ,
1

Hello,
any idea how to convert the following transfer function:

julia> res
(-9.770075179266223e-8(s^4) - 2.543931256784638e-6(s^5) - 3.198336516372632e-5(s^6) - 0.00019354918819587642(s^7) - 0.000453357561824716(s^8)) / ((0.00013932867142612727 + 0.0124908808235244s + s^3 + 0.2134619608161619(s^2))^3)
julia> typeof(res)
Num

into a transfer function with the type TransferFunction ?

2

I mean, I can convert it manually:

using ControlSystemsBase
num = [0.0007205235056679737, 0.00029812769938757924, 4.8098814990661017e-5, 3.763125715852026e-6, 1.4315223278917721e-7,0,0,0,0]
den = [0.206882701981523, 0, 0, 1, 0, 0.011977460042917651, 0.0001334537369270392]
tf(num, den)

But it should also be possible automatically somehow, I thinkā€¦

3

One (untestested) idea is to create a function

f=build_function(res, [s])

f = eval(f[1])

And then evaluate with the transfer function ā€œsā€

S=tf([1.0, 0.0],[1.0])

H=f(S)

1 Like
4

SymbolicControlSystems.jl contains a number of functions for converting to and from sympy, but the support for Symbolics.jl is not that well developed. However, it does look like the following works

julia> using Symbolics

julia> @variables s;

julia> (s+1) / (s+2);

julia> using SymbolicControlSystems

julia> G = (s+1) / (s+2)
(1 + s) / (2 + s)

julia> tf(G)
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s + 1.0
----------
1.0s + 2.0

Continuous-time transfer function model
2 Likes
5

This works:

using Symbolics, ControlSystemsBase

@variables s

import SymbolicControlSystems

G = (s+1) / (s+2)
tf(G)
6

Isnā€™t this exactly what the previous post did?

7

Well, you forgot to import ControlSystemBase, and you wrote using SymbolicControlSystems, not import SymbolicControlSystemsā€¦

Using results in a warning because s is alreay defined.

But I have another problem: After executing import SymbolicControlSystems my linear algebra code fails.

This line of my code then fails:

 G1 = (C*inv(M)*B + D)

with

ERROR: MethodError: no method matching oneunit(::Type{Any})

Stacktrace:
 [1] oneunit(#unused#::Type{Any})
   @ Base .\missing.jl:106
 [2] inv(A::Matrix{Any})
   @ LinearAlgebra C:\Users\uwefechner\.julia\juliaup\julia-1.9.3+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\dense.jl:910

Is SymbolicControlSystems doing some type piracy?

I thought that import is safe and should not impact already defined functionsā€¦

8

I donā€™t think so, in particular, we do not define any methods for neither inv nor oneunit. Maybe SymPy or Latexify does?

9

I will checkā€¦

10

I use an external script now:

# script to convert symbolic expression into a TransferFunction
using Symbolics, ControlSystemsBase, JLD2
import SymbolicControlSystems

expr = load("data/H_num.jld2")["H_num"]
H_tf = tf(expr)
jldsave("data/H_tf.jld2"; H_tf)

and execute it with:

jldsave("data/H_num.jld2"; H_num)
# execute conversion script ...
mycommand = `julia --project -e "include(\"src/convert_to_tf.jl\")"`
run(mycommand)
H_tf = load("data/H_tf.jld2")["H_tf"]

I am open for suggestions how to do this in a more elegant way, but at least this approach worksā€¦

11

What is M, it looks like itā€™s a Matrix of Any in your stack trace. What is M if you do not load SymbolicControlSystems?

12

A matrix of symbolic expressions:

julia> M
2Ɨ2 Matrix{Num}:
 s + (-Pge(t) - Pgc(t) + 0.5A*(U^3)*b*Ī“*Ļ) / (J*(Ļ‰(t)^2))           0
      -(0.5A*R*(U^2)*a*Ī“*Ļ - 0.5A*R*a*(Uest(t)^2)*Ī“est*Ļ)*Ku  s - (-1.5A*(b*Uest(t) + R*a*Ļ‰(t))*Uest(t)*Ī“est*Ļ + 0.5A*R*a*Ļ‰(t)*Uest(t)*Ī“est*Ļ)*Ku

And the same if I load SymbolicControlSystemsā€¦ Only matrix operations failā€¦

Update:
At some point in time the type of the matrix changes to 2Ɨ2 Matrix{Any}, even though each element
is still of type Numā€¦

13

Much more elegant and much faster:

using Symbolics, ControlSystemsBase

@variables R a b c d Ī“ Ļ U J Q s

input = 2s + 4U*a*s^3 + Q*(J*s+s)/(4Ī“*U*s + 2a)

function split_s(input)
    inp = simplify(input, expand=true)
    vcat([Num(Symbolics.coeff(inp, s^i)) for i in Symbolics.degree(inp, s):-1:1],
    substitute(inp, s => 0))
end

function symbolic2tf(expr::Num)
    numerator, denominator = Symbolics.arguments(Symbolics.value(simplify(expr)))
    num = split_s(numerator)
    den = split_s(denominator)
    num, den
    tf(num, den)
end

res = symbolic2tf(input)

@baggepinnen What do you think about this solution?

14

The use of Symbolics.arguments is only correct if / is the top-level operator, you should probably check for that

15

Thatā€™s why I simplify the expression firstā€¦

16

Thatā€™s not sufficient, your function produces the wrong result for this PID controller

@variables s
1 + s + 1/s

julia> symbolic2tf(1 + s + 1/s)
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Num}}
1
--
1s

Continuous-time transfer function model
17

Thanks for testing!