pidsim1.m

状态反馈matlab仿真程序

clf,
num=1,den=[1 3 2],tao=1,
w=logspace(-2,2,200)',
pause;
% plot Nyquist graph
clf;
[r,i]=nyquist(num,den,w),
rd=real((r+j*i).*exp(-j*w*tao)),
id=imag((r+j*i).*exp(-j*w*tao)),
plot(rd,id,'--'),
pause;
% slove the magnitude and phase
%[mag,phase]=nichols(num,den,w),
%mag,phase,
%pause;
%magdb=20*log10(mag),
%[gm,pm,wpc,wgc]=margin(magdb,phase,w),grid,
%gm,pm,wpc,wgc,
%pause;
% System Identification
A=abs((rd(109,1)+rd(110,1))/2), w0=(w(109)+w(110))/2,
t=[0:0.1:10],
[y,x]=step(num,den,t);
K=y(100,1),
T=(sqrt(K*K/A/A-1))/w0,ti=(pi-atan(w0*T))/w0,
num1=K, den1=[T 1],
printsys(num1,den1),
[r1,i1]=nyquist(num1,den1,w);
%step(num1,den1),
rd1=real((r1+j*i1).*exp(-j*w*ti)),
id1=imag((r1+j*i1).*exp(-j*w*ti)),
hold on,
plot(rd1,id1,':'),
pause;
%Simulition by SIMULINK
%k2=(1.41*0.6174*T-T-ti)/(K*ti),
%k1=(0.6174^2*T/ti-1)/K-k2,
%Kp=(2.15*T*ti-K*k1-K*k2-1)/K,
%Td=(1.75*T*ti-T-ti-K*k2*ti)/(K*Kp),
%Kd=Kp*Td,
%Ti=(K*Kp*ti^2)/(T*0.6174^3),
%Ki=Kp/Ti,
num2=0.3099*0.6174^2,
den2=[1 1.41*0.6174 0.6174^2],
[r2,d2]=nyquist(num2,den2,w),
plot(r2,d2),grid,