# 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 $t_0$: 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 $t=4$ 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

ControlSystems.SimulatorMethod
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")
source
Base.stepMethod
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.

source
ControlSystemsBase.impulseMethod
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)

source
ControlSystemsBase.lsimMethod
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 ommitted, 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 inteface, 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 continuos 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.

Note: 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)
source
ControlSystemsBase.stepinfoMethod
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)
source
ControlSystemsBase.LsimWorkspaceMethod
LsimWorkspace(sys::AbstractStateSpace, N::Int)
LsimWorkspace(sys::AbstractStateSpace, u::AbstractMatrix)
LsimWorkspace{T}(ny, nu, nx, N)

Genereate 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.

source
ControlSystemsBase.StepInfoType
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.
source
ControlSystemsBase.bodeMethod
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.

source
ControlSystemsBase.bodemag!Method
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ω).

source
ControlSystemsBase.evalfrMethod
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.

source
ControlSystemsBase.freqrespMethod
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.

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

Genereate 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ω)
source
ControlSystemsBase.TransferFunctionMethod

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
source