help me smith predictor and Ziegler Nichols

My question is to help me to write program to find optimal PID parameter by using x = fminimax(@myfun,x0,Ainq,Binq,Aeq,Beq,LB,UB); Kp=x(1); Ki=x(2); Kd=x(3); N=x(4); Smith predictor method And Ziegler Nichols method I want to synthesis study behaviors of system and do comparison between smith predictor method and Ziegler Nichols method Study margin of stability Study robustness of the method

clear all;clc;close all; s=tf('s'); % l'opèrateur de Laplace

%---------> Transfer function of model system <----------------% %****** km *****% %****** Gm(s)= ------------ exp(-rhom*s) ******% %****** (1+taum*s)^2 *****% %--------------------------------------------------------------%

    km=1;rhom=5.57;taum=1.64;
    Gm0=km/((0.999+taum*s)^2);
    delay=exp(-rhom*s);
    Gm=Gm0*delay;
    Ym=step(Gm); plot(Ym);

%---------> Transfer function of real system <----------------% %****** kr *****% %****** Gr(s)= ----------------- exp(-rhor*s) ******% %****** (0.999+taur*s)^5 *****% %--------------------------------------------------------------%

    kr=1;rhor=4;taur=1;
    Gr=(kr*exp(-rhor*s))/((0.999+taur*s)^5);
    Yr=step(Gr); plot(Yr);

%---------> optimisation<--------------------------------------% %****** 1+alpha*s *****% %****** | W1*Sd | W1=L1 ---------- *****% %****** Min | | 1 alpha*s *****% %****** | W3*Sc | W3=L2 ---------- *****% %****** 1+beta*s *****% %****** fonctions de sensibilité *****% %****** 1 GK *****% %****** Sd= ------, Sc= ------, Sd+Sc=1; *****% %****** 1+GK 1+GK *****% %****** W1 et W3 donner par l'utilisateur *****% %****** x=[Kp; Ki; Kd; N] *****% %--------------------------------------------------------------%

    L1=0.9;    L2=0.9;    alpha=12;   beta=12;
    Ainq=[]; Binq=[];
    Aeq=[];  Beq=[];
    LB=[0.01;  0.01;  0.0001;     99];
    UB=[   100;     100;      100;   100]; 
    x0=[0.01;  0.01;  0.0001;      1];
    x = fminimax(@myfun,x0,Ainq,Binq,Aeq,Beq,LB,UB);
    Kp=x(1); Ki=x(2); Kd=x(3); N=x(4);

w=logspace(-3,3,300);% la plage fréquentielle allant de 10^(-3) jusqu'à 10^(3), on a 50 fréquences for i=1:length(w) j=sqrt(-1); p=j*w(i); G=(km*exp(-rhom*p))/(1+taum*p)^2; C0=Kp+(Ki/p)+(Kd*p)/(1+N*p); GK=G*C0; Sd=1/(1+GK); Sc=GK/(1+GK); W1=L1*(1+alpha*p)/(alpha*p); W3=L2*(1+beta*p); f(i)=max(svd([W1*Sd;W3*Sc]));

% valeurs singulières maximales

G_svd(i)=max(svd(G)); C0_svd(i)=max(svd(C0));

Sd_svd(i)=max(svd(Sd)); Sc_svd(i)=max(svd(Sc)); W1_svd(i)=max(svd(W1));W1_inv(i)=1/W1_svd(i); W3_svd(i)=max(svd(W3));W3_inv(i)=1/W3_svd(i);

end

%---------> PID Controller;(Kp,Ki,Kd,N)<-----------------------% %****** Ki Kd *****% %****** K(s)= Kp + ----- + ------- *****% %****** s (1+Ns) *****% %****** 1+N*s: filtrage de l'action dérivèe *****% %****** x=[ Kp Ki Kd N] *****% %--------------------------------------------------------------% % Kp=0.592060348482518; Ki=0.209359822894674; Kd=0; N=100;

    K=Kp+(Ki/s)+(Kd*s)/(1+N*s);

%---------> System step responce <-----------------------------% %****** *****% %****** *****% %****** *****% %--------------------------------------------------------------%

    [ak,bk,ck,dk]=ssdata(K);
    [am0,bm0,cm0,dm0]=ssdata(Gm0);
    [am,bm,cm,dm]=ssdata(Gm);
    [ar,br,cr,dr]=ssdata(Gr);

    [abo,bbo,cbo,dbo]=series(ak,bk,ck,dk,ar,br,cr,dr);
    [abf,bbf,cbf,dbf]=cloop(abo,bbo,cbo,dbo);

    ym=step(abf,bbf,cbf,dbf,1,0:100);

%---------> figure <-------------------------------------------% %****** *****% %****** *****% %****** *****% %--------------------------------------------------------------%

plot(ym); grid on figure,loglog(w,Sd_svd,'r');hold on;loglog(w,W1_inv,'b');grid;title('bleu: inverse de W1 rouge: sensibilité Sd'); figure,loglog(w,Sc_svd,'r');hold on;loglog(w,W3_inv,'b');grid;title('bleu: inverse de W3 rouge: sensibilité Sc'); figure,plot(ym); grid