Tutorial forControl System Toolbox for MATLABby
This publication can be downloaded and copied freely, but reference to the source is appreciated. Contents:
1 Introduction 1 IntroductionThis text gives an easy guide to Control System Toolbox version 5 which assumes the installation of MATLAB version 6. 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. It is assumed that you have the basic knowledge about dynamic systems and control theory. It is also assumed that you have basic skills in MATLAB. For Norwegian users I may suggest Lær MATLAB trinn for trinn (http://techteach.no/publications/matlab). And if you want to simulate your (control) systems, you may consider using SIMULINK which is based on MATLAB. Tip for Norwegian users: Lær SIMULINK trinn for trinn (http://techteach.no/publications/simulink). 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 modelsWe will start with continuous-time models, and then take discrete-time models. Continuous-time transfer functionsThe 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
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 modelsThe 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 functionsDiscrete-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)/(z^{2}+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 modelsDiscrete-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 subsystemsTo 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-propertiesSeveral 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 modelsFrom state-space model to transfer functionThe 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 modelsFrom discrete-time to continuous-timeA 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: 'zoh' for zero order hold. 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 modelsUse 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 toolsAs 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 zerosWe 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) s^{2 }+ 2z w _{0}s + w _{0}^{2} Here is one example:
MATLAB's response is
3.3 Frequency responseThe 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 responseThe 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-approximationPade-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. 3.6 The LTI-viewerLTI-viewer makes it easy to perform simulations, frequency response plots etc. for one or several of the LTI-models which exist in the working directory. The LTI-viewer is easy to use. Just execute
and try it. 4 Controller design tools4.1 Root locusGiven 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 placementSuppose 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)=A_{d}x(k)+B_{d}u(k) and the controller is u(k)=-G_{d}x(k) Then, the closed-loop system has the state-space model x(k+1)=(A_{d}-B_{d}G_{d})x(k) "place" calculates G_{d} so that the eigenvalues of (A_{d}-B_{d}G_{d}) (the closed-loop system) are as specified in the vector E_{d}. 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 dx_{e}/dt=Ax_{e}+Bu+Ke, y_{e}=Cx_{e}+Du where the residual e is given by e=y-y_{e}. 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 controlThe 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)=A_{d}x(k)+B_{d}u(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(A_{d},B_{d},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 (A_{d}-B_{d}G). 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 Overview over functions in Control System ToolboxVersion 4.0.1 07-Mar-1997 Creation of LTI models ss - Create a state-space model. Data extraction ssdata - Extract state-space matrices. Model characteristics class - Model type ('ss', 'zpk', or 'tf'). Conversions ss - Conversion to state space. Overloaded arithmetic operations + and - - Add and subtract LTI systems (parallel connection). Model dynamicspole, eig - System poles. State-space models rss,drss - Random stable state-space models. Time response step - Step response. Frequency responsebode - Bode plot of the frequency response. System interconnections append - Group LTI systems by appending inputs and outputs. Classical design tools rlocus - Evans root locus. LQG design tools lqr,dlqr - Linear-quadratic (LQ) state-feedback regulator. Matrix equation solvers lyap - Solve continuous Lyapunov equations. Demonstrations ctrldemo - Introduction to the Control System Toolbox. |