-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 0.5 1 1.5 2 2.5 3 3.5 Figure 1: Rawstep response data Solution 8.8.15 The rst step is to get the data into MATLAB. The MATLAB statements load mtr6stepc t=mtr6stepc(1:500,1);; y=mtr6stepc(1:500,2);; y=y-y(99);; t=t-t(99);; y=y(99:500);; t=t(99:500);; y(1) = 0;; t(1) = 0;; load the data and makesome adjustments so that the response is zero when t =0.Theunadjusted and adjusted time responses are shown in Figures 1 and 2, respectively. The MATLAB statements y(150:402) = mean(y(150:402));; 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 3 3.5 Figure 2: Adjsuted step response data y(145:149) = mean(y(145:149));; y(140:144) = mean(y(140:144));; y(135:139) = mean(y(135:139));; y(130:134) = mean(y(130:134));; y(125:129) = mean(y(125:129));; y(120:124) = mean(y(120:124));; y(115:119) = mean(y(115:119));; y(110:114) = mean(y(110:114));; y(105:109) = mean(y(105:109));; y(100:104) = mean(y(100:104));; y(95:99) = mean(y(95:99));; y(90:94) = mean(y(90:94));; y(85:89) = mean(y(85:89));; y(80:84) = mean(y(80:84));; y(75:79) = mean(y(75:79));; y(70:74) = mean(y(70:74));; y(65:69) = mean(y(65:69));; y(60:64) = mean(y(60:64));; 2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 3 3.5 Figure 3: Smoothed step response y(55:59) = mean(y(55:59));; y(50:54) = mean(y(50:54));; plot(t,y) print -deps sr8815a.eps representacrude smoothing that results in the step response shown in Fig- ure 3. We are nowinaposition to apply the identi cation techniques discussed in Chapter 3. As a rst step, the MATLAB statments y76 =y(76) y86 = y(86) y96 = y(96) y106 = y(106) A=mean(y(150:402));; yy= [y(76) y(86) y(96) y(106)];; tt = [t(76) t(86) t(96) t(106)];; plot(t,y,'b.',tt,yy,'bd') print -deps sr8815d.eps 3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 3 3.5 Figure 4: Smoothed step response with selected data points superimposed produces the response shown in Figure 4. Since the data(set) is noisy,we use the MATLAB statments y(96) = y(96) + 0.06;; y(106) = y(106)-0.03;; yy= [y(76) y(86) y(96) y(106)];; plot(t,y,'b.',tt,yy,'bd') print -deps sr8815c.eps to makesome adjustnments an produce the plot shown in Figure 5. Wenow use the MATLAB statements s=-3 B30 = (yy - max(y))./exp(s*tt) s=-3.1 B31 = (yy - max(y))./exp(s*tt) s=-3.2 B32 = (yy - max(y))./exp(s*tt) s=-3.3 B33 = (yy - max(y))./exp(s*tt) 4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 3 3.5 Figure 5: Smoothed step response with corrected data points superimposed Btab = [B30;;B31;;B32;;B33] to calculate the required values of B whichwestoreinthe array\Btab." Those values are: Btab = -3.77166038351623 -3.08530856324902 -2.84487754983542 -2.82324801283584 -4.06541299056572 -3.35902812243020 -3.12839482070992 -3.13580998697516 -4.38204427315173 -3.65703128097996 -3.44016710132597 -3.48297571793418 -4.72333611773837 -3.98147240886751 -3.78301025391670 -3.86857618991797 A brief perusal of these numbers shows that are best estimate of the rst pole location is essentially the same as that found in the chapter, p 1 =3:2 or 3.3 with just a slightbiastoward 3.3. We use the MATLAB statements B32m = mean(B32(2:4)) B33m = mean(B33(2:4)) B=B33m s=-3.3 y2 = y-(max(y) + B*exp(s*t));; 5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 Figure 6: Determining p 2 w=log(y2);; plot(t(1:50),w(1:50)) print -deps 8815lny2.eps to generate the plot shown in Figure 6 Note that 3.25656 was increased to 3.2566 to avoid ln(0). The initial slope is roughly ;14:28. We then use the MATLAB statments T=linspace(0,4,402);; Kplant = A*s*s1 g=zpk([],[s s1],Kplant) [yhat T] = step(g,T);; plot(t,y,'k-',T,yhat,'k--') print -deps sr8815d.eps to compare the recorded response to that of our model whichis G(s)= 142:299 (s +3:3)(s+14:28) 6 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1.5 2 2.5 3 3.5 Figure 7: Plot of ^y(t)versusy(t) The comparison is shown in Figure 7. As can be seen the t is very good. The gain is determined from the expression K p 1 p 2 = A;; whichinthe presentcasebecomes K =3:1263:314:28 = 142:299: However, Figure 8 shows the armature voltage for the recorded step re- sponse. Thus, we need to divide K by27toachieve the correct gain. Finally G(s)= 5:456 (s +3:3)(s++14:28) : This is essentially the result wegotinthis chapter, so wehavenotimproved the model very much, if at all. 7 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 5 10 15 20 25 30 Figure 8: Recorded armature voltage The entire MATLAB program used to obtain this transfer function is % %load the data % load mtr6stepc t=mtr6stepc(1:500,1);; y=mtr6stepc(1:500,2);; plot(t,y) print -deps sr8815a.eps y=y-y(99);; t=t-t(99);; y=y(99:500);; t=t(99:500);; y(1) = 0;; t(1) = 0;; plot(t,y) print -deps sr8815b.eps y(150:402) = mean(y(150:402));; 8 y(145:149) = mean(y(145:149));; y(140:144) = mean(y(140:144));; y(135:139) = mean(y(135:139));; y(130:134) = mean(y(130:134));; y(125:129) = mean(y(125:129));; y(120:124) = mean(y(120:124));; y(115:119) = mean(y(115:119));; y(110:114) = mean(y(110:114));; y(105:109) = mean(y(105:109));; y(100:104) = mean(y(100:104));; y(95:99) = mean(y(95:99));; y(90:94) = mean(y(90:94));; y(85:89) = mean(y(85:89));; y(80:84) = mean(y(80:84));; y(75:79) = mean(y(75:79));; y(70:74) = mean(y(70:74));; y(65:69) = mean(y(65:69));; y(60:64) = mean(y(60:64));; y(55:59) = mean(y(55:59));; y(50:54) = mean(y(50:54));; plot(t,y) print -deps sr8815c.eps y76 =y(76) y86 = y(86) y96 = y(96) y106 = y(106) A=mean(y(150:402)) yy= [y(76) y(86) y(96) y(106)];; tt = [t(76) t(86) t(96) t(106)];; plot(t,y,'b.',tt,yy,'bd') print -deps sr8815d.eps y(96) = y(96) + 0.06;; y(106) = y(106)-0.03;; yy= [y(76) y(86) y(96) y(106)];; plot(t,y,'b.',tt,yy,'bd') print -deps sr8815c.eps s=-3 B30 = (yy - max(y))./exp(s*tt) s=-3.1 B31 = (yy - max(y))./exp(s*tt) 9 s=-3.2 B32 = (yy - max(y))./exp(s*tt) s=-3.3 B33 = (yy - max(y))./exp(s*tt) Btab = [B30;;B31;;B32;;B33] B32m = mean(B32(2:4)) B33m = mean(B33(2:4)) B=B33m s=-3.3 y2 = y-(max(y) + B*exp(s*t));; w=log(y2);; plot(t(1:50),w(1:50)) print -deps 8815lny2.eps s1 = -14.28 T=linspace(0,4,402);; Kplant = A*s*s1 g=zpk([],[s s1],Kplant) [yhat T] = step(g,T);; plot(t,y,'k-',T,yhat,'k--') print -deps sr8815d.eps K=Kplant/27 load mtr6av tav = mtr6av(1:500,1);; av = mtr6av(1:500,2);; plot(tav,av) print -deps 881av.eps 10