基于MATLAB对控制课程的辅助应用分析.doc

1、基于MATLAB对控制课程的辅助应用分析摘要：早期的控制系统设计可以由纸笔等工具容易地计算出来，如 Ziegler 与 Nichols 于1942年提出的 PID经验公式就可以十分容易地设计出来。随着控制理论的迅速发展，光利用纸笔以及计算器等简单的运算工具难以达到预期的效果，加之在计算机领域取得了迅速的发展，于是很自然地出现了控制系统的计算机辅助设计 (computer-aided control system design , CACSD)方法。随着MATLAB 语言出现以来，就深受控制领域学生和研究者的欢迎，已经成为控制界最流行、最有影响的通用计算机语言， MATLAB 作为主要程序设计语

2、言来介绍控制系统计算机辅助设计的算法，可以使学者将主要精力集中在控制系统理论和方法上，而不是将主要精力花费在没有太大价值的底层重复性机械性劳动上，这样可以对控制系统计算机辅助设计技术有较好的整体了解，避免“只见树木，不见森林”的认识偏差，提高控制器设计的效率和可靠性。子曰:“工欲善其事，必先利其器”。跟踪国际最先进的CACSD软件环境及发展，以当前国际上最流行的 CACSD软件环境MATLAB为基本出发点来系统地介绍控制系统计算机辅助设计技术及软件实现，从而大大提高 CACSD算法研究与实际应用的效率和可靠性。关键词：Matlab 经典控制 现代控制 数学模型 传递函数 根轨迹 稳定性 劳斯判

3、据 奈氏图 伯德图 能控性 能观性 状态空间模型 正文：控制理论的主要内容主要分为两大部分：经典控制理论和现代控制理论。但是将这两者总结道一起又可以详尽的分为以下6个分类：1 数学模型经典（时域法）现代（频域法）理论基础建立在以1. 常微分方程稳定性理论2. Fourier变换为基础的根轨迹和奈奎斯特判据理论之上1. 常微分方程稳定性理论2. 状态空间分析3. 泛函分析、微分几何等现代数学分支数学模型传递函数(研究系统外部特性,属于外部描述,不完全描述。)状态空间表达式(深入系统内部,是内部描述,完全描述。)适用对象仅适用于: 单输入单输出线性定常集总参数可推广至: 多输入多输出非线性时变分布

4、参数传递函数表示方法MATLAB语言表达Num=b1,b2,，bm;Den=a1,a2,an;G=tf(num,den);G=ss(A,B,C,D)2. 稳定性分析及时域分析2.1对控制系统的性能的要求，主要是稳定性、暂态性能和稳态性能几个方面。系通过分析是系统设计的基础，而大部分的系统设计方法都是在系统稳定性的基础上发展起来的。常见的线性定常系统的时域分析方法有劳斯判据、赫尔维兹稳定判据等。线性定常连续系统稳定的充要条件是系统的全部特征根或闭环极点都具有负实部，或者说位于复平面左半部。劳斯判据不仅能够判别系统是否稳定，而且能够确定有多少正实部根，也能够具体确定对称于远点的特征根。线性系统的动

5、态性能取决于系统的闭环极点和零点的分布。虽然可以通过计算机直接求解得到闭环极点，但不能看到系统闭环极点随着系统参数变化的情况。但是通过根轨迹法，能够根据系统的开环零、极点分布，用图解的方法画出系统闭环极点随着系统参数变化的轨迹。 虽然用手工精确绘制系统更轨迹是非常困难的，但是Matlab中专门提供了绘制根轨迹有关的函数。r,k=rlocus(num,den)和r,k=rlocus(num,den,k)的功能是绘制根轨迹图。r,k=rlocus(num,den)是绘制部分的根轨迹。如果要以给定的参数范围绘制根轨迹，则执行命令k,poles=rlofind(num,den,p)。 k,poles=

6、rlofind(num,den,)和k,poles=rlofind(num,den,p)的功能是确定根轨迹上poles处的根轨迹放大系数的值。例：键入：n=1,3; d=1,13,54,82,60,0; rlocus(n,d)可得图1所示根轨迹图：图13.频域分析 当系统是高阶系统时，系统微分方程的求解时很困难的；另外，系统的时间响应没有明确反映出系统响应与系统结构、参数之间的关系，一旦系统很难满足要求，就很难确定如何去调整系统的结构和参数。频域分析法克服了时域的不足。根据系统的频率特性，可以直观地分析系统的稳定性，并且很容易的二和系统的结构、参数联系起来，因此可以根据系统的频率特性选择系统的

7、结构和参数，使之满足控制要求。实验中常用的频率分析方法是伯德图和奈氏图。因此，可以使用Matlab绘制系统的博得图、奈氏图，并确定系统的相位裕度和幅值裕度。 用Matlab很容易精确的绘制奈氏图和伯德图，但是只能绘制其中的一部分，不能绘制出的全貌，所以不能用来分析系统的稳定性。3.1用Matlab绘制奈氏图 绘制奈氏图的Matlab命令是nyquist(num,den),若要指定频率时，可用函数nyquist(num,den，)。还有以下两种形式：re,im,w= nyquist(num,den)re,im,w= nyquist(num,den,w)通过以上两种形式的调用，可计算频率特性的实部

8、和虚部，但不能产生奈氏图，还需调用plot（re,im）函数才能得到。例如：开环传递函数为用Matalb键入：G=tf(75*1,2,conv(1,10,1,3,2,5),nyquist(G)得到如图2所示：图23.2用Matlab绘制伯德图 绘制伯德图可用命令bode(num,den).如果需要给出频率的频率范围，可调用之林w=logspace(a,b,n),频率的采样点可在一定范围内产生n个十进制对数分点的等距离点。 如果需要制定幅值mag和phase范围，这执行命令mag,phase,w= bode(num,den)，Matlab在频率响应范围内能够自动选取频率值绘图。例：用Matlab

9、绘制伯德图 键入命令：G=tf(2000*1.5,conv(1,2,0,1,4,100),bode(G) 输出结果如下： Transfer function: 3000-s4 + 6 s3 + 108 s2 + 200 s 图3若键入G=tf(2000*1.5,conv(1,2,0,1,4,100);kg,r=margin(G)可以得到G的幅值域度kg = 0.8296和幅值穿越频率r =-14.40624.结构问题关于结构问题主要涉及现代控制理论，它包含了系统的能控能观性、以及对应的能控规范型和能观规范型、系统的能控分解和能观分解、系统的最小实现。线性系统的能控性和能观性事基于状态方程的控制

10、理论的基础.关于线性系统空间结构涉及的两个主要问题就是系统的状态空间模型的结构性分解以及传递函数阵与能控性/能观性的关系。而由系统的传递函数建立状态空间模型这类问题称为系统实现问题,而求得的状态空间模型称为相应的传递函数的一个实现。两种系统实现方法-能控/能观规范形实现。4.1线性系统的可控性判定 基于 MATLAB 的判定方法：rank(T) 构造可控性判定矩阵: 例：离散状态方程的可控性MATLAB 求解判定矩阵判定矩阵构造方法4.1.1线性系统的能观性判定 4.1.1.1判定矩阵 4.1.1.2 MATLAB 求解 4.1.2可控标准型和能观标准型 使用下面两个转换函数时线输入如下函数s

11、scanform（）function Gs=sscanform(G,type)switch type case ctrl G=tf(G); Gs=; G.num1=G.num1/G.den1(1); % 传递函数归一化 G.den1=G.den1/G.den1(1); d=G.num1(1); G1=G; G1.ioDelay=0; G1=G1-d; num=G1.num1; den=G1.den1; n=length(G.den1)-1; A=zeros(n-1,1) eye(n-1); -den(end:-1:2); B=zeros(n-1,1);1; C=num(end:-1:2); D

12、=d; Gs=ss(A,B,C,D,Ts,G.Ts,ioDelay,G.ioDelay); case obsv Gc=sscanform(G,ctrl); Gs=ss(Gc.a,Gc.c,Gc.b,Gc.d,Ts,G.Ts,ioDelay,G.ioDelay); otherwise error(Only options ctrl and obsv are applicable.)end可控标准型 Gs=sscanform(G,obsv)例如：对求解可控标准型 键入命令：num=6 0 2 8 10; den=2 0 6 4 8;G=tf(num,den); Gs=sscanform(G,obs

13、v)标准型：可观测标准型 4.1.2.1Kalman 规范分解子空间 示意图 4 图44.1.2.2最小实现例：多变量模型 是否能最小实现？MATLAB求解A=-6,-1.5,2,4,9.5; -6,-2.5,2,5,12.5; -5,0.25,-0.5,3.5,9.75; -1, 0.5, 0, -1, 1.5; -2, -1, 1, 2, 3; B=6,4; 5,5; 3,4; 0,2; 3,1; D=zeros(2);C=2,0.75,-0.5,-1.5,-2.75; 0,-1.25,1.5,1.5,2.25;G=ss(A,B,C,D); G1=minreal(G) 可得如下：a = x

14、1 x2 x3 x1 -2.169 -1.967 0.1868 x2 0.1554 -0.1811 -0.2789 x3 0.1403 1.613 -1.65b = u1 u2 x1 -7.907 -5.5 x2 3.228 2.506 x3 2.461 5.046c = x1 x2 x3 y1 -0.8326 -0.1501 -0.04028 y2 -0.3842 0.2764 0.4348d = u1 u2 y1 0 0 y2 0 05系统设计 系统设计主要是基于输入输出模型和基于状态空间模型的计。由此可以使用Matlab设计出相应的算法直接设计，并进行整个系统的仿真分析。5.1基于输入输

15、出模型的设计方法5.1.1.超前滞后校正器设计方法串联超前滞后校正器，如图5 图5 超前滞后校正器1 ;1：超前滞后校正器零极点分配图与伯德图，如图6所示图65.1.2超前滞后校正器的设计方法基于剪切频率和相位裕度的设计方法式中，为校正器的增益根据超前滞后校正器的涉及规则，可以写出相应的Matlab语言的超前滞后校正器的设计函数leadladc()，内容如下function Gc=leadlagc(G,Wc,Gam_c,Kv,key)G=tf(G); Gai,Pha=bode(G,Wc);Phi_c=sin(Gam_c-Pha-180)*pi/180);den=G.den1; a=den(le

16、ngth(den):-1:1);ii=find(abs(a)0, a=a(ii(1)+1); else, a=a(1); end;alpha=sqrt(1-Phi_c)/(1+Phi_c); Zc=alpha*Wc; Pc=Wc/alpha; Kc=sqrt(Wc*Wc+Pc*Pc)/(Wc*Wc+Zc*Zc)/Gai; K1=G_n*Kc*alpha/a;if nargin=4, key=1; if Phi_c0, key=2; else, if K1Kv, key=3; end, endendswitch keycase 1, Gc=tf(1 Zc*Kc,1 Pc);case 2 Kc=1

17、/Gai; K1=G_n*Kc/a; Gc=tf(1 0.1*Wc,1 K1*Gcn(2)/Kv); case 3 Zc2=Wc*0.1; Pc2=K1*Zc2/Kv; Gcn=Kc*conv(1 Zc,1,Zc2); Gcd=conv(1 Pc,1,Pc2); Gc=tf(Gcn,Gcd);end例：， 设计超前滞后校正器。键入命令：s=tf(s); G=4*(s+1)*(s+0.5)/s/(s+0.1)/(s+2)/(s+10)/(s+20)wc=20; f1=figure; f2=figure; % 打开两个图形窗口for gam=20:10:90 Gc=leadlagc(G,wc,ga

18、m,1000,3); figure(f1); step(feedback(G*Gc,1),1); hold on figure(f2); bode(Gc*G); hold on;endGc1=zpk(leadlagc(G,20,60,1000,3) % 设计超前滞后校正器Gc2=zpk(leadlagc(G,20,60,1000,1) % 设计超前校正器figure;step(feedback(G*Gc1,1),-,feedback(G*Gc2,1),:,1) % 绘制闭环响应结果如下：Transfer function: 4 s2 + 6 s + 2-s5 + 32.1 s4 + 263.2

19、 s3 + 426 s2 + 40 s其中为超前滞后校正器，为超前校正器，如图7，图8，图9所示。如果要知道相位裕度为60，不同剪切频率对应的单位阶跃响应和伯德图，键入如下命令即可，如图10、11、12、13所示。 图7 图8图9 图10 图11 图12 图13gam=60; f1=figure; f2=figure; % 打开两个图形窗口for wc=5:5:30 Gc=leadlagc(G,wc,gam,1000,3); a,b,c,d=margin(Gc*G); figure(f1); step(feedback(G*Gc,1),3); hold on figure(f2); bode(

20、Gc*G); hold on;end5.2基于状态空间模型的设计方法 使用Matlab辅助设计状态空间模型主要是对状态空间模型控制器的设计。5.2.1状态反馈控制结构框图如图14 所示 图14 将u(t)=v(t)-Kx(t)代入开环系统的状态方程模型，则在状态反馈矩阵K下，系统的闭环状态方程模型可以写成如果系统(A,B)完全可控，则选择合适的K矩阵，可以将闭环系统矩A-BK的特征值配置到任意地方。对于开环状态空间模型（A，B，C，D）来说，在状态反馈向量K下，闭环系统的状态方程可以写成（A-BK，B，C-DK，D）。就极点配置来说，使用以下算法可以实现，算法如下：Ackermann 算法 ：

21、K=acker（A,B,p）;鲁棒极点配置算法：K=place（A，B，p）;place( ) 函数不适用于含有多重期望极点的问题acker( ) 函数可以求解配置多重极点的问题例： ， 键入命令：A=0 1 0 0 ; 0 0 -1 0; 0 0 0 1; 0 0 5 0; B=0 1 ; 0 -1; 0 0 ; 0 0; p=-0.1; -0.2; -0.5+0.2i; -0.5-0.2i; K=place(A,B,p) rank(ctrb(A,B)5.2.2观测器设计及基于观测调节器设计 如图15所示图15原系统的（A，C）为完全可观测，则状态观测器数学模型的状态空间表示为：方程解析解为

22、:稳定，输入.m函数，如下function xh,x,t=simobsv(G,L)y,t,x=step(G); G=ss(G); A=G.a; B=G.b; C=G.c; D=G.d; y1,xh1=step(A-L*C),(B-L*D),C,D,1,t);y2,xh2=lsim(A-L*C),L,C,D,y,t); xh=xh1+xh2;G为对象的状态方程的状态对象模型，L为观测器向量。例：。键入命令：P=-10;-10;-10;-10; L=acker(A,C,P); L % 设计新观测器xh,x,t=simobsv(ss(A,B,C,D),L);plot(t,x,t,xh,:); set

23、(gca,XLim,0,30,YLim,-0.5,4)figure结果如图16所示：图16如果极点位于10键入命令：P=-10;-10;-10;-10; L=acker(A,C,P); L % 设计新观测器xh,x,t=simobsv(ss(A,B,C,D),L);plot(t,x,t,xh,:); set(gca,XLim,0,30,YLim,-0.5,4)可得结果：ans = -421.1634 260.7255 33.2946 -20.8091带有观测器的状态反馈控制结构图如图17所示：图17求取控制器和反馈环节的Matlab函数function Gc,H=obsvsf(G,K,L)H=

24、ss(G.a-L*G.c, L, K, 0);Gc=ss(G.a-G.b*K-L*G.c+L*G.d*K, G.b, -K, 1);如果参考输入信号r(t )=0,控制结构Gc(s)化简为,c=reg(G,K,L)则有Gc=reg(G,K,L)例：， 选择加权矩阵Q=diag(0.01,0.01,2,3),R=1,设计LQ最有控制器。键入命令：A=0,2,0,0; 0,-0.1,8,0; 0,0,-10,16; 0,0,0,-20;B=0;0;0;0.3953; C=0.09882,0.1976,0,0; D=0;Q=diag(0.01,0.01,2,3); R=1; % 输入加权矩阵K=lq

25、r(A,B,Q,R), step(ss(A-B*K,B,C,D) % 设计rm LQ 最优控制器结果如下：K =0.1000 0.9429 0.7663 0.6387如图18所示：设观测器极点均位于-5，用极点配置发设计极点观测器：键入命令：P=-5;-5;-5;-5; G=ss(A,B,C,D); L=acker(A,C, P); % 设计观测器Gc,H=obsvsf(G,K,L); % 设计控制器step(ss(A-B*K,B,C,D),feedback(G*Gc,H) % 比较基于观测器的控制器与状态反馈结果如图19 所示：图18图19基于观测器下控制器闭环系统的最小实现模型：键入命令：

26、zpk(minreal(feedback(G*Gc,H) % 最小实现模型zpk(minreal(ss(A-B*K,B,C,D) % 这个结果和上式忽略两个一阶零点几乎一致结果如下：Zero/pole/gain: 9.9982 (s+1)-(s+20.01) (s+10.01) (s2 + 0.3341s + 0.05052)Zero/pole/gain: 9.9982 (s+1)-(s+20.01) (s+10.01) (s2 + 0.3341s + 0.05052)参考文献：1 薛定宇 控制系统计算机辅助设计MATLAB语言与应用（第二版）.北京：清华大学出版社2 王万良 自动控制原理 .北京：高等教育出版社3 刘豹，唐万生 现代控制理论（第三版）.北京：机械工业出版社