% programa 17.1, Luyben, pag. 600
% "BLT tuning for Wood and Berry column"
% modif. E.Quintero.M 17nov'98
%-------------------------------------------------------
clear
format short e
% process transfer function data
kp=[12.8 -18.9; 6.6 -19.4]; xx=[];
n=2; D=[1 3; 7 3];
tau=[16.7 21; 10.9 14.4];
% BLT
flagp=1;
flagm=1;
loop=0;
fblt=10;
df=2;
% find ziegler-nichols settings
dw=10^(1/80);
[ku,wu]=znt(n,kp,tau,D);
while(1)
   for k=1:n
      kc(k)=ku(k)/(2.2*fblt);
      reset(k)=fblt*2*pi/(wu(k)*1.2);
   end
   dbmax=-200;
   w=0.01;
   loop=loop+1;
   if(loop>30) error('BLT loop'); end
   while(w<10)
%     evaluate process transfer function matrix
      s=i*w;
      for j=1:n
         for k=1:n
            gm(j,k)=kp(j,k)*exp(-D(j,k)*s)/(tau(j,k)*s+1);
         end
      end
%     evaluate diagonal PI controller matrix
      b=zeros(n);
      for j=1:n
         b(j,j)=kc(j)*(reset(j)*s+1)/(reset(j)*s); 
      end
      q=gm*b;
      char=diag(ones(n,1),0)+q;
      xdeter=det(char);
%     form w function=det(1+gb)-1
      xdeter=xdeter-1;
      cltf=xdeter/(1+xdeter);
      dbcl=20*log10(abs(cltf));
      if(dbcl>dbmax)
         dbmax=dbcl;
         wmax=w;
      end
      w=w*dw;
   end
   xx=[xx;fblt,dbmax,wmax];
%  dbmax specification is 2*n
   if(abs(dbmax-2*n)<0.01) break; end
%  use interval halving to change fblt
   if((dbmax-2*n)<=0) 
      if(flagm<0) df=df/2; end
      fblt=fblt-df;
      if(fblt>0.5) flagp=-1; end
   else
      if(flagp<0) df=df/2; end
      fblt=fblt+df;
      flagm=-1;
   end
end
xx,kc,reset

