-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 1: Rawstep response data Solution 8.8.3mtr5cy The rst step is to get the data into MATLAB. The MATLAB statements load mtr5stepc t=mtr5stepc(1:500,1);; y=mtr5stepc(1:500,2);; plot(t,y) print -deps sr883mtr5cya.eps 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+0.06;; y=y(45:500);; t=t(45:500);; y(1) = 0;; y(2) = 0.05;; 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 2: Adjusted step response y(3) = 0.1;; y(4) = 0.15;; y(5) = 0.2;; t(1) = 0;; plot(t,y) print -deps sr883mtr5cyb.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 883mtr5cylny1.eps produces the response shown in Figure 3. Note that the maximum resposne was increased to0.001 to avoid ln(0). 2 0 0.5 1 1.5 2 2.5 3 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3: Plot of ln(1;y(t) The slope is roughly ;3:85=1:90625 = 2:01976: So there is a pole around s = ;2:0. Our next goal is to nd B using the MATLAB statements: s=-2.01967 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 883mtr5cyB.eps B=mean(Bv(55:120)) 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 =3:7121 and B = ;5:3425: 3 0 0.5 1 1.5 -8 -7.5 -7 -6.5 -6 -5.5 -5 -4.5 -4 -3.5 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);(3:7121; 5:3425e ;2:012t : 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 -2.0. We conclude 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)=;(3:7121;5:3425) = 1:6304: and place the second pole at s = ;10, on the assumption that it is, relatively much farther to the left in the s plane. Weknowthe pole is farther out, so 4 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Figure 5: Plot of ln[y 2 (t)] we start with this guess. Then our estimate of y(t)is ^y(t)=3:7121; 5:3425e ;2:012t +1:6304e ;10t : The MATLAB statements K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,4.5,455);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr5cyd.eps generate the comparison of y(t)and^y(t) shown in Figure 6. As can be seen the estimated response is o near t =0.However, the overall t is not that bad. The culprits are probably our estimates of B and C.Rather than try to re ne our estimates of B and C,whichdon't really interest us anyway, weuse some of the conveientMATLAB commands to nish up the analysis. If we use the MATLAB statments K=A*abs(s)*abs(s1) 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 6: Plot of ^y(t)versusy(t) g=zpk([],[s s1],K) T=linspace(0,4.5,455);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr5cyd.eps we get the responses shown in Figure 7, and the t is very prettygood.If wemove the pole to s = ;2:2, weget the response in Figure 8 The last task is to nd K.Wehave K p 1 p 2 =4:0696;; whichinthe presentcasebecomes K =3:7121 2:210 = 59:9778: However, Figure 9 shows the armature voltage for the recorded step re- sponse. Thus, we need to divide K by32toachieve the correct gain. Finally G(s)= 2:5521 (s +2:2)(s+10) : 6 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 7: Plot of ^y(t)versusy(t) load mtr5stepc t=mtr5stepc(1:500,1);; y=mtr5stepc(1:500,2);; plot(t,y) print -deps sr883mtr5cya.eps t=t+0.06;; y=y(45:500);; t=t(45:500);; y(1) = 0;; y(2) = 0.05;; y(3) = 0.1;; y(4) = 0.15;; y(5) = 0.2;; t(1) = 0;; plot(t,y) print -deps sr883mtr5cyb.eps maxy = max(y);; y1 = (max(y) + 0.0001) -y;; 7 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 8: Plot of ^y(t)versusy(t)fors = ;2:2 w=log(y1);; plot(t(1:275),w(1:275));; print -deps 883mtr5cylny1.eps s=-2.01967 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 883mtr5cyB.eps B=mean(Bv(55:120)) A=mean(y(250:455)) y2 = y-A +B*exp(s*t);; w2 = log(y2);; plot(t(1:150),w2(1:150)) print -deps 883mtr5cylny2.eps C=-(A + B) s1 = -8 yhat = A + B*exp(s*t) + C*exp(s1*t);; plot(t,y,'k-',t,yhat,'k--') print -deps sr883mtr5cyc.eps 8 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 10 15 20 25 30 35 Figure 9: Recorded armature voltage pause K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,4.5,455);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr5cyd.eps s=-2.2 s1 = -10 K=A*abs(s)*abs(s1) g=zpk([],[s s1],K) T=linspace(0,4.5,455);; [yhat1 T] = step(g,T);; plot(t,y,'k-',T,yhat1,'k--') print -deps sr883mtr5cye.eps pause load mtr5avc tav = mtr5avc(1:500,1);; 9 av = mtr5avc(1:500,2);; plot(tav,av) print -deps 883mtr5cyav.eps 10