%PID
close all;
clear all;
%u-T1
num1=-1.4261;
den1=conv([10 1],[30 1]);
Midsys1=tf(num1,den1);
% Kp2=0.7;
%  Kp2=-0.252;
 Kp2=-1.03;
sys1=Kp2*Midsys1/(1+Kp2*Midsys1);
%T1-T2
num2=0.3181;
den2=conv([30 1],[90 1]);
sys2=tf(num2,den2);
%T-T2
num3=0.8022;
den3=conv([90 1],[90 1]);
sys3=tf(num3,den3);
%T-u
sys=sys1*sys2*sys3;
%sample time
Ts=40;
sysd=c2d(sys,Ts,'zoh');
set(sysd,'variable','z^-1');
A=sysd.den{1};   
B=sysd.num{1};
n_a=length(A)-1;
n_b=length(B)-1;
kp=3.0;
ki=kp*Ts/219;
kd=kp*20/Ts;
r=1;
N=100;
%simulation step 
S_Time=300;
y_actual1=0;
u=0;
e=0;
deltau=0;
for i=1:S_Time
    time(i)=i*Ts;
    r(i)=1;
    e(i)=r(i)-y_actual1(max(1,i-1));
    u(i)=u(max(1,i-1))+kp*(e(i)-e(max(1,i-1)))+ki*e(i)+kd*(e(i)-2*e(max(1,i-1))+e(max(1,i-1)));
     y_actual1(i)=B*u(max(1,[i-1:-1:i-1-n_b]))'-A(2:end)*y_actual1(max(1,[i-1:-1:i-n_a]))';
end
hold on 
grid on
plot(time,r,'r--',time,y_actual1,'b');
title('PID control');
xlabel('Time (sec)');