mod_reduction.m

很优良的PID控制器设计仿真程序与模型,经过严格检验

%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);
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);
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,'Td',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.Td=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,'Td',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.Td; 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      
   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 v2==2 & v3==0, 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 v2==2 & v3==0, 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 v2>=2, 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];