%mod_reduction is used to administrate the reduction of the plant model.
%
%   mod_reduction(nTask,arg1)
%
%If there is no input argument supplied, then reduced model comparison dialog box 
%will appear which allows the user to compare the behaviours of different reduced 
%models.
%
%If there is only one input argument supplied, then the parameters of the reduction 
%process should be entered, and it will chose properiate algorithm to do reduction.
%
%If there are two input arguments supplied, then 
%    
%   arg1=0 --- setting reduction parameters
%   arg1=1,2 --- perform reduction
%   arg1=3 --- compare reduction
%
%Available lists of functions under this module
%
%  do_reduction -- manipulate the model reduction procedures
%  pade_red -- implements the Pade approximation technique
%  routh_red -- implements the Routh approximation technique
%  dom_mod_red -- implements the dominent mode technique
%  ffpade_red -- implements the FF Pade approximation technique
%  bal_real_red -- implements the balanced rr\ealization technique
%  robust_red -- implements the robust reduction technique
%  opt_reduction -- implements the optimal reduction technique
%  tf_monic -- performs monic transformation to transfer functions
%  tf_revs -- write the TF polynomials in reversed order
%  cmp_reduction -- compare different reduced order models.
%  modcmp_win -- creates a dialog box for model comparison
%  modred_win -- creates a dialog box for reduction parameters
%  opt_func -- defines the optimal criterion for model reduction
%  r_intgl -- evaluate the H_2 norm using recursive algorithm
%  tf_derv -- obtain the first derivative of a TF
%  vec_plus -- adds two polynomials of different lengths
%

%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 G_out=mod_reduction(nTask,arg1,arg2,arg3,arg4,arg5,arg6,arg7)

switch nargin
case 0, modcmp_win;
case 1, 
   g_main=findobj('Tag','CtrlLABMain');
   uu0=get(g_main,'UserData'); g1=get(uu0{1}(1),'UserData');
   if length(g1)>0, modred_win(nTask);
   else, 
      warndlg('Plant model dose not exist, please enter it first!','Model reduction failed'); 
   end   
case 2   
   switch arg1
   case 0, do_reduction;
   case 1
      uu=get(gcf,'UserData'); nTask=get(uu{1}(1),'Value');
      switch nTask
      case {1,3,4}, extra_funs(4,2,'Visible',1:4,5:10);
      case 6, extra_funs(4,2,'Visible',[1:4,7:10],5:6);
      otherwise, extra_funs(4,2,'Visible',5:6,[1:4,7:10]);   
      end
   case 2
      uu=get(gcf,'UserData'); nTask=get(uu{1}(1),'Value');
      switch get(uu{2}(10),'Value')
      case 1, avec(4:5)=[0,0];
      case 2, avec(4:5)=[0,1];
      case 3, avec(4:5)=[1,0]; 
      end; 
      set(uu{2}(8),'String',mat2str(avec));  
   case 3, cmp_reduction;
   case 4
      uu0=get(findobj('Tag','CtrlLABMain'),'UserData');
      g4=get(uu0{1}(4),'UserData'); ModPars=g4{3}; 
      m_list=findobj(gcf,'Tag','MList'); m_list=m_list(end:-1:1);
      keyMod=extra_funs(5,m_list,'Value');
      if length(keyMod)==0, keyMod=n_max; end
      nn=ModPars(keyMod,1); nd=ModPars(keyMod,2); 
      num=ModPars(keyMod,7:7+nd); den=ModPars(keyMod,8+nd:8+2*nd);
      Tau=ModPars(keyMod,9+2*nd); g4{4}=tf(num,den,'Td',Tau);
      set(uu0{1}(4),'UserData',g4);
      proc_model(1,5);   
   end
case {3,4,5,6,7}   
   G_Sys=nTask; 
   nn=arg1(1); nd=arg1(2); avec=arg2; id=arg3(1); ic=arg3(2);
   nPade=arg4(1); Pdly=arg4(2);
   G_out=opt_reduction(G_Sys,nn,nd,avec,id,ic,[nPade,Pdly]);
case 8, G_out=opt_func(nTask,arg1,arg2,arg3,arg4,arg5,arg6,arg7);
end

%------------------------------------------------------------------------------
%do_reduction is the function which calls the corresponding routines to perform
%model reduction.
%------------------------------------------------------------------------------
function do_reduction()
uu=get(gcf,'UserData'); nTask=get(uu{1}(1),'Value'); Pdly=get(uu{1}(2),'Value');
if any([1,3,4,6]==nTask),
   nn=eval(get(uu{2}(2),'String')); nd=eval(get(uu{2}(4),'String'));
   if nTask==6,
      avec=eval(get(uu{2}(8),'String')); ic=get(uu{2}(10),'Value');
   else, ic=0; id=0; end;
else, 
   ic=0; id=0; nd=eval(get(uu{2}(6),'String')); nn=nd-1; 
   if nTask==5, nn=nd; end
end

uu0=get(findobj('Tag','CtrlLABMain'),'UserData');
g1=get(uu0{1}(1),'UserData'); G1=tf(g1{2}); g4=get(uu0{1}(4),'UserData');
Td=g4{1}; nPade=g4{2}; G1.Td=Td; 

if nTask==6
   if abs(avec(2))<1e-5, id=2;
   elseif abs(avec(1))<1e-5, if abs(avec(3))<1e-5, id=1; else, id=3; end
   else, if abs(avec(3))<1e-5, id=4; else, id=5; end
   end
end

%check whether such a reduced order model exists
ModPars=g4{3}; [n,m]=size(ModPars); key1=0;
for i=1:n,
   ii=(ModPars(i,1:6)==[nn,nd,nTask,Pdly id ic]);
   if sum(ii)==6, key1=1; break; end
end

key=1; 
if key1==1
   %if it exists, then load and display it directly
   nn=ModPars(i,1); nd=ModPars(i,2);
   num=ModPars(i,7:7+nd); den=ModPars(i,8+nd:8+2*nd);
   Td=ModPars(i,9+2*nd); G_red=tf(num,den,'Td',Td); 
else
   %if not exist, then perform reduction
   switch nTask
   case 1, G_red=pade_red(G1,nn,nd,[nPade,Pdly]);
   case 2, G_red=routh_red(G1,nd,[nPade,Pdly]);
   case 3, [G_red,key]=dom_mod_red(G1,nn,nd,[nPade,Pdly]);
   case 4, G_red=ffpade_red(G1,nn,nd,[nPade,Pdly]);
   case 5, G_red=bal_real_red(G1,nd,[nPade,Pdly]);
   case 6, 
      g_wait=extra_funs(12,'Performing optimal reduction');
      pause(0.01)
      G_red=opt_reduction(G1,nn,nd,avec,id,ic,[nPade,Pdly]);
      close(g_wait);
   case {9,10}
      if exist('balmr')~=2,
         key=0;
         errordlg('Robust Control Toolbox is not available!','Warning: Reduction failed!');
      else,
         G_red=robust_red(G1,nd,nTask,[nPade,Pdly]);
      end
   otherwise,
      key=0;
      errordlg('Algorithm not available yet!','Warning: Reduction failed!');
   end
   if key==1, 
      %if reduction is success, then write the results back to the system
      g4=get(uu0{1}(4),'UserData'); ModPars=g4{3}; [n1,n2]=size(ModPars);
      [num,den,Td]=tfdata(G_red,'v'); PVec=[nn,nd,nTask,Pdly,id,ic,num,den,Td];
      n3=length(PVec); 
      ModPars=[ModPars,zeros(n1,n3-n2); PVec, zeros(1,n2-n3)];
      % reserve only the last 7 models
      if n1==7, ModPars=ModPars(2:8,:); end
      g4{3}=ModPars; 
   end
end
if key==1
   g4{4}=G_red; set(uu0{1}(4),'UserData',g4); close(gcf); proc_model(1,5); 
   uu=get(gcf,'UserData'); set(uu(6),'Visible','on');
end

%-----------------------------------------------------------------------------
%pade_red implements the Pade approximation technique
%
%   G_red=pade_red(G_Sys,nn,nd,arg1)
%where G_Sys is the original plant model transfer function which possibily have 
%delay terms.  (nn,nd) is the expected order of num and den, respectively. arg1 
%contains delay parameters -- [nPade, Pdly], where nPade is the order of Pade 
%approximation, and Pdly is where delay terms are used in the reduced model.
%G_red is the reduced model which can contain delays.
%------------------------------------------------------------------------------
function G_red=pade_red(G_Sys,nn,nd,arg1)
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
n=length(den)-1; k=length(num)-1; d0=den(n+1); L=nn+nd;
num0=[zeros(1,L+n-k-1) num]/d0; den0=[zeros(1,L) den]/d0;
c=zeros(1,nn+nd); c(1)=num0(n+L);
for i=2:nn+nd+1, c(i)=num0(n+1+L-i)-sum(den0(n-i+2+L:n+L).*c(1:i-1)); end

w=-c(nn+2:nn+nd+1)'; vv=[c(nn+1:-1:1)'; zeros(nd-1-nn,1)];
W=rot90(hankel(c(nn+nd:-1:nn+1),vv)); V=rot90(hankel(c(nn:-1:1)));
x=[1 (W\w)']; dred=x(nd+1:-1:1)/x(nd+1); y=[c(1) x(2:nn+1)*V'+c(2:nn+1)]; nred=y(nn+1:-1:1)/x(nd+1);

if Pdly==0, Td=0; end
G_red=tf(nred,dred,'Td',Td);

%----------------------------------------------------------
%routh_red implements the Routh approximation technique
%
%   G_red=routh_red(G_Sys,nd,arg1)
%descriptions to the I/O arguments can be found in pade_red
%----------------------------------------------------------
function G_red = routh_red(G_Sys,nr,arg1)
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
n0=length(den); n1=length(num); a1=den(n0:-1:1); b1=[num(n1:-1:1) zeros(1,n0-n1-1)];
for k=1:n0-1, 
   k1=k+2; alpha(k)=a1(k)/a1(k+1); beta(k)=b1(k)/a1(k+1);
   for i=k1:2:n0-1, a1(i)=a1(i)-alpha(k)*a1(i+1); b1(i)=b1(i)-beta(k)*a1(i+1); end, 
end
nn=[]; dd=[1]; nn1=beta(1); dd1=[alpha(1), 1]; nred=nn1; dred=dd1;
for i=2:nr,
   nred=[alpha(i)*nn1, beta(i)]; dred=[alpha(i)*dd1, 0]; n0=length(dd); n1=length(dred);
   nred=nred+[zeros(1,n1-n0),nn]; dred=dred+[zeros(1,n1-n0),dd];
   nn=nn1; dd=dd1; nn1=nred; dd1=dred;
end
nred=nred(1:nr); [nred,dred]=tf_revs(nred,dred); [nred,dred]=tf_monic(nred,dred,1);

if Pdly==0, Td=0; end
G_red=tf(nred,dred,'Td',Td);

%----------------------------------------------------------
%dom_mod_red implements the dominent mode technique
%
%   G_red=dom_mod_red(G_Sys,nn,nd,arg1)
%descriptions to the I/O arguments can be found in pade_red
%----------------------------------------------------------
function [G_red,key]=dom_mod_red(G_Sys,nn,nd,arg1)
nPade=arg1(1); Pdly=arg1(2); key=1; [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
n0=length(den); n1=length(num); j=sqrt(-1); r=roots(den); rr=abs(real(r)); ri=imag(r);
[rr,i]=sort(rr); dmode(1:nd)=-rr(1:nd)+ri(i(1:nd))*j; dred=[1 -dmode(1)];
for i=2:nd, dred=conv(dred,[1 -dmode(i)]); end
d0=dred(nd+1); dred=dred(nd+1:-1:1)/d0; num1=[zeros(1,n0-n1) num]; c(1)=num(n1)/den(n0);
for i=2:nd, c0=0;
   for j=1:i-1, 
      i0=n0-j; if i0>0, c0=c0+den(i0)*c(i-j); end
   end
   i0=n0+1-i; if i0>0, c(i)=(num1(i0)-c0)/den(n0);
   else, c(i)=-c0/den(n0);
end, end
v0=dred(2:nn+1); 
for i=1:nn, for j=1:nn, if i+1-j>0, v(i,j)=c(i+1-j); end, end, end
nred=[c(1) v0*v'];
for i=2:nn+1, nred(i)=nred(i)+c(i); end
[nred,dred]=tf_revs(nred,dred); [nred,dred]=tf_monic(nred,dred,1);
if (abs(sum(imag(nred)))+abs(sum(imag(dred)))<10000*eps), 
   key=0; 
   errordlg('Cannot find the specified dominent modes!','Warning: Reduction failed!');
else,
   if Pdly==0, Td=0; end
   G_red=tf(nred,dred,'Td',Td);
end

%------------------------------------------------------------------------
%ffpade_red implements the frequency fitting Pade approximation technique
%
%   G_red=ffpade_red(G_Sys,nn,nd,arg1)
%descriptions to the I/O arguments can be found in pade_red
%------------------------------------------------------------------------
function G_red=ffpade_red(G_Sys,nn,nd,arg1)
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
w=logspace(-1,3); [m,p]=bode(num,den,w); graf_tool(1);
semilogx(w,20*log10(m)); [xx,yy]=ginput(nd-1);
np=nd-1;xx=[1;4]; wp=xx'; [mp,pp]=bode(num,den,wp); pp=pp*pi/180;

n0=length(den); n1=length(num); num1=[zeros(1,n0-n1) num]; c(1)=num(n1)/den(n0);
aa=den(n0:-1:1); bb=num1(n0:-1:1); [bb,aa]=tf_monic(bb,aa);
for i=2:nn+nd+1,
   for j=1:n0-1, cc(j)=bb(1)*aa(j+1)-bb(j+1); end
   bb=[cc zeros(1,n0-length(cc))]; cc=[]; c(i)=(-1)^(i-1)*bb(1);  
end
w11=eye(nn+1); w21=zeros(nd,nn+1); w12=zeros(1,nd); temp=w12;
for i=1:nn, temp=[-c(i), temp(1:nd-1)]; w12=[w12; temp]; end
w22=[c(nn+1), -temp(1:nd-1)]; temp=w22;
for i=1:nd-1, temp=[c(nn+1+i), temp(1:nd-1)]; w22=[w22; temp]; end
W=[w11,w12; w21,w22]; V=[c(1:nn+1)'; -c(nn+2:nn+nd+1)'];
%  ** replacing 2*np equation by Frequency Fit
for i=1:np,
   i1=nd+nn+2-2*i; i2=i1+1; W(i1,1)=-sin(pp(i)); W(i2,1)=cos(pp(i)); xp=wp(i)*wp(i);
   if nn~=0, W(i1,2)=wp(i)*cos(pp(i)); W(i2,2)=wp(i)*sin(pp(i)); end
   for j=3:nn+1, W(i1,j)=-xp*W(i1,j-2); W(i2,j)=-xp*W(i2,j-2); end
   W(i1,nn+2)=-wp(i)*mp(i); W(i1,nn+3)=0; W(i2,nn+2)=0; W(i2,nn+3)=xp*mp(i);
   for j=nn+4:nn+nd+1, W(i1,j)=-xp*W(i1,j-2); W(i2,j)=xp*W(i2,j-2); end
   V(i1)=0; V(i2)=mp(i);
end
theta=inv(W)*V; nred=theta(nn+1:-1:1)'/theta(nd+nn+1); dred=[theta(nd+nn+1:-1:nn+2)' 1]/theta(nd+nn+1);

if Pdly==0, Td=0; end
G_red=tf(nred,dred,'Td',Td);

%----------------------------------------------------------
%bal_real_red implements the balanced realization technique
%
%   G_red=bal_real_red(G_Sys,nr,arg1)
%descriptions to the I/O arguments can be found in pade_red
%----------------------------------------------------------
function G_red=bal_real_red(G_Sys,nr,arg1)
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
G_Sys=ss(tf(num,den)); G_red=minreal(G_Sys); [G_red,g]=balreal(G_red);
nn=length(g); aa=nn:-1:1; [g,ii]=sort(g); aa=aa(ii); elim=aa(nr+1:nn);
G_red=modred(G_red,elim); 
if Pdly==0, Td=0; end
G_red=tf(G_red); G_red.Td=Td; 

%--------------------------------------------------------------------------------
%robust_red implements the robust reduction technique
%
%   robust_red(G_Sys,nd,nTask,arg1)
%where nTask=1 for optimal Hankel approximation, while 2 for Schur approximation. 
%Other descriptions to the I/O arguments can be found in pade_red
%--------------------------------------------------------------------------------
function G_red=robust_red(G_Sys,nd,nTask,arg1);
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); G_Sys=G_Sys*tf(nP,dP); end
ss_G=ss(G_Sys); ii=find(real(eig(ss_G))>=0);
if nTask==9, [am,bm,cm,dm,totbnd,hsv]=balmr(ss_G.a,ss_G.b,ss_G.c,ss_G.d,1,nd);
else, [am,bm,cm,dm,totbnd,hsv]=schmr(ss_G.a,ss_G.b,ss_G.c,ss_G.d,1,nd); end
G_red=tf(ss(am,bm,cm,dm)); 
num=G_red.num{1}; ii=find(abs(num)>eps); G_red.num{1}=num(ii(1):end);
if Pdly==0, Td=0; end
G_red.Td=Td;

%----------------------------------------------------------------
%opt_reduction implements the optimal reduction technique
%
%   G_red=opt_reduction(G_Sys,nn,nd,avec,id,ic,arg1)