第第6 6章章 仿真软件仿真软件 —MATLAB —MATLAB基础与应用基础与应用 6.4 线性系统分析与设计 MATLAB的控制工具箱是的控制工具箱是MATLAB最早最早的工具箱之一，也是控制系统的计算机辅助设计的工具箱之一，也是控制系统的计算机辅助设计中最为流行的设计工具控制工具箱适用于线性中最为流行的设计工具控制工具箱适用于线性时不变系统时不变系统(LTI)，可实现线性系统时域或频域，可实现线性系统时域或频域的分析、设计和建模可处理连续系统，也可处的分析、设计和建模可处理连续系统，也可处理离散系统；可使用经典或现代的技术理离散系统；可使用经典或现代的技术线性系统的描述线性系统的描述MATLAB只处理只处理矩阵矩阵这一种数学形式，各种这一种数学形式，各种控制系统的描述必须使用矩阵来表达控制系统的描述必须使用矩阵来表达 一、一、连续系统的数学描述连续系统的数学描述1 1、、状态空间描述法状态空间描述法在在MATLAB中，这个系统写为中，这个系统写为A、、B、、C、、D四个矩阵的形式即可，当然矩阵维数要匹配四个矩阵的形式即可，当然矩阵维数要匹配。
也可用也可用SYS = SS(A,B,C,D) 建立建立ss模型，模型，SYS = SS(A,B,C,D,Ts) 建立离散建立离散ss模型 %控制系统模型的描述方式控制系统模型的描述方式a=[1 2;3 4];b=[0;1];c=[1 1];d=1;f=ss(a,b,c,d) a = x1 x2 x1 1 2 x2 3 4c = x1 x2 y1 1 1Continuous-time model.b = u1 x1 0 x2 1d = u1 y1 1 f1=ss(a,b,c,d,0.1)a = x1 x2 x1 1 2 x2 3 4c = x1 x2 y1 1 1Discrete-time model. b = u1 x1 0 x2 1d = u1 y1 1 2 2、传递函数描述法、传递函数描述法 传递函数使用分子、分母的多项式表示，即传递函数使用分子、分母的多项式表示，即传递函数使用分子、分母的多项式表示，即传递函数使用分子、分母的多项式表示，即numnum和和和和denden两个向量。
两个向量两个向量两个向量同样可用同样可用同样可用同样可用SYS = TF(NUM,DEN)SYS = TF(NUM,DEN)建立建立建立建立tftftftf模型 num=[1 2 3];den=[2 2 3 4];yy=tf(num,den)Transfer function: s^2 + 2 s + 3-----------------------2 s^3 + 2 s^2 + 3 s + 4 3 3、零极点描述法、零极点描述法在在MATLAB中，这种形式使用增益中，这种形式使用增益k、分子、分子零点向量零点向量z、分母极点向量、分母极点向量p表示注意：根据注意：根据MATLAB的约定，多项式的根的约定，多项式的根（零极点）存在列向量中，行向量中存多项式（零极点）存在列向量中，行向量中存多项式的系数这里，除了的系数这里，除了系数系数k使用行向量使用行向量外，外，z和和p使用列向量使用列向量同样可用同样可用 SYS = ZPK(Z,P,K)建立建立zpk模型 Zero/pole/gain: 2 (s-1) (s-2)-----------------(s-3) (s-5) (s-7)z=[1;2];p=[3;5;7];k=2;sys=zpk(z,p,k) 4 4、部分分式描述法、部分分式描述法 在传递函数没有相同极点时与部分分式在传递函数没有相同极点时与部分分式相互转换：相互转换： [r,p,k]=residue(num,den) [num,den]=residue(r,p,k)二、离散系统的描述二、离散系统的描述在在MATLAB中，离散系统表达方式与连续系统中，离散系统表达方式与连续系统相同，也可以为状态空间、多项式传递函数、相同，也可以为状态空间、多项式传递函数、零极点增益以及部分分式形式。
零极点增益以及部分分式形式 三、闭环系统的表达三、闭环系统的表达 以上已经给出开环系统的模型表达有时需要系以上已经给出开环系统的模型表达有时需要系以上已经给出开环系统的模型表达有时需要系以上已经给出开环系统的模型表达有时需要系统的闭环模型，统的闭环模型，统的闭环模型，统的闭环模型，MATLABMATLAB提供了一组这样的函数提供了一组这样的函数提供了一组这样的函数提供了一组这样的函数: :• •feedbackfeedback 反馈连接反馈连接反馈连接反馈连接 SYS = feedback(SYS1,SYS2,sign) SYS = feedback(SYS1,SYS2,sign) [A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign[A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign) )[num,den]=feedback(num1,den1,num2,den2,sign) [num,den]=feedback(num1,den1,num2,den2,sign) sign=+(-)1sign=+(-)1反馈极性反馈极性反馈极性反馈极性, ,缺省缺省缺省缺省-1-1为负反馈为负反馈为负反馈为负反馈• •parallelparallel系统并联系统并联系统并联系统并联• •seriesseries系统串联系统串联系统串联系统串联 FEEDBACK Feedback connection of two LTI models. SYS = FEEDBACK(SYS1,SYS2) computes an LTI model SYS for the closed-loop feedback system u --->O---->[ SYS1 ]----+---> y | | y = SYS * u +-----[ SYS2 ]<---+ Negative feedback is assumed and the resulting system SYS maps u to y. To apply positive feedback, use the syntax SYS = FEEDBACK(SYS1,SYS2,+1). num1=[1 1];den1=[1 5 6];sys1=tf(num1,den1);sys2=tf(1,1);sysb=feedback(sys1,sys2)[numb,denb]=feedback(num1,den1,1,1)Transfer function: s + 1-------------s^2 + 6 s + 7 numb = 0 1 1 denb = 1 6 7 6.4.2 6.4.2 模型之间的转换模型之间的转换一、线性系统模型之间的转换一、线性系统模型之间的转换ssss—状态空间、状态空间、tftf—传递函数、传递函数、zpzp—零极点：零极点：[num,den]=ss2tf(a,b,c,d,iu) 状态空间到传状态空间到传函函 [z,p,k]=ss2zp(a,b,c,d,iu) 状态空间到零极状态空间到零极 [a,b,c,d]=tf2ss(num,den) 传函到状态空传函到状态空间间 [z,p,k]=tf2zp(num,den) 传函到零极传函到零极 [a,b,c,d]=zp2ss(z,p,k) 零极到状态空零极到状态空间间 [num,den]=zp2tf(z,p,k) 零极到传函零极到传函 [r,p,k]=residue(num,den) 传函到部分分传函到部分分式式 [num,den]=residue(r,p,k) 部分分式到传函部分分式到传函 用法举例：用法举例：已知系统状态空间模型为：已知系统状态空间模型为：%转换为传递函数模型：转换为传递函数模型：A=[0 1; -1 -2]; B=[0;1]; C=[1,3]; D=[1];[num,den]=ss2tf(A,B,C,D,1)%iu用来指定第用来指定第n个输入，当只有一个输入时可忽略。
个输入，当只有一个输入时可忽略》》num=1 5 2; den=1 2 1;%转换为零极点增益模型：转换为零极点增益模型：[z,p,k]=ss2zp(A,B,C,D)》》z= -4.5616 p= -1 k=1 -0.4384 -1tf(num,den) s^2 + 5 s + 2-------------s^2 + 2 s + 1zpk(z,p,k) (s+0.4384) (s+4.562)-------------------- (s+1)^2 已知一个单输入三输出系统的传递函数模型为：已知一个单输入三输出系统的传递函数模型为：%转换为状态空间模型转换为状态空间模型num=[0 0 -2;0 -1 -5;1 2 0];den=[1 6 11 6];[A,B,C,D]=tf2ss(num,den)》》A= -6 -11 -6 B= 1 C= 0 0 -2 D= 0 1 0 0 0 0 -1 -5 0 0 1 0 0 1 2 0 0 1 1、连续系统的离散化、连续系统的离散化格式：格式：[Ad,Bd]=c2d(A,B,ts) [Ad,Bd,Cd,Dd]=c2dm(A,B,C,D,ts,’method’)， [numz,denz]=c2dm(num,den,ts,’method’) 1)c2d 1)c2d使用零阶保持器离散化，只有状态空间形式使用零阶保持器离散化，只有状态空间形式 2)c2dm2)c2dm既有状态空间形式，又有传递函数形式；既有状态空间形式，又有传递函数形式； 3) 3)参数参数tsts是采样周期是采样周期T T；； 4)method4)method指指定定转转换换方方式式，，其其中中“zoh”“zoh”表表示示采采用用零零阶阶保保持持器器；；“foh”“foh”表表示示采采用用三三角角形形近近似似；；“tustin”“tustin”表表示示采采用用双双线线性性变变换换；；“prewarp”“prewarp”表表示示采采用用指指定定转转折折频频率率的的双双线性变换线性变换; ; 系统默认为零阶保持器法。
系统默认为零阶保持器法二、连续系统与离散系统之间的转换二、连续系统与离散系统之间的转换 [Ad,Bd]=c2d(A,B,0.1)[Ad1,Bd1,Cd1,Dd1]=c2dm(A,B,C,D,0.1,'tustin')[numd,dend]=c2dm(num,den,0.1,'tustin')A=[0 1; -1 -2]; B=[0;1]; C=[1,3]; D=[1];[num,den]=ss2tf(A,B,C,D)s^2 + 5 s + 2-------------s^2 + 2 s + 1BdA = 0 1 B = 0 -1 -2 1 C = 1 3 D = 1-------------------------------------- c2d 连续连续ss到离散到离散ss d2c 离散离散ss到连续到连续ss c2dm 连续到离散连续到离散(可用不同方法可用不同方法)d2cm离散到连续离散到连续 (可用不同方法可用不同方法) 6.4.3 6.4.3 时间响应分析时间响应分析一、单位阶跃响应一、单位阶跃响应1.1.对对状态空间状态空间描述的系统描述的系统y=step(A,B,C,D,iu,T) 求求ss系统对系统对iu个输入的在个输入的在T时间向量下的单时间向量下的单位阶跃响应，位阶跃响应，y是列向量。
是列向量[y,x]=step(A,B,C,D,iu,T) 同上同上,并返回状态并返回状态x的变化的变化step(A,B,C,D,iu,T) 直接画出单位阶跃响应曲线直接画出单位阶跃响应曲线 2.2.对传递函数描述的系统对传递函数描述的系统y=step(num,den,T) 求求tf系统在系统在T时间向量下时间向量下的单位阶跃响应，的单位阶跃响应，y是列向量是列向量step(num,den,T) 直接画出单位阶跃响应曲线直接画出单位阶跃响应曲线3.3.对各种形式描述的系统对各种形式描述的系统step(sys1,’s1’,sys2,’s2’,…) sys(i)可是各种形式描述的可是各种形式描述的LTI系统，时间系统，时间是自动给出的，也可用时间向量和给出最终时是自动给出的，也可用时间向量和给出最终时间两种方法确定，同样可直接画出响应曲线间两种方法确定，同样可直接画出响应曲线 %单位阶跃响应t=0:0.1:100;num1=[0 2];den1=[10 1];wn=1.5;z=0.2;A=[0 1;-wn^2 -2*z*wn];B=[0;wn^2];C=[1 0];D=0;z=[-1 -2]';p=[-3;-4;-5];k=10; %z,p is column vectorsys1=tf(num1,den1)sys2=ss(A,B,C,D)sys3=zpk(z,p,k)subplot(3,1,1)step(sys1,'r',sys2,'b',sys3,'c')% time range and number of points are chosen automatically.subplot(3,1,2)step(sys1,t,‘b') %uses the user-supplied time vector T subplot(3,1,3)step(sys2,40,'b') % the step response from t=0 to the final time Transfer function: 2--------10 s + 1Zero/pole/gain: 10 (s+1) (s+2)-----------------(s+3) (s+4) (s+5)a = x1 x2 x1 0 1b = u1 x1 0c = x1 x2 y1 1 0d = u1 y1 0 Continuous-time model. %阻尼系数对典型二阶系统的阶跃响应的影响a=[0 0.3 0.7 1 3];b=0.2*ones(1,5);for i=1:5 num=b(i)^2; den=[1 2*a(i)*b(i) b(i).^2]; y(:,i)=step(num,den,t);Endplot(t,y(:,1)','r',t,y(:,2)','b',t,y(:,3)','m',t,y(:,4)','g',t,y(:,5)','k')legend('a=0','a=0.3','a=0.7','a=1','a=3') t=0:0.1:40;sys1=zpk([-1 -0.5],[-0.1 -1.5 -2],2);%0型系统开环传函sys2=zpk([-1 -0.5],[0 -0.1 -3 -4],2);%1型系统开环传函sysf=tf(1,1);%反馈环节sys11=feedback(sys1,sysf);%0型系统闭环环传函sys22=feedback(sys2,sysf);%1型系统闭环环传函subplot(2,1,1)step(sys11,'r',sys22,'b')%对单位阶跃，0型有误差，1型无误差subplot(2,1,2)syss=tf(1,[1 0]);%1/ssys111=series(sys11,syss);%没有直接的斜坡响应函数， % y=G/s^2=1/s*(G/s)=step(G/s)sys222=series(sys22,syss);y1=step(sys111,t);y2=step(sys222,t);h=plot(t,y1,'r',t,y2,'g',t,t,'k');%对单位斜坡，0型误差无穷，1型有误差 sys1 2 (s+1) (s+0.5) --------------------- (s+0.1) (s+1.5) (s+2) sys2 2 (s+1) (s+0.5)---------------------s (s+0.1) (s+1.5) (s+2) 二、脉冲响应1.1.对状态空间描述的系统对状态空间描述的系统y=impulse(A,B,C,D,iu,T) 求求ss系统对系统对iu个输个输入的在入的在T时间向量下的单位阶跃响应，时间向量下的单位阶跃响应，y是列向是列向量。
量[y,x]=impulse(A,B,C,D,iu,T) 同上同上,并返回状并返回状态态x的变化的变化2.对传递函数描述的系统对传递函数描述的系统y=impulse (num,den,T) 求求tf系统在系统在T时间向时间向量下的单位阶跃响应，量下的单位阶跃响应，y是列向量是列向量 3.3.对各种形式描述的系统对各种形式描述的系统直接画出响应曲线直接画出响应曲线impulse(sys1,’s1’,sys2,’s2’,…) sys(i)可以是各种形式描述的可以是各种形式描述的LTI系统，时间系统，时间是自动给出的，也可以用时间向量和给出最终时是自动给出的，也可以用时间向量和给出最终时间两种方法确定间两种方法确定4.对离散系统对离散系统 dimpulse--------------------------------------t=0:0.1:5;dimpulse(numd,dend,t); t=0:0.1:10; wn=1.5; z=0.2;num1=[0 2]; den1=[10 1];A=[0 1; -wn^2 -2*z*wn];B=[0;wn^2]; C=[1 0]; D=0;z=[-1 -2]'; p=[-3;-4;-5]; k=10;sys1 = tf(num1,den1)sys2 = ss(A,B,C,D)sys3 = zpk(z,p,k)subplot(2,1,1)y1=impulse(num1,den1,t);y2=impulse(A,B,C,D,1,t);y3=impulse(sys3,t);plot(t,y1,'r',t,y2,'b',t,y3,'k')subplot(2,1,2)impulse(sys1,‘r’,sys2,‘b’,sys3,‘m',5)Transfer function: 2--------10 s + 1a = x1 x2 x1 0 1 b = u1 x1 0 c = x1 x2 y1 1 0 d = u1 y1 0 Zero/pole/gain: 10 (s+1) (s+2)-----------------(s+3) (s+4) (s+5) 三、对任意输入的响应三、对任意输入的响应1.1.1.1.对状态空间描述的系统对状态空间描述的系统对状态空间描述的系统对状态空间描述的系统y=lsim(A,B,C,D,U,T)y=lsim(A,B,C,D,U,T) 求求求求ssss系统对输入序列系统对输入序列系统对输入序列系统对输入序列U U在在在在T T时间向量下的响应。
时间向量下的响应时间向量下的响应时间向量下的响应[y,x]=lsim (A,B,C,D,U,T)[y,x]=lsim (A,B,C,D,U,T) 同上同上同上同上, ,并返回状态并返回状态并返回状态并返回状态x x的变化的变化的变化的变化2.2.对传递函数描述的系统对传递函数描述的系统对传递函数描述的系统对传递函数描述的系统y=lsim (num,den,U,T)y=lsim (num,den,U,T) 求求求求tf tf系统对任意输入系统对任意输入系统对任意输入系统对任意输入U U的响应的响应的响应的响应3.3.3.3.对各种形式描述的系统对各种形式描述的系统对各种形式描述的系统对各种形式描述的系统 还可直接画出响应曲线还可直接画出响应曲线•lsim(sys1,’s1’,sys2,’s2’,…,U,T) sys(i)可以是各种形式描述的可以是各种形式描述的LTI系统系统 a=[0 1 0;0 0 1;0 -12 -7];b=[0 0 1]';c=[2 3 1];d=[0];t=0:0.1:10;u=sin(t);y=lsim(a,b,c,d,u,t);plot(t,y,'r',t,u,'b--');gtext('- 输出y')gtext('--输入u')t=0:0.01:30;sys1=zpk([-1 -0.5],[-0.1 -1.5 -2],2);%0型系统开环传函sys2=zpk([-1 -0.5],[0 -0.1 -3 -4],2);%1型系统开环传函sysf=tf(1,1);%反馈环节sys1b=feedback(sys1,sysf);%0型系统闭环环传函sys2b=feedback(sys2,sysf);%1型系统闭环环传函lsim(sys1b,'r',sys2b,'b',t,t); 6.4.4 6.4.4 频率响应分析频率响应分析 频率响应是指系统对正弦输入信号的稳态响应，从 频率响应是指系统对正弦输入信号的稳态响应，从 频率响应是指系统对正弦输入信号的稳态响应，从 频率响应是指系统对正弦输入信号的稳态响应，从频率响应中可以得出带宽、增益、转折频率、闭环稳定频率响应中可以得出带宽、增益、转折频率、闭环稳定频率响应中可以得出带宽、增益、转折频率、闭环稳定频率响应中可以得出带宽、增益、转折频率、闭环稳定性等系统特征。
性等系统特征性等系统特征性等系统特征 频率特性是指系统在正弦信号作用下，稳态输出与 频率特性是指系统在正弦信号作用下，稳态输出与 频率特性是指系统在正弦信号作用下，稳态输出与 频率特性是指系统在正弦信号作用下，稳态输出与输入之比对频率的关系特性频率特性函数与传递函数输入之比对频率的关系特性频率特性函数与传递函数输入之比对频率的关系特性频率特性函数与传递函数输入之比对频率的关系特性频率特性函数与传递函数有直接的关系，记为：有直接的关系，记为：有直接的关系，记为：有直接的关系，记为： 一、一、BodeBode图图bode(num,den,w)bode(num,den,w) 直接绘出直接绘出直接绘出直接绘出bodebode图 包括了对数幅频特性图和对数相频特性图横包括了对数幅频特性图和对数相频特性图横坐标为频率坐标为频率w，采用对数分度，单位为弧度，采用对数分度，单位为弧度/秒；纵秒；纵坐标均匀分度，分别为幅值函数坐标均匀分度，分别为幅值函数20lgA(w)，以，以dB表表示；相角，以度表示示；相角，以度表示[m,p]= bode(num,den,w)[m,p]= bode(num,den,w) 求出幅值和相角，可以再用半对数坐标画图求出幅值和相角，可以再用半对数坐标画图求出幅值和相角，可以再用半对数坐标画图求出幅值和相角，可以再用半对数坐标画图subplot(2,1,1),semilogx(w,20*log10(m))subplot(2,1,2),semilogx(w,p) num=[1];den=[1 1 1];w=logspace(-1,2);figure(1)bode(num,den,w)gridfigure(2)[m,p]=bode(num,den,w);subplot(2,1,1)semilogx(w,20*log10(m))gridsubplot(2,1,2)semilogx(w,p)grid 二二.Nyquist.Nyquist曲线曲线 对于频率特性函数对于频率特性函数对于频率特性函数对于频率特性函数G(jw)G(jw)，给出，给出，给出，给出ww从负无穷到正无穷从负无穷到正无穷从负无穷到正无穷从负无穷到正无穷的一系列数值，分别求出的一系列数值，分别求出的一系列数值，分别求出的一系列数值，分别求出Im(G(jw))Im(G(jw))和和和和Re(G(jw))Re(G(jw))。
以以Re(G(jw)) Re(G(jw)) 为横坐标，为横坐标，为横坐标，为横坐标， Im(G(jw)) Im(G(jw)) 为纵坐标绘制成为极为纵坐标绘制成为极为纵坐标绘制成为极为纵坐标绘制成为极坐标频率特性图坐标频率特性图坐标频率特性图坐标频率特性图 Nyquist(num,den,{wmin,wmax}) 直接绘出直接绘出直接绘出直接绘出NyquistNyquist曲线，曲线，曲线，曲线， {wmin,wmax} {wmin,wmax}是频率范围，是频率范围，是频率范围，是频率范围，缺省自动给定缺省自动给定缺省自动给定缺省自动给定ww从负无穷到正无穷从负无穷到正无穷从负无穷到正无穷从负无穷到正无穷[real,imag]=[real,imag]= Nyquist(num,den,{wmin,wmax}) Nyquist(num,den,{wmin,wmax}) 得到实部、虚部，用得到实部、虚部，用得到实部、虚部，用得到实部、虚部，用plotplot绘制绘制w从零到正无穷变化部分从零到正无穷变化部分 另外，以上两种可用另外，以上两种可用另外，以上两种可用另外，以上两种可用ssss描述描述描述描述(A,B,C,D)(A,B,C,D)；或几种方法；或几种方法；或几种方法；或几种方法建立的建立的建立的建立的syssys所描述的模型来作图。
所描述的模型来作图所描述的模型来作图所描述的模型来作图 %典型环节wmin=0.01;wmax=100;sys1=tf([2],[1])%比例环节sys2=tf([2],[10 1])%惯性环节sys21=tf([2],[-10 1])%不稳定惯性环节sys3=tf([1],[1 0])%积分环节sys4=tf([1 0],[1])%微分环节wn=0.1;a=[0.1 0.3 1];for i=1:3sys5(i)=tf(wn^2,[1 2*a(i)*wn wn^2]) %不同阻尼比的二阶振荡环节end 2--------10 s + 1 2 - 2--------10 s - 1 1 --- s s ----------------------- ----------------------- -----------------------figure(1)%bode图subplot(2,2,1),bode(sys1)subplot(2,2,2),bode(sys2,'r',sys21,'b')subplot(2,2,3),bode(sys3,'r',sys4,'b')subplot(2,2,4)bode(sys5(1),'r',sys5(2),'g',sys5(3),'b') figure(2)%Nyquist图subplot(2,2,1)nyquist(sys1,{wmin,wmax})subplot(2,2,2)[r1,m1,w1]=nyquist([2],[10 1],{wmin,wmax});[r2,m2,w2]=nyquist([2],[-10 1],{wmin,wmax});plot(r1,m1,'r',r2,m2,'b')subplot(2,2,3)[r3,m3,w3]=nyquist([1],[1 0],{wmin,wmax});[r4,m4,w4]=nyquist([1 0],[1],{wmin,wmax});plot(r3,m3,'r',r4,m4,'b')axis([-2 2 -2 2])subplot(2,2,4)nyquist(sys5(1),'r',sys5(2),'g',sys5(3),'b',{wmin,wmax}) margin函数求函数求幅值和相位稳定裕度幅值和相位稳定裕度 marginmargin函数可以从频率响应数据中计算出幅值裕度、函数可以从频率响应数据中计算出幅值裕度、函数可以从频率响应数据中计算出幅值裕度、函数可以从频率响应数据中计算出幅值裕度、相角裕度以及对应的频率。
幅值裕度和相角裕度是针对相角裕度以及对应的频率幅值裕度和相角裕度是针对相角裕度以及对应的频率幅值裕度和相角裕度是针对相角裕度以及对应的频率幅值裕度和相角裕度是针对开环开环开环开环SISOSISO系统而言，它指示出系统闭环时的相对稳定系统而言，它指示出系统闭环时的相对稳定系统而言，它指示出系统闭环时的相对稳定系统而言，它指示出系统闭环时的相对稳定性margin(num,den)margin(num,den) 当不带输出变量引用时，当不带输出变量引用时，当不带输出变量引用时，当不带输出变量引用时，marginmargin可在当前图形窗可在当前图形窗可在当前图形窗可在当前图形窗口中绘制出带有裕量及相应频率显示的口中绘制出带有裕量及相应频率显示的口中绘制出带有裕量及相应频率显示的口中绘制出带有裕量及相应频率显示的BodeBode图，其中图，其中图，其中图，其中幅值裕度以分贝为单位幅值裕度以分贝为单位幅值裕度以分贝为单位幅值裕度以分贝为单位[Gm,Pm,Wcg,Wcp]=margin(num,den)[Gm,Pm,Wcg,Wcp]=margin(num,den) 计算出系统计算出系统计算出系统计算出系统wcpwcp，，，，GmGmGmGm幅值裕度，幅值裕度，幅值裕度，幅值裕度，PmPmPmPm相位裕度，相位裕度，相位裕度，相位裕度，及相应的相角交界频率及相应的相角交界频率及相应的相角交界频率及相应的相角交界频率WcgWcgWcgWcg、、、、截止频率截止频率截止频率截止频率WcpWcpWcpWcp 而不直接绘出而不直接绘出而不直接绘出而不直接绘出BodeBode图曲线。
图曲线 num=[1];den=[1 1 1];[z,p,k]=tf2zp(num,den);sys1=zpk(z,p,k)[a,b,c,d]=tf2ss(num,den);num2=[1];den2=[1 1 1 0];%带积分环节[z2,p2,k2]=tf2zp(num2,den2);sys2=zpk(z2,p2,k2)w=0.1:0.1:100;figure(1)subplot(2,2,1)nyquist(num,den)subplot(2,2,2)[r,img]=nyquist(num,den,w);plot(r,img)[g,p,wc,wp]=margin(num,den) %幅值和相位的稳定裕度subplot(2,2,3)nyquist(sys1,w);subplot(2,2,4)nyquist(a,b,c,d,1,w);g = Infp = 90wc = Infwp figure(2)subplot(2,1,1)nyquist(sys2,w)subplot(2,1,2)[rr,mm,w2]=nyquist(num2,den2,w);plot(rr,mm) 三、时间延迟系统三、时间延迟系统带延迟的系统相当于通常的系统后串接一个带延迟的系统相当于通常的系统后串接一个纯延迟系统：纯延迟系统： 对纯延迟环节的处理可用有理近似法（相当对纯延迟环节的处理可用有理近似法（相当于幂级数展开方式）：于幂级数展开方式）：• •[num,den]= pade(T,n)或或 [A,B,C,D]= pade(T,n) T为延迟时间常数，为延迟时间常数，n为要求拟合的阶数（阶为要求拟合的阶数（阶数越高，拟合精度越高），这样相当于用线性模数越高，拟合精度越高），这样相当于用线性模型（传递函数或状态方程形式）来近似代替纯延型（传递函数或状态方程形式）来近似代替纯延迟环节。
迟环节 4阶阶Transfer function:--------------------------------------tt=5;[nt1,dt1]=pade(tt,4);[nt2,dt2]=pade(tt,20);subplot(2,1,1),step(nt1,dt1,[0:0.1:100])subplot(2,1,2),step(nt2,dt2,100)20阶阶Transfer function: s^20 - 84 s^19 + 3511 s^18 - 9.691e004 s^17 + 1.977e006 s^16 - 3.163e007 s^15 + 4.112e008 s^14 - 4.441e009 s^13 + 4.041e010 s^12 - 3.125e011 s^11 + 2.063e012 s^10 - 1.163e013 s^9+ 5.581e013 s^8 - 2.267e014 s^7 + 7.706e014 s^6 - 2.158e015 s^5 + 4.855e015 s^4 - 8.453e015 s^3 + 1.071e016 s^2 - 8.791e015 s + 3.517e015 -----------------------------------------------------------------------------------------------s^20 + 84 s^19 + 3511 s^18 + 9.691e004 s^17 + 1.977e006 s^16 + 3.163e007 s^15 + 4.112e008 s^14+ 4.441e009 s^13 + 4.041e010 s^12 + 3.125e011 s^11 + 2.063e012 s^10 + 1.163e013 s^9+ 5.581e013 s^8 + 2.267e014 s^7 + 7.706e014 s^6 + 2.158e015 s^5 + 4.855e015 s^4+ 8.453e015 s^3 + 1.071e016 s^2 + 8.791e015 s + 3.517e015 20阶传函有理近似阶传函有理近似4阶传函有理近似阶传函有理近似实际纯滞后实际纯滞后 阶跃响应阶跃响应 num=[3 2 1];den=[1 5 7 2][nt2,dt2]=pade(tt,20);[num1,den1]=series(num,den,nt2,dt2)subplot(2,1,1)step(num,den,40)subplot(2,1,2)step(num1,den1,40) 6.4.5 6.4.5 根轨迹分析根轨迹分析 根轨迹是指，当 根轨迹是指，当开环系统开环系统某一参数从零变某一参数从零变到无穷大时，到无穷大时，闭环系统特征方程的根闭环系统特征方程的根在在s平面上平面上的轨迹。
一般来说，这一参数选作开环系统的的轨迹一般来说，这一参数选作开环系统的增益增益K，而在无零极点对消时，闭环系统特征方，而在无零极点对消时，闭环系统特征方程的根就是闭环传递函数的极点程的根就是闭环传递函数的极点 根轨迹分析方法是分析和设计线性定常控根轨迹分析方法是分析和设计线性定常控制系统的图解方法，使用十分简便利用它可制系统的图解方法，使用十分简便利用它可以对系统各种性能：稳定性、稳态性能、动态以对系统各种性能：稳定性、稳态性能、动态性能进行分析性能进行分析 一、常规根轨迹一、常规根轨迹rlocus(num,den) 直接绘出根轨迹直接绘出根轨迹[r,k]= rlocus(num,den) 返回增益和对应的返回增益和对应的 复极点复极点[r]= rlocus(num,den,k) 返回给定增益对应的返回给定增益对应的 复极点复极点注意：注意：这里输入的是这里输入的是开环开环传递函数，得到的是传递函数，得到的是闭环闭环系统极点随增益变化的轨迹。
系统极点随增益变化的轨迹 num=[1];den=conv([1 3 0],[1 2 2]);subplot(2,1,1)rlocus(num,den)subplot(2,1,2)sys1=tf(num,den);rlocus(sys1,'r')[r,k]=rlocus(num,den);rr=rlocus(num,den,100)rr 1 1、零极点图绘制、零极点图绘制、零极点图绘制、零极点图绘制 Ø[p,z]=[p,z]=pzmap(a,b,c,dpzmap(a,b,c,d) ) 返回状态空间描述系统的极 返回状态空间描述系统的极 返回状态空间描述系统的极 返回状态空间描述系统的极点矢量和零点矢量，而不在屏幕上绘制出零极点图点矢量和零点矢量，而不在屏幕上绘制出零极点图点矢量和零点矢量，而不在屏幕上绘制出零极点图点矢量和零点矢量，而不在屏幕上绘制出零极点图Ø[p,z]=[p,z]=pzmap(num,denpzmap(num,den) ) 返回传递函数描述系统的极 返回传递函数描述系统的极 返回传递函数描述系统的极 返回传递函数描述系统的极点矢量和零点矢量，而不在屏幕上绘制出零极点图。
点矢量和零点矢量，而不在屏幕上绘制出零极点图点矢量和零点矢量，而不在屏幕上绘制出零极点图点矢量和零点矢量，而不在屏幕上绘制出零极点图Øpzmap(a,b,c,dpzmap(a,b,c,d) )或或或或pzmap(num,denpzmap(num,den) ) 不带输出参数项， 不带输出参数项， 不带输出参数项， 不带输出参数项，则直接在则直接在则直接在则直接在s s复平面上绘制出系统对应的零极点位置，极点复平面上绘制出系统对应的零极点位置，极点复平面上绘制出系统对应的零极点位置，极点复平面上绘制出系统对应的零极点位置，极点用用用用× ×表示，零点用表示，零点用表示，零点用表示，零点用o o表示Øpzmap(p,zpzmap(p,z) ) 根据系统已知的零极点列向量或行向量 根据系统已知的零极点列向量或行向量 根据系统已知的零极点列向量或行向量 根据系统已知的零极点列向量或行向量直接在直接在直接在直接在s s复平面上绘制出对应的零极点位置，极点用复平面上绘制出对应的零极点位置，极点用复平面上绘制出对应的零极点位置，极点用复平面上绘制出对应的零极点位置，极点用× ×表示，表示，表示，表示，零点用零点用零点用零点用o o表示。
表示其它命令：其它命令：其它命令：其它命令： z=[1;2];p=[3;5;7];k=2;sys=zpk(z,p,k)[pp,zz]=pzmap(sys)pzmap(sys,'r')Zero/pole/gain: 2 (s-1) (s-2)-----------------(s-3) (s-5) (s-7)pp = 3 5 7zz = 1 2 ２、２、rlocfind()函数函数找出给定的一组根（闭环极点）对应的根轨迹增益找出给定的一组根（闭环极点）对应的根轨迹增益找出给定的一组根（闭环极点）对应的根轨迹增益找出给定的一组根（闭环极点）对应的根轨迹增益q[k,p]=[k,p]=rlocfind(a,b,c,drlocfind(a,b,c,d) )或者或者或者或者[k,p]=[k,p]=rlocfind(num,denrlocfind(num,den) )q它要求在屏幕上先已经绘制好有关的根轨迹图然后，它要求在屏幕上先已经绘制好有关的根轨迹图然后，它要求在屏幕上先已经绘制好有关的根轨迹图然后，它要求在屏幕上先已经绘制好有关的根轨迹图然后，此命令将产生一个光标以用来选择希望的闭环极点。
命令此命令将产生一个光标以用来选择希望的闭环极点命令此命令将产生一个光标以用来选择希望的闭环极点命令此命令将产生一个光标以用来选择希望的闭环极点命令执行结果：执行结果：执行结果：执行结果：k k为对应选择点处根轨迹开环增益；为对应选择点处根轨迹开环增益；为对应选择点处根轨迹开环增益；为对应选择点处根轨迹开环增益；p p为此点处为此点处为此点处为此点处的系统闭环特征根的系统闭环特征根的系统闭环特征根的系统闭环特征根 q不带输出参数项不带输出参数项不带输出参数项不带输出参数项[k,p][k,p]时，同样可以执行，只是此时只时，同样可以执行，只是此时只时，同样可以执行，只是此时只时，同样可以执行，只是此时只将将将将k k的值返回到缺省变量的值返回到缺省变量的值返回到缺省变量的值返回到缺省变量ansans中 二、广义根轨迹 将变化的参数经变形后放到原来增益所处的将变化的参数经变形后放到原来增益所处的位置，从而得到等效传函，再画根轨迹位置，从而得到等效传函，再画根轨迹 １、开环零点变化的根轨迹１、开环零点变化的根轨迹２、开环极点变化的根轨迹２、开环极点变化的根轨迹３、零度根轨迹３、零度根轨迹 对非最小相位系统，如含正反馈环节的系统。
对非最小相位系统，如含正反馈环节的系统绘制时，取分子多项式的负值即可绘制时，取分子多项式的负值即可 广义根轨迹举例：求下列系统随开环零点参数广义根轨迹举例：求下列系统随开环零点参数TaTa变化的根轨迹变化的根轨迹给定开环传递函数给定开环传递函数解：系统闭环传递函数解：系统闭环传递函数闭环特征方程闭环特征方程同除以同除以5 5同除前同除前3 3项项得等效闭环特征方程得等效闭环特征方程按照按照 确定开环零极点位置就绘制出随确定开环零极点位置就绘制出随TaTa变化的根轨迹变化的根轨迹 g=sym('5*(1+Ta*s)/(s*(5*s+1))');%开环传函g1=symadd('1',g);gg=symdiv(g,g1);%闭环传函[ng,dg]=numden(gg);dx=solve(dg,‘Ta');[nn,dn]=numden(symmul('-1',dx));denn=sym2poly(nn);numn=sym2poly(dn);rlocus(numn,denn)闭环特征方程是闭环特征方程是注意分子分母颠倒注意分子分母颠倒 广义根轨迹举例：求下列系统随广义根轨迹举例：求下列系统随开环极点开环极点参数参数 P1 P1 变化的根轨迹变化的根轨迹给定开环传递函数给定开环传递函数解：系统闭环传递函数解：系统闭环传递函数闭环特征方程闭环特征方程等效闭环特征方程等效闭环特征方程按照按照 确定开环零极点位置就绘制出随确定开环零极点位置就绘制出随P1P1变化的根轨迹变化的根轨迹求解过程是否和按开求解过程是否和按开环零点变化的广义根轨迹一致环零点变化的广义根轨迹一致 brainstormbrainstorm function [mn,md]=mydxtf(G,X)%g原开环传函，x新变量，均为符号表达式G1=symadd('1',G);GG=simple(symdiv(G,G1));[nng,dng]=numden(GG);dxg=solve(dng,X);[mnn,mdn]=numden(symmul('-1',dxg));md=sym2poly(mnn);mn=sym2poly(mdn);等效传递函数变换的函数等效传递函数变换的函数figure(1)%随开环零点随开环零点ta变化的根轨迹变化的根轨迹g=sym('5*(1+Ta*s)/(s*(5*s+1))');%开环传函 [n1,d1]=mydxtf(g,'Ta');rlocus(n1,d1)figure(2)%随开环极点随开环极点p1变化的根轨迹变化的根轨迹gp=sym('1/(s*(s+2)*(s+p1))');[n2,d2]=mydxtf(gp,'p1');rlocus(n2,d2) 现代控制理论分析设计现代控制理论分析设计一、状态反馈设计一、状态反馈设计一、状态反馈设计一、状态反馈设计 线性系统模型为：线性系统模型为：线性系统模型为：线性系统模型为： 利用状态反馈利用状态反馈利用状态反馈利用状态反馈U(t)= U(t)= - - - - KX(t)+r(t) KX(t)+r(t) 其中其中其中其中K K为状态反馈为状态反馈为状态反馈为状态反馈矩阵，矩阵，矩阵，矩阵， r(t) r(t)为参考输入，为参考输入，为参考输入，为参考输入， C Au Byx kr使得系统闭环方程为：使得系统闭环方程为：使得系统闭环方程为：使得系统闭环方程为： 几个常用命令几个常用命令几个常用命令几个常用命令•M=ctrb(A,B) 系统的能控矩阵系统的能控矩阵 M=[B AB A2B … An-1B]•rank(M) 得到矩阵的秩，得到矩阵的秩，M的秩为的秩为n，则能控，则能控•N=obsv(A,C) 求取系统的能观矩阵求取系统的能观矩阵N=[C CA CA2 … CAn-1]T，， N的秩为的秩为n，则能观，则能观a=[-1 -2 1;1 0 1;-1 0 1];b=[1;1;1];c=[1 0 1];d=0;m=ctrb(a,b); %能控性矩阵n=rank(m); if n= =length(a) disp('The system state can be controlled!')else disp('The system state can not be controlled!')endm =1 -2 -2 1 2 -2 1 0 2n = 3The system state can be controlled!N=obsv(a,c)nn=rank(N) %能观性判断if nn==length(a) disp('The system state can be observed!')else disp('The system state can not be observed!')end N = 1 0 1 -2 -2 2 -2 4 -2 nn = 3 The system state can be observed! 二、单输入系统的极点配置•k=acker(A,B,p) 对于期望极点对于期望极点p，求出状态反馈增益阵，求出状态反馈增益阵ka=[0 1 0;0 -1 1;0 0 -2];b=[0;0;1];c=[10 0 0];d=0;%系统矩阵pp=[-2 -1+j -1-j];%期望极点[z,p,k]=ss2zp(a,b,c,d);g0=zpk(z,p,k)m=ctrb(a,b)n=rank(m) %能控性判断if n==length(a) disp('The system state can be controlled!') k=acker(a,b,pp)else disp('The system state can not be controlled!')end10-------------s (s+1) (s+2) m = 0 0 1 0 1 -3 1 -2 4n = 3The system state can be controlled!k = 4 3 1 %设计后的效果比较[z2,p2,k2]=ss2zp(a-b*k,b,c-d*k,d);g=zpk(z2,p2,k2)figure(1)subplot(2,1,1)pzmap(g0)subplot(2,1,2)pzmap(g)figure(2)subplot(2,1,1)step(g0,'b',10)subplot(2,1,2)step(g,'r',10) 10---------------------(s+2) (s^2 + 2s + 2) 三、单输出系统的观测器设计三、单输出系统的观测器设计•若下面系统完全可观，若下面系统完全可观，•则全阶观测器方程为则全阶观测器方程为 x C Au By C A B G- • •全阶观测器方程为全阶观测器方程为• • 适当设计观测器的输出误差反馈阵适当设计观测器的输出误差反馈阵G，就可，就可改变观测器的极点，达到期望的闭环性能。
改变观测器的极点，达到期望的闭环性能• •设计设计G时采用的时采用的MATLAB命令仍是命令仍是acker，但，但要注意调用时，系统矩阵要进行转置要注意调用时，系统矩阵要进行转置• •G=(acker(A’,C’,p))’ a=[1 0;0 0];b=[1;1];c=[2 -1];d=0;p=[-10 -10];N=obsv(a,c);%能观性判断nn=rank(N)if nn==length(a) g=acker(a',c',p) disp('The system state can be observed!')else disp('The system can not be observed!')endnn = 2The system state can be observed! 。