Linear Systems and Control – EXPERIMENT 1

ECE356S – Linear Systems and Control –
EXPERIMENT 1
Introduction to MATLAB and Simulink for Systems Analysis
1 Purpose
The purpose of this experiment is to introduce you to MATLAB and Simulink as tools for systems analysis.
MATLAB is used primarily for numerical computations, while Simulink provides a block diagram style
representation of systems for analysis, simulation, and design.
2 Preparation
Before you go to your lab session, read through carefully the description of this laboratory provided in
the following pages. Useful terms and commands are highlighted in bold font, while laboratory work is
written in bold italics. Identify the parts you have to do in the lab.
3 Linear Algebra
MATLAB is heavily based on linear algebra for doing its computations. In MATLAB, matrices are entered
row by row. The end of a row is indicated by a “;” or a carriage return. For example,
A = [1 2 3; 4 5 6; 7 8 9]
produces the matrix
A =


1 2 3
4 5 6
7 8 9


Individual rows and columns can be extracted from A. For example, A(:, 1) and A(1, :) give the first
column and row of A, respectively.
I. Inverse of a Matrix and Solution of Linear Equations:
MATLAB provides a command to compute the inverse of an n × n matrix A, namely inv(A).
However, if the objective is to solve a linear system of equations Ax = b, there is a better alternative
than using x = inv(A) ∗ b, namely x = A\b. The MATLAB operation “\” is called mldivide or
matrix left divide. This is illustrated in the following calculation.
(a) Save the following code to a file and name it inv matrix.m. This is a script file which you can
1
run in MATLAB.
n = 500;
Q = orth(randn(n,n));
d = logspace(0,-10,n);
A = Q*diag(d)*Q’;
x = randn(n,1);
b = A*x;
tic, y = inv(A)*b; toc
err = norm(y-x)
res = norm(A*y-b)
pause
tic, z = A\b; toc
err = norm(z-x)
res = norm(A*z-b)
Consult the MATLAB documentation to understand what each line means. Then
run “inv matrix” in MATLAB and note the outputs.
(b) Similarly, there is an operation called mrdivide or matrix right divide, A/B, which has the
effect of AB−1 when B is invertible. Let A = randn(4, 4). Compute the inverse of A using
inv, mldivide, and mrdivide. Show your results to your TA. Explain how you would
compare the accuracy of the three methods. For square matrices with low dimension, all
three calculations should give comparable accuracy.
II. Eigenvalues and Eigenvectors
Of great importance in the analysis of linear systems are the notions of eigenvalues and eigenvectors.
A scalar λ is an eigenvalue and a non-zero vector x is an eigenvector of a linear transformation A if
Ax = λx. If A is a matrix, then eigenvalues and eigenvectors are only defined if A is square.
(a) Let A be given by
A =


7 2 −3
4 6 −4
5 2 −1


We will use this A from item (a) to (e) in this part. Use the MATLAB command [V, D] =
eig(A) to determine the eigenvalues and eigenvectors of A.
(b) The command eig(A) gives eigenvectors which are normalized to have norm 1. If you compute the eigenvectors by hand, you are more likely to come up with eigenvectors which are not
normalized. For the matrix in II(a), you can use the command eig(A,’nobalance’)
to produce more “readily recognizable” eigenvectors. Recall that eigenvectors are determined up to a scalar multiple. From the results of eig(A,’nobalance’), determine 3
eigenvectors of A whose entries are all integers.
(c) For manual computation of eigenvalues, usually you determine the characteristic polynomial of
A given by det(sI − A). Then you find the roots of the characteristic polynomial, i.e., find
solutions of the characteristic equation det(sI − A) = 0. Determine the characteristic
polynomial of A using the command poly(A), and determine the eigenvalues by
applying the command roots to the resulting polynomial. Compare the answer with
those from (a).
(d) The eigenvalue-eigenvector equations can be written as one matrix equation
AV = V D,
2
with D a diagonal matrix consisting of the eigenvalues of A. From the results of IV(a),
determine norm(AV-VD). Show more significant digits in MATLAB using the command format long. From this determine the exact values of the eigenvalues and
eigenvectors (up to a scalar multiple). Now verify that AV − V D = 0.
(e) Recall that if the eigenvalues of A are distinct, the eigenvectors are linearly independent. In
that case, the V matrix computed using eig is invertible and that V
−1AV = D. This process is
called diagonalizing A. Check that the matrix in II(a) can be diagonalized.
(f) When A has repeated eigenvalues, it may not be possible to diagonalize A. Suppose A has λ
as a repeated eigenvalue, then det(λI − A) = 0 and the number of times λ repeats as roots is
called its algebraic multiplicity. Thus the following matrix
A =

1 1
0 1
has 1 as its only eigenvalue with algebraic multiplicity 2. Determine by hand calculation
the eigenvector corresponding to the eigenvalue 1, and verify there is only one
independent eigenvector. Try using eig on this A. Does the resulting [V,D] satisfy
AV = V D? Is V invertible?
(g) When A cannot be diagonalized, one can transform it to a form called the Jordan (canonical)
form. The MATLAB command is jordan. For the matrix
A =


0 4 3
0 20 16
0 −25 −20


find its Jordan form.
4 Ordinary Differential Equations and Transfer Functions
Control systems studied in this course are modelled by in the time domain by state equation of the form
x˙ = Ax + Bu (1)
y = Cx + Du, (2)
or by input/output models of the form
y
(n)
(t) + a1y
(n−1)(t) + a2y
(n−2)(t) + · · · + any(t) = b0u
(m)
(t) + · · · + bmu(t), with m ≤ n
In the state equation, the transfer matrix from u to y is given by
G(s) = C(sI − A)
−1B + D.
In this course, we focus on single-input single-output (SISO) systems, i.e., u and y are scalar-valued. In this
case, the transfer function G(s) is a scalar-valued proper rational function. In the case of an input/output
model, the transfer function G(s) is given by
G(s) = b0s
m + b1s
m−1 + · · · + bm
s
n + a1s
n−1 + · · · + an
MATLAB provides tools to solve such linear systems. We illustrate some of these tools using a second
order system.
3
Consider a second order system described by the differential equation
y¨ + 2 ˙y + 4y = 4u
The transfer function from the input u to the output y for this system is given by
G(s) = 4
s
2 + 2s + 4
(3)
To model this as a state equation, let
x(t) = 
y(t)
y˙(t)

Then the state x satisfies the differential equation
x˙ =

0 1
−4 −2

x +

0
4

u (4)
y =

1 0
x (5)
MATLAB provides a data object which conveniently describes the state space system (1)-(2). It is created
by the command sys=ss(A,B,C,D).
(a) Create the sys object corresponding to the second order system described by (4)-(5).
(b) One can study the response of the second order system due to particular inputs such as a step (to get
the step response), or the response due to nonzero initial conditions with no input, or generally with
nonzero initial conditions and inputs. Here, determine the step response using the command
[Y,T,X]=step(sys), where sys is the object you created using the ss command. Use
plot(T,X) to plot the trajectories of both states.
(c) To determine the response due only to initial conditions, it would be convenient to create a second
sys object using sys init=ss(A,0,C,0) (0 is a matrix of appropriate dimension). Use the command
[Y,T,X]=initial(sys init,x0) to determine the response to initial conditions when x0 =

0
1

. Again plot the state responses.
(d) In general, the system may have both nonzero initial conditions as well as inputs. For example, the
commands
t=0:0.01:20;
u=sin(t);
define an input vector u whose components are given by the values of sin(t) at various times specified
by the time vector t. Use [Y,t,X]=lsim(sys,u,t,x0) to produce the response of the second
order system due to the initial condition and sinusoidal input defined above. Plot the
state trajectories.
(e) If you are given a state space description, the command ss2tf generates the numerator and denominator polynomials of the transfer function associated with the state space system. Use the ss2tf
command, followed by the tf command to produce the transfer function for the systems
described by (4)-(5). Verify that it is the same as (3).
4