# Matlab tutorial

Document Sample

```					                     Matlab tutorial for Z-domain analysis

1. Overview of tutorial
-    Quick review on the chapter 4, 5, 6, 8, and 9
-    Specifying LTI system with transfer function and zero-pole-gain model
-    Building system from subsystem
-    Inferring system property
-    Viewing system response
-    P Controller Design
-    PI Controller Design
-    PD Controller Design
-    Controller Design using sisotool

2. Quick review of chapter 4, 5, 6, 8, and 9
See the power point slides

3. Specifying LTI System with Transfer Function and Zero-pole-gain
model
Commands to cover: tf, tfdata, zpk, zpkdata, celldisp

* specifying the LTI system with transfer function
 Gz= tf(0.47, [1 -0.43], -1);

* Get numerator and denominator from transfer function
 [num, den] = tfdata(Gz);

* Use celldisp to view the contents of cell array
 celldisp(num);
 celldisp(den);

* converting TF to zero-pole-gain model object
 Hz=zpk(Gz);

* back to transfer function
 Gz=tf(Hz)

1
* Get zero, pole, and gain from zero-pole-gain object
 [z, p, k] = zpkdata(Hz)
 [z, p, k] = zpkdata(Gz)

* specifying the system with zero-pole-gain model
 Hz = zpk(z, p, k);

* Using zero and pole
 pole(Gz);
 zero(Gz);

4. Building system from subsystem
Commands to cover: feedback and sisotool

* Concatenation

U(z)             G 1 (z)                G 2 (z)             Y(z)

G(z)

 Gz = tf(0.47, [1 -0.43], -1);
 Hz = tf(0.05, [1 -0.95], -1);
 Gz = Gz*Hz

* Summation

G1(z)

+
U(z)                                                 Y(z)
+
G2(z)

G(z)

2
 Gz= tf(0.47, [1 -0.43], -1);
 Gz = Gz+Hz

* Feedback Composition

U(z)                                       Y(z)
+                  G1(z)
-

G2(z)

G(z)

 Gz = tf(0.47, [1 -0.43], -1);
 Fz = feedback(Gz, Hz)

 from the textbook (pg 117)
FR(z) =            G(z)K(z)           feedback(K(z)*G(z), Hz)
1+H(z)G(z)K(z)
FD(z) =            G(z)                feedback(K(z)*G(z), Hz) / K(z)
1+H(z)G(z)K(z)
FN(z) =               1                feedback(K(z)*G(z), Hz) / K(z)*G(z)
(1+H(z)G(z)K(z)

Exercise 1)

U(z)                                   G(z)=          Y(z)
Kp =
+                               0.47
-                2.0       z-0.43

H(z) =
0.05
z-0.95

G(z)

 Kp=2.0;
 Gz=tf(0.47, [1 -0.43], -1);
 Hz = tf(0.05, [1 -0.95], -1);
 Fz = feedback(Kp*Gz, Hz)

3
Exercise 2)

PID Controller

KI =                z
2.0                z-1

R(z)           E(z)                                        U(z)   G(z)=    Y(z)
KP =                  +
+                                       +                   0.47
-                 0.5                   +
z-0.43

KD =               z-1
0.9                 z

H(z) =
0.05
z-0.95

 Ki=2.0
 Kp=0.5
 Kd=0.9
 tf1=tf([1 0],[1 -1], -1)
 tf2=tf([1 -1],[1 0], -1)
 Fz = feedback((Ki*tf1+Kp+Kd*tf2)*Gz,Hz)

sisotool: GUI-based LTI feedback system builder
 sisotool

5. Inferring System Properties
Commands to cover: dcgain, isct(isdt), pzmap, dsort
* Steady state gain = G(1)
 Gz
 dcgain(Gz)
 0.47/(1-0.43)
* Knowing whether the system is continuous or discrete
 Gz = tf(0.47, [1 -0.43], -1);
 isdt(Gz)
 isct(Gz)
 Gz=tf(0.47, [1 -0.43]);
 isdt(Gz)

4
 isct(Gz)
* Viewing pole-zero map
 Gz=zpk([0.3 0.1], [0.2+2j 0.2-2j], 1);
 pzmap(Gz);
* Sorting the discrete-time poles by magnitude
 Gz=zpk(0, [0.15 0.1+2j 0.1-2j], 1)
 p=pole(Gz)
 dsort(p)

6. Viewing system response to various signals
Commands to cover: impulse, step, lsim, ltiview

* Impulse signal
 Gz = tf(0.47, [1 -0.43], -1);
 impulse(Gz)
 Hz=tf(-0.014, [1 -0.59], -1)
 impulse(Hz)
 [y, t]=impulse(Gz)
 stem(t, y, 'filled')

* Step signal
 step(Gz)
 step(10*Gz)
 [y, t]=step(Gz, 20)
 stem(t, y, 'filled')

For other signals, use lsim(transfer function, input signal, time).

* Ramp signal
 t=0:10:200
 ramp=0:1:20
 lsim(Gz,ramp,t)

* Sinusoid signal
 t=0:20
 sinusoid = sin(t*pi/3)

5
 [y time]=lsim(Gz, sinusoid, t)
 stem(time, y, ‘filled’);

ltiview: Multi-purpose LTI system responses viewer
 ltiview(Gz)
 ltiview(Gz, Hz)
 ltiview(‘impulse’,Gz,Hz);
 ltiview(‘pzmap’,Gz,Hz);

Exercise)
Plot the system response for the following feedback system (Assume that reference
value is 10 and disturbance of 20 is introduced after 10 time units)

D(z) = 20 (after 10
time units)

R(z)=10             E(z)               U(z)                G(z)=     Y(z)
Kp =                +
+                                 +             0.47
-              2.0                         z-0.43

1) Get the system response with respect to the reference value
 FRz=feedback(Kp*Gz,1)
 [yr,t] = step(10*FRz,30)

2) Get the system response with respect to the disturbance
 FDz = FRz/Kp;
% (why? See page 3 of tutorial)
 d=zeros(1,30);
 for i=10:length(d)
       d(i)=20;
 end;
 [yd] = lsim(FDz, d, t);

6
3) Plot sum of two response
    stem(t, yr+yd, ‘filled’);

Now we need to add moving average filter in our system. There are three possible
location of the filter: before controller, after target system, and as a transducer.
Compare the response at different locations of filter. Transfer function of filter is (1-c)/z-
c where c=0.95. (See figure 8.36 in page 290).

 c=0.95;
 Hz=tf(1-c, [1 -c], -1);

1) before controller
 hold on
 Kp=Kp*Hz;
 FRz=feedback(Kp*Gz,1);
 [yr, t] = step(10*FRz,30);
 FDz=FRz/Kp;
 [yd]=lsim(FDz,d,t);
 stem(t, yr+yd,'filled',’r’)

2) after target system
 Kp=2.0;
 Gz=Gz*Hz;
 FRz=feedback(Kp*Gz,1);
 [yr, t] = step(10*FRz,30);
 FDz=FRz/Kp;
 [yd]=lsim(FDz,d,t);
 stem(t, yr+yd,'filled',’g’)

3) as a transducer
 Gz=tf(0.47, [1 -0.43], -1);
 FRz=feedback(Kp*Gz,Hz);
 [yr, t] = step(10*FRz,30);
 FDz=FRz/Kp;
 [yd]=lsim(FDz,d,t);
 stem(t, yr+yd,'filled','b')

7
7. P-Controller Design

D(z)

R(z) =                                               G(z)=
E(z)            U(z)           +               Y(z)
10       +              Kp = ?              +        0.47
-                                      z-0.43

H(z)=
0.05
z-0.95

Four properties of Controller: Stability, Accuracy, Settling time, Overshoot
Goal of design is to find the value of P-controller gain Kp.

1) Stability : Use rlocus to find the range of KP where poles are within unit circle
 c=0.95
 Hz=tf(1-c, [1 -c], -1);
 Gz=tf(0.47, [1 -0.43], -1);
 rlocus(Gz*Hz)
2) Accuracy: We want to minimize ess = rss-yss. Find ess with varying KP
 Kp=0:1:25;
 for i=1:length(Kp)
       FRz = feedback(Kp(i)*Gz, Hz);
 end;
3) Settling time: KS depends on the magnitude of the largest closed-loop pole
KS = -4 / log(r), where r is magnitude of the largest closed-loop pole
 for i=1:length(Kp)
       FRz=feedback(Kp(i)*Gz,Hz);
       mag_larg_pole = max(abs(pole(FRz)));
       Ks(i) = -4/log(mag_larg_pole);
 end;
4) Maximum Overshoot
Mp ≈ 0        real dominant pole and >= 0
|p| real dominant pole and < 0

8
rpi/|Ɵ| dominant poles p1, p2=re +/- jƟ

for dominant poles, c+/- dj,

r = sqrt(c2+d2), Ɵ=tan-1(d/c)

 for i=1:length(Kp)
      FRz=feedback(Kp(i)*Gz,Hz);
      pole_vector = pole(FRz)
      [r, index]=max(abs(pole_vector));
% if dominant pole is complex, then r become sqrt(c2+d2);
      theta=angle(pole_vector(index));
% theta =0 if real dominant pole r > 0
% theta = pi if real dominant pole r < 0
% theta = tan-1(d/c) if complex dominant pole(c+/- dj)
      if(abs(theta)<0.00)
          mp(i)=0;
      else if(abs(theta - pi)<0.001)
          mp(i)=r;
      else
          mp(i) = r^(pi/theta);
      end; end;
 end;
4) plot the vectors, ess, Ks, and mp, choose appropriate value for Kp
 hold on;
 plot(Kp, ess, ‘r’);
 plot(Kp, Ks, ‘g’);
 plot(Kp, mp, ‘b’);
   axis([0 20 0 25]);

8. PI-Controller Design

9
D(z)

R(z) =                                                              G(z)=
E(z)      (K P +K I )z – K P        U(z)       +            Y(z)
10      +                                                  +        0.47
-                   z-1                                z-0.43

Goal is to find the values for KP and KI
Zero: KP/( KP +KI)
Overall gain: KP +KI

1) Pick the location of zero with respect to the open-loop pole (0.43). In this case,
we pick one to the left, and the other to the right of pole. Draw root locus with
chosen zero. See the value of possible closed-loop pole.
 zero_value = 0.2;
 Kpi = tf([1 –(zero_value)], [1 -1], -1);
 rlocus(Kpi*Gz);
 zero_value = 0.7;
 Kpi = tf([1 –(zero_value)], [1 -1], -1);
 rlocus(Kpi*Gz);
2) Now, we know the zero at the left of open-loop pole is better in terms of fast
transient response. Pick some values of the zero within the range (< 0.43). In
this case, we choose 0.1, 0.2, and 0.4. Now for each of the zero location,
compute expected settling time and overshoot as the overall gain varies (note:
PID controller has ess=0).
 zero_value=[0.1 0.2 0.4];
 overall_gain=0.1:0.5:6;
 Ks=[]; mp=[];
 for i=1:length(zero_value)
      for j=1:length(overall_gain)
              Kpi=tf([overall_gain(j) -(overall_gain(j)*zero_value(i))], [1 -1],-1);
              FRz=feedback(Kpi*Gz, 1);
              pole_vector=pole(FRz)
              [r, index]=max(abs(pole_vector));
              Ks(i,j)=-4/log(r);
              theta=angle(pole_vector(index));

10
                 if(abs(theta)<0.00) mp(i,j)=0;
                 else if(abs(theta-pi)<0.001) mp(i,j)=r;
                 else
                       mp(i,j)=r^(pi/theta);
                 end;
                 end;
        end;
 end;
3) Plot the result and find the appropriate values of zero and overall gain
(compare the results with figures in page 310 of textbook)
 plot(overall_gain,Ks(1,:),'r')
 hold on
 plot(overall_gain,Ks(2,:),'g')
 plot(overall_gain,Ks(3,:),'b')
 plot(overall_gain,mp(1,:),'--r')
 plot(overall_gain,mp(2,:),'--g')
 plot(overall_gain,mp(3,:),'--b')
4) Solve the equation to get the value of KP and KI(Zero: KP/( KP +KI, Overall gain:
KP +KI)

9. PD-Controller Design

D(z)

R(z) =                                                               G(z)=
E(z)       (KP+KD)z – KD       U(z)       +                  Y(z)
10        +                                             +             1
-                   z                             Z2-1.3z+0.49

zero = KD/( KP +KD)
overall gain = KP +KD

1) Pick the candidate location of PD-controller’s zero and the range of overall gain.
Note that zero must be positive (KP>0 and KD>0).
 Gz=tf(1, [1 -1.3 0.49],-1)

11
 rlocus(Gz)
 Kpd=tf([1 -0.3],[1 0],-1);
 rlocus(Kpd*Gz)
 Kpd=tf([1 -0.8],[1 0],-1);
 rlocus(Kpd*Gz)
 zero_value=0.1:0.3:1.0
 overall_gain=0.1:0.1:0.9
2) For each zero location, plot the settling time and maximum overshoot.
 Ks=[]; mp=[];
 for i=1:length(zero_value)
      for j=1:length(overall_gain)
         Kpd = tf([overall_gain(j) -(zero_value(i)*overall_gain(j))], [1 0], -1);
         FRz=feedback(Kpd*Gz, 1);
         pole_vector=pole(FRz);
         [r, index]=max(abs(pole_vector));
         Ks(i,j)=-4/log(r);
         theta=angle(pole_vector(index));
         if(abs(theta)<0.00) mp(i,j)=0;
         else if(abs(theta-pi)<0.001) mp(i,j)=r;
         else
             mp(i,j)=r^(pi/theta);
         end;end;
      end;
 end;
 plot(overall_gain,Ks(1,:),'r')
 hold on
 plot(overall_gain,Ks(2,:),'g')
 plot(overall_gain,Ks(3,:),'b')
 plot(overall_gain,Ks(4,:),'c')
 plot(overall_gain,mp(1,:),'--r')
 plot(overall_gain,mp(2,:),'--g')
 plot(overall_gain,mp(3,:),'--b')
 plot(overall_gain,mp(4,:),'--c')
3) Also plot the actual value of maximum overshoot and compare it with estimated
value
 for i=1:length(zero_value)

12
       for j=1:length(overall_gain)
             Kpd = tf([overall_gain(j) -(zero_value(i)*overall_gain(j))], [1 0], -1);
             FRz=feedback(Kpd*Gz, 1);
             [y,t]=step(FRz);
             mp_actual(i,j)=(max(y) - dcgain(FRz))/dcgain(FRz);
       end;
 end;
 hold on;
 plot(overall_gain,mp(1,:),'r')
 plot(overall_gain,mp(2,:),'g')
 plot(overall_gain,mp(3,:),'b')
 plot(overall_gain,mp(4,:),'c')
 plot(overall_gain,mp_actual(1,:),'--r')
 plot(overall_gain,mp_actual(2,:),'--g')
 plot(overall_gain,mp_actual(3,:),'--b')
 plot(overall_gain,mp_actual(4,:),'--c')
 axis([0.1 0.7 0 1])
5) See the system response to the reference (=10) and disturbance (=20 after 10
time units)
(We chose zero=0.7, overall gain0.2, and thus KP=0.06, KD=0.14)
 Kp=0.06; Kd=0.14;
 Kpd = tf([(Kp+Kd) -(Kd)], [1 0], -1);
 FRz=feedback(Kpd*Gz, 1);
 [yr, t] = step(10*FRz,30);
 FDz=FRz/Kpd;
 [yd]=lsim(FDz,d,t);
 stem(t, yr+yd,'filled')

10.Controller design using sisotool
   sisotool(Gz)

1) Import the controller transfer function, and transducer(if exists)
2) Change compensator format to zero-pole-gain model
3) Open step response viewer
4) In root locus editor, move zero and poles or set the current compensator(gain)
value

13
5) Find the best location of zero and gain

14

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 0 posted: 5/28/2013 language: English pages: 14
tang shuming