Simulation

we wish simulate 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)$.

simulate takes in the output $Y(s)$ and, essentially, inverts $Y(s)$ into the time domain to arrive at $y(t)$. in each case, simulate returns an array of times ($t_i$'s) in t and corresponding output values ($y_i$'s where $y_i=y(t_i)$ in y. we can then plot the solution or interpolate these arrays to obtain the value of $y(t)$ at a particular time $t=\tilde{t}$.

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

t, y = simulate(Y, 50.0) # simulate until t = 50

we can then plot the y array versus the t array:

viz_response(t, y, 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

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

viz_response(t, y, 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

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

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

the nb_time_points argument gives a t array with higher resolution for plotting the response.

$y(t)$ at an arbitrary time $\tilde{t}$

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=\tilde{t}$, we can call interpolate to linearly interpolate these data points.

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
t, y = simulate(g / s, 10.0) # unit step response
y_at_τ = interpolate(t, y, τ) # 0.63

under the hood

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

detailed docs

Controlz.simulate — Function

t, y = 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.
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

t::Array{Float64, 1}: array of times: $t_i$'s
y::Array{Float64, 1}: array of $y$ values at corresponding times in t: $y_i$'s, where $y_i=y(t_i)$.

Example

One can simulate the first order step response as, given the Laplace transform of the output, Y:

g = 4 / (3 * s + 1) # first-order transfer function
U = 1 / s #  unit step input
Y = g / s # output
t, y = simulate(Y, 12.0)

source

Controlz.interpolate — Function

y_at_t̃ = interpolate(t, y, t̃)

given an array of times $t_i$ in t and corresponding values of $y(t)$, $y_i=y(t_i)$ in the array y, interpolate to approximate the function $y(t)$ at a new time t̃, i.e. $y(\tilde{t})$.

the output of simulate is an array of times t and corresponding output values y. interpolate is useful for obtaining the solution at a particular time t̃ that is not necessarily present in the array t.

an error is thrown if maximum(t) < t̃ < minimum(t) (extrapolation) or if the t array is not sorted.

Arguments

t::Array{Float64, 1}: array of times $t_i$. these must be sorted.
y::Array{Float64, 1}: an array of corresponding $y$-values $y_i=y(t_i)$
t̃::Float64: the new time at which we wish to know $y$. i.e. we wish to know $y(\tilde{t})$.

Returns

y_at_t̃::Float64: the value of $y$ when $t$ is t̃, $y(\tilde{t})$, according to linear interpolation.

Example

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

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

source