Tutorial for

Control System Toolbox for MATLAB

by

Finn Haugen


February 23, 2018

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
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 Control System Toolbox


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:

  • 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

Alternatively, you can define the symbol ‘s’ as the Laplace transform variable, and then write transfer functions as you would do in a document:

s=tf(‘s’);%sufficient to execute once in each MATLAB session
Hd1=2*s/(4*s+1)*exp(-1*s)

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]

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.


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 To get an overview over Control System Toolbox

help control

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.