Solution 3.8.6.11 The rst step is to get the data into MATLAB. The followng commands do that. EDU>load tfid EDU>size(tfid) ans = 500 2 EDU> the rst column is the time, and the second column is the step response. Weusethe following commands to put the data in a form wecan plot the step response versus time. EDU>t = tfid(1:500,1);; EDU>y = tfid(1:500,2);; EDU>plot(t,y) EDU>print -deps sr38611a.eps EDU> The time response is shown in Figure`1 As can be seen the time response, whichwas taken o the oscilloscope shows negativetimebeforethe step is applied and also the system does not start exactly at zero. The downward trend just before t =0iswherethe switchwas thrown in the analog circuit to initiate the step. So wehave some choices to make. We need to clean the data up a bit to apply the identi cation technique described in this chapter. First of all, weadd0.12 seconds to the t vector so that t(75) becomes 0. Then wesety(75) equal to zero and t(75) equal to zero. Then wetruncate the t and y vectors to the entries in locations 72 through vehundred. The commands are: t=t+0.12;; t(75) = 0;; y(75) = 0;; t=t(75:500);; y=y(75:500);; plot(t,y);; print -deps sr38611b.eps;; 1 -2 0 2 4 6 8 10 -0.5 0 0.5 1 1.5 2 2.5 Figure 1: Rawstep response data 0 1 2 3 4 5 6 7 8 9 0 0.5 1 1.5 2 2.5 Figure 2: Adjusted step response 2 The adjusted response is shown in Figure 2. The rationale for the adjust- ments is as follows. The response was generated bycreating a second order comepensator on the Feedback and Control box(FAC). The system response at steady state is almost exactly four in reponse to a 1 V input. Thus, all wewanttodoisgetzero to be the pointatwhich the switchonthe FAC was thrown to initiate the step and the response goes down close to zero. This turns out to be approximately t = ;1:8s.Wearenowinaposition to apply the identi cation techniques discussed in Chapter 3. As a rst step, wecreate the function y 1 = A;y(t) and then w(t)=ln(y(t));; as shown bytheMATLAB dialogue: EDU>max(y) ans = 2.0100 EDU>y1 = 2.02 -y;; EDU>w =log(y1);; EDU>plot(t,w) EDU>print -deps lny1.eps EDU> Note that 2.01 was increased to 2.02 to avoid ln(0). The result is shown in Figure 3. The slope is roughly ; 3:5 2:1 = ;1:67 So there is a pole around s = ;2. Our next goal is to nd B using the MATLAB dialogue: B=(y-2.01)./exp(-1.67*t);; EDU>plot(t(25:200),B(25:200)) EDU>print-deps Ba.eps EDU>mean(B(50:75)) ans = 3 0 1 2 3 4 5 6 7 8 9 -5 -4 -3 -2 -1 0 1 Figure 3: Plot of ln(2:02;y(t) 4 0 0.5 1 1.5 2 2.5 3 3.5 4 -25 -20 -15 -10 -5 0 Figure 4: B versus t, p 1 =0:42 -2.2735 EDU>mean(B(50:100)) ans = -2.3649 EDU> Note that wehave used the MATLAB command: ./ to do an elementbyelement division using two matrices. The plotofB versus time is shown in Figure 4. Wechoose B =2:36: 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -9 -8 -7 -6 -5 -4 -3 -2 -1 Figure 5: Plot of ln[y 2 (t)] Wecomputed the mean in the region just before steady state, namely for 1 <t<2 s: Having found B we can nowattempt to nd p 2 .Todosoweform the function y 2 (t)=y(t);(A + Be ;p 1 t = Ce ;p 2 t : In the presentcase y 2 (t)=y(t);2:01+ 2:36e ;1:67t : We then take the natural logarithm of y 2 and plot it versus time. as shown in Figure 5. It should be clear from the gure that wecannot identify the second pole. The response due to the second pole dies out too quickly.To haveanychance wewould need to sample faster. Wegoahead and compute C = ;(A + B)=;(2:01;2:36) = 0:35;; and let p 2 =10. Weknow p 2 is fairly large so this guess is, if anything conservative. Then our estimate of y(t)is ^y(t)=2:01;2:36e ;0:42t +0:35e ;10t : 6 0 1 2 3 4 5 6 7 8 9 -0.5 0 0.5 1 1.5 2 2.5 Figure 6: Plot of ^y(t)versusy(t) A comparison of y(t) and ^y(t)isshown in Figure 6. As can be seen the t is very good, and wecould stickwith this estimate of the transfer function. There is alternativeapproach. Wehaveavery good estimate of p 1 . Suppose wepickavalue for p 2 , nd the gain K,and then build and test the transfer function using MATLAB. For instance, weletp 2 = 20. Then we calculate the mean steady state value: EDU>size(y) ans = 426 1 EDU>mean(y(350:426)) ans = 1.9913 EDU> 7 0 1 2 3 4 5 6 7 8 9 0 0.5 1 1.5 2 2.5 Figure 7: Plot of ^y(t)versusy(t) using alternativemethod p 2 =10 Then K p 1 p 2 =1:9913;; or K =1:67101:9913 = 33:25: Then weusethe MATLAB dialogue: EDU>g =zpk([],[-1.67 -10],33.25) Zero/pole/gain: 33.25 --------------- (s+1.67) (s+10) EDU>[yhat1,T]=step(g,t);; EDU>plot(t,y,T,yhat1) EDU> to generate the plot of y versus ^y shown in Figure 7. The t is prettygood but using the MATLAB dialogue: 8 0 1 2 3 4 5 6 7 8 9 0 0.5 1 1.5 2 2.5 Figure 8: Plot of ^y(t)versusy(t) using alternativemethod p 2 =20 EDU>g =zpk([],[-1.67 -20],66.509) Zero/pole/gain: 66.509 --------------- (s+1.67) (s+20) EDU>[yhat2,T]=step(g,t);; EDU>plot(t,y,T,yhat2) EDU> whichgenerates the plot shown in Figure 8. and the t is perfect. The only y in the ointmentisthat using the MATLAB dialogue: EDU>g =zpk([],[-1.67 -25],83.137) Zero/pole/gain: 83.137 --------------- 9 0 1 2 3 4 5 6 7 8 9 0 0.5 1 1.5 2 2.5 Figure 9: Plot of ^y(t)versusy(t) using alternativemethod p 2 =25 (s+1.67) (s+25) EDU>[yhat3,T]=step(g,t);; EDU>plot(t,y,T,yhat3) EDU> wegenerate the plot shown in Figure 9, whichisjust as good a t. Thus, wede nitely know that p 2  20: The most conservativemodelwould then be G(s)= 66:51 (s +1:67)(s+20) : Wecould also experimentwiththe value of p 1 , but the model wehave found would be adequate for designing control systems. 10