%PID_DESIGN is the function used in CtrlLAB for PID controller design.  In this 
%function, different tuning algorithms are implemented, and details can be referred 
%from the lecture notes.
%
%Available Algorithms:
%
%   1 for Ziegler-Nichols tuning formula
%   2-4 for optimum setting, with ISE, ISTE and IST2E
%   5 for optimum gain-phase margin setting
%   7 for refined Ziegler-Nichols tuning formula
%   8 for Cohen-Coon tuning formula
%   9 for Chien-Hoones-Reswick tuning algorithm
%  10 for modified Ziegler-Nichols tuning
%  11 for PID with internal model control structure
%
%bPIDType variable for different PID structures
%
%   1 for P, 2 for PI, 3 for normal PID
%   4 for PID with D in feedback
%   5 for refined Ziegler-Nichols tuning structure
%
%bFodIdent variable for first order delay approximation
%
%   1 for optimal reduction based method
%   2 for frequency response based method
%   3 for transfer function based method
%
%Available lists of functions under pid_design()
%
%   design_pid --- to administrate PID tuning algorithms
%   pid_dat_win -- allows the user to set PID parameters manually
%   pid_set --- obtain the PID controller from the dialog box
%   proc_pid --- process PID parameters and set them to windows
%   ziegler_nic --- implements Ziegler-Nichols tuning formula
%   cohen_pid --- implements Cohen-Coon tuning formula
%   opt_pid --- implements optimum PID controller setting
%   rziegler_nic --- implements refined Ziegler-Nichols algorithm
%   otherpid -- administrates PID with other algorithms
%   chr_pid --- implements Chien-Hoones-Reswick tuning algorithm
%   imc_pid --- implements internal model control PID structure
%   display_pid -- displays the PID controller obtained
%   get_fod --- obtain the first-order plus delay approximation 
%   tf_derv --- obtain the first derivative of transfer function
%

%Designed by Professor Dingyu Xue
%School of Information Science and Engineering, Northeastern University
%Shenyang 110006, P R China
%Email: xue_dy@hotmail.com
%
%This module is for CtrllAB 3.0, (c) 1996-1999
%Last modified 5 October, 1999
%-------------------------------------------------------------------------

function pid_design(nTask,arg1)
switch nargin
case 0,    pid_dat_win;
   %enter the PID controller in a manual way, i.e., enter 
   %the PID parameters in the fialog box provided.
case 1,    design_pid(nTask);
   %for certain algorithms, if there is no extra parameters expected, the PID 
   %controller can be design directly.
case 2,    otherpid(nTask,arg1);   
   %if there are extra parameters expected in the design, then the parameters 
   %provided in the dialog box can be used in design.
end   

%---------------------------------------------------------------------
%design_pid is used to design PID using the selected tuning algorithm.
%---------------------------------------------------------------------
function design_pid(nTask)
g_main=figure(findobj('Tag','CtrlLABMain')); uu=get(g_main,'UserData');
g1=get(uu{1}(1),'UserData'); G1=tf(g1{2}); z=roots(G1.den{1}); G0=G1;
bFodIdent=extra_funs(5,5,'Checked',29:31); bPIDType=extra_funs(5,5,'Checked',23:26);
if length(find(real(z)>=0))>0
   %display error message for unstable plants 
   warndlg('Sorry, PID not available now for unstable plant','Warning: Design failed!'), 
   return
else
   %get gain/phase margins of the original plant
   g1=get(uu{1}(4),'UserData'); Td=g1{1}; nPade=g1{2};
   if Td>0, [nn,dd]=pade(Td,nPade); G1=G1*tf(nn,dd); end
   [Kc,pp,wg,wp]=margin(G1); den=G1.den{1};
   if any(2:5==nTask) & abs(den(end))<eps
      warndlg('There exists pole at s=0, PID controller not found!','Warning: Design failed!'), 
      return;   
   end
end

Tc=2*pi/wg; 
if any([2:4,7:11]==nTask)|(any([1,5]==nTask)&isfinite(Kc))
   %get the first order plus delay approximation to the plant
   switch bFodIdent
   case 1,
      extra_funs(12,'Optimal first-order delay model fitting ...');
      pause(0.0002);
      G_Sys=G1; G_Sys.Td=Td;
      G_red=mod_reduction(G_Sys,[0,1],[1,0,0,0,0],[2,1],[nPade,1]);
      [nc,dc,L]=tfdata(G_red,'v'); k=nc(2)/dc(2); T=dc(1)/dc(2); 
      delete(findobj('Name','Please Wait'));
   case 2, [k,L,T]=get_fod(G1); %frequency response based fitting
   case 3, [k,L,T]=get_fod(G1,1); %transfer function based fitting
   end
   NN=uu{7}(4);
   
   %design PID controller using different algorithms
   switch nTask
   case 1
      if isfinite(Kc), [Gc_Sys,H_Sys,Kp,Ti,Td]=ziegler_nic(bPIDType,[Kc,Tc,NN]);
      else, [Gc_Sys,H_Sys,Kp,Ti,Td]=ziegler_nic(bPIDType,[k,L,T,NN]); end
   case {8,9,10,11}
      if nTask==10, vars=[nTask,Kc,Tc,NN];
      else, vars=[nTask,k,L,T,NN,Kc,Tc]; end
      if nTask==8, [Gc_Sys,H_Sys,Kp,Ti,Td]=cohen_pid(bPIDType,[k,L,T,NN]);
      else, otherpid(vars); %allow an extra dialog box to rea specifications of design
      end
   case {2,3,4,7}
      if bPIDType==1 & bFodIdent==1, bPIDType=2; end
      if nTask==7, 
         [Gc_Sys,H_Sys,Kp,Ti,Td,beta]=rziegler_nic([k,L,T,NN,Kc,Tc]);
      else
         [Gc_Sys,H_Sys,Kp,Ti,Td]=opt_pid(bPIDType,1,[k,L,T,NN,nTask-1]);
      end
   case 5
      kappa=dcgain(G1)*Kc;
      if bPIDType==1, bPIDType=2; end
      [Gc_Sys,H_Sys,Kp,Ti,Td]=opt_pid(bPIDType,1,[Kc,Tc,kappa,NN]);
   end
else,
   warndlg('Sorry, PID not found for the plant','Warning: Design failed!'), 
   return;   
end

%if PID controller is design, call PID display & process function
if nTask==7, proc_pid(Gc_Sys,H_Sys,Kp,Ti,Td,NN,beta);
elseif any([1:5,8]==nTask), proc_pid(Gc_Sys,H_Sys,Kp,Ti,Td,NN); end

%-------------------------------------------------------------------------------
%pid_proc is used to save PID controller information and display the controller.
%-------------------------------------------------------------------------------
function proc_pid(Gc_Sys,H_Sys,Kp,Ti,Td,Nf,beta)

%write PID controller to the main window objects
h_main=findobj('Tag','CtrlLABMain'); uu=get(h_main,'UserData');
[nc,dc]=tfdata(Gc_Sys,'v'); set(uu{1}(2),'UserData',{1,Gc_Sys,mat2str(nc),mat2str(dc)});
if length(H_Sys)>0, 
   [nh,dh]=tfdata(H_Sys,'v'); set(uu{1}(3),'UserData',{1,H_Sys,mat2str(nh),mat2str(dh)});
else, set(uu{1}(3),'UserData',[]); end
uu{7}(1:3)=[Kp,Ti,Td]; if nargin==7, uu{7}(5)=beta; else, uu{7}=uu{7}(1:4); end
%display PID controller
set(h_main,'UserData',uu); display_pid;

%-------------------------------------------------------------------------------
%pid_dat_win is used to specify the parameters of PID controller manually from 
%the user.  The data will be passed in the interface in the 7-th UserData of the 
%main window.
%-------------------------------------------------------------------------------
function pid_dat_win

%read in PID parameters manually using the dialog box
h_win=findobj('Name','Enter PID parameters');
if length(h_win)==0
   uu=get(findobj('Tag','CtrlLABMain'),'UserData'); keyD=0;
   if strcmp(get(uu{5}(25),'Checked'),'on'), keyD=1; end
   if uu{7}(1)==0, K_p=1; T_i=Inf; T_d=0; N_f=10;
   else, K_p=uu{7}(1); T_i=uu{7}(2); T_d=uu{7}(3); N_f=uu{7}(4); end   
   h_win=figure('Units','normalized','Position',[0.24875 0.332 0.35 0.233],...
      'NumberTitle','off','Name','Enter PID Parameters','Tag','CtrlLABExtras',...
      'MenuBar','none','Color',0.8*[1,1,1],'Resize','off');
   extra_funs(1); display_str(0.05,0.83,'Proportional K_p',[0,0,0],'on',9);
   aText(1)=uicontrol('Style','Edit','String',num2str(K_p),...
      'Units','normalized','Position',[.39 .78 .25 .13],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   display_str(0.05,0.65,'Integral T_i',[0,0,0],'on',9);
   aText(2)=uicontrol('Style','Edit','String',num2str(T_i),...
      'Units','normalized','Position',[0.39 .60 .25 .13],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   display_str(0.05,0.47,'Derivative T_d',[0,0,0],'on',9);
   aText(3)=uicontrol('Style','Edit','String',num2str(T_d),...
      'Units','normalized','Position',[.39 .42 .25 .13],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   display_str(0.05,0.29,'Filter Constant N_f',[0,0,0],'on',9);
   aText(4)=uicontrol('Style','Edit','String',num2str(N_f),...
      'Units','normalized','Position',[.39 .24 .25 .13],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   [v,d]=version; v1=eval(v(1)); v2=eval(v(3)); v3=eval(v(5));
   if v2==2 & v3==0, strCheck='ToggleButton'; else, strCheck='CheckBox'; end 
   aText(5)=uicontrol('Style',strCheck,'String','D in Feedback','Value',keyD,...
      'Units','normalized','Position',[0.15 0.07 0.45 0.14],'BackgroundColor',0.8*[1,1,1]);
   hOK=uicontrol('Style','Pushbutton','String','Apply',...
      'Units','normalized','Position',[0.76 0.68 0.2 0.17],'Callback','pid_design(2);');
   hCancel=uicontrol('Style','Pushbutton','String','Cancel',...
      'Units','normalized','Position',[0.76 0.48 0.2 0.17],'Callback','close(gcf);');
   set(gcf,'UserData',aText);
else, figure(h_win); end

%----------------------------------------------------------------------------
%pid_set is used to obtain the PID controller parameters from the dialog box.
%----------------------------------------------------------------------------
function pid_set()
%get PID parameters from the dialog box
uu0=get(gcf,'UserData');
K_p=eval(get(uu0(1),'String')); T_i=eval(get(uu0(2),'String'));
T_d=eval(get(uu0(3),'String')); N_f=eval(get(uu0(4),'String'));
keyD=get(uu0(5),'Value'); close(gcf);

figure(findobj('Tag','CtrlLABMain'));
if keyD==1, extra_funs(4,5,'Checked',26,23:25);
else, extra_funs(4,5,'Checked',25,[23:24,26]); end   
proc_pid(Gc_Sys,H_Sys,K_p,T_i,T_d,N_f);

%------------------------------------------------------------------------------
%ziegler_nic is used to design PID controllers using the Ziegler-Nichols algorithm.
%------------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td]=ziegler_nic(key,vars)
Ti=0; Td=0;
switch length(vars)
case 3, 
   K=vars(1); Tc=vars(2); N=vars(3);
   switch key
   case 1, Kp=0.5*K; 
   case 2, Kp=0.4*K; Ti=0.8*Tc; 
   case {3,4}, Kp=0.6*K; Ti=0.5*Tc; Td=0.12*Tc; 
   end
case 4, 
   K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T; 
   switch key
   case 1, key==1,  Kp=1/a; 
   case 2, Kp=0.9/a; Ti=3.33*L; 
   case {3,4}, Kp=1.2/a; Ti=2*L; Td=L/2; 
   end
case 5, 
   K=vars(1); Tc=vars(2); rb=vars(3); pb=pi*vars(4)/180; 
   N=vars(5); Kp=K*rb*cos(pb); 
   if key==2, Ti=-Tc/(2*pi*tan(pb)); 
   elseif key==3 | key==4, Ti=Tc*(1+sin(pb))/(pi*cos(pb)); Td=Ti/4; end
end
switch key
case 1, Gc_Sys=tf(Kp,1); H_Sys=[];
case 2, Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=[];
case 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=[];
case 4
   d0=sqrt(Ti*(Ti-4*Td)); Ti0=Ti; 
   Kp=0.5*(Ti+d0)*Kp/Ti; Ti=0.5*(Ti+d0); Td=Ti0-Ti; nn=Kp*[Ti,1]; dd=[Ti,0]; 
   nH=[(1+Kp/N)*Ti*Td, Kp*(Ti+Td/N), Kp]; dH=Kp*conv([Ti,1],[Td/N,1]);
   Gc_Sys=tf(nn,dd); H_Sys=tf(nH,dH);
end

%------------------------------------------------------------------------
%cohen_pid is use to design PID controller using the Cohen-Coon algorithm.
%------------------------------------------------------------------------
function [Gc_Sys,H_Sys,Kp,Ti,Td]=cohen_pid(key,vars)
Ti=0; Td=0; K=vars(1); L=vars(2); T=vars(3); N=vars(4); a=K*L/T; tau=L/(L+T); 
switch key
case 1, Gc_Sys=tf((1+0.35*tau/(1-tau))/a,1); H_Sys=[];
case 2, 
   Kp=0.9*(1+0.92*tau/(1-tau))/a; Ti=(3.3-3*tau)*L/(1+1.2*tau);
   Gc_Sys=tf(Kp*[Ti,1],[Ti,0]); H_Sys=[];
case {3,4}
   Kp=1.35*(1+0.18*tau/(1-tau))/a; Ti=(2.5-2*tau)*L/(1-0.39*tau); Td=0.37*(1-tau)*L/(1-0.81*tau); 
   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