1、控制系统MATLAB计算与仿真,控制系统MATLAB仿真基础 1 控制系统数学模型 2 系统时域分析 3 系统频域分析 4 系统状态空间分析 5 控制系统的MATLAB仿真及计算 6 控制系统设计,MATLAB概述,Matlab=Matrix Laboratory：矩阵实验室Matlab 语言的首创人是 C. Moler到九十年代初期，在国际上 30 几个数学类科技应用软件中， Matlab 在数值计算方面独占鳌头，而 Mathematica 和 Maple 则分居符号计算软件的前两名。,MATLAB的主要功能,数值计算功能 符号运算功能 绘图功能 编程功能 丰富的工具箱(ToolBox) S
2、imulink动态仿真集成环境,命令窗口,当前工作目录,当前工作空间,输入命令的历史记录,命令 提示符,MATLAB的语言规则,Matlab 区分大小写，它的命令和函数全是小写的一行可以输入几个命令，用分号“;” 或逗号“,”隔开 续行符： “”（三个点），如果语句很长，可用续行符将一个语句写成多行。续行符的前最好留一个空格。注释符： “%”，其后面的内容为注释，对 Matlab 的计算不产生任何影响,MATLAB的变量,Matlab中变量名是以字母开头，可以由字母、数字或下划线组成who，whos，clc，clear系统预定义变量：eps, pi, Inf, NaN, i, j,MATLAB
3、数据的输入,直接输入：同一行中各元素之间用“空格”或“,” (英文状态下) 分开；行与行之间用 “;” 或“回车”分开a=1 2 ; 3, 4 冒号 “:” 运算符：初值 : 步长 : 终值 a=1:5， b=0:pi/4:pi 由向量或小矩阵生成：X=a ; b交互式输入：inputn=input(Please input n: ),MATLAB的基本语句结构,赋值语句 直接赋值语句：赋值变量=赋值表达式 （赋值表达式后无分号“;”: 命令窗直接显示赋值结果赋值表达式后加分号“;”：命令窗不显示结果，变量保留在workspace)函数调用语句：返回变量列表=函数名(输入变量列表),MATLA
4、B的基本语句结构,循环语句 for语句：for 变量=向量，循环体，end；s=0;for i=1:2:100s=s+i;end while语句：while (条件式)，循环体，end；s=0;i=1;while(i=100)s=s+i;end,MATLAB的基本语句结构,跳转语句 if语句：if (条件1)语句组1;elseif(条件2)语句组2;else语句组2;end,MATLAB的基本语句结构,跳转语句 switch语句：case (表达式1)语句组1;case (表达式2)语句组2;, otherwise语句组n;end,MATLAB的M-文件,M-脚本文件M-函数function
5、返回变量列表=函数名 (输入变量列表)function a=dbness(n)a=2*n;end,MATLAB的图形可视化,plot函数:,clc; x=-5:0.02:5; m,n=size(x); %返回矩阵x的维数 for a=1:1:nif x(a)=-1.1y(a)=x(a);elsey(a)=1.1;end end plot(x,y);,1 控制系统的数学描述与建模,控制系统的数学模型在控制系统的研究中有着相当重要的地位，要对系统进行仿真处理，首先应当知道系统的数学模型，然后才可以对系统进行模拟。同样，如果知道了系统的模型，才可以在此基础上设计一个合适的控制器，使得系统响应达到预期
6、的效果，从而符合工程实际的需要。,在线性系统理论中，一般常用的数学模型形式有：传递函数模型（系统的外部模型）、状态方程模型（系统的内部模型）、零极点增益模型和部分分式模型等。这些模型之间都有着内在的联系，可以相互进行转换。,按系统性能分：线性系统和非线性系统；连续系统和离散系统；定常系统和时变系统；确定系统和不确定系统。 1、线性连续系统：用线性微分方程式来描述，如果微分方程的系数为常数，则为定常系统；如果系数随时间而变化，则为时变系统。今后我们所讨论的系统主要以线性定常连续系统为主。 2、线性定常离散系统：离散系统指系统的某处或多处的信号为脉冲序列或数码形式。这类系统用差分方程来描述。 3、
7、非线性系统：系统中有一个元部件的输入输出特性为非线性的系统。,一 系统的分类,微分方程是控制系统模型的基础，一般来讲，利用机械学、电学、力学等物理规律，便可以得到控制系统的动态方程，这些方程对于线性定常连续系统而言是一种常系数的线性微分方程。 如果已知输入量及变量的初始条件，对微分方程进行求解，就可以得到系统输出量的表达式，并由此对系统进行性能分析。 MATLAB提供了ode23、ode45等微分方程的数值解法函数，不仅适用于线性定常系统，也适用于非线性及时变系统。,二 线性定常连续系统的微分方程模型,控制系统MATLAB仿真基础,1.数学模型种类与转换 MATLAB中用四种数学模型表示控制系
8、统：传递函 数模型、零极点增益模型、状态空间模型、动态 结构图（SIMULINK中使用）。 实际应用中，根据要解决问题的需要，往往要进 行不同种类模型之间的转换，MATLAB控制系统工 具箱中为此提供了相应的命令函数。,对线性定常系统，式中s的系数均为常数，且a1不等于零，这时系统在MATLAB中可以方便地由分子和分母系数构成的两个向量唯一地确定出来，这两个向量分别用num和den表示。num=b1,b2,bm,bm+1den=a1,a2,an,an+1注意：它们都是按s的降幂进行排列的。,三 传递函数描述,1、连续系统的传递函数模型 连续系统的传递函数如下：,（1）传递函数模型 MATLAB
9、中用函数命令tf（）来建立控制系统的传 递函数模型。 函数命令的调用格式为： sys=tf（num，den） 函数返回的变量为连续系统的传递函数模型，函数 输入参量num与den分别为系统的分子与分母多项式 系数向量。 sys=tf（num，den，Ts） 函数返回的变量为离散系统的传递函数模型，函数 输入参量num与den分别为系统的分子与分母多项式 系数向量，Ts为采样时间。 MATLAB中还用函数printsys（）来输出系统模型。,零极点模型实际上是传递函数模型的另一种表现形式，其原理是分别对原系统传递函数的分子、分母进行分解因式处理，以获得系统的零点和极点的表示形式。,在MATLAB
10、中零极点增益模型用z,p,K矢量组表示。即： z=z1,z2,zm p=p1,p2,.,pn K=k 函数tf2zp()可以用来求传递函数的零极点和增益。,2、零极点增益模型,K为系统增益，zi为零点，pj为极点,（2）零极点增益模型 MATLAB中用函数命令zpk（）来建立控制系统的零 极点增益模型。 函数命令的调用格式为： sys=zpk（z，p，k） 其中的z、p、k分别代表系统零点、极点、增益向量， 函数返回连续系统零极点模型。 sys=zpk（z，p，k，Ts） 其中的z、p、k分别代表系统零点、极点、增益向量， 函数返回离散系统零极点模型，Ts为采样时间。,控制系统常用到并联系统，
11、这时就要对系统函数进行分解，使其表现为一些基本控制单元的和的形式。 函数r,p,k=residue(b,a)对两个多项式的比进行部分展开，以及把传函分解为微分单元的形式。 向量b和a是按s的降幂排列的多项式系数。部分分式展开后，余数返回到向量r，极点返回到列向量p，常数项返回到k。 b,a=residue(r,p,k)可以将部分分式转化为多项式比p(s)/q(s)。,3、部分分式展开,举例：传递函数描述 1）num=12,24,0,20;den=2 4 6 2 2;2）借助多项式乘法函数conv来处理： num=4*conv(1,2,conv(1,6,6,1,6,6); den=conv(1,
12、0,conv(1,1,conv(1,1,conv(1,1, 1,3,2,5);,num=1,11,30,0; den=1,9,45,87,50;z,p,k=tf2zp(num,den),结果表达式：,零极点增益模型：,num=2,0,9,1; den=1,1,4,4; r,p,k=residue(num,den),结果表达式：,部分分式展开：,（4）状态空间模型 MATLAB中用函数命令ss（）来建立控制系统的零极 点增益模型。 函数命令的调用格式为： sys=ss（a，b，c，d） 其中的a、b、c、d分别代表系统状态矩阵、控制矩阵、 输出矩阵、直接传输矩阵，函数返回连续系统状态空 间模型。
13、 sys=ss（a，b，c，d，Ts） 其中的a、b、c、d分别代表系统状态矩阵、控制矩阵、 输出矩阵、直接传输矩阵，函数返回离散系统状态空 间模型，Ts为采样时间。,模 型 转 换,例 1用Matlab表示传递函数为 的系统解：num=2 3；den=4 3 2 1；sys=tf（num，den） 执行结果： Transfer function:2 s + 3 - 4 s3 + 3 s2 + 2 s + 1,例 2 用Matlab表示传递函数为 的系统。 解：z=-3；p=0 -1 -2；k=2；sys=zpk（z，p，k） 执行结果： Zero/pole/gain:2 (s+3) - s
14、(s+1) (s+2),例 3已知某控制系统的传递函数为 ，求Matlab描述的传递函数模型及零极点增益模型。解：num=0 0 1；den=1 3 2；sys1=tf（num，den）；z p k= tf2zp（num，den）；sys2=zpk（z,p,k）执行结果： Transfer function:1 - s2 + 3 s + 2,Zero/pole/gain:1 - (s+2) (s+1),2. 方框图模型化简在MATLAB中描述系统的模型形式不仅仅拘泥于 数学表达式，还有应用在SIMULINK仿真环境中的动 态方框图形式。只要按照一定的规则画出系统模型图， 然后用 实际系统的数据
15、进行设置，就可以对其实现仿 真。 模型化简包括：环节串联化简；环节并联化简；反馈环节化简。,（1）环节串联化简 多个环节相串联的连接形式是控制系统最基本的组成 结构形式之一。在MATLAB中用函数命令series（）来实现串联化简。 其调用格式为： sys=series（sys1，sys2） 说明： 也可以简单地通过命令sys=sys1sys2实现 。,（2）环节并联化简 多个环节相并联的连接形式也是控制系统最基本的 组成结构形式之一。在MATLAB中用函数命令parallel（）来实现串联化简。 其调用格式为： sys= parallel （sys1，sys2） 说明： 也可以简单地通过命令
16、sys=sys1+sys2实现 。,（3）反馈环节化简 反馈连接结构是控制系统动态方框图动中常见连接 形式，为此MATLAB提供了相应的函数命令实现反 馈化简。在MATLAB中用函数命令feedback（）来实现串联化简。 其调用格式为： sys= feedback （sys1，sys2，sign） 其中sign缺省值为-1，即表示负反馈；若为正反馈，则sign=1。,二 控制系统的时域分析,一个动态系统的性能常用典型输入作用下的响应来描述。响应是指零初始值条件下某种典型的输入函数作用下对象的响应，控制系统常用的输入函数为单位阶跃函数和脉冲激励函数（即冲激函数）。在MATLAB的控制系统工具箱
17、中提供了求取这两种输入下系统响应的函数。,1、时域分析的一般方法,求取系统单位阶跃响应：step() 求取系统的冲激响应：impulse(),1、step()函数的用法,y=step(num,den,t)：其中num和den分别为系统传递函数描述中的分子和分母多项式系数，t为选定的仿真时间向量，一般可以由t=0:step:end等步长地产生出来。该函数返回值y为系统在仿真时刻各个输出所组成的矩阵。,y,x,t=step(A,B,C,D,iu)：其中A,B,C,D为系统的状态空间描述矩阵，iu用来指明输入变量的序号。x为系统返回的状态轨迹。,如果对具体的响应值不感兴趣，而只想绘制系统的阶跃响应曲
18、线，可调用以下的格式： step(num,den)；step(num,den,t)；step(A,B,C,D,iu,t)；step(A,B,C,D,iu)；,线性系统的稳态值可以通过函数dcgain()来求取，其调用格式为：dc=dcgain(num,den)或dc=dcgain(a,b,c,d),y,x,t=step(num,den)：此时时间向量t由系统模型的特性自动生成, 状态变量x返回为空矩阵。,exp4_3.m（传递函数模型），exp4_3_.m（状态空间模型）,2、impulse()函数的用法,求取脉冲激励响应的调用方法与step()函数基本一致。 y=impulse(num,de
19、n,t)；y,x,t=impulse(num,den)；y,x,t=impulse(A,B,C,D,iu,t) impulse(num,den)；impulse(num,den,t) impulse(A,B,C,D,iu)；impulse(A,B,C,D,iu,t),仿真时间t的选择： 对于典型二阶系统根据其响应时间的估算公式 可以确定。 对于高阶系统往往其响应时间很难估计，一般采用试探的方法，把t选大一些，看看响应曲线的结果，最后再确定其合适的仿真时间。 一般来说，先不指定仿真时间，由MATLAB自己确定，然后根据结果，最后确定合适的仿真时间。 在指定仿真时间时，步长的不同会影响到输出曲线的
20、光滑程度，一般不易取太大。 例exp4_6_.m,2、常用时域分析函数,时间响应探究系统对输入和扰动在时域内的瞬态行为，系统特征如：上升时间、调节时间、超调量和稳态误差都能从时间响应上反映出来。MATLAB除了提供前面介绍的对系统阶跃响应、冲激响应等进行仿真的函数外，还提供了大量对控制系统进行时域分析的函数，如： covar：连续系统对白噪声的方差响应 initial：连续系统的零输入响应 lsim：连续系统对任意输入的响应 对于离散系统只需在连续系统对应函数前加d就可以，如dstep，dimpulse等。 它们的调用格式与step、impulse类似，可以通过help命令来察看自学。,时域响
21、应仿真 时域响应MATLAB仿真的方法有两种： 一种是在MATLAB的函数指令方式下进行时域仿真； 另一种是在SIMULINK环境下的菜单方式的时域仿真。 （1）命令方式下的时域仿真 控制系统工具箱中提供了一系列关于时域响应求取的 函数命令。 阶跃响应：step（sys） step（sys，t） dstep（a，b，c，d） dstep（num，den） 脉冲响应： impulse（sys） impulse （sys，t） dimpulse （a，b，c，d） dimpulse （num，den）,例 5 用MATLAB绘制典型二阶系统的单位阶跃响 应曲线。 解：程序实现如下： c=0 3 7
22、 10 20 40; k=25; t=linspace(0,10,100); num=k; for i=1:6den=1 c(i) 25;sys=tf(num,den);y(:,i)=step(sys,t); end plot(t,y(:,1:6) ),例 6 用MATLAB绘制三阶系统的单位阶跃响应曲线。 解：程序实现如下：rangek=0.12 0.24 0.48 0.6 3 180; rangep=0.12 0.24 0.48 0.6 3 180; t=linspace(0,10,100); for j=1:6 num=rangek(j); den=1 1.2+rangep(j) 1+1
23、.2*rangep(j) rangep(j); sys=tf(num,den); y(:,j)=step(sys,t); end plot(t,y(:,1:6) gtext(=0.2);gtext(=0.4);gtext(=0.8);gtext(=1);gtext(=5);gtext(=); hold on； plot(t,1),三 控制系统的频域分析,频率响应是指系统对正弦输入信号的稳态响应，从频率响应中可以得出带宽、增益、转折频率、闭环稳定性等系统特征。 频率特性是指系统在正弦信号作用下，稳态输出与输入之比对频率的关系特性。频率特性函数与传递函数有直接的关系，记为：,1、频域分析的一般方法
24、,求取系统对数频率特性图（波特图）：bode() 求取系统奈奎斯特图（幅相曲线图或极坐标图）：nyquist(),频域分析法是应用频率特性研究控制系统的一种典型方法。采用这种方法可直观地表达出系统的频率特性，分析方法比较简单，物理概念比较明确，对于诸如防止结构谐振、抑制噪声、改善系统稳定性和暂态性能等问题，都可以从系统的频率特性上明确地看出其物理实质和解决途经。通常将频率特性用曲线的形式进行表示，包括对数频率特性曲线和幅相频率特性曲线简称幅相曲线，MATLAB提供了绘制这两种曲线的函数。,1、对数频率特性图（波特图） exp4_10.m exp4_10_.m,对数频率特性图包括了对数幅频特性图
25、和对数相频特性图。横坐标为频率w，采用对数分度，单位为弧度/秒；纵坐标均匀分度，分别为幅值函数20lgA(w)，以dB表示；相角，以度表示。,MATLAB提供了函数bode()来绘制系统的波特图，其用法如下：,bode(a,b,c,d)：自动绘制出系统的一组Bode图，它们是针对连续状态空间系统a,b,c,d的每个输入的Bode图。其中频率范围由函数自动选取，而且在响应快速变化的位置会自动采用更多取样点。 bode(a,b,c,d,iu)：可得到从系统第iu个输入到所有输出的波特图。 bode(num,den)：可绘制出以连续时间多项式传递函数表示的系统的波特图。 bode(a,b,c,d,i
26、u,w)或bode(num,den,w)：可利用指定的角频率矢量绘制出系统的波特图。 当带输出变量mag,pha,w或mag,pha引用函数时，可得到系统波特图相应的幅值mag、相角pha及角频率点w矢量或只是返回幅值与相角。相角以度为单位，幅值可转换为分贝单位：magdb=20log10(mag),2、奈奎斯特图（幅相频率特性图） exp4_11.m exp4_11_.m,对于频率特性函数G(jw)，给出w从负无穷到正无穷的一系列数值，分别求出Im(G(jw)和Re(G(jw)。以Re(G(jw) 为横坐标， Im(G(jw) 为纵坐标绘制成为极坐标频率特性图。,MATLAB提供了函数nyq
27、uist()来绘制系统的极坐标图，其用法如下：,nyquist(a,b,c,d)：绘制出系统的一组Nyquist曲线，每条曲线相应于连续状态空间系统a,b,c,d的输入/输出组合对。其中频率范围由函数自动选取，而且在响应快速变化的位置会自动采用更多取样点。 nyquist(a,b,c,d,iu)：可得到从系统第iu个输入到所有输出的极坐标图。 nyquist(num,den)：可绘制出以连续时间多项式传递函数表示的系统的极坐标图。 nyquist(a,b,c,d,iu,w)或nyquist(num,den,w)：可利用指定的角频率矢量绘制出系统的极坐标图。 当不带返回参数时，直接在屏幕上绘制出
28、系统的极坐标图（图上用箭头表示w的变化方向，负无穷到正无穷） 。当带输出变量re,im,w引用函数时，可得到系统频率特性函数的实部re和虚部im及角频率点w矢量（为正的部分）。可以用plot(re,im)绘制出对应w从负无穷到零变化的部分。,2、常用频域分析函数,MATLAB除了提供前面介绍的基本频域分析函数外，还提供了大量在工程实际中广泛应用的库函数，由这些函数可以求得系统的各种频率响应曲线和 特征值。如：,margin：求幅值裕度和相角裕度及对应的转折频率 freqs：模拟滤波器特性 nichols：求连续系统的尼科尔斯频率响应曲线（即对数幅相曲线） ngrid：尼科尔斯方格图,margi
29、n()函数 exp4_12.m exp4_12_.m,margin函数可以从频率响应数据中计算出幅值裕度、相角裕度以及对应的频率。幅值裕度和相角裕度是针对开环SISO系统而言，它指示出系统闭环时的相对稳定性。当不带输出变量引用时，margin可在当前图形窗口中绘制出带有裕量及相应频率显示的Bode图，其中幅值裕度以分贝为单位。,幅值裕度是在相角为-180度处使开环增益为1的增益量，如在-180度相频处的开环增益为g，则幅值裕度为1/g；若用分贝值表示幅值裕度，则等于：-20*log10(g)。类似地，相角裕度是当开环增益为1.0时，相应的相角与180度角的和。 margin(mag,phase
30、,w)：由bode指令得到的幅值mag（不是以dB为单位） 、相角phase及角频率w矢量绘制出带有裕量及相应频率显示的bode图。 margin(num,den) ：可计算出连续系统传递函数表示的幅值裕度和相角裕度并绘制相应波特图。类似，margin(a,b,c,d)可以计算出连续状态空间系统表示的幅值裕度和相角裕度并绘制相应波特图。 gm,pm,wcg,wcp=margin(mag,phase,w)：由幅值mag（不是以dB为单位） 、相角phase及角频率w矢量计算出系统幅值裕度和相角裕度及相应的相角交界频率wcg、截止频率wcp，而不直接绘出Bode图曲线。,freqs()函数 exp
31、4_13.m,freqs用于计算由矢量a和b构成的模拟滤波器H(s)=B(s)/A(s)的幅频响应。,h=freqs(b,a,w)用于计算模拟滤波器的幅频响应，其中实矢量w用于指定频率值，返回值h为一个复数行向量，要得到幅值必须对它取绝对值，即求模。 h,w=freqs(b,a)自动设定200个频率点来计算频率响应，这200个频率值记录在w中。 h,w=freqs(b,a,n)设定n个频率点计算频率响应。 不带输出变量的freqs函数，将在当前图形窗口中绘制出幅频和相频曲线，其中幅相曲线对纵坐标与横坐标均为对数分度。,（2）伯德图 MATLAB中绘制系统的伯德图使用命令函数bode（） 实现。
32、其调用格式为： bode（sys） 得到关于连续系统的伯德图； dbode （num，den，Ts） dbode（a，b，c，d，Ts） 得到关于离散系统的伯德图。 （3）系统相对稳定性分析 控制理论中用幅值裕度和相角裕度评价系统相对稳定性， MATLAB中提供命令函数margin（）从频域响应数据 中计算出幅值裕度和相角裕度及其对应的角频率。 gm,pm,wcp,wcg=margin(mag,phase,w) 其中，wcg为剪切频率 , pm为相角裕度, wcp为穿越频率, gm为幅值裕度Kg,系统频域分析仿真,频域分析法是控制系统研究中应用频率特性来分析系 统性能：稳定性、快速性和稳态精度
33、的方法。此法不 必直接求解系统的微分方程，而是间接地运用系统开 环频率特性曲线分析闭环系统的响应。 MATLAB中提供了相应绘制频率特性曲线的命令函数。 （1）nyquist曲线图 调用格式为： nyquist（sys） 得到连续系统的极坐标图。 dnyquist（num，den，Ts） dnyquist（a，b，c，d，Ts） 得到关于离散系统的极坐标图。,例 10 已知系统开环传递函数为试利用MATLAB画出系统的奈奎斯特图。 解： num=0,0,4, den=3,7,2,0, nyquist(num,den), grid on, title(Nyquist plot G（s）=4/（3
34、s3+7s2+2s)）,例 11 已知系统的开环传递函数为试利用MATLAB绘制系统的奈奎斯特图 解： 其程序：num=0 0 500； den=1 1 100 0； nyquist(num,den)；,例 12 已知系统的开环传递函数为试利用MATLAB绘制系统的伯德图。 解： 其程序：num=0 0 500; den=1 1 100 0; bode(num,den); grid on;,例 13 已知一单位反馈系统开环传递函数为 G(s)=2/(s3+6s2+5s)，试绘制bode图并计算系统频 域性能指标。 解： 程序如下： num=0 0 0 2; den=1 6 5 0;sys=tf(num,den);mag,phase,w=bode(sys);gm,pm,wcp,wcg=margin(mag,phase,w) 执行结果： gm =15.0000 pm =65.1454 wcp =2.2361 wcg =0.3731margin(sys),
