SymbolicControlSystems.jl

An interface between ControlSystems.jl and SymPy.jl
Author JuliaControl
Popularity
11 Stars
Updated Last
2 Years Ago
Started In
October 2020

SymbolicControlSystems

Utilities for

  • Working with ControlSystems.jl types with SymPy.jl symbols as coefficients.
  • Generation of C-code for filtering with LTI systems.

This package exports the names s,z for the Laplace and Z-transform variables. These can be used to build symbolic transfer functions.

Usage examples

julia> using ControlSystems, SymbolicControlSystems

julia> @vars w T d # Define (SymPy) symbolic variables
(w, T, d)

julia> h = 0.01;

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    222
Ts  + s (2Tdw + 1) + s⎝Tw  + 2dw⎠ + w 

julia> ex = w^2 / (s^2 + 2*d*w*s + w^2) # Define symbolic expression
         2       
        w        
─────────────────
           2    2
2dsw + 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 model

Get 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}\$"

Code generation

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.

A usage example follows

using ControlSystems, SymbolicControlSystems

@vars w T d # Define symbolic variables
h        = 0.01
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(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(Gd_, u); # Filter using Julia
@test norm(y-y_)/norm(y_) < 1e-10
plot([u y y_], lab=["u" "y c-code" "y julia"]) |> display

NOTE: 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.