第3章 控制系统特性分析的MATLAB的实现 3.1 控制系统稳定性分析的MATLAB实现 一.直接求根判定系统稳定性 求解控制系统闭环特征方程的根并判断所有根的实部是否小于零，在MATLAB里这是很容易用函数roots( )实现的。 【例3.1】已知系统开环传递函数为： ，试判别系统的稳定性。 【解】根据题意，利用roots( )函数给出以下MATLAB程序段：   k=100;z=[-2];p=[0,-1,-20]; [n1,d1]=zp2tf(z,p,k)； G=tf(n1,d1)； p=n1+d1； roots(p) %运行结果：n1 =[ 0 0 100 200]，d1 =[ 1 21 20 0]； %运行结果： Transfer function: 100 s + 200 ------------------- s^3 + 21 s^2 + 20 s ans =   -12.8990 -5.0000 -3.1010 %运行结果：p =[ 1 21 120 200]

【例3.2】已知系统的动态结构图如图3.1-1所示，试对系统闭环判别其稳定性。 【解】根据题意，利用roots( )函数给出以下MATLAB程序段： n1=[ 10];d1=[1 1 0];s1=tf(n1,d1)； n1=[2 0];d1=[ 0 1];s2=tf(n1,d1)； s12=feedback(s1,s2)； %运行结果：Transfer function: 10 ------ s^2 + s %运行结果：Transfer function: 2 s %运行结果：Transfer function: 10 ---------- s^2 + 21 s

n3=[ 1 1];d3=[1 0];s3=tf(n3,d3)； sys1=s12*s3； sys=feedback(sys1,1) roots(sys.den{1}) %运行结果：Transfer function: s + 1 ----- s %运行结果：Transfer function: 10 s + 10 ------------ s^3 + 21 s^2 %运行结果：Transfer function: 10 s + 10 ------------------------ s^3 + 21 s^2 + 10 s + 10 ans =   -20.5368 -0.2316 + 0.6582i -0.2316 - 0.6582i

补充知识：MATLAB基础Ⅱ 一.M文件 所谓M文件，就是用户把要实现的命令写在一个以.m为扩展名的文件中。与在命令窗口中输入命令行方式比，M文件的有的是可调试、可重复使用。M文件分为程序文件与函数(function)文件两大类。 1.M程序文件 在MATLAB桌面建立新的M文件或打开已经建立好的M文件（点击桌面操作环境画面中File→New→m File或File→Open→所需文件）等操作，或直接在命令窗口键入“edit”，或直接点击图标 ，都能进入图3.1-2所示的MATLAB编辑和调试操作环境的窗口画面。

如【例3.2】可以用M文件方式建立，文件名为：zhang3L2.m 图3.1-2 MATLAB编辑和调试操作环境的窗口画面 如【例3.2】可以用M文件方式建立，文件名为：zhang3L2.m

2.函数式M文件 在函数文件的第一行必须是以关键字“function”开始的函数说明语句。下面是一个只有两行的函数文件的例子。 Function c=myfile(a,b) C=sqrt((a^2)+(b^2)); >>a=3; >>b=4; >>c=myfile(a,b); c= 5 >>

二.MATLAB控制系统工具箱函数 1.主要模型建立函数表

举例： 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)； [A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,inp1,out1)； [num，den]=feedback(num1，den1，num2，den2)； [num，den]=feedback(num1，den1，num2，den2，sign)； [A，B，C，D]=feedback(A1，B1，C1，D1，A2，B2，C2，D2，sign)可将两个系统按反馈方式连接，系统1的所有输出连接到系统2的输入，系统2的所有输出连接到系统1的输入；sign是用于反馈连接的符号，默认为负号。 [num，den]=feedback(num1，den1，num2，den2，sign)；用于传递函数形式所表示的系统，sign是用于反馈连接的符号，默认为负号。

[A，B，C，D]=feedback(A1，B1，C1，D1，A2，B2，C2，D2，inp1，out1)将系统1的指定输出（out1）连接到系统2的输入，系统2的输出连接到系统1的指定输入（inp1），以次构成闭环系统，如图所示。

在【例3.2】就用到了该函数。

2.主要模型变换函数表

3.主要模型特性函数表

4.主要时域响应函数表

《线控》例2-11已知系统矩阵 ，输入矩阵 ， ， 且 ，单输入u(t)为单位阶跃函数，试求系统的状态响应和输出响应。 解：直接用MATLAB的lsim命令求解 A=[0 1;-2 -3]; B=[0 1]'; C=[1 0]; D=0; x0=[0 0.5]; [y,x]=lsim(A,B,C,D,1+t*0,t,x0); plot(x);

5.主要频域响应函数表

6.主要的根轨迹函数表

二．用根轨迹法判断系统稳定性及其举例 【例3.3】已知一个单位负反馈系统开环传递函数为G(s)，试在系统闭环的根轨迹图上选择一点，求出该点的增益k及其系统的闭环极点位置，并判断在该点系统闭环的稳定性。 【解】根据题目要求，调用函数命令rlocfind( )的程序如下： num=[1 3];den=conv(conv(conv([1 0],[1 5]),[1 6]),[1 2 2])； sys=tf(num,den)； rlocus(sys) %运行结果： den = 1 13 54 82 60 0 图3.1-2 %运行结果：Transfer function: s + 3 -------------------------------------------------- s^5 + 13 s^4 + 54 s^3 + 82 s^2 + 60 s

[k,poles]=rlocfind(sys) selected_point = -0.6498 - 0.0585i k = 7.2349 poles = -5.8248 -5.2067 -0.6585 + 0.8155i -0.6585 - 0.8155i -0.6514 Select a point in the graphics window 求给定一组根的系统根轨迹增益函数rlocfind( ) 函数命令调用格式： [k,poles]=rlocfind(sys) 函数命令使用说明：

【例3.4】续【例3.3】，试计算当k在33—37范围内时系统的闭环极点位置，并判断系统闭环的稳定性。 【解】根据题目要求，用函数命令编写MATLAB程序如下： num=[1 3];den=conv(conv(conv([1 0],[1 5]),[1 6]),[1 2 2]); cpole=rlocus(num,den,[33:1:37])； %运行结果： cpole = Columns 1 through 4   -5.5745 + 0.6697i -5.5745 - 0.6697i -1.7990 -0.0260 + 1.3210i -5.5768 + 0.6850i -5.5768 - 0.6850i -1.8154 -0.0155 + 1.3340i -5.5791 + 0.7001i -5.5791 - 0.7001i -1.8313 -0.0052 + 1.3467i -5.5815 + 0.7147i -5.5815 - 0.7147i -1.8466 0.0048 + 1.3591i -5.5838 + 0.7291i -5.5838 - 0.7291i -1.8615 0.0146 + 1.3712i Column 5 -0.0260 - 1.3210i -0.0155 - 1.3340i -0.0052 - 1.3467i 0.0048 - 1.3591i 0.0146 - 1.3712i

range=[33:1:37]'； [range,cpole] ans = Columns 1 through 4   Columns 1 through 4 33.0000 -5.5745 + 0.6697i -5.5745 - 0.6697i -1.7990 34.0000 -5.5768 + 0.6850i -5.5768 - 0.6850i -1.8154 35.0000 -5.5791 + 0.7001i -5.5791 - 0.7001i -1.8313 36.0000 -5.5815 + 0.7147i -5.5815 - 0.7147i -1.8466 37.0000 -5.5838 + 0.7291i -5.5838 - 0.7291i -1.8615 Columns 5 through 6 -0.0260 + 1.3210i -0.0260 - 1.3210i -0.0155 + 1.3340i -0.0155 - 1.3340i -0.0052 + 1.3467i -0.0052 - 1.3467i 0.0048 + 1.3591i 0.0048 - 1.3591i 0.0146 + 1.3712i 0.0146 - 1.3712i

rlocus( )函数命令用来绘制SISO系统的根轨迹图。给定前向通道传递函数G(s)，反馈补偿为k*F(s)的受控对象，其闭环传递函数为： 函数命令使用说明： rlocus( )函数命令用来绘制SISO系统的根轨迹图。给定前向通道传递函数G(s)，反馈补偿为k*F(s)的受控对象，其闭环传递函数为： 函数命令调用格式： [r,k]=rlocus(a,b,c,d) [r,k]=rlocus(sys) 可以用以下程序来校核当k=35与k=36时闭环系统的阶跃给定响应曲线： n1=[1 3];d1=conv(conv(conv([1 0],[1 5]),[1 6]),[1 2 2]); for k=35:36 n=k*n1; s1=tf(n,d1); G=feedback(s1,1); step(G);hold on end gtext('k=35'),gtext('k=36') 图3.1-.3

step(sys) 函数命令调用格式如下： 三. 用Bode图法判断系统稳定性及其举例 求连续系统单位阶跃响应的函数step( ) 图3.1-4 函数命令调用格式如下： step(sys) 三. 用Bode图法判断系统稳定性及其举例 【例3.5】已知两个单位负反馈系统开环传递函数分别为： 用Bode图法判断系统闭环的稳定性。 【解】根据要求，对第一个系统用margin( )函数命令给出如下程序： num=[0 0 0 2.7];den=[1 5 4 0];s1=tf(num,den); [Gm,Pm,Wcp,Wcg]=margin(s1) Gm = 7.4074 Pm = 51.7320 Wcp = 2 Wcg = 0.5783 margin(s1)

[Gm,Pm,Wcp,Wcg]=margin(sys) margin(sys) 函数命令调用格式： [Gm,Pm,Wcp,Wcg]=margin(sys) margin(sys) 同样对第二个系统给出程序如下： n2=[0 0 0 2.7];d2=[1 5 -4 0];s2=tf(n2,d2); margin(s2) [Gm,Pm,Wcp,Wcg]=margin(s2) 图3.1-5 Gm = 1 Pm =

四. 李雅普诺夫(Lyapunov)稳定性判据 复习： 《现代控制理论》 [定理4.8] 线性系统渐进稳定的判别方法 设线性定常连续系统为 则平衡状态 为大范围渐进稳定的充要条件是：对任意给定的一个正定实对称矩阵Q，比必存在一个正定的实对称矩阵P，且满足李雅普诺夫方程 并且 是系统的李雅普诺夫函数。

在《线代控制理论》的理论计算中，要判系统（3.1-1）是否稳定,需要做以下工作: ① 由式(3.1-4) 求P，并验证P是正定的； ② 由式(3.1-3) 求V(x)，并判验证V(x) 是正定的； ③ 结论：系统是稳定的。 在MATLAB中，Lyapunov方程（3.1-2）可以由控制系统工具箱中提供的lyap( )函数求解，该函数的调用格式： V=lyap(A,W)。 [例3.6] 已知系统的状态方程为： 试分析系统的稳定性。

3.2 控制系统稳态误差分析的MATLAB实现 一.控制系统稳态误差分析的有关概念 1.稳态误差ess 2.自动控制系统的型别 3.控制系统的稳态误差系数 (1)稳态位置误差系数kp (2) 稳态速度误差系数kv (3) 稳态加速度误差系数ka

二. 控制系统稳态误差计算举例 【例3.7】已知一个单位负反馈系统开环零极点增益模型为： 试绘制出该系统的单位斜坡响应曲线并求单位斜坡响应稳态误差。 【解】（1）对系统判稳 程序： k=6;z=-0.5;p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k)； s=tf(n1,d1)； sys=feedback(s,1)； roots(sys.den{1}) %运行结果： n1 = 0 0 6 3 d1 = 1 1 -2 0 %运行结果： Transfer function: 6 s + 3 --------------- s^3 + s^2 - 2 s ans =   -0.1084 + 1.9541i -0.1084 - 1.9541i -0.7832 % Transfer function: 6 s + 3 ------------------- s^3 + s^2 + 4 s + 3

subplot(121),plot(t,y),grid； （2）求系统单位阶跃给定响应与稳态误差 程序： k=6;z=-0.5;p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k); s=tf(n1,d1); sys=feedback(s,1); step(sys)； subplot(121),plot(t,y),grid； subplot(122),ess=1-y； plot(t,ess),grid 图4.6 系统阶跃响应与阶跃误差响应曲线 t=[0:0.1:30]'; y=step(sys,t); %运行结果：y = 0.0294 ： 1.0224 1.0322 ess(length(ess))   ans = -0.0322

求单位阶跃响应的稳态误差函数esst( ) [ess]=esst(sys,t) % MATLAB PROGRAM esst.m   function [ess]=esst(sys,t) y=step(sys,t); subplot(121),plot(t,y),grid subplot(122),es=1-y; plot(t,es),grid ess=es(length(es))

k=6;z=-0.5;p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k); s=tf(n1,d1); 如果运行以下调用esst.m函数的程序，会得到同样的结果。 k=6;z=-0.5;p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k); s=tf(n1,d1); sys=feedback(s,1); t=[0:0.1:30]' [ess]=esst(sys,t); ans = -0.0322

3.3 用MATLAB分析系统能控能观测性的方法 一.能控能观性判别 系统能观测 利用这两个函数很容易判定系统的能控性和能观测性。

[例3.9] 选自《现代控制理论》 已知系统的状态空间表达式为： 用MATLAB函数分析系统能控性与能观测性。

MATLAB中，提供了可将连续或离散系统状态空间表达式化为对角线或约旦标准型的jordan()函数，该函数的格式为： [v，j]＝jordan(A)。 [例3.10] 已知系统的状态空间表达式为： 用MATLAB函数将系统化为对角线或约旦标准型，并分析系统能控性与能观测性。

[例3.11] 已知系统的状态空间表达式为： 用MATLAB函数将系统化为对角线或约旦标准型，并分析系统能控性与能观测性。

二.系统的结构分解 复习： 定理[3.15] 设一个不能控系统的状态方程为: 能控子空间 不能控子空间

MATLAB中，提供了将系统进行能控与不能控分解的函数ctrbf（ ），该函数的调研格式为： [AB，BB，CB，T，K]=ctrbf(A,B,C) 同样，设一个不能观测系统的状态方程为: 能观测子空间 不能观测子空间

MATLAB中，提供了将系统进行能观测与不能观测分解的函数obsvf（ ），该函数的调研格式为： [AB，BB，CB，T，K]=obsvf(A,B,C) [例3.12] 已知系统的状态空间表达式为： 对系统进行能控与不能控分解和能观测与不能观测分解。

程序运行结果： >> n = 2 AB = 1.0000 -0.0000 -0.0000 2.1213 4.0000 -1.2247 1.7321 2.4495 0 BB = 0.0000 -0.0000 -1.4142 CB = -1.6330 -0.5774 0   T = -0.8165 0.4082 -0.4082 0.5774 0.5774 -0.5774 0 -0.7071 -0.7071  K =  1 1 0 m = 2 AB1 = 1.0000 -2.8868 -3.5355 -0.0000 2.0000 1.2247 -0.0000 0.8165 2.0000  BB1 = 1.2247 0.7071  CB1 = 0 0 -1.7321  T1 = -0.4082 0.4082 0.8165 0.7071 0.7071 0.0000 -0.5774 0.5774 -0.5774 K1 = 1 1 0 >>

系统按能控与不能控分解为新的状态空间表达式： 系统按能观测与不能观测分解为新的状态空间表达式：

P=floplr(eye(size(sys,’order’)))。 令： 则：

结果： AB2 = 0 2.4495 1.7321 -1.2247 4.0000 2.1213 -0.0000 -0.0000 1.0000   BB2 = -1.4142 -0.0000 0.0000 CB2 = 0 -0.5774 -1.6330