📄 sys_analysis.m
字号:
%sys_analysis is the function used in CtrlLAB for administrating system analysis tasks.
%
%Available Analysis Tasks are:
%
% 1 for Bode diagrams
% 2 for Nyquist plots
% 3 for Nichols charts
% 4 for inverse Nyquist plots
% 5 for root loci
% 6 for step responses
% 7 for impulse responses
% 8 for stability margins
% 9 for H2 and Hinf norms evaluation
% 10 for time moments
% 11 for Markov parameters
% 12 for analytical solutions
%
%Available lists of functions under this module
%
% make_model -- prepares models to be analysed.
% plt_range_dat -- updates the plot range parameters
% refine_plot -- refines the plots according to the preset format
% dly_freq -- evaluates the frequency response of systems with delay
% time_resp -- evaluates the time domain response of the given system
% plt_common -- extracts common setting in system analysis
% make_model -- prepares models for further analysis tasks
% tim_moments -- evaluates the time moments of the system
% markov_paras -- evaluates the Markov parameters of the system
% asol_disp -- displays the analytical solutions of the LTI systems
% p_frac -- performs trianglar manipulations on the residue() function
% sim_pars --
% bode_asymp --
% rloc_asymp --
%
%Designed by Professor Dingyu Xue
%School of Information Science and Engineering, Northeastern University
%Shenyang 110006, P R China
%Email: xue_dy@hotmail.com
%
%This module is for CtrllAB 3.0, (c) 1996-1999
%Last modified 5 October, 1999
%---------------------------------------------------------------------------------
function sys_analysis(nTask,arg1,arg2)
switch nargin
case 0
%call simulation parameters dialog box
sim_pars;
case 1
%bring the main window to front
figure(findobj('Tag','CtrlLABMain'));
uu=get(gcf,'UserData'); g1=get(uu{1}(1),'UserData');
if length(g1)==0,
warndlg('Plant model dose not exist, please enter it first!','Matrix Processor Warning:');
else
uu{1}(7)=nTask; set(gcf,'UserData',uu);
g_comp=extra_funs(5,4,'Checked',26:27,[]);
if any([1:7,14:16]==nTask), [x1,x2,n_points,k_extra]=plt_range_dat(nTask); end
switch nTask
case {1,2,3,4,14,15}
if length(x1)==0, [x,y,w,arg4]=dly_freq(k_extra(2),g_comp(1),nTask);
else,
w=logspace(log10(x1),log10(x2),n_points);
[x,y,w0,arg4]=dly_freq(k_extra(2),g_comp(1),nTask,w);
end
if length(g_comp)>1, [x1,y1,w0,arg5]=dly_freq(k_extra(2),g_comp(2),nTask,w); end
switch nTask
case {1,14}, graf_tool(1);
if length(g_comp)>1,
subplot(211), semilogx(w,20*log10(x),w,20*log10(x1));
if nTask==14, line(arg5(1,:),arg5(2,:)); end
ylabel('Magnitude (dB)'), title('Bode Diagram')
subplot(212), semilogx(w,y,w,y1); ylabel('Phase (degrees)')
else
subplot(211), semilogx(w,20*log10(x));
if nTask==14, line(arg4(1,:),arg4(2,:)); end
ylabel('Magnitude (dB)'), title('Bode Diagram')
subplot(212), semilogx(w,y); ylabel('Phase (degrees)')
end
title('Frequenct \omega (Radian)')
case {2,4,15},
if nTask<=7, graf_tool(nTask); else, graf_tool(2); end
p_tmp=pi/180*y; H=x.*(cos(p_tmp)+sqrt(-1)*sin(p_tmp)); H=[H,conj(H)];
if length(g_comp)>1,
p_tmp=pi/180*y1; H1=x1.*(cos(p_tmp)+sqrt(-1)*sin(p_tmp)); H1=[H1,conj(H1)];
if any([2,15]==nTask), h=plot(real(H),imag(H),real(H1),imag(H1));
else, h=plot(real(1./H),imag(1./H),real(1./H1),imag(1./H1)); end
else
if any([2,15]==nTask), h=plot(real(H),imag(H));
else, h=plot(real(1./H),imag(1./H)); end
end
set(h,'UserData',w); xlabel('Real Axis'); ylabel('Imaginary Axis'),
switch nTask
case 2, title('Nyquist Plot')
case 4, title('Inverse Nyquist Plot');
case 15, title('ATAN Nyquist Plot')
end
case 3, graf_tool(3);
bHold=ishold;
if ~bHold, delete(extra_funs(2)); ngrid('new'); end
if length(g_comp)>1, h=plot(y,20*log10(x),y1,20*log10(x1));
else, h=plot(y,20*log10(x)); end
set(h,'UserData',w); xlabel('Phase (degree)')
ylabel('Magnitude (dB)'), title('Nichols Chart')
if ~bHold, hold off; end
end
case {5,16},
%for root locus analysis, only the open-loop model is meaningful.
[G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_reltn]=plt_common(1,g_comp(1));
graf_tool(5);
if length(x1)==0, [r,k]=rlocus(G_forward);
else,
if ~isfinite(x2), x2=1000; end
if abs(x1)<eps, x1=0.001; end
k=logspace(log10(x1),log10(x2),n_points); r=rlocus(G_forward,k);
end
if length(g_comp)>1,
r1=rlocus(G_Sys,k); h=plot(real(r'),imag(r'),real(r1'),imag(r1'));
[z,p,k]=zpkdata(G_Sys,'v');
if length(z)>0, h=line(real(z),imag(z),'Marker','o','LineStyle','none'); end
if length(p)>0, h=line(real(p),imag(p),'Marker','x','LineStyle','none'); end
else, h=plot(real(r'),imag(r')); end
[z,p,k]=zpkdata(G_forward,'v');
if length(z)>0, h=line(real(z),imag(z),'Marker','o','LineStyle','none'); end
if length(p)>0, h=line(real(p),imag(p),'Marker','x','LineStyle','none'); end
set(h,'UserData',k); xlabel('Real Axis'); ylabel('Imaginary Axis'), title('Root Locus')
if nTask==16, rloc_asymp(G_forward); if length(g_comp)>1, rloc_asymp(G_Sys); end, end
case {6,7},
if length(x1)==0
[y,t]=time_resp(k_extra,g_comp(1),nTask);
if length(y)==0; return; end
else
t=x1:(x2-x1)/n_points:x2; y=time_resp(k_extra,g_comp(1),nTask,t);
end
if length(g_comp)>1, y1=time_resp(k_extra,g_comp(2),nTask,t); end
graf_tool(nTask);
if length(g_comp)>1, plot(t,y,t,y1); else, plot(t,y); end
xlabel('Time (Sec)'); ylabel('Response Magnitude');
if nTask==6, title('Step Response');
else, title('Impulse Response'); end
case 8, %gain phase margins
figure(findobj('Tag','CtrlLABMain')); uu0=get(gcf,'UserData');
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 Td>0, [nP,dP]=pade(Td,nPade); GP=tf(nP,dP); else, GP=1; end
G1=G_Sys*H_Sys*GP; [Gm,Pm,Wcg,Wcp]=margin(G1);
if length(g_comp)>1, G2=G1*Gc_Sys; [Gm1,Pm1,Wcg1,Wcp1]=margin(G2); end
display_str;
display_str(0.1,0.8,['Gain Margin: ',display_str(Gm),' at \omega=' display_str(Wcg)],[1,0,0]);
display_str(0.1,0.6,['Phase Margin: ',display_str(Pm) '^o at \omega=' display_str(Wcp)],[1,0,0]);
if length(g_comp)>1,
display_str(0.1,0.4,['Gain Margin: ',display_str(Gm1),' at \omega=' display_str(Wcg1)],[0,0,1]);
display_str(0.1,0.2,['Phase Margin: ',display_str(Pm1) '^o at \omega=' display_str(Wcp1)],[0,0,1]);
end
case {9,10,11,12}
[G_Sys1,g1]=make_model(g_comp(1));
if length(g_comp)>1, [G_Sys2,g2]=make_model(g_comp(2)); end
switch nTask
case 9, %H_2 & H_inf norm evaluation
H2=norm(g1,2); Hinf=norm(g1,inf); display_str;
display_str(0.1,0.8,['H_2 Norm=',display_str(H2),', H_\infty Norm=' display_str(Hinf)],[1,0,0]);
if length(g_comp)>1,
H2=norm(g2,2); Hinf=norm(g2,inf);
display_str(0.1,0.5,['H_2 Norm=',display_str(H2),', H_\infty Norm=' display_str(Hinf)],[0,0,1]);
end
case 10, % time moments
M1=tim_moments(g1,7); display_str;
display_str(0.1,0.8,'Time Moments:',[1,0,0]); display_str(0.1,0.65,M1,[1,0,0]);
if length(g_comp)>1,
M2=tim_moments(g2,7);
display_str(0.1,0.4,'Time Moments:',[0,0,1]); display_str(0.1,0.25,M2,[0,0,1]);
end
case 11, %Markov parameters
M1=markov_paras(g1,7); display_str;
display_str(0.1,0.8,'Markov Parameters:',[1,0,0]); display_str(0.1,0.65,M1,[1,0,0]);
if length(g_comp)>1,
M2=markov_paras(g2,7);
display_str(0.1,0.4,'Markov Parameters:',[0,0,1]); display_str(0.1,0.25,M2,[0,0,1]);
end
case 12, %analytic solutions
key=get(uu{4}(7),'Checked'); display_str;
xL(1)=asol_disp(g1,key,0.1,0.8,[1,0,0]);
if length(g_comp)>1, xL(2)=asol_disp(g2,key,0.1,0.3,[0,0,1]); end
end
case 13
[G_Sys1,g1]=make_model(g_comp(1),[]);
if length(g_comp)>1, [G_Sys2,g2]=make_model(g_comp(2),[]); end
inv_T=1/norm(g1,inf); t1=1-inv_T; t2=1+inv_T; t3=-2*asin(0.5*inv_T)*180/pi; display_str;
display_str(0.1,0.8,['Guaranteed Gain Margin: [' num2str(t1) ',' num2str(t2) ']'],[1,0,0]);
display_str(0.1,0.6,['Guaranteed Phase Margin: [' num2str(t3) '^o,' num2str(-t3) '^o]'],[1,0,0]);
if length(g_comp)>1,
inv_T=1/norm(g2,inf); t1=1-inv_T; t2=1+inv_T; t3=-2*asin(0.5*inv_T)*180/pi;
display_str(0.1,0.4,['Guaranteed Gain Margin: [' num2str(t1) ',' num2str(t2) ']'],[0,0,1]);
display_str(0.1,0.2,['Guaranteed Phase Margin: [' num2str(t3) '^o,' num2str(-t3) '^o]'],[0,0,1]);
end
end
if nTask<=7, refine_plot(nTask); end
end
case 2
uu=get(gcf,'UserData'); g_graf=gcf;
switch arg1
case 1
kk=get(uu{3}(1),'Value'); set(uu{3}(2),'Value',~kk);
set(uu{1}(1:6),'Visible',extra_funs(6,~kk)); set(uu{1}(7:8),'Visible',extra_funs(6,kk));
case 2
kk=get(uu{3}(2),'Value'); set(uu{3}(1),'Value',~kk);
set(uu{1}(1:6),'Visible',extra_funs(6,kk)); set(uu{1}(7:8),'Visible',extra_funs(6,~kk));
case 4
keyAlg=get(uu{4},'Value'); ww(6)=num2str(get(uu{1}(9),'String'));
bFixedStep=get(uu{3}(1),'Value');
if bFixedStep==0,
ww(2)=get(uu{1}(2),'String'); eval(['tmin=',ww(2),';']);
ww(3)=get(uu{1}(4),'String'); eval(['tmax=',ww(3),';']);
ww(4)=get(uu{1}(6),'String'); eval(['terr=',ww(4),';']);
else,
ww(5)=get(uu{1}(8),'String'); eval(['tmin=',ww(5),';']);
end
ww(7)=get(uu{1}(10),'String'); ww(8)=get(uu{1}(11),'String');
ww(9)=get(uu{1}(12),'String'); ww(1)=get(uu{1}(13),'String');
uu{6}=ww; set(g_graf,'UserData',uu);
end
case 3
[mag,pha,w]=bode(arg1);
subplot(211), semilogx(w,20*log10(mag(:)));
ylabel('Magnitude (dB)'), title('Bode Diagram')
subplot(212), semilogx(w,pha(:)); ylabel('Phase (degrees)')
title('Frequenct \omega (Radian)')
end
%-------------------------------------------------------
%plt_range_dat is used to set the plot range parameters.
%-------------------------------------------------------
function [x1,x2,n_points,k_extra]=plt_range_dat(nTask)
%output argument initialization
key=0;
x1=[]; x2=[]; n_points=[]; k_extra=zeros(1,2); k_match=0;
figure(findobj('Tag','CtrlLABMain')); uu=get(gcf,'UserData');
if length(uu)>=10,
u_data=uu{10}; if length(u_data)>0, k_match=u_data(5); end
end
switch k_match
case {1,2,3,4}
%set frequency response data
if any([1:4]==nTask), key=1; end
case 5,
%set root locus data
if nTask==5, key=1; end
case {6,7}
%set time-domain response data
if any([6:7]==nTask), key=1; end
end
if key==1
x1=u_data(1); x2=u_data(2); n_points=u_data(3); k_extra=u_data([4,6]);
end
%-----------------------------------------------------------------------
%refine_plot is used to refine the plots according to the preset format.
%-----------------------------------------------------------------------
function refine_plot(nTask)
it=get(gca,'Title'); pp=get(it,'Position'); xL=get(gca,'XLim'); yL=get(gca,'YLim');
set(it,'Position',[pp(1), 1.006*yL(2)-0.006*yL(1),pp(3)])
if nTask==1
ii=extra_funs(2); axes(ii(2)); it=get(gca,'Title'); pp=get(it,'Position');
yL=get(gca,'YLim'); set(it,'Position',[pp(1), 1.006*yL(2)-0.006*yL(1),pp(3)])
else,
ix=get(gca,'xLabel'); x0=0.85*xL(2)+0.15*xL(1);
y0=0.08*yL(2)+0.92*yL(1); set(ix,'Position',[x0,y0,0]);
end
%update plot properties according to the settings
graf_tool(2,1); graf_tool(2,2); graf_tool(2,4);
%-------------------------------------------------------------------------------------
%dly_freq is a function used to evaluate the frequency response of a system containing
%a pure time delay element exp(-tdly*s) over the frequency range of W.
%
% [mag,pha,w,arg4]=dly_freq(k_Pade,nCtrl,nTask,w)
%
%where
% mag,pha -- the magnitude and phase over w
% w -- the frequency vector
% nCtrl -- whether the compensated system are required
% nTask -- the system analysis task code
% k_Pade -- whether Pade approximation is used for delay
%-------------------------------------------------------------------------------------
function [mag,pha,w,arg4]=dly_freq(k_Pade,nCtrl,nTask,w)
[G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_reltn]=plt_common(k_Pade,nCtrl);
arg4=[];
if nargin==3
%selecting the frequency vector w in an automatic way
[num,den]=tfdata(G_Sys,'v'); w=freqint2(num,den,30);
end
if g_loop==1,
%processing open-loop
if k_Pade==1 & Td>0, G_forward=G_forward*GP;
elseif Td>0, G_forward.Td=Td; end
G_overall=G_forward; [mag,pha]=bode(G_forward,w); mag=mag(:); pha=pha(:);
else
%processing closed-loop
if Td>0 & k_Pade==0
Gf=G_forward; Gf.Td=Td;
[mag,pha]=bode(Gf,w); mag=mag(:); pha=pi*pha(:)/180;
Hf=mag.*exp(-sqrt(-1)*pha);
if isa(H_Sys,'double'), H=H_Sys; else, H=nyquist(H_Sys,w); end
switch g_reltn
case 1, HH=Hf./(1+H.*Hf);
case {2,3}, HH=ones(size(Hf))./(1+H.*Hf); if g_reltn==3, HH=1-HH; end
end
mag=abs(HH); pha=phase(HH')'*180/pi;
else
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -