a closed-loop transfer function that relates an output Y and an input U in a feedback loop.

the resulting closed-loop transfer function is:

 Y      top
--- = --------
 U    1 + g_ol

example

g_ol = 4 / (s + 1) * 2 / (s + 2)
top = 5 / (s + 4)
g = ClosedLoopTransferFunction(top, g_ol)
# output
closed-loop transfer function.
      top
    -------
    1 + g_ol

  top =
    5.0
-----------
1.0*s + 4.0

  g_ol =
         8.0
---------------------
1.0*s^2 + 3.0*s + 2.0

attributes

• top::TransferFunction: numerator
• g_ol::TransferFunction: open-loop transfer function

pc = PController(Kc)

Construct a Proportional (P) controller by specifying the controller gain defined under the following transfer function representation:

\[g_c(s)=K_c\]

Arguments

• Kc::Float64: controller gain

Example

pc = PController(1.0) # specify P controller gain
gc = TransferFunction(pc) # construct transfer function with this P-controller gain

pic = PIController(Kc, τI)

Construct a Proportional-Integral (PI) controller by specifying the controller gain and integral time constant defined under the following transfer function representation:

\[g_c(s)=K_c \left[1+\frac{1}{\tau_I s}\right]\]

Arguments

• Kc::Float64: controller gain
• τI::Float64: integral time constant

Example

pic = PIController(1.0, 3.0) # specify PI controller params
gc = TransferFunction(pic) # construct transfer function with these PI-controller params

pidc = PIDController(Kc, τI, τD, α=0.0)

Construct a Proportional-Integral-Derivative (PID) controller by specifying the controller gain, integral time constant, derivative time constant, and derivative filter defined under the following transfer function representation:

\[g_c(s)=K_c \left[1+\frac{1}{\tau_I s}+\tau_D s \frac{1}{\alpha \tau_D s + 1}\right]\]

Arguments

• Kc::Float64: controller gain
• τI::Float64: integral time constant
• τD::Float64: derivative time constant
• α::Float64: derivative filter

Example

pidc = PIDController(1.0, 3.0, 0.1) # specify PID controller params
gc = TransferFunction(pidc) # construct transfer function with these PID-controller params

tf = TransferFunction([1, 2], [3, 5, 8])
tf = TransferFunction([1, 2], [3, 5, 8], 3.0)

construct a transfer function representing a linear, time-invariant system.

example

to construct the transfer function

\[G(s) = \frac{4e^{-2.2s}}{2s+1}\]

in Julia:

tf = TransferFunction([4], [2, 1], 2.2)
# output
    4.0
----------- e^(-2.2*s)
2.0*s + 1.0

attributes

• numerator::Polynomial{Float64, :s}: the polynomial in the numerator of the transfer function
• denominator::Polynomial{Float64, :s}: the polynomial in the denominator of the transfer function
• time_delay::Float64: the associated time delay

tf_time_delay = exp(- 3 * s)

Conveniently create a time delay by exp(- θ * s).

Example

g = 1 / (s + 1) * exp(-2.0 * s) # introduce time delay of 2.0

axs = bode_plot(tf, log10_ω_min=-4.0, log10_ω_max=4.0, nb_pts=300)

draw the Bode plot of a transfer function tf to visualize its frequency response. returns the two axes of the plot for further tuning via matplotlib commands.

adjust the range of frequencies that the Bode plot presents with log10_ω_min and log10_ω_max.

increase the resolution of the Bode plot with nb_pts.

returns a CairoMakie.jl Figure object for further modification.

p = characteristic_polynomial(g_ol)

Determine the characteristic polynomial associated with open loop transfer function g_ol.

The characteristic polynomial is $1+g_{ol}(s)$. The roots of the characteristic polynomial determine the character of the response of the closed loop system to bounded inputs.

Arguments

• g_ol::TransferFunction: open loop transfer function

Returns

a polynomial of type Polynomial

Example

g_ol = 4 / (s + 3) / (s + 2) / (s + 1)
characteristic_polynomial(g_ol)
# output
Polynomial(10.0 + 11.0*s + 6.0*s^2 + 1.0*s^3)

ξ = damping_coefficient(g)

compute the damping coefficient ξ of an order (0, 2) transfer function.

order (0, 2) representation:

\[g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}\]

returns

ξ::Float64: the damping coefficient

examples

g = 1.0 / (8 * s^2 + 0.8 * s + 2)
damping_coefficient(g)
# output
0.1

evaluate(tf, z)

evaluate a TransferFunction, tf, at a particular number z (could be complex).

example

tf = TransferFunction([1], [3, 1])
evaluate(tf, 1.0)
# output
0.25

g = first_order_system(K, τ)

construct a first-order transfer function with gain K and time constant τ:

\[g(s)=\frac{K}{\tau s+1}\]

example

K = 1.0
τ = 3.0
g = first_order_system(K, τ)
# output
    1.0
-----------
3.0*s + 1.0

returns

• g::TransferFunction: the first order transfer function. well, (0, 1) order.

margins = gain_phase_margins(g_ol, ω_c_guess=0.001, ω_g_guess=0.001)

compute critical frequency (radians / time), gain crossover frequency (radians / time), gain margin, and phase margin (radians) of a closed loop, given its closed loop transfer function g_ol::TransferFunction.

if ωc or ωg is not found (i.e. if either are NaN), but the bode_plot clearly shows a critical/gain crossover frequency, adjust ω_c_guess or ω_g_guess to find the root.

Example

g_ol = 2 * exp(-s) / (5 * s + 1)
margins = gain_phase_margins(g_ol)
margins.ω_c # critical freq. (radians / time)
margins.ω_g # gain crossover freq. (radians / time)
margins.gain_margin # gain margin
margins.phase_margin # phase margin (radians)

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
U = 1 / s                     # input, U(s)
Y = g * U                     # output, Y(s)
data = simulate(g / s, 10.0)  # output, y(t)
y_at_τ = interpolate(data, τ) # ≈ y(τ)
# output
0.6320802858877126

mk_gif(data, title="", xlabel="time, t", 
             ylabel="output, y(t)",
             savename="response")

make a .gif of the process response. data is a data frame with two columns, :t and :output, likely returned from simulate. accepts same arguments as viz_response. ImageMagick must be installed to create the .gif. the .gif is saved as a file savename.

Arguments

• data::DataFrame: data frame of time series data, containing a :t column for times and :output column for the outputs.
• title::String: title of plot
• xlabel::String: x-label
• ylabel::String: y-label
• savename::String: filename to save as a .gif. .gif extension automatically appended if not provided.

nyquist_diagram(tf, nb_pts=500, ω_max=10.0, savename=nothing)

plot the Nyquist diagram for a transfer function tf to visualize its frequency response. s=-1 is plotted as a red +. nb_pts changes the resolution. ω_max gives maximum frequency considered.

returns a CairoMakie.jl Figure object for further modification.

tf = pole_zero_cancellation(tf, verbose=false, digits=8)

find (pole, zero) pairs such that pole = zero and return a new transfer function with those pairs cancelled. this is achieved by comparing the poles and zeros with isapprox, with poles and zeros rounded to digits digits (also applies to reconstruction).

arguments

• tf::TransferFunction: the transfer function
• verbose::Bool=false: print off which poles, zeros are cancelled.
• digits::Int: number of digits to round poles and zeros to, for (i) cancelling and (ii) reconstruction.

example

pole_zero_cancellation(s * (s - 1) / (s * (s + 1)))
# output
1.0*s - 1.0
-----------
1.0*s + 1.0

Convert a polynomial to a string for printing.

proper(tf)

Return true if transfer function tf is proper and false otherwise.

root_locus(g_ol, max_mag_Kc=10.0, nb_pts=500, savename=nothing, legend_pos=:rt)

visualize the root locus plot of an open-loop transfer function g_ol.

Arguments

• g_ol::TransferFunction: the open-loop transfer function of the closed loop system
• max_mag_Kc::Float64=10.0: the maximum magnitude by which the gain of g_ol is scaled in order to see the roots traversing the plane
• nb_pts::Int=500: the number of gains to explore. increase for higher resolution.
• legend_pos::Symbol: Makie command for where to place legend

returns

a CairoMakie.jl Figure object for further modification.

g = second_order_system(K, τ, ξ)

construct a second-order transfer function with gain K, time constant τ, and damping coefficient ξ:

\[g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}\]

example

K = 1.0
τ = 2.0
ξ = 0.1
g = second_order_system(K, τ, ξ)
# output
         1.0
---------------------
4.0*s^2 + 0.4*s + 1.0

returns

• g::TransferFunction: the second order transfer function. well, (0, 2) order.

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::Union{TransferFunction, ClosedLoopTransferFunction}: 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::Union{Int64, 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].

examples

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ᵢ)
first(data, 5)           # show the first 5 rows of the data frame
# output
5×2 DataFrame
 Row │ t          output     
     │ Float64    Float64    
─────┼───────────────────────
   1 │ -0.6       0.0
   2 │ -1.0e-5    0.0
   3 │  1.0e-5    1.33333e-5
   4 │  0.123721  0.161606
   5 │  0.247432  0.316671

strictly_proper(tf)

Return true if transfer function tf is strictly proper and false otherwise.

o = system_order(tf::TransferFunction)

return the order of the numerator and denominator of the transfer function tf.

use pole_zero_cancellation first if you wish to cancel poles and zeros that are equal before determining the order.

returns

o::Tuple{Int, Int}: (order of numerator, order of denominator)

examples

g = 1 / (s + 1)
system_order(g)
# output 
(0, 1)

g = (s + 1) / ((s + 2) * (s + 3))
system_order(g)
# output
(1, 2)

dx –– = A * x(t) + B * u(t) + C * x(t - ϕ) dt y(t) = D * x(t) we compute the matrices A, B, C, D here. only works if strictly proper.

A, B, C = ss_to_tf(tf)

this formulation is based on the notes here: http://web.mit.edu/2.14/www/Handouts/StateSpace.pdf

     Y(s)    bₙ sⁿ + ... + b₁s + b₀

G(s) = ––- = ––––––––––––- U(s) aₙ sⁿ + ... + a₁s + a₀

then we can formulate the system described by transfer function G(s) as a state space representation known as the controllable canonical form:

dx –– = A * x(t) + B * u(t) dt

y(t) = C * x(t) + bₙ/aₙ * u(t)

we compute the matrices A, B, C here.

τ = time_constant(g)

compute the time constant τ of an order (0, 1) or order (0, 2) transfer function.

order (0, 1) representation:

\[g(s)=\frac{K}{\tau s+1}\]

order (0, 2) representation:

\[g(s)=\frac{K}{\tau^2 s^2 + 2\tau \xi s +1}\]

returns

τ::Float64: the time constant.

examples

g = 4 / (6 * s + 2)
time_constant(g)
# output 
3.0

g = 1.0 / (8 * s^2 + 0.8 * s + 2)
time_constant(g) 
# output
2.0

viz_poles_and_zeros(g, savename=nothing, title::String="poles and zeros")

plot the zeros and poles of the transfer function g in the complex plane.

returns a CairoMakie.jl Figure object for further modification.

viz_response(data, 
             title="", xlabel="time, t", 
             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.

Arguments

• data::DataFrame: data frame of time series data, containing a :t column for times and :output column for the outputs.
• title::String: title of plot
• xlabel::String: x-label
• ylabel::String: y-label
• savename::Union{Nothing, String}: filename to save as a figure in .png format.

Returns

CairoMakie.jl Figure object. this will display in a Pluto.jl notebook.

CairoMakie.jl commands can be invoked after viz_response to make further changes to the figure panel by e.g.:

fig = viz_response(data)
ax = current_axis(fig)
ax.xlabel = "new xlabel"
xlims!(ax, 0, 15)

Example

g = 4 / (4 * s ^ 2 + 0.8 * s + 1)
U = 1 / s
Y = g * U
data = simulate(Y, 50.0)
fig = viz_response(data)

K = zero_frequency_gain(tf)

compute the (signed) zero frequency gain of a transfer function $g(s)$, which is:

\[K := \lim_{s\rightarrow 0} G(s)\]

the zero-frequency gain "represents the ratio of the steady state value of the output with respect to a step input" source

example

g = 5 / (3 * s + 1)
K = zero_frequency_gain(g)
# output
5.0

arguments

• tf::TransferFunction: the transfer function

returns

• K::Float64: the zero-frequency gain of the transfer function

z, p, gain = zeros_poles_gain(tf)

Compute the zeros, poles, and zero-frequency gain of a transfer function.

• the zeros are the zeros of the numerator of the transfer function.
• the poles are the zeros of the denominator of the transfer function.
• the zero-frequency gain is the transfer function evaluated at $s=0$

# compute the zeros, poles, and k-factor of a transfer function
z, p, k = zeros_poles_k(tf)
# construct a transfer function from its zeros, poles, and k-factor
tf = zeros_poles_k(z, p, k, time_delay=0.0)

the representation of a transfer function in this context is:

\[g(s)=k\dfrac{\Pi_j (s-z_j)}{\Pi_j (s-p_j)}\]

where $z_j$ is zero $j$, $p_j$ is pole $j$, and $k$ is a constant factor (not equal to the zero-frequency gain) that uniquely specifies the transfer function.

• the zeros are the zeros of the numerator of the transfer function.
• the poles are the zeros of the denominator of the transfer function.

tf = zpk_form(tf)

write transfer function tf in zeros, poles, k-factor form:

\[g(s)=k\dfrac{\Pi_j (s-z_j)}{\Pi_j (s-p_j)}\]

where $z_j$ is zero $j$, $p_j$ is pole $j$, and $k$ is a constant factor (not equal to the zero-frequency gain) that uniquely specifies the transfer function.

this is achieved by multiplying by 1.0 in a fancy way such that the highest power of $s$ in the denominator is associated with a coefficient of $1$.

Example

g = 8.0 / (2 * s^2 + 3 * s + 4)
g_zpk = zpk_form(g)
# output
         4.0
---------------------
1.0*s^2 + 1.5*s + 2.0