Matlab tutorial

Document Sample
Matlab tutorial Powered By Docstoc
					                     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
     rss=10;
     Kp=0:1:25;
     for i=1:length(Kp)
           FRz = feedback(Kp(i)*Gz, Hz);
           ess(i)=rss - rss*dcgain(FRz);
     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 tang shuming
About