iden_c21h.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;

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);

h=6;
d=h/K+0.05;
Tal=5*tao1+5; 

P=1; Pmax=1;

while (P<=Pmax)
   
% 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(1:(length(rd1)-91)),id1(1:(length(rd1)-91)),'r-'),
xlabel('Re'); ylabel('Im'); Title('Nyquist Graph');

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

% 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
w01=2*pi/T0;
a=abs(ym(n));
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];
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(1:(length(rd2)-100)),id2(1:(length(rd2)-100)),'b--'),

% 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.4367);
A(P)=A1*dA;

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(1:(length(rd2)-91)),id2(1:(length(rd2)-91)),'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' 
%h=h+0.1;
P=P+1;

end

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