-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 1: Rawstep response data Solution 8.8.3mtr4y The rst step is to get the data into MATLAB. The MATLAB statements load mtr4step size(mtr4step) t=mtr4step(1:500,1);; y=mtr4step(1:500,2);; EDU>t = tfid(1:500,1);; EDU>y = tfid(1:500,2);; EDU>plot(t,y) EDU>print -deps sr883mtr4ya.eps EDU> seperates the data into a column vector of times and a column vector of responses, and then plots and saves the step response. 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+.002;; 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 2: Adjusted step response y=y(49:500);; t=t(49:500);; y(1) = 0;; y(2) = 0.13;; t(1) = 0;; plot(t,y) print -deps sr883mtr4yb.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 883mtr4ylny1.eps produces the response shown in Figure 3. Note that the maximum resposne was increased to0.001 to avoid ln(0). 2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 -10 -8 -6 -4 -2 0 2 Figure 3: Plot of ln(1;y(t) The slope is roughly ;5=0:175 = 28:57: So there is a pole around s = ;29. Our next goal is to nd B using the MATLAB statements: s=-4/0.15 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 883mtr4yB.eps B=mean(Bv(90:140)) A=mean(y(250:450)) Note that wehave used the MATLAB command: ./ to do an elementbyelement division using two matrices. The plotofB versus time is shown in Figure 4. The values obtained are A =4:2719 and B = ;14:9691: 3 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -25 -20 -15 -10 -5 0 Figure 4: Comparison of adjusted response and y 1 (t) These values havetheright relationship, at least, 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:2719; 14:9691e ;28:5714t : 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 -25. Weconclude that the second pole is either very far to the left, or very close to the rst pole. We will haveto\ sh" a bit here. So wegoaheadand nd C = ;(A + B)=;(4:2719;14:9691) = 10:6973: and place the second pole at s = ;40. Weknow the pole is farther out, so we start with this guess. Then our estimate of y(t)is ^y(t)=4:2719; 14:9691e ;27t +10:6972e ;40t : 4 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 Figure 5: Plot of ln[y 2 (t)] A comparison of y(t) and ^y(t)isshown in Figure 6. As can be seen the t is prettygood. Wenow use some of the MATLAB commands to nish up the analysis. If we use the MATLAB statments K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,0.4,449);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr4yd.eps wegettheresponses shown in Figure 7, and the t is still prettygood, which means our estimates of B and C are reasonable. The last taskisto ndK.Wehave K p 1 p 2 =4:2719;; whichinthe presentcasebecomes K =4:271928:57 40 = 4882:1: 5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Figure 6: Plot of ^y(t)versusy(t) However, Figure 8 shows the armature voltage for the recorded step re- sponse. Thus, we need to divide K by32toachieve the correct gain. Finally G(s)= 135:6147 (s +28:57)(s++40) : 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 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 mtr4step t=mtr4step(1:500,1);; y=mtr4step(1:500,2);; plot(t,y) print -deps sr883mtr4ya.eps t=t+.002;; y=y(49:500);; t=t(49:500);; y(1) = 0;; y(2) = 0.13 t(1) = 0;; plot(t,y) print -deps sr883mtr4yb.eps maxy = max(y) y1 = (max(y) + 0.0001) -y;; w=log(y1);; plot(t(1:300),w(1:300));; 7 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 5 10 15 20 25 30 35 40 Figure 8: Recorded armature voltage print -deps 883mtr4ylny1.eps s=-5/0.175 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 883mtr4yB.eps B=mean(Bv(90:140)) A=mean(y(250:450)) y2 = y-A +B*exp(s*t);; w2 = log(y2);; plot(t(1:150),w2(1:150)) print -deps 883mtr4ylny2.eps C=-(A + B) s1 = -40 yhat = A + B*exp(s*t) + C*exp(s1*t);; plot(t,y,'k-',t,yhat,'k--') print -deps sr883mtr4yc.eps K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) 8 T=linspace(0,0.4,449);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr4yd.eps load mtr4av tav = mtr4av(1:500,1);; av = mtr4av(1:500,2);; plot(tav,av) print -deps 883mtr4yav.eps 9