% My Heater Transfer Function
close all;clear all;clc
% Water Heater Volume
r = 0.07; % radius
h = 0.1;  % height  
V = pi*(r^2)*h; % volume

% Thermal Capacitance of The Liquid
C = V*(4.2*10^6); % for Water: row*ch=4.2*10^6

% Thermal Resistance of The Wall
l = 0.0012; % thickness
k = 0.1; % conductivity of the material, PP (Polypropylene): 0.1-0.22
R = l/k;

% Dead time process
L = 30;

% Constant Provide by the Manufacture
K = 300;
T = C*R;
% Transfer function
num = K*R;
den = [T 1];
H = tf(num,den,'InputDelay',L);
H = pade(H,10);
% H = exp(-L*s)*K*R/(T*s+1);

% Maclaurin Series pole
p3 = (((C*R*(L^2))/2)+((L^3)/6))/(K*R);
p2 = (C*R*L+((L^2)/2))/(K*R);
p1 = (C*R+L)/(K*R);
p0 = 1/(K*R);

% PD feedback
chi = 1.0412;
Beta2 = (chi+((chi^2)/2))/(1+chi)^2;
Beta3 = 0.08;
sigma = (Beta2*p3)/(Beta3*p2);
f0 = (p2/(Beta2*(sigma^2)))-p0;
f1 = ((p0+f0)*sigma)-p1;
Kf = f0; 
Tf = f1/f0; 
kappa = 0.1;

% Controller Parameter Setting
Kc = f0+p0;
Tc = sigma/(1+chi);
Lc = chi*Tc;

% Tuning parameters
lamda = 1; % respond speed parameter
alfa = 1-((1-lamda)^2)*exp(-Lc/Tc); % distubance regulation parameter

s = tf('s');
num1 = 1+lamda*Tc*s;
den1 = 1+alfa*Tc*s;
SP = num1/den1;

num2 = (1+Tc*s)*(1+alfa*Tc*s);
den2 = (1+lamda*Tc*s)^2;
Q = num2/den2;

M = tf(1,[Tc 1],'InputDelay',Lc);
M = pade(M,10);

num3 = Kf*(1+Tf*s);
den3 = 1+kappa*Tf*s;
PD = num3/den3;

G = feedback(H,num3);
QM = feedback(Q,-M);
Co = QM*Kc;
A = Co*G;
B = feedback(A,1);
MD = SP*B;

% PI Controller Parameters
% Kp = 0.2178;
% Ki = 0.003866;
% PIController = pid(Kp,Ki);
% Gpid = feedback(H*PIController,1);
% Show Response Analysis
t = 0:10:300;
figure
step(t,MD);
title('MD PID vs PI Controller')
grid on
% hold on
% step(t,Gpid);
% legend('MD PID','PI');
% hold off