仿真5线性系统分析设计课件

第五章线性系统分析与设计

MATLAB的控制工具箱是MATLAB最早的工具箱之一，也是控制系统的计算机辅助设计中最为流行的设计工具。控制工具箱适用于线性时不变系统(LTI)，可实现线性系统时域或频域的分析、设计和建模。可处理连续系统，也可处理离散系统；可使用经典或现代的技术。5.1线性系统的描述

MATLAB只处理矩阵这一种数学形式，各种控制系统的描述必须使用矩阵来表达。

MATLAB中主要采用传递函数和状态空间表达式。MATLAB中传递函数的描述方法

传递函数有三种常用形式：（1）一般形式（2）零极点增益形式

（3）部分分式形式（1）传递函数的一般形式

传递函数用分子、分母多项式表示，即num和den两个向量num=[bmbm-1…b1b0],den=[1an-1…a1a0]还可用SYS=TF(NUM,DEN)建立tf对象模型。

num=[123];den=[2234];yy=tf(num,den)Transferfunction:s^2+2s+3-----------------------2s^3+2s^2+3s+4（3）部分分式描述法

在传递函数没有相同极点时与部分分式相互转换：

[r,p,k]=residue(num,den)%部分分式展开 [num,den]=residue(r,p,k)%部分分式拟合状态空间表达式对SISO，A是n*n维系统矩阵，B是n*1维输入列向量，C是1*n维输出行向量，D是1*1维的直接传递矩阵。从仿真的角度来看，有时，仅仅实现系统输入与输出之间的关系是不够的，还必须实现模型内部变量，即状态变量，因此仿真要求采用系统内部模型，可采用状态空间表达式。状态方程输出方程

在MATLAB中，这个系统写为A、B、C、D四个矩阵的形式即可，当然矩阵维数要匹配。

也可用SYS=SS(A,B,C,D)建立ss模型，

SYS=SS(A,B,C,D,Ts)建立离散ss模型。典型的反馈控制系统结构图基本环节通常由各种联接关系来构成复杂系统：

串联series

并联parallel

反馈feedbackcloop（2）并联连接[A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2)传递函数形式[num,den]=parallel(num1,den1,num2,den2)（3）反馈连接[A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign)[num,den]=feedback(num1,den1,num2,den2,sign)

[A,B,C,D]=cloop(A1,B1,C1,D1,sign)单位反馈[num,den]=cloop(num1,den1,sign)sign反馈极性，正反馈1，负反馈-1或缺省，MATLAB中线性系统模型之间的转换

ss—状态空间、tf—传递函数、zp—零极点：[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)部分分式到传函1)c2d使用零阶保持器离散化，只有状态空间形式

2)c2dm既有状态空间形式，又有传递函数形式；3)参数ts是采样周期T；

4)method指定转换方式，其中“zoh”表示采用零阶保持器；“foh”表示采用三角形近似；“tustin”表示采用双线性变换；“prewarp”表示采用指定转折频率的双线性变换;系统默认为零阶保持器法。MATLAB中连续系统与离散系统之间的转换c2d连续ss到离散ssd2c离散ss到连续ssc2dm连续到离散(可用不同方法)d2cm离散到连续(可用不同方法)5.2时间响应分析step()连续系统的单位阶跃响应

dstep()离散系统的单位阶跃响应

impulse()连续系统的单位脉冲响应

dimpulse()离散系统的单位脉冲响应

lsim()连续系统的任意输入响应

dlsim()离散系统的任意输入响应

initial()连续系统的零输入响应

dinitial()离散系统的零输入响应

5.2时间响应分析一、单位阶跃响应1.对状态空间描述的系统y=step(A,B,C,D,iu,T)

求ss系统对iu个输入的在T时间向量下的单位阶跃响应，y是列向量。[y,t,x]=step(A,B,C,D,iu,T)同上,并返回状态x的变化step(A,B,C,D,iu,T)

直接画出单位阶跃响应曲线没有返回参数，直接画出单位阶跃响应曲线返回参数指定时间序列总结单位阶跃响应step调用时是否带输出变量

不带输出变量，直接绘出响应曲线step(sys1)带输出变量，得到响应值[y,t,x]=step(sys1)是否指定仿真持续时间不指定，缺省时根据系统的零极点自动给出指定：两种形式终止时间或时间向量序列

step(sys1,Tfinal)step(sys1,0:dt:Tfinal)模型形式各种sysobject如由tfsszpk建立传递函数numden状态空间表达式

ABCD%阻尼系数对典型二阶系统的阶跃响应的影响a=[00.30.713];b=0.2;T=0:100;fori=1:5num=b^2;den=[12*a(i)*b(i)b^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,[10]);%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型有误差sys12(s+1)(s+0.5)---------------------(s+0.1)(s+1.5)(s+2)sys22(s+1)(s+0.5)---------------------s(s+0.1)(s+1.5)(s+2)离散时间系统例如二阶离散时间系统num=[2-3.41.5],den=[1-1.60.8]dstep(num,den)3.对各种形式描述的系统直接画出响应曲线

impulse(sys1,’s1’,sys2,’s2’,…)

sys(i)可以是各种形式描述的LTI系统，时间是自动给出的，也可以用时间向量和给出最终时间两种方法确定。4.对离散系统dimpulset=0:0.1:10;wn=1.5;z=0.2;num1=[02];den1=[101];A=[01;-wn^2-2*z*wn];B=[0;wn^2];C=[10];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,'k',5)Transferfunction:2--------10s+1a=x1x2x101x2-2.25-0.6b=u1x10x22.25c=x1x2y110d=u1y10Zero/pole/gain:10(s+1)(s+2)-----------------(s+3)(s+4)(s+5)四、initial(sys,x0)零输入响应（自由运动）线性定常系统带初始状态的阶跃响应——全响应step函数自由运动强制运动initial函数根据线性定常系统的叠加性，由step和initial函数得到的结果进行叠加实现[y1,t1,x1]=step(sys,T)[y2,t2,x2]=initial(sys,x0,T)subplot(3,1,3)plot(t1,x1+x2)legend('x1','x2')根据线性定常系统的叠加性，由step和initial函数得到的结果进行叠加实现[y1,t1,x1]=step(sys,T)[y2,t2,x2]=initial(sys,x0,T)subplot(3,1,3)plot(t1,y1+y2)还可由lsim函数得到[y,t,x]=lsim(sys,u,t,x0)subplot(3,1,3)tt=0:0.1:20;nn=length(tt);u(1:nn)=1;lsim(sys,u,tt,x0)用符号工具箱求解a=[-0.5572-0.7814;0.78140];b=[1;0];c=[1.96916.4493];d=0x0=[1;0],T=20;symstE=expm(a*t)x0t=E*x0h=0.1;T=20;tt=0:h:T;fori=1:length(tt)t=tt(i);xx0(:,i)=eval(x0t);end

plot(tt,xx0)legend('x1','x2')5.3频域分析法当系统G(s)的输入为正弦函数系统的稳态响应为：对于稳定的线性定常系统，在谐波输入作用下，其系统输出的稳态分量仍然是与输入同频率的谐波函数，而幅值和相位的变化是频率的函数，且与系统数学模型相关。相位之差为系统的相频特性为此，定义谐波输入下，系统的响应中与输入同频率的谐波分量与谐波输入的幅植之比为系统的幅频特性为系统的频率特性。对于稳定的线性系统，当系统的输入为频率为w的正弦信号，其稳态输出也是频率为w的正弦信号。输出与输入的幅值比为系统的幅频特性，输出和输入的相位差为系统的相频特性。产生频率向量频率向量由logspace()函数来构成，调用格式：w=logspace(m,n,npts)生成以10为底的指数向量（10m

~10n），点数由npts选定。w=logspace(-2,2,10)w=0.010.02780.07740.21540.59951.66814.641612.915535.9381100类似地，在线性空间也有生成向量的函数

linspace(-2,2,10)ans=-2-1.5556-1.1111-0.6667-0.22220.22220.66671.11111.55562

linspace(0,1,11)ans=00.100.200.300.400.500.600.700.800.901.00x=0:0.1:1x=00.100.200.300.400.500.600.700.800.901.00一、Bode图伯德图即对数频率特性图，由对数幅频特性图和对数相频特性图组成，是工程中广泛使用的一种图示方法。对系统的频率特性取以10为底的对数，得

令为系统的对数幅频特性，也称为系统的增益，单位为贝尔（B）或分贝(dB)，且有1B=20dB。同时为了和对数幅频特性相一致，也称为对数相频特性（注意没有取对数）。对数幅频特性图和对数相频特性图均以为横坐标分度。

对数幅频特性图的纵坐标按线性分度，单位是分贝(db)。线性分度，单位为度或弧度对数相频特性图的纵坐标按由此构成的坐标系均为半对数坐标系。用MATLAB绘制Bode图bode(num,den,w)没有返回参数,直接绘出bode图。

包括了对数幅频特性图和对数相频特性图。横坐标为频率w，采用对数分度，单位为弧度/秒；纵坐标均匀分度，分别为幅值函数20lgA(w)，以dB表示；相角，以度表示。[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=[111];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)grida=[0.20.40.60.813];b=6*ones(1,6);w=logspace(0,2)fori=1:6num=b(i)^2;den=[12*a(i)*b(i)b(i).^2];[m(:,i)p(:,i)]=bode(num,den,w);endmm=20*log(m');pp=p';subplot(2,1,1)semilogx(w,mm(1,:),‘r’,w,mm(2,:),‘b’,w,mm(3,:),‘m’,w,mm(4,:),‘g’,…w,mm(5,:),'k',w,mm(6,:),'c')legend('a=0.2','a=0.4','a=0.6','a=0.8','a=1','a=3')subplot(2,1,2)semilogx(w,pp(1,:),'r',w,pp(2,:),'b',w,pp(3,:),'m',w,pp(4,:),'g',w,pp(5,:),'k',w,pp(6,:),'c')legend('a=0.2','a=0.4','a=0.6','a=0.8','a=1','a=3')二.Nyquist曲线奈奎斯特图即极坐标图，主要用于对闭环系统稳定性。又称为幅相频率特性图或幅相曲线。对于任意给定的频率，频率特性值为复数。若将频率特性表示为实部和虚部的形式，即则实部为实频特性，虚部为虚频特性。MATLAB中绘制Nyquist曲线

对于频率特性函数G(jw)，给出w从负无穷到正无穷的一系列数值，分别求出Im(G(jw))和Re(G(jw))。以Re(G(jw))为横坐标，Im(G(jw))为纵坐标绘制成为极坐标频率特性图。

Nyquist(num,den,{wmin,wmax})

直接绘出Nyquist曲线，{wmin,wmax}是频率范围，缺省自动给定w从负无穷到正无穷[real,imag]=

Nyquist(num,den,{wmin,wmax})得到实部、虚部，用plot绘制w从零到正无穷变化部分。 另外，以上两种可用ss描述(A,B,C,D)；或几种方法建立的sys所描述的模型来作图。

%典型环节wmin=0.01;wmax=100;sys1=tf([2],[1])%比例环节sys2=tf([2],[101])%惯性环节sys21=tf([2],[-101])%不稳定惯性环节sys3=tf([1],[10])%积分环节sys4=tf([10],[1])%微分环节wn=0.1;a=[0.10.31];fori=1:3sys5(i)=tf(wn^2,[12*a(i)*wnwn^2])%不同阻尼比的二阶振荡环节end2--------10s+12-2--------10s-11---ss0.01-----------------------s^2+0.02s+0.010.01-----------------------s^2+0.06s+0.010.01-----------------------s^2+0.2s+0.01margin函数求幅值和相位稳定裕度

margin函数可以从频率响应数据中计算出幅值裕度、相角裕度以及对应的频率。幅值裕度和相角裕度是针对开环SISO系统而言，它指出系统闭环时的相对稳定性。margin(num,den)当不带输出变量引用时，margin可在当前图形窗口中绘制出带有裕量及相应频率显示的Bode图，其中幅值裕度以分贝为单位。[Gm,Pm,Wcg,Wcp]=margin(num,den)

计算出系统wcp，Gm幅值裕度，Pm相位裕度，及相应的相角交界频率Wcg、截止频率Wcp而不直接绘出Bode图曲线。

num=[1];den=[111];[z,p,k]=tf2zp(num,den);sys1=zpk(z,p,k)[a,b,c,d]=tf2ss(num,den);num2=[1];den2=[1110];%带积分环节[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=1.0000figure(2)subplot(2,1,1)nyquist(sys2,w)%axis([-200-105])subplot(2,1,2)[rr,mm,w2]=nyquist(num2,den2,w);plot(rr,mm)figure(3)margin(num,den)例某反馈系统的开环传递函数为其中，要求判断闭环系统是否稳定。析：正实部极点数=0第二步：作出系统的开环极坐标图：第一步：把开环比例系数增大为100，要求重新判断闭环系统的稳定性。研究图可知，当从连续增大到时，曲线不包围点，符合Nyquist稳定判据的要求，所以闭环系统是稳定的！k1=20;num=k1;den=conv(conv([101],[21]),[0.21]);sys1=tf(num,den)k2=100;num2=k2;sys2=tf(num2,den)

figure(1)subplot(2,1,1)nyquist(sys1)subplot(2,1,2)nyquist(sys2)研究图可知，当从连续增大到时，曲线不包围点，符合Nyquist稳定判据的要求，所以闭环系统是稳定的！K=20正实部极点数仍为P=01)当从连续增大到时，的曲线顺时针方向包围点2圈，即．由于，闭环系统不稳定。2)根据知道闭环系统具有两个右半平面极点K=100k1=20;num1=k1;den=conv(conv([101],[21]),[0.21]);k2=100;num2=k2;numf=1;denf=1;[numb1,denb1]=feedback(num1,den,numf,denf)[numb2,denb2]=feedback(num2,den,numf,denf)sysb1=tf(numb1,denb1)sysb2=tf(numb2,denb2)p1=roots(denb1)p2=roots(denb2)时间延迟系统

带延迟的系统相当于通常的系统后串接一个纯延迟系统：

对纯延迟环节的处理可用有理近似法（相当于幂级数展开方式）：[num,den]=pade(T,n)或[A,B,C,D]=pade(T,n)

T为延迟时间常数，n为要求拟合的阶数（阶数越高，拟合精度越高），这样相当于用线性模型（传递函数或状态方程形式）来近似代替纯延迟环节。

4阶Transferfunction:s^4-4s^3+7.2s^2-6.72s+2.688--------------------------------------s^4+4s^3+7.2s^2+6.72s+2.688tt=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阶Transferfunction:

s^20-84s^19+3511s^18-9.691e004s^17+1.977e006s^16-3.163e007s^15+4.112e008s^14-4.441e009s^13+4.041e010s^12-3.125e011s^11+2.063e012s^10-1.163e013s^9+5.581e013s^8-2.267e014s^7+7.706e014s^6-2.158e015s^5+4.855e015s^4-8.453e015s^3+1.071e016s^2-8.791e015s+3.517e015-----------------------------------------------------------------------------------------------s^20+84s^19+3511s^18+9.691e004s^17+1.977e006s^16+3.163e007s^15+4.112e008s^14+4.441e009s^13+4.041e010s^12+3.125e011s^11+2.063e012s^10+1.163e013s^9+5.581e013s^8+2.267e014s^7+7.706e014s^6+2.158e015s^5+4.855e015s^4+8.453e015s^3+1.071e016s^2+8.791e015s+3.517e01520阶传函有理近似4阶传函有理近似实际纯滞后阶跃响应num=[321];den=[1572][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)5.4rootlocus根轨迹根轨迹是指，当开环系统某一参数从零变到无穷大时，闭环系统特征方程的根在s平面上的轨迹。一般来说，这一参数选作开环系统的增益K，而在无零极点对消时，闭环系统特征方程的根就是闭环传递函数的极点——常规根轨迹

根轨迹分析方法是分析和设计线性定常控制系统的图解方法，使用十分简便。利用它可以对系统各种性能：稳定性、稳态性能、动态性能进行分析。闭环特征方程即根轨迹方程为G(s)H(s)=–1系统的闭环传递函数:C（s）R（s）－G(s)H(s)

模条件角条件常规根轨迹rlocus(num,den)

直接绘出根轨迹（缺省增益范围）rlocus(num,den,k)

直接绘出根轨迹（制定增益范围）[r,k]=rlocus(num,den)

返回增益和对应的复极点[r]=rlocus(num,den,k)

返回给定增益对应的复极点

注意：这里输入的是开环传递函数，得到的是闭环系统极点随增益变化的轨迹。

num=[1];den=conv([130],[122]);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=-3.5492+2.1028i-3.5492-2.1028i1.0492+2.1852i1.0492-2.1852i1、零极点图绘制[p,z]=pzmap(a,b,c,d)

返回状态空间描述系统的极点矢量和零点矢量，而不在屏幕上绘制出零极点图。[p,z]=pzmap(num,den)返回传递函数描述系统的极点矢量和零点矢量，而不在屏幕上绘制出零极点图。pzmap(a,b,c,d)或pzmap(num,den)

不带输出参数项，则直接在s复平面上绘制出系统对应的零极点位置，极点用×表示，零点用o表示。pzmap(p,z)根据系统已知的零极点列向量或行向量直接在s复平面上绘制出对应的零极点位置，极点用×表示，零点用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=357zz=12２、rlocfind()函数找出给定的一组根（闭环极点）对应的根轨迹增益[k,p]=rlocfind(a,b,c,d)[k,p]=rlocfind(num,den)[k,p]=rlocfind(sys)

它要求在屏幕上先已经绘制好有关的根轨迹图。然后，此命令将产生一个光标以用来选择希望的闭环极点。命令执行结果：k为对应选择点处根轨迹开环增益；p为此点处的系统闭环特征根。

不带输出参数项[k,p]时，同样可以执行，只是此时只将k的值返回到缺省变量ans中。G0=tf([240],[157])rlocus(G0)[k,p]=rlocfind(G0)Selectapointinthegraphicswindowselected_point=-1.6600+0.5568ik=0.6474p=-1.6537+0.5620i-1.6537-0.5620i3、sgrid()函数在根轨迹平面上绘制出阻尼比和等固有频率的网格。

sgrid

sgrid(z,wn)G0=tf([240],[157])figure(1)subplot(1,2,1)rlocus(G0)grid

subplot(1,2,2)rlocus(G0)sgrid(0:0.1:1,0:pi/10:pi)ImβRe(b)—衰减系数—阻尼振荡角频率—阻尼角广义根轨迹

将变化的参数经变形后放到原来增益所处的位置，从而得到等效传函，再画根轨迹。１、开环零点变化的根轨迹２、开环极点变化的根轨迹３、零度根轨迹 对非最小相位系统，如含正反馈环节的系统。绘制时，取分子多项式的负值即可。

广义根轨迹举例：求下列系统随开环零点参数Ta变化的根轨迹给定开环传递函数解：系统闭环传递函数闭环特征方程同除以5同除前3项得等效闭环特征方程按照确定开环零极点位置就绘制出随Ta变化的根轨迹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变化的根轨迹求解过程是否和按开环零点变化的广义根轨迹一致

brainstormfunction[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)5.5现代控制理论分析设计一、利用MATLAB求状态方程的解1、矩阵指数函数利用符号工具箱（SymbolicToolbox）和expm（）函数求解矩阵指数函数。

例1,求矩阵指数函数2、状态方程的解

线性系统满足叠加原理，把系统同时在初始状态和输入作用下的状态运动分解为：由初始状态和输入分别单独作用所产生的运动的叠加。A=[01;-2-3];B=[01]';u=1;x0=[00]';symstt1;%t1=taoE=expm(A*t)E2=subs(E,'t','t-t1');f=E2*B*u;%symmul(symmul(E2,B),u)JF=int(f,'t1','0','t')x=E*x0+JF%symadd(symmul(E,x0),JF)E=[-exp(-2*t)+2*exp(-t),exp(-t)-exp(-2*t)][-2*exp(-t)+2*exp(-2*t),2*exp(-2*t)-exp(-t)]JF=1/2-exp(-t)+1/2*exp(-2*t)exp(-t)-exp(-2*t)x=1/2-exp(-t)+1/2*exp(-2*t)exp(-t)-exp(-2*t)t=0:0.1:10;x1=eval(x(1));x2=eval(x(2))；figure(1)plot(t,x1,'r',t,x2,'k')legend('x1','x2')A=[01;-2-3];B=[01]';C=[10];D=0;sys=ss(A,B,C,D);[y,t,x]=step(sys);subplot(1,3,1)plot(t,x(:,1),'--',t,x(:,2))xlabel('Time(sec)')legend('x1(t)','x2(t)')subplot(1,3,2)plot(x(:,1),x(:,2))xlabel('x1');ylabel('x2')subplot(1,3,3)step(A,B,C,D)？直接用时域响应分析函数stepA=[01;-2-3];B=[01]';x0=[11]';symstt1;u=t;E=expm(A*t)E2=subs(E,'t','t-t1');f=E2*B*u;JF=int(f,'t1','0','t')x=E*x0+JF

t=0:0.1:10;x1=-2*exp(-2*t)+3*exp(-t)+1/2*t+1/2*t.*exp(-2*t)-t.*exp(-t);x2=-3*exp(-t)+4*exp(-2*t)-t.*(-exp(-t)+exp(-2*t));figure(1)plot(t,x1,'r',t,x2,'k')legend('x1','x2')E=[-exp(-2*t)+2*exp(-t),exp(-t)-exp(-2*t)][-2*exp(-t)+2*exp(-2*t),2*exp(-2*t)-exp(-t)]JF=1/2*t+1/2*t*exp(-2*t)-t*exp(-t)-t*(-exp(-t)+exp(-2*t))x=-2*exp(-2*t)+3*exp(-t)+1/2*t+1/2*t*exp(-2*t)-t*exp(-t)-3*exp(-t)+4*exp(-2*t)-t*(-exp(-t)+exp(-2*t))A=[01;-2-3];B=[01]';x0=[11]';C=[10];D=0;sys=ss(A,B,C,D);[y1,t1,x1]=lsim(sys,t,t);[y2,t2,x2]=initial(sys,x0,t);x=x1+x2;y=y1+y2;plot(t,x(:,1),'r',t,x(:,2),'k')legend('x1','x2')

？直接用时域响应分析函数step二、MATLAB中系统的能控性和能观性分析1、能控性和能观性的判定

在MATLAB中利用以下几个函数，可以直接从矩阵A、B、C中判断系统的能控性和能观性。

M=ctrb(A,B)

求取系统的能控性判别矩阵M=[BABA2B…An-1B]

N=obsv(A,C)

求取系统的能观性判别矩阵N=[CCACA2…CAn-1]T

rank()得到矩阵的秩2、按系统的能控性和能观性进行分解

（1）按能控性进行分解若系统是不完全能控的，则通过相似变换可进行结构分解，分解为能控子系统和不能控子系统。MATLAB提供分解函数[Ac,Bc,Cc,T,K]=ctrbf(A,B,C)T为相似变换阵，K是长度为n的矢量，其元素是各个块的秩，sum(K)得出系统能控部分的秩。判断是否完全能控，否则按能控性进行分解例：系统为

（2）按能观性进行分解若系统是不完全能观的，则通过相似变换可进行结构分解，分解为能观子系统和不能观子系统。MATLAB提供分解函数[Ao,Bo,Co,T,K]=obsvf(A,B,C)T为相似变换阵，K是长度为n的矢量，其元素是各个块的秩，sum(K)得出系统能观部分的秩。判断是否完全能观，否则按能观性进行分解例：系统为

三、MATLAB中系统的稳定性分析对于线性定常系统而言，一个连续时间系统的所有极点都具有负实部，则系统稳定。一个离散时间系统的所有极点都位于单位园内，则系统稳定。在MATLAB中，易于求出系统的极点。

判断系统稳定性，并给出不稳定极点例：系统为