Linear systems

3. LINEAR SYSTEMS

Although a great deal of effort has gone into finding general theories for non-linear systems, such a
theory does not yet exist and even theories applied to classes of non-linear systems can be very
difficult to use. For the very special sub-set of non-linear systems called linear systems, we do
have a well-developed theory. It is for this reason and because many systems have regions of
operation which are rather linear that we develop our control design understanding for linear
systems. We do not discard the non-linear modelling, as it is useful for accurate simulation of the
design that comes out of a simplified, linear model. Simulation has become a very accessible tool
for control engineering with the advent of fast, cheap computing systems. In addition to the
benefits for simulation, the non-linear tools that are available can to some extent be used to check
the outcome of engineering design.

3.1. Linearisation
The linearisation of eq(3.1) is just a vector version of the linearisation of a function of many
variables. If we have some operating condition, x 0 , u 0 then the linearisation of eq(3.1) is

f f
x - x0 =
& & (x - x 0 ) + (u - u 0 ) (3.1a)
x x , u u x , u
0 0 0 0


g g
y - y0 = (x - x 0 ) + (u - u 0 ) (3.1b)
x x ,u u x ,u
0 0 0 0


The partial derivative of a vector with respect to a vector is the Jacobian and is evaluated as
follows,

f1 f1
L L
x1 x 2
f
f 2 x
O

= (3.2)
x 1
M O

M f n

x n
We define
f
A( t ) = is the [n×n] state transition matrix
x x
0 ,u0

f
B( t ) = is the [n×m] input matrix
u x , u
0 0

g
C( t ) = is the [p×n] output matrix
x x , u
0 0

g
D( t ) = is the[p×m] throughput matrix
u x , u
0 0



Systems & Simulation ­ ENEL3SS 16
SS1_3.DOC
and get

x = A( t ) x + B( t ) u
&
y = C( t ) x + D( t ) u

Often the operating condition is a "steady state" (i.e. x 0 = 0 ) and we rename the state variables
&
and inputs and outputs to get rid of offsets (which must not be forgotten in practical applications).
We then get the familiar linear state space equations,

x = A( t ) x + B( t ) u
&
(3.3)
y = C( t ) x + D( t ) u

Eq(3.3) is linear in the sense that the variables, x, u and y appear only as linear combinations of
themselves and parameters. There are no x1x2 or exp(x1) terms. (Bilinear systems allow terms in
uixj).



x
& = A x + B u




y = C x + D u



In still more specialised applications, the Jacobians do not vary with time (i.e. are constant) and we
get the time invariant linear model equations,

x = Ax + B u
&
(3.4)
y = Cx + D u


Examples
The magnetic levitation control problem is of considerable practical engineering and academic
interest and has been investigated by a number of authors using different approaches (for example,
Yang and Tateishi, 2001; Lin and Tho, 1998, Hurley and Wölfe, 1997).

Section 2 presents the model and data used in the study. Section 3 presents an overview of the
LTIE method, assuming some familiarity with the QFT method for single-input, single-output
systems. Section 4 presents the design details and measured results for the tutorial example.


MAGNETIC LEVITATOR MODEL
Assuming a linear electrical circuit, the dynamic equations of the magnetic levitator illustrated in
Fig. 1 are given by,




Systems & Simulation ­ ENEL3SS 17
SS1_3.DOC
v - g + f (i, x) / m
d
x = v (E.1)
dt i (V (t ) - iR ) / L

where, v is the velocity, x is the position, i is the coil current, V is the applied coil voltage, R is the
coil resistance, L is the coil inductance, m is the mass of the object to be levitated, and f(i,x) is the
levitating force developed as a function of current and distance.

ELECTRO-
MAGNET


SENSOR

-x
EARTH




Fig. 1 ­ Magnetic Levitator



The electromagnet shown in Fig. 1 was designed and built as part of an undergraduate project. Our
set-up levitates a hollow steel moon (mass, mm=0.236 kg, diameter, m=150mm) or earth (mass,
me=0.333kg, diameter, e=180mm). The coil inductance is L[80, 100] mH and R[3, 3.5] ,
depending on operating conditions.

The sensor shown in Fig. 1 is an inductive proximity sensor that was built by our workshop. It
consists of a spiral inductor etched on a printed circuit board that is excited by a tuned circuit at
around 1 MHz. The proximity of the steel ball detunes the circuit, reducing the sensed voltage. The
distance is evaluated from calibration data via a look-up table. (An undergraduate project has built
and commissioned an optical sensor based on a linear array of charge-coupled devices.)

The magnetic force is often approximated by,
f (i, x) = k (i x )2 (E.2)
The measured force to current and position is shown in Fig. 2, and k=0.0017 roughly fits the
measured data shown. With the QFT design method, there is freedom to use mix of theoretical and
empirical relationships for the non-linear plant model.

If the magnetic force satisfies eq(2) or is available in empirical form, a simple non-linear approach
to control the levitation height, x, would be to use say a linear controller operating on the error to
generate a force set-point as its output (eq(1) is linear if the magnetic force can be controlled). The
required current is calculated from the height measurement and force demand. Finally a current
controller is used to achieve the required current ­ see Fig. 3. The scheme is a cascaded (inner and
outer) controller with ratio control on the inner loop and it is appealing and works. Students should
be alert to the possibility of using heuristic schemes. There are many successful applications gain
scheduling (ratio control) in the process industry when the dynamics are process dependent in a
benign way. Examples include mixing vessels where the lag depends on the flow-rate and systems
Systems & Simulation ­ ENEL3SS 18
SS1_3.DOC
with variable but measurable transport lag). This is not the main subject of this paper and will not
be pursued further

2 2
Local linearisation of eq(1) and eq(2), around positive i0 and negative x0, with k i0 x 0 = mg for
steady state, gives some insight into the behaviour of the system and has often been used in design
studies,

v 0 - 2 g x 0 2 g i 0 v 0
d
x = 1 0 0 x + 0 V (t ) (E.3)
dt
i 0 0 R L i 1 / L



This system has two real eigenvalues with the same magnitude located on either side of the
imaginary axis and a third associated with the electrical circuit,
1,2,3 = ± - 2 g x 0 , - R / L (E.4)
1,2 become very fast as the distance to the magnet becomes small. The gain from current to
magnetic force becomes small with increasing x0.

An obvious problem with using a single local linearisation for control design is that it is not valid
for tracking or if the levitated object is perturbed too far away from the operating point.

FORCE VS DISTANCE
18

16

14

12 I=[3.4:-0.2:2.4]
FORCE [N]




10

8

6

4 mearth g
mmoon g
2

0
-75 -70 -65 -60 -55 -50 -45 -40 -35 -30
DISTANCE [MM]

Fig. 2 ­ Force to current and distance
Controller
PID sqrt(u/k) V(t) I(t)
Step Honeywell 1
Input PID Controller fcn PID
Sum Product Iref L.s+R
position -> force Saturation PID Current Voltage
Electrical
controller Saturation
System

m*g
weight
current actual 1 v 1 x
limit 1/m
force s s
Sum2 1/m acc vel Distance
2-D Look-Up
Table disturbance


Physical System Model




Systems & Simulation ­ ENEL3SS 19
SS1_3.DOC
3.2 Solution of the linear state space equations

We will now take some care to be clear about the solution of the state differential equations.


3.2.1. Solution of a first order, time invariant, linear differential equation
Consider the equation,

x = ax+ bu
& (3.5)

with initial condition, x( t 0 ) = x 0 and u(t) given for all t0.
The general solution is the sum of the forced response and the natural response,

x( t ) = e ( 0 ) x 0 + e a ( t - ) b u ( ) d
a t- t t
(3.6)
t0


We can formally obtain this result by a number of means and should quickly revise the
possibilities.

1) Differentiate eq(3.5) to obtain eq(3.4): The first term is easy,
dt
(
d a( t- t 0 )
e )
x 0 = a e at x 0 . The
second part is more tricky if we want to be mathematically accurate:

d t a ( t - )
b u( ) d = b u( t ) + t a ea ( t - ) b u( ) d
t
dt
t0 e 0

= b u( t ) + a ( x( t ) - e a ( t - t 0 ) x 0 )
d t t
dt t 0
(Leibnitz' rule for a "smooth" function, f(t,), states: f ( t, ) d = f ( t, t ) + f ( t, ) d )
t 0 t

We see that eq(3.6) is a solution of eq(3.5) and invoke the uniqueness of the solution of linear
equations to show that it is indeed the only solution.
0 .5
(As an aside, a non-linear differential equation does not
always have a unique solution. It may have no solution, 0 .4

bi-furcations, chaotic solutions etc. For example, 0 .3
x = x , x0=0 has a bifurcation at t=1 as x=0 t is a
x, xdot




&
0 .2 xd o t


( )
2
solution but x = t-1 ,
2
t >1 is also a solution.)




2) Use Laplace transforms. To use the Laplace transform, we must be sure that the signals, u(t) and
2
x(t) (which we don't yet know) are Laplace transformable. (E.g. x = e t does not have a Laplace
transform. Why?) This is actually a mild condition for practical systems. Applying Laplace
transform to eq(3.5) gives,

sX(s ) - x 0 = a X(s) + b U(s )


Systems & Simulation ­ ENEL3SS 20
SS1_3.DOC
1
X( s ) = ( x 0 + b U(s)) (3.7)
s-a

Recalling that a product in the Laplace domain gives a convolution in the time domain, we will
obtain eq(3.5) on inverse transformation.

We notice that the behaviour of the solution depends very significantly on the value of "a" and that
(at least for u(t)=0) a>0 implies that the solution becomes unbounded as t. We called "a" a pole
(or a singularity) of eq(3.6).


3.2.2. The matrix case

Consider the state differential equation with u(t)=0.

x = A x , x( t 0 ) = x 0
& (3.8)
If x(t) is sufficiently smooth, we can approximate it by a power (Taylor) series around t0,

x( t ) = x t + x t ( t - t 0 ) + && t
& x (t - t 0 ) 2 2 ! + L +x ( n ) (t - t 0 ) n n !+ L (3.9)
0 0 0 t0


By successively differentiating eq(3.8) we get,

&& = A x = A 2 x , x ( n ) = A n x
x &

and substituting into eq(3.9) we have,


(
x ( t ) = I + A( t - t 0 ) + A 2 ( t - t 0 )
2
2 ! + L + A n (t - t 0 )
n
n !+L x0 )
By analogy with the series expansion for ea, we define the matrix exponent

eA =
& Ak k!
k =0
and the transition matrix,

( t ) = e
& At
= Ak tk k!
k =0
The unforced (natural) response is then

x( t ) = e A( t- t0 ) x 0 = ( t - t 0 ) x 0


Properties of the transition matrix:
(1) e0t = eA0 = I
d At
(2) e = Ae At = e At A
dt



Systems & Simulation ­ ENEL3SS 21
SS1_3.DOC
(Order of multiplication may not generally be reversed with matrix operations - A B B A in
general.)

(3) (e ) At -1
= e -At

(The inverse of the transition matrix always exists.)
If A is regular (i.e. non-singular),

( ) ( )
t A t
(5) 0 e d = A -1 e A = A -1 e At - I = e At - I A -1
0

If A is singular, the integral exists but must be found by means of a power series.
(6) e A( t + ) = e At e A but e A + B e A e B generally


The complete solution of the state differential equation is

x( t ) = e ( 0 ) x 0 + e ( ) B u( ) d
A t-t t A t -
(3.10)
t0
Eq(3.10) can be obtained via the Laplace transform. We assume that t0=0 for simplicity and take
Laplace transforms of eq(3.8):

s X ( s ) - x 0 = A X( s ) + B U ( s )
(3.11)
X( s ) = ( s I - A )
-1
(x 0 + B U(s))
Inverse Laplace transform of eq(3.11) will give eq(3.10). This also gives us a way of calculating
the transition matrix as

( t ) = e At = 1
{(sI - A) } -1
(3.12)




3.3. Conversion between state-space and transfer function descriptions

We are often interested in the behaviour of the system output as a response to the input and then
would like to find the transfer function

-1
P(s) = C (sI - A) B+ D (3.13)

We can go backwards and forwards between state-space and input-output descriptions. The
transformation is not unique and several useful canonic forms have been defined:




Systems & Simulation ­ ENEL3SS 22
SS1_3.DOC
3.3.1. Observer form state space description







L Y - b4 s + a n -1Y - b n -1U s + L s + a1Y - b1U s + a 0Y - b 0 U
123
4 nU
X
4444 244444
1 n 4 3
X n -1
4444444444 244444444444
O4
1 3
X1


s X1 = - a 0 Xn + ( b 0 - a 0 b n ) U
s X2 = X1 - a1 X n + ( b 1 - a1 b n ) U
M
sXn = X n -1 - a n -1 X n + ( b n -1 - a n -1 b n ) U
Y = Xn + bn U

0 -a 0 b0 - a0bn
1 0 - a1 b - a1 b n
1
&= O
x x + M u

0 O 0 -a n-2 b n - 2 - a n-2 b n
1 - a n -1 b n -1
- a n -1 b n
(3.14)
y = [ 0 L L 0 1] x + b n u



3.3.2. Controller form state space description

sn s n -1 s 1
Y = bn U + b n -1 U + L + b1 U + b0 U
D (s )
123 D (s )
123 D (s )
123 D (s )
123
sX n X n =sX n -1 X 2 =sX1 X1

D(s) = s n + a n -1s n -1 + L + a 1s + a 0




Systems & Simulation ­ ENEL3SS 23
SS1_3.DOC
sX1 = X 2
sX 2 = X 3
M (3.15)
sX n -1 = X n
sX n = -a0 X1 - a1 X 2 -L- an -1 X n + U
Y = ( b 0 - a 0 b n )X 1 + ( b 1 - a 1 b n )X 2 + L + ( b n -1 - a n -1 b n )X n + b n U
etc.




Systems & Simulation ­ ENEL3SS 24
SS1_3.DOC
3.3.3. Similarity transforms

The state variables of the state space description,
x = Ax + Bu
&
(3.16)
y = Cx + Du

can be transformed by the invertible similarity transform, T
z = Tx

to get

z = A' z + B'u
&
(3.17)
y = C'z + Du

with, A' = TAT-1, B' = TB, C' = CT-1
and with identical input-output behaviour (i.e. transfer function).




Systems & Simulation ­ ENEL3SS 25
SS1_3.DOC

4. INPUT-OUTPUT DESCRIPTIONS - TRANSFER FUNCTION OF LINEAR
SYSTEMS

Note: Matlab has a number of utilities for doing routine calculations used in this chapter. In
addition to being able to do these calculations by hand for small problems, you should be able to
use the following.

Command Function
residue Partial-fraction expansion or residue computation
ss2tf State-space to transfer function conversion
tf2ss Transfer function to state-space conversion
eig Eigenvalues and eigenvectors.
roots Polynomial roots.
step Step response.
bode Bode plot (frequency response)
logspace Logarithmically spaced vector

4.1. Transfer Functions
Very often, we are interested in the external behaviour of a system and the state variables are not
important to the problem on hand. Obviously the external behaviour cannot depend on the choice
of state variables (for example the outflow from a parallel sided tank under gravity with a square-
root flow law does not depend on whether we use head, volume or mass as the state variable). We
also saw that for linear systems there are canonic forms with the same external behaviour. We
define the "transfer function", P(s) of a system as the relationship between the Laplace transform of
the input signal and the Laplace transform of the output signal assuming zero initial conditions.

Y(s) = P(s) U(s) |zero i.c.'s (4.1)

For single input, single output (SISO) systems

P(s)=Y(s)/U(s) |zero i.c.'s (4.2)

For any state space description,
A B
S= ,
C D
-1
P(s) = C (sI - A) B + D (4.3)
It is important to note that we can apply the definition of the transfer function literally by actually
taking an input signal (with sufficient energy), calculating the corresponding output signal and then
dividing the Laplace transforms. Alternatively, (and more usually) we apply the Laplace transform
to the describing differential equation.




Systems & Simulation ­ ENEL3SS 25
SS1_4.doc
Example: Consider the system,

y + 3y = 2 u , y ( 0) = 0
& &

which has transfer function,

2
P( s) =
s+ 3

If we find the response to a bounded ramp,

u(t) = t (t) - (t-1) (t-1)
1 - e- s
U( s) =
s2

Y(s) = P(s) U(s) =
(
2 1 - e -s )
(s + 3)s2
( ) 1
= 2 1 - e -s 2 - +

1 1
9s 9(s + 3)

3s
2
( (
y( t ) = 3(t ( t ) - ( t - 1)( t - 1) ) - (t ( t ) - ( t - 1)( t - 1) ) + e -3t ( t ) - e -3( t -1) ( t - 1)
9
))
Bounded ramp input response to bounded ramp
0.7
1

0.9 0.6

0.8
0.5
0.7

0.6 0.4
u(t)




y(t)




0.5
0.3
0.4

0.3 0.2

0.2
0.1
0.1

0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t t

» t=linspace(0,2); % produces 100 equally spaced points
» u=t-(t-1).*(t>=1);
» y=2/9*(3*t-1+exp(-t*3));
» y(51:100)=y(51:100)-y(1:50); % for t>=1
» plot(t,u)
» title('bounded ramp'), ylabel('u(t)'), xlabel('t')
» plot(t,y)
» title('response to bounded ramp'), ylabel('y(t)'), xlabel('t')




Systems & Simulation ­ ENEL3SS 26
SS1_4.doc
4.2. Frequency response and Bode plots
In control design we will often use the so-called frequency response. If a stable (to be defined
below) linear system is excited with a steady sinusoidal input, the output will also be sinusoidal
after initial transients have died out.

2
Example: The system above was P( s) = . Take an input signal u(t) = A sin(t)
s+ 3
A
U(s) = 2 and zero initial conditions (for definition of the transfer function to be valid).
s + 2

The output is
Y(s) = P(s) U(s)
2 A
=
s + 3 s2 + 2
a bs + c
= + 2
s + 3 s + 2
2 A 1 -s + 3
= 2 s+3
+ 2
9+ s + 2

2A 2A
y( t ) = e -3t + ( - cos(t ) + 3 sin(t ))
9 + 2 9 + 2
2A -3t 2A
=
2
e + (sin(t - arctan( / 3)))
9+ 9 + 2

(recall X sin + Y cos = X2 + Y2 sin( + arctan(Y / X)) )

The first term is the transient response, and the second (which is of interest at present) is the
harmonic response.

The harmonic response is a sinusoid which has an amplitude A|P(j)| and a phase Arg(P(j))
where
j = -1
2 2
P( j ) = = e - j arctan( / 3)
j + 3 2
+9

We could have obtained the same result using phasor methods (from an Electricity course) or using
the Fourier transform (which is an integral defined over all time and therefore would not have had
the problem with the initial conditions.

The transfer function calculated with s=j is a very useful tool in understanding system behaviour.




Systems & Simulation ­ ENEL3SS 27
SS1_4.doc
4.2.1. Bode plots

In the 1940's (before PC's) all calculations were done by hand. Bode worked at Bell Laboratory in
the USA on problems related to amplifier performance (before transistors and op amps!). We will
discover that there are very useful relationships between the frequency domain and time domain
behaviour and the Bode plot is a simple way of representing behaviour in the frequency domain. A
Bode plot displays |P(j)| in dB (defined below) and Arg(P(j)) as a function of frequency using
semi logarithmic paper so that large frequency range can be represented. (Incidentally the
frequency of notes on a piano is logarithmically spaced.)

Definition - deciBell (dB)

x dB = 20 log10 x

( )
= 10 log10 x2
(We remember the first formula but the second shows that the 20 is not arbitrary: 10 = deci and the
power 2 indicates that the dB is a measure of power (i2R, V2/R etc).)

Discussion
Any polynomial can be factored into a product of first and second order terms

m
b h sh m
h =1(s + z h ) n
r
i.e. P(s) = h=0
n
=
n
= (s +hp +k
a ks k k =1(s + p k ) k =1 k)

k =0
Where s = -pk are the "poles" and s = -zh are the "zeros". If there is a second order term, say in the
denominator, there is a complex conjugate pair of poles, pr = p*r+1

This means that we only need to be able to draw first and second order factors and add since
multiplication becomes addition in the complex logarithmic plane (dB and degrees)
2
Example: Consider P( s) = above.
s+ 3
2
P ( j ) = Arg( P( j )) = - arctan( / 3)
2 + 9

Considering |P(j)|,

<<3 |P(j)| 2/3
|P(j)|dB = 20 log10(2/3) = -3.52dB
>>3 |P(j)| 2/
|P(j)|dB = 20 log10(2) - 20 log10()
3 |P(j)| 2/ 3 2
|P(j)|dB = 20 log10(2/3)-20 log10(2) = (-3.52-3.01)dB
= -6.53dB



Considering Arg(P(j)),
Systems & Simulation ­ ENEL3SS 28
SS1_4.doc
<<3 Arg(P(j) = -arctan(/3) 0°
>>3 Arg(P(j) = -arctan(/3) -90°
3 Arg(P(j) = -arctan(/3) -45°

P (s )= 2/(s + 3)
0


-5


-10


-15
|P| [dB]




-20


-25


-30


-35
-1 0 1 2
10 10 10 10
freq [rad/s ]
P (s )= 2/(s + 3)
0




-20
Arg(P) [deg]




-40




-60




Figure - First order example with gain of 2/3 and corner frequency of 3 rad/s




Systems & Simulation ­ ENEL3SS 29
SS1_4.doc
1
FIRST ORDER BODE PLOTS P(s) =
s + 1
|1/(s/wn+1)|
0




-5
magnitude [dB]




-10




-15




-20
-1 0 1
10 10 10
normalised frequency w/wn

Arg(1/(s/wn+1))



-10


-20


-30
angle [deg]




-40

-50




Systems & Simulation ­ ENEL3SS 30
SS1_4.doc
SECOND ORDER BODE PLOTS
(
1 (s / n )2 + 2s / n s + 1 )dB
=0.1
10
=0.2
=0.3
=0.4
0 =0.5
=0.7
=1.0
magnitude [dB]




-10



-20



-30



-40 -1 0 1
10 10 10
normalised frequency /n

((
Arg 1 (s / n )2 + 2s / n s + 1 ))
=0.1
-20 =0.2
=0.3
=0.4
-40 =0.5
=0.7
-60 =1.0
angle [deg]




-80

-100

-120

-140

-160

-180 -1 0 1
10 10 10
normalised frequency /n

Systems & Simulation ­ ENEL3SS 31
SS1_4.doc