Simulation

we wish to simulate the response $y(t)$ (output) of a linear, time-invariant system, characterized by a transfer function $g(s)$, to an input $u(s)$.

pass the output $Y(s)$ in the frequency domain into the function simulate to invert it into the time domain to obtain $y(t)$. simulate returns a time series data frame (a DataFrame, see DataFrames.jl docs), with a :t column for the times and an :output column for $y(t)$. we can then plot the time series data, interpolate the data to obtain the value of $y(t)$ at a particular time $t=\Tau$, etc. we provide examples below.

response of an underdamped second-order system to a unit step input

g = 4 / (4 * s ^ 2 + 0.8 * s + 1) # second order transfer function, underdamped

U = 1 / s # unit step input
Y = g * U # system output

data = simulate(Y, 50.0) # simulate until t = 50, returns DataFrame
data[:, :t]      # array of times, tᵢ's
data[:, :output] # array of outputs, yᵢ's ≈ y(tᵢ)'s

we can then plot the time series via:

viz_response(data, plot_title="SO underdamped step response")

response of a first-order plus time delay system to a unit step input

K = 2.0 # gain
τ = 4.0 # time constant
θ = 1.5 # time delay
g = K * exp(-θ * s) / (τ * s + 1) # FOPTD transfer function

U = 1 / s # step input
Y = g * U

data = simulate(Y, 15.0) # simulate until t = 15

viz_response(data, plot_title="FOPTD step response")

inverse Laplace transform

to emphasize that our simulate function takes a function of the complex frequency s and inverts it into the time domain, consider the Laplace transform of $t \cos(at)$:

\[\mathcal{L}[t \cos(at)] = \dfrac{s^2-a^2}{(s^2+a^2)^2}.\]

we can numerically invert $\dfrac{s^2-a^2}{(s^2+a^2)^2}$ as follows:

a = π
U = (s^2 - a^2) / (s^2 + a^2) ^ 2

data = simulate(U, 8.0, nb_time_points=300) # simulate until t=8, use 300 time points for high resolution

viz_response(data, plot_title=L"inverting an input $U(s)$", plot_ylabel=L"$u(t)$")

the nb_time_points argument allows us to return a time series with a higher resolution in time. if the plot of the response appears jagged, likely you need to increase nb_time_points.

$y(t)$ at an arbitrary time $\Tau$

the simulate function returns an array of times $t_i$'s and corresponding $y_i=y(t_i)$'s. if we wish to know $y(t)$ at a particular time $t=\Tau$, we can call interpolate to linearly interpolate the time series data.

for example, to obtain the output of a first-order system with time constant $\tau$ in response to a unit step input at $t=\tau$:

τ = 3.45
g = 1 / (τ * s + 1) # FO system
data = simulate(g / s, 10.0)  # unit step response
y_at_τ = interpolate(data, τ) # 0.63 ≈ 1 - 1/e, as we expect

under the hood

under the hood, simulate converts the system passed to it into a state space ODE (a system of ODEs) in the time domain and uses DifferentialEquations.jl (see here) to numerically solve the resulting ODE.

detailed docs

Controlz.simulate — Function

data = simulate(Y, final_time, nb_time_points=100) # invert Y(s)

simulate the output $y(t)$ of an LTI system, given the Laplace transform of the output, $Y(s)$, Y.

in other words, simulate inverts an expression in the frequency domain into the time domain.

arguments

Y::TransferFunction: the Laplace transform of the output $y(t)$. usually formed by $g(s)U(s)$, where $U(s)$ is the Laplace transform of the input and $g(s)$ is the transfer function governing the dynamics of the system.
final_time::Tuple{Float64, Float64}: the duration over which to simulate the output of the LTI system, starting at time zero.
nb_time_points::Int=100: the number of time points at which to save the solution $y(t)$

two time points preceding $t=0$ are included to illustrate that it is assumed $y(t)=0$ for $t<0$.

returns

data::DataFrame: data frame containing two columns: :t for time $t$ and :output for $y(t)$. each row corresponds to a $(t_i, y(t_i))$ pair. i.e., row i of the :t column is time $i$, $t_i$, and row i of the :output column is $y_i=y(t_i)$. access the columns by data[:, :t] and data[:, :output].

example

simulate the first order step response to a step, given the Laplace transform of the output, Y:

g = 4 / (3 * s + 1)      # first-order transfer function g(s)
U = 1 / s                # unit step input U(s)
Y = g / s                # output Y(s)
data = simulate(Y, 12.0) # time series data frame
data[:, :t]              # array of time points tᵢ
data[:, :output]         # array of corresponding outputs y(tᵢ)

source

Controlz.interpolate — Function

y_at_Τ = interpolate(data, Τ)

interpolate a data frame containing a time series characterizing $y(t)$, with :t and :output columns, whose rows are $(t_i, y(t_i))$ pairs. interpolate the data to approximate the function $y(t)$ at a new time Τ, i.e. $y(\Tau)$.

simulate returns such a data frame, thus interpolate is useful for obtaining the solution at a particular time Τ that is not necessarily present in the :t column of data.

arguments

data::DataFrame: a data frame of a time series, containing a :t column for times and :output column for outputs.
Τ::Float64: the new time at which we wish to know $y$. i.e. we wish to know $y(\Tau)$.

Returns

y_at_Τ::Float64: the value of $y$ when $t$ is equal to Τ, $y(\Tau)$, according to linear interpolation.

example

the unit step response of a first-order process with time constant $\tau$ is $\approx 63\%$ of the final value when $t=\tau$.

τ = 3.45
g = 1 / (τ * s + 1)           # FO system
data = simulate(g / s, 10.0)  # unit step response
y_at_τ = interpolate(data, τ) # 0.63

source

Controlz.viz_response — Function

viz_response(data, 
             plot_title="", plot_xlabel="time, t", 
             plot_ylabel="output, y(t)",
             savename=nothing)

plot data[:, :output] vs. data[:, :t] to visualize the response of a system to an input. typically the data frame, data, is returned from simulate.

note that PyPlot.jl (matplotlib) commands can be invoked after viz_response to make further changes to the figure panel. e.g. xlim([0, 1]) can be applied after viz_response.

Arguments

data::DataFrame: data frame of time series data, containing a :t column for times and :output column for the outputs.
plot_title::Union{String, LaTeXString}: title of plot
plot_xlabel::Union{String, LaTeXString}: x-label
plot_ylabel::Union{String, LaTeXString}: y-label
savename::Union{Nothing, String}: filename to save as a figure in .png format (dpi 250).

Example

g = 4 / (4 * s ^ 2 + 0.8 * s + 1)
u = 1 / s
data = simulate(g * u, (0.0, 50.0))
viz_response(data)

source