Creating Systems

This page illustrates how to create system models such as transfer functions and statespace models. This topic is also treated in the introductory video below:

Transfer Functions

tf - Rational Representation

The basic syntax for creating a transfer function is tf

tf(num, den)     # Continuous-time system
tf(num, den, Ts) # Discrete-time system

where num and den are the polynomial coefficients of the numerator and denominator of the polynomial and Ts, if provided, is the sample time for a discrete-time system.

Example:

tf([1.0],[1,2,1])

# output

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
        1.0
-------------------
1.0s^2 + 2.0s + 1.0

Continuous-time transfer function model

The transfer functions created using this method will be of type TransferFunction{SisoRational}. For more general expressions, it is often more convenient to define s = tf("s"):

Example:

julia> s = tf("s")

TransferFunction{Continuous,ControlSystems.SisoRational{Int64}}
s
-
1

Continuous-time transfer function model

This allows us to use s to define transfer-functions:

julia> (s-1)*(s^2 + s + 1)/(s^2 + 3s + 2)/(s+1)

TransferFunction{Continuous,ControlSystems.SisoRational{Int64}}
       s^3 - 1
---------------------
s^3 + 4*s^2 + 5*s + 2

Continuous-time transfer function model

zpk - Pole-Zero-Gain Representation

Sometimes it's better to represent the transfer function by its poles, zeros and gain, this can be done using the function zpk

zpk(zeros, poles, gain)     # Continuous-time system
zpk(zeros, poles, gain, Ts) # Discrete-time system

where zeros and poles are Vectors of the zeros and poles for the system and gain is a gain coefficient.

Example

zpk([-1.0,1], [-5, -10], 2)

# output

TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, Float64}}
   (1.0s + 1.0)(1.0s - 1.0)
2.0-------------------------
   (1.0s + 5.0)(1.0s + 10.0)

Continuous-time transfer function model

The transfer functions created using this method will be of type TransferFunction{SisoZpk}.

State-Space Systems

A state-space system

\[\begin{aligned} \dot{x} &= Ax + Bu \\ y &= Cx + Du \end{aligned}\]

in continuous time, or

\[\begin{aligned} x_{t+T_s} &= Ax_t + Bu_t \\ y_t &= Cx_t + Du_t \end{aligned}\]

in discrete time, is created using

ss(A,B,C,D)    # Continuous-time system
ss(A,B,C,D,Ts) # Discrete-time system

and they behave similarly to transfer functions.

The ss constructor allows you to

• Pass 0 instead of a $D$ matrix, and an appropriately sized zero matrix is created automatically.
• Pass I instead of a $C$ matrix, and an appropriately sized identity matrix is created automatically. The UniformScaling operator I lives in the LinearAlgebra standard library which must be loaded first.

State-space systems with heterogeneous matrix types are also available, which can be used to create systems with static or sized matrices, e.g.,

using ControlSystemsBase, StaticArrays
sys = ss([-5 0; 0 -5],[2; 2],[3 3],[0])
HeteroStateSpace(sys, to_sized)
HeteroStateSpace(sys, to_static)

StaticStateSpace{Continuous, 1, 1, 2, StaticArraysCore.SMatrix{2, 2, Int64, 4}, StaticArraysCore.SMatrix{2, 1, Int64, 2}, StaticArraysCore.SMatrix{1, 2, Int64, 2}, StaticArraysCore.SMatrix{1, 1, Int64, 1}}
A = 
 -5   0
  0  -5
B = 
 2
 2
C = 
 3  3
D = 
 0

Continuous-time state-space model

Notice the different matrix types used.

To associate names with states, inputs and outputs, see named_ss from RobustAndOptimalControl.jl.

Converting between types

It is sometime useful to convert one representation to another. This is possible using the constructors tf, zpk, ss, for example

tf(zpk([-1], [1], 2, 0.1))

# output

TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Int64}}
2z + 2
------
z - 1

Sample Time: 0.1 (seconds)
Discrete-time transfer function model

Converting between continuous and discrete time

A continuous-time system represents differential equations or a transfer function in the Laplace domain, while a discrete-time system represents difference equations or a transfer function in the Z-domain.

The functions c2d and d2c implement sampling/discretization of continuous-time systems and the inverse mapping from discrete-time to continuous-time systems.

Delay Systems

The constructor delay creates a pure delay, which may be connected to a system by multiplication:

delay(1.2)               # Pure delay or 1.2s
tf(1, [1, 1])*delay(1.2) # Input delay
delay(1.2)*tf(1, [1, 1]) # Output delay

Delayed systems can also be created using

s = tf("s")
L = 1.2 # Delay time
tf(1, [1, 1]) * exp(-L*s)

Padé approximations of delays can be created using pade. Models with delays can be discretized using c2d, currently, only delays that are integer multiples of the sample time are supported. Pure fractional delays can be approximately discretized using the function thiran.

A tutorial on delay systems is available here:

Nonlinear Systems

See Nonlinear functionality.

Simplifying systems

A statespace system with a non-minimal realization, or a transfer function with overlapping zeros and poles, may be simplified using the function minreal. Systems that are structurally singular, i.e., that contains outputs that can not be reached from the inputs based on analysis of the structure of the zeros in the system matrices only, can be simplified with the function sminreal.

Examples:

julia> using ControlSystemsBase
julia> G = tf([1, 1], [1, 1])TransferFunction{Continuous, ControlSystemsBase.SisoRational{Int64}}
s + 1
-----
s + 1

Continuous-time transfer function model
julia> minreal(G) # Performs pole-zero cancellationTransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0
---
1.0

Continuous-time transfer function model
julia> P = tf(1, [1, 1]) |> ssStateSpace{Continuous, Float64}
A =
 -1.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model
julia> G = P / (1 + P) # this creates a non-minimal realization, use feedback(P) insteadStateSpace{Continuous, Float64}
A =
 -2.0
B =
 1.4142135623730954
C =
 0.7071067811865476
D =
 0.0

Continuous-time state-space model
julia> feedback(P) # Creates a minimal realization directlyStateSpace{Continuous, Float64}
A =
 -2.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model
julia> Gmin = minreal(G) # this simplifies the realization to a minimal realizationStateSpace{Continuous, Float64}
A =
 -2.0
B =
 1.4142135623730954
C =
 0.7071067811865476
D =
 0.0

Continuous-time state-space model
julia> norm(Gmin - feedback(P), Inf) # No difference1.110224134848181e-16
julia> bodeplot([G, Gmin, feedback(P)]) # They are all identicalPlot{Plots.GRBackend() n=6}

Multiplying systems

Two systems can be connected in series by multiplication

using ControlSystemsBase
P1 = ss(-1,1,1,0)
P2 = ss(-2,1,1,0)
P2*P1

StateSpace{Continuous, Int64}
A = 
 -2   1
  0  -1
B = 
 0
 1
C = 
 1  0
D = 
 0

Continuous-time state-space model

The state of the resulting system is the concatenation of the states of the two systems, starting with the left/first operand (P2 above).

If the input dimension of P2 does not match the output dimension of P1, an error is thrown. If one of the systems is SISO and the other is MIMO, broadcasted multiplication will expand the SISO system to match the input or output dimension of the MIMO system, e.g.,

Pmimo = ssrand(2,2,1)
Psiso = ss(-2,1,1,0)
# Psiso * Pmimo # error
Psiso .* Pmimo ≈ [Psiso 0; 0 Psiso] * Pmimo # Broadcasted multiplication expands SISO into diagonal system

true

Broadcasted multiplication between a system and an array is only allowed for diagonal arrays

using LinearAlgebra
Psiso .* I(2)

StateSpace{Continuous, Int64}
A = 
 -2   0
  0  -2
B = 
 1  0
 0  1
C = 
 1  0
 0  1
D = 
 0  0
 0  0

Continuous-time state-space model

Adding systems

Two systems can be connected in parallel by addition

P12 = P1 + P2

StateSpace{Continuous, Int64}
A = 
 -1   0
  0  -2
B = 
 1
 1
C = 
 1  1
D = 
 0

Continuous-time state-space model

The state of the resulting system is the concatenation of the states of the two systems, starting with the left/first operand (P1 above).

MIMO systems and arrays of systems

Concatenation of systems creates MIMO systems, which is different from an array of systems. For example

using ControlSystemsBase
P = ss(-1,1,1,0)
P_MIMO = [P 2P]

StateSpace{Continuous, Int64}
A = 
 -1   0
  0  -1
B = 
 1  0
 0  1
C = 
 1  2
D = 
 0  0

Continuous-time state-space model

is a 1×2 MISO system, not a 1×2 array.

From SISO to MIMO

SISO systems do not multiply MIMO systems directly, i.e.,

using Test
siso = ss(-1,1,1,0)
mimo = ssrand(2,2,2)
@test_throws DimensionMismatch siso * mimo

Test Passed
      Thrown: DimensionMismatch

To multiply siso with each output channel of mimo in the example above, use broadcasting:

siso .* mimo

StateSpace{Continuous, Float64}
A = 
 -1.0   0.0  -0.814115710635757    2.05611954192475
  0.0  -1.0  -0.309671182140287    0.8213151945608457
  0.0   0.0  -0.4615064500828379  -2.3246447006361297
  0.0   0.0   1.214054438269907    0.40767445019951687
B = 
 -0.4465259344260694   0.8557846776610339
 -1.074561991927812   -1.898978885310723
  2.7521659903912368  -0.2210180097529923
 -1.1205600313664303   1.2754295013740533
C = 
 1.0  0.0  -0.0  0.0
 0.0  1.0  -0.0  0.0
D = 
 -0.0  0.0
 -0.0  0.0

Continuous-time state-space model

This is equivalent to first expanding the SISO system into a diagonal system

using LinearAlgebra
(siso .* I(2)) * mimo

StateSpace{Continuous, Float64}
A = 
 -1.0   0.0  -0.814115710635757    2.05611954192475
  0.0  -1.0  -0.309671182140287    0.8213151945608457
  0.0   0.0  -0.4615064500828379  -2.3246447006361297
  0.0   0.0   1.214054438269907    0.40767445019951687
B = 
 -0.4465259344260694   0.8557846776610339
 -1.074561991927812   -1.898978885310723
  2.7521659903912368  -0.2210180097529923
 -1.1205600313664303   1.2754295013740533
C = 
 1.0  0.0  -0.0  0.0
 0.0  1.0  -0.0  0.0
D = 
 -0.0  0.0
 -0.0  0.0

Continuous-time state-space model

Converting an array of systems to a MIMO system

Diagonal MIMO systems can be created from a vector of systems using append

P1 = ssrand(1,1,1)
P2 = ssrand(1,1,1)
append(P1, P2)

StateSpace{Continuous, Float64}
A = 
 -0.6071761154690204   0.0
  0.0                 -0.22966324841010835
B = 
 -0.31855392586848136   0.0
  0.0                  -0.3665955381662163
C = 
 0.5776702975692982   0.0
 0.0                 -0.48894241370063657
D = 
 0.7519721567302424   0.0
 0.0                 -0.12899827293248062

Continuous-time state-space model

More general arrays of systems can be converted to a MIMO system using array2mimo.

sys_array = fill(P, 2, 2) # Creates an array of systems
mimo_sys = array2mimo(sys_array)

StateSpace{Continuous, Int64}
A = 
 -1   0   0   0
  0  -1   0   0
  0   0  -1   0
  0   0   0  -1
B = 
 1  0
 0  1
 1  0
 0  1
C = 
 1  1  0  0
 0  0  1  1
D = 
 0  0
 0  0

Continuous-time state-space model

Converting MIMO system to an array of systems

This conversion is not explicitly supported, but is easy enough to accomplish with standard Julia code, for example:

P = ssrand(2,3,1) # A random 2×3 MIMO system
sys_array = getindex.(Ref(P), 1:P.ny, (1:P.nu)')

2×3 Matrix{StateSpace{Continuous, Float64}}:
 StateSpace{Continuous, Float64}
A = 
 -0.08691007291327735
B = 
 -0.5806774863172849
C = 
 0.5737119614126563
D = 
 1.0096819268508332

Continuous-time state-space model  …  StateSpace{Continuous, Float64}
A = 
 -0.08691007291327735
B = 
 1.324518848959261
C = 
 0.5737119614126563
D = 
 1.0462125553246147

Continuous-time state-space model
 StateSpace{Continuous, Float64}
A = 
 -0.08691007291327735
B = 
 -0.5806774863172849
C = 
 0.012172225030040618
D = 
 0.86166223315889

Continuous-time state-space model     StateSpace{Continuous, Float64}
A = 
 -0.08691007291327735
B = 
 1.324518848959261
C = 
 0.012172225030040618
D = 
 0.26750937193882823

Continuous-time state-space model

Creating arrays with different types of systems

When calling hcat/vcat, Julia automatically tries to promote the types to the smallest common supertype, this means that creating an array with one continuous and one discrete-time system fails

P_cont = ssrand(2,3,1)
P_disc = ssrand(2,3,1, Ts=1)
@test_throws ErrorException [P_cont, P_disc] # ERROR: Sampling time mismatch

Test Passed
      Thrown: ErrorException

You can explicitly tell Julia that you want a particular supertype, e.g,

StateSpace[P_cont, P_disc]

2-element Vector{StateSpace}:
 StateSpace{Continuous, Float64}
A = 
 -0.06682087565482409
B = 
 1.4946819804967477  -1.7753277533475647  0.06871071737381562
C = 
 1.1058131091255579
 1.4554082338585987
D = 
 -0.3690864860323875   0.38877368689346825  -0.5381864459037952
  0.9468404821940324  -0.18773341014198228  -0.37455239368384247

Continuous-time state-space model
 StateSpace{Discrete{Int64}, Float64}
A = 
 0.9
B = 
 -0.8947598765146301  2.824305648773669  1.8515305564795363
C = 
 -0.9997315444062458
 -0.2682635474824611
D = 
  0.7792688578785669  -0.4184642057104451  -1.3735905782913203
 -1.320275850416764   -0.5373239855237092   0.6193200509609528

Sample Time: 1 (seconds)
Discrete-time state-space model

The type StateSpace is abstract, since the type parameters are not specified.

Demo systems

The module ControlSystemsBase.DemoSystems contains a number of demo systems demonstrating different kinds of dynamics.

From block diagrams to code

This section lists a number of block diagrams, and indicates the corresponding transfer functions and how they are built in code.

The function feedback(G1, G2) can be thought of like this: the first argument G1 is the system that appears directly between the input and the output (the forward path), while the second argument G2 (defaults to 1 if omitted) contains all other systems that appear in the closed loop (the feedback path). The feedback is assumed to be negative, unless the argument pos_feedback = true is passed (lft is an exception, which due to convention defaults to positive feedback). This means that feedback(G, 1) results in unit negative feedback, while feedback(G, -1) or feedback(G, 1, pos_feedback = true) results in unit positive feedback.

The returned closed-loop system will have a state vector comprised of the state of G1 followed by the state of G2.

Closed-loop system from reference to output

r   ┌─────┐     ┌─────┐
───►│     │  u  │     │ y
    │  C  ├────►│  P  ├─┬─►
 -┌►│     │     │     │ │
  │ └─────┘     └─────┘ │
  │                     │
  └─────────────────────┘

\[Y = \dfrac{PC}{I+PC}R\]

Code: feedback(P*C) or equivalently comp_sensitivity(P, C). Here, the system $PC$ appears directly between the input $r$ and the output $y$, and the feedback loop is negative identity. We thus call feedback(P*C) = feedback(P*C, 1)

d     ┌───┐   y
───+─►│ P ├─┬───►
  -▲  └───┘ │
   │        │
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

\[Y = \dfrac{P}{I+PC}D = PSD\]

Code: feedback(P, C) or equivalently G_PS(P, C). Here, only $P$ appears directly between $d$ and $y$, while $C$ appears first in the feedback loop. We thus call feedback(P, C)

Sensitivity function at plant input

d    e┌───┐   
───+─►│ P ├─┬───►
  -▲  └───┘ │
   │        │
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

\[E = \dfrac{1}{I+CP}D = SD\]

Code: feedback(1, C*P) or equivalently input_sensitivity(P, C). Here, there are no systems directly between the input and the output, we thus call feedback(1, C*P). Note the order in C*P, which is important for MIMO systems. This computes the sensitivity function at the plant input. It's more common to analyze the sensitivity function at the plant output, illustrated below (for SISO systems they are equivalent).

Sensitivity function at plant output

      ┌───┐   
───+─►│ P ├─+◄── e
  -▲  └───┘ │
   │        │y
   │  ┌───┐ │
   └──┤ C │◄┘
      └───┘

\[Y = \dfrac{1}{I+PC}E = SE\]

Code: feedback(1, P*C) or equivalently output_sensitivity(P, C). Note the reverse order in $PC$ compared to the input sensitivity function above.

Reference $r$ and input disturbance $d$ to output $y$ and control signal $u$. This example forms the transfer function matrix with $r$ and $d$ as inputs, and $y$ and $u$ as outputs.

              d
     ┌─────┐  │  ┌─────┐
r    │     │u ▼  │     │ y
──+─►│  C  ├──+─►│  P  ├─┬─►
  ▲  │     │     │     │ │
 -│  └─────┘     └─────┘ │
  │                      │
  └──────────────────────┘

\[\begin{bmatrix} y \\ u \end{bmatrix} = \begin{bmatrix} \dfrac{PC}{I + PC} & \dfrac{C}{I + PC} \\ \dfrac{P}{I + PC} & \dfrac{-PC}{I + PC} \end{bmatrix} \begin{bmatrix} r \\ d \end{bmatrix}\]

Code: feedback(C, P, W2=:, Z2=:, Zperm=[(1:P.ny).+P.nu; 1:P.nu]) # y,u from r,d. Here, we have reversed the order of P and C to get the correct sign of the control signal. We also make use of the keyword arguments W2 and Z2 to specify that we want to include the inputs and outputs of P as external inputs and outputs, and Zperm to specify the order of the outputs ($y$ before $u$).

Linear fractional transformation

     ┌─────────┐
z◄───┤         │◄────w
     │    P    │
y┌───┤         │◄───┐u
 │   └─────────┘    │
 │                  │
 │      ┌───┐       │
 │      │   │       │
 └─────►│ K ├───────┘
        │   │
        └───┘

\[Z = \operatorname{lft}{(P, K)} W\]

Code: lft(P, K)

      z1          z2
      ▲  ┌─────┐  ▲      ┌─────┐
      │  │     │  │      │     │
w1──+─┴─►│  C  ├──┴───+─►│  P  ├─┐
    │    │     │      │  │     │ │
    │    └─────┘      │  └─────┘ │
    │                 w2         │
    └────────────────────────────┘

The transfer function from $w_1, w_2$ to $z_1, z_2$ contains all the transfer functions that are commonly called "gang of four" (see also gangoffour).

\[\begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} I \\ C \end{bmatrix} (I + PC)^{-1} \begin{bmatrix} I & P \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \end{bmatrix}\]

Code: This function requires the package RobustAndOptimalControl.jl.

RobustAndOptimalControl.extended_gangoffour(P, C, pos=true)
# For SISO P
S  = G[1, 1]
PS = G[1, 2]
CS = G[2, 1]
T  = G[2, 2]

# For MIMO P
S  = G[1:P.ny,     1:P.nu]
PS = G[1:P.ny,     P.nu+1:end]
CS = G[P.ny+1:end, 1:P.nu]
T  = G[P.ny+1:end, P.nu+1:end]

See also

• output_sensitivity
• input_sensitivity
• output_comp_sensitivity
• input_comp_sensitivity
• G_PS
• G_CS
• gangoffour)
• gangoffourplot)

This diagram is more complicated and forms several connections, including both feedforward and feedback connections. A code file that goes through how to form such complicated connections using named signals is linked below. This example uses the package RobustAndOptimalControl.jl.

                 yF
              ┌────────────────────────────────┐
              │                                │
    ┌───────┐ │  ┌───────┐ yR   ┌─────────┐    │    ┌───────┐
uF  │       │ │  │       ├──────►         │ yC │  uP│       │    yP
────►   F   ├─┴──►   R   │      │    C    ├────+────►   P   ├────┬────►
    │       │    │       │   ┌──►         │         │       │    │
    └───────┘    └───────┘   │- └─────────┘         └───────┘    │
                             │                                   │
                             └───────────────────────────────────┘

See code example complicated_feedback.jl.

Filter design

Filters can be designed using DSP.jl. This results in filter objects with types from the DSP package, which can be converted to transfer functions using tf from ControlSystemsBase.

using DSP, ControlSystemsBase, Plots

fs = 100
df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2))
G = tf(df, 1/fs) # Sample time must be provided in the conversion to get the correct frequency scale in the Bode plot
bodeplot(G, xscale=:identity, yscale=:identity, hz=true)

See also

• ControlSystemsBase.seriesform