Time and Frequency response analysis

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 response of the system sys to a unit step at time t = 0. 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. A SimResul can be plotted by plot(result), 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 response of the system sys to an impulse at time t = 0. For continous-time systems, the impulse is a unit Dirac impulse. For discrete-time systems, the impulse lasts one sample and has magnitude 1/Ts. 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. A SimResul can be plotted by plot(result), or destructured as y, t, x = result.

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

See also lsim.

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. For continuous-time systems, keyword arguments are forwarded to the ODE solver. By default, the option dtmax = t[2]-t[1] is used to prevent the solver from stepping over discontinuities in u(x, t). This prevents the solver from taking too large steps, but may also slow down the simulation when u is smooth. To disable this behavior, set dtmax = Inf.

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.

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)

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