Solution 8.8.9.a
The full MATLAB program used to nd the transfer function is:
%
%Find the length of the array y(t) Note:this the step response. It has to be
%pulled into Matlab before this program is run.
%
load p3id1;;
y=p3id1;;
topm = size(y);;
top = topm(1,1);;
%
%Initialize counter.
%
asum = 0.0;;
%
%Puu data into arrays
%
t=p3id1(1:top,1);;
y=p3id1(1:top,2);;
plot(t,y)
pause
print -deps sr889a.a.eps
%
%
%
k0 = y(top);;
%
f0 = k0-y;;
%
%Reset counter
%
k=2;;
%
%Integrate f0 = k0 - y(t). and form the function y1(t).
%
while ( k <= top)
l= k-1;;
asum = asum + ( ( (f0(k)+f0(l))/2 ) *( t(k) - t(l) ) );;
1
y1(k) = asum;;
k=k+1;;
end
%
%Save data for plotting.
%
save time.dat t -ascii
save y.dat y -ascii
save y1.dat y1 -ascii
save f0.dat f0 -ascii
%
%SetK1 = to max(y1(t)), that is its final steady state value.
%
k1 = asum;;
%
%form function:f1 = k1 - y1(t).
%
f1 = k1-y1;;
%
%Save function for later plotting
%
save f1.dat f1 -ascii
%
%reset summer
%
asum = 0.0;;
%
%Reset counter
%
k=2;;
%
%Integrate f1(k) = k1 - y1(t) to get y2(t)
%
while ( k <= top)
l= k-1;;
asum = asum + ( ( (f1(k)+f1(l))/2 ) *( t(k) - t(l) ) );;
y2(k) = asum;;
k=k+1;;
end
%
2
%Setk2 equal to maximum value of y2(t).
%
k2 = asum;;
%
%Save function y2(t) for later plotting.
%
save y2.data y2 -ascii
%
%
%form function: f2 = k2 - y2(t).
%
f2 = k2-y2;;
%
%Save function for later plotting
%
save f2.data f2 -ascii
%
%reset summer
%
asum = 0.0;;
%
%Reset counter
%
k=2;;
%
%Integrate f2 = k2 - y2(t) to get y3(t)
%
while ( k <= top)
l= k-1;;
asum = asum + ( ( (f2(k)+f2(l))/2 ) *( t(k) - t(l) ) );;
y3(k) = asum;;
k=k+1;;
end
%
%SetK3 equal to maximum value of y3(t).
%
k3 = asum;;
A=[k0,0,0;;k1,-k0,0;;k2,-k1,k0]
b=[k1;;k2;;k3]
x=inv(A)*b
3
0 1 2 3 4 5 6 7 8 9
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 1: Recorded step response
a1 = x(1)
a2 = x(2)
a3 = x(3)
polyd=[a3 a2 a1 1]
roots(polyd)
a1 = k1/k0
a2 = -(k2-a1*k1)/k0
a3 = (k3 - a1*k2 + a2*k1)/k0
polyd = [a3a2a11]
p=roots(polyd)
T=linspace(0,4,200);;
g=zpk([],[p(1) p(2) p(3)],k0/a3)
[yhat,T] = step(g,T);;
plot(t,y,'k-',T,yhat,'k--')
print -deps sr889a.b.eps
The step response is shownin Figure 1. The step response is verysmooth.
The program, when executed, yields
EDU>sm889a
4
A=
1.0000 0 0
1.1329 -1.0000 0
1.1443 -1.1329 1.0000
b=
1.1329
1.1443
1.1338
x=
1.1329
0.1391
-0.0050
a1 =
1.1329
a2 =
0.1391
a3 =
-0.0050
polyd =
-0.0050 0.1391 1.1329 1.0000
5
ans =
34.7021
-5.7178
-1.0134
a1 =
1.1329
a2 =
0.1391
a3 =
-0.0050
polyd =
-0.0050 0.1391 1.1329 1.0000
p=
34.7021
-5.7178
-1.0134
Zero/pole/gain:
-201.0743
----------------------------
(s-34.7) (s+5.718) (s+1.013)
6
EDU>
The algorithm does not converge. The only hope would be to try a more
accurate integration scheme than trapezoidal integration.
7