%原結構（單位反饋）的修改程序
%文章階躍信號用于k>=500時過程kp,tol變大.tp變小。
clear all;
close all;
a=readfis('m500nochange');

%被控對象模型：一階純滯后
ts=1;
kp1=1;
Tp1=30;
tol1=30;
sys1=tf([kp1],[Tp1,1],'inputdelay',tol1);
dsys1=c2d(sys1,ts,'zoh');
[num1,den1]=tfdata(dsys1,'v');
rin=1.0;
n1=(tol1/ts);
zz=0;

% signal initialisation
xee(1)=0;
xec(1)=0;
xm(1)=0;
e2(1)=0;
ec(1)=0;
ee(1)=0;
u(1)=0;
r(1)=1;
lanmda0=30;
p=1000;
ym=zeros([1,p]);
y=zeros([1,p]);
for k=2:1:p
%timef(k)=k*ts;
   r(k)=rin;
    zero(k)=zz;
   if k>=500
        kp=kp1*1.4;
       Tp=Tp1*0.6;
        tol=tol1*1.4;
    
    else
        kp=kp1;
       Tp=Tp1;
       tol=tol1;
   end
   n0=(tol/ts);
    sys=tf([kp],[Tp,1],'inputdelay',tol);
    dsys=c2d(sys,ts,'zoh');
    [num,den]=tfdata(dsys,'v');
   xm(k)=-den1(2)*xm(k-1)+num1(2)*u(k-1);
       if k>n0
          y(k)=-den(2)*y(k-1)+num(2)*u(k-n0);
      end
      if k>n1
      ym(k)=-den1(2)*ym(k-1)+num1(2)*u(k-n1);
      end
  %No precise Model:PI+Smith
   xee(k)=r(k)-y(k);
    xec(k)=xee(k)-xee(k-1);
    ee(k)=y(k)-ym(k);
    ec(k)=ee(k)-ee(k-1);
    e1(k)=r(k)-ee(k);
    e2(k)=e1(k)-xm(k);
   alf1=1;
    alf2=1;
    
   %alf1= 1- 0.7*exp(-1*xee(k)^2);
  % alf2= 1- 0.7*exp(-1*xec(k)^2); 
     

   w=[xee(k)/alf1,xec(k)/alf2];
    lanmda_pid=2*evalfis(w,a);
   %lanmda_pid=0;
    lanmda(k)=lanmda0+lanmda_pid; 
    u(k)=u(k-1)+[Tp1/(kp1*lanmda(k))]*e2(k)-[Tp1/(kp1*lanmda(k))]*e2(k-1)+[1/(kp1*lanmda(k))]*e2(k-1);  
end
%displaying the results
figure(1)
plot(r,':')
%axis([3 length(r-10) 0 1.5]);
hold on 
plot(y(2:length(y)));

xlabel('time(s)');
ylabel('r,y');
%figure(2)
%plot(uo,':')
%hold on 
%stairs(u(2:length(u)));
%xlabel('time(s)');
%ylabel('u');
%figure(3)
%plot(e1(2:length(e1)))
%xlabel('time(s)');
%ylabel('e1');
figure(4)
plot(lanmda,'r');
xlabel('time(s)');
ylabel('lanmda');
%figure(5)
%plot(ee(2:length(ee)))
%xlabel('time(s)');
%ylabel('model error');
figure(6)
plot(zero,':')
hold on
plot(xee(2:length(xee)))
xlabel('time(s)');
ylabel('error');