% Paper 2021 Vrancic/Huba: High-Order Filtered PID Controller Tuning Based
% on Magnitude Optimum; Mathematics 2021, 9, 1340, 24 pages;
% plant equation (18a):

Vs1a=1;
num1a=1;
den1a=conv([1 1],conv([1 1],conv([1 1],[1 1])));
Gplt=tf(num1a,den1a);

% Filter in Regler für Ableitungen (D-Terme)
TFRe=0.1;    % Berücksichtigung deutlicher Reglerdynamik: TFRe > TFNi
TFNi=0.001;  % Nur Berücksichtigung der Realisierung, Simulink typ. N=1000
     % left: paper with TFRe            right: simulink standard TFN=0.001
den_ad1=[TFRe  1];                 den_ad1N=[TFNi  1];
den_ad2=conv(den_ad1,[TFRe  1]);   den_ad2N=conv(den_ad1N,[TFNi  1]);
den_ad3=conv(den_ad2,[TFRe  1]);   den_ad3N=conv(den_ad2N,[TFNi  1]);
den1a1=conv(den_ad1,den1a);        den1a1N=conv(den_ad1N,den1a);
den1a2=conv(den_ad2,den1a);        den1a2N=conv(den_ad2N,den1a);
den1a3=conv(den_ad3,den1a);        den1a3N=conv(den_ad3N,den1a);

% PID  undelayed [VR,TN,T0,TV]
Ru3N=bo_u_gen02(13,Vs1a,num1a,den1a), % Berechnung Nenner ohne TFRe: den1a
Ru3zaeN=[Ru3N(4)*Ru3N(2)  Ru3N(2)  1],%           Ru3(3) = TR0 = integrator
Ru3nenN=conv(den_ad1N,[Ru3N(3)  0]),  % ABER Realisierung Nenner mit: TFNi
Ru3=bo_u_gen02(13,Vs1a,num1a,den1a1), % Berechnung Nenner mit TFRe: den1a1
Ru3zae=[Ru3(4)*Ru3(2)  Ru3(2)  1],              % Ru3(3) = TR0 = integrator
Ru3nen=conv(den_ad1,[Ru3(3)  0]),     % UND Realisierung Nenner mit: TFRe
Cu3=tf(Ru3zae,Ru3nen); Cu3N=tf(Ru3zaeN,Ru3nenN);

% PID2  undelayed [VR,TN,T0,TV,TQ]
Ru4N=bo_u_gen02(14,Vs1a,num1a,den1a), % Berechnung Nenner ohne TFRe: den1a
Ru4zaeN=[Ru4N(5)*Ru4N(2)  Ru4N(4)*Ru4N(2)  Ru4N(2)  1],         %Ru4(3)=TR0
format long
Ru4nenN=conv(den_ad2N,[Ru4N(3)  0]),  % ABER Realisierung Nenner mit: TFNi
format short
Ru4=bo_u_gen02(14,Vs1a,num1a,den1a2), % Berechnung Nenner mit TFRe: den1a2
Ru4zae=[Ru4(5)*Ru4(2)  Ru4(4)*Ru4(2)  Ru4(2)  1],
Ru4nen=conv(den_ad2,[Ru4(3)  0]),     % UND Realisierung Nenner mit: TFRe
Cu4=tf(Ru4zae,Ru4nen); Cu4N=tf(Ru4zaeN,Ru4nenN);

% PID3 undelayed [VR,TN,T0,TV,TQ,TT]
% Ru5N=bo_u_gen02(15,Vs1a,num1a,den1a), % KEINE LÖSUNG ohne TFRe!!!
% Ru5N=bo_u_gen02(15,Vs1a,num1a,den1a3N),
% auch den1a3N liefet keine brauchbare Lösung; Verstärkung zu hoch
Ru5=bo_u_gen02(15,Vs1a,num1a,den1a3), % Berechnung Nenner mit TFRe: den1a3
Ru5zae=[Ru5(6)*Ru5(2)  Ru5(5)*Ru5(2)  Ru5(4)*Ru5(2)  Ru5(2)  1],
Ru5nen=conv(den_ad3,[Ru5(3)  0]),     % UND Realisierung Nenner mit: TFRe
Cu5=tf(Ru5zae,Ru5nen);   %    Cu5N=tf(Ru5zaeN,Ru5nenN); NICHT berechenbar

% PID  delayed [VR,TN,T0,TV]
Rd3N=bo_d_gen02(13,Vs1a,num1a,den1a), % Berechnung Nenner ohne TFRe: den1a
Rd3zaeN=[Rd3N(4)*Rd3N(2)  Rd3N(2)  1],         % Ru3(3) = TR0 = integrator
format long
Rd3nenN=conv(den_ad1N,[Rd3N(3)  0]),  % ABER Realisierung Nenner mit: TFNi    
format short
Fd3N=Rd3zaeN,                         % MO: Filternenner=Reglerzähler
Rd3=bo_d_gen02(13,Vs1a,num1a,den1a1), % Berechnung Nenner mit TFRe: den1a1
Rd3zae=[Rd3(4)*Rd3(2)  Rd3(2)  1],             % Ru3(3) = TR0 = integrator
Rd3nen=conv(den_ad1,[Rd3(3)  0]),     % UND Realisierung Nenner mit: TFRe
Fd3=Rd3zae,
Cd3=tf(Rd3zae,Rd3nen); Cd3N=tf(Rd3zaeN,Rd3nenN);
Fd3=tf(1,Rd3zae);      Fd3N=tf(1,Rd3zaeN);
Fd3Opt=bo_prefi_op(Gplt,Cu3,Cd3,Fd3,2),     % optimized pre-filters
Fd3O4dig=tf([ 0.9468 0.5649 0.4179 ],[ 1 0.8616 0.4179 ]),  % 5 digits only
Fd3NOpt=bo_prefi_op(Gplt,Cu3N,Cd3N,Fd3N,2),
Fd3NO4dig=tf([ 0.9551 0.5806 0.4487 ],[ 1 0.889 0.4487 ]),  % 5 digits only

% PID2 delayed  [VR,TN,TR0,TV,TQ]
Rd4N=bo_d_gen02(14,Vs1a,num1a,den1a), % Berechnung Nenner ohne TFRe: den1a
Rd4zaeN=[Rd4N(5)*Rd4N(2)   Rd4N(4)*Rd4N(2)  Rd4N(2)  1],
format long
Rd4nenN=conv(den_ad2N,[Rd4N(3)  0]),  % ABER Realisierung Nenner mit: TFNi
format short
Fd4N=Rd4zaeN,
Rd4=bo_d_gen02(14,Vs1a,num1a,den1a2), % Berechnung Nenner mit TFRe: den1a2
Rd4zae=[Rd4(5)*Rd4(2)   Rd4(4)*Rd4(2)  Rd4(2)  1],
Rd4nen=conv(den_ad2,[Rd4(3)  0]),     % UND Realisierung Nenner mit: TFRe
Fd4=Rd4zae,
Cd4=tf(Rd4zae,Rd4nen); Cd4N=tf(Rd4zaeN,Rd4nenN);
Fd4=tf(1,Rd4zae);      Fd4N=tf(1,Rd4zaeN);
Fd4Opt=bo_prefi_op(Gplt,Cu4,Cd4,Fd4,3),     % optimized pre-filters
Fd4O4dig=tf([ 0.924 1.962 2.032 1.626 ],[ 1 2.542 3.169 1.626 ]),% 5 digits
Fd4NOpt=bo_prefi_op(Gplt,Cu4N,Cd4N,Fd4N,3),
Fd4NO4dig=tf([0.9976 2.576 2.853 2.745],[1 3.138 4.416 2.745]),  % 5 digits

% PID3 delayed  [VR,TN,TR0,TV,TQ,TT]
%optim=optimset('MaxFunEvals',80000,'MaxIter',80000);
%Rd5=bo_d_gen02(15,Vs1a,num1a,den1a,[],[],[36 36 0.05 0.38 0.2 0.2],optim),
                                      % KEINE LÖSUNG ohne Reglerdynamik !!!
% !!! Um fsolve zu Lösung zu bringen ist TFR nötig !!!
% OHNE TFR findet fsolve keine Lösung
%optim=optimset('MaxFunEvals',80000,'MaxIter',80000)
%Rd5N=bo_d_gen02(15,Vs1a,num1a,den1a3N,[],[],[],optim), % TFNi ist zu klein

TFS=0.06  % oder doch 0.05 ??? Größerer Wert als TFNi nötig
den_ad1S=[TFS  1];                 den_ad2S=conv(den_ad1S,[TFS  1]);
den_ad3S=conv(den_ad2S,[TFS  1]);  den1a3S=conv(den_ad3S,den1a);
optim=optimset('MaxFunEvals',80000,'MaxIter',80000);
Rd5N=bo_d_gen02(15,Vs1a,num1a,den1a3S,[],[],[],optim) 
format long % KEINE Berechnung Nenner ohne TFRe (den1a)
            % ABER Berechnung Nenner mit TFS (den1a3S): TFNi < TFS < TFRe
Rd5zaeN=[ Rd5N(6)*Rd5N(2)   Rd5N(5)*Rd5N(2)   Rd5N(4)*Rd5N(2)  Rd5N(2)  1],

%Rd5nenN=conv(den_ad3N,[Rd5N(3)  0]), % ABER Realisierung Nenner mit: TFNi
Rd5nenS=conv(den_ad3S,[Rd5N(3)  0]), % UND Realisierung Nenner mit: TFS

format short
FNed5=Rd5zaeN,
Rd5=bo_d_gen02(15,Vs1a,num1a,den1a3), % Berechnung Nenner mit TFRe: den1a3
Rd5zae=[ Rd5(6)*Rd5(2)   Rd5(5)*Rd5(2)   Rd5(4)*Rd5(2)  Rd5(2)  1],
format long
Rd5nen=conv(den_ad3,[Rd5(3)  0]),     % UND Realisierung Nenner mit: TFRe
format short
Fd5=Rd5zae,
Cd5=tf(Rd5zae,Rd5nen); Cd5N=tf(Rd5zaeN,Rd5nenS);%Achtung hier "S" statt "N"
Fd5=tf(1,Rd5zae);      Fd5N=tf(1,Rd5zaeN);
Fd5Opt=bo_prefi_op(Gplt,Cu5,Cd5,Fd5,4),     % optimized pre-filters
     % Fd5NOpt=bo_prefi_op(Gplt,Cu5N,Cd5N,Fd5N,4), % Cu5N nicht berechenbar
Fd5O4dig=tf([0.8791 4.087 14.07 20.36 29.61],[1 5.913 22.48 34.56 29.61]),
                                                            % 5 digits only
disp('"Files ex1_sim_Vgl_XXX.slx" will be called - compare parameters')
ex1_sim_Vgl_und; sim('ex1_sim_Vgl_und',40),
ex1_sim_Vgl_del; sim('ex1_sim_Vgl_del',40),
ex1_sim_Vgl_del_FO; sim('ex1_sim_Vgl_del_FO',40),
ex1_sim_Vgl_del_FO_4d; sim('ex1_sim_Vgl_del_FO_4d',40),
disp('Repeat bo_d_gen call at MATLAB promt, if search algorithm fails')
