Solution 8.8.3.mtr1cy The full MATLAB program used to nd the transfer function is: load mtr1stepc t=mtr1stepc(1:500,1);; y=mtr1stepc(1:500,2);; plot(t,y) print -deps sr883mtr1cya,eps y=y(51:500);; t=t(51:500);; y(1) = 0;; t(1) = 0;; plot(t,y) print -deps sr883mtr1cyb.eps maxy = max(y) y1 = (max(y) + 0.0001) -y;; w=log(y1);; plot(t(1:200),w(1:200));; print -deps 883mtr1clny1.eps s1 = -1.5/2 s2 = -4.2/1.75 s=s2 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 8813mtr1cyB.eps B=mean(Bv(75:100)) A=mean(y(200:450)) y2 = y-A +B*exp(s*t);; w2 = log(y2);; plot(t(1:150),w2(1:150)) print -deps 883mtr1cyln2.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 sr883mtr1cyc.eps yhat1 =A+B*exp(s*t) + C*exp(-2*t);; plot(t,y,'k-',t,yhat1,'k--') print -deps sr883mtr1cyd.eps 1 yhat2 =A+B*exp(s*t) + C*exp(-4*t);; plot(t,y,'k-',t,yhat2,'k--') print -deps sr883mtr1cye.eps K=max(y)*abs(s)*abs(s1) load mtr1av tav = mtr1av(1:500,1);; av = mtr1av(1:500,2);; plot(tav,av) print -deps 883mtr1yav.eps pp1 = 2.4 pp2 = 1.75 KK= A*pp1*pp2 Kplant = KK/24 g=zpk([],[-pp1 -pp2],KK);; [yhat3, T] = step(g,t);; plot(t,y,'k-',T,yhat3,'k--') print -deps sr883mtr1cyf.eps The rst step is to get the data into MATLAB. That is done by the MAT- LAB statements: load mtr1stepc t=mtr1stepc(1:500,1);; y=mtr1stepc(1:500,2);; plot(t,y) print -deps sr883mtr1cya,eps The plot of the raw date is shown in Figure1. As can be seen, the plot does not start at t =0and the initial voltage is not zero. The MATLAB statements: y=y(51:500);; t=t(51:500);; y(1) = 0;; t(1) = 0;; plot(t,y) print -deps sr883mtr1cyb.eps seperate the data into a column vector of times and a column vector of responses, and adjusts the plot so that it begins with a voltage of zero at t =0.The adjusted plot is shown in Figure2 We are nowinaposition to apply the identi cation techniques discussed in Chapter 3. As a rst step, the MATLAB statements 2 -1 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 Figure 1: Rawstep response data 3 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 Figure 2: Adjusted step response 4 0 0.5 1 1.5 2 2.5 3 3.5 4 -5 -4 -3 -2 -1 0 1 2 Figure 3: Plot of ln(1 ; y(t) maxy = max(y) y1 = (max(y) + 0.0001) -y;; w=log(y1);; plot(t(1:200),w(1:200));; print -deps 883mtr1clny1.eps produces the response shown in Figure 3. Note that the maximum value of y =3:0917 was increased to 0.0001 to avoid ln(0). Wesee that there are two distinct slopes s 1 = ;1:5 2 = ;0:75;; and s 2 = ;3:2 1:75 = ;2:4: Wechoose s 2 because we try to measure for the slowmodeastheresponse approaches steady state. Thus weassume there is a pole near s = ;2:4. Our next goal is to nd B using the MATLAB statements: s1 = -1.5/2 5 0 0.5 1 1.5 2 2.5 3 -300 -250 -200 -150 -100 -50 0 Figure 4: Comparison of adjusted response and y 1 (t)=1; e ;24:9t s2 = -4.2/1.75 s=s2 Bv = (y - max(y))./exp(s*t);; plot(t(1:150),Bv(1:150)) print -deps 8813mtr1cyB.eps B=mean(Bv(75:100)) Note that wehave used the MATLAB command: ./ to do an elementbyelement division using two matrices. The plotofB versus time is shown in Figure 4. Wesee that it is going to be dicult to estimate B using our strategy of looking for a regin where B is reasonably constant. A =5:0417 and B = ;56:07;; and weknow that jBj > jAj 6 0 0.5 1 1.5 2 2.5 3 -4 -3 -2 -1 0 1 2 3 4 5 Figure 5: Plot of ln[y 2 (t)] 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) ; (5:0417; 56:07e ;p 1 t ): We then take the natural logarithm of y 2 and plot it versus time. as shown in Figure 5. The inital slope is on the order of ;2, so wereally can't gain any information about the second pole, if we believe the poles arewidely separated. The other possibilityisthat the poles are not widely seperated. Wewill address that issue as well. If we assume the poles are widely seperated, then will just havetoguess at the location of the second pole. We can nd C = ;(A + B)=;(325656; 3:8992) = 51:025: Then our estimate of y(t)is ^y(t)=5:0417; 56:07e ;2:4t +51:025e ;10t ;; 7 0 1 2 3 4 5 6 7 8 9 -25 -20 -15 -10 -5 0 5 10 Figure 6: Plot of ^y(t)versusy(t) where wehave simply guessed at the second pole location. In the case of the motor whose transfer function weare trying to nd, we believe that the poles are widely separated. Acomparison of y(t) and ^y(t)isshown in Figure 6. As can be seen the t is very bad. Placing the second pole at s = ;2ors = ;4does not help, as shown byFigure 7 and Figure 7. Wehavetoconclude that our value of B is wayo .Wenowcaneither try to restimate B,whichdoesn't look too promising, or use the pole at s = ;2:4 that wearefairly sure of and use MATLAB to nd the other. After some experimentation, we arriveat G(s)= 0:8823 (s +1:75)(s+2:4) : As can be seen, from Figure 9, the t is is better prettygood.The MATLAB code to nd this transfer function is pp1 = 2.4 pp2 = 1.75 KK= A*pp1*pp2 8 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 Figure 7: Trial system for s 1 = ;2:4, s 2 = ;2, B =56,C =51 9 0 1 2 3 4 5 6 7 8 9 -8 -6 -4 -2 0 2 4 6 Figure 8: Trial system for s 1 = ;2:4, s 2 = ;4, B =56,C =51 Kplant = KK/24 g=zpk([],[-pp1 -pp2],KK);; [yhat3, T] = step(g,t);; plot(t,y,'k-',T,yhat3,'k--') print -deps sr883mtr1cyf.eps Note that wehaveset K = A  2:4 1:75 24 =0:8823 to accountfor the roughly 24 V that is applied to the system. 10 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 Figure 9: Step response for alternativemodelwith poles at s = ;1:75 and ;2:4 11