Utilities for
- Working with ControlSystems.jl types with SymPy.jl and Symbolics.jl symbols as coefficients.
- Generation of C-code for filtering with LTI systems.
This package exports the names s,z of type SymPy.Sym for the Laplace and Z-transform variables. These can be used to build symbolic transfer functions. To build symbolic transfer functions with Symbolics.jl symbols, create s or z using using Symbolics; @variables s.
Users typically want to install both ControlSystemsBase and SymbolicControlSystems. ControlSystemsBase contains the basic control-systems functionality, like system types etc., that is used when working with SymbolicControlSystems.
using Pkg
Pkg.add(["ControlSystemsBase", "SymbolicControlSystems"])julia> using ControlSystemsBase, SymbolicControlSystems
julia> @vars w T d # Define (SymPy) symbolic variables
(w, T, d)
julia> h = 0.01; # Sample time
julia> G = tf([w^2], [1, 2*d*w, w^2]) * tf(1, [T, 1])
TransferFunction{Continuous, SisoRational{Sym}}
w^2
-----------------------------------------------
T*s^3 + 2*T*d*w + 1*s^2 + T*w^2 + 2*d*w*s + w^2
Continuous-time transfer function model
julia> Gd = tustin(G, h); # Discretize
julia> Sym(G) # Convert a TransferFunction to symbolic expression
2
w
───────────────────────────────────────────────
3 2 ⎛ 2 ⎞ 2
T⋅s + s ⋅(2⋅T⋅d⋅w + 1) + s⋅⎝T⋅w + 2⋅d⋅w⎠ + w
julia> ex = w^2 / (s^2 + 2*d*w*s + w^2) # Define symbolic expression
2
w
─────────────────
2 2
2⋅d⋅s⋅w + s + w
julia> tf(ex) # Convert symbolic expression to TransferFunction
TransferFunction{Continuous, SisoRational{Sym}}
w^2
---------------------
1*s^2 + 2*d*w*s + w^2
Continuous-time transfer function model
julia> # Replace symbols with numbers
T_, d_, w_ = 0.03, 0.2, 2.0 # Define system parameters
(0.03, 0.2, 2.0)
julia> Gd_ = sym2num(Gd, h, Pair.((T, d, w), (T_, d_, w_))...)
TransferFunction{Discrete{Float64}, SisoRational{Float64}}
1.4227382019434605e-5*z^3 + 4.2682146058303814e-5*z^2 + 4.2682146058303814e-5*z + 1.4227382019434605e-5
-------------------------------------------------------------------------------------------------------
1.0*z^3 - 2.705920013658287*z^2 + 2.414628594192383*z - 0.7085947614779404
Sample Time: 0.01 (seconds)
Discrete-time transfer function modelGet a Latex-string
julia> latextf(G)
"\$\\dfrac{1.0w^2}{0.003s^3 + s^2(0.006dw + 1.0) + s(2.0dw + 0.003w^2) + 1.0w^2}\$"The function code = SymbolicControlSystems.ccode(G::LTISystem) returns a string with C-code for filtering of a signal through the linear system G. All symbolic variables present in G will be expected as inputs to the generated function. The transfer-function state is handled by the C concept of static variables, i.e., a variable that remembers it's value since the last function invocation. The signature of the generated function transfer_function expects all input arguments in alphabetical order, except for the input u which always comes first.
Code generation for systems with multiple inputs and outputs (MIMO) is only handled for statespace systems, call ss(G) to convert a transfer function to a statespace system.
A usage example follows
using ControlSystemsBase, SymbolicControlSystems
@vars w T d # Define symbolic variables
h = 0.01 # Sample time
G = tf([w^2], [1, 2*d*w, w^2]) * tf(1, [T, 1])
Gd = tustin(G, h) # Discretize
code = SymbolicControlSystems.ccode(Gd, cse=true)
path = mktempdir()
filename = joinpath(path, "code.c")
outname = joinpath(path, "test.so")
write(joinpath(path, filename), code)
run(`gcc $filename -lm -shared -o $outname`)
## Test that the C-code generates the same output as lsim in Julia
function c_lsim(u, T, d, w)
Libc.Libdl.dlopen(outname) do lib
fn = Libc.Libdl.dlsym(lib, :transfer_function)
map(u) do u
@ccall $(fn)(u::Float64, T::Float64, d::Float64, w::Float64)::Float64
end
end
end
u = randn(1,100); # Random input signal
T_, d_, w_ = 0.03, 0.2, 2.0 # Define system parameters
y = c_lsim( u, T_, d_, w_); # Filter u through the C-function filter
Gd_ = sym2num(Gd, h, Pair.((T, d, w), (T_, d_, w_))...) # Replace symbols with numeric constants
y_,_ = lsim(ss(Gd_), u); # Filter using Julia
using Plots, LinearAlgebra, Test
@test norm(y-y_)/norm(y_) < 1e-10
plot(u', lab="u", layout=2)
plot!([y; y_]', lab=["y c-code" "y julia"], sp=2, linestyle=[:solid :dash]) |> displayNOTE: Numerical accuracy
The usual caveats for transfer-function filtering applies. High-order transfer functions might cause numerical problems. Consider either filtering through many smaller transfer function in series, or convert the system into a well-balanced statespace system and generate code for this instead. See lecture notes slide 45 and onwards as well as the ControlSystems docs on numerical accuracy.. The function
ControlSystems.ssconverts a transfer function to a statespace system and performs automatic balancing.
The following example writes C-code that interpolates between two linear systems.
The interpolation vector t defines the interpolation points.
The system in the example is a double-mass-spring damper, where the inertia of the load is allowed to vary:
function double_mass_model(Jl) # Inertia load
Jm = 1 # Inertia motor
k = 100 # Spring constant
c0 = 1 # Dampings
c1 = 1
c2 = 1
A = [
0.0 1 0 0
-k/Jm -(c1 + c0)/Jm k/Jm c1/Jm
0 0 0 1
k/Jl c1/Jl -k/Jl -(c1 + c2)/Jl
]
B = [0, 1/Jm, 0, 0]
C = [1 0 0 0]
ss(A,B,C,0)
end
t = [1, 5] # The different inertias in the interpolation
sys = [c2d(double_mass_model(inertia), 0.01) for inertia in t]
SymbolicControlSystems.print_c_array(stdout, sys, t, "mass_spring_damper")This will print C-code functions for the interpolation of each of the system matrices. See the docstring for print_c_array for more customization options.