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) # Continuoustime system
tf(num, den, Ts) # Discretetime 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 discretetime system.
Example:
tf([1.0],[1,2,1])
# output
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0

1.0s^2 + 2.0s + 1.0
Continuoustime 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
Continuoustime transfer function model
This allows us to use s
to define transferfunctions:
julia> (s1)*(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
Continuoustime transfer function model
zpk  PoleZeroGain 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) # Continuoustime system
zpk(zeros, poles, gain, Ts) # Discretetime 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)
Continuoustime transfer function model
The transfer functions created using this method will be of type TransferFunction{SisoZpk}
.
StateSpace Systems
A statespace system is created using
ss(A,B,C,D) # Continuoustime system
ss(A,B,C,D,Ts) # Discretetime 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. TheUniformScaling
operatorI
lives in theLinearAlgebra
standard library which must be loaded first.
Statespace 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
Continuoustime statespace 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)
Discretetime transfer function model
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
.
A tutorial on delay systems is available here:
Nonlinear Systems
Simplifying systems
A statespace system with a nonminimal 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 Continuoustime transfer function model
julia> minreal(G) # Performs polezero cancellation
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}} 1.0  1.0 Continuoustime transfer function model
julia> P = tf(1, [1, 1]) > ss
StateSpace{Continuous, Float64} A = 1.0 B = 1.0 C = 1.0 D = 0.0 Continuoustime statespace model
julia> G = P / (1 + P) # this creates a nonminimal realization, use feedback(P) instead
StateSpace{Continuous, Float64} A = 1.0 1.0 0.0 2.0 B = 1.0 1.0 C = 1.0 0.0 D = 0.0 Continuoustime statespace model
julia> feedback(P) # Creates a minimal realization directly
StateSpace{Continuous, Float64} A = 2.0 B = 1.0 C = 1.0 D = 0.0 Continuoustime statespace model
julia> Gmin = minreal(G) # this simplifies the realization to a minimal realization
StateSpace{Continuous, Float64} A = 1.9999999999999993 B = 1.4142135623730951 C = 0.7071067811865472 D = 0.0 Continuoustime statespace model
julia> norm(Gmin  feedback(P), Inf) # No difference
1.110224134848181e16
julia> bodeplot([G, Gmin, feedback(P)]) # They are all identical
Plot{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
Continuoustime statespace model
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
Continuoustime statespace model
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
Continuoustime statespace 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.6301251733700558 0.5440945787509249
0.0 1.0 0.3161165857934258 1.4905406575307747
0.0 0.0 0.08137138135051014 0.7663871194310197
0.0 0.0 0.26732907412161555 2.2761941421919145
B =
0.648183459949459 1.529966154571083
0.955613886564255 0.14067524622311925
0.7839948365963444 0.5643820396125516
2.4152849650581976 1.4732683418250079
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
Continuoustime statespace 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.6301251733700558 0.5440945787509249
0.0 1.0 0.3161165857934258 1.4905406575307747
0.0 0.0 0.08137138135051014 0.7663871194310197
0.0 0.0 0.26732907412161555 2.2761941421919145
B =
0.648183459949459 1.529966154571083
0.955613886564255 0.14067524622311925
0.7839948365963444 0.5643820396125516
2.4152849650581976 1.4732683418250079
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
Continuoustime statespace 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.5905598444774278 0.0
0.0 0.011119569569505924
B =
0.10452978590150008 0.0
0.0 0.4291173451623655
C =
1.241840508516736 0.0
0.0 0.2927855281016749
D =
0.7142973869184909 0.0
0.0 1.2554602788066522
Continuoustime statespace 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
Continuoustime statespace 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.4805382685533392
B =
0.26175783427038896
C =
1.1137397390799957
D =
0.06176235368453573
Continuoustime statespace model … StateSpace{Continuous, Float64}
A =
0.4805382685533392
B =
0.39236691435753857
C =
1.1137397390799957
D =
0.4601995570692816
Continuoustime statespace model
StateSpace{Continuous, Float64}
A =
0.4805382685533392
B =
0.26175783427038896
C =
0.3139465635247748
D =
0.8216046518638513
Continuoustime statespace model StateSpace{Continuous, Float64}
A =
0.4805382685533392
B =
0.39236691435753857
C =
0.3139465635247748
D =
0.4058248425045739
Continuoustime statespace 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 discretetime 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]
2element Vector{StateSpace}:
StateSpace{Continuous, Float64}
A =
0.4340395229534122
B =
1.3663836936700429 0.6028726959153275 0.8968729363401868
C =
0.39641883859610055
0.6255045953296506
D =
0.0580210103680844 0.9530488237065597 0.9958826414184869
1.4539123276582613 1.2488468106142518 1.2284194600208673
Continuoustime statespace model
StateSpace{Discrete{Int64}, Float64}
A =
0.9
B =
1.3172537467090075 0.6175595715259448 1.020432606637834
C =
1.8843334579545388
0.3025716881225607
D =
0.5458313302517364 0.5984941369468755 1.8199819249179046
0.3331445651825328 0.33261536987797535 0.6229213866695263
Sample Time: 1 (seconds)
Discretetime statespace 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.
Closedloop 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