iden_c21i.m

状态反馈matlab仿真程序

%-----------------------------------------------------------------------
%    Model Identification Using Relay Feed-back with delay width 'h'
%-----------------------------------------------------------------------
% Note: This program is available only after the Chap21.mdl is executed 
%                      which output the 'y' value

% Give the value according to the parameters in Chap21e.mdl
clear all; clf;
 
num1=[1];den1=[1 3 2];tao1=10;
w=logspace(-2,2,200)';
[r1,i1]=nyquist(num1,den1,w);
tt=[0:0.1:50];
z=step(num1,den1,tt);
K=z(length(z),1);

rf=10;
L=rf/K;
h=0.001;
d=0.05*rf+h*1.01;
Tal=20*tao1+35;

P=1; Pmax=1;

while (P<=Pmax)
   
% Calculate the Bode Graph of the Real Object
rd1=real((r1+j*i1).*exp(-j*w*tao1));
id1=imag((r1+j*i1).*exp(-j*w*tao1));

% Calculate the relay system output by sim the model of 'chap21e.mdl'
[t,xx,y]=sim('chap22d',Tal);
y=y-rf;
dm=max(y)/2;

% Identification using relay feedback output response
m=0;
dy=y(2:length(y))-y(1:(length(y)-1));
for i=2:1:(length(y)-1)
   if (y(i-1)*y(i)<0 & abs(dy(i))>0) 
      m=m+1;
      tm(m)=i;
   end
end
m0=1;
tm0(m0)=tm(1);
for (i=2:1:m)
   if ((tm(i)-tm(i-1))>15)
      m0=m0+1;
      tm0(m0)=tm(i);
   end
end
m=m0;
tm=tm0;

T0=Tal*(tm(m)-tm(m-2))/(length(y)-1);  % Tal must be same as the time in Chap21
w01=2*pi/T0;
a=(max(y(tm(m-2):tm(m)))-min(y(tm(m-2):tm(m))))/2;
A1=(pi*a)/(4*d);            % 'd' is variable

w0(P)=w01;
A(P)=A1;
T=sqrt((K/A(P))^2-1)/w0(P),
tao2=(pi-atan(w0(P)*T))/w0(P),
num2=K;
den2=[T 1];
[r2,i2]=nyquist(num2,den2,w);
rd2=real((r2+j*i2).*exp(-j*w*tao2));
id2=imag((r2+j*i2).*exp(-j*w*tao2));

% To correct the value of 'A' and 'w0'
dw=0.981/((-0.078*h+1)^0.1);    
w0(P)=w01*dw;
dA=3.2287/(tao1^(-0.6667)+2.4567);
A(P)=A1*dA;

T=sqrt((K/A(P))^2-1)/w0(P),
tao3=(pi-atan(w0(P)*T))/w0(P),
num3=K;
den3=[T 1];

%Calculate the Bode Graph of the Identified Model
[r3,i3]=nyquist(num3,den3,w);
rd3=real((r3+j*i3).*exp(-j*w*tao3));
id3=imag((r3+j*i3).*exp(-j*w*tao3));

P=P+1;

end

% Calculate the coefficient fo the state feed back
N=1;
l0=1;
for (i=1:1:N)
   l(i)=tao2^i/fac(i);
end

% k(1)=-1.0 ~ 1.5
% k(2)=-2 ~ 100
figure(1)
k=[-1 -0.1];
for (i=1:1:5)
   [t,xx,y]=sim('t21a1',50);
   plot(t,y),
   hold on
   k(1)=k(1)+0.5;
end

figure(2)
k=[0.1 -1];
for (i=1:1:5)
   [t,xx,y]=sim('t21a1',50);
   plot(t,y),
   hold on
   k(2)=k(2)+1;
end

k=[-0.1344 2.0343];
b2=l(1)*T;
b1=l(1)+K*k(2)*l(1)+T*l0;
b0=l0+K*k(2)*l0+K*k(1);
num4=[K/b2];
den4=[1 b1/b2 b0/b2];
printsys(num4,den4);
roots(den4);
[rd4,id4]=nyquist(num4,den4,w);

% Plot the Nyquist graph
figure(3)
xlabel('Re'); ylabel('Im'); Title('Nyquist Graph');
plot(rd1(1:(length(rd1)-91)),id1(1:(length(rd1)-91)),'r-'),
hold on
plot(rd2(1:(length(rd2)-100)),id2(1:(length(rd2)-100)),'b--'),
plot(rd3(1:(length(rd3)-91)),id3(1:(length(rd3)-91)),'b:'),
plot(rd4(1:(length(rd4)-71)),id4(1:(length(rd4)-71))),