iden_c21e.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
clf; clear all;
d=10;
h=0.15;
dm=d/4;

w=logspace(-2,2,200)';
num1=1;den1=[1 3 2];
[r1,i1]=nyquist(num1,den1,w);
tt=[0:0.1:10];
z=step(num1,den1,tt);
K=z(100,1);

Tal=20; tao1=1,
P=1; Pmax=50;

while (P<=Pmax)
   
tao1, 
% Plot 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));
for i=2:length(id1)
   if (id1(i-1)*id1(i)<0 & rd1(i)<0)
      i0=i;
      break;
   end
end
A0(P)=-rd1(i0);
wc0(P)=w(i0);
%plot(rd1,id1,'r-'),
%xlabel('Re'); ylabel('Im'); Title('Nyquist Graph');

% Calculate the relay system output by sim the model of 'chap21e.mdl'
[t,xx,y]=sim('chap21e',Tal);

% Identification using relay feedback output response
m=0; n=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
   if (dy(i-1)*dy(i)<0 & abs(y(i))>dm)  % 'dm' is cariable with 'd'
      n=n+1;
      ym(n)=y(i);
   end
end

T0=Tal*(tm(m)-tm(m-2))/(length(y)-1);  % Tal must be same as the time in Chap21
w0(P)=2*pi/T0;
a=abs(ym(n));
A1=(pi*a)/(4*d);            % 'd' is variable
%dh=(pi-asin(h/a))/pi;    % To correct the value of 'A'
%A=A1/dh;
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];

% Plot the Bode Graph of the Identified Model
%hold on;
%[r2,i2]=nyquist(num2,den2,w);
%rd2=real((r2+j*i2).*exp(-j*w*tao2));
%id2=imag((r2+j*i2).*exp(-j*w*tao2));
%plot(rd2,id2,'b:'),

% Plot the graph of the output value
%figure(2);
plot(t,y); grid,
xlabel('t'); ylabel('y'); Title('Output Oscillation Curve');

% Calculate the error of A and w0
err_A(P)=abs((A(P)-A0(P))/A0(P));
err_w0(P)=(w0(P)-wc0(P))/wc0(P);
beta_A(P)=A0(P)/A(P);
beta_w0(P)=wc0(P)/w0(P);

% Revise Time-delay Constant 'tao1' and Sim time 'Tal' 
tao1=tao1+1;
Tal=Tal+8;
P=P+1;

end

plot(t,y); grid,
xlabel('t'); ylabel('y'); Title('Output Oscillation Curve');

% Plot Error Curve of 'A' and 'w0'
figure(2);
plot(1:1:Pmax,err_A,'b*',1:1:Pmax,err_w0,'r.');
xlabel('Time delay constant'); ylabel('Err');
Title('The Error Curves of A and w0');
figure(3);
plot(1:1:Pmax,beta_A,'b*',1:1:Pmax,beta_w0,'r.');
xlabel('Time delay constant'); ylabel('Revising Coefficient');
Title('Revising Coefficient of A and w0');