第2章 应用MATLAB建立控制系统模型

第2章应用MATLAB建立控制系统模型MATLAB是国际控制界目前使用最广的工具软件，几乎所有的控制理论与应用分支中都 有MATLAB工具箱。本节结合前面所学自控理论的基本内容，采用控制系统工具箱(Control Systems Toolbox)和仿真环境(Simulink),学习 MATLAB 的应用。一、MATLAB在控制系统微分方程中的应用111物理系统的微分方程利用机械学、电学、流体力学和热力学等的物理规律，我们可以得到物理系统的动态 方程。它们通常用常系数线性微分方程来描述。1.1.2数值解通过拉普拉斯变换和反变换，可得到线性时不变方程的解析解，也可用状态转移矩阵 (t)求解。这些分析方法通常只限于常系数的线性微分方程。解析解是精确的，然而通常寻 找解析解是困难的，甚至是不可能的。而数值分析方法直接在时域里求解微分方程，不仅适 用于线性时不变方程，也适用于非线性以及时变微分方程。MATLAB提供了两个求微分方程数值解的函数，它们采用龙格-库塔(Runge-kutta)法。 Ode23和ode45分别表示采用2阶和4阶龙格-库塔公式，后者具有更高的精度。n阶微分方程必须化为n个首1的一阶微分方程组，且放入M-文件中，以便返回方程 状态变量的导数，下面的例子介绍这些函数的用法。例1.1对图1-1的机械系统，已知三个量一一拉力、摩擦力、以及弹簧力都影响质量 M的加速度。型-1利用牛顿运动定理，建立系统的力平衡方程式d2xdxM + B + Kx 二 f(t)d t dt2人dx令 x 二 x,x 二 ，有12 dtdx1 = xdt 2设质量M=1kg，摩擦系数B=5N/m/sec，弹簧常数K=25N/m。在t=0时刻，施加25N的 拉力。上述方程及已知量在M-文件mechsys.m中定义如下：function xdot=mechsys(t, x);F=25;M=1;B=5;K=25;xdot=x(2);1/M*(F-B*x(2)-K*x(1);F面的M-文件使用ode23对系统在零初始条件下进行仿真:t0=0;tfinal=3;x0=0,0; tol=0.001;trace=0;%时间间隔03秒%零初始条件%精度%如果非零，则打印出每一步的计算值t, x=ode23(mechsys,t0,tfinal,x0,tol,trace) subplot(211),plot(t, x);title (Time response of mechanical translational system5)xlabel (Time-sec)text (2,1.2, displacement)text (2,.2,veloclty)d=x(:,1);v=x(:,2); subplot(212),plot(d,v);title (velocity versus displacement)xlabel (displacement)ylabel (velocity)00.511.522.53Time-sec3210-1图1-2例1.2 电路图见图1-3。R=1.4Q, L=2H, C=0.32F。初始状态：电感电流为零，电容 电压为0.5Vt = 0时刻接入1伏的电压。求0vtv15秒时i(t)和v(t)的值，并且画出电流与 电容电压的关系曲线。根据基尔霍夫定律RLC 有：subplot(111) 仿真结果如图Ve=1vVo (t)Ri + Ldi + v 二 v dt 0x 二 v ,x 二 i1 0 2X 二 (v - x - Rx )2 L s 12以上方程式在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对系统进行仿真。 t0=0;tfinal=15;* 时间间隔 015 秒x0=0.5,0;%初值条件 x1=0.5,x2=0tol=0.001;%精度trace=0;%若非零值则打印出每一步的计算值t, x=ode23(electsys,t0,tfinal,x0,tol,trace); subplot(211), plot(t, x);title(Time response of a RLC series circuit) xlabel(Time-sec)text(8,1.05,Capacitor voltage),text(8,0.05,current) vc=x(:,1),i=x(:,2);subplot(212),plot(vc, i);title(current versus capacitor voltage) xlabel(Capacitor Voltage) ylabel(current)subplot(111)仿真结果见图1-4。Time-secCapacitor Voltage 图1-4113非线性系统在变量的一定范围内大多数物理系统是线性的。然而，随着变量范围的延伸，所有的系统最终都是非线性的。对于非线性系统，叠加原理不再成立Ode23和ode45可以求解非线 性微分方程，见例1-3。mg图1-5例1.3图1-5的单摆，一长L米的无重量绳悬挂重为W=mg千克 的物体，阻尼系数为B千克/米/秒。通常，该系统用线性微分方程近似描 述，实际上它是非线性的。如果8为绳的倾角，物体的速度将为L0，那么离心力为F 二一Wsin 0 BL0根据牛顿定理有F 二 mL 0联结上两式，得mL 0 + BL 0 + Wsin 0 二 0令 x =0,x =0 得1 2. B W .x =x 一 sinx2 m 2 mL 1%状态导数上面方程式在M-文件pendulum.m中定义如下: function xdot=pendulum(t, x);W=2; L=.6; B=0.02; g=9.81; m=W/g; xdot=x(2); -B/m*x(2)-W/(m*L)*sin(x(1);%时间间隔05秒%初始条件X=1,x2=0%精度，不打印每一步的计算值下面的M-文件对系统进行仿真。t0=0; tfinal=5;x0=1,0;tol=0.0001; trace=0;t, x=ode23(pendulum,t0,tfinal,x0,tol,trace);subplot(211),plot(t, x)title(Time response of a rigid pendulum),xlabel(Time-sec)text(3.2,3.5,Velocity),text(3.2,1.2,Angle-Rad) th=x(:,1);w=x(:,2);subplot(212),plot(th, w)title(Phase plane plot of pendulum) xlabel(Position-Rad),ylabel(Angular velocity) subplot(111)仿真结果见图1-6。-44-1-0.8-0.6-0.4-0.200.20.40.60.81Position-Rad00.511.522.533.544.55Time-secPhase plane plot of pendulum2 0 2-y m fie LrKuan-4图1-6二、用MATLAB建立控制系统模型MATLAB是基于数学模型分析控制系统性能的功能强大的软件工具。由于微分方程与传 递函数的模型参数完全一致，在MATLAB中只需要用传递函数描述，这里介绍用MATLAB及 Simulink建立不同形式的传递函数模型，进行不同形式传递函数之间的转换。1 有理函数模型线性系统的传递函数模型可一般地表示为:b sm + b sm-1 hf b s + b(1)G (s)二 T2mm+1sn + a sn-i ff a s + a1n-1n将系统的分子和分母多项式的系数按降幕的方式以向量的形式输入给两个变量num和 den，就可以轻易地将传递函数模型输入到MATLAB环境中。命令格式为：(2)num = b , b，,b , b ；12m m+1den = 1, a , a，,a , a ；1 2n-1n在MATLAB控制系统工具箱中，定义了 tf()函数，它可由传递函数分子分母给出的变 量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。该函数的调用格式为：例2.1G=tf(num, den);(4)一个简单的传递函数模型：s + 5G (s)二s 4 + 2 s 3 + 3s 2 + 4s + 5可以由下面的命令输入到MATLAB工作空间中去。 num=1, 5;den=1, 2, 3, 4, 5;G二tf(num, den)运行结果：Transfer function:s 4 + 2s 3 + 3s 2 + 4s + 5这时对象G可以用来描述给定的传递函数模型，作为其它函数调用的变量。 例2.2 一个稍微复杂一些的传递函数模型：G (s) =6( s + 5)(s 2 + 3s + 1)2( s + 6)该传递函数模型可以通过下面的语句输入到MATLAB工作空间。 num=6*1, 5;den=conv(conv(1, 3, 1, 1, 3, 1), 1, 6); tf( num,den)运行结果Transfer function:6 s + 30s 5 + 12 s 4 + 47 s 3 + 72 s 2 + 37 s + 6其中conv ()函数(标准的MATLAB函数)用来计算两个向量的卷积，多项式乘法当然也 可以用这个函数来计算。该函数允许任意地多层嵌套，从而表示复杂的计算。2. 零极点模型线性系统的传递函数还可以写成极点的形式:(s + z )(s + z )(s + z )12m(s + p )(s + p )(s + p )12n将系统增益、零点和极点以向量的形式输入给三个变量KGain、Z和P,就可以将系统的零 极点模型输入到MATLAB工作空间中，命令格式为：KGain = K;Z = z ; Z Z ;(7)12mP 二P ； P ； ； P ；(8)12n在MATLAB控制工具箱中，定义了 zpk()函数，由它可通过以上三个MATLAB变量构造出 零极点对象，用于简单地表述零极点模型。该函数的调用格式为：G=zpk(Z,P,KGain)(9)例2.3某系统的零极点模型为：G( _ 6(s +1.9294)(s + 0.0353 土 0.9287 j)s _ (s + 0.9567 1.2272j)(s 0.0433 0.6412j)该模型可以由下面的语句输入到MATLAB工作空间中。 KGain=6; z=-1.9294;-0.0353+0.9287j;-0.0353-0.9287j;p=-0.9567+1.2272j,-0.9567-1.2272j,0.0433+0.6412j,0.0433-0.6412j; G=zpk(z， p， KGain)运行结果：Zero/pole/gain:6 (s+1.929) (s2 + 0.0706s + 0.8637)(s2 - 0.0866s + 0.413) (s2 + 1.913s + 2.421)注意：对于单变量系统，其零极点均是用列向量来表示的，故Z、P向量中各项均用分 号(；)隔开。3. 反馈系统结构图模型设反馈系统结构图如图5所示。控制系统工具箱中提供 了 feedback()函数,用来求 取反馈连接下总的系统模图5反馈系统结构图型，该函数调用格式如下：G=feedback(Gl,G2,sign);(10)其中变量sign用来表示正反馈或负反馈结构，若sign=-l表示负反馈系统的模型，若 省略sign变量，则仍将表示负反馈结构。G1和G2分别表示前向模型和反馈模型的LTI(线 性时不变)对象。例2.4若反馈系统图5中的两个传递函数分别为：1(S + 1)21S+1则反馈系统的传递函数可由下列的MATLAB命令得出 G1=tf(l,l,2,l);G2=tf(l,l,l);G=feedback(G1,G2)运行结果：Transfer function:s 3 + 3 s 2 + 3 s + 2 若采用正反馈连接结构输入命令 G=feedback(Gl,G2,l) 则得出如下结果：Transfer function:s + 1s 3 + 3 s 2 + 3 ss 3 + 7 s 2 + 24s + 241 s _ s4 +10s3 + 35s2 + 50s + 24G2(s)二10 s + 510.01s +1例2.5若反馈系统为更复杂的结构如图6所示。其中则闭环系统的传递函数可以由下面的MATLAB命令得出: G1=tf(l,7,24,24,l,10,35,50,24);G2=tf(10,5,l,0);H=tf(l,0.01,l);G_a=feedback(G1*G2,H)得到结果：Transfer function:0.1 s5 + 10.75 s4 + 77.75 s3 + 278.6 s2 + 361.2 s + 1200.01 s 6 + 1.1 s 5 + 20.35 s 4 + 110.5 s 3 + 325.2 s 2 + 384 s + 1204. 有理分式模型与零极点模型的转换有了传递函数的有理分式模型之后，求取零极点模型就不是一件困难的事情了。在控制 系统工具箱中，可以由zpk()函数立即将给定的LTI对象G转换成等效的零极点对象G1。该 函数的调用格式为：G1=zpk(G)(11)例2.6给定系统传递函数为：6.8s 2 + 61.2 s + 95.2 G (s)=s 4 + 7.5s 3 + 22s 2 + 19.5s对应的零极点格式可由下面的命令得出 num二6.8, 61.2, 95.2; den=1, 7.5, 22, 19.5, 0; G=tf( num,den);G1=zpk(G)显示结果：Zero/pole/gain:6.8 (s+7) (s+2)s (s+1.5) (s2 + 6s + 13)可见，在系统的零极点模型中若出现复数值，则在显示时将以二阶因子的形式表示相应 的共轭复数对。同样，对于给定的零极点模型，也可以直接由MATLAB语句立即得出等效传递函数模型。 调用格式为：G1=tf(G)(12)例2.7给定零极点模型：G (s) = 6.8(s + 2)( s + 7)s(s + 3 土 j2)(s +1.5)可以用下面的MATLAB命令立即得出其等效的传递函数模型。输入程序的过程中要注意 大小写。 Z=-2,-7;P=0,-3-2j,-3+2j,-1.5;K=6.8;G=zpk(Z,P,K);G1=tf(G)结果显示：Transfer function:6.8 s2 + 61.2 s + 95.2s 4 + 7.5 s 3 + 22 s 2 + 19.5 s5. Simulink建模方法在一些实际应用中，如果系统的结构过于复杂，不适合用前面介绍的方法建模。在这种 情况下，功能完善的Simulink程序可以用来建立新的数学模型。Simulink是由Math Works 软件公司1990年为MATLAB提供的新的控制系统模型图形输入仿真工具。它具有两个显著的 功能：Simul(仿真)与Link(连接)，亦即可以利用鼠标在模型窗口上“画”出所需的控制系 统模型。然后利用SIMULINK提供的功能来对系统进行仿真或线性化分析。与MATLAB中逐行 输入命令相比，这样输入更容易，分析更直观。下面简单介绍SIMULINK建立系统模型的基本步骤:(1)SIMULINK的启动：在MATLAB命令窗口的工具栏中单击按钮门或者在命令提示符下键入simulink命令，回车后即可启动Simulink程序。启动后软件自动打开Simullink 模型库窗口，如图7所示。这一模型库中含有许多子模型库，如Sources(输入源模块库)、 Sinks(输出显示模块库)、Nonlinear(非线性环节)等。若想建立一个控制系统结构框图，则 应该选择File| New菜单中的Model选项，或选择工具栏上new Model 按钮，打开一个 空白的模型编辑窗口如图8所示。(2) 画出系统的各个模块：打开相应的子模块库，选择所需要的元素，用鼠标左键点中后拖图7 simulink模型库HQE3Ilii L住i 争1 诞Xctila HaljIjiwdy图8模型编辑窗口到模型编辑窗口的合适位置。(3) 给出各个模块参数：由于选中的各个模块只包含默认的模型参数，如默认的传递函数模 型为1/(s+1)的简单格式，必须通过修改得到实际的模块参数。要修改模块的参数，可以用鼠标双击该模块图标，则会出现一个相应对话框，提示用户修改模块参数。（4）画出连接线：当所有的模块都画出来之后，可以再画出模块间所需要的连线，构成完整 的系统。模块间连线的画法很简单，只需要用鼠标点按起始模块的输出端（三角符号）, 再拖动鼠标，到终止模块的输入端释放鼠标键，系统会自动地在两个模块间画出带箭头 的连线。若需要从连线中引出节点，可在鼠标点击起始节点时按住Ctrl键，再将鼠标 拖动到目的模块。（5）指定输入和输出端子：在Simulink下允许有两类输入输出信号，第一类是仿真信号， 可从source（输入源模块库）图标中取出相应的输入信号端子，从Sink（输出显示模块库） 图标中取出相应输出端子即可。第二类是要提取系统线性模型，则需打开Connection（连 接模块库）图标，从中选取相应的输入输出端子。例2.8典型二阶系统的结构图如图9所示。用SIMULINK对系统进行仿真分析。图9典型二阶系统结构图按前面步骤，启动simulink并打开一个空白的模型编辑窗口。（1）画出所需模块，并给出正确的参数： 在sources子模块库中选中阶跃输入（step）图标，将其拖入编辑窗口，并用鼠标 左键双击该图标，打开参数设定的对话框，将参数step time（阶跃时刻）设为0。 在Math（数学）子模块库中选中加法器（sum）图标，拖到编辑窗口中，并双击该图 标将参数List of signs （符号列表）设为|+-（表示输入为正，反馈为负）。 在continuous（连续）子模块库中、选积分器（Integrator）和传递函数（Transfer Fcn）图标拖到编辑窗口中，并将传递函数分子（Numerator）改为900，分母（Denominator） 改为1, 9。 在sinks（输出）子模块库中选择scope（示波器）和0ut1（输出端口模块）图标并将之 拖到编辑窗口中。（3）将画出的所有模块按图9用鼠标连接起来，构成一个原系统的框图描述如图10所示。（4）选择仿真算法和仿真控制参数，启动仿真过程。 在编辑窗口中点击Simulation|Simulationparameters菜单，会出现一个参数对话框， 在solver模板中设置响应的仿真范围StartTime（开始时间）和StopTime（终止时间）， 仿真步长范围Maxinum step size（最大步长）和Mininum step size（最小步长）。对于 本例，StopTime可设置为2。最后点击Simulation|Start菜单或点击相应的热键启动 仿真。双击示波器，在弹出的图形上会“实时地”显示出仿真结果。输出结果如图11 所示。在命令窗口中键入whos命令，会发现工作空间中增加了两个变量一一tout和yout,这 是因为Simulink中的Outl模块自动将结果写到了 MATLAB的工作空间中。利用MATLAB命 令plot(tout,yout),可将结果绘制出来，如图12所示。比较11和12,可以发现这两种输图10二阶系统的simulink实现 出结果是完全一致的。图11仿真结果示波器显示图12 MATLAB命令得出的系统响应曲线三、MATLAB在控制系统模型转换中的应用3.1传递函数向状态空间描述的转换下面一组函数允许控制系统模型的不同表示形式之间可以互相转换。(1) 状态空间模型到传递函数模型的转换。num, den=ss2tf (a, b, c, d, u)(2) 状态空间模型到零-极点增益模型的转换。z, p, k= sstzp (a, b, c, d, u)(3) 传递函数模型到状态空间模型的转换。a, b, c, d = tf2ss (num, den)(4) 传递函数模型到零-极点增益模型的转换。z, p, k = tf2zp (num, den)(5) 零一极点增益模型到状态空间模型的转换。a, b, c, d= zp2ss (z, p， k)(6) 零-极点增益模型到传递函数模型的转换。num, den= zp2tf(z, p, k)(7) 传递函数模型到部分分式模型的转换。z, p, k = residue (num, den)(8) 部分分式模型到传递函数模型的转换。num, den= residue (z, p, k)(9) 连续系统模型到离散系统模型的转换。ad, bd=c2d (a, b, Ts)(10) 离散系统模型到连续系统模型的转换。a, b=d2c (ab, bd, Ts)控制系统工具箱包含一组模型转换的函数A,B,C,D=tf2ss(num, den)将传递函数转换为 状态空间描述。例3.1求下面传递函数的状态空间描述S2+ 7S +2 S3+9S2+26S+24 num=1 7 2; den=1 9 2624; A, B, C, D=tf2ss(num, den)状态方程各矩阵如下：- 9 - 26 - 24_A =1000 1 0丁B =0D=003.2状态空间描述向传递函数的转换已知状态方程和输出方程x = Ax + Buy=Cx+Du采用拉普拉斯变换有Y(s)=C(sI-A)-1Bu(s)+Du(s)则G(s) = C(sI- A) -1 B+ DU(s)函数ss2tf(A,B,C,D,i)是将状态空间描述转换为对第一个输入的传递函数。 num,den=ss2tf(A,B,C,D,i)是将状态空间描述化为分子、分母多项式形式的传递函数。 z,p,k=ss2zp(A,B,C,D,i)将状态空间描述化为零极点形式表示的传递函数。例3.3 一个系统的状态空间描述如下一 x 一01 0x1011x=0 0 1x+022x-1 -2 -3x01-3111-311y=100x求传递函数G(s)=Y(s)/U(s)A=010; 0 001; -1-2 -3;B=10; 0; 0;C=10 0;D=0;num,den=ss2tf(A,B,C,D,1)z,p,k=ss2zp(A,B,C,D,1) 其中，ss2tf(A,B,C,D,1)中“1”表示对第一个输入。传递函数的分子、分母多项式系数如下：num=010.000030.000020.0000den=1.00003.00002.00001.0000传递函数的零、极点如下：z=-1-2p=-0.3376+0.5623i-0.3376-0.5623i-2.3247k=10因而传递函数为、10(s2 + 3s + 2)G(s)二_ s3 + 3s2 + 2s +110(s + l)(s+ 2)-(s+ 0.3376 - 0.5623i)(s + 0.3376 + 0.5623i)(s + 2.3247)3.3由方框图求状态空间描述和传递函数控制系统工具箱中提供了函数A,B,C,D=connect(a, b, c, q, iu, iy)。将方框图描述转换成 状态空间描述和传递函数。其中q矩阵规定了各框之间的连接关系。其每一行的第一个元素 是框号，其余的元素依次是于该框连接的框号，iu,iy分别表示输入，输出施加的框号。例34将图1-7由框图表示的系统转换成状态空间描述和传递函数。n1=1; d=l; n2=0.5; d2=1; n3=4; d3=1 4;n4=1; d4=1 2; n5=1; d5=1 3; n6=2; d6=1;n7=5; d7=1; n8=1; d8=1;nblocks=8; blkbuildq=10000%q矩阵表示框图的结构。21-6-7-8如第2个框于第1个框按320001的关系连接，于第6.7.843000个框按-1关系连接，依次类推。54000630007400085000;%输入施加于第1个框上%由第8个框输出iu=1;iy=8；A, B, C, D=connect(a, b, c, d, q, iu, iy)num, den=ss2tf(A,B,C,D,1)%转换成传递函数结果为A=-8.0-2.5-0.50.4-2.0001.0-3.0B=0.500C=0 0 1D=0num=0 0 0 2den=1.013.056.080.0C(s) _2R(s)S3 + 13s2 + 56s+ 80