1 / 56

第 六 章

第 六 章. MATLAB 在控制系统中的应用. §6 .1 控制系统的数学描述与建模. 一、微分方程 常微分方程是控制系统模型的基础。 MATLAB 提供了两个求解微分方程数值解的函数，龙格 — 库塔法 ode23 、 ode45 （采用 2 阶和 4 阶龙格 — 库塔公式）。 例：电路图如图所示， R=1.4 Ω,L=2H,C=0.32F 。初始状态：电感电流为零，电容电压为 0.5V 。 t=0 时接入 1V 的电压。求 0<t<15s 时 i(t) 和 v(t) 的值，并且画出电流和电容电压的关系曲线。 根据电流定律有：.

odessa
Download Presentation

第 六 章

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 第 六 章 MATLAB在控制系统中的应用

  2. §6.1 控制系统的数学描述与建模 一、微分方程 常微分方程是控制系统模型的基础。MATLAB提供了两个求解微分方程数值解的函数，龙格—库塔法ode23、ode45（采用2阶和4阶龙格—库塔公式）。 例：电路图如图所示，R=1.4Ω,L=2H,C=0.32F。初始状态：电感电流为零，电容电压为0.5V。t=0时接入1V的电压。求0<t<15s时i(t)和v(t)的值，并且画出电流和电容电压的关系曲线。 根据电流定律有： 1.4 Ω 2H 令： x1=v0, x2=i 得： 0.32F Vs=1V

  3. 方程的M文件定义如下： ELECTSYS.M function xdot=electsys(t,x) v=1; R=1.4; L=2; C=0.32; xdot=[x(2)/C;1/L*(V-x(1)-R*x(2))]; 下面的M文件使用ode23对系统进行仿真。 SIMU.M t0=0; tfinal=15; %时间间隔为0~15 x0=[0.5,0]; %初始化 tol=0.001; %精度 trace=0 %若非零则打印出每一步的计算值 [t,x]=ode23('electsys',t0,tfinal,x0,tol,trace); subplot(211)

  4. plot(t,x) title('Time Response of a RLC series ciruit'); xlabel('Time-sec'); text(8,1.05,'Capatitor Voltage'); text(8,0.05,'current'); vc=x(:,1); i=x(:,2); subplot(212); plot(vc,i) title('current versus capatitor voltage'); xlabel('Capacitor Voltage'); ylabel('current'); subplot(111);

  5. 二、传递函数 1. 多项式和特征多项式的根 root(P) poly(r) poly2sym(p) 例：求多项式s6+9s5+31.25s4+61.25s3+67.25s2+14.75s+15的根 多项式系数以降幂次序排列在一个行向量中，用roots( )函数求根。 P=[1 9 31.25 61.25 67.25 14.75 15] r=roots(P) 例：多项式的根为-1，-2，-3+j4，-3-j4。求多项式方程。 r=[-1 –2 –3+j4 –3-j4] p=poly(r) poly2sym(p) 例：求下列矩阵的特征方程的根。 A=[0 –1 –1 ；-6 –11 6；-6 –11 5] p=poly(A); r=roots(p); eig(A); poly2sym(p)

  6. 2. 传递函数的零点和极点 tf2zp( ) 例：求下列传递函数的零点、极点和增益。 程序如下： num=[1 11 30 0]; den=[1 9 45 87 50]; [z,p,k]=tf2zp(unm,den) 例：已知系统的零点为-6，-5，0，极点为-3+4j，-3-4j，-2，-1，增益 为1，求其传递函数。 程序如下： z=[-6 –5 0]’; k=1; p=[-3+4*j,-3-4*j,-2,-1]’; [num,den]=zp2tf(z,p,k)

  7. 3. 部分分式展开 [r,p,k]=residue(b,a) 把 分解为微分单元和的形式，即： 例：对F(s)进行部分分式展开。 程序如下： b=[2 0 9 1]; a=[1 1 4 4]; [r,p,k]=residue(b,a) 三、状态空间描述 1. 将微分方程化为状态方程 2. 矩阵的对角化

  8. 四、模型的转换与连接 1. 数学模型的转换（见下表： ）

  9. 对部分常用函数用法举例如下： (1) ss2tf—实现由状态空间描述到传递函数形式的转换。 [num,den]=ss2tf(A,B,C,D,iu) 例：将下列状态空间描述化为函数模型 程序如下： a=[0 1 –1;-6 –11 6;-6 –11 5]; b=[0 0 1]’; c=[1 0 0]; d=0 [num,den]=ss2tf(a,b,c,d) (2) ss2zp—由状态空间形式转化为零极点增益形式。 [z,p,k]=ss2zp(A,B,C,D,iu) 转换后的传函：

  10. 例：已知二输入二输出系统的状态空间模型为：例：已知二输入二输出系统的状态空间模型为： 求其零极点模型。 程序如下： a=[0 1 –1;-6 –11 6;-6 –11 5]; b=[0 0 ;0 1;1 0]’; c=[1 0 0;0 1 0]; d=[2 0;0 2]; [num,den]=ss2zp(a,b,c,d,1); [num,den]=ss2zp(a,b,c,d,2); (3) tf2ss—由传递函数形式转化为状态空间形式。 [A,B,C,D]=tf2ss(num,den) 例：已知系统传递函数如下所示： 求其状态空间模型。

  11. 程序如下： num=[0 0 –2;0 –1 –5; 1 2 0]; den=[1 6 11 6]; [A,B,C,D]=tf2ss(num,den) 2. 系统模型的连接（见下表： ）

  12. 对部分函数的用法举例如下： （1）append( )—将两个状态空间描述的动态特性连接在一起。 [A,B,C,D]=append(A1,B1,C1,D1,A2,B2,C2,D2) (2) series( )—将两个状态空间形式表示的系统进行串联。 [A,B,C,D]=series(A1,B1,C1,D1,A2,B2,C2,D2) (3) parallel( )—两个状态空间系统的并联。 [A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2) [num,den]=parallel(num1,den1,num2,den2) (4) feedback( )—用于两个系统的反馈连接。 [A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2) [A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign) [num,den]=feedback(num1,den1,num2,den2) [num,den]=parallel(num1,den1,num2,den2,sign)

  13. u y 例：两个子系统如下所示： g(s) - h(s) 按右图方式连接，求闭环系统的传递函数。 程序如下： numg=[2 5 1]; deng=[1 2 3]; numh=[5 10]; denh=[1 10]; [num,den]=feedback(numg,deng,numh,denh) 五、控制系统的模型属性 MATLAB提供了许多用来分析控制系统模型属性的函数，例如常 用的可控性，可观测性、阻尼系数、自然频率等，如下表：

  14. 对部分函数用法举例如下： （1）gram( ),dgram( )—其可控、可观测的格来姆(Gram)矩阵。 Gc=gram(A,B,’c’) Gc=dgram(A,B,’c’) Go= gram(A,B,’o’) Gc=dgram(A,B,’o’) （2）ctrb( ),obsv( )—求系统的可控性和可观测性。 Ac=ctrb(A,B) Ao=obsv(A,C) 例：一线性系统如下所示： 当α分别取-1，0，1时，判断系统的可控性和可观测性，并求出相 应的状态方程。 程序如下面的M文件所示： CANDO.M %判断系统在不同参数下的可控性和可观测性

  15. for alpha=[-1 0 2] alpha num=[1 alpha]; den=[1 -15 14 0]; [A,B,C,D]=tf2ss(num,den) Ac=ctrb(A,B); rAc=rank(Ac) Ao=obsv(A,C); rAo=rank(Ao) end (3) damp( ),ddamp( )—求系统的阻尼系数和自然频率。其用法如下： [Wn,z]=damp(A) mag=ddamp(A） %得到A的特征值向量。 [mag,Wn,z]=ddamp(A,Ts) %Ts为抽样时间间隔。

  16. damp( ),ddamp( )函数用于计算自然频率和阻尼系数，当不带输出变量时，可在屏幕上显示出特征值表、阻尼系数和自然频率。A变量可以取三种形式，如下表： 例1：求系统的特征值、阻尼系数和自然频率。 A=[1 2 –1 2;2 6 3 0;4 7 –8 –5;7 2 1 6] damp(A)

  17. 例2：某离散系统如下所示： 求其特征值、幅值、阻尼系数和自然频率。 程序如下： num=[3 –5.1 6]; den=[4 3.2 –1.9 8]; ddamp(den,0.01) （4）printsys( )—显示或打印线性系统，其用法如下： printsys(A,B,C,D) printsys(A,B,C,D,ulabels,ylabels,xlabels) printsys(unm,den,’s’) printsys(unm,den,’z’)

  18. 例： A=[1 2 –1 2;2 6 3 0;4 7 –8 –5;7 2 1 6]; B=[-1 0 0 1]’; C=[-2 5 6 1]; D=8; printsys(A,B,C,D) or: printsys(A,B,C,D,’current’,’voltage’,’Vr I1 Vc I’) or: [num,den]=ss2tf(A,B,C,D) printsys(num,den) (5) dcgain( ),ddcgain( )—计算控制系统的稳态增益，其用法如下： k=dcgain(A,B,C,D); k= ddcgain(A,B,C,D); k= dcgain(num,den); k= ddcgain(num,den); 例： A=[1 2 –1 2;2 6 3 0;4 7 –8 –5;7 2 1 6]; B=[-1 0 0 1]’; C=[-2 5 6 1]; D=8; k= dcgain(A,B,C,D)

  19. 六、控制系统常用数学方程求解 1）are( )——求解里卡第(ATx+xB=-C)方程，其用法为： x=are(A,B,C) 例：求解某系统的代数里卡第方程的解。 A=[1 2 –1 2;2 6 3 0;4 7 –8 –5;7 2 1 6]; B=eye(size(A)); C=[1 2 3 0;2 4 4 6;3 4 1 5;0 6 5 0]; x=are(A,B,C) 2）lyap( ),lyap2( ),dlyap( )——求解李亚普诺夫方程。其用法如下： x=lyap(A,B,C) 可求解连续系统的一般形式的李亚普诺夫方程： Ax+xB=-C (A,B,C应当有适当的维数) x=lyap(A,C)可求解特殊形式的李亚普诺夫方程： Ax+xAT=-C x=dlyap(A,B,C)用于求解离散系统的李亚普诺夫方程：

  20. lyap2( )与lyap( )的功能相同，但采用了特征值分解技术求解李亚普 诺夫方程，因此一般来讲，运算速度比较快。 例：求解上例系统中的李亚普诺夫方程。 A=[1 2 –1 2;2 6 3 0;4 7 –8 –5;7 2 1 6]; B=eye(size(A)); C=[1 2 3 0;2 4 4 6;3 4 1 5;0 6 5 0]; x=lyap(A,B,C)

  21. §6.2 控制系统的分析方法 一、控制系统的稳定性分析 由控制理论的一般规律可知：对线性系统而言，如果一个连续系统所有的极点都位于s平面的左半平面，则该系统为一个稳定系统。对离散系统而言，如果一个系统的全部极点都位于z平面的单位圆内部，则该系统是一个稳定系统。 从控制理论又可知，所谓最小相位系统首先是指一个稳定的系统，同时对于连续系统而言，系统的所有零点都位于s平面的左半平面，即零点实部小于零；对于一个离散系统而言，系统的所有零点都位于z平面的单位圆内部。 所以，只要知道了系统的模型，不论哪种形式的数学模型，都可以很方便的由MATLAB求出系统的零点、极点，从而判断系统的稳定以及判断系统是否为最小相位系统。 例1：已知某控制系统的模型为：

  22. 要求判断系统的稳定性以及系统是否为最小相位系统。要求判断系统的稳定性以及系统是否为最小相位系统。 程序如下： function sstable(A,B,C,D) % 判断系统的稳定性以及系统是否为最小相位系统 [z,p,k]=ss2zp(A,B,C,D); %求出系统的零点、极点 ii=find(real(z)>0); n1=length(ii); jj=find(real(p)>0); n2=length(jj); if (n1>0) disp('The System is Unstable.');

  23. else disp('The System is stable.'); if (n2>0) disp('The System is a Nonminimal Phase One.'); else disp('The System is a Minimal Phase One.'); end end 在程序中使用了find( )函数，find函数用法：ii=find(条件) 例2：已知系统模型如下所示： 判断系统的稳定性，以及系统是否为最小相位系统。 程序如下：

  24. function sstable(A,B,C,D) % 判断系统(A,B,C,D)的稳定性以及系统是否为最小相位系统 if(nargin==2) [z,p,k]=tf2zp(A,B) %求出系统的零点、极点 elseif(nargin==4) [z,p,k]=ss2zp(A,B,C,D); %求出系统的零点、极点 end ii=find(real(z)>0); n1=length(ii); jj=find(real(p)>0); n2=length(jj); if (n1>0) disp('The System is Unstable.');

  25. else disp('The System is stable.'); if (n2>0) disp(‘But the System is a Nonminimal Phase One.'); else disp('The System is a Minimal Phase One.'); end end 二、控制系统的时域分析 1. 时域分析的一般方法 一个动态系统的性质常用典型输入下的响应来描述，响应是指零初 始条件下某种典型的输入作用下对象的响应，控制系统常用的输入 为单位阶跃函数和脉冲激励函数。在MATLAB的控制系统工具箱 中提供了求取两种输入下系统典型响应的函数step( )和impulse( )。 step( )—求取系统的阶跃响应。其用法：

  26. [y,x]=step(num,den,t) [y,x]=step(A,B,C,D,iu,t) t为选定的仿真时间，y为系统在仿真时刻各个输出所组成的矩阵，x为自动选择的状态变量的时间响应数据。 绘制系统解阶跃响应曲线，可由如下格式完成： step(num,den,t) or step(num,den) step(A,B,C,D,iu,t) or step(A,B,C,D,iu) 例1：已知系统的开环传递函数如下所示： 求出该系统在单位反馈下的阶跃响应曲线。 首先求出系统的闭环传递函数模型：

  27. 对系统进行仿真，执行下面的M文件： EXP4.M num0=20; den0=[1 8 36 40 0]; t=1:0.1:10; numc=num0; denc=[zeros(1,length(den0)-length(num0)),num0]+den0; numc denc [y,T,x]=step(numc,denc,t); plot(t,y,t,x); title('The Step Response') xlabel('Time_Sce'); text(4,1.5,'The Output'); text(5,6.5,'The State')

  28. impulse( )—求取脉冲响应的函数，用法与step( )函数基本一致。 例子见EXP4.M。 2. 常用的时域分析函数 （见下表）

  29. 部分函数的用法举例如下： (1) initial( )—求连续系统的零输入响应。用法如下： [y,x,t]=iniial(A,B,C,D,x0) [y,x,t]=iniial(A,B,C,D,x0,t) 例：某三阶系统如下所示： 当初始状态x0=[1 0 0]T时，求该系统的零输入响应。 执行程序如下： A=[1 –1 0.5;2 –2 0.3;1 –4 –0.1]; B=[0 0 1]’; C=[0 0 1]; D=0;

  30. x0=[1 0 0]’; t=0:0.1:20; initial(A,B,C,D,x0,t) title(‘The Initial Condition Response’) (2) dinitial( )—求离散系统的零输入响应。用法如下： [y,x,t]=diniial(A,B,C,D,x0) [y,x,t]=diniial(A,B,C,D,x0,n) (3) lsim( )—对任意输入的连续系统进行仿真，用法如下： [y,x]=lsim(A,B,C,D,u,t) [y,x]=lsim(A,B,C,D,u,t,x0) [y,x]=lsim(num,den,u,t) 例：已知某系统如下所示

  31. 求输入为正弦波的输出响应。 程序如下： num=[1 6.8 13.85 8.05] den=[1 11.2 46.4 88.4 77.4 25.2] t=0:0.2*pi:2*pi u=sin(t) lsim(num,den,u,t) 3. 时域分析应用实例 例1：典型二阶系统如下所示： 其中ωn为自然频率（无阻尼振荡频率），ξ为阻尼系数，要求绘出当 ωn=4， ξ分别为0.1,0.2,…,1.0,2.0时系统的单位阶跃响应。

  32. 程序如下： (SysStepko.M) wn=4; kosai=[0.1:0.1:1,2]; figure(1) hold on for i=kosai num=wn.^2; den=[1,2*i*wn,wn.^2] step(num,den) end title('The Step Response of Two Order System'); 在过阻尼和临界阻尼响应曲线中，临界阻尼响应具有最短的上升时间，响应速度最快；在欠阻尼（0< ξ<1)响应曲线中，阻尼系数越小，超调量越大，上升时间越短。一般取ξ=0.4~0.8。

  33. 例2：对例1中的典型二阶系统，绘出当ξ=0.6时， ωn取2,4,6,8,10,12时的单位阶跃响应。 程序如下： (SysStepwn.M) w=2:2:12; kosai=0.6 figure(1) hold on for wn=w num=wn.^2; den=[1,2*kosai*wn,wn.^2] step(num,den) end title('The Step Response of Two Order System');

  34. 例3：一伺服系统的方框图如图所示，求d和e的值，使系统的阶跃响应满足：（1）超调量不大于40%，（2）峰值时间为0.8s。例3：一伺服系统的方框图如图所示，求d和e的值，使系统的阶跃响应满足：（1）超调量不大于40%，（2）峰值时间为0.8s。 由控制理论可知，对于二阶系统， 计算超调量和峰值时间的公式如下： R(s) C(s) 1+es 由此得： 伺服系统的传递函数为： 二阶系统的传递函数标准形式为： 比较二式得：

  35. 程序如下： （SysStepde.M) os=40; tmax=0.8; z=log(100/os)/sqrt(pi^2+(log(100/os))^2); wn=pi/(tmax*sqrt(1-z^2)); num=wn^2; den=[1 2*z*wn wn^2]; t=0:0.02:4; c=step(num,den,t); plot(t,c); xlabel('Time-Sec'); ylabel('y(t)'); d=wn^2 e=(2*z*wn-1)/d title('The Step Response of Two Order System'); grid;

  36. 三、控制系统的频域分析 1. 频域分析的一般方法 1）对数频率特性图（波特图） 对数频率特性图包括了对数幅频特性图和相频特性图。画对数频 率特性图时横坐标的单位为十倍频程（lgw）。 MATLAB提供了函数bode( )来绘制波特图，其用法如下： [mag,phase]=bode(num,den,w) 例1：典型二阶系统如下所示： 绘制出当ξ取不同值时的波特图。 取ωn=5， ξ取[0.1:0.2:2]时的波特图。 程序如下：(见Sysbodeko.M)

  37. wn=5; kosai=[0.1:0.2:2]; w=logspace(-1,1,100) %生成对数横坐标，区间为0.1~10 num=[wn.^2]; for ii=kosai den=[1 2*ii*wn wn.^2]; [mag,pha,w1]=bode(num,den,w); subplot(2,1,1); hold on semilogx(w1,mag); subplot(2,1,2); hold on semilogx(w1,pha); end subplot(2,1,1);

  38. grid on title('Bode Plot'); xlabel('Frequency(rad/sec)'); ylabel('Gain dB'); text(5.5,4.5,'0.1'); subplot(2,1,2); grid on xlabel('Frequency(rad/sec)'); ylabel('Phase deg'); text(4,-20,'0.1'); text(2.5,-90,'2.0'); 2） 极坐标图（奈奎斯特图） 对于频率特性函数G(jω)，给出ω从- 到+ 的一系列数值，分别求 出Im(G(jω))和Re(G(jω))。以Re(G(jω))为横坐标， Im(G(jω))为纵坐 标绘制成极坐标频率特性图。

  39. MATLAB提供了函数nyquist来绘制系统的奈奎斯特图。 [re,im,w]=nyquist(num,den,w) [re,im,w]=nyquist(A,B,C,D) 例1：已知系统的传递函数如下所示： 求当K分别取1300和5200时，系统的极坐标频率特性图。 程序如下：(Sysnyquist.M) k1=1300; k2=5200; w=8:1:80; num1=[k1];num2=[k2]; den=[1 52 100 0]; %subplot(2,1,1); figure(1)

  40. nyquist(num1,den,w); axis([-1.0,0,-0.04,0.04]); grid; %subplot(2,1,2); figure(2) nyquist(num2,den,w) grid; 3）频率响应 MATLAB提供了频率响应函数freqs( ),其用法如下： y=freqs(num,den,w) 例1：系统的闭环函数如下所示： 要求画出系统的幅频特性。 程序如下：

  41. %Sysfreqs.M num=4; den=[1 2 4]; w=0:0.01:3 g=freqs(num,den,w); mag=abs(g); plot(w,mag); xlabel('Frequency -rad/s'); ylabel('Magnitude'); grid; axis([0 3 0.5 1.2]) title('幅频特性')

  42. 2. 常用频域分析函数 (见下表：)

  43. 频域分析应用实例 • 例1: 已知某开环系统如下所示： 要求（1）绘制系统的奈奎斯特曲线，判断闭环系统的稳定性，求 出系统的阶跃响应。（2）给系统增加一个开环极点p=2，求此时 的奈奎斯特曲线，判断此时闭环系统的稳定性，并求出系统的阶 跃响应。 程序如下： （SysExp1.M) %SysExp1.M k=26; z=[];p1=[1 -6];p2=[1 -6 2]; [num1,den1]=zp2tf(z,p1,k); [num2,den2]=zp2tf(z,p2,k); subplot(2,2,1); nyquist(num1,den1);

  44. subplot(2,2,2) hold on [numc1,denc1]=cloop(num1,den1); impulse(numc1,denc1) step(numc1,denc1) subplot(2,2,3) nyquist(num2,den2) title('Nyquist Diagrams'); subplot(2,2,4) hold on t=0:0.1:3 [numc2,denc2]=cloop(num2,den2); impulse(numc2,denc2) step(numc2,denc2,t); axis([0,3,-100,20]);

  45. 例2：线性时不变系统如下所示： 要求绘制系统的波特图和奈奎斯特图，判断系统的稳定性，如果系统稳定，求出系统稳定裕度，并绘制出系统的单位脉冲响应以验证判断绪论。 程序如下： （SysExp2.M） %SysExp2.M a=[-0.6 -1.044 0 0;1.044 0 0 0;0 0.9578 -0.7 -0.3162;0 0 0.3162 0]; b=[1;0;0;0]; c=[0 0 0 0.3162];d=0; w=logspace(-1,1,100); %生成对数横坐标，区间为0.1~10

  46. figure(1);bode(a,b,c,d); figure(2); nyquist(a,b,c,d); figure(3); [ac,bc,cc,dc]=cloop(a,b,c,d); margin(a,b,c,d); figure(4); impulse(ac,bc,cc,dc);

  47. 例3：一个系统如下所示： 要求（1）找出系统的主导极点；（2）求出系统的低阶模型；（3）比较原系统与低阶模型系统的阶跃响应和频率响应。 deno=[1 36 205 750]; r=roots(deno) 求得：r=-30.0000 -3.0000+4.0000i -3.0000-4.0000i 因而传递函数可以写成如下形式： 因为极点s=-30远离原点。其对系统影响很小，可以忽略。因而系统的主导极点应该为-3±4i，于是系统的近似传递函数如下所示：

  48. 对两个系统的性能进行比较。程序如下： %SysExp3.M num1=25;den1=[1 6 25]; %近似的二阶系统 num2=750;den2=[1 36 205 750]; %原三阶系统 figure(1); hold on; bode(num1,den1,'r'); bode(num2,den2,'g'); %比较频率响应 figure(2); hold on step(num1,den1,'r'); step(num2,den2,'g'); %比较单位阶跃响应

  49. 四、根轨迹分析方法 MATLAB绘制轨迹的函数为rlocus(num,den,K)，num和den分别是系 统开环传递函数的分子和分母多项式系统，K为开环增益，K的范围 可以指定，若不给定K的范围，则隐含K从0变至+∞，该函数可以精 确地绘制出根轨迹。常用的根轨迹函数如下表：

More Related