Constructing systems

append(systems::StateSpace...), append(systems::TransferFunction...)

Append systems in block diagonal form

sysd = c2d(sys::AbstractStateSpace{<:Continuous}, Ts, method=:zoh; w_prewarp=0)
Gd = c2d(G::TransferFunction{<:Continuous}, Ts, method=:zoh)

Convert the continuous-time system sys into a discrete-time system with sample time Ts, using the specified method (:zoh, :foh, :fwdeuler or :tustin).

method = :tustin performs a bilinear transform with prewarp frequency w_prewarp.

• w_prewarp: Frequency (rad/s) for pre-warping when using the Tustin method, has no effect for other methods.

See also c2d_x0map

Extended help

ZoH sampling is exact for linear systems with piece-wise constant inputs (step invariant), i.e., the solution obtained using lsim is not approximative (modulu machine precision). ZoH sampling is commonly used to discretize continuous-time plant models that are to be controlled using a discrete-time controller.

FoH sampling is exact for linear systems with piece-wise linear inputs (ramp invariant), this is a good choice for simulation of systems with smooth continuous inputs.

To approximate the behavior of a continuous-time system well in the frequency domain, the :tustin (trapezoidal / bilinear) method may be most appropriate. In this case, the pre-warping argument can be used to ensure that the frequency response of the discrete-time system matches the continuous-time system at a given frequency. The tustin transformation alters the meaning of the state components, while ZoH and FoH preserve the meaning of the state components. The Tustin method is commonly used to discretize a continuous-time controller.

The forward-Euler method generally requires the sample time to be very small relative to the time constants of the system, and its use is generally discouraged.

Classical rules-of-thumb for selecting the sample time for control design dictate that Ts should be chosen as $0.2 ≤ ωgc⋅Ts ≤ 0.6$ where $ωgc$ is the gain-crossover frequency (rad/s).

Qd     = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Ts;             opt=:o)
Qd, Rd = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Rc::Matrix, Ts; opt=:o)
Qd     = c2d(sys::StateSpace{Discrete},   Qc::Matrix;                 opt=:o)
Qd, Rd = c2d(sys::StateSpace{Discrete},   Qc::Matrix, Rc::Matrix;     opt=:o)

Sample a continuous-time covariance or LQR cost matrix to fit the provided discrete-time system.

If opt = :o (default), the matrix is assumed to be a covariance matrix. The measurement covariance R may also be provided. If opt = :c, the matrix is instead assumed to be a cost matrix for an LQR problem.

The method used comes from theorem 5 in the reference below.

Ref: "Discrete-time Solutions to the Continuous-time Differential Lyapunov Equation With Applications to Kalman Filtering", Patrik Axelsson and Fredrik Gustafsson

On singular covariance matrices: The traditional double integrator with covariance matrix Q = diagm([0,σ²]) can not be sampled with this method. Instead, the input matrix ("Cholesky factor") of Q must be manually kept track of, e.g., the noise of variance σ² enters like N = [0, 1] which is sampled using ZoH and becomes Nd = [1/2 Ts^2; Ts] which results in the covariance matrix σ² * Nd * Nd'.

Example:

The following example designs a continuous-time LQR controller for a resonant system. This is simulated with OrdinaryDiffEq to allow the ODE integrator to also integrate the continuous-time LQR cost (the cost is added as an additional state variable). We then discretize both the system and the cost matrices and simulate the same thing. The discretization of an LQR contorller in this way is sometimes refered to as lqrd.

using ControlSystemsBase, LinearAlgebra, OrdinaryDiffEq, Test
sysc = DemoSystems.resonant()
x0 = ones(sysc.nx)
Qc = [1 0.01; 0.01 2] # Continuous-time cost matrix for the state
Rc = I(1)             # Continuous-time cost matrix for the input

L = lqr(sysc, Qc, Rc)
dynamics = function (xc, p, t)
    x = xc[1:sysc.nx]
    u = -L*x
    dx = sysc.A*x + sysc.B*u
    dc = dot(x, Qc, x) + dot(u, Rc, u)
    return [dx; dc]
end
prob = ODEProblem(dynamics, [x0; 0], (0.0, 10.0))
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
cc = sol.u[end][end] # Continuous-time cost

# Discrete-time version
Ts = 0.01 
sysd = c2d(sysc, Ts)
Qd, Rd = c2d(sysd, Qc, Rc, opt=:c)
Ld = lqr(sysd, Qd, Rd)
sold = lsim(sysd, (x, t) -> -Ld*x, 0:Ts:10, x0 = x0)
function cost(x, u, Q, R)
    dot(x, Q, x) + dot(u, R, u)
end
cd = cost(sold.x, sold.u, Qd, Rd) # Discrete-time cost
@test cc ≈ cd rtol=0.01           # These should be similar

c2d(G::DelayLtiSystem, Ts, method=:zoh)

feedback(sys)
feedback(sys1, sys2)

For a general LTI-system, feedback forms the negative feedback interconnection

>-+ sys1 +-->
  |      |
 (-)sys2 +

If no second system is given, negative identity feedback is assumed

feedback(sys1::AbstractStateSpace, sys2::AbstractStateSpace;
         U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[],
         Wperm=:, Zperm=:, pos_feedback::Bool=false)

Basic use feedback(sys1, sys2) forms the (negative) feedback interconnection

           ┌──────────────┐
◄──────────┤     sys1     │◄──── Σ ◄──────
    │      │              │      │
    │      └──────────────┘      -1
    │                            |
    │      ┌──────────────┐      │
    └─────►│     sys2     ├──────┘
           │              │
           └──────────────┘

If no second system sys2 is given, negative identity feedback (sys2 = 1) is assumed. The returned closed-loop system will have a state vector comprised of the state of sys1 followed by the state of sys2.

Advanced use feedback also supports more flexible use according to the figure below

              ┌──────────────┐
      z1◄─────┤     sys1     │◄──────w1
 ┌─── y1◄─────┤              │◄──────u1 ◄─┐
 │            └──────────────┘            │
 │                                        α
 │            ┌──────────────┐            │
 └──► u2─────►│     sys2     ├───────►y2──┘
      w2─────►│              ├───────►z2
              └──────────────┘

• U1, W1 specify the indices of the input signals of sys1 corresponding to u1 and w1. W1 contains the indices of the inputs of sys1 that are included among the inputs to the returned system, i.e., external inputs.
• Y1, Z1 specify the indices of the output signals of sys1 corresponding to y1 and z1. Z1 contains the indices of the outputs ofsys1` that are included among the outputs of the returned system, i.e., external outputs.
• U2, W2, Y2, Z2 specify the corresponding signals of sys2. W2 contains the indices of the inputs ofsys2that are included among the inputs to the returned system, i.e., external inputs.Z2contains the indices of the outputs ofsys2` that are included among the outputs of the returned system, i.e., external outputs.

Specify Wperm and Zperm to reorder the inputs (corresponding to [w1; w2]) and outputs (corresponding to [z1; z2]) in the resulting statespace model.

Negative feedback (α = -1) is the default. Specify pos_feedback=true for positive feedback (α = 1).

See also lft, starprod, sensitivity, input_sensitivity, output_sensitivity, comp_sensitivity, input_comp_sensitivity, output_comp_sensitivity, G_PS, G_CS.

The manual section From block diagrams to code contains higher-level instructions on how to use this function. See also RobustAndOptimalControl.jl: Connections using named signals for a higher-level interface.

See Zhou, Doyle, Glover (1996) for similar (somewhat less symmetric) formulas.

feedback2dof(P,R,S,T)
feedback2dof(B,A,R,S,T)

• Return BT/(AR+ST) where B and A are the numerator and denominator polynomials of P respectively
• Return BT/(AR+ST)

feedback2dof(P, C, F)

Return the transfer function P(F+C)/(1+PC) which is the closed-loop system with process P, controller C and feedforward filter F from reference to control signal (by-passing C).

         +-------+
         |       |
   +----->   F   +----+
   |     |       |    |
   |     +-------+    |
   |     +-------+    |    +-------+
r  |  -  |       |    |    |       |    y
+--+----->   C   +----+---->   P   +---+-->
      |  |       |         |       |   |
      |  +-------+         +-------+   |
      |                                |
      +--------------------------------+

minreal(tf::TransferFunction, eps=sqrt(eps()))

Create a minimal representation of each transfer function in tf by cancelling poles and zeros will promote system to an appropriate numeric type

minreal(sys::StateSpace; fast=false, balance=true, kwargs...)

Minimal realisation algorithm from P. Van Dooreen, The generalized eigenstructure problem in linear system theory, IEEE Transactions on Automatic Control

For information about the options, see ?ControlSystemsBase.MatrixPencils.lsminreal

See also sminreal, which is both numerically exact and substantially faster than minreal, but with a much more limited potential in removing non-minimal dynamics.

parallel(sys1::LTISystem, sys2::LTISystem)

Connect systems in parallel, equivalent to sys2+sys1

series(sys1::LTISystem, sys2::LTISystem)

Connect systems in series, equivalent to sys2*sys1

sminreal(sys)

Compute the structurally minimal realization of the state-space system sys. A structurally minimal realization is one where only states that can be determined to be uncontrollable and unobservable based on the location of 0s in sys are removed.

Systems with numerical noise in the coefficients, e.g., noise on the order of eps require truncation to zero to be affected by structural simplification, e.g.,

trunc_zero!(A) = A[abs.(A) .< 10eps(maximum(abs, A))] .= 0
trunc_zero!(sys.A); trunc_zero!(sys.B); trunc_zero!(sys.C)
sminreal(sys)

In contrast to minreal, which performs pole-zero cancellation using linear-algebra operations, has an 𝑂(nₓ^3) complexity and is subject to numerical tolerances, sminreal is computationally very cheap and numerically exact (operates on integers). However, the ability of sminreal to reduce the order of the model is much less powerful.

See also minreal.

sys = ss(A, B, C, D)      # Continuous
sys = ss(A, B, C, D, Ts)  # Discrete
sys = ss(P::TransferFunction; balance=true, minimal=false) # Convert TransferFunction to StateSpace

Create a state-space model sys::StateSpace{TE, T} with matrix element type T and TE is Continuous or <:Discrete.

This is a continuous-time model if Ts is omitted. Otherwise, this is a discrete-time model with sampling period Ts.

D may be specified as 0 in which case a zero matrix of appropriate size is constructed automatically. sys = ss(D [, Ts]) specifies a static gain matrix D.

When a transfer function P is converted to a state-space model, the user may choose whether to balance the state-space model (default=true) and/or to return a minimal realization (default=false). This conversion has more keyword argument options, see ?ControlSystemsBase.MatrixPencils.rm2ls for additional details.

To associate names with state variables, inputs and outputs, see named_ss.

sys = tf(num, den[, Ts])
sys = tf(gain[, Ts])

Create as a fraction of polynomials:

• sys::TransferFunction{SisoRational{T,TR}} = numerator/denominator

where T is the type of the coefficients in the polynomial.

• num: the coefficients of the numerator polynomial. Either scalar or vector to create SISO systems

or an array of vectors to create MIMO system.

• den: the coefficients of the denominator polynomial. Either vector to create SISO systems

or an array of vectors to create MIMO system.

• Ts: Sample time if discrete time system.

The polynomial coefficients are ordered starting from the highest order term.

Other uses:

• tf(sys): Convert sys to tf form.
• tf("s"), tf("z"): Create the continuous-time transfer function s, or the discrete-time transfer function z.
• numpoly(sys), denpoly(sys): Get the numerator and denominator polynomials of sys as a matrix of vectors, where the outer matrix is of size n_output × n_inputs.

See also: zpk, ss.

zpk(gain[, Ts])
zpk(num, den, k[, Ts])
zpk(sys)

Create transfer function on zero pole gain form. The numerator and denominator are represented by their poles and zeros.

• sys::TransferFunction{SisoZpk{T,TR}} = k*numerator/denominator

where T is the type of k and TR the type of the zeros/poles, usually Float64 and Complex{Float64}.

• num: the roots of the numerator polynomial. Either scalar or vector to create SISO systems

or an array of vectors to create MIMO system.

• den: the roots of the denominator polynomial. Either vector to create SISO systems

or an array of vectors to create MIMO system.

• k: The gain of the system. Obs, this is not the same as dcgain.
• Ts: Sample time if discrete time system.

Other uses:

• zpk(sys): Convert sys to zpk form.
• zpk("s"): Create the transferfunction s.

delay(tau)
delay(tau, Ts)
delay(T::Type{<:Number}, tau)
delay(T::Type{<:Number}, tau, Ts)

Create a pure time delay of length τ of type T.

The type T defaults to promote_type(Float64, typeof(tau)).

If Ts is given, the delay is discretized with sampling time Ts and a discrete-time StateSpace object is returned.

Example:

Create a LTI system with an input delay of L

L = 1
tf(1, [1, 1])*delay(L)
s = tf("s")
tf(1, [1, 1])*exp(-s*L) # Equivalent to the version above

pade(τ::Real, N::Int)

Compute the Nth order Padé approximation of a time-delay of length τ.

See also thiran for discretization of delays.

pade(τ::Real, N_num::Int, N_den::Int)

Compute the Padé approximation of a time-delay of length τ with N_num and N_den degrees in the numerator and denominator, respectively.

pade(G::DelayLtiSystem, N)

Approximate all time-delays in G by Padé approximations of degree N.

thiran(τ::Real, Ts)

Discretize a potentially fractional delay $τ$ as a Thiran all-pass filter with sample time Ts.

The Thiran all-pass filter gives an a maximally flat group delay.

If $τ$ is an integer multiple of $Ts$, the Thiran all-pass filter reduces to $z^{-τ/Ts}$.

Ref: T. I. Laakso, V. Valimaki, M. Karjalainen and U. K. Laine, "Splitting the unit delay [FIR/all pass filters design]," in IEEE Signal Processing Magazine, vol. 13, no. 1, 1996.

A, B, C, D = ssdata(sys)

Outputs the statespace matrices of sys. The matrices are not copies: no new memory is allocated, but modifying the matrices in-place will change the behavior of sys.

Gs, k = seriesform(G::TransferFunction{Discrete})

Convert a transfer function G to a vector of second-order transfer functions and a scalar gain k, the product of which equals G.