Solution 3.8.6.10 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 sr38610a.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, we add 0.2534 seconds to the t vector so that t(74) becomes 0. Then wesety(74) equal to zero and t(74) equal to zero. Then wetruncate the t and y vectors to the entries in locations 72 through vehundred. The commands are: t=t+0.2534;; t(74) = 0;; y(74) = 0;; t=t(74:500);; y=y(74:500);; plot(t,y);; print -deps sr38610b.eps;; 1 -5 0 5 10 15 20 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Figure 1: Rawstep response data 0 2 4 6 8 10 12 14 16 18 0 0.2 0.4 0.6 0.8 1 1.2 1.4 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 = 1.0280 EDU>y1 = 1.029-y;; EDU>w = log(y1);; EDU>plot(t,w) EDU>print -deps lny1.eps EDU> Note that 1.028 was increased to 1.029 to avoid ln(0). The result is shown in Figure 3. The slope is roughly ; 1:25 3:1 = ;0:4 So there is a pole around s = ;0:5. Our next goal is to nd B using the MATLAB dialogue: EDU>B = (y-1.028)./exp(-0.42*t);; EDU>plot(B(150:200)) EDU>mean(B(150:200)) ans = 3 0 2 4 6 8 10 12 14 16 18 -7 -6 -5 -4 -3 -2 -1 0 1 Figure 3: Plot of ln(1:029; y(t) -1.1175 EDU>print -deps B.eps 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. Note wehave increased p 1 slightly to make jBj >A: Wealso computed the mean in the region just before steady state, namely for 6 <t<8 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 : 4 5.5 6 6.5 7 7.5 8 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 Figure 4: B versus t, p 1 =0:42 In the presentcase y 2 (t)=y(t);1:028+ 1:1175e ;0:42t : 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)=;(1:028; 1:1175) = 0:0895;; and let p 2 =10. Then our estimate of y(t)is ^y(t)=1:028;1:1175e ;0:42t +0:0895e ;10t : A comparison of y(t) and ^y(t)isshown in Figure 6. As can be seen the t is prettygood. 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: 5 0 0.5 1 1.5 2 2.5 3 3.5 4 -4.5 -4 -3.5 -3 -2.5 -2 Figure 5: Plot of ln[y 2 (t)] 0 2 4 6 8 10 12 14 16 18 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Figure 6: Plot of ^y(t)versusy(t) 6 EDU>size(t) ans = 427 1 EDU>mean(y(375:427)) ans = 1.0164 EDU> Then K p 1 p 2 =1:0164;; or K =0:43201:0164 = 8:74: Then weusethe MATLAB dialogue: EDU>size(t) ans = 427 1 EDU>mean(y(375:427)) ans = 1.0164 EDU>g = zpk([],[-0.43 -20],8.74) Zero/pole/gain: 8.74 --------------- (s+0.43) (s+20) 7 0 2 4 6 8 10 12 14 16 18 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Figure 7: Plot of ^y(t)versusy(t) using alternativemethod p 2 =20 EDU>[yhat,T]=step(g,t);; EDU>plot(t,y,T,yhat) EDU>print -deps yhat1.eps EDU> to generate the plot of y versus ^y shown in Figure 7. The only problem with this is given bytheMATLAB dialogue: EDU>g = zpk([],[-0.43 -50],21.853) Zero/pole/gain: 21.853 --------------- (s+0.43) (s+50) EDU>[yhat,T]=step(g,t);; EDU>plot(t,y,T,yhat) EDU> which generates the plot shown in Figure 8. There is not a lot to choose 8 0 2 4 6 8 10 12 14 16 18 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Figure 8: Plot of ^y(t)versusy(t) using alternativemethod p 2 =50 here. The t is marginally better for p 2 =50Ifwegothe other direction, assuming values of p 2 less than 20, we nd that the t is worse. Thus, we de nitely knowthat p 2 > 20: The most conservativemodelwould then be G(s)= 8:74 (s+0:43)(s+20) : We could also experimentwiththe value of p 1 , but the model would be adequate for designing control systems. 9