% Paper 2009 Vrancic et al. Improving disturbance rejection of PID 
% controllers by means of the magnitude optimum method;
% ISA Transactions 2009, 10 pages; plants 1 till 4, plant equations (8):
Vs=1;   % same for all 4 plants
num1=1;   num2=[1  1];   num3=[-2  1];  num4=1;  % numerators
den1=conv([9  1],[1  1]);                        % denominators
den2=conv([8  1],conv([1.5  1],[1.5  1]));
den3=conv([4  1],[4  1]);
     den4h=conv([1  1],conv([1  1],conv([1  1],[1  1])));
den4=conv(den4h,den4h);
Td4=2;  % dead time
GPl1=tf(Vs*num1,den1,'Variable','p');  % transfer functions
GPl2=tf(Vs*num2,den2,'Variable','p');
GPl3=tf(Vs*num3,den3,'Variable','p'); 
GPl4=tf(Vs*num4,den4,'Variable','p','IODelay',Td4);

% TF: Filter inside controller <--> Simulink: TF=0.001 typically
TF=0.5;   fipo=[TF   1];                      % order 1, paper (9)

den1f=conv(fipo,den1);
Ru1=bo_u_gen02(13,Vs,num1,den1f), %PID undel [VR,TN,T0,TV,TQ,TT]
% for simulink:
Ru1zae=[ Ru1(4)*Ru1(2)  Ru1(2)  1],
Ru1nen=conv(fipo,[Ru1(3)  0]),    % Ru1(3) = TR0 = integrator
GCu1=tf(Ru1zae,Ru1nen);           % transfer function

den2f=conv(fipo,den2);
Ru2=bo_u_gen02(13,Vs,num2,den2f), %PID undel [VR,TN,T0,TV,TQ,TT]
% for simulink:
Ru2zae=[ Ru2(4)*Ru2(2)  Ru2(2)  1],
Ru2nen=conv(fipo,[Ru2(3)  0]),    % Ru2(3) = TR0 = integrator
GCu2=tf(Ru2zae,Ru2nen);           % transfer function

den3f=conv(fipo,den3);
Ru3=bo_u_gen02(13,Vs,num3,den3f), %PID undel [VR,TN,T0,TV,TQ,TT]
% for simulink:
Ru3zae=[ Ru3(4)*Ru3(2)  Ru3(2)  1],
Ru3nen=conv(fipo,[Ru3(3)  0]),    % Ru3(3) = TR0 = integrator
GCu3=tf(Ru3zae,Ru3nen);           % transfer function

den4f=conv(fipo,den4);
Ru4=bo_u_gen02(13,Vs,num4,den4f,Td4,23), %PID undel [VR,TN,T0,TV,TQ,TT]
% for simulink:
Ru4zae=[ Ru4(4)*Ru4(2)  Ru4(2)  1],
Ru4nen=conv(fipo,[Ru4(3)  0]),    % Ru4(3) = TR0 = integrator
GCu4=tf(Ru4zae,Ru4nen);           % transfer function
% for comparison paper table 2:
Ru_Kx(1,:)=[ Ru1(1)  Ru1(1)/Ru1(2)  Ru1(1)*Ru1(4) ];
Ru_Kx(2,:)=[ Ru2(1)  Ru2(1)/Ru2(2)  Ru2(1)*Ru2(4) ];
Ru_Kx(3,:)=[ Ru3(1)  Ru3(1)/Ru3(2)  Ru3(1)*Ru3(4) ];
Ru_Kx(4,:)=[ Ru4(1)  Ru4(1)/Ru4(2)  Ru4(1)*Ru4(4) ],

Rd1=bo_d_gen(13,Vs,num1,den1f), % PID delayed   [VR,TN,TR0,TV,TQ]
% for simulink:
Fd1=[ Rd1(4)*Rd1(2)  Rd1(2)  1],
Rd1zae=[Rd1(4)*Rd1(2)  Rd1(2)  1],
Rd1nen=conv(fipo,[Rd1(3)  0]),    % Rd1(3) = TR0 = integrator
GCd1=tf(Rd1zae,Rd1nen);           % transfer function
GFd1=tf(1,Fd1);                   % pre-filter transfer function

Rd2=bo_d_gen(13,Vs,num2,den2f), % PID delayed   [VR,TN,TR0,TV,TQ]
% for simulink:
Fd2=[ Rd2(4)*Rd2(2)  Rd2(2)  1],
Rd2zae=[Rd2(4)*Rd2(2)  Rd2(2)  1],
Rd2nen=conv(fipo,[Rd2(3)  0]),    % Rd2(3) = TR0 = integrator
GCd2=tf(Rd2zae,Rd2nen);           % transfer function
GFd2=tf(1,Fd2);                   % pre-filter transfer function

Rd3=bo_d_gen(13,Vs,num3,den3f), % PID delayed   [VR,TN,TR0,TV,TQ]
% for simulink:
Fd3=[ Rd3(4)*Rd3(2)  Rd3(2)  1],
Rd3zae=[Rd3(4)*Rd3(2)  Rd3(2)  1],
Rd3nen=conv(fipo,[Rd3(3)  0]),    % Rd3(3) = TR0 = integrator
GCd3=tf(Rd3zae,Rd3nen);           % transfer function
GFd3=tf(1,Fd3);                   % pre-filter transfer function

Rd4=bo_d_gen(13,Vs,num4,den4f,Td4,39), % PID delayed   [VR,TN,TR0,TV,TQ]
% for simulink:
Fd4=[ Rd4(4)*Rd4(2)  Rd4(2)  1],
Rd4zae=[Rd4(4)*Rd4(2)  Rd4(2)  1],
Rd4nen=conv(fipo,[Rd4(3)  0]),    % Rd4(3) = TR0 = integrator
GCd4=tf(Rd4zae,Rd4nen);           % transfer function
GFd4=tf(1,Fd4);                   % pre-filter transfer function
% for comparison paper table 3:
Rd_Kx(1,:)=[ Rd1(1)  Rd1(1)/Rd1(2)  Rd1(1)*Rd1(4) ];
Rd_Kx(2,:)=[ Rd2(1)  Rd2(1)/Rd2(2)  Rd2(1)*Rd2(4) ];
Rd_Kx(3,:)=[ Rd3(1)  Rd3(1)/Rd3(2)  Rd3(1)*Rd3(4) ];
Rd_Kx(4,:)=[ Rd4(1)  Rd4(1)/Rd4(2)  Rd4(1)*Rd4(4) ],
% pre-filter optimization, order 4
OptPreFi1 = bo_prefi_op(GPl1,GCu1,GCd1,GFd1,4),
OptPreFi2 = bo_prefi_op(GPl2,GCu2,GCd2,GFd2,4),
OptPreFi3 = bo_prefi_op(GPl3,GCu3,GCd3,GFd3,4),
OptPreFi4 = bo_prefi_op(GPl4,GCu4,GCd4,GFd4,4),


