mod_reduction.m

非常不错的一本资料

%----------------------------------------------------------
%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,'InputDelay',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,'InputDelay',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.InputDelay=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.InputDelay=Td;

%----------------------------------------------------------------
%opt_reduction implements the optimal reduction technique
%
%   G_red=opt_reduction(G_Sys,nn,nd,avec,id,ic,arg1)
%where 
%   avec -- the weighting vector for optimal reduction
%   (id,ic) -- criterion weighting and performance weighting
%Other descriptions to the I/O arguments can be found in pade_red
%----------------------------------------------------------------
function G_red=opt_reduction(G_Sys,nn,nd,avec,id,ic,arg1)
nPade=arg1(1); Pdly=arg1(2); [num,den,Td]=tfdata(G_Sys,'v');
if num(1)==0, key=0; num=num(2:end); else, key=1; end
if abs(Td)>0 & Pdly==0, [nP,dP]=pade(Td,nPade); num=conv(num,nP); den=conv(den,dP); end
G_routh=routh_red(G_Sys,nd,[3,0]); beta=G_routh.num{1}(nd+1-nn:nd+1); alph=G_routh.den{1};
n0=1; if id==1 | key==1, n0=0; end
x=[beta(1:nn+1-n0),alph(2:nd+1)];
Tau=1.5*Td; if abs(Tau)<1e-5, Tau=0.5; end
if Pdly==1, x=[x,Tau]; end
y=opt_func(x,num,den,[nn,nd,id,n0,ic,key],Td,nPade,Pdly,avec);
[v,d]=version; v1=eval(v(1)); v2=eval(v(3));
if (v1==5 & v2==3) | v1==6
   fopts=optimset; fopts.TolX=1e-4; fopts.TolFun=1e-4; fopts.MaxIter=500;
   x=fminsearch('mod_reduction',x,fopts,num,den,[nn,nd,id,n0,ic,key],Td,nPade,Pdly,avec);
else
   x=fmins('mod_reduction',x,[0,1e-4,1e-4,zeros(1,10),500],[],num,den,[nn,nd,id,n0,ic,key],Td,nPade,Pdly,avec);
end
alph=[1,x(nn+2-n0:nn+1-n0+nd)]; beta=x(1:nn+1);
if id~=1 & key==0, beta(nn+1)=alph(end)*num(end)/den(end); end

if Pdly==1, Tau=x(end)+Td; else, Tau=0; end
G_red=tf(beta,alph,'InputDelay',Tau); 

%-------------------------------------------------------------------------------
%tf_monic is a function which can be used to perform monic transformation to the
%transfer function of a control system.
%
%  [NN, DD]=tf_monic(NUM, DEN)
%where
%   NN, DD are the num and den after monic transformation
%   NUM, DEN are the num and den of original system
%-------------------------------------------------------------------------------
function [nn,dd]=tf_monic(num,den,key)
if nargin==2
   nb=length(den); nn=num./den(1); i=length(den):-1:1; dd(i)=den(i)/den(1); 
   nn=[zeros(1,length(dd)-length(nn)) nn];
elseif key==1, nn=num/den(1); dd=den/den(1);
else, nn=num/den(end); dd=den/den(end); end

%--------------------------------------------------------------------------
%tf_revs is a function which can be used to transform the transfer function
%                 B(1) s^M + ... + B(M+1) 
%       ^G(s) = -----------------------------
%                 A(1) s^N + ... + A(N+1)
%into its equivelent form of
%               B(1) + B(2) s + ... + B(M+1) s^M
%      G(s) = ------------------------------------
%              A(1) + A(2) s + ... + A(N+1) s^N
%or vise versa.
%    [N, D] = tf_revs(NUM, DEN)
%     N, D are the num and DEN of G(s) 
%     NUM, DEN are the num and den of ^G(s)
%---------------------------------------------------------------------------
function [n,d]=tf_revs(num,den)
d=den(end:-1:1); n=num(end:-1:1);

%----------------------------------------------------------------
%cmp_reduction is used to compare different reduced order models.
%----------------------------------------------------------------
function cmp_reduction()
g0=gcf;
g_Slnk=findobj(g0,'Tag','SLnk'); g_list=findobj(g0,'Tag','List');
nTask=get(g_list,'Value'); b_Simulink=get(g_Slnk,'Value');
uu0=get(findobj('Tag','CtrlLABMain'),'UserData');
if b_Simulink==1
   g1=get(uu0{1}(1),'UserData');   
   if g1{1}==4, sim_model=g1{3};
   else, b_Simulink=0; set(g_Slnk,'Value',0); end   
end
g4=get(uu0{1}(4),'UserData'); ModPars=g4{3};
g1=get(uu0{1}(1),'UserData'); G_Sys=g1{2}; G_Sys.InputDelay=g4{1};
ii=[]; 
m_list=findobj(g0,'Tag','MList'); m_list=m_list(end:-1:1);
for i=1:length(m_list), if get(m_list(i),'Value')==1, ii=[ii,i]; end, end
lstColorMap=[0         0    1.0000
             0    0.5000         0
        1.0000         0         0
             0    0.7500    0.7500
        0.7500         0    0.7500
        0.7500    0.7500         0
        0.2500    0.2500    0.2500
        0.5000    0.2500    0.7500];
switch nTask     
case 3, graf_tool(3); delete(gca); hold off; ngrid; 
case 4, tr_cmd='step';
case 5, tr_cmd='impulse';   
end
for i=0:length(ii)
   if i==0, G_red=G_Sys; 
   else
      i_now=ii(i); nn=ModPars(i_now,1); nd=ModPars(i_now,2); 
      num=ModPars(i_now,7:7+nd); den=ModPars(i_now,8+nd:8+2*nd);
      Tau=ModPars(i_now,9+2*nd); G_red=tf(num,den,'InputDelay',Tau); 
   end
   switch nTask
   case 1,
      graf_tool(1); [mag,phi,w]=bode(G_red);
      subplot(211), h(1)=semilogx(w,20*log10([mag(:)]));
      subplot(212), h(2)=semilogx(w,[phi(:)]);
   case 2
      graf_tool(2); [re,im,w]=nyquist(G_red); h=plot([re(:)],[im(:)]);
   case 3
      graf_tool(3); [mag,phi,w]=nichols(G_red); h=plot([phi(:)],20*log10([mag(:)]));
   case {4,5}
      graf_tool(nTask+2); eval(['[y,t]=' tr_cmd '(G_red); ']); 
      if i==0 & b_Simulink==1
         Td=G_Sys.InputDelay; t_end=t(end);
         if nTask==4, ut=[0,1; t_end,1]; 
         else, ut=[0,1000; 0.001,1000; t_end,0]; end  
         [t,x,y]=sim(sim_model,[0,t_end],[],ut); ii=find(t>=Td); 
         y=[zeros(length(t)-length(ii),1); y]; y=y(1:length(t)); h=plot(t,y);
      else, eval(['[y,t]=' tr_cmd '(G_red); ']); h=plot(t,y); end