More free stuff from TechTeach:


Tutorial for

Control System Toolbox for MATLAB

by

Finn Haugen


October 11, 2003

This publication can be downloaded and copied freely, but reference to the source is appreciated.


Contents:

1 Introduction
2 Creating and handling LTI-models (Linear Time-Invariant)
2.1 Creating models
2.2 Extracting subsystems
2.3 Setting and accessing LTI-properties
2.4 Conversion between state-space models and transfer functions
2.5 Conversion between continuous-time and discrete-time models
2.6 Combining models (series, parallel, and feedback)
3 Model analysis tools
3.1 Simulation (time response)
3.2 Poles, eigenvalues, and zeros
3.3 Frequency response
3.4 Stability analysis based on frequency response
3.5 Pade-approximations
3.6 The LTI-viewer
4 Controller design tools
4.1 Root locus
4.2 Pole placement
4.3 Optimal control
4.4 Kalman estimator (or filter)
5 Overview over functions in Control System Toolbox


1 Introduction

This 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:

  • transfer functions
  • state-space models

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)

K0=2; T0=4;
num0=[K0,0]; den0=[T0,1];
H0=tf(num0,den0);

To display H0 we execute

H0

MATLAB's response is

Transfer function:
2 s
-------
4 s + 1

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):

Tdelay=1;
Hd=tf(num0,den0,'InputDelay',Tdelay)

To display Hd we execute

Hd

MATLAB's response is

Transfer function:

2 s
-------
4 s + 1

Input delay: 1

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]

K1=2; T1=4; num1=[K1,0]; den1=[T1,1];
K2=3; T2=5; num2=[K2]; den2=[T2,1];
Tdelay1=1;
Tdelay2=2;
H1=tf(num1,den1,'InputDelay',Tdelay1);
H2=tf(num2,den2,'InputDelay',Tdelay2);
H=[H1,H2];

To display H we execute

H

MATLAB's response is

Transfer function from input 1 to output:

                     2 s
exp(-1*s) * -------
                  4 s + 1

                     3
exp(-1*s) * -------
                  5 s + 1

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)
H2(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:

[numHd,denHd,Ts]=tfdata(Hd)

MATLAB's response is

numHd=
[1x2 double]
denHd=
[1x2 double]
Ts=
0

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:

inputno=1;
outputno=1;
numpolHd=numHd{outputno,inputno}
denpolHd=denHd{outputno,inputno}

MATLAB's response is

numpolHd=
2 0
denpolHd=
4 1

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

A=[0,1;-4,-2]; B=[0;2]; C=[1,0]; D=[0];
ss1=ss(A,B,C,D);

To display ss1 we execute

ss1

MATLAB's response is


a=

x1 x2

x1 0 1.00000

x2 -4.00000 -2.00000


b=

u1

x1 0

x2 2.00000


c=

x1 x2

y1 1.00000 0


d=

u1

y1 0

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

[A1,B1,C1,D1,Ts]=ssdata(ss1)

MATLAB's response is


A1=

0 1

-4 -2


B1=

0

2


C1=

1 0


D1=

0


Ts=

0

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

numdis=[1,0.5]; dendis=[1,1.5,2]; Tsdis=0.4;
Hdis=tf(numdis,dendis,Tsdis);

To display Hdis we execute

Hdis

MATLAB's response is

Transfer function:
z + 0.5
---------------
z^2 + 1.5 z + 2
Sampling time: 0.4

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

numdsp=[1,0.5];
dendsp=[1,1.5,2];
Tsdsp=0.4;
Hdsp=filt(numdsp,dendsp,Tsdsp);

To display Hdsp we execute

Hdsp

MATLAB's response is

 

Transfer function:
1 + 0.5 z^-1
---------------------
1 + 1.5 z^-1 + 2 z^-2
Sampling time: 0.4

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:

outputno=1; inputno=2;
Hsub=H(outputno,inputno)

MATLAB's response is

Transfer function from input 2 to output:

3
-------
5 s + 1

Input delay: 2

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.

Name

Description

Applies to LTI-object of type

Ts

Sample Time

All

Td

Input delays

All (only continuous systems)

Inputname

Input names

All

Outputname

Output names

All

Notes

Notes on model history

All

Userdata

Additional data

All

Num

Numerator

TF

Den

Denominators

TF

Variable

Transfer function variable (e.g. 's' or 'z')

TF, ZPK

Z

Zeros

ZPK

P

Poles

ZPK

K

Gains

ZPK

A

Matrix A

SS

B

Matrix B

SS

C

Matrix C

SS

D

Matrix D

SS

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:

get(H)

MATLAB's response is

num={1x2 cell}

den={1x2 cell}

Variable='s'

Ts=0

Td=[1 2]

InputName={2x1 cell}

OutputName={''}

Notes={}

UserData=[]

We can access specific LTI-properties. Let us as an example assign the value of the property Td of H to the variable Tdvec:

Tdvec=get(H,'Td')

MATLAB's response is

Tdvec=

1 2

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:

set(H,'InputName',{'Voltage';'Flow'})

These input names are now used when displaying information about H:

H

MATLAB's response is

Transfer function from input "Voltage" to output:

2 s
-------
4 s + 1

Transfer function from input "Flow" to output:

3
-------
5 s + 1

Input delays (listed by channel): 1 2


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)

H(s)=C(sI-A)-1B + D (continuous-time case)

Hd(z)=Cd(zI-Ad)-1Bd + Dd (discrete-time case)

In the following example a continuous-time state-space model is converted to a corresponding (continuous-time) transfer function:

A= -1; B=2; C=3; D=4;
sscont1=ss(A,B,C,D);
Hcont1=tf(sscont1)

MATLAB's response is

Transfer function:

4 s + 10
--------
s + 1

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):

Ad=0.5; Bd=0.2; Cd=0.3; Dd=0.4; Ts=0.1;
ssdis1=ss(Ad,Bd,Cd,Dd,Ts);
Hdis1=tf(ssdis1)

MATLAB's response is

Transfer function:

0.4 z - 0.14
------------
z - 0.5

Sampling time: 0.1

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:

'zoh' for zero order hold.
'foh' for first order hold.
'tustin' for Tustin's method of bilinear transformation.
'prewarp' for Tustin's method with frequency prewarping.
'matched' for the matched pole-zero method.

As an example, let us discretize a continuous-time model (transfer function) which contains a time-delay, assuming zero order hold:

K=2; T=4; num=[K,0]; den=[T,1]; Tdelay=1;
Hcont4=tf(num,den,'Td',Tdelay);
Ts=0.2;
Hdis4=c2d(Hcont4,Ts,'zoh')

MATLAB's response is


Transfer function:

0.5 z - 0.5
----------------
z^6 - 0.9512 z^5

Sampling time: 0.2

(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)

Hb(s)=3/(s+4)

First, let us create the above two systems:

k1=1; k2=2; k3=3; k4=4;
Ha=tf(k1,[1,k2]);
Hb=tf(k3,[1,k4]);

Let us connect Ha(s) and Hb(s) in series, which means that Hser(s) becomes Ha(s)*Hb(s):

Hser=series(Ha,Hb)

MATLAB's response is

Transfer function:
3
-------------
s^2 + 6 s + 8

Below, Ha(s) and Hb(s) are connected in parallel, which means that Hpar(s)=Ha(s)+Hb(s).

Hpar=parallel(Ha,Hb)

MATLAB's response is

Transfer function:
4 s + 10
-------------
s^2 + 6 s + 8

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:

feedbsign=-1; Hfeedb=feedback(Ha,Hb,feedbsign)

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

Transfer function:
s + 4
--------------
s^2 + 6 s + 11

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

n1=[1,0.5]; d1=[1,2,4];
H1c=tf(n1,d1)

MATLAB's response is


Transfer function:

s + 0.5
-------------
s^2 + 2 s + 4

Then, let H1d(z) be the zero-order hold discretized version of H1c(s):

Ts=0.2;
H1d=c2d(H1c,Ts,'zoh')

MATLAB's response is


Transfer function:

0.1692 z - 0.1529

---------------------

z^2 - 1.54 z + 0.6703

Sampling time: 0.2


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:

Tf=10; step(H1c,H1d,Tf)

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:

Tf=10; [y,t]=step(H1d,Tf);

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:

Tf=10; impulse(H1c,Tf)

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:

Tp=0.5; Tf=10; Tstep=0.01; [u,t]=gensig('sin',Tp,Tf,Tstep);
lsim(H1c,u,t)

 


The resulting plot is not shown here in this text.

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

p=pole(H1c)

MATLAB's response is

p=

-1.0000 + 1.7321i

-1.0000 - 1.7321i

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):

z=tzero(H1c)

MATLAB's response is

z=

-0.5000

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:

pzmap(H1c)

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:

damp(H1c)

MATLAB's response is

Eigenvalue Damping Freq. (rad/s)

-1.00e+000 + 1.73e+000i 5.00e-001 2.00e+000

-1.00e+000 - 1.73e+000i 5.00e-001 2.00e+000


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:

bode(H1c,H1d)

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):

n1=[1]; d1=[1,2,4]; H1=tf(n1,d1);
[mag,phase,w]=bode(H1);
semilogx(w,mag(:));
set(gca,'xtick',[.1,.2,.5,1,2,5,10]);
set(gca,'ytick',[0:0.05:0.4]);
grid


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:

K=0.19; Tdel=4.17;
Hloop=tf(K,[1,0],'InputDelay',Tdel);
[Gm,Pm,Wcg,Wcp]=margin(Hloop)

MATLAB's response is

Gm =
1.9826
Pm =
44.6046
Wcg =
0.3767
Wcp =
0.1900

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:

nichols(H1c),ngrid

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):

Tdelay=2; N=5;
[numpade,denpade]=pade(Tdelay,N);
tfpade=tf(numpade,denpade)

MATLAB's response is

Transfer function:

-s^5 + 15 s^4 - 105 s^3 + 420 s^2 - 945 s + 945
--------------------------------------------------
s^5 + 15 s^4 + 105 s^3 + 420 s^2 + 945 s + 945

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".

tfinal=5; [ypade,time]=step(tfpade,tfinal);
T0=2; yideal=stepfun(time,T0);
plot(time,[ypade,yideal])

The plot will not be shown in this text, but it illustrates the approximation effectively.


3.6 The LTI-viewer

LTI-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

Ltiview

and try it.


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):

L0=tf(1,[1 2 1 0])

MATLAB's response is

Transfer function:

1
--------------
s^3 + 2 s^2 + s

 

Now, to plot the poles of the closed-loop system which have the loop transfer function of k*L0(s), we execute:

rlocus(L0),sgrid

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:

[k0,pole0]=rlocfind(L0)

"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:

selected_point=

-0.0023 + 0.9912i

k0=

1.9606

pole0=

-1.9921

-0.0040 + 0.9921i

-0.0040 - 0.9921i

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:

E=-4;
A=0;B=1;
g1=place(A,B,E)


MATLAB's response is

place: ndigits=15

Warning: Pole locations are more than 10% in error.

g1=

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:

Ad=0.85;Bd=0.15;
Q=1;R=0.2;
[G,S,E]=dlqr(Ad,Bd,Q,R)

MATLAB's response is

G=

1.1797

S=

2.3370

E=

0.6731

"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)
=Ax(k|k) + Bu(k) (Eq. 2)

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:

A=0.85;B=0.15;C=1;D=1;Ts=0.1;
G=1;H=0;
s1=ss(A,[B,G],C,[D,H],Ts);
Q=1; R=0.2;
[Kest,L,P,M,Z]=kalman(s1,Q,R)

MATLAB's response is

a=

x1_e

x1_e 0.12853

b=

u1 u2

x1_e -0.57147 0.72147

c=

x1_e

y1_e 0.15121

x1_e 0.15121

d=

u1 u2

y1_e 0.15121 0.84879

x1_e -0.84879 0.84879

Sampling time: 0.1

Discrete-time system.

L=

0.7215

P=

1.1226

M=

0.8488

Z=

0.1698


5 Overview over functions in Control System Toolbox

Version 4.0.1 07-Mar-1997

Creation of LTI models

ss - Create a state-space model.
zpk - Create a zero/pole/gain model.
tf - Create a transfer function model.
dss - Specify a descriptor state-space model.
filt - Specify a digital filter.
set - Set/modify properties of LTI models.
ltiprops - Detailed help for available LTI properties.

Data extraction

ssdata - Extract state-space matrices.
zpkdata - Extract zero/pole/gain data.
tfdata - Extract numerator(s) and denominator(s).
dssdata - Descriptor version of SSDATA.
get - Access values of LTI model properties.

Model characteristics

class - Model type ('ss', 'zpk', or 'tf').
size - Input/output dimensions.
isempty - True for empty LTI models.
isct - True for continuous-time models.
isdt - True for discrete-time models.
isproper - True for proper LTI models.
issiso - True for single-input/single-output systems.
isa - Test if LTI model is of given type.

Conversions

ss - Conversion to state space.
zpk - Conversion to zero/pole/gain.
tf - Conversion to transfer function.
c2d - Continuous to discrete conversion.
d2c - Discrete to continuous conversion.
d2d - Resample discrete system or add input delay(s).

Overloaded arithmetic operations

+ and - - Add and subtract LTI systems (parallel connection).
* - Multiplication of LTI systems (series connection).
\ - Left divide -- sys1\sys2 means inv(sys1)*sys2.
/ - Right divide -- sys1/sys2 means sys1*inv(sys2).
' - Pertransposition.
.' - Transposition of input/output map.
[..] - Horizontal/vertical concatenation of LTI systems.
inv - Inverse of an LTI system.

Model dynamics

pole, eig - System poles.
tzero - System transmission zeros.
pzmap - Pole-zero map.
dcgain - D.C. (low frequency) gain.
norm - Norms of LTI systems.
covar - Covariance of response to white noise.
damp - Natural frequency and damping of system poles.
esort - Sort continuous poles by real part.
dsort - Sort discrete poles by magnitude.
pade - Pade approximation of time delays.

State-space models

rss,drss - Random stable state-space models.
ss2ss - State coordinate transformation.
canon - State-space canonical forms.
ctrb, obsv - Controllability and observability matrices.
gram - Controllability and observability gramians.
ssbal - Diagonal balancing of state-space realizations.
balreal - Gramian-based input/output balancing.
modred - Model state reduction.
minreal - Minimal realization and pole/zero cancellation.
augstate - Augment output by appending states.

Time response

step - Step response.
impulse - Impulse response.
initial - Response of state-space system with given initial state.
lsim - Response to arbitrary inputs.
ltiview - Response analysis GUI.
gensig - Generate input signal for LSIM.
stepfun - Generate unit-step input.

Frequency response

bode - Bode plot of the frequency response.
sigma - Singular value frequency plot.
nyquist - Nyquist plot.
nichols - Nichols chart.
ltiview - Response analysis GUI.
evalfr - Evaluate frequency response at given frequency.
freqresp - Frequency response over a frequency grid.
margin - Gain and phase margins

System interconnections

append - Group LTI systems by appending inputs and outputs.
parallel - Generalized parallel connection (see also overloaded +).
series - Generalized series connection (see also overloaded *).
feedback - Feedback connection of two systems.
star - Redheffer star product (LFT interconnections).
connect - Derive state-space model from block diagram description.

Classical design tools

rlocus - Evans root locus.
rlocfind - Interactive root locus gain determination.
acker - SISO pole placement.
place - MIMO pole placement.
estim - Form estimator given estimator gain.
reg - Form regulator given state-feedback and estimator gains.

LQG design tools

lqr,dlqr - Linear-quadratic (LQ) state-feedback regulator.
lqry - LQ regulator with output weighting.
lqrd - Discrete LQ regulator for continuous plant.
kalman - Kalman estimator.
kalmd - Discrete Kalman estimator for continuous plant.
lqgreg - Form LQG regulator given LQ gain and Kalman estimator.

Matrix equation solvers

lyap - Solve continuous Lyapunov equations.
dlyap - Solve discrete Lyapunov equations.
care - Solve continuous algebraic Riccati equations.
dare - Solve discrete algebraic Riccati equations.

Demonstrations

ctrldemo - Introduction to the Control System Toolbox.
jetdemo - Classical design of jet transport yaw damper.
diskdemo - Digital design of hard-disk-drive controller.
milldemo - SISO and MIMO LQG control of steel rolling mill.
kalmdemo - Kalman filter design and simulation.