控制系统仿真应用

MATLAB在控制系统仿真中的应用从本章开始，正式进入到有关MATLAB在控制系统仿真的领域中去。我们知道，对控制系统进行仿真，首先要建立系统的数学模型，这是计算机仿真的基础。我们下面先介绍控制系统数学模型的描述方法。第一节 控制系统数学模型的基本描述方法在控制系统仿真中，主要用4种形式的数学模型：传递函数.、零极点模型、结构图形式和状态方程模型。这些模型之间存在着内在的等效关系。在不同的场合下可能使用的模型形式要求不同，需要了解模型之间的转换方法。这一节主要介绍它们的MATLAB实现。一控制系统的传递函数描述：（在MATLAB中称tf模型）对系统的微分方程在零初始条件下做拉氏变换，则可得系统的传递函数(SISO系统)：G(s)=对线性时不变（线性定常）系统（LTI）来说，a、b均为常数a10。前面讲过多项式的表示方法。这里分子分母都为多项式，可将分子分母分别表示出来。即用分子分母的系数构成两个向量，唯一的确定出来： num=b,b,bden=a,a,a注意：构成分子，分母向量按降幂排列的顺序。这只是一个简单的特例，很多时候，传递函数的分子、分母均为多项式相乘的形式，如：不能直接写出，可借助多项式乘法运算函数conv( )来处理，以便获得分子、分母多项式向量。conv( )函数的调用方式为：c=conv(a,b)a,b各表示一个多项式，c表示a和b的乘积多项式，此函数允许嵌套使用：c=conv(a,conv(b,e)，则上例的G(s)可用下面的语句来输入：num=4*conv(1,2,conv(1, 6, 6,1, 6, 6);den=conv(1,0,conv(1,1,conv(1,1,conv(1,1,1,3,2,5);还可以进一步地编写一个convs( )函数来一次性地求出若干个（十个）多项式的连乘积：function a=convs(a1,a2,a3,a4,a5,a6,a7,a,8,a9,a10)a=a1;for i=2:nargin %在MATLAB中它也是一个固定变量，表示在函数调用时实际输入变量的个数。 eval(a=conv(a,a int2str(i) );) %eval( )函数执行MATLAB语句构成的字符串。end这里定义的函数，可以一次执行十个多项式的乘积。如上例可以写成：convs(1,0,1,1,1,1,1,1,1,3,2,5)注意：写完函数后一定要按函数名来保存该文件。相应地，对离散时间系统的动态模型，用脉冲传递函数描述，其输入方法与此类似，不多举例。如果是MIMO系统，则用传递矩阵描述，如： 则可表示为：num= 1,1;1den=1 2 2 ;1 1 二控制系统的零极点模型：（zpk模型）是传递函数的另一种表现形式，对SISO系统：G(s)=将零点、极点及K值输入即可建立零极点模型。 z=-z,-z ,-z p=-p,-p ,-p k=k对于给出的传函来说，分别对分子分母作因式分解，则可以得出系统的零极点模型。这可以通过求出分子，分母多项式的根来实现。MATLAB中提供了多项式求根的函数，roots ( )调用格式： z=roots(a) 其中：z 各个根所构成的向量 a 多项式系数向量系统增益k即为原传递函数分子的最高项系数与分母最高项系数的比值。分别求出分子分母多项式的根，即可得到系统的零，极点模型。对于多输入多输出系统，应分别对每个输入求出系统的零极点模型最后才可以获得整个系统的零极点模型（为矩阵形式）。两种模型之间的转换函数：z ,p ,k=tf2zp(num , den);num , den =zp2tf(z ,p ,k)三. 控制系统状态方程模型：( ss 模型)LTI系统的状态方程： = Ax + Bu Y = Cx + Du只要将A，B，C，D几个矩阵输入进去即可。对于离散系统来说，也与上面类似。MATLAB还提供了由系统状态方程转换为系统零极点模型及传递函数模型的函数。num , den=ss2tf(A ,B ,C ,D , iu ) %iu表示输入的序号（对多输入系统）z, p ,k=ss2zp(A ,B ,C ,D , iu ) %z ,p ,k 表示对第iu个输入信号的传递函数的零极点。A , B ,C ,D= tf2ss (num , den)A , B ,C ,D=zp2ss(z ,p ,k)四 在MATLAB控制系统工具箱中的LTI对象：为了避免对一个系统采用多个分离变量进行描述，新版本的控制系统工具箱，将LTI系统的各种描述封装成一个对象，即用一个变量来描述。在控制系统工具箱中，有以上讲述的三种对象，即ss 对象，tf对象和zpk对象。每种系统模型的生成和模型间的转换均可以通过一个函数来实现。下面介绍这些函数：1 dss( )函数：生成系统的状态空间模型。格式： sys= dss (a ,b ,c ,d ,e ) 该调用生成连续系统的状态空间模型：= = E为非奇异阵（E阵奇异，为奇异系统）。sys = (a ,b , c, d ,e ,Ts)生成离散系统的状态空间模型。= Ts为采样周期2filt()函数：生成DSP（数字信号处理）形式的离散传递函数： 格式：sys=filt(num ,den) - 生成离散传递函数模型 sys=filt(num ,den ，Ts)-定义模型的采样周期注：DSP（数字信号处理）形式的离散传递函数如下面的形式：3ss()函数：生成状态空间模型，或者将传递函数及零极点模型转换成状态空间模型。 格式：1）sys=ss(a,b,c,d)-生成连续系统的状态空间模型。形式为： y =Cx + Du 2）sys=ss(a,b,c,d,Ts)-生成离散状态空间模型x(n+1)=Ax(n)+Bu(n)y(n)=Cx(n)+Du(n) 3）sys_ss=ss(sys)-将任意的LTI对象sys转换成状态空间模型。4tf()函数：生成传递函数模型，或将零极点模型及状态空间模型转换成传递函数模型。 格式为：sys=tf(num,den)-生成连续时间系统传函模型: G(s)=sys=tf(num ,den, Ts)-生成离散时间系统传函。tfsys=tf(sys)-将任意的LTI对象转换成传递函数模型。5zpk()函数：生成零极点模型或者将其他模型转化成零极点模型。 格式：sys=zpk(z ,p ,k) - 连续系统的零极点增益模型。sys=zpk(z ,p ,k,Ts)- 离散时间系统的零极点增益模型。zsys=zpk(sys) - 将任意LTI对象转换成零极点增益模型。五. 控制系统的结构图描述及转换： 由于复杂的结构图建模可以用SIMUKINK。在这里我们仅介绍简单的结构图描述，也可称为系统建模。 典型连接：1串联连接：21 =子系统1： =Ax+B u 子系统2： =Ax+Bu y= cx+Du y= Cx+Du串联后系统模型为： =+ y= +D2D1u1 2并联连接： 1 +2并联后系统模型： =+u y=+(+)u3负反馈连接： 1 -2 为了简化,设，则反馈后：=+=+(r-) =+r-(+) =+r- = 所以，状态方程为： 在MATLAB工具箱中，提供了子系统的连接处理函数：1） series()函数：系统串联实现格式：sys=series(,)A,B,C,D=series(A,B,C,D, A,B,C,D)2） parallel()函数：系统并联实现。格式：sys=parallel(,)3） feedback()函数：系统反馈连接格式：sys=feedback(,) sys=feedback(,sign) sign定义反馈形式：正反馈，sign=+1；负反馈，sign=-1。 4） append函数：多个LTI系统的组合格式:sys=append(,)对N个子系统组合： 设各个子系统的模型为：(s), (s), (s)则合成后子系统为分块对角阵： 如子系统为状态空间模型，则组合后： 第二节 控制系统的稳定性分析在计算机仿真中，控制系统稳定性的分析方法有： 求闭环特征方程的根； 化为零极点模型，看极点是否在s右半平面； 对状态空间形式（闭环），求A阵的特征值eig(A)； 用Lyapunov方程求解。 在自控原理中，我们均采用间接方法劳斯判据，奈氏判据等，但由于在MATLAB中很容易地求解多项式方程。因此，我们可以直接求出特征方程的根，可用多项式求根函数roots()。 例1：判断如下系统的稳定性：G(s)= -可编程如下：numg=1; deng=1 1 2 23; numf=1; denf=1; (num,den)= feedback(numg,deng,numf,denf,-1); roots(den)例2：判断状态空间形式描述的系统的稳定性，两个输入。A=0 1 1；-6 11 6；-6 11 5；B=0 0 1；C=1 0 0； D=0；,=ss2zp(A,B,C,D,1),=ss2zp(A,B,C,D,2)ii=find(real()0),=length(ii);%find（ ）：找出数组中满足实部值大于0的所有元素的下标，并将结果返回到数组中。ii1= find(real()0), =length(ii1);if (+0),disp(system is Unstable)else disp(system is Stable);end一般当系统中含有非线性环节时，即不能用上述方法，一般需求解李雅普诺夫方程：P+PA=-Q 求出P为正定对称矩阵，则系统是稳定的。李雅普诺夫（Lyapunov）方程求解函数Lyap()（采用Schur方法编写的），其调用格式：x=lyap(A,B,C)该函数可求解：AX+XB=-C方程如要求解：P+PA=-Q 则可写成：p=lyap(A,Q)注：判断稳定性的方法还可确定是否为最小相位系统（没有右半平面零极点的系统）。第三节 控制系统的根轨迹在控制系统工具箱中提供了系统根轨迹绘制及分析函数，还有一个GUI分析工具。下面介绍一下这些根轨迹的绘制及分析函数：1 rlocfind():计算给定根的根轨迹增益 格式：k,poles=rlocfind(sys) k,poles= rlocfind(sys,p) 在LTI对象的根轨迹图中显示出十字光标，当用户选择其中一点时，其相应的增益由k记录，与增益相关的所有极点记录poles中。若要使用该函数，必须首先在当前窗口上绘制系统的根轨迹。 定义要得到增益的根矢量P，即事先给出极点。除了显示出该根对应的增益以外，还显示出该增益对应的其它根。2 rlocus() 函数：功能为求系统根轨迹。 格式：rlocus(sys)- 计算SISO开环LTI对象的根轨迹，增益自动选取。 rlocus(sys,k)-显式设置增益。 r,k= rlocus(sys)- 返回系统的增益k（向量）和闭环极点r（向量），即对应于增益的闭环极点为r(i)。3 sgrid()函数：绘制连续时间系统根轨迹和零极点图中的阻尼系数和自然频率栅格线。（在根轨迹图中加和网格）格式：sgrid-在连续系统根轨迹或零极点图中绘制栅格线。栅格线由等阻尼系数和等自然频率线构成，阻尼系数步长为0.1范围从0到1，自然频率步长为1弧度/秒，范围从0到10。在绘制前当前窗口必须包含根轨迹图。 sgrid(z,)-自己定义想要绘制的阻尼系数向量z和自然频率向量的范围。 注：在画出的根轨迹图上面点击，可以给出对应于所点击点的极点值、对应增益、阻尼比和超调量等。4 zgrid()函数：绘制离散时间系统根轨迹和零极点图中的栅格。格式：zgrid- 在离散系统的根轨迹图或零极点图上绘制出栅格线。 zgrid(z，)-显式定义想要绘制的z向量和向量。第四节 系统的频率响应一频率响应的计算：若系统的传递函数为：G(s)= 则频率特性函数 G（j）=若系统由状态空间模型来描述，则系统的频率特性为：G(j)=CjI+AB+D频率响应的求法：1）对传函形式可用：GWI=polyval(num,sqrt(-1)*w)./polyval(den,sqrt(-1)*w) Polyval（）：多项式求值函数 sqrt(-1)*w：构造j向量2）对状态空间模型形式可用：evalfr()和freqresp()函数来求evalft()格式：frsp=evalfr(sys，f) f：复数频率点，f值以复数给出功能：计算系统在单个频率点的频率响应freqresp()功能：计算系统频率响应 格式：H=freqresp(sys,) H：是三维向量 ：实数频率向量，可指定一些实数频率点如：计算下述系统在频率1，10，100处的频率响应。H(s)=num=0,1;1,-1,1;den=1,1 1;1 2,1;p=tf(num,den)w=1 10 100;H=freqresp(p,w)MATLAB将返回：H(:,:,1)= 对第1个频点 H(:,:,2)= 对第2个频点 H(:,:,3)= 对第3个频点二频率响应曲线的绘制： 提供了3种频率响应曲线的绘制函数：Bode图绘制，Nyquist 曲线及Nichols曲线的绘制，它们可用于SISO或MIMO的连续时间和离散时间系统。1. Bode图绘制：bode()函数调用格式：bode(sys): bode(num,den)bode(A,B,C,D,iu) -可自动地选择一个合适的频率范围。bode(sys，)-给出频率范围，这里频率范围一般由=logspace(a,b,n)给出。logspace(a,b,n)：表示在 到之间的 n个点。这两种格式可直接画出规范化的图形（有网格，有说明）。mag,phase,=bode(sys)或m,p=bode(sys) 这种格式只计算Bode图的幅值向量和相位向量，不画出图形。要在此基础上画图，可用：subplot(211);semilogx(,20*log10(m) %对数幅频曲线subplot(212);semilogx(,p) %对数相频曲线 bode(,)bode(,)这两种格式可在一个图形窗口同时绘多个系统的bode图。2. Nyquist曲线的绘制：nyquist()函数格式：nyquist(sys) nyquist(sys,) nyquist(,)nyquist(,) re,im,=nyquist(sys) re频率响应实部 im虚部3. Nichols(sys)曲线：nichols()函数格式：nichols(sys).nichols(sys,). nichols(,)nichols(,)mag,phase,=nichols(sys,) = nichols(sys)mag：幅值 phase：相位四计算控制系统的稳定裕度：margin( )函数利用该函数可以直接求出系统的幅值裕度和相角裕度。调用格式为：，= margin(num,den) = margin(A,B,C,D) = margin(sys)- 幅值裕度；-相位裕度；-幅值裕度处对应的频率；-相位裕度处对应的频率。 ，=margin(mag ,phase,)-根据幅值、相位数据计算稳定裕度。 margin(sys)-在当前图形窗口中绘制出系统裕度的Bode图。 例：对本节内容的综合使用程序： A=0 1 1；-6 11 6；-6 11 5；B=0 0 1；C=1 0 0； D=0；=logspace(-1,1)；m,p=bode(A,B,C,D,1,)subplot(211),semilogx(,20*log10(m)；subplot(212),semilogx(,p)；X,Y=nyquist(A,B,C,D,1,);plot(X,Y)，=margin(A,B,C,D)第五节 线性系统时间响应分析 系统仿真的实质就是对描述系统的数学模型进行求解，一般的求解方法为微分方程的数值解法（龙格库塔法ode23( ),ode45( )等），这是系统仿真的基本方法。 但我们没有必要直接采用上面的方法，因为工具箱已提供了直接求各种输入下系统响应的函数。 一般控制系统的输入有：阶跃，斜坡，加速度及脉冲输入等。在工具箱主要提供了如下三种函数：1step( )函数-求系统阶跃响应 其调用格式为：1） 2）-定义仿真时间这两种格式直接画出响应曲线 3）y,x=step(sys)-只计算仿真数据，不画图。=step(sys,t) y -各个仿真时刻的输出向量。x -自动选择的状态向量的时间响应数据。4）step(,)同时仿真多个系统。2impulse( )函数：求取系统的脉冲响应格式与step函数相同。3lsim( )函数：求系统的任意输入下的仿真格式：1）lism(sys ,u ,t) u-给定输入构成的列向量，它的元素个数应与t的元素个数相一致。这种形式还可以写成：lsim(num ,den ,u ,t)lsim(A ,B ,C ,D ,iu ,u ,t)其中：iu-输入变量序号2）lsim(sys ,u ,t,)-定义初始状态3）y,x=lsim(sys ,u ,t ,)-不画出图形只算数据4）lsim(,t)-同时仿真多个LTI对象例：求下面系统在阶跃信号为0.011(t)时系统的响应。num=20； den=1 8 36 40 20;t=0:0.1:20;=length(t);u=ones(,1)*0.01;y，x=lsim(num ,den ,u ,t)plot(t,y)第六节 延迟系统的仿真 带有时间延迟的控制系统传函模型 G(s)= =T为延迟时间常数。纯时间滞后环节可以由有理函数来近似，1892年法国数字家Pade曾提出了一种著名的有理近似方法，后人将它命名为Pade近似法。其表达式为：=求这样的近似有理函数，控制系统工具箱提供了一个函数pade( )。 pade( )：计算pade近似的系数。 它的调用格式为： num ,den= pade(T ,n) - 返回近似模型 T-延迟时间常数 n-要求拟和的阶数或：A，B，C，D=pade(T ,n)-返回等效状态方程模型拟和的阶数越高，精度也越高。另，建立延迟系统的数学模型还有两种方法：1） 用函数TF、ZPK或SS建模并同时给模型时间延迟属性“IoDelayMatrix、InputDelay、OutputDelay”赋值。2） 先用函数TF、ZPK或SS建模，再用set命令给模型时间延迟属性赋值。例：已知系统传递函数为用MATLAB生成系统模型程序及运行结果如下：%Create system modelMun=2;Den=1 1;sys=tf(num,den,InputDelay,0.1)Transfer function:exp(-0.1*s)*或：%Create system modelMun=2;Den=1 1;sys=tf(num,den)set(sys, InputDelay,0.1)第七节 状态空间法仿真实现在控制系统工具箱中，有各种状态空间法实现的函数： 1ss2ss( )：功能是进行相似变换如有状态方程：引入非奇异变换矩阵T，使变换后的状态变量(t)=Tx(t) 则变换后系统的状态方程： 其中： = 该变换为相似变换。该函数的格式为： sysT=ss2ss(sys,T)-完成状态空间中LTI对象sys的相似变换,T为变换矩阵2canon( )函数：功能是状态空间的正则实现 该函数计算连续时间或者离散时间LTI对象的正则实现，它支持两种正则形式：模态矩阵(modal matrix)和伴随矩阵(companion matrix)形式。格式：1）csys=canon(sys ,type) type可以为以下字符串：modal-计算模态矩阵正则实现 companion-计算伴随矩阵正则实现sys-如果系统不是以状态空间模型给出的，要首先进行模型转换。 2）csys ,T=canon(sys ,type) T-返回相似变换矩阵T，当给出的LTI对象不是ss对象时T为是空矩阵。 注：模态实现要求系统的A矩阵为可对角化的。3ctrb()函数：功能是计算可控性矩阵 Co=B AB B 格式：Co=ctrb(A,B) Co=ctrb(sys) 这里sys为状态空间对象如果rank(Co)=n 则系统可控4 ctrbf()函数：系统的可控与不可控分解格式：Abar ,Bbar ,Cbar ,T,k=ctrbf(A,B,C) T-相似变换阵 k-是长度为n的一个矢量，其元素为各个块的矩阵Abar ,Bbar ,Cbar-对应于转换后系统的（A，B，C）,通过它们看出可控与不可控部分。 5. obsv() 函数：计算可观性矩阵格式：ob=obsv(A,C) ob=obsv(sys)如果rank(ob)=n ,则系统可观。6. obsvf():系统的可观与不可观分解格式：Abar ,Bbar ,Cbar ,T,k=obsvf(A,B,C)与ctrbf()函数相似。7. minreal():状态空间的最小实现格式：sysr=minreal(sys)-实现LTI对象的最小实现，即删除状态空间模型中的不可观和不可控状态，或者对消零，极点模型中相同的零极点对。sysr具有最小的阶数。如可给零极点模型 =用该函数后可返回：8模型降阶函数： 格式：1）降阶后系统rsys具有与原系统相匹配的直流增益(DC gain) elim-给出需要删除的状态的索引 2）rsys=modred(sys ,elim ,del) 仅仅简单删除要删除的状态，不能保证匹配的直流增益。9SISO系统的极点配置函数：acker()格式：k=acker(a ,b ,p)- 利用Ackermann公式计算反馈增益矩阵k,使u=-kx的单输入系统具有指定的闭环极点p 。 a,b-p=eig(a-bk) 10MIMO极点配置函数：place() 格式：k=place(A,B,p)-p=eig(A-Bk) k ,prec ,message=place(A,B,p) prec-系统闭环实际极点与希望极点的接近程度。 message-如果系统闭环的实际极点偏离希望极点10%以上，将给出警告信息。 place 函数还可用来计算估计器增益：L=place(A,C,p)11.求系统状态估计器函数：estim() 格式：est=estim(sys,L)-生成状态和输出估计器 L- 估计器增益矩阵，由函数形成或由滤波函数生成。 sys- 状态空间模型 est=estim(sys ,L ,sensors ,known)-生成更为普遍的状态和输入估计器 sensors-指定可观输出； known-指定已知输入。12系统调节器生成函数： reg( ) 格式：-生成给定状态估计增益L及状态反馈增益阵k下的系统的调节器或补偿器。sys为状态空间模型. -生成更为普遍的LTI系统的状态和输出调节器 -可测输出； -已知输入；-控制输入最优控制仿真 LQ和LQG： LQ-线性二次型最优控制问题 LQG-线性二次型最优控制问题（对系统中含有系统噪声和量测噪声的系统的最优控制问题）。1 函数：连续系统的LQ调节器设计 格式：= 或 = -状态反馈阵 -方程的解 -闭环系统特征值 Q，R-性能指标 2函数：离散系统的LQ调节器设计 格式：k,s,e=slqr(a,b,Q,R) k,s,e=slqr(a,b,Q,R,N) 3. 函数：系统的LQ调节器设计 格式：k,s,e=lqry(sys,Q,R) k,s,e=lqry(sys ,Q ,R ，N)Q,R-性能指标： 4lqrd()：连续时间系统的离散LQ调节器设计 格式：，s ,e =lqrd(A,B,Q,R,) LQG控制问题可以分解为两个问题求解：1） 利用滤波理论，从状态中得到最优估计。2） 利用估计状态代替，设计满足二次型性能指标的LQ控制律。 58