mod_reduction.m

非常不错的一本资料

%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: [e-mail]
%
%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,'InputDelay',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.InputDelay=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,'InputDelay',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,'InputDelay',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,'InputDelay',Td);