mod_reduction.m

非常不错的一本资料

   end
   set(h,'Color',lstColorMap(i+1,:));   
   if i==0
      h_axes=extra_funs(2); for j=1:length(h_axes), axes(h_axes(j)); hold on; end   
   end
end   
for i=1:length(h_axes), axes(h_axes(i)); hold off; end   

%--------------------------------------------------------
%modcmp_win creates a dialog box for model comparison
%--------------------------------------------------------
function modcmp_win()
g_cmp=findobj('Name','Reduced Model Comparison');
if length(g_cmp)==0
   g_cmp=figure('Units','normalized','Position',[0.18625 0.248 0.425 0.367],...
      'NumberTitle','off','Name','Reduced Model Comparison','Tag','CtrlLABExtras',...
      'MenuBar','none','Color',0.8*[1,1,1],'Resize','off');
   extra_funs(1); hh=extra_funs(10,[0.05,0.05],[0.56,0.9]);
   uicontrol('Style','Text','String','Select a Model',...
      'Units','normalized','Position',[0.07,0.87,0.30,0.06],'BackgroundColor',0.8*[1,1,1]);
   uicontrol('Style','PushButton','String','Compare',...
      'Units','normalized','Position',[0.80,0.81,0.18,0.1],'CallBack','mod_reduction(0,3); ');
   uicontrol('Style','PushButton','String','Show',...
      'Units','normalized','Position',[0.60,0.81,0.18,0.1],'CallBack','mod_reduction(0,4);');
   uicontrol('Style','PushButton','String','Close',...
      'Units','normalized','Position',[0.60,0.68,0.18,0.1],'CallBack','close(gcf);');
   uicontrol('Style','PushButton','String','Help',...
      'Units','normalized','Position',[0.80,0.68,0.18,0.1],'CallBack','clab_help(26);');
   uicontrol('Style','PushButton','String','Refresh List',...
      'Units','normalized','Position',[0.65,0.55,0.28,0.1],'CallBack','mod_reduction;');
   add_modlist;
   extra_funs(10,[0.60,0.05],[0.97,0.48]);
   uicontrol('Style','Text','String','Assign Specifications',...
      'Units','normalized','Position',[0.62,0.45,0.31,0.06],'BackgroundColor',0.8*[1,1,1]);
   str=str2mat('Bode diagram','Nyquist plot','Nichols chart',...
        'Step response','Impulse response','Criterion');
   uicontrol('Tag','List','Style','Popupmenu','String',str,'Value',4,...
      'Units','normalized','Position',[0.63,0.32,0.30,0.07],'BackgroundColor',[1,1,1]);
   uu=get(findobj('Tag','CtrlLABMain'),'UserData');
   g1=get(uu{1}(1),'UserData');
   [v,d]=version; v1=eval(v(1)); v2=eval(v(3)); v3=eval(v(5));
   if (v1==5 & v2==2 & v3==0) | v1==6, strCheck='ToggleButton'; else, strCheck='CheckBox'; end 
   uicontrol('Tag','SLnk','Style',strCheck,'String','Use SIMULINK','Value',0,...
      'Units','normalized','Position',[0.66,0.14,0.26,0.08],...
      'BackgroundColor',0.8*[1,1,1],'Enable',extra_funs(6,g1{1}==4));
else, 
   h_win=figure(g_cmp); delete(findobj(gcf,'Tag','MList')); add_modlist;
end

%--------------------------------------------------------
%modred_win creates a dialog box for reduction parameters
%--------------------------------------------------------
function modred_win(nTask)
g_red=findobj('Name','Model Reduction Parameters');
if length(g_red)==0
   avec=[1,0,0,0,0];
   g_red=figure('Units','normalized','Position',[0.18625 0.248 0.425 0.283],...
      'NumberTitle','off','Name','Model Reduction Parameters','Tag','CtrlLABExtras',...
      'MenuBar','none','Color',0.8*[1,1,1],'Resize','off');
   extra_funs(1); kDly=0; 
   display_str(0.04,0.9,'Select a Method',[0,0,0],'on',9);
   str=str2mat('Pade approximation','Routh''s method','Dominant modes',...
        'FF-Pade method','Balanced realisation','Optimal reduction',...
        'Modal reduction','Aggregation','Optimal Hankel Approxation','Schur reduction');
   h_alg(1)=uicontrol('Style','Popupmenu','String',str,'Value',nTask,...
      'Units','normalized','Position',[0.1,0.72,0.45,0.11],...
      'BackgroundColor',[1,1,1],'CallBack','mod_reduction(0,1);');
   [v,d]=version; v1=eval(v(1)); v2=eval(v(3)); v3=eval(v(5));
   if (v1==5 & v2==2 & v3==0) | v1==6, strCheck='ToggleButton'; else, strCheck='CheckBox'; end 
   h_alg(2)=uicontrol('Style',strCheck,'String','With Delay',...
      'Units','normalized','Position',[0.1,0.55,0.30,0.11],...
      'BackgroundColor',0.8*[1,1,1],'CallBack','extra_funs(3,gco,''Value'');');
   extra_funs(10,[0.0386, 0.0509],[0.55,0.4281]);
   uicontrol('Style','Text','String','Expected Reduction Order',...
      'Units','normalized','Position',[0.08,0.36,0.38,0.10],'Backgroundcolor',0.8*[1,1,1]);
   [xL,h_red(1)]=display_str(0.08,0.31,'Numerator Order',[0,0,0],'on',9);
   h_red(2)=uicontrol('Style','Edit','String','1',...
      'Units','normalized','Position',[0.4,0.25,0.1,0.11],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   [xL,h_red(3)]=display_str(0.08,0.16,'Denominator Order',[0,0,0],'on',9);
   h_red(4)=uicontrol('Style','Edit','String','2',...
      'Units','normalized','Position',[0.4,0.11,0.1,0.11],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1]);
   [xL,h_red(5)]=display_str(0.08,0.29,'Reduction Order',[0,0,0],'off',9);
   h_red(6)=uicontrol('Style','Edit','String','2',...
      'Units','normalized','Position',[0.4,0.24,0.1,0.11],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1],'Visible','off');
   [xL,h_red(7)]=display_str(0.58,0.24,'Weighting Vector',[0,0,0],'off',9);
   h_red(8)=uicontrol('Style','Edit','String','[1,0,0,0,0]',...
      'Units','normalized','Position',[0.60,0.08,0.35,0.11],...
      'HorizontalAlignment','left','BackgroundColor',[1,1,1],'Visible','off');
   [xL,h_red(9)]=display_str(0.58,0.49,'Criterion',[0,0,0],'off',9);
   h_red(10)=uicontrol('Style','PopupMenu','Value',1,'Visible','off',...
      'String','ISE criterion |TISE criterion |Exponential',...
      'Units','normalized','Position',[0.60,0.32,0.34,0.11],...
      'BackgroundColor',[1,1,1],'CallBack','mod_reduction(0,2);');
   set(gcf,'UserData',{h_alg,h_red});
   uicontrol('Style','Pushbutton','String','Reduce',...
      'Units','normalized','Position',[0.80,0.83,0.17,0.14],'CallBack','mod_reduction(0,0);');
   uicontrol('Style','Pushbutton','String','Cancel',...
      'Units','normalized','Position',[0.80,0.67,0.17,0.14],'CallBack','close(gcf);');
   uicontrol('Style','Pushbutton','String','Help',...
      'Units','normalized','Position',[0.80,0.51,0.17,0.14],'Callback','clab_help(4);');
else, figure(g_red); end
mod_reduction(0,1);

%------------------------------------------------------------------------------
%opt_func defines the optimal criterion for model reduction of control systems.
%
%  y=opt_func(x,b,a,nvec,Tdly,nPad,Pdly,avec)
%where
%  y -- the actual criterion calculated
%  x -- the vector of parameters to be optimized
%  (b,a) -- the transfer function of the original system
%  nvec -- the vector contains information for reduction
%  Tdly -- delay constant of the original system
%  nPad -- order of Pade approximation
%  Pdly -- key for whether delay is expected in reduction
%  avec -- the weighting vector for optimal reduction
%------------------------------------------------------------------------------
function y=opt_func(x,b,a,nvec,Tdly,nPad,Pdly,avec)
ff0=1e10;

nn=nvec(1); nd=nvec(2); id=nvec(3); n0=nvec(4); ic=nvec(5); key=nvec(6);
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)*b(end)/a(end); end

a0=conv(a,alph); b0=conv(b,alph); c0=conv(beta,a);
if Pdly==1,
   tau=x(end); [nP,dP]=pade(tau,nPad); a0=conv(a0,dP); b0=conv(b0,dP); c0=conv(c0,nP); 
end

d0=vec_plus(c0,-b0); c0=a0; n_max=length(d0);
switch id
case 2, d0=[0,d0(1:n_max)];
case 3, c0=conv(c0,[a3,1]); d0=[0,d0(1:n_max)];
case 4, d0=[0 conv(d0,[a2,a1])];
case 5, c0=conv(a0,[a3,a1]); d0=[0 conv([a2+a1*a3,a1],d0)];
end
ct=c0; dt=d0;

if ic==2, for i=1:avec(5), [b0,a0]=tf_derv(d0,c0); c0=a0; d0=b0; end
elseif ic==3, c0=polyexpd(c0,avec(4)); d0=polyexpd(d0,avec(4)); end

[y,ierr]=r_intgl(c0,d0);
if ic==4, [y1,ierr]=r_intgl(ct,dt); y=y+y1; end
if ierr==1, y=10*ff0; else, ff0=y; end

%--------------------------------------------------------------------------------
%r_intgl is used to evaluate the H_2 norm of a transfer function model, using the 
%Astrom recursive algorithm.
%
%   [v,ierr]=r_intgl(a,b)
%where
%   (a,b) --- den and num of the transfer function
%   v --- the H_2 norm
%   ierr --- 1 if there are error or unstable
%--------------------------------------------------------------------------------
function [v,ierr]=r_intgl(a,b)
ierr=0; v=0; 
n=length(a); nb=length(b);
if n<nb,
   warndlg('System not proper!','Warning: Reduction failed!');
   ierr=1; return
elseif n==nb
   if abs(b(1))>eps,
      warndlg('System not strictly proper!','Warning: Reduction failed!');
      ierr=1; return
   else, a1=a; b1=b(2:nb); end
else, a1=a; b1=[zeros(1,n-nb-1) b]; end
if (a1(1)<=eps), ierr=1;  
else
   for k=1:n-1
      if (a1(k+1)<=eps), ierr=1; return
      else,
         alpha=a1(k)/a1(k+1); beta=b1(k)/a1(k+1); v=v+beta*beta/alpha; k1=k+2;
         for i=k1:2:n-1, a1(i)=a1(i)-alpha*a1(i+1); b1(i)=b1(i)-beta*a1(i+1); end
   end, end,
   v=0.5*v;
end

%----------------------------------------------------------------------------
%tf_derv is used to obtain the first derivative of a given transfer function.
%
%   [e,f]=tf_derv(b,a)
%where 
%  (b,a) --- the num, den of the original TF
%  (e,f) --- the num, den of the derivative of the TF
%----------------------------------------------------------------------------
function [e,f]=tf_derv(b,a)
f=conv(a,a); e1=conv((length(b)-1:-1:1).*b(1:end-1),a);
e2=conv((length(a)-1:-1:1).*a(1:end-1),b); e=vec_plus(e1,-e2);

%----------------------------------------------------------------------------
%add_modlist is used to add the obtained reduced order model to the existing
%model list so that they can be compared.
%----------------------------------------------------------------------------
function add_modlist()

%define the reduction name abbreviate form
MRedName=str2mat('Pade approximation','Routh approximation','Dominant mode',...
   'FF-Pade approximation','Balance Realization','Optimal reductition',...
   'Modal reduction','Aggregation','Optimal Hankel','Optimal Schur');
uu=get(findobj('Tag','CtrlLABMain'),'UserData');
g4=get(uu{1}(4),'UserData'); ModPars=g4{3}; [nModels,m]=size(ModPars); chkModels=[];
[v,d]=version; v1=eval(v(1)); v2=eval(v(3));
if (v1==5 & v2>=2) | v1==6, strCheck='ToggleButton'; else, strCheck='CheckBox'; end 

%set a check box for each of the reduced models
X0=0.1; Y0=0.75;
for i=1:nModels
   str=[MRedName(ModPars(i,3),:),int2str(ModPars(i,1)),'/' int2str(ModPars(i,2))];
   uicontrol('Tag','MList','Style',strCheck,'String',str,'Value',0,...
      'Units','normalized','Position',[X0,Y0,0.45,0.09],'BackgroundColor',0.8*[1,1,1]);
   Y0=Y0-0.10;
end

%----------------------------------------------------------------------------
%vec_plus is used to add two polynomials together
%    
%        e=vec_plus(e1,e2)
%----------------------------------------------------------------------------
function e=vec_plus(e1,e2)
L0=length(e1)-length(e2); e=[zeros(1,-L0) e1]+[zeros(1,L0) e2];