Solution 8.8.3mtr2y The rst step is to get the data into MATLAB. The MATLAB statements load mtr6step size(mtr6step) t=mtr6step(1:500,1);; y=mtr6step(1:500,2);; plot(t,y) print -deps sr881a.eps seperates the data into a column vector of times and a column vector of responses, and then plots and saves the step response EDU>t = tfid(1:500,1);; EDU>y = tfid(1:500,2);; EDU>plot(t,y) EDU>print -deps sr883mtr2ya.eps EDU> shown in Figure`1 As can be seen the time response, whichwas taken o the oscilloscope, does not begin exactly at t =0,andnot quite at y =0.The MATLAB program t=t+.004;; y=y(49:500);; t=t(49:500);; y(1) = 0;; t(1) = 0;; plot(t,y) print -deps sr883mtr2yb.eps produces the adjusted reponse shonwinFigure 2. We are nowinaposition to apply the identi cation techniques discussed in Chapter 3. As a rst step, the MATLAB statements maxy = max(y) y1 = (max(y) + 0.0001) -y;; w=log(y1);; plot(t(1:300),w(1:300));; print -deps 883mtr2ylny1.eps 1 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 1: Rawstep response data produces the response shown in Figure 3. Note that the maximum resposne was increased to0.001 to avoid ln(0). The slope is roughly ; [1:5; (;1:2)] 0:3 = ;9: So there is a pole around s = ;9. Our next goal is to nd B using the MATLAB statements: s=-(1.5 +1.2)/0.3 Bv = (y - max(y))./exp(s*t);; plot(t(1:125),Bv(1:125)) print -deps 8813mtr2yB.eps B=mean(Bv(25:75)) A=mean(y(200:450)) Note that wehave used the MATLAB command: ./ 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 2: Adjusted step response to do an elementbyelement division using two matrices. The plotofB versus time is shown in Figure 4. The values obtained are A =4:3343 and B =4:5960: These values to be about right B because weknow jBj > jAj 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); (4:3343;4:596e ;9t : We then take the natural logarithm of y 2 and plot it versus time. as shown in Figure 5. This plot has a slope of about -10. Weconclude that the second pole is either very far to the left, or very close to the rst pole. slope to be Wewillhaveto\ sh" a bit here. So wegoahead and nd C = ;(A + B)=;(4:3343;4:596) = 0:2617: 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -4 -3 -2 -1 0 1 2 Figure 3: Plot of ln(1;y(t) 0 0.05 0.1 0.15 0.2 0.25 -5.6 -5.4 -5.2 -5 -4.8 -4.6 -4.4 -4.2 -4 -3.8 Figure 4: Comparison of adjusted response and y 1 (t)=1;e ;24:9t 4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 -1 -0.5 0 0.5 1 1.5 2 2.5 Figure 5: Plot of ln[y 2 (t)] and place the second pole at s = ;10. Then our estimate of y(t)is ^y(t)=4:3343; 4:596e ;9t +0:2617e ;10t : A comparison of y(t) and ^y(t)isshown in Figure 6. As can be seen the t is prettygood. However, if weusethe MATLAB statments K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,0.9,449);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr2yd.eps we get the responses shown in Figure 7, and the t is not as good. The culprit maybeour choices of B and C.Rather than go back and try to recompute B and C,we just sh around abit. If we put the poles at s = ;12 and s = ;30, weget the responses shown in Figure 8. The last taskisto ndK.Wehave K p 1 p 2 =4:3343;; 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 6: Plot of ^y(t)versusy(t) whichinthe presentcasebecomes K =4:33343 12 30 = 1560:34: However, Figure 9 shows the armature voltage for the recorded step re- sponse. Thus, we need to divide K by36toachieve the correct gain. Finally G(s)= 43:34 (s +12)(s++30) : 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 7: Plot of ^y(t)versusy(t) The entire MATLAB program that does this analysis is load mtr2step t=mtr2step(1:500,1);; y=mtr2step(1:500,2);; plot(t,y) print -deps sr883mtr2ya.eps t=t+.004;; y=y(49:500);; t=t(49:500);; y(1) = 0;; t(1) = 0;; plot(t,y) print -deps sr883mtr2yb.eps maxy = max(y) y1 = (max(y) + 0.0001) -y;; w=log(y1);; plot(t(1:300),w(1:300));; print -deps 883mtr2ylny1.eps 7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 8: Plot of ^y(t)versusy(t) s=-(1.5 +1.2)/0.3 Bv = (y - max(y))./exp(s*t);; plot(t(1:125),Bv(1:125)) print -deps 8813mtr2yB.eps B=mean(Bv(25:75)) A=mean(y(200:450)) y2 = y-A +B*exp(s*t);; w2 = log(y2);; plot(t(1:150),w2(1:150)) print -deps 883mtr2ylny2.eps C=-(A + B) s1 = -10 yhat = A + B*exp(s*t) + C*exp(s1*t);; plot(t,y,'k-',t,yhat,'k--') print -deps sr883mtr2yc.eps K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,0.9,449);; 8 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 5 10 15 20 25 30 35 40 Figure 9: Recorded armature voltage [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr2yd.eps s=-12 s1 = -30 K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,0.9,449);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr2ye.eps pause load mtr2av tav = mtr2av(1:500,1);; av = mtr2av(1:500,2);; plot(tav,av) print -deps 883mtr2yav.eps 9