《MATLAB语言与控制系统仿真课程设计六条腿步行器系统的校正与优化.doc》由会员分享，可在线阅读，更多相关《MATLAB语言与控制系统仿真课程设计六条腿步行器系统的校正与优化.doc（30页珍藏版）》请在三一办公上搜索。
1、 2008 级MATLAB语言与控制系统仿真课程大作业姓 名 * 学 号 * 所在院系 电气与电子工程 班 级 * 日 期 2011年1月12日 同组人员 作业评分 评阅人 六条腿步行器系统的校正与优化一、简介如图1所示，本报告的研究对象是一个六条腿自动行走的机器系统。由于该系统的特定用途，它需要准确而快速的响应，对带宽的要求严格，在控制系统中具有典型的代表性。因而本报告选择该系统进行分析，分别用三种方法来进行校正和优化：时域法进行串联校正；频率法进行超前校正；状态空间法进行极点配置。以达到举一反三，活学活用的效果。图1，六条腿步行器系统的实物图该系统的控制对象是六条腿的位移，通过输入信号来控
2、制系统的位移，其传递函数为，因为系统要求准确迅速和抗干扰能力强，故有如下性能指标要求：1、1型系统，阶跃响应的稳态误差为零；2、通频带带宽要大于1Hz;3、阶跃响应的超调量要小于15% ；4、响应要快，调整时间在1秒以内。首先用MATLAB描述步行器系统模型，并观察该系统的阶跃响应和频谱图。%描述系统模型，传递函数形式num0=1;den0=1 2 10 0; num,den=cloop(num0,den0);t=0:0.001:50;step(num,den,t) %求闭环模型阶跃响应图2,原系统闭环阶跃响应由响应图可知系统阻尼比大，响应很慢，调节时间38.5s。%求闭环系统Bode图sys
3、=tf(num,den)w=logspace(-1,2,400);bode(sys,w)margin(sys)图2,原系统闭环幅频相频曲线由图可以看出，系统的带宽不足。其通频带带宽可由如下程序算出：w=logspace(-1,2,400);mag,phase,w=bode(num,den,w);y,L=max(mag);wblist=find(mag0.707);wb=w(wblist(1)得结果 wb =0.1035，远小于6.28。为此，我们主要的目标是改善阻尼比和带宽。二、正文1、时域法串联校正a、设计方法指标1要求系统是一阶系统即可，指标2要求系统的自然频率比较大，由书中所给公式有，且
4、系统的阻尼比要小于0.707,这样能得到比较理想的带宽，指标3要求超调量PO15%，则要求阻尼比。最后，指标4要求系统有尽可能小的时间，为此在满足带宽的条件下，我们可以寻找到合适的阻尼比系数，程序如下：%寻找合适的阻尼比，阴尼比以damp表示damp=0.5:0.001:0.707;wB=6.28;Wn=wB./(-1.1961*damp+1.8508);Ts=4./(damp.*Wn);plot(damp,Ts)nb=find(damp=0.7);wn=Wn(nb)得结果wn =6.1962，带宽满足要求，由输出的Ts随阻尼比变化图可知，选择，能够得到比较快的响应速度。（注：教材中所选阻尼比
5、是0.52,是因为分析有误，错误之处为“in general , selecting larger results in a larger settling time”1）图3,系统调节时间Ts随阻尼比变化的关系本小节采用教材所用方案，意在对其优化，我们所加的串联校正器为。之所以采用此种PID控制器的变形，是因为，它有四个示知数，串联入系统中后，系统是四阶的，配置所需极点时正好能列出四个方程，于是可以得出合理的解，否则参数比方程多或者参数比方程少，都会给求解带来困难，而且常常是无解。在MATLAB中，可用以下形式表示控制器: numc=k*1 a b;denc=0 1 c;同教材，所需配置的闭
6、环极点特征方程式为；选取10,编程可得校正器的具体参数值：%求当阻尼比为0.7时的系统阶跃响应damp=0.7;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alph=10;%求解出四元方程解c=2*damp*wn*(1+alph)-2k=wn2*(1+alph*damp2*(4+alph)-10-2*ca=(2*alph*damp*wn3*(1+damp2*alph)-10*c)/kb=alph2*damp2*wn4/k得到结果为：c = 93.4210 k =2.4753e+003 a =7.5609 b =29.1786 接着我们得到校正后闭环系统的阶跃响应图。
7、%建系统闭环模型numc=k*1 a b;denc=0 1 c;num0=1;den0=1 2 10 0;num,den=series(numc,denc,num0,den0);n,d=cloop(num,den);t=0:0.002:1.5;step(n,d,t)得到结果响应图：图4,阻尼比等于0.7时的闭环系统阶跃响应验证频率特性是否满足条件，求带宽：%求其带宽w=logspace(-1,2,400);mag,phase,w=bode(n,d,w);%得到其幅频图semilogx(w,mag),gridylabel(Mag),xlabel(Frequency)y,L=max(mag);wb
8、list=find(mag0.707);wb=w(wblist(1)Mp=ywr=w(L)得到结果为wb =40.6464 Mp =1.2165 wr =10 带宽为40多，满足要求。其幅频图如下：图5, 阻尼比等于0.7时的闭环系统幅频图b、结果分析与探讨与教材所选阻尼比相比较可得，系统响应的超调量相差不大，但当阻尼比等于0.7时，系统的响应速度要快很多。图6,教材中阻尼比等于0.52时的闭环系统阶跃响应但我们可以发现一个问题，超调量与我们的理论值不相符，比我们的理论大一些。：“一个不能忽略的零点对系统的影响是使超调量加大，响应速度加快，这是由于零点具有微分的作用；一个不能忽略的极点对系统的
9、影响是使超调量减小，调节时间增加，这是由于极点的滤波作用（或称为阻尼作用）。2”由此，我们可知，是串联的校正器引入了零点的原因，使超调量大于理论值。可以用仿真来验证这一结论的正确性，另一比较对象是只包含该系统极点，不含零点的闭环系统。%求不引入零点的闭环系统的响应，用作比较damp=0.7;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alph=10;%配置的极点模型dspole=conv(1 2*damp*wn wn2,conv(1 alph*damp*wn,1 alph*damp*wn);kdc=dcgain(1,dspole); %求直流增益，让稳态误差为零t=
10、0:0.002:1.5;step(1/kdc,dspole,t)图7,不含零点的闭环系统响应图 由图7所示，不含零点时其超调量为4.49%，符合理论值，且反应速度变慢，可知前述结论是正确的。c、优化设计以前述理论为基础，我们可以对改变系统的极点来优化校正方案，具体表现为：闭环极点特征方程式中，通过改变的大小，可以改变闭环系统的极点位置，由于所添加的极点有阻尼作用，可以通过控制添加极点的位置，来达到改善超调量及响应时间的目的。令分别等于2，4，8，10，15，20，25，30，35，40，45，50，求出超调量PO、调节时间Ts及带宽wb和谐振峰值Mp，列成表以供分析比较用。程序如下：%子函数部
11、分%解四元方程得到校正器参数再得到系统闭环模型function n,d=model(alph); %定义子函数global damp wn %全局变量阻尼比和自然频率%求校正器参数c=2*damp*wn*(1+alph)-2;k=wn2*(1+alph*damp2*(4+alph)-10-2*c;a=(2*alph*damp*wn3*(1+damp2*alph)-10*c)/k;b=alph2*damp2*wn4/k;%求系统闭环模型numc=k*1 a b;denc=0 1 c;num0=1;den0=1 2 10 0;num,den=series(numc,denc,num0,den0);
12、n,d=cloop(num,den);%主函数部分%当变化时求系统的超调量PO、调节时间Ts及带宽wb和谐振峰值Mpglobal damp wn %定义全局变量阻尼比和自然频率，并赋值，供子函数调用damp=0.52;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alph=2 4 8 10 15 20 25 30 35 40 45 50; %的取值for i=1:1:12 %循环求值t=0:0.0001:4;n,d=model(alph(i); %调用子函数得到系统闭环模型y=step(n,d,t); %求系统闭环模型的阶跃响应数据PO(i)=100*(max(y)-
13、1); %求超调量%求调节时间TsL=find(abs(y-1)0.02); Ts(i)=t(L(length(L); w=logspace(-1,2,400);%求系统闭环模型的幅频特性mag,phase,w=bode(n,d,w);y,L=max(mag);wblist=find(mag0.707); %求带宽wb(i)=w(wblist(1);Mp(i)=y; %求谐振峰值wr(i)=w(L); %求谐振频率endPO,Ts,wb,Mp,wrfigure(1) %图1得到PO、TsAX1,H11,H12=plotyy(alph,PO,alph,Ts);gridset(get(AX1(1)
14、,Ylabel),String,PO)set(get(AX1(2),Ylabel),String,Ts)set(H11,Marker,*,LineStyl,none)set(H12,Marker,d,LineStyl,none)figure(2) %图2得到Mp、wbAX2,H21,H22=plotyy(alph,wb,alph,Mp);gridset(get(AX2(1),Ylabel),String,wb)set(get(AX2(2),Ylabel),String,Mp)set(H21,Marker,*,LineStyl,none)set(H22,Marker,d,LineStyl,non
15、e)以下为输出图：图8,PO与Ts随变化时的取值图9,wb与Mp随变化时的取值由图分析可得，当增大时，调节时间减小，超调减小，且带宽增加，性能会有改善，所以选一些的值，会有比较优良的性能。同时我们也可看出，当大于40以后，对超调量与调节时间的影响已经并不大了，增大作用不大，且太大时，校正器增益比较大，不易实现。所以，选择，能得到比较理想的性能。取，适当修改前面的程序即可得到其性能参数，在此不再赘述。其阶跃响应图为：图10, 时的闭环系统阶跃响应超调量为10%，调节时间为0.41秒。其校正器参数为c =109.6146，k =3.1868e+003，a =5.2387，b =23.1521。幅频
16、特性参数为：Mp =1.1391，wr =6.8326，wb =41.3563。可见性能很优越。如果实际系统还需改进的话，可以进一步减小的取值范围进行精确计算查找。2、频域法串联校正由前面的PID控制器串联校正法，我们可以观察到控制器的开环增益非常大，往往要加到近千倍的增益，这在实际系统中是很不容易实现的。例如：我们增加第5条性能指标要求校正环节的开环增益不能超过250此时，我们不能仅仅依靠增加开环系统的增益来达到增加带宽的目的，本小节将运用超前校正来解决这一问题。因为超前校正能够使系统的相位超前，故引入超前校正系统能够在增加系统开环带宽的同时也增加系统的相位裕度（相位裕度与阻尼比直接相关，相
17、位裕度加大，表示阻尼比系数也加大，意味着超调减小），能够满足系统的要求。a、设计方法超前校正环节的通用表达式为，本系统的开环带宽比较小，所以要先增加开环增益k以增大带宽，然后再进行相位超前校正。编程序先找到原系统在6.28处的增益，然后以此为依据进行补偿，过程如下%增加系统开环增益，使带宽增加sys0=tf(1,1 2 10 0);w=logspace(-1,2,400);mag0,phase0,w=bode(sys0,w);wB=6.28;mk=spline(w,mag0,wB); %插值得到穿越频率为6.28处的幅值k=1/mk %确定校正环节的增益值sysk=tf(k*sys0);mag
18、k,phasek=bode(sysk,w);bode(sysk,w) %得到幅相频特性图，找到相角裕度margin(sysk)图11, 增加开环增益后的Bode图程序得结果k =200.9966，校正环节的开环增益在要求的范围内。且由输出Bode图可以观察到系统的相位裕度为-66.9度，不稳定。由于超调量要求比较高，我们为达到比较高的阻尼比，需将目标相位裕度取为70度，由此可见，一个超前校正装置无法满足要求，三个串联的方案才可行。其校正环节应为。编程进行校正：%超前校正过程dsPm=70; %目标相位裕度Pmk=-67; %增加开环增益后的开环系统相位裕度Pmcha=dsPm-Pmk+10Pm
19、a=Pmcha*pi/180/3; %每个超前环节应该校正的相位裕度量alph=(1+sin(Pma)/(1-sin(Pma); MdB=20*log10(magk);am=-30*log10(alph);wgc=spline(MdB,w,am); %插值，寻找到实际的穿越频率点T=1/(wgc*sqrt(alph); num2,den2=series(alph*T 1,T 1,alph*T 1,T 1);num3,den3=series(num2,den2,alph*T 1,T 1); %三超前环节组成的校正器模型参数num,den=series(k,1 2 10 0,num3,den3);
20、n,d=cloop(num,den)figure(1)margin(num,den) %校正后系统开环Bode图figure(2) step(n,d) %校正后闭环系统的阶跃响应图得到图形如下：图12, 校正后系统开环Bode图图13, 校正后闭环系统的阶跃响应图运用程序得到各环节的Bode图，了解各环节的作用。%系统各环节的Bode图sys=tf(num,den);sysc=tf(num3,den3);%原系统sys0，增加开环增益k后系统sysk，校正环节sysc，加校正后的整体系统sysbode(sys0,b,sysk,g,sysc,r,sys,k)legend(sys0,sysk,sy
21、sc,sys)grid 图14, 校正后系统各环节的Bode图b、结果分析与探讨由结果图形可知，增加超前校正环节后，系统的开环带宽及相度都得以增大，系统的超调量得到改善。校正的目标相位裕度为70度，对应的阻尼比接近0.7,超调量应该比较小，但是由阶跃响应图反映系统的超调量是11.7%，比设想的要大。原因是校正环节仍然引入了增加超调量的极点（验证方法与时域校正法中的验证环节相同，在此不必再说），校正环节仍然无法避免不了此种影响的出现。而且，超前校正系统的形式单一，并不能得到理想的极点配置以得到好的性能。可见超前校正还是有诸多限制和不足的。3、状态空间法极点任意配置前两种方法都能改善系统的性能，但
22、都不能任意配置所需的极点，因而常常会有一些局限。运用状态空间法进行极点任意配置能够很好地解决这一问题（前提是系统可控）。a、设计方法运用状态反馈进行极点配置，本身并不会增加或减小原系统的阶数，因而本题中系统的理想极点不能再设定为前面所述的了，原系统是三阶的，为了方便研究，仍旧取一对主导极点，其特征方程为，还有一个相对远离虚轴的极点，特征方程为，在此可取，所需配置的系统极点特征方程为，其性能指标与时域法中的分析差别不大，因为这个极点影响作用不大。用MATLAB描述如下：%求取所需配置的系统极点多项式damp=0.52;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alp
23、h=10;dspole=conv(1 2*damp*wn wn2,1 alph*damp*wn)得到结果为dspole=1.0000 31.8899 167.3631 694.0791进行极点配置需前先判断可控性：%判断系统的可控性num0=1;den0=1 2 10 0;A,B,C,D = tf2ss(num0,den0)n=length(A)%CO=B A*B A2*B;%CO=ctrb(A,B);rCO=rank(CO)if rCO=ndisp(系统可控)elseif rCOndisp(系统不可控)end由结果可知此系统是可控的，即可进行极点任意配置。JA=poly(A); %矩阵A的特
24、征多项式a1=JA(2);a2=JA(3);a3=JA(4);W=a2 a1 1;a1 1 0;1 0 0;P=CO*W; %求取变换矩阵aa1=dspole(2);aa2=dspole(3);aa3=dspole(4); %期望极点特征多项式与原系统极点特征多项式相减K=aa3-a3 aa2-a2 aa1-a1*(inv(P) %求取反馈矩阵得到结果为状态反馈矩阵为K =29.8899 157.3631 694.0791num,den=ss2tf(A-B*K,B,C,D)Kpot=dcgain(num,den); %求取直流增益，确定系统应该增加的频率figure(1) %图1为极点配置后系
25、统阶跃响应图step(A-B*K,B/Kpot,C,D,1,t)figure(2) %图2为极点配置后闭环系统的Bode图numnew,dennew=ss2tf(A-B*K,B/Kpot,C,D);bode(numnew,dennew)得到极点配置后系统的传递函数为：num = 0 -0.0000 -0.0000 1.0000den = 1.0000 31.8899 167.3631 694.0791图15,极点配置后系统阶跃响应图图16,极点配置后闭环系统阶跃Bode图b、结果分析与探讨由结果可知理想极点位置的特征多项多为：dspole=1.0000 31.8899 167.3631 694
26、.0791所配置极点传递函数为：num = 0 -0.0000 -0.0000 1.0000den = 1.0000 31.8899 167.3631 694.0791两相比较可发现是一致的，并且极点配置法不会引入不必要的零点，因而系统的性能参数与理论值是相符合的（从图中可知超调量为14.5%，调节时间为1.56秒，与理论值一致）。我们可以任意配置所需极点，使系统的性能指标达到我们所期望的值，这是极点任意配置的优越之处。但是，状态反馈在应用中也有很多限制，实际系统中有很多系统都是组装好封装的，不能得到所设定的状态，或者是状态在物理电路中不能测量或测量的成本太高，这些都使得状态反馈法并不能尽如人
27、意。c、优化设计由前结果分析，可知，系统的状态往往不能在物理电路中测量出来或者是测量的成本太高，此时，就有必要构造一个观察器了。构造系统观察器，先得判断系统的可观性：%判断系统的可观性num0=1;den0=1 2 10 0;A,B,C,D = tf2ss(num0,den0)n=length(A)OB=obsv(A,C);rOB=rank(OB)；if rOB=ndisp(系统可控)elseif rOBndisp(系统不可控)可得结果，系统是可控的。由于观察器要求响应速度比原系统要快，一般观察器的极点距虚轴的位置是原系统极点的五倍以上。JA=poly(A); %矩阵A的特征多项式roots(
28、JA)得原系统的开环极点为ans = 0 、-1.0000 + 3.0000i 、-1.0000 - 3.0000i；因而，我们可取观察器的极点为-5+5j ， -5-5j ，-10。配置过程如下：%求取观察器增益矩阵a1=JA(2);a2=JA(3);a3=JA(4);W=a2 a1 1;a1 1 0;1 0 0;P=W*OB; %求取变换矩阵 J=-5+5j 0 0;0 -5-5j 0;0 0 -10; %期望极点JJ=poly(J)aa1=JJ(2);aa2=JJ(3);aa3=JJ(4); Ke=(inv(P)*aa3-a3;aa2-a2;aa1-a1 %求取反馈矩阵得到结果为Ke =
29、112.0000；104.0000；18.0000。由于观察器的具体实现涉及到具体的物理电路实现，已经超出了本课题所研究的范围，故在此不必深入探讨，只就自动控制原理在MATLAB上的方法实现作一介绍即可。三、总结本课题运用时域分析法进行串联校正，频率分析法进行超前校正，状态空间分析法进行极点配置，将所学的自动控制理论知识比较充分地运用了出来，有利于提高自控知识的活学活用水平。本次课题针对实际系统的分析研究，让我对各种控制校正方法有了一个比较清晰全面的认识，了解到了在工程实际中各种方法的优点及缺陷，有益于完成将来的实际工程工作。同时，本次课题研究也让我意识到，MATLAB作为了一款强大的工程应用
30、软件，在处理实际工程问题中发挥了重要的作用，本次研究只是一个比较肤浅的基础尝试，将来还必须有更深入的研究，MATLAB的魅力才能更突出地显示出来。致谢和参考文献感谢张蓉老师平时的授课指导，让我学到了很多关于MATLAB的基础知识，并启发我将理论运用于实际，让我的能力有了比较大的提高。感谢华中科技大学王永骥、王敏老师编著的自动控制原理一书，书中的自动控制的原理和方法给了我很多启发，本文中状态空间的计算方法就是根据该书的相关知识而来。参考文献：1Robert H.Bishop .现代控制系统分析与设计应用MATLAB和Simulink .北京.清华大学出版社 20032 王永骥，王金城，王敏.自动
31、控制原理.北京：化学工业出版社 2007附录1、简介分析中所用到的MATLAB代码：%描述系统模型，传递函数形式num0=1;den0=1 2 10 0; num,den=cloop(num0,den0);t=0:0.001:50;step(num,den,t) %求闭环模型阶跃响应%求闭环系统Bode图sys=tf(num,den)w=logspace(-1,2,400);bode(sys,w)margin(sys)%计算通频带带宽w=logspace(-1,2,400);mag,phase,w=bode(num,den,w);y,L=max(mag);wblist=find(mag0.70
32、7);wb=w(wblist(1)2、时域分析法校正环节中所用到的MATLAB代码：%寻找合适的阻尼比，阴尼比以damp表示damp=0.5:0.001:0.707;wB=6.28;Wn=wB./(-1.1961*damp+1.8508);Ts=4./(damp.*Wn);plot(damp,Ts)nb=find(damp=0.7);wn=Wn(nb)numc=k*1 a b;denc=0 1 c; %串联控制器的描述%求当阻尼比为0.7时的系统阶跃响应damp=0.7;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alph=10;%求解出四元方程解c=2*damp*
33、wn*(1+alph)-2k=wn2*(1+alph*damp2*(4+alph)-10-2*ca=(2*alph*damp*wn3*(1+damp2*alph)-10*c)/kb=alph2*damp2*wn4/k%建系统闭环模型numc=k*1 a b;denc=0 1 c;num0=1;den0=1 2 10 0;num,den=series(numc,denc,num0,den0);n,d=cloop(num,den);t=0:0.002:1.5;step(n,d,t)%求其带宽w=logspace(-1,2,400);mag,phase,w=bode(n,d,w);%得到其幅频图se
34、milogx(w,mag),gridylabel(Mag),xlabel(Frequency)y,L=max(mag);wblist=find(mag0.02); Ts(i)=t(L(length(L); w=logspace(-1,2,400);%求系统闭环模型的幅频特性mag,phase,w=bode(n,d,w);y,L=max(mag);wblist=find(mag0.707); %求带宽wb(i)=w(wblist(1);Mp(i)=y; %求谐振峰值wr(i)=w(L); %求谐振频率endPO,Ts,wb,Mp,wrfigure(1) %图1得到PO、TsAX1,H11,H12
35、=plotyy(alph,PO,alph,Ts);gridset(get(AX1(1),Ylabel),String,PO)set(get(AX1(2),Ylabel),String,Ts)set(H11,Marker,*,LineStyl,none)set(H12,Marker,d,LineStyl,none)figure(2) %图2得到Mp、wbAX2,H21,H22=plotyy(alph,wb,alph,Mp);gridset(get(AX2(1),Ylabel),String,wb)set(get(AX2(2),Ylabel),String,Mp)set(H21,Marker,*,
36、LineStyl,none)set(H22,Marker,d,LineStyl,none)3、频率法超前校正环节中所用到的MATLAB代码：%增加系统开环增益，使带宽增加sys0=tf(1,1 2 10 0);w=logspace(-1,2,400);mag0,phase0,w=bode(sys0,w);wB=6.28;mk=spline(w,mag0,wB); %插值得到穿越频率为6.28处的幅值k=1/mk %确定校正环节的增益值sysk=tf(k*sys0);magk,phasek=bode(sysk,w);bode(sysk,w) %得到幅相频特性图，找到相角裕度margin(sysk
37、)%超前校正过程dsPm=70; %目标相位裕度Pmk=-67; %增加开环增益后的开环系统相位裕度Pmcha=dsPm-Pmk+10Pma=Pmcha*pi/180/3; %每个超前环节应该校正的相位裕度量alph=(1+sin(Pma)/(1-sin(Pma); MdB=20*log10(magk);am=-30*log10(alph);wgc=spline(MdB,w,am); %插值，寻找到实际的穿越频率点T=1/(wgc*sqrt(alph); num2,den2=series(alph*T 1,T 1,alph*T 1,T 1);num3,den3=series(num2,den2
38、,alph*T 1,T 1); %三超前环节组成的校正器模型参数num,den=series(k,1 2 10 0,num3,den3);n,d=cloop(num,den)figure(1)margin(num,den) %校正后系统开环Bode图figure(2) step(n,d) %校正后闭环系统的阶跃响应图%系统各环节的Bode图sys=tf(num,den);sysc=tf(num3,den3);%原系统sys0，增加开环增益k后系统sysk，校正环节sysc，加校正后的整体系统sysbode(sys0,b,sysk,g,sysc,r,sys,k)legend(sys0,sysk,
39、sysc,sys)grid 4、状态空间分析法极点配置环节中所使用到的MATLAB环节：%求取所需配置的系统极点多项式damp=0.52;wB=6.28;wn=wB/(-1.1961*damp+1.8508);alph=10;dspole=conv(1 2*damp*wn wn2,1 alph*damp*wn)%判断系统的可控性num0=1;den0=1 2 10 0;A,B,C,D = tf2ss(num0,den0)n=length(A)%CO=B A*B A2*B;%CO=ctrb(A,B);rCO=rank(CO)if rCO=ndisp(系统可控)elseif rCOndisp(系统不可控)endJA=poly(A); %矩阵A的特征多项式a1=JA(2);a2=JA(3);a3=JA(4);W=a2 a1 1;a1 1 0;1 0 0;P=CO*W; %求取变换矩阵aa1=dspole(2);aa2=dspole(3);aa3=dspole(4); %期望极点特征多项式与原系统极点特征多项式相减K=aa3-a3 aa2-a2 aa1-a1*(inv(P) %求取反馈矩阵num,den=ss2tf(A-B*K,B,C,D)Kpot=dcgain(num,den); %求取直流增益，确定系统应该增
