％仿真程序1
%big delay plant with conventional pid controller
clear all;
close all;
ts=20;
%delay plant0
sys=tf([1],[20,1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');
%delay plant1
sys1=tf([1],[60,1],'inputdelay',80);
dsys1=c2d(sys1,ts,'zoh');
[num1,den1]=tfdata(dsys1,'v');
%delay plant2
sys2=tf([1],[160,1],'inputdelay',80);
dsys2=c2d(sys2,ts,'zoh');
[num2,den2]=tfdata(dsys2,'v');
u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
y_1=0;y_2=0;y_3=0;
error_1=0;error_2=0;
ei=0;
u_11=0;u_21=0;u_31=0;u_41=0;u_51=0;
y_11=0;y_21=0;y_31=0;
error_11=0;error_21=0;
ei1=0;u_12=0;u_22=0;u_32=0;u_42=0;u_52=0;
y_12=0;y_22=0;y_32=0;
error_12=0;error_22=0;
ei2=0;

for k=1:1:200
time(k)=k*ts;
%delay plant
yout(k)=-den(2)*y_1+num(2)*u_5;   
yout1(k)=-den1(2)*y_11+num1(2)*u_51;   
yout2(k)=-den2(2)*y_12+num2(2)*u_52;    
rin(k)=1;
error(k)=rin(k)-yout(k);
ei=ei+error(k)*ts;
error1(k)=rin(k)-yout1(k);
ei1=ei1+error1(k)*ts;
error2(k)=rin(k)-yout2(k);
ei2=ei2+error2(k)*ts;

kp1=0.80;ki1=0.005;kd1=3.0;kp=0.59;ki=0.0052;kd=3.0;
u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+ki*ei;
u1(k)=kp1*error1(k)+kd1*(error1(k)-error_11)/ts+ki1*ei1;
u2(k)=kp1*error2(k)+kd1*(error2(k)-error_12)/ts+ki1*ei2;

if u(k)>=10
       u(k)=10;
end
if u(k)<=-10
       u(k)=-10;
end
if u1(k)>=10
       u1(k)=10;
end
if u1(k)<=-10
       u1(k)=-10;
end
if u2(k)>=10
       u2(k)=10;
end
if u2(k)<=-10
       u2(k)=-10;
end
u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=yout(k);
u_51=u_41;u_41=u_31;u_31=u_21;u_21=u_11;u_11=u1(k);
y_31=y_21;y_21=y_11;y_11=yout1(k);
u_52=u_42;u_42=u_32;u_32=u_22;u_22=u_12;u_12=u2(k);
y_32=y_22;y_22=y_12;y_12=yout2(k);
error_2=error_1;error_1=error(k);error_21=error_11;error_11=error1(k);
error_22=error_12;error_12=error2(k);
end

figure(1);
plot(time,rin,'b',time,yout,'r',time,yout1,'g',time,yout2,':');grid;
xlabel('time(s)');ylabel('yout&yout1&yout2');

figure(2);
plot(time,u,'r',time,u1,'g',time,u2,':');grid;
xlabel('time(s)');ylabel('u&u1&u2');




％仿真程序2
% big delay plant with integration separation pid controller
clear all;
close all;
ts=20;
%delay plant
sys=tf([1],[60,1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');

u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
y_1=0;y_2=0;y_3=0;
error_1=0;error_2=0;
ei=0;
u_11=0;u_21=0;u_31=0;u_41=0;u_51=0;
y_11=0;y_21=0;y_31=0;
error_11=0;error_21=0;
ei1=0;

for k=1:1:200
time(k)=k*ts;
%delay plant
yout(k)=-den(2)*y_1+num(2)*u_5; 
yout1(k)=-den(2)*y_11+num(2)*u_51;   
rin(k)=1;
error(k)=rin(k)-yout(k);
ei=ei+error(k)*ts;
error1(k)=rin(k)-yout1(k);
ei1=ei1+error1(k)*ts;
if abs(error1(k))>=0.8&abs(error1(k))<=1
        beta=0.35;
elseif abs(error1(k))>=0.5&abs(error1(k))<=0.8
        beta=0.65; 
elseif abs(error1(k))>=0.3abs(error1(k))<=0.5 
        beta=0.95;
else 
        beta=1.0;
end
kp=0.80;
ki=0.005
kd=3.0;
u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+ki*ei;
u1(k)=kp*error1(k)+kd*(error1(k)-error_11)/ts+beta*ki*ei1;

if u(k)>=10
       u(k)=10;
end
if u(k)<=-10
       u(k)=-10;
end
if u1(k)>=10
       u1(k)=10;
end
if u1(k)<=-10
       u1(k)=-10;
end

u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);
y_3=y_2;y_2=y_1;y_1=yout(k);
u_51=u_41;u_41=u_31;u_31=u_21;u_21=u_11;u_11=u1(k);
y_31=y_21;y_21=y_11;y_11=yout1(k);
error_2=error_1;error_1=error(k);
error_21=error_11;error_11=error1(k);
end
figure(1);
plot(time,rin,'k',time,yout,'r',time,yout1,':');grid;
xlabel('time(s)');ylabel('rin,yout&yout1');
figure(2);
plot(time,u,'r',time,u1,':');grid;
xlabel('time(s)');ylabel('u&u1');


















％仿真程序3
%delay control with dalin algorithm 
clear all;
close all;
ts=20;
%plant
l=1;
if l==1
Kp=1;
Tp=60;
tol=80;
else if l==2
Kp=1*1.05;
Tp=60*1.05;
tol=80*1.10
%plant
sys1=tf([Kp],[Tp 1],'inputdelay',tol);
dsys1=c2d(sys1,ts,'zoh');
[num1,den1]=tfdata(dsys1,'v');
%ideal closed loop 
sys2=tf([1],[20 1],'inputdelay',80);
dsys2=c2d(sys2,ts,'zoh');
%design dalin controller
dsys=1/dsys1*dsys2/(1-dsys2);
[num,den]=tfdata(dsys,'v');

u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;u_6=0;
y_1=0;y_2=0;y_3=0;
error_1=0;error_2=0;
ei=0;
for k=1:1:200
    time(k)=k*ts;
    rin(k)=1;
    yout(k)=-den1(2)*y_1+num1(2)*u_5;
    error(k)=rin(k)-yout(k);
    
    m=1;
    if m==1;       u(k)=(-den(6)*u_1-den(7)*u_2-den(10)*u_5-den(11)*u_6+num(1)*error(k)+num(2)*error_1+num(3)*error_2)/den(5);
    elseif m==2
           ei=ei+error(k)*ts;
           u(k)=0.80*error(k)+3.0*(error(k)-error_1)/ts+0.005*ei;
    end
   u_6=u_5;u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;
   if k==100
   u_1=u(k)+1.0;
   else
   u_1=u(k);
   end
   y_3=y_2;y_2=y_1;y_1=yout(k);
   
   error_2=error_1;
   error_1=error(k);
end
plot(time,rin,'b',time,yout,'r');grid;
xlabel('time(s)');
ylabel('rin,yout');






























％仿真程序5
%big delay pid with smith algorithm
clear all;
close all;
ts=20;
%delay plant0
sys=tf([1],[20,1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');
%delay plant1
sys1=tf([1],[60,1],'inputdelay',80);
dsys1=c2d(sys1,ts,'zoh');
[num1,den1]=tfdata(dsys1,'v');
%delay plant2
sys2=tf([1],[60,1],'inputdelay',80);
dsys2=c2d(sys2,ts,'zoh');
[num2,den2]=tfdata(dsys2,'v');
u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; y_1=0;
u_11=0;u_21=0;u_31=0;u_41=0;u_51=0;y_11=0;
u_12=0;u_22=0;u_32=0;u_42=0;u_52=0;y_12=0;
e2=0;e2_1=0;ei=0;
e21=0;e2_11=0;ei1=0;
e22=0;e2_12=0;ei2=0;
xm_1=0; xm_11=0; xm_12=0;
for k=1:1:200
    time(k)=k*ts;
    rin(k)=1;
    xm(k)=-den(2)*xm_1+num(2)*u_1;
    yout(k)=-den(2)*y_1+num(2)*u_5;
    xm1(k)=-den1(2)*xm_11+num1(2)*u_11;
    yout1(k)=-den1(2)*y_11+num1(2)*u_51;
    xm2(k)=-den2(2)*xm_12+num2(2)*u_12;
    yout2(k)=-den2(2)*y_12+num2(2)*u_52;
     
e2(k)=rin(k)-xm(k);
ei=ei+e2(k)*ts;
u(k)=0.50*e2(k)+0.035*ei;
e21(k)=rin(k)-xm1(k);
ei1=ei1+e21(k)*ts;
u1(k)=0.50*e21(k)+0.010*ei1;
e22(k)=rin(k)-xm2(k);
ei2=ei2+e22(k)*ts;
u2(k)=0.50*e22(k)+0.0045*ei2;
xm_1=xm(k);xm_11=xm1(k);xm_12=xm2(k);
u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);
u_51=u_41;u_41=u_31;u_31=u_21;u_21=u_11;u_11=u1(k);
u_52=u_42;u_42=u_32;u_32=u_22;u_22=u_12;u_12=u2(k);
y_1=yout(k);y_11=yout1(k);y_12=yout2(k);
end
plot(time,rin,'b',time,yout,'r',time,yout1,'g',time,yout2,':');
grid;
xlabel('time(s)');
ylabel('rin,yout');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%big delay pid with smith algorithm
clear all;
close all;
ts=20;
%delay plant
sys=tf([1],[60 1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');
u_1=0.0;u_2=0.0;u_3=0.0;u_4=0.0;u_5=0.0;
e2=0; e2_1=0;ei=0;xm_1=0;y_1=0;
for k=1:1:200
    time(k)=k*ts;
    rin(k)=1;
    xm(k)=-den(2)*xm_1+num(2)*u_1;
    yout(k)=-den(2)*y_1+num(2)*u_5;
    e2(k)=rin(k)-xm(k);
ei=ei+e2(k)*ts;
u(k)=0.50*e2(k)+0.010*ei;
xm_1=xm(k);
u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1; 
if k==100
       u_1=u(k)+1.0;
    else
       u_1=u(k);
end
y_1=yout(k);
end
plot(time,rin,'b',time,yout,'r');
grid;
xlabel('time(s)');
ylabel('rin,yout');



%仿真程序6
%fuzzy tuning pid control