Chapter 3 MATLAB Results

The collection of MATLAB statements and screen display for:

Figure 3.1a

>> %Figure 3.1a. Effect of increasing time constant

>> t=0:.2:12;
>> G=tf(1,[2 1]) %time constant = 2, pole at -1/2

Transfer function:
   1
-------
2 s + 1

>> y1=step(G,t);

>> G=tf(1,[1 1]) %time constant = 1, pole at -1

Transfer function:
  1
-----
s + 1

>> y2=step(G,t);

>> G=tf(1,[1/2 1]) %time constant = 1/2, pole at -2

Transfer function:
    1
---------
0.5 s + 1

>> y3=step(G,t);

>> G=tf(1,[1/4 1]) %time constant = 1/4, pole at -4

Transfer function:
    1
----------
0.25 s + 1

>> y4=step(G,t);

>> plot(t,y1,t,y2,t,y3,t,y4)
>> legend('pole at -1/2','pole at -1','pole at -2','pole at -4')

Figure 3.1b

>> %Effect of changing the steady state gain
>> t=0:0.2:10;
>> G=tf(0.5,[1.5 1]) %Steady state gain = 0.5

Transfer function:
   0.5
---------
1.5 s + 1

>> y1=step(G,t);

>> G=tf(1,[1.5 1]) %Steady state gain = 1

Transfer function:
    1
---------
1.5 s + 1

>> y2=step(G,t);

>> G=tf(2,[1.5 1]) %Steady state gain = 2

Transfer function:
    2
---------
1.5 s + 1

>> y3=step(G,t);
>> plot(t,y1,t,y2,t,y3)
>> legend('K = 0.5','K = 1','K = 2')

Figure 3.2 with different damping ratios

>> %Effect of damping ratio
>> t=0:0.5:30;
>> G=tf(1,[1 0.4 1])  %damping ratio = 0.2

Transfer function:
       1
---------------
s^2 + 0.4 s + 1

>> damp(G)

       Eigenvalue          Damping     Freq. (rad/s)
 -2.00e-01 + 9.80e-01i     2.00e-01       1.00e+00
 -2.00e-01 - 9.80e-01i     2.00e-01       1.00e+00
>> %The Eigenvalue is the pole of G

>> y1=step(G,t);   %Save this step response for plotting later

>> %Repeat calculation wtih other damping ratios
>> G=tf(1,[1 0.8 1])   %damping ratio = 0.4

Transfer function:
       1
---------------
s^2 + 0.8 s + 1

>> damp(G)

       Eigenvalue          Damping     Freq. (rad/s)
 -4.00e-01 + 9.17e-01i     4.00e-01       1.00e+00
 -4.00e-01 - 9.17e-01i     4.00e-01       1.00e+00

>> y2=step(G,t);

>> G=tf(1,[1 2 1]) %Critically damped

Transfer function:
      1
-------------
s^2 + 2 s + 1

>> pole(G)

ans =
    -1
    -1

>> y3=step(G,t);
>> G=tf(1,[1 4 1]) %Overdamped

Transfer function:
      1
-------------
s^2 + 4 s + 1

>> pole(G)

ans =
   -3.7321
   -0.2679

>> y4=step(G,t);
>> plot(t,y1,t,y2,t,y3,t,y4)
>> legend('zeta = 0.2','zeta = 0.4','zeta = 1','zeta = 2')

Chapter 3 Review Problems

>> %Review problem 9
>> td=3;
>> P1=tf([-td/2 1],[td/2 1]);
>>
>> %Plotting deviates from text so we can superimpose the step input
>> t=0:.2:10;
>> y=step(P1,t);
>> plot(t,y) %Note how the response starts from negative values
>> hold
Current plot held
>> plot([0 3 3 9],[0 0 1 1]) %Plot a step input with a dead time of 3
>> title('Compare a first order Pade approx with a real time delay')

>> t=0:0.1:50;
>> taup=10;
>> G1=tf(1,[taup 1]);
>>
>> %y1 is first order with Pade approx of dead time
>> y1=step(G1*P1,t);
>> %y2 is just the function. Will be plotted with a shifted time axis
>> y2=step(G1,t);
>> t2=t+td;
>> plot(t,y1, t2,y2);
>> hold
Current plot held
>> plot([0 20],[0 0]);
>> legend('Using Pade approx','Theoretical time delay')
>> %Note how the Pade approx has a dip at the beginning

Example 3.2

>> %Review Problem 10 -- Example 3.2
>> q=3;
>> p2=conv([1 1],[3 1]);
>> p4=conv(conv(p2,[0.5 1]),[0.1 1]);
>> roots(p4) %check the poles

ans =

  -10.0000
   -2.0000
   -1.0000
   -0.3333

>> G2=tf(q,p2); % Second order reduced model
>> G4=tf(q,p4); %The original full order function

>> t=0:0.2:18;
>>
>> %Plotting here deviates from text for a better image
>> y4=step(G4,t);
>> td=0.1+0.5;
>> P1=tf([-td/2 1],[td/2 1]);
>> y2=step(P1*G2,t); %Second order with dead time approx.
>> plot(t,y4,t,y2)
>> %We'll do this again below with the first order approximation
>>
>> %We now do the more common first order with dead time approximation
>> td=1+0.1+0.5;
>> P1=tf([-td/2 1],[td/2 1]);
>> G1=tf(q,[3 1]); % First order approximation
>> y1=step(P1*G1,t);
>> plot(t,y4,t,y2,t,y1)
>> legend('Full model','2nd order with dead time','1st order with dead time')



>> %The second order with dead time approximation is quite good.
>> %The first order with dead time is very poor especially near the
>> %beginning. Nevertheless, the first order with dead time is the typical
>> %approximation used in empirical controller tuning relations


>> %Review Problem 11 - observe the effects of multicapacity models
>> tau=3;
>> G=tf(1,[tau 1]);
>> step(G); %First order function
>> hold
Current plot held
>> step(G*G); %Second order function
>> step(G*G*G) %Third oder function
>> step(G*G*G*G) %Fourth order function
>> step(G*G*G*G*G) %Fifth order function


>> %As the order increases, the shape of the step response resembles more of a
>> %first order with dead time approximation. Can we do it here? Let's try.
>>
>> %Repeat the fifth order calculation to get a time axis
>> [y5,t]=step(G*G*G*G*G); 
>>
>> %For the first order with dead time approximation, we'll try with
>> %a time constant of 4*tau = 12, and a dead time of td = tau = 3.
>> %No particular good reason. See comments with plot below.
>> G1=tf(1,[12 1])

Transfer function:
   1
--------
12 s + 1

>> y1=step(G1,t);  %Using a time shift axis to do the first order
>> t1=t+3;         % with dead time plot
>> plot(t,y5, t1,y1)

Note: The first order with dead time approximation does a poor job here. Why? That's because we don't have a dominant time constant in this example. The full order function is with a pole at -1/tau repeated five times. We, of course, could do a data fitting exercise and find a better choice of time constant and dead time for the approximation. But the result will never match the shape of the fifth order function without a dominant pole. That's one reason why we casually picked 12 and 3 as the values for this illustration.
 

>> %Review Problem 12
>> tz=3;
>> G=tf([tz 1],[5 1]);
>> step(G);

>> %But we can do better and recreate Fig. 3.4
>> G=tf([4 1],[2 1])

Transfer function:
4 s + 1
-------
2 s + 1

>> step(G)
>> hold
Current plot held

>> G=tf([3 1],[2 1])

Transfer function:
3 s + 1
-------
2 s + 1

>> step(G)

>> G=tf([2 1],[2 1]) %Case of pole-zero cancellation

Transfer function:
2 s + 1
-------
2 s + 1

>> step(G)   %This is the horizontal line at y=1

>> G=tf([-1 1],[2 1])

Transfer function:
-s + 1
-------
2 s + 1

>> step(G)

>> G=tf([-2 1],[2 1])

Transfer function:
-2 s + 1
--------
2 s + 1

>> step(G)

>> G=tf([-4 1],[2 1])

Transfer function:
-4 s + 1
--------
2 s + 1

>> step(G)

>> %Review Problem 13 - parallel functions
>> k1=3;
>> k2= -1; %k2 is negative
>> tau1=2; tau2=0.5;

>> %To repeat calculations with different values as in the text,
>> %it is best to put the following  statements in an M-file.
>> %Or copy-and-paste from a scratch file.
>> k=k1+k2;
>> tz=(k1*tau2+k2*tau1)/k;
>> G=tf(k*[tz 1], conv([tau1 1],[tau2 1]) );
>> step(G);
>> -1/tz          %Calculate the zero
ans =
     4

Note: The statements that we used above are clumsy because we want
to see the value of the zero. If we just want to do the plotting,
a cleaner approach is listed below. This way, we get to see the
curves of the individual contribution too.

>> G1=tf(k1, [tau1 1]);
>> G2=tf(k2, [tau2 1]);
>> G=G1+G2;
>> step(G,'-', G1,'--',G2,'--') 
>>     %Style: '-' means solid, '--' means dashed