E58 - Controls and Matlab

Contents

Define a transfer functions

mySys1=tf([1 30],[1 8 15])    %overdamped
mySys2=tf([16],[1 1 16])     %underdamped
mySys1 =
      s + 30
  --------------
  s^2 + 8 s + 15
Continuous-time transfer function.

mySys2 =
       16
  ------------
  s^2 + s + 16
Continuous-time transfer function.

Step response

plot step response (and find characteristics)

step(mySys1)

plot multiple step responses

step(mySys1,mySys2)

Calculate values for a step function

[y,t]=step(mySys1);

Calculate values for a step function at chosen points

t=0:0.01:5;
y=step(mySys1,t);
plot(t,y);

Frequency Response

Plot frequency response

bode(mySys1)

Calculate values for freq response

[mag, phase, freq]=bode(mySys1);
m=mag(:);       %Find magnitude vector
p=phase(:);     %Find phase vector

Calculate freq response from 0.1 to 1000 rad/sec.

freq=logspace(-1,3);
[mag, phase]=bode(mySys1,freq);
m=mag(:);       %Find magnitude vector
p=phase(:);     %Find phase vector

Transfer function manipulation

Two transfer functions in series

SeriesSys=series(mySys1,mySys2)
SeriesSys=mySys1*mySys2
SeriesSys =
              16 s + 480
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 143 s + 240
Continuous-time transfer function.

SeriesSys =
              16 s + 480
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 143 s + 240
Continuous-time transfer function.

Two transfer functions in parallel

ParallelSys=parallel(mySys1,mySys2)
ParallelSys=mySys1+mySys2
ParallelSys =
      s^3 + 47 s^2 + 174 s + 720
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 143 s + 240
Continuous-time transfer function.

ParallelSys =
      s^3 + 47 s^2 + 174 s + 720
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 143 s + 240
Continuous-time transfer function.

Feedback loop (negative feedback is assumed)

FeedbackSys=feedback(mySys1,mySys2)
FeedbackSys =
      s^3 + 31 s^2 + 46 s + 480
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 159 s + 720
Continuous-time transfer function.

Feedback using algebra (need to cancel equal poles and zeros (minreal))

FeedbackSys=mySys1/(1+mySys1*mySys2)
FeedbackSys=minreal(FeedbackSys)
FeedbackSys =
       s^5 + 39 s^4 + 309 s^3 + 1313 s^2 + 4530 s + 7200
  ------------------------------------------------------------
  s^6 + 17 s^5 + 126 s^4 + 606 s^3 + 2577 s^2 + 8145 s + 10800
Continuous-time transfer function.

FeedbackSys =
      s^3 + 31 s^2 + 46 s + 480
  ----------------------------------
  s^4 + 9 s^3 + 39 s^2 + 159 s + 720
Continuous-time transfer function.

Analysis of system

Get numerator and denominator

[n,d]=tfdata(FeedbackSys,'v')
n=FeedbackSys.num{1}
d=FeedbackSys.den{1}
n =
         0    1.0000   31.0000   46.0000  480.0000
d =
    1.0000    9.0000   39.0000  159.0000  720.0000
n =
         0    1.0000   31.0000   46.0000  480.0000
d =
    1.0000    9.0000   39.0000  159.0000  720.0000

Poles of transfer function

roots(d)
damp(d)
damp(FeedbackSys)
ans =
  -5.2914 + 2.7255i
  -5.2914 - 2.7255i
   0.7914 + 4.4381i
   0.7914 - 4.4381i
                                                 
       Eigenvalue          Damping     Frequency                                             
 -5.29e+00 + 2.73e+00i     8.89e-01     5.95e+00  
 -5.29e+00 - 2.73e+00i     8.89e-01     5.95e+00  
  7.91e-01 + 4.44e+00i    -1.76e-01     4.51e+00  
  7.91e-01 - 4.44e+00i    -1.76e-01     4.51e+00  
(Frequencies expressed in rad/TimeUnit)
                                          
       Eigenvalue          Damping     Frequency                                                    
  7.91e-01 + 4.44e+00i    -1.76e-01     4.51e+00  
  7.91e-01 - 4.44e+00i    -1.76e-01     4.51e+00  
 -5.29e+00 + 2.73e+00i     8.89e-01     5.95e+00  
 -5.29e+00 - 2.73e+00i     8.89e-01     5.95e+00  
(Frequencies expressed in rad/seconds)

General Characterization

ltiview

Inverse Laplace Transform. Partial fraction expansion

[n1,d1]=tfdata(mySys1,'v');
[r,p,k]=residue(n1,d1)
r =
  -12.5000
   13.5000
p =
    -5
    -3
k =
     []

Symbolic manipulation

syms s
n1s=poly2sym(n1,'s')   % Polynomial in 's'
d1s=poly2sym(d1,'s')
H = simple(n1s/d1s)  % Transfer function
n1s =
s + 30
d1s =
s^2 + 8*s + 15
H =
(s + 30)/(s^2 + 8*s + 15)

Inverse laplace of transfer function

ilaplace(H)
ans =
(27*exp(-3*t))/2 - (25*exp(-5*t))/2

Step response and inverse laplace

StpRspS=(1/s)*H
StpRspT=ilaplace(StpRspS)
StpRspS =
(s + 30)/(s*(s^2 + 8*s + 15))
StpRspT =
(5*exp(-5*t))/2 - (9*exp(-3*t))/2 + 2

Plot the function

t=0:0.01:5;
y=eval(vectorize(StpRspT));
plot(t,y)

Manipulation of transfer functions

syms G H k
G=1/(s^2+s+16);
H=k*s+1
myTF=G/(1+G*H)
H =
k*s + 1
myTF =
1/(((k*s + 1)/(s^2 + s + 16) + 1)*(s^2 + s + 16))

Simplify tranfer function

k=10;
T = simple(myTF)
T = eval(T)
T =
1/(s + k*s + s^2 + 17)
T =
1/(s^2 + 11*s + 17)

Pretty

pretty(T)
 
        1 
  -------------- 
   2 
  s  + 11 s + 17

LaTex

latex(T)
ans =
\frac{1}{s^2 + 11\, s + 17}

Get numerator and denominator

[ns, ds] = numden(T) % Get numerator and denominator
ns =
1
ds =
s^2 + 11*s + 17

Convert back to numeric

H = tf(sym2poly(ns),sym2poly(simple(ds)));

Simulink

Take System from Simulink sys=linmod(‘MySys’); [n,d]=ss2tf(sys.a, sys.b, sys.c, sys.d)

Help

doc symbolic
doc control