
%Written by Iosif Pappas in Spring 2016.
%The script approximates a high-order process model to a second-order process with time delay.
%In addition, it calculates the ideal PID SIMC parameters which satisfy the desired Ms values. Finally the resulting PID settings can be sent to a Simulink file to calculate the IAE.

tic
%Set the Laplace variable
s=tf('s');


%Original Model G(s) = ((20*s+1)*(14*s+1)*exp(-theta*s))/((42*s+1)*(5*s+1)*(2*s+1)*(s+1))

%process time constants
To_1  = 20;
To_2  = 14;
tau_1 =42;
tau_2 =5;
tau_3 =2;
tau_4 =1;

%counters
p=0
j=0

%specify the desired Ms value
for ms = 1.29:0.01:2
p=0
j=j+1

for tc = 0:0.01:5
p=p+1;

%Half-rule for (14*s+1)/(5*s+1)
if (To_2>=tau_2) && (tau_2>=tc)
g1(p,1) = To_2/tau_2;
elseif (To_2>=tc) && (tc>=tau_2)
g1(p,1) = To_2/tc;
elseif (tc>=To_2) && (To_2>=tau_2)
g1(p,1)= 1;
else
g1(p,1)=1
end

%Half-rule for (20*s+1)/((42*s+1)
if (tau_1>=To_1) && (To_1>=5*tc);
h = To_1/tau_1;
H(p,1) = (g1(p,1)*h*exp(-1*s))/((2*s+1)*(s+1))
q=0;
k=1;
tau1=2;
tau2=1;
kp=g1(p,1)*h;
else
tau = min(tau_1,5*tc);
q = tau-To_1;
v = tau/tau_1;

if (q>2)
w = q;
H(p,1) = (g1(p,1)*v*exp(-1.5*s))/((w*s+1)*(2.5*s+1));
k=1.5;
tau1=w;
tau2=2.5;
kp=g1(p,1)*v;
elseif (q<2)&&(q>1)
H(p,1) = (g1(p,1)*v*exp(-1.5*s))/((2*s+1)*((q+0.5)*s+1));
k=1.5;
tau1=2;
tau2=q+0.5;
kp=g1(p,1)*v;
elseif (q<1)
H(p,1) = (g1(p,1)*v*exp(-(1+q/2)*s))/((2*s+1)*((1+q/2)*s+1));
k=1+q/2;
tau1=2;
tau2=1+q/2;
kp=g1(p,1)*v;
end

end

% end
G(p,1) = H(p,1);

kc = ((tau1)/(kp*(tc+k)));
taud = tau2;
%set integral time
if tau1 < (4*(tc+k));
ti=tau1;
else
ti = 4*(tc+k);
end
f = 1 + (taud/ti);
kc = kc*f;
ti = ti*f;
taud = taud/f;

%Calculate the Ms
[G_n G_d] = tfdata(G(p,1),'v');
delay = totaldelay(G(p,1));
w = logspace(-10,10,10^4);
GG_no = conv([kc*ti*taud kc*ti kc],G_n);
GG_de = conv([ti 0],G_d);
GG_w1 = polyval(GG_no,w*1i);
GG_w2 = polyval(GG_de,w*1i);
delay_w = exp(-delay*w*1i);
G_G = GG_w1./GG_w2.*delay_w;
S = 1./abs(1+G_G);
MS_S = max(S);

%save the values which satisfy the requirements
%if Ms = 1.59 save the controller parameters
if (abs(ms-MS_S)<1e-4)

break;

%print the results
end
b(j,1) = kc;
b(j,2) = ti;
b(j,3) = taud;
b(j,4) = tc;
b(j,5) = k;
b(j,6) = MS_S;
b(j,7) = ms;
b(j,8) = p;

end
end

%If needed, calculate the IAE in Simulink
for p=1:1:71
SimOut = sim('Test');
h(p,1) = output1;
h(p,2) = output2;
end

toc

