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 systemwhere 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 modelThe 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 modelThis 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 modelzpk - 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 systemwhere 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 modelThe 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 systemand they behave similarly to transfer functions.
The ss constructor allows you to
- Pass
0instead of a $D$ matrix, and an appropriately sized zero matrix is created automatically. - Pass
Iinstead of a $C$ matrix, and an appropriately sized identity matrix is created automatically. TheUniformScalingoperatorIlives in theLinearAlgebrastandard 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 modelNotice 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 modelConverting 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 delayDelayed 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
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 ControlSystemsBasejulia> G = tf([1, 1], [1, 1])TransferFunction{Continuous, ControlSystemsBase.SisoRational{Int64}} s + 1 ----- s + 1 Continuous-time transfer function modeljulia> minreal(G) # Performs pole-zero cancellationTransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}} 1.0 --- 1.0 Continuous-time transfer function modeljulia> P = tf(1, [1, 1]) |> ssStateSpace{Continuous, Float64} A = -1.0 B = 1.0 C = 1.0 D = 0.0 Continuous-time state-space modeljulia> 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 modeljulia> 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 modeljulia> 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 modeljulia> norm(Gmin - feedback(P), Inf) # No difference1.110224134848181e-16julia> 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*P1StateSpace{Continuous, Int64}
A =
-2 1
0 -1
B =
0
1
C =
1 0
D =
0
Continuous-time state-space modelThe 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 systemtrueBroadcasted 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 modelAdding systems
Two systems can be connected in parallel by addition
P12 = P1 + P2StateSpace{Continuous, Int64}
A =
-1 0
0 -2
B =
1
1
C =
1 1
D =
0
Continuous-time state-space modelThe 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 modelis 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 * mimoTest Passed
Thrown: DimensionMismatchTo multiply siso with each output channel of mimo in the example above, use broadcasting:
siso .* mimoStateSpace{Continuous, Float64}
A =
-1.0 0.0 0.37385441724439283 0.862322246448323
0.0 -1.0 -0.8016667112399992 0.7809536785619572
0.0 0.0 -0.08038106672173129 -0.5651317320682765
0.0 0.0 0.15132389642526553 -2.1567059726692532
B =
1.6868266819077053 0.2187308117561154
-0.4105371256956794 -0.7755470761027035
0.5129689411118268 1.773376525151089
0.4737326192557296 1.0012940225824272
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 modelThis is equivalent to first expanding the SISO system into a diagonal system
using LinearAlgebra
(siso .* I(2)) * mimoStateSpace{Continuous, Float64}
A =
-1.0 0.0 0.37385441724439283 0.862322246448323
0.0 -1.0 -0.8016667112399992 0.7809536785619572
0.0 0.0 -0.08038106672173129 -0.5651317320682765
0.0 0.0 0.15132389642526553 -2.1567059726692532
B =
1.6868266819077053 0.2187308117561154
-0.4105371256956794 -0.7755470761027035
0.5129689411118268 1.773376525151089
0.4737326192557296 1.0012940225824272
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 modelConverting 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.09430462988506938 0.0
0.0 -0.20324073901028386
B =
0.49287732786556393 0.0
0.0 -0.4152000481271276
C =
0.11208860790311481 0.0
0.0 1.3192671136140668
D =
-0.048152251418552054 0.0
0.0 -0.9664117560187861
Continuous-time state-space modelMore 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 modelConverting 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.014362973735765672
B =
-1.3066679333573978
C =
0.7049464775809035
D =
-0.7390333218358129
Continuous-time state-space model … StateSpace{Continuous, Float64}
A =
-0.014362973735765672
B =
-0.6157945702241067
C =
0.7049464775809035
D =
0.6398259271248952
Continuous-time state-space model
StateSpace{Continuous, Float64}
A =
-0.014362973735765672
B =
-1.3066679333573978
C =
-0.15477688963221442
D =
-1.441325417634307
Continuous-time state-space model StateSpace{Continuous, Float64}
A =
-0.014362973735765672
B =
-0.6157945702241067
C =
-0.15477688963221442
D =
-0.8130181521737997
Continuous-time state-space modelCreating 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 mismatchTest Passed
Thrown: ErrorExceptionYou 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.009707866809683469
B =
0.35393569584688434 -0.8039793259760136 0.4836185371255346
C =
1.0554397550023749
-0.6344743099373784
D =
0.5315314397048061 0.34793751219901997 0.45155618547215215
-0.013621797134377526 -0.44016042099445496 1.025037208029557
Continuous-time state-space model
StateSpace{Discrete{Int64}, Float64}
A =
0.9
B =
0.13704601517425163 -1.6654678590428913 -2.0159596192845703
C =
1.9454096808860768
-1.2113794275654202
D =
-0.701048425713861 0.8943900646326443 0.2427115330691682
-0.7297849044613821 0.04232175432745541 1.4699799434815548
Sample Time: 1 (seconds)
Discrete-time state-space modelThe 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_sensitivityinput_sensitivityoutput_comp_sensitivityinput_comp_sensitivityG_PSG_CSgangoffour)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