Engee documentation

Time and Frequency response analysis

Frequency response

Frequency responses are calculated using freqresp, bode, sigma and nyquist. Frequency-response plots are obtained using bodeplot, sigmaplot, nyquistplot, marginplot and nicholsplot.

Any TransferFunction can be evaluated at a point using F(s), F(omega, true), F(z, false)

  • F(s) evaluates the continuous-time transfer function F at s.

  • F(omega,true) evaluates the discrete-time transfer function F at exp(i*Ts*omega)

  • F(z,false) evaluates the discrete-time transfer function F at z

A video demonstrating frequency-response analysis in ControlSystems.jl is available below.

Time response (simulation)

Simulation with arbitrary inputs is primarily handled by the function lsim, with step and impulse serving as convenience functions to simulate responses to particular inputs.

The function lsim can take an input vector u containing a sampled input trajectory, or an input function taking the state and time as arguments, u(x,t). This function can be used to easily simulate, e.g., ramp responses or saturated state-feedback control etc. See the docstring of lsim for more details.

For more extensive nonlinear simulation capabilities, see the notes on ModelingToolkit and DifferentialEquations under The wider Julia ecosystem for control.

Example step response:

The following simulates a step response of a second-order system and plots the result.

using ControlSystemsBase, Plots
G = tf(1, [1, 1, 1])
res = step(G, 20) # Simulate 20 seconds step response
plot(res)

Using the function stepinfo, we can compute characteristics of a step response:

si = stepinfo(res)
StepInfo:
Initial value:     0.000
Final value:       1.000
Step size:         1.000
Peak:              1.163
Peak time:         3.652 s
Overshoot:         16.30 %
Undershoot:         0.00 %
Settling time:     8.134 s
Rise time:         1.660 s

We can also plot the StepInfo object

plot(si)

Example lsim:

The function lsim can take the control input as either

  1. An array of equidistantly sampled values, in this case the argument u is expected to have the shape nu × n_time

  2. A function of the state and time u(x,t). This form allows simulation of state feedback, a step response at time : u(x, t) = amplitude * (t > t0), or a ramp response: u(x, t) = t etc.

The example below simulates state feedback with a step disturbance at by providing the function u(x,t) = -L*x .+ (t > 4) to lsim:

using ControlSystems
using LinearAlgebra: I
using Plots

A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x .+ (t > 4) # State feedback + step disturbance
t  = 0:0.1:12
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]"); vline!([4], lab="Step disturbance", l=(:black, :dash, 0.5))

A video demonstrating time-domain simulation in ControlSystems.jl is available below.

Docstrings

res = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), abstol=1e-6, reltol=1e-6, force_dtmin=true, kwargs...)

Simulate system sys, over time t, using input signal u, with initial state x0, using method alg .

Arguments:

t: Has to be an AbstractVector with equidistant time samples (t[i] - t[i-1] constant) u: Function to determine control signal ut at a time t, on any of the following forms:

  • u: Function to determine control signal uₜ at a time t, on any of the following forms:

    • A constant Number or Vector, interpreted as a constant input.

    • Function u(x, t) that takes the internal state and time, note, the state representation for delay systems is not the same as for rational systems.

    • In-place function u(uₜ, x, t). (Slightly more efficient)

alg, abstol, reltol and kwargs...: are sent to DelayDiffEq.solve.

This methods sets force_dtmin=true by default to handle the discontinuity implied by, e.g., step inputs. This may lead to the solver taking a long time to solve ill-conditioned problems rather than exiting early with a warning.

Returns an instance of SimResult which can be plotted directly or destructured into y, t, x, u = res.

lsim(sys::HammersteinWienerSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=Tsit5(), abstol=1e-6, reltol=1e-6, kwargs...)

Simulate system sys, over time t, using input signal u, with initial state x0, using method alg .

Arguments:

  • t: Has to be an AbstractVector with equidistant time samples (t[i] - t[i-1] constant)

  • u: Function to determine control signal uₜ at a time t, on any of the following forms: Can be a constant Number or Vector, interpreted as uₜ .= u , or Function uₜ .= u(x, t), or In-place function u(uₜ, x, t). (Slightly more efficient)

  • alg, abstol, reltol and kwargs...: are sent to OrdinaryDiffEq.solve.

Returns an instance of SimResult.

Simulator

Fields:

P::StateSpace
f = (x,p,t) -> x
y = (x,t)   -> y
Simulator(P::StateSpace, u = (x,t) -> 0)

Used to simulate continuous-time systems. See function ?solve for additional info.

Usage:

using OrdinaryDiffEq, Plots
dt             = 0.1
tfinal         = 20
t              = 0:dt:tfinal
P              = ss(tf(1,[2,1])^2)
K              = 5
reference(x,t) = [1.]
s              = Simulator(P, reference)
x0             = [0.,0]
tspan          = (0.0,tfinal)
sol            = solve(s, x0, tspan, Tsit5())
plot(t, s.y(sol, t)[:], lab="Open loop step response")
y, t, x = step(sys[, tfinal])
y, t, x = step(sys[, t])

Calculate the step response of system sys. If the final time tfinal or time vector t is not provided, one is calculated based on the system pole locations. The return value is a structure of type SimResult that can be plotted or destructured as y, t, x = result.

y has size (ny, length(t), nu), x has size (nx, length(t), nu)

See also stepinfo and lsim.

y, t, x = impulse(sys[, tfinal])
y, t, x = impulse(sys[, t])

Calculate the impulse response of system sys. If the final time tfinal or time vector t is not provided, one is calculated based on the system pole locations. The return value is a structure of type SimResult that can be plotted or destructured as y, t, x = result.

y has size (ny, length(t), nu), x has size (nx, length(t), nu)

res = lsim!(ws::LsimWorkspace, sys::AbstractStateSpace{<:Discrete}, u, [t]; x0)

In-place version of lsim that takes a workspace object created by calling LsimWorkspace. Notice, if u is a function, res.u === ws.u. If u is an array, res.u === u.

result = lsim(sys, u[, t]; x0, method])
result = lsim(sys, u::Function, t; x0, method)

Calculate the time response of system sys to input u. If x0 is omitted, a zero vector is used.

The result structure contains the fields y, t, x, u and can be destructured automatically by iteration, e.g.,

y, t, x, u = result

result::SimResult can also be plotted directly:

plot(result, plotu=true, plotx=false)

y, x, u have time in the second dimension. Initial state x0 defaults to zero.

Continuous-time systems are simulated using an ODE solver if u is a function (requires using ControlSystems). If u is an array, the system is discretized (with method=:zoh by default) before simulation. For a lower-level interface, see ?Simulator and ?solve

u can be a function or a matrix of precalculated control signals and must have dimensions (nu, length(t)). If u is a function, then u(x,i) (for discrete systems) or u(x,t) (for continuous ones) is called to calculate the control signal at every iteration (time instance used by solver). This can be used to provide a control law such as state feedback u(x,t) = -L*x calculated by lqr. To simulate a unit step at t=t₀, use (x,t)-> t ≥ t₀, for a ramp, use (x,t)-> t, for a step at t=5, use (x,t)-> (t >= 5) etc.

The function u will be called once before simulating to verify that it returns an array of the correct dimensions. This can cause problems if u is stateful. You can disable this check by passing check_u = false.

For maximum performance, see function lsim!, available for discrete-time systems only.

Usage example:

using ControlSystems
using LinearAlgebra: I
using Plots

A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x # Form control law
t  = 0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")

# Alternative way of plotting
res = lsim(sys,u,t,x0=x0)
plot(res)
stepinfo(res::SimResult; y0 = nothing, yf = nothing, settling_th = 0.02, risetime_th = (0.1, 0.9))

Compute the step response characteristics for a simulation result. The following information is computed and stored in a StepInfo struct:

  • y0: The initial value of the response

  • yf: The final value of the response

  • stepsize: The size of the step

  • peak: The peak value of the response

  • peaktime: The time at which the peak occurs

  • overshoot: The percentage overshoot of the response

  • undershoot: The percentage undershoot of the response. If the step response never reaches below the initial value, the undershoot is zero.

  • settlingtime: The time at which the response settles within settling_th of the final value

  • settlingtimeind: The index at which the response settles within settling_th of the final value

  • risetime: The time at which the response rises from risetime_th[1] to risetime_th[2] of the final value

Arguments:

  • res: The result from a simulation using step (or lsim)

  • y0: The initial value, if not provided, the first value of the response is used.

  • yf: The final value, if not provided, the last value of the response is used. The simulation must have reached steady-state for an automatically computed value to make sense. If the simulation has not reached steady state, you may provide the final value manually.

  • settling_th: The threshold for computing the settling time. The settling time is the time at which the response settles within settling_th of the final value.

  • risetime_th: The lower and upper threshold for computing the rise time. The rise time is the time at which the response rises from risetime_th[1] to risetime_th[2] of the final value.

Example:

G = tf([1], [1, 1, 1])
res = step(G, 15)
si = stepinfo(res)
plot(si)
LsimWorkspace(sys::AbstractStateSpace, N::Int)
LsimWorkspace(sys::AbstractStateSpace, u::AbstractMatrix)
LsimWorkspace{T}(ny, nu, nx, N)

Generate a workspace object for use with the in-place function lsim!. sys is the discrete-time system to be simulated and N is the number of time steps, alternatively, the input u can be provided instead of N. Note: for threaded applications, create one workspace object per thread.

StepInfo

Computed using stepinfo

Fields:

  • y0: The initial value of the step response.

  • yf: The final value of the step response.

  • stepsize: The size of the step.

  • peak: The peak value of the step response.

  • peaktime: The time at which the peak occurs.

  • overshoot: The overshoot of the step response.

  • settlingtime: The time at which the step response has settled to within settling_th of the final value.

  • settlingtimeind::Int: The index at which the step response has settled to within settling_th of the final value.

  • risetime: The time at which the response rises from risetime_th[1] to risetime_th[2] of the final value

  • i10::Int: The index at which the response reaches risetime_th[1]

  • i90::Int: The index at which the response reaches risetime_th[2]

  • res::SimResult{SR}: The simulation result used to compute the step response characteristics.

  • settling_th: The threshold used to compute settlingtime and settlingtimeind.

  • risetime_th: The thresholds used to compute risetime, i10, and i90.

mag, phase, w = bode(sys[, w]; unwrap=true)

Compute the magnitude and phase parts of the frequency response of system sys at frequencies w. See also bodeplot

mag and phase has size (ny, nu, length(w)). If unwrap is true (default), the function unwrap! will be applied to the phase angles. This procedure is costly and can be avoided if the unwrapping is not required.

For higher performance, see the function bodemag! that computes the magnitude only.

mag = bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)

Compute the Bode magnitude operating in-place on an instance of BodemagWorkspace. Note that the returned magnitude array is aliased with ws.mag. The output array mag is ∈ 𝐑(ny, nu, nω).

evalfr(sys, x)

Evaluate the transfer function of the LTI system sys at the complex number s=x (continuous-time) or z=x (discrete-time).

For many values of x, use freqresp instead.

freqresp!(R::Array{T, 3}, sys::LTISystem, w_vec::AbstractVector{<:Real})

In-place version of freqresp that takes a pre-allocated array R of size (ny, nu, nw)`

sys_fr = freqresp(sys, w)

Evaluate the frequency response of a linear system

w -> C*((iw*im*I - A)^-1)*B + D

of system sys over the frequency vector w.

re, im, w = nyquist(sys[, w])

Compute the real and imaginary parts of the frequency response of system sys at frequencies w. See also nyquistplot

re and im has size (ny, nu, length(w))

sv, w = sigma(sys[, w])

Compute the singular values sv of the frequency response of system sys at frequencies w. See also sigmaplot

sv has size (min(ny, nu), length(w))

BodemagWorkspace(sys::LTISystem, N::Int)
BodemagWorkspace(sys::LTISystem, ω::AbstractVector)
BodemagWorkspace(R::Array{Complex{T}, 3}, mag::Array{T, 3})
BodemagWorkspace{T}(ny, nu, N)

Generate a workspace object for use with the in-place function bodemag!. N is the number of frequency points, alternatively, the input ω can be provided instead of N. Note: for threaded applications, create one workspace object per thread.

Arguments:

  • mag: The output array ∈ 𝐑(ny, nu, nω)

  • R: Frequency-response array ∈ 𝐂(ny, nu, nω)

F(s), F(omega, true), F(z, false)

Notation for frequency response evaluation.

  • F(s) evaluates the continuous-time transfer function F at s.

  • F(omega,true) evaluates the discrete-time transfer function F at exp(imTsomega)

  • F(z,false) evaluates the discrete-time transfer function F at z