sys_analysis.m

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

         if g_reltn==3, G_overall=1-G_overall; end
      end  
      [mag,pha]=bode(G_overall,w); mag=mag(:); pha=pha(:);
   end      
end
switch nTask
case 14, [xpos,ypos]=bode_asymp(G_overall,w); arg4=[xpos; ypos];
case 15, mag=2*atan(mag)/pi;
end

%---------------------------------------------------------------------------
%time_resp is used to evaluate the time domain response of the given system.
%---------------------------------------------------------------------------
function [y,t]=time_resp(k_Extras,nCtrl,nTask,t)

y=[]; if nargout==2, t=[]; end
[G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_sigs]=plt_common(k_Extras,nCtrl);
if length(G_forward)==0, return; end

if nTask==6, plt_cmd='step'; else, plt_cmd='impulse'; end
if nargin==3, %selecting the time vector t in an automatic way
   eval(['[a,t]=' plt_cmd '(G_Sys);']) 
end

%check whether there is delay and whether Pade appr is needed

if k_Extras(1)==1
   str=str2mat('This function is not available yet!  Please check off',...
      'the SIMULINK selction and try again!');
   warndlg(str,'Cannot use SIMULINK!'); 
else   
   if g_loop==1
      %processing open-loop
      if k_Extras(2)==1 & Td>0, G_forward=G_forward*GP; 
      elseif Td>0, G_forward.Td=Td; end
      G_overall=G_forward;
   else   
      %processing closed-loop
      G_forward=G_forward*GP;
      switch g_sigs
      case 1, %closed-loop error signal
         G_overall=feedback(1,G_forward*H_Sys*GP); 
      case 2, %evaluate the actuating signal
         G_overall=feedback(Gc_Sys,G_forward*H_Sys*GP); 
      case 3, %closed loop I/O
         G_overall=feedback(G_forward,H_Sys*GP); G_overall.Td=Td;
      end  
   end
   if nargin==3, eval(['[y,t]=' plt_cmd '(G_overall);']);   
   else, eval(['y=' plt_cmd '(G_overall,t);']); end
end

%---------------------------------------------------------------------------
%plt_common is used to extract common setting in system analysis.
%
%  [G_Sys,Gc_Sys,H_Sys,Td,nPade,g_loop,g_sel]=plt_common(key)
%where 
%  key = 1 for time response, else for frequency response
%  G_Sys, Gc_Sys, H_Sys, Td and nPade are models extracted
%  g_loop -- the loop information, i.e., closed-loop via open-loop
%  g_sel -- selection signal, for freq response, select the transfer 
%          function, else select the signal.
%---------------------------------------------------------------------------
function [G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_sel]=plt_common(k_Pade,nCtrl)
if length(k_Pade)>1, if k_Pade(1)==1; 
   %direct inclusion if SIMULINK blocks are not enabled now   
   G_forward=[]; G_Sys=[]; Gc_Sys=[]; H_Sys=[]; Td=[]; 
   nPade=[]; GP=[]; g_loop=[]; g_sel=[];
   str=str2mat('This function is not available yet!  Please check off',...
      'the SIMULINK selction and try again!');
   warndlg(str,'Cannot use SIMULINK!'); 
   return; 
end, end

figure(findobj('Tag','CtrlLABMain')); uu0=get(gcf,'UserData');
g_loop=extra_funs(5,4,'Checked',24:25); 

%use I/O relationship/y(t) only
if g_loop==1, extra_funs(4,4,'Checked',30:31,[28:29,32:33]); end   

if length(k_Pade)>1, %set signal selection
   k_Pade=k_Pade(2); g_sel=extra_funs(5,4,'Checked',28:30);
else, g_sel=extra_funs(5,4,'Checked',31:33); %set TF relationsship
end

%get all the sub-models
g1=get(uu0{1}(1),'UserData'); G_Sys=g1{2}; 
g2=get(uu0{1}(2),'UserData'); if length(g2)>0, Gc_Sys=g2{2}; else, Gc_Sys=1; end
g3=get(uu0{1}(3),'UserData'); if length(g3)>0, H_Sys=g3{2}; else, H_Sys=1; end
g4=get(uu0{1}(4),'UserData'); Td=g4{1}; nPade=g4{2};

if k_Pade==1 & Td>0, [nP,dP]=pade(Td,nPade); GP=tf(nP,dP);
else, GP=1; end

G_forward=G_Sys; if nCtrl==1, G_forward=Gc_Sys*G_Sys; end

%----------------------------------------------------------------------
%make_model prepares models for further analysis tasks
%----------------------------------------------------------------------
function [G_overall,G_pade]=make_model(nCtrl,arg1)
[G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_reltn]=plt_common(1,nCtrl);
if nargin==2, g_loop=2; r_reltn=3; end
if g_loop==1
   %processing open-loop
   if Td>0, G_forward.Td=Td; end
   G_overall=G_forward; 
else   
   %processing closed-loop
   G_forward=G_forward*GP;
   switch g_reltn
   case 1, G_overall=feedback(G_forward,H_Sys*GP); %closed loop I/O
   case {2,3}, G_overall=feedback(1,G_forward*H_Sys*GP); %closed-loop sensitivity & comp sens
      %evaluate complementary sensitivity function
      if g_reltn==3, G_overall=1-G_overall; end
   end  
end
G_pade=G_overall; G_pade.Td=0;
if Td>0, [nP,dP]=pade(Td,nPade); G_pade=G_pade*tf(nP,dP); end

%----------------------------------------------------------------------
%tim_moments evaluates the time moments of the specified system
%----------------------------------------------------------------------
function M=tim_moments(G,k)
G=ss(G); C=G.c; B=G.b; iA=inv(G.a); iA1=iA; M=[];
for i=1:k, 
   mm=-C*iA1*B; sgn='+'; if mm<0, sgn=''; end
   M=[M, sgn display_str(mm) 's^{' int2str(i) '}'];
   iA1=iA*iA1; 
end
M=[M '+o(s^{' int2str(k+1) '})'];

%----------------------------------------------------------------------
%markov_paras evaluates the Markov parameters of the specified system
%----------------------------------------------------------------------
function M=markov_paras(G,k)
G=ss(G); A=G.a; B=G.b; C=G.c; D=G.d; 
A1=A; M=display_str(C*B+D);
for i=1:k-1, 
   mm=C*A1*B; sgn='+'; if mm<0, sgn=''; end
   M=[M, sgn display_str(mm) 's^{-' int2str(i) '}'];
   A1=A*A1; 
end
M=[M '+o(s^{-' int2str(k) '})'];

%--------------------------------------------------------------------------------
%asol_disp displays in styled format the analytical solutions of the LTI systems
%--------------------------------------------------------------------------------
function [xL,t_vec]=asol_disp(sys_mod,key,xL1,yL1,nCol);
sys_mod=tf(sys_mod); key0=1; 
if strcmp(key,'on'), sys_mod.den={[sys_mod.den{1},0]}; end
if nargin==4, 
   t_vec=xL1; key0=0; if length(t_vec)==0, t_vec=0:.01:1; end; xL=zeros(size(t_vec)); 
else, xL=xL1; end

[r,p,k]=p_frac(sys_mod); Str_s=[]; pp=1.2345678; kk=0;
for i=1:length(r)
   g_i=imag(p(i)); ss1=[]; if real(r(i))>0 & i~=1, ss1='+'; end,
   if g_i>=0
      ss2='-'; if real(p(i))>0, ss2=''; end, ss_x=display_str(r(i));
      if abs(p(i)-pp)<1e-5*abs(pp), kk=kk+1; else, kk=0; pp=p(i); end   
      if ~strcmp(ss_x,'0'), 
         if (strcmp(ss_x,'1')& abs(real(p(i)))>1e-12 &i==length(r)), ss_x=''; 
         elseif strcmp(ss_x,'-1'), ss_x='-'; end
         Str_s=[Str_s,ss1,ss_x]; 
         if abs(real(p(i)))>1e-12, 
            if kk==0, ss0=''; elseif kk==1, ss0='t'; else, ss0=['t^{' int2str(kk) '}']; end   
            Str_s=[Str_s,ss0,'e^{' ss2, display_str(abs(real(p(i))),[]), 't}'];
         end
         if g_i>1e-10
            ss3='+'; if r(i+1)<0, ss3=''; end
            Str_s=[Str_s,'sin(' display_str(g_i,[]),'t',ss3,num2str(180*r(i+1)/pi),'^o)'];
         end, 
      end   
   end
end

%show analytical solutions
xL=display_str(xL1,yL1,Str_s,nCol); 

%--------------------------------------------------------------------------------
%p_frac performs extra trianglar manipulations on the results of residue function
%--------------------------------------------------------------------------------
function [R,P,K]=p_frac(G)
[num,den]=tfdata(G,'v'); [R,P,K]=residue(num,den); 
for i=1:length(R),
   if imag(P(i))>eps
      a=real(R(i)); b=imag(R(i)); R(i)=-2*sqrt(a^2+b^2); R(i+1)=atan2(-a,b);
   end
end
R=real(R);

%--------------------------------------------------------------------------------
%sim_pars creates an extra menu item for system simulation parameters using 
%SIMULINK program.
%--------------------------------------------------------------------------------
function sim_pars()
g_graf=gcf; g_SimPars=findobj('Tag','CtrlLABSimPars'); 
if length(g_SimPars)==0
   uu=get(g_graf,'UserData'); ww=uu{6};
   g_SimPars=figure('Units','normalized','Position',[0.1425 0.132 0.5 0.41],...
      'MenuBar','none','Color',0.8*[1,1,1],'Resize','off','Tag','CtrlLABSimPars',...
      'NumberTitle','off','Name','Simulation Parameters Setting ...','Resize','off');
   extra_funs(1); extra_funs(10,[0.035,0.58],[0.58,0.92]);
   uicontrol('Style','Text','Unit','normalized','Position',[0.055,0.85,0.3,0.08],...
      'String','Simulation Algorithm','BackgroundColor',0.8*[1,1,1]);
   hSimAlg=uicontrol('Style','Popupmenu','Unit','normalized','Position',[0.115,0.75,0.35,0.08],...
      'BackgroundColor',[1,1,1],'Value',1,...
      'String','ode45 (defaults) | ode45 | ode23 | ode113 | ode15s | ode23s');
   extra_funs(10,[0.035,0.05],[0.58,0.50]);
   uicontrol('Style','Text','Unit','normalized','Position',[0.065,0.46,0.33,0.07],...
      'String','Simulation Parameters','BackgroundColor',0.8*[1,1,1]);
   display_str(0.075,0.40,'Terminate Time',[0,0,0],'on',9);
   h(9)=uicontrol('Style','Edit','Unit','normalized','Position',[0.315 0.35 0.15 0.07],...
      'String',num2str(ww(6)),'BackgroundColor',[1,1,1]);
   [xL,h(1)]=display_str(0.075,0.31,'Min Step',[0,0,0],'on',9);
   [xL,h(3)]=display_str(0.075,0.21,'Max Step',[0,0,0],'on',9);
   h(2)=uicontrol('Style','Edit','Unit','normalized','Position',[0.315 0.26 0.15 0.07],...
      'BackgroundColor',[1,1,1],'String',num2str(ww(2)));
   h(4)=uicontrol('Style','Edit','Unit','normalized','Position',[0.315 0.17 0.15 0.07],...
      'BackgroundColor',[1,1,1],'String',num2str(ww(3)));
   [xL,h(5)]=display_str(0.075,0.12,'Tolerant Error',[0,0,0],'on',9);
   h(6)=uicontrol('Style','Edit','Unit','normalized','Position',[0.315 0.08 0.15 0.07],...
      'BackgroundColor',[1,1,1],'String',num2str(ww(4)));
   [xL,h(7)]=display_str(0.075,0.24,'Step Size',[0,0,0],'on',9);
   h(8)=uicontrol('Style','Edit','Unit','normalized','Position',[0.315 0.22 0.12 0.07],...
      'BackgroundColor',[1,1,1],'String',num2str(ww(5)));
   [v,d]=version; v1=eval(v(1)); v2=eval(v(3)); v3=eval(v(5));
   if v2==2 & v3==0, strRadio='ToggleButton'; strCheck='ToggleButton';
   else, strRadio='RadioButton'; strCheck='CheckBox'; end 
   rdFixed(1)=uicontrol('Style',strRadio,'Unit','normalized','Position',[0.055,0.63,0.19,0.08],...
      'String','Fixed Step','BackgroundColor',0.8*[1,1,1],'Value',0,'CallBack','sys_analysis(1,1);');
   rdFixed(2)=uicontrol('Style',strRadio,'Unit','normalized','Position',[0.285,0.63,0.25,0.08],...
      'String','Variable Step','BackgroundColor',0.8*[1,1,1],'Value',1,'CallBack','sys_analysis(1,2);');
   rdFixed(3)=uicontrol('Style',strCheck,'Unit','normalized','Position',[0.615,0.73,0.28,0.08],...
      'String','Show Linearised','BackgroundColor',0.8*[1,1,1],'Value',1);
   extra_funs(10,[0.61,0.054],[0.975,0.68]);
   uicontrol('Style','Text','Unit','normalized','Position',[0.635 0.62 0.19 0.08],...
      'BackgroundColor',0.8*[1,1,1],'HorizontalAlignment','left','String',' Input Signal');
   hh(1)=uicontrol('Style','PushButton','Unit','normalized','Position',[0.625,0.52,0.1,0.1],...
      'String','Sine','BackgroundColor',0.8*[1,1,1],'CallBack','sys_analysis(1,5);');
   hh(2)=uicontrol('Style','PushButton','Unit','normalized','Position',[0.74,0.52,0.1,0.1],...
      'String','Sqr','BackgroundColor',0.8*[1,1,1],'CallBack','sys_analysis(1,6);');
   hh(3)=uicontrol('Style','PushButton','Unit','normalized','Position',[0.855,0.52,0.1,0.1],...
      'String','Saw','BackgroundColor',0.8*[1,1,1],'CallBack','sys_analysis(1,7)');
   display_str(0.65,0.44,'Peak Value',[0,0,0],'on',9);
   h(10)=uicontrol('Style','Edit','Unit','normalized','Position',[0.84 0.40 0.1 0.07],...
      'String',num2str(ww(7)),'BackgroundColor',[1,1,1]);
   display_str(0.65,0.34,'Peak Range',[0,0,0],'on',9);
   h(11)=uicontrol('Style','Edit','Unit','normalized','Position',[0.84 0.30 0.10 0.07],...
      'String',num2str(ww(8)),'BackgroundColor',[1,1,1]);
   display_str(0.65,0.24,'Frequency',[0,0,0],'on',9);
   h(12)=uicontrol('Style','Edit','Unit','normalized','Position',[0.84 0.20 0.1 0.07],...
      'String',num2str(ww(9)),'BackgroundColor',[1,1,1]);
   display_str(0.65,0.14,'Freq Range',[0,0,0],'on',9);
   h(13)=uicontrol('Style','Edit','Unit','normalized','Position',[0.84 0.10 0.1 0.07],...
      'String',num2str(ww(1)),'BackgroundColor',[1,1,1]);

%%  Execute and Close Window
   uicontrol('Style','PushButton','Unit','normalized','Position',[0.65,0.85,0.14,0.09],...
      'String','Change','BackgroundColor',0.8*[1,1,1],'CallBack','sys_analysis(1,10);');
   uicontrol('Style','PushButton','Unit','normalized','Position',[0.81,0.85,0.14,0.09],...
      'String','Default','BackgroundColor',0.8*[1,1,1],'CallBack','delete(gcf)');
   set(gcf,'UserData',{h,hh,rdFixed,hSimAlg});
   sys_analysis(1,1);
else
   figure(g_SimPars);   
end

%--------------------------------------------------------------------------------
%bode_asymp gets the information for Bode asymptotics
%--------------------------------------------------------------------------------
function [wpos,ypos]=bode_asymp(G,w)
[zer,pol,gain]=zpkdata(G,'v'); wpos=[]; pos1=[];
for i=1:length(zer);
   if isreal(zer(i)), wpos=[wpos, abs(zer(i))]; pos1=[pos1,20];
   else, if imag(zer(i))>0, wpos=[wpos, abs(zer(i))]; pos1=[pos1,40]; end
end, end
for i=1:length(pol);
   if isreal(pol(i)), wpos=[wpos, abs(pol(i))]; pos1=[pos1,-20];
   else, if imag(pol(i))>0, wpos=[wpos, abs(pol(i))]; pos1=[pos1,-40]; end
end, end
wpos=[wpos w(2) w(end)]; pos1=[pos1,0,0];
[wpos,ii]=sort(wpos); pos1=pos1(ii); ii=find(abs(wpos)<eps); kslp=0; w_start=1000*eps; 
if length(ii)>0, 
   kslp=sum(pos1(ii)); ii=(ii(end)+1):length(wpos); wpos=wpos(ii); pos1=pos1(ii); 
end
while 1
   [ypos1,pp]=bode(G,w_start);
   if isinf(ypos1), w_start=w_start*10; else, break; end
end
wpos=[w_start wpos]; ypos(1)=20*log10(ypos1); pos1=[kslp pos1];
for i=2:length(wpos)
   kslp=sum(pos1(1:i-1)); ypos(i)=ypos(i-1)+kslp*log10(wpos(i)/wpos(i-1));      
end
ii=find(wpos>=w(2) & wpos<=w(end)), wpos=wpos(ii); ypos=ypos(ii);

%--------------------------------------------------------------------------------
%rloc_asymp gets the information for root locus asymptotics
%--------------------------------------------------------------------------------
function [xpos,ypos]=rloc_asymp(G)
[zer,pol,gain]=zpkdata(G,'v');
ii=find(abs(zer)<1e10); 
if length(ii)>0, zer=zer(ii); end
nExcess=length(pol)-length(zer); 
if nExcess>0
   pp=(sum(real(pol))-sum(real(zer)))/nExcess; deltaP=pi/nExcess; 
   xx=get(gca,'Xlim'); yy=get(gca,'YLim'); xx1=(xx(1)-pp)*tan(deltaP); 
end
xpos=[pp*ones(1,nExcess); zeros(1,nExcess)]; ypos=zeros(2,nExcess); 
for i=1:nExcess
   PAngle=(2*i-1)*deltaP; Kslp=tan(PAngle); 
   if (pi/2>=PAngle & PAngle>=0), xP=xx(2); yP=yy(2);
   elseif (pi>=PAngle&PAngle>pi/2), xP=xx(1); yP=yy(2);
   elseif (3*pi/2>=PAngle&PAngle>pi), xP=xx(1); yP=yy(1);
   else, xP=xx(2); yP=yy(1); end
   xx1=xP; yy1=Kslp*(xx1-pp); 
   if yy1>yy(2) | yy1<yy(1), yy1=yP; xx1=yy1/Kslp+pp; end
   xpos(2,i)=xx1; ypos(2,i)=yy1;
end
if nargout==0, line(xpos,ypos); end