d0=sqrt(Ti*(Ti-4*Td)); Ti0=Ti; Kp=0.5*(Ti+d0)*Kp/Ti; Ti=0.5*(Ti+d0); Td=Ti0-Ti;
      Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp]; 
      dH=Kp*conv([Ti,1],[Td/N,1]); H_Sys=tf(nH,dH);
   end
case 5
   Kp=1.24*(1+0.13*tau/(1-tau))/a; Td=(0.27-0.36*tau)*L/(1-0.87*tau); 
   dd=[Td/N,1]; nn=Kp*[Td*(N+1)/N, 1]; Gc_Sys=tf(nn,dd); H_Sys=[];
end

%-----------------------------------------------------------------------------
%opt_pid is used to design PID controller using the optimal setting algorithm.
%-----------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td]=opt_pid(key,typ,vars)
Ti=0; Td=0; N=vars(4); 
if length(vars)==5, k=vars(1); L=vars(2); T=vars(3); iC=vars(5); tt=0;
else, Kc=vars(1); Tc=vars(2); kappa=(3); tt=1; end
if tt==0
   switch key
   case 2,
      PIDtab=[0.980, 0.712, 0.569, 1.072, 0.786, 0.628;
             -0.892,-0.921,-0.951,-0.560,-0.559,-0.583;
              0.690, 0.968, 1.023, 0.648, 0.883, 1.007;
             -0.155,-0.247,-0.179,-0.114,-0.158,-0.167];
   case 3
      PIDtab=[1.048, 1.042, 0.968, 1.154, 1.142, 1.061; 
             -0.897,-0.897,-0.904,-0.567,-0.579,-0.583;
              1.195, 0.987, 0.977, 1.047, 0.919, 0.892;
             -0.368,-0.238,-0.253,-0.220,-0.172,-0.165;
              0.489, 0.385, 0.316, 0.490, 0.384, 0.315;
              0.888, 0.906, 0.892, 0.708, 0.839, 0.832];
   case 4
      PIDtab=[1.260, 1.053, 0.942, 1.295, 1.120, 1.001;
             -0.887,-0.930,-0.933,-0.619,-0.625,-0.624;
              0.701, 0.736, 0.770, 0.661, 0.720, 0.754;
             -0.147,-0.126,-0.130,-0.110,-0.114,-0.116;
              0.375, 0.349, 0.308, 0.378, 0.350, 0.308;
              0.886, 0.907, 0.897, 0.756, 0.811, 0.813];
   end
   ii=0; if (L/T>1) ii=3; end; tt=L/T; 
   a1=PIDtab(1,ii+iC); b1=PIDtab(2,ii+iC); a2=PIDtab(3,ii+iC); b2=PIDtab(4,ii+iC); 
   Kp=a1/k*tt^b1; Ti=T/(a2+b2*tt); 
   if key==3| key==4, a3=PIDtab(5,ii+iC); b3=PIDtab(6,ii+iC); Td=a3*T*tt^b3; end
else
   switch key
   case 2, Kp=0.361*Kc; Ti=0.083*(1.935*kappa+1)*Tc;   
   case 3, Kp=0.509*Kc; Td=0.125*Tc; Ti=0.051*(3.302*kappa+1)*Tc;
   case 4,
      Kp=(4.437*kappa-1.587)/(8.024*kappa-1.435)*Kc; Ti=0.037*(5.89*kappa+1)*Tc;  Td=0.112*Tc;
   end
end
switch key
case 2, nn=Kp*[Ti,1]; dd=[Ti,0];
case 3
   nn=[Kp*Ti*Td*(N+1)/N, Kp*(Ti+Td/N), Kp];
   dd=Ti*[Td/N,1,0]; Gc_Sys=tf(nn,dd); H_Sys=[];
case 4
   nn=Kp*[Ti,1]; dd=[Ti,0]; dH=Kp*conv([Ti,1],[Td/N,1]);
   nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp]; Gc_Sys=tf(nn,dd); H_Sys=tf(nH,dH);
end

%--------------------------------------------------------------------------------------
%rziegler_nic is used to design PID controller using the refined Ziegler-Nichols algorithm.
%--------------------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td,beta]=rziegler_nic(vars)
Ti=0; Td=0; K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T; 
Kp=1.2/a; Ti=2*L; Td=L/2; Kc=vars(5); Tc=vars(6); kappa=Kc*K; tau=L/T;
if (kappa > 2.25 & kappa<15) | (tau>0.16 & tau<0.57)
   beta=(15-kappa)/(15+kappa);
elseif (kappa<2.25 & kappa>1.5) | (tau<0.96 & tau>0.57)
   mu=4*jappa/9; beta=8*(mu-1)/17; Ti=0.5*mu*Tc;
elseif (kappa>1.2 & kappa<1.5), 
   Kp=5*(12+kappa)/(6*(15+14*kappa));  Ti=0.2*(4*kappa/15+1); beta=1;
end

dH=conv([Ti*beta,1],[Td/N,1]); nH=[Ti*Td*beta*(N+2-beta)/N, Ti+Td/N, 1];
Gc_Sys=tf(Kp*[beta*Ti,1],[Ti,0]); H_Sys=tf(nH,dH);

%---------------------------------------------------------------------------------
%otherpid is used to give an extra dialog box to allow the PID setting using other
%algorithms.
%---------------------------------------------------------------------------------
function otherpid(arg1,arg2)

if nargin==1
   nTask=arg1(1); g_win=findobj('Name','Enter PID Specifications');
   if length(g_win)==0
      g_win=figure('Units','normalized','Position',[0.24875 0.332 0.35 0.167],...
         'NumberTitle','off','Name','Enter PID Specifications','Tag','CtrlLABExtras',...
         'MenuBar','none','Color',0.8*[1,1,1],'Resize','off');
      extra_funs(1);
      switch nTask
      case 9
         [v,d]=version; v1=eval(v(1)); v2=eval(v(3)); v3=eval(v(5));
         if v2==2 & v3==0, strRadio='ToggleButton'; else, strRadio='RadioButton'; end 
         [xL,aText(1)]=display_str(0.05,0.85,'Overshoot Requrements',[0,0,0],'on',9);
         aText(2)=uicontrol('Style',strRadio,'String','No Overshoot (0%)',...
            'Units','normalized','Position',[0.1 0.55 0.45 0.19],...
            'Value',1,'BackgroundColor',0.8*[1,1,1],'Callback','extra_funs(4,1,''Value'',2,3);');
         aText(3)=uicontrol('Style',strRadio,'String','Less than 20%',...
            'Units','normalized','Position',[0.1 0.30 0.45 0.19],...
            'Value',0,'BackgroundColor',0.8*[1,1,1],'Callback','extra_funs(4,1,''Value'',3,2);');
      case 10
         [xL,aText(1)]=display_str(0.05,0.75,'Value of r_b',[0,0,0],'on',9);
         aText(2)=uicontrol('Style','Edit','String','0.45',...
            'Units','normalized','Position',[0.40 0.71 0.25 0.18],...
            'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
         [xL,aText(3)]=display_str(0.05,0.45,'Value of \phi_b',[0,0,0],'on',9);
         aText(4)=uicontrol('Style','Edit','String','45',...
            'Units','normalized','Position',[0.40 0.35 0.25 0.18],...
            'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
      case 11
         [xL,aText(1)]=display_str(0.05,0.75,'Filter Constant T_f',[0,0,0],'on',9);
         aText(2)=uicontrol('Style','Edit','String','10',...
            'Units','normalized','Position',[0.42 0.71 0.2 0.18],...
            'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
      end
      uicontrol('Style','Pushbutton','String','Design',...
         'Units','normalized','Position',[0.76 0.75 0.2 0.2],'Callback','pid_design(6,1);');
      uicontrol('Style','Pushbutton','String','Cancel',...
         'Units','normalized','Position',[0.76 0.45 0.2 0.2],'Callback','delete(gcf);');
      uicontrol('Style','Pushbutton','String','Help',...
         'Units','normalized','Position',[0.76 0.15 0.2 0.2],'Callback','pid_design(6,2);');
      set(gcf,'UserData',{aText,arg1}); 
   else, figure(g_win); end
else
   uu=get(gcf,'UserData'); g_pid=gcf; vars=uu{2}; nTask=vars(1);
   g_main=figure(findobj('Tag','CtrlLABMain'));
   uu0=get(g_main,'UserData'); g1=get(uu0{1}(1),'UserData'); G1=tf(g1{2});
   bFodIdent=extra_funs(5,5,'Checked',29:31); bPIDType=extra_funs(5,5,'Checked',23:26);
   bOvshoot=~get(uu{1}(2),'Value');
   switch arg2
   case 1
      switch nTask
      case 9
         k=vars(2); L=vars(3); T=vars(4); NN=vars(5); 
         [Gc_Sys,H_Sys,Kp,Ti,Td]=chr_pid(bPIDType,1,[k,L,T,NN,bOvshoot]);
      case 10,
         rb=eval(get(uu{1}(2),'String')); phib=eval(get(uu{1}(4),'String'));
         Kc=vars(2); Tc=vars(3); NN=vars(4); 
         [Gc_Sys,H_Sys,Kp,Ti,Td]=ziegler_nic(bPIDType,[Kc,Tc,rb,phib,NN]);
      case 11,
         k=vars(2); L=vars(3); T=vars(4); NN=vars(5); Tf=eval(get(uu{1}(2),'String'));
         [Gc_Sys,H_Sys,Kp,Ti,Td]=imc_pid(bPIDType,[k,L,T,NN,Tf]);
      end;
      close(g_pid);
      proc_pid(Gc_Sys,H_Sys,Kp,Ti,Td,NN);      
   case 2
      switch nTask
      case {9,10,11}, clab_help(14);
      case {2,3,4,5}, clab_help(15);
      end      
   end  
end   

%------------------------------------------------------------------------------
%chr_pid is used to design PID using the Chien-Hoones-Reswick tuning algorithm.
%------------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td]=chr_pid(key,typ,vars)
Ti=0; Td=0; K=vars(1); L=vars(2); T=vars(3); N=vars(4); 
a=K*L/T; ovshoot=vars(5);
if typ==1, TT=T; else TT=L; typ=2;  end
if ovshoot==0, 
   KK=[0.3,0.35,1.2,0.6,1,0.5; 0.3,0.6,4,0.95,2.4,0.42];
else, 
   KK=[0.7,0.6,1,0.95,1.4,0.47; 0.7,0.7,2.3,1.2,2,0.42]; 
end
switch key 
case 1, Kp=KK(typ,1)/a; Gc_Sys=tf(Kp,1); H_Sys=[];
case 2, Kp=KK(typ,2)/a; Ti=KK(typ,3)*TT; Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=[]; 
otherwise
   Kp=KK(typ,4)/a; Ti=KK(typ,5)*TT; Td=KK(typ,6)*L; 
   if key==3
      dd=Ti*[Td/N,1,0]; nn=[Kp*Ti*Td*(N+1)/N, Kp*(Ti+Td/N), Kp]; Gc_Sys=tf(nn,dd); H_Sys=[];
   elseif key==4
      d0=sqrt(Ti*(Ti-4*Td)); Ti0=Ti; Kp=0.5*(Ti+d0)*Kp/Ti; Ti=0.5*(Ti+d0); Td=Ti0-Ti;
      nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp]; dH=Kp*conv([Ti,1],[Td/N,1]);
      Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=tf(nH,dH);
   end
end

%-------------------------------------------------------------------------
%imc_pid is used to design PID using the internal model control structure.
%-------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td]=imc_pid(key,vars)
Ti=0; Td=0; K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T; Tf=vars(5);
switch key
case 2, 
   Kp=T/(K*(L+Tf)); Ti=T; Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=[]; 
otherwise
   Kp=(L/2+T)/(K*(L+Tf)); Ti=T+L/2; Td=L*T/(L+2*T); 
   if key==3
      dd=Ti*[Td/N,1,0]; nn=[Kp*Ti*Td*(N+1)/N, Kp*(Ti+Td/N), Kp]; Gc_Sys=tf(nn,dd); H_Sys=[];
   elseif key==4
      d0=sqrt(Ti*(Ti-4*Td)); Ti0=Ti; Kp=0.5*(Ti+d0)*Kp/Ti; Ti=0.5*(Ti+d0); Td=Ti0-Ti;
      nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp]; dH=Kp*conv([Ti,1],[Td/N,1]);
      Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=tf(nH,dH); 
   end
end

%-----------------------------------------------------------------------------
%display_pid is used to display the designed PID controller in the PID format.
%-----------------------------------------------------------------------------
function display_pid()
g_main=findobj('Tag','CtrlLABMain'); uu=get(g_main,'UserData'); figure(g_main);
bPIDType=extra_funs(5,5,'Checked',23:26); set(uu{4}(26),'Enable','on','Checked','on');
Kp=uu{7}(1); Ti=uu{7}(2); Td=uu{7}(3); N_f=uu{7}(4);
if length(uu{7})==5, beta=uu{7}(5); bPIDType=5; end
%clear information display window to get ready for diaplay
display_str;
%according to PID controller type to display the PID controller
xL=0.1; yL=0.5; display_str(xL,0.75,'PID Controller');
xL=display_str(xL,yL,['G_c(s) = ' display_str(Kp)]);
if bPIDType>1
   %display integral control component
   xL=display_str(xL,yL,' (1+ '); vv=[display_str(Ti) 's']; xx=0.12;
   if length(find(vv=='^'))==0, xx=0.08; end
   xL2=display_str(xL,yL-xx,vv); h=line([xL,xL2],[yL,yL]); set(h,'Color',[0,0,0]);
   display_str(0.5*xL2+0.5*xL,yL+0.065,'1'); xL=xL2; 
   switch bPIDType
   case 2, xL=display_str(xL,yL,' )'); %PI control
   case 3, xL=display_str(xL,yL,[' +' display_str(Td) 's )']); %Normal PID control
   case 4, xL=display_str(xL,yL,['),   F(s)= ', display_str(Td) 's']); %PID control with D in feedback
   case 5, %Refined Ziegler-Nichols structure
      xL=display_str(xL,yL,'),  '); xL=display_str(0.1,yL-0.3,['F(s)=P+I+',display_str(Td) 's']);
      xL=display_str(xL,yL-0.3,[',   \beta=' num2str(beta)]);
   end
end
ctrllab(0,2);

%----------------------------------------------------------------------------------
%get_fod is used to obtain the first-order plus delay approximation of the original 
%model using different algorithms.
%----------------------------------------------------------------------------------
function [K,L,T]=get_fod(G1,method)

K=real(dcgain(G1));
if nargin==1
   [Kc,Pm,wc,wcp]=margin(G1); ikey=0; L=1.6*pi/(3*wc); T=0.5*Kc*K*L; 
   if isfinite(Kc), 
       x0=[L;T]; 
       while ikey==0, 
          ww1=wc*x0(1); ww2=wc*x0(2);
          FF=[K*Kc*(cos(ww1)-ww2*sin(ww1))+1+ww2^2; sin(ww1)+ww2*cos(ww1)];
          J=[-K*Kc*wc*sin(ww1)-K*Kc*wc*ww2*cos(ww1),-K*Kc*wc*sin(ww1)+2*wc*ww2;
              wc*cos(ww1)-wc*ww2*sin(ww1), wc*cos(ww1)]; 
          x1=x0-inv(J)*FF; 
          if norm(x1-x0)<1e-8, ikey=1; else, x0=x1; end
       end
       L=x0(1); T=x0(2);
   end
elseif nargin==2 & method==1
   [nn,dd]=tfdata(G1,'v'); [n1,d1]=tf_derv(nn,dd); [n2,d2]=tf_derv(n1,d1);
   K1=n1(end)/d1(end); K2=n2(end)/d2(end); Tar=-K1/K; T=sqrt(K2/K-Tar^2); L=Tar-T;
end
L=real(L); T=real(T);

%----------------------------------------------------------------------------
%tf_derv is used to obtain the first derivative of a given transfer function.
%----------------------------------------------------------------------------
function [e,f]=tf_derv(b,a)
f=conv(a,a);
e1=conv((length(b)-1:-1:1).*b(1:end-1),a); e2=conv((length(a)-1:-1:1).*a(1:end-1),b);
L0=length(e1)-length(e2); e=[zeros(1,-L0) e1]-[zeros(1,L0) e2];