Tutorial for
Control System Toolbox for MATLAB
by
This publication can be downloaded and copied
freely, but reference to the source is appreciated. Contents:
1 Introduction 1 Introduction
This text gives an easy
guide to Control System Toolbox. The text is
self-instructive: You are asked to perform a number of simpe tasks through
which you will learn to master this toolbox, and the expected responses are
shown in the text. In the following, CST
will be used as an abbreviation for Control System Toolbox. 2 Creating and handling LTI-models (Linear
Time-Invariant)
The models supported by
CST are continuous-time models and discrete-time models of the following
types:
The CST also supports
zero-pole-gain models, but since we will in practice hardly use that model
format, it is not described in this text. If you need to use it, you can get
information by executing "help zpk". The various functions
will be introduced via simple examples. 2.1 Creating models
We will start with
continuous-time models, and then take discrete-time models. Continuous-time transfer functions
The function
"tf" creates transfer functions. "tf" needs two
MATLAB-vectors, one containing the coefficients of the numerator
polynomial - taken in descending orders - and one for the denominator
polynomial of the transfer function. As an example we will create the
transfer function H0(s)=K0*s/(T0*s+1)
To display H0 we execute
MATLAB's response is
An LTI-model can contain
a time-delay, as in the following model: Hd(s)=[K0*s/(T0*s+1)]e^(-Tdelay*s) This system is created by
(except from the time-delay term, Hd is equal to H0 defined above):
To display Hd we execute
MATLAB's response is
Alternatively, you can define
the symbol ‘s’ as the Laplace transform variable, and then write transfer
functions as you would do in a document:
MATLAB’s response is Hd1 = 2 s exp(-1*s) * ------- 4 s + 1 Continuous-time transfer
function. The expressions below
create a two-input one-output MIMO (Multiple-Input Multple-Output)
system of the form y(s)=H1(s) u1(s) + H2(s) u2(s) The transfer function can
be expressed as H=[H1,H2]
To display H we execute
MATLAB's response is
Note 1: If the system has
no time-delay, it is sufficient to write tf(num1,den1). Note 2: If the transfer
function has the form H(s)=[H1(s) which may correspond to a
one-input two-output system, the CST requires that the time-delays are equal.
Suppose we have created a
transfer function, and we want to retrieve detailed information
about the model. The function "tfdate" retreives the information,
which is stored in cell arrays (both in SISO-cases and in MIMO-cases).
As an example, the following expressions retreive the polynomials containing
the coefficients of the numerator and denominator polynomials of the (SISO-)
transfer function Hd(s) created earlier:
MATLAB's response is
Here, Ts is the sampling
time which is zero for continuous-time models but non-zero for discrete-time
models. The response above is not
very descriptive since we do not see any numerical values for the numerator
and denominator coefficients. Suppose we want to retreive these coefficients
as vectors. This can be done as follows:
MATLAB's response is
Note: The use of inputno
and outputno above can be applied also for MIMO-systems. For example,
inputno=3 and outputno=2 corresponds to the transfer function from input no.
3 to output no. 2 in a MIMO-system. Note: The time-delay of
Hd(s) was ignored in the data retreival above. Continuous-time state-space models
The function
"ss" creates linear state-space models having the form dx/dt=Ax + Bu y=Cx + Du as illustrated in the
following example: Given the following state-space model: dx1/dt=x2 dx2/dt=-4*x1 -2*x2 + 2*u y=x1 This model is created by
To display ss1 we execute
MATLAB's response is
Above, a, b, c, and d,
which are default names, are used only for displaying purposes, and they do
in fact not exist in the workspace. To retreive the system parameters and
give them specific names, we can execute
MATLAB's response is
Discrete-time transfer functions
Discrete-time
(z-)transfer functions and state-space models are created in the same way as
for continuous-time models. As an example, to create
the transfer function Hdis(z)=(z+0.5)/(z2+1.5*z+2) with sampling time 0.4
(time-units), we execute
To display Hdis we
execute
MATLAB's response is
In digital signal
processing (DSP) it is common to use z-1 in the transfer functions
(and not explicitly z). The function "filt" supports this
DSP-format, and "filt" requires two vectors (one for the numerator
and one for the denominator) containing the coefficients for descending
orders of z-1. As an example, to create the transfer function Hdsp(z)=(1+0.5*z-1)/(1+1.5*z-1+2* z-2) with sampling time 0.4
(time-units), we execute
To display Hdsp we
execute
MATLAB's response is
Note: Hdsp(s) is not
equal to Hdis(s), although the coefficient vectors are equal! Discrete-time state-space models
Discrete-time state-space
models are assumed to have the form x(k+1)=Ad*x(k) + Bd*u(k) y(k)=Cd*x(k) + Dd*u(k) where Ad, Bd, Cd, and Dd
are matrices. x(k) is the state-variable, u(k) is the input variable, and y
is the output variable (these variables are generally vectors having one or
more elements). k is the discrete-time or time-index. Discrete-time
state-space models are created using the syntax ssdis=ss(Ad, Bd, Cd, Dd, Ts) where Ts is the
sampling-time. 2.2 Extracting subsystems
To create a new system
(say, sys1) which is a subsystem of an existing system (say, sys), we use the
syntax sys1=sys(outputno,inputno) As an example, let us
extract the transfer function Hsub from input 2 (u2) to the output (y) from
the MIMO-transfer function H defined earlier in this chapter:
MATLAB's response is
In fact, we can use
ordinary matrix-like accessing, like for example
"subsys=sys(1,2:4)" which selects inputs 2 to 4 and output 1, that is,
three subsystems altogether. 2.3 Setting and accessing LTI-properties
Several properties
can be attached to the LTI-objects (or LTI-models). The table belos gives an
overview over the properties. To get detailed information about legal values
of the properties, execute "help ltiprops". In the table below, NU,
NY, and NX refer to the number of inputs, outputs, and states (when
applicable). TF refers to transfer functions, SS to state-space models, and
ZPK to zero-pole-gain models.
To see the values of the
properties of a model, we use "get". Let us as an example look at
the properties of the MIMO-transfer function H defined earlier:
MATLAB's response is
We can access specific
LTI-properties. Let us as an example assign the value of the property Td of H
to the variable Tdvec:
MATLAB's response is
To set certain
LTI-properties, we use "set". Let us as an example set the property
InputName of H to {'Voltage';'Flow'}, which is a cell array:
These input names are now
used when displaying information about H:
MATLAB's response is
2.4 Conversion between transfer functions and
state-space models
From state-space model to transfer function
The function
"tf" converts state-space models to a corresponding transfer
function according to the formula(s)
In the following example
a continuous-time state-space model is converted to a corresponding
(continuous-time) transfer function:
MATLAB's response is
And in the following
example a discrete-time state-space model is converted to a
corresponding (discrete-time) transfer function (the sampling-time is 0.1):
MATLAB's response is
From transfer function to state-space model
"ss" converts a
transfer function to a state-space model. But this convertion is not unique,
since there are an infinite number of state-space models which have the same
transfer functions. Hence, "ss" chooses one such realization,
called the controller-canonical form. 2.5 Conversion between continuous-time and
discrete-time models
From discrete-time to continuous-time
A continuous-time
LTI-model can be discretized using the function "c2d", using the
syntax sysdis=c2d(syscont,Ts,method) Ts is the sampling time.
method is a string telling which discretizing method to be used, and the
following are available:
As an example, let us
discretize a continuous-time model (transfer function) which contains a
time-delay, assuming zero order hold:
MATLAB's response is
(The high order of the
denominator is due to the time-delay of the continuous-time system.) From discrete-time models to continuous-time models
Use the function "d2c". 2.6 Combining models (in series, parallel, and
feedback)
A model can be thought of
as a block with inputs and outputs and containing a transfer function or a
state-space model inside it. The functions "series",
"parallel", and "feedback" can be used to perform basic
block-diagram manipulation, as demonstrated in the following. These functions
apply to continuous-time models as well as discrete-time systems. As an example, assume
given the following transfer functions: Ha(s)=1/(s+2) First, let us create the
above two systems:
Let us connect Ha(s) and
Hb(s) in series, which means that Hser(s) becomes Ha(s)*Hb(s):
MATLAB's response is
Below, Ha(s) and Hb(s)
are connected in parallel, which means that Hpar(s)=Ha(s)+Hb(s).
MATLAB's response is
Finally, Ha(s) and Hb(s)
are connected in a feedback loop with Ha(s) being in the forward path
and Hb(s) being in the feedback path:
Above, feedbsign is the
sign of the feedback, and we set it to either -1 (in which case it can be
omitted) or 1. The transfer function Hfeedb(s) is given by Hfeedb(s)=Ha(s)/[1-feedbsign*Ha(s)*Hb(s)] MATLAB's response to the
expressions above is
3 Model analysis tools
As a starting point for
this chapter, let us define a continuous-time transfer function, H1c(s), and
a discrete-time transfer function, H1d(z) which will be the discretized
version of H1c(s). We take H1c(s) as H1c(s)=(s+0.5)/(s^2+2s+4) and we create H1c(s) with
MATLAB's response is
Then, let H1d(z) be the
zero-order hold discretized version of H1c(s):
MATLAB's response is
3.1 Simulation (time response)
We can simulate the step
response of one or several LTI-models with "step". The input
step is then a unit step (the step height is equal to one). Let us simulate
both H1c and H1d from t=0 to Tfinal=10 time units:
The resulting plot is not
shown here in this text. We can store the
simulated time response, y, and the time vector, t, by using return (or left)
arguments:
In this case, no plot
will be shown on the screen. Note that only single systems can be simulated
when you use left-hand arguments, as above, that is
"[y,t]=step(H1c,H1d,Tf);" will not work. "impulse"
simulates the impulse response:
The resulting plot is not
shown here in this text. "initial"
simulates the response from the initial state of an LTI-model. "lsim" is a
more general simulation function as it accepts any input signal, not just
step or impulse. To generate the input signal, we may use "gensig",
which produces a sine wave, or a square wave, or periodic pulses. Let us as
an example simulate H1c with a sinusoidal input, u, having period Tp=0.5,
final time Tf=10, and time step Tstep=0.01:
Above, "gensig"
was used to generate the input signal. Other signal generators are
"square", "sawtooth", and "stepfun". 3.2 Poles, eigenvalues, and zeros
We can calculate the
poles of an LTI-model with "pole", as in
MATLAB's response is
For a state-space model
we can analogously calculate the eigenvalues with "eig". "tzero"
calulates zeros (transmission zeros for MIMO-models). Here is an example
based on H1c(s):
MATLAB's response is
With
"[pol,zer]=pzmap" we can calculate and plot both the poles and the
zeros. Omitting the return-parameters just plots the poles and the zeros in
the complex plane:
The resulting plot is not
shown here in this text. Poles and eigenvalues
will always be real numbers and/or complex conjugate numbers.
"damp" calculates the relative damping, z , and the natural (or udamped) resonance
frequency, w 0 for such complex poles/eigenvalues. It is assumed that the poles
or eigenvalues are given as the roots of (the characteristic equation) s2 + 2z w 0s + w 02 Here is one example:
MATLAB's response is
3.3 Frequency response
The frequency response of
an LTI-model can be calculated with "bode". Here, "bode"
shows the Bode-plots of both H1c and H1d, in the same diagram:
The resulting plot is not
shown here in this text. Using the syntax [mag,phase,w]=bode(sys) the magnitude and the
phase (in degrees), and the frequency vector (in rad/s) used by
"bode", are returned. With "bode(sys,w)" the Bode-plot
will be drawn using our own frequency vector, w. The frequency response
can be drawn in a Nichols-diagram with "nichols", and in a
Nyquist-diagram with "nyquist" (examples follows in the next
section). Unfortunately, the
Bode-plots are often drawn with very few grid lines, and it may be very
difficult to read off precise values from it. To have more grids, you can do
as in the following example (in the example only the amplitude plot is drawn,
but similar expressions can be used for the phase plot):
3.4 Stability analysis based on frequency response
The function "margin"
calculates the stability margins and the cross-over frequencies for a
feedback system (continuous-time or discrete-time). In the following example
"margin" draws a Bode-plot and writes the stability margins and the
cross-over frequencies in the diagram. The loop transfer function of the
feedback system is Hloop(s)=(K/s)*e-Tdel*s where K=0.19, and
Tdel=4.17 (time-delay). Let us first create the LTI-object Hloop:
MATLAB's response is Gm = Here, Gm is gain margin
in dB, Pm is phase margin in degrees, Wcg is cross-over frequency where phase
is -180 degrees, and Wcp is cross-over frequency where gain is 1 = 0 dB. If you omit the left-hand
arguments of "margin", that is "margin(Hloop)", a Bode
plot of Hloop is presented, with the values of the four above mentioned
left-hand arguments of "margin" indicated in the plot. Stability analysis can
also be performed using a Nyquist-plot or a Nichols-plot. Let us try Nichols:
Above, ngrid adds a
Nichols grid to the plot. The resulting Nichols-plot is not shown in this
text. However, you will see that the interesting details around the critical
point (gain=1=0dB, and phase=- 180 degrees) do not appear clearly. Try
"zoom on" (executed at the command line) to see the details! 3.5 Pade-approximation
Pade-approximations are
rational transfer functions which approximates the transfer function of a
time-delay. Pade-approximations can be calculated with "pade". Let
us as an example develop a 5'th order approximation to a time-delay of 2
seconds (or time-units), and represent it as an LTI-model (or -object):
MATLAB's response is
Now, let us (for fun)
compare the step response of tfpade with the ideal step response of the
time-delay (of 2 sec.). To generate a time-delayed step, we will use
"stepfun".
The plot will not be
shown in this text, but it illustrates the approximation effectively. 4 Controller design tools
4.1 Root locus
Given a closed-loop
system with negative feedback and loop transfer function k*L0(s), where k is
a constant. [k*L0(s) is then the product of the individual transfer functions
in the loop.] The function "rlocus" calculates the poles of the
closed-loop system as a function of k. These poles are the root locus of the
closed-loop system. The poles are the roots of the characteristic equation
1+k*L0(s). Here is one example: First, let us create the transfer function
L0(s):
MATLAB's response is
Now, to plot the poles of
the closed-loop system which have the loop transfer function of k*L0(s), we
execute:
The resulting figure
showing the root locus is not shown here in this text. Note that
"sgrid" plots a grid of z and w 0-curves.
("zgrid" does the same in the discrete-time case.) MATLAB chooses
automatically the proper range of k for which the root locus is drawn.
However, we can use our own values of k by executing "rlocus(L0,k0)".
Of coourse, we must have created k0 (a vector) in advance. Once the root locus is
drawn (using "rlocus"), we can choose a pole location in the root
locus plot, and ask MATLAB to calculate the corresponding value of k. This is
done with "rlocfind". As an example, let us find (the approximate)
value of k which causes one of the closed-loop poles to be placed on the
imaginary axis (the system is then on the stability limit). (It can be shown
that k0=2 produces closed-loop poles at -2, ± i.). We execute:
"rlocfind" puts
a crosshair cursor on the root locus plot. Let us, with the cursor, select
the pole at +i (on the imaginary axis). As a result, MATLAB returns the both
the selected point and the corresponding value of k:
In fact, our selected
point (-0.0023 + 0.9912i) is not exactly on the root locus, but MATLAB finds
the value of k0 corresponding to the nearest point, and calculates the
corresponding poles (pole0). 4.2 Pole placement
Suppose the system to be controlled has the (continuous-time-) state-space model dx/dt=Ax+Bu. The system is to be controlled by state-feedback: u=-Gx where G is a matrix of gains. Then the closed-loop system has the state-space model dx/dt=Ax+B(-Gx)=(A-BG)x. The function "place" calculates G so that the eigenvalues of (A-BG) (the closed-loop system) are as specified. Here is one simple example: Assume the system to be controlled is given by dx/dt=u (which is an integrator, from u to x). The closed-loop eigenvalue shall be E=-4:
(The answer g1=4 is
correct despite the warning message.) "place" works
similarly for discrete-time systems: Suppose the system to be
controlled has the model x(k+1)=Adx(k)+Bdu(k) and the controller is u(k)=-Gdx(k) Then, the closed-loop
system has the state-space model x(k+1)=(Ad-BdGd)x(k) "place"
calculates Gd so that the eigenvalues of (Ad-BdGd)
(the closed-loop system) are as specified in the vector Ed. Thus,
we can execute "Gd=place(Ad,Bd,Ed)". "place" can be
used to calculate state-estimator gains, too: Suppose given a system
with state-space model dx/dt=Ax+Bu, y=Cx+Du and that the states x are
to be estimated by the estimator dxe/dt=Axe+Bu+Ke, ye=Cxe+Du where the residual e is
given by e=y-ye. It can be shown that the state estimation error z
is given by dz/dt=(A-KC)z. The eigenvalues of this error model are given by
E=eig(A-KC) º eig(A'-C'K'). (Here, A' means the
transpose of A.) Thus, we can calculate the estimator-gain by executing
"K1=place(A',C',Ee); K=K1' ". 4.3 Optimal control
The function
"dlqr" performs linear-quadratic regulator design for discrete-time
systems. This means that the controller gain G in the feedback u(k)=-Gx(k) for the controlled system
x(k+1)=Adx(k)+Bdu(k) is calculated so that the
cost function J=Sum {x'Qx + u'Ru + 2*x'Nu} is minimized. The syntax
is [K,S,E]=dlqr(Ad,Bd,Q,R,N) The matrix N is set to zero if omitted. Above, S is the solution of the Riccati equation, E contains the closed-loop eigenvalues which are the eigenvalues of (Ad-BdG). Here is one example:
MATLAB's response is
"lqr" works analogously for continuous-time systems. "dlqry" minimizes the cost function J=Sum {y'Qy + u'Ru}in which only the output vector y, and not the whole state-vector x, is weighted. "dlqy" works analogously for continuous-time systems. 4.4 Kalman estimator (or -filter)
"kalman"
designs both continuous-time and discrete-time Kalman-filters. The discrete-time
Kalman-filter is more commonly used. In this case, it is assumed that the
system is given by a discrete-time state-space model in the form of an
LTI-object. The model is x(k+1)=Ax(k) + Bu(k) + Gw(k) y(k)=Cx(k) + Du(k) + Hw(k) + v(k) where w is the process
noise, and v is the measurement noise. Note that the LTI-object must be
created by the expression sys=ss(A,[B,G],C,[D,H]) Q, R, and N are the
white-noise covariances as follows: E{ww'}=Q, E{vv'}=R, E{wv'}=N The Kalman-estimator is x(k+1|k)=Ax(k|k-1) + Bu(k) + L[y(k) - y(k|k)] (Eq. 1) x(k|k)=x(k|k-1) + M[y(k) - y(k|k)] (Eq. 3) where y(k|k)=Cx(k|k) + Du(k) In a "two-step"
Kalman-estimator the aposteriori estimate or the "measurement update
estimate" - which is the estimate we usually need - is given by Eq. 3,
while the apriori estimate or the "time update" estimate is given
by Eq. 2. In a "one-step" Kalman-estimator Eq. 1 gives the estimate. The syntax of
"kalman" is [Kest,L,P,M,Z]=kalman(sys,Q,R,N) Here, N may be omitted if
it is zero. Kest is the Kalman estimator in the form of an LTI-object having
[u;yv] as inputs (yv is the measurements), and [yest,xest] as outputs. M is
the steady-state innovation gain or estimator gain for a "two-step"
Kalman-estimator. L is the steady-state innovation gain or estimator gain for
a "one-step" Kalman-estimator. P, and Z are the steady-state error
covariances: P=E{[x - x(k|k-1)][x - x(k|k-1)]'} (Covariance of apriori estimate
error) Z=E{[x - x(k|k)][x - x(k|k)]'} (Covariance of aposteriori estimate
error) Here is one example:
MATLAB's response is
5 To get an overview over Control System Toolbox
The result may be like
this: Control System Toolbox Version 10.3 (R2017b) 24-Jul-2017 Table of Contents ----------------- Modeling - Create and combine models of linear
systems Analysis - Analyze systems in the time or
frequency domain Design - Design and tune control systems Plotting - Create and customize response plots General - Preferences, linear algebra
utilities Apps - Control System Toolbox apps Examples - Control System Toolbox examples For example, by clicking
on the hypertext Analysis, you will get a list of functions for analysis of
dynamic systems, but details are not shown here. |
