当前位置:首页 >> 信息与通信 >>

控制系统计算机仿真与CAD-05-8


控制系统计算机仿真与CAD
陈 在 平

自动化与能源工程学院 天津理工大学
1

学时安排
总学时:54(其中实验8学时)

2

第一章 绪言
第一节 控制系统计算机仿真和辅助设计概述
在进行自动控制系统分析与设计综合工作过程中 除了需要进行理

论分析之外还要对系统的特性进行实 验研究。 控制系统计算机仿真与辅助设计是目前对复杂控 制系统进行分析设计重要手段与方法。

3

? 控制系统计算机仿真的概念
仿真的概念 计算机仿真的概念

?

控制系统计算机仿真是以控制系统的数学模 型为基础,借助计算机对控制系统的动静态过 程进行实验研究(仿真分析)。 其特点:将实际系统的运动规律用数学表 达式描述(通常是一组常微分方程或差分方 程),然后利用计算机来求解以达到对系统进 行分析研究的目的。
4

计算机辅助设计的概念
?
? ? 计算机辅助设计(CAD) 利用计算机高速而精确的计算能力大容量存储 和处理数据的能力 加快设计进程、缩短设计周期、提高设计质量。 控制系统计算机辅助设计(CS-CAD) 根据给定的系统性能指标的要求借助计算机进 行控制系统的分析设计。 计算机仿真实际就是计算机辅助设计工作的一 部分,以精确模型为验证对象,检验系统的设 计效果。
5

?

计算机仿真模型
? 模型的概念 ? 模型是人们认识事物的一种手段或工具。
模型的分类 (1)物理模型 (2)数学模型 (3)仿真模型

6

控制系统计算机仿真的三要素及基本步骤
1、三要素 (1)系统(2)模型(3)计算机
? ? ? ? ? ? ? ? ? ? ¨? ? §? ? ? ? ? ? ·? ? ? ? ? ? ? ? ? ? ? ? ú ? ? ·? ? ? ¨? ? ? ? ? ? ? ·? á ? ? ?

? 1.1 ? ? ? ·? ? ? ? ? ? ? ? ? ? ú ? ? ? ? × × ? ?

7

2、基本步骤 ? 包括三个基本的内容:建模 结果分析 ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ± ? ? ? ? ¨? ? ? à ? ò ±? ? · ? ? ¤? ·? ? ? ? ? ? ? · ? ? ·? ? ? ? ? ·? ? ? ? ? ? ? ? ? ? ? ? ? ? ·? ? ? ? ? á ? ? ? ? ?
? ? 1.2 ? ? ? ·? ? ? ? ? ? ? ú ? ? ? ò ÷?

仿真实验

· ?

8

第二节

仿真的分类

一、按模型分类
1、物理仿真:采用物理模型,有实物介入 ? 具有效果逼真,精度高等优点,但造价高或耗 时长,具有实时性、在线的特点。 2、数学仿真:采用数学模型 ? 在计算机上进行,具有非实时性、离线的特点, 经济、快速、实用。

9

二、按计算机类型分类
1、模拟仿真:采用数学模型,在模拟计算机上 进行的实验研究。(50年代) ? 描述连续物理系统的动态过程比较自然、逼 真,具有速度快、失真小、结果可靠的优点, 但精度较低,对计算机控制系统的仿真较困 难,自动化程度低。 ? 模拟计算机的核心是运算部分,由 “模拟运 算放大器”为主要部件所构成。 2、数字仿真:采用数学模型,在数字计算机上 借助于数值计算进行的仿真实验。(60年代)
10

? 计算与仿真的精度较高。 ? 对计算机控制系统的仿真比较方便。

? 计算速度比较低,在一定程度上影响到仿真 的可信度。 ? 数字仿真需要用高级程序语言编写求解系统 模型及结果输出的程序。

11

3、混合仿真:结合了模拟仿真与数字仿真。

12

第三节 仿真技术的应用与发展
一、仿真技术在工程中的应用
1、航空与航天工业 ? 飞行器设计中的三级仿真体系 ? 飞行员及宇航员训练用飞行仿真模拟器。 2、电力工业 ? 电力系统动态模型实验 ? 电站操作人员培训模拟系统。

13

3、原子能工业 ? 模拟核反应堆 ? 核电站仿真器

4、石油、化工及冶金工业
5、非工程领域 ? 医学 ? 社会学 ? 宏观经济与商业策略
14

二、应用仿真技术的意义 1、经济 ? 大型、复杂系统直接实验是十分昂贵的, 仿真实验仅需其成本的1/10~1/5,而且 可以重复使用。 2、安全 ? 某些系统(如载人飞行器、核电装置 等),直接实验往往会有很大的危险, 而仿真实验可以有效降低危险程度。
15

3、快捷 ? 提高设计效率 4、具有优化设计和预测的特殊功能 ? 对一些真实系统进行优化设计非常困难时, 仿真可以发挥特殊的优化设计功能。 ? 在非工程系统中(如社会、管理、经济等系 统),直接实验几乎不可能,而通过可以获 得对系统的某种超前认识。
16

三、仿真技术的发展趋势
1、硬件方面: 基于多CPU并行处理技术大大增强数字仿真 的实时性。 2、应用软件方面: 直接面向用户的数字仿真软件不断推陈出新, 各种专家系统与智能化技术更加深入。 3、分布式数字仿真: 充分利用网络技术,协调合作,投资少,效果 好。
17

4、虚拟现实技术: 综合了计算机图形技术、多媒体技术、传感器 技术、显示技术以及仿真技术等多学科,使人 置身于真实环境之中。 5、在线实时仿真: 随着仿真环境的不断丰富,仿真开始向仿真-控 制一体化过渡,向进行实时仿真发展。

18

第四节 计算机仿真软件
一、仿真软件的发展
1、程序编程阶段 ? 所有问题都用高级语言编写。 2、程序软件包阶段 ? 出现了“应用子程序库”。 3、交互式语言阶段(仿真语言) ? 仿真语言可用一条指令实现某种功能 4、模型化图形组态阶段 ? 符合基于模型图形化的描述。
19

二、几种仿真软件
1、PSPICE、ORCAD:通用的电子电路仿真软 件,适合于元件级仿真。 2、SYSTEM VIEW:系统级的电路动态仿真 软件 3、MATLAB:具有强大的数值计算能力,包含 各种工具箱 4、SIMULINK:是MATLAB附带的基于模型化 图形组态的动态仿真环境。

20

第五节 MATLAB语言基础
使用MATLAB的窗口环境 一、MATLAB语言的显著特点 1、具有强大的矩阵运算能力:Matrix Laboratory(矩 阵实验室),使得矩阵运算非常简单。 2、是一种演算式语言 ? MATLAB的基本数据单元是既不需要指定维数,也不 需要说明数据类型的矩阵(向量和标量为矩阵的特 例),而且数学表达式和运算规则与通常的习惯相同。 因此MATLAB语言编程简单,使用方便。 例 exp2_1.m

21

二、MATLAB命令窗口
1、启动MATLAB命令窗口 ? 计算机安装好MATLAB之后,双击MATLAB图 标,就可以进入命令窗口,此时意味着系统处 于准备接受命令的状态,可以在命令窗口中直 接输入命令语句。 ? MATLAB语句形式 》变量=表达式; 通过等于符号将表达式的值赋予变量。当 键入回车键时,该语句被执行。 如果希望结果不被显示,则只要在语句之 后加上一个分号(;)即可。但它依然被赋值 并在MATLAB工作空间中分配了内存。 22

2、命令行编辑器 (1)方向键和控制键可以编辑修改已输入的命 令 :回调上一行命令 :回调下一行命 令 (2)命令窗口的分页输出 ? more off:不允许分页 more on:允许分页 ? more(n):指定每页输出的行数 ? 回车前进一行,空格键显示下一页,q结束当 前显示。 (3)多行命令(…) ? 如果命令语句超过一行或者太长希望分行输入, 则可以使用多行命令继续输入。 23 ? S=1-12+13+4+…

三、变量和数值显示格式
1、变量 (1)变量的命名:变量的名字必须以字母开头 (不能超过19个字符),之后可以是任意字母、 数字或下划线;变量名称区分字母的大小写; 变量中不能包含有标点符号。 (2)一些特殊的变量 ans:用于结果的缺省变量名 i、j:虚数单位 pi:圆周率 nargin:函数的输入变量个数 eps:计算机的最小数 nargout:函数的输出变量 个数 inf:无穷大 realmin:最小正实数 realmax:最大正实数 nan:不定量 24 flops:浮点运算数

(3)变量操作 ? 在命令窗口中,同时存储着输入的命令和创建 的所有变量值,它们可以在任何需要的时候被 调用。如要察看变量a的值,只需要在命令窗 口中输入变量的名称即可:》a 2、数值显示格式 ? 任何MATLAB的语句的执行结果都可以在屏幕 上显示,同时赋值给指定的变量,没有指定变 量时,赋值给一个特殊的变量ans,数据的显 示格式由format命令控制。 ? format只是影响结果的显示,不影响其计算与 存储;MATLAB总是以双字长浮点数(双精度) 来执行所有的运算。
25

? 如果结果为整数,则显示没有小数;如果结果不 是整数,则输出形式有: format (short):短格式(5位定点数)99.1253 format long:长格式(15位定点数 99.12345678900000 format short e:短格式e方式 9.9123e+001 format long e:长格式e方式 9.912345678900000e+001 format bank:2位十进制 99.12 format hex:十六进制格式
26

四、简单的数学运算(例exp2_2.m)
1、常用的数学运算符 ? +,—,*(乘),/(左除),\(右除),^(幂) ? 在运算式中,MATLAB通常不需要考虑空格;多条 命令可以放在一行中,它们之间需要用分号隔开; 逗号告诉MATLAB显示结果,而分号则禁止结果显 示。 2、常用数学函数 abs,sin,cos,tan,asin,acos,atan,sqrt,exp,imag,real,si gn, log,log10,conj(共扼复数)等

五、MATLAB的工作空间
1、MATLAB的工作空间包含了一组可以在命令窗口中 调整(调用)的参数 ? who:显示当前工作空间中所有变量的一个简单列表

27

? whos:则列出变量的大小、数据格式等详细信 息 ? clear :清除工作空间中所有的变量 ? clear 变量名:清除指定的变量 2、保存和装入workspace (1)save filename variables ? 将变量列表variables所列出的变量保存到磁盘 文件filename中 ? Variables所表示的变量列表中,不能用逗号, 各个不同的变量之间只能用空格来分隔。 ? 未列出variables时,表示将当前工作空间中所 有变量都保持到磁盘文件中。 ? 缺省的磁盘文件扩展名为“.mat”,可以使用 “-”定义不同的存储格式(ASCII、V4等) 28

(2)load filename variables ? 将以前用save命令保存的变量variables从磁盘 文件中调入MATLAB工作空间。 ? 用load 命令调入的变量,其名称为用save命令 保存时的名称,取值也一样。 ? Variables所表示的变量列表中,不能用逗号, 各个不同的变量之间只能用空格来分隔。 ? 未列出variables时,表示将磁盘文件中的所有 变量都调入工作空间。 3、退出工作空间 ? quit 或 exit

29

六、文件管理
? 文件管理的命令,包括列文件名、显示或删除 文件、显示或改变当前目录等。(what、dir、 type、delete、cd、which) ? what:显示当前目录下所有与matlab相关的文 件及它们的路径。 ? dir:显示当前目录下所有的文件 ? which:显示某个文件的路径 ? cd path:由当前目录进入path目录 ? cd ..:返回上一级目录 ? cd:显示当前目录 ? type filename:在命令窗口中显示文件 filename 30 ? delete filename:删除文件filename

七、使用帮助
1、help命令,在命令窗口中显示 ? MATLBA的所有函数都是以逻辑群组方式进行 组织的,而MATLAB的目录结构就是以这些群 组方式来编排的。 ? help matfun :矩阵函数-数值线性代数 ? help general:通用命令 ? help graphics:通用图形函数 ? help elfun:基本的数学函数 ? help elmat:基本矩阵和矩阵操作 ? help datafun:数据分析和傅立叶变换函数 ? help ops:操作符和特殊字符
31

? help polyfun:多项式和内插函数 ? help lang:语言结构和调试 ? help strfun:字符串函数 ? help control:控制系统工具箱函数 2、helpwin:帮助窗口 3、helpdesk:帮助桌面,浏览器模式 4、lookfor命令:返回包含指定关键词的那些项 5、demo:打开示例窗口

32

MATLAB矩阵运算及多项式处理
一、矩阵的生成 1、在命令窗口中输入
》a=1; b=2; c=3; 》x=[5 b c; a*b a+c c/b] x= 5.000 2.000 3.000 2.000 4.000 1.500 》y=[2,4, 5 3 6 8] y= 245 368

? 矩阵生成不但可以使用纯数字(含复数),也可以使 用变量(或者说采用一个表达式)。矩阵的元素直接 排列在方括号内,行与行之间用分号隔开,每行内的 元素使用空格或逗号隔开。大的矩阵可以用分行输入, 回车键代表分号。
33

2、语句生成 (1)用线性等间距生成向量矩阵(start:step:end) 》a=[1:2:10] 其中start为起始值,step为步长, a= end为终止值。当步长为1时可省略 1 3 5 7 9 step参数;另外step也可以取负数。 (2)a=linspace(n1,n2,n) ?在线性空间上,行矢量的值从n1到n2,数据 个数为n,缺省n为100。 》a=linspace(1,10,10) a= 1 2 3 4 5 6 7 8 9 10
34

(3)a=logspace(n1,n2,n) ?在对数空间上,行矢量的值从10n1到10n2, 数据个数为n,缺省n为50。这个指令为建立 对数频域轴坐标提供了方便。 》a=logspace(1,3,3) a= 10 100 1000 (4)一些常用的特殊矩阵 单位矩阵:eye(m,n); eye(m) 零矩阵:zeros(m,n); zeros(m) 一矩阵:ones(m,n); ones(m) 对角矩阵:对角元素向量 V=[a1,a2,…,an] A=diag(V) 随机矩阵:rand(m,n)产生一个m×n的均匀分别的随机35 矩阵

》eye(2,3) ans= 100 010 》zeros(2,3) ans= 000 000 》ones(2,3) ans= 111 111 》V=[5 7 2]; A=diag(V) A= 500 070 002

》eye(2) ans= 10 01 》zeros(2) ans= 00 00 》ones(2) ans= 11 11

如果已知A为方阵,则 V=diag(A)可以提取A的对角元 素构成向量V。

36

二、矩阵的运算
1、转置:对于实矩阵用(’)符号或(.’)求转置结
果是一样的;然而对于含复数的矩阵,则(’)将同 时对复数进行共轭处理,而 (.’)则只是将其排列形 式进行转置。
》a=[1 2 3;4 5 6]' 》a=[1 2 3;4 5 6].' 》b=[1+2i 2-7i]' b= 4 5 1.0000 - 2.0000i 2.0000 + 7.0000i 》b=[1+2i 2-7i].' b= 1.0000 + 2.0000i

a=
1 2 4 5

a=
1 2

3

6

3

6

2.0000 - 7.0000i

37

2、四则运算与幂运算
+ ;- ;* ;\和/ ; ^;.* ;.\ ;./;.^ 如:a=[1 2;3 4];b=[ 3 5; 5 9] ? 》c=a+b d=a-b 》c= d= ? 4 7 -2 -3 8 13 -2 -5 》a*b=[13 23; 29 51] 》a/b=[-0.50 0.50;3.50 –1.50] 》a\b=[-1 -1;2 3] 》a^3=[37 54; 81 118] 》a.*b=[3 10;15 36] ? 》a./b=[0.33 0.40;0.60 0.44] 》a.\b=[3.00 2.50;1.67 2.25] 》a.^3= [1 8; 27 64]

只有维数相同的矩阵才能进 行加减运算。 注意只有当两个矩阵中前一 个矩阵的列数和后一个矩阵 的行数相同时,才可以进行 乘法运算。a\b运算等效于 求a*x=b的解;而a/b等效于 求x*b=a的解。只有方阵才 可以求幂。 点运算是两个维数相同矩阵 对应元素之间的运算,在有 的教材中也定义为数组运算。
38

3、逆矩阵与行列式计 算 求逆:inv(A); 求行列式:det(A) 要求矩阵必须为方阵 》a=[1 2 3; 4 5 6; 2 3 5];
》b=inv(a)

4、了解矩阵超越函数 ?在MATLAB中exp、sqrt等 命令也可以作用到矩阵上, 但这种运算是定义在矩阵 的单个元素上的,即分别 对矩阵的每一个元素进行 计算。 ?超越数学函数可以在函数 后加上m而成为矩阵的超 越函数,例如: expm,sqrtm。矩阵的超越 函数要求运算矩阵为方阵。
39

b=
-2.3333 0.3333 1.0000 2.6667 0.3333 -2.0000

-0.6667 -0.3333 1.0000
》det(a) ans =

-3

三、矩阵的操作
1、矩阵下标 ?MATLAB通过确认矩阵下标,可以对矩阵进 行插入子块,提取子块和重排子块的操作。

? ? ? ?

A(m,n):提取第m行,第n列元素 A(:,n):提取第n列元素 A(m,:):提取第m行元素 A(m1:m2,n1:n2):提取第m1行到第m2 行和第n1列到 第n2列的所有元素(提 取子块)。
40

? A(:):得到一个长列矢量,该矢量的元素按矩 阵的列进行排列。 ? 矩阵扩展:如果在原矩阵中一个不存在的地址 位置上设定一个数(赋值),则该矩阵会自动 扩展行列数,并在该位置上添加这个数,而且 在其他没有指定的位置补零。 ? 消除子块:如果将矩阵的子块赋值为空矩阵[ ], 则相当于消除了相应的矩阵子块。

41

2、矩阵的大小 ? [m,n]=size(A,x):返回矩阵的行列数m与n,当x=1,则 只返回行数m,当x=2,则只返回列数n。 ? length(A)=max(size(A)):返回行数或列数的最大值。 ? rank(A):求矩阵的秩
》a=[1 2 3;3 4 5]; 》length(a) 》rank(a) ans = 2

》[m,n]=size(a)
m= 2

ans =
3 》max(size(a))

n=
3

ans =
3

3、了解矩阵操作函数:flipud;fliplr; rot90

42

四、多项式处理
(1)多项式的建立与表示方法 ? 在MATLAB中,多项式使用降幂系数的行向量表 x 4 ? 12 x 3 ? 0 x 2 ? 25 x ? 116 示,如:多项式 表示为:p=[1 -12 0 25 116],使用函数roots可 以求出多项式等于0的根,根用列向量表示。 若已知多项式等于0的根,函数poly可以求出相 应多项式。 r=roots(p) p=poly(r) r= p= 1 -12 -0 25 116 11.7473 2.7028 -1.2251 + 1.4672i -1.2251 - 1.4672i

43

(2)多项式的运算

?相乘conv a=[1 2 3] ; b=[1 2] c=conv(a,b)=1 4 7 6 conv指令可以嵌套使用,如 conv(conv(a,b),c) ?相除deconv [q,r]=deconv(c,b) q=1 2 3 %商多项式 r=0 0 0 %余多项式
44

? 求多项式的微分多项式polyder polyder(a)=2 2 ? 求多项式函数值polyval(p,n):将值n代入多项 式求解。polyval(a,2)=11

45

(3)*多项式的拟合 ? 多项式拟合又称为曲线拟合,其目的就是在众 多的样本点中进行拟合,找出满足样本点分布 的多项式。这在分析实验数据,将实验数据做 解析描述时非常有用。 ? 命令格式:p=polyfit(x,y,n),其中x和y为样本 点向量,n为所求多项式的阶数,p为求出的多 项式。 ? 例exp2_15.m (4)*多项式插值 ? 多项式插值是指根据给定的有限个样本点,产 生另外的估计点以达到数据更为平滑的效果。 该技巧在信号处理与图像处理上应用广泛。
46

? 所用指令有一维的interp1、二维的interp2、三 维的interp3。这些指令分别有不同的方法 (method),设计者可以根据需要选择适当的 方法,以满足系统属性的要求。Help polyfun 可以得到更详细的内容。 y=interp1(xs,ys,x,’method’) ? 在有限样本点向量xs与ys中,插值产生向量x 和y,所用方法定义在method中,有4种选择: ? nearest:执行速度最快,输出结果为直角转折 ? linear:默认值,在样本点上斜率变化很大 ? spline:最花时间,但输出结果也最平滑 ? cubic:最占内存,输出结果与spline差不多 例exp2_16.m
47

五*、MATLAB数据处 理

求矩阵A的奇异值及分解矩阵,满 1、矩阵分解 足U*S*V’=A,其中U、V矩阵为正 (1)奇异值分解 交矩阵(U*U’=I),S矩阵为对角 [U,S,V]=svd(A) 矩阵,它的对角元素即A矩阵的奇 异值。 [u,s,v]=svd(a) 例:a = u= 9 8 0.7705 -0.6375 6 8 0.6375 0.7705 s= 可以验证: 15.5765 0 u*u’=I 0 1.5408 v*v’=I v= u*s*v’=a 0.6907 -0.7231 48 0.7231 0.6907

(2)特征值分解 求矩阵A的特征向量V及特征 [V,D]=eig(A) 值D,满足A*V=V*D。其中D 例: a = 的对角线元素为特征值,V的 9 8 列为对应的特征向量。如果 6 8 D=eig(A)则只返回特征值。 [v,d]=eig(a) v= 可以验证:A*V=V*D 0.7787 -0.7320 0.6274 0.6813 d= 15.4462 0 0 1.5538
49

(3)正交分解 将矩阵A做正交化分解,使得 [Q,R]=qr(A) Q*R=A,其中Q为正交矩阵 例: a = (其范数为1,指令 9 8 norm(Q)=1),R为对角化的 6 8 上三角矩阵。 [q,r]=qr(a) q= norm(q) q*r -0.8321 -0.5547 ans = ans = -0.5547 0.8321 1 9.0000 r= 8.0000 -10.8167 -11.0940 6.0000 0 2.2188 8.0000
50

(4)三角分解 [L,U]=lu(A) ? 将A做对角线分解,使得A=L*U,其中L为下三角 矩阵,U为上三角矩阵。 ? 注意:L实际上是一个“心理上”的下三角矩 阵,它事实上是一个置换矩阵P的逆矩阵与一 个真正下三角矩阵L1(其对角线元素为1)的 乘积。 [L1,U1,P]=lu(A) 例:a=[1 2 3;4 5 6;7 8 9] 比较: [l1,u1,p]=lu(a) [l,u]=lu(a)
51

l1 = 1.00 0 0.14 1.00 0.57 0.50 u1 = 7.00 8.00 0 0.86 0 0 p= 0 0 1 1 0 0 0 1 0

0 0 1.00

9.00 1.71 0.00

l= 0.14 1.00 0.57 0.50 1.00 0 u= 7.00 8.00 0 0.86 0 0 可以验证:

0 1.00 0 9.00 1.71 0.00

u1=u,inv(p)*l1=l
a=l*u

p*a=l1*u1

52

2*、数据分析 (1)绘制函数图形:fplot() (2)求极值:fmin,fmins (3)求零点:寻找一维函数的过零点fzero() (4)频谱分析(fft):y=FFT(x);unwrap(); abs;angle画出幅频和相频曲线 ( 5 ) 了 解 数 据 分 析 函 数 : max,min,mean,sum,prod等 (6)了解积分运算:trap2,quad,quad8 3*、常微分方程数值解 [t,x]=ode23(‘xfun’,t0,tf,x0,tol) [t,x]=ode45(‘xfun’,t0,tf,x0,tol)
53

MATLAB绘图功能简介
? MATLAB提供了丰富的绘图功能 help graph2d可得到所有画二维图形的命令 help graph3d可得到所有画三维图形的命令 下面介绍常用的二维图形命令
1、基本的绘图命令 plot(x1,y1,option1,x2,y2,option2,…) x1,y1给出的数据分别为x,y轴坐标值,option1 为选项参数,以逐点连折线的方式绘制1个二维 图形;同时类似地绘制第二个二维图形,……等。
54

这是plot命令的完全格式,在实际应用中 可以根据需要进行简化。比如: plot(x,y);plot(x,y,option) 选项参数option定义了图形曲线的颜色、 线型及标示符号,它由一对单引号括起来。 例exp2_3.m exp2_3_.m

55

2、选择图像 figure(1);figure(2);…;figure(n) 打开不同的图形窗口,以便绘制不同的图形。
3、grid on:在所画出的图形坐标中加入栅格 grid off:除去图形坐标中的栅格 4、hold on:把当前图形保持在屏幕上不变, 同时 允许在这个坐标内绘制另外一个 图形。 hold off:使新图覆盖旧的图形 例exp2_4.m exp2_4_.m
56

5、设定轴的范围 axis([xmin xmax ymin ymax]) 例exp2_5.m axis(‘equal’):将x坐标轴和y坐标轴的单位刻度 大小调整为一样。
6、文字标示 ? text(x,y,’字符串’) 在图形的指定坐标位置(x,y)处,标示单引号括起来 的字符串。 ? gtext(‘字符串’) 利用鼠标在图形的某一位置标示字符串。
57

? title(‘字符串’) 在所画图形的最上端显示说明该图形标题的字符 串。 ? xlabel(‘字符串’),ylabel(‘字符串’) 设置x,y坐标轴的名称。 ? 输入特殊的文字需要用反斜杠(\)开头。 7、legend(‘字符串1’,‘字符串2’,…,‘字符串n’) ? 在屏幕上开启一个小视窗,然后依据绘图命令 的先后次序,用对应的字符串区分图形上的线。 例exp2_5_.m

58

8、subplot(m n k):分割图形显示窗口 m:上下分割个数,n:左右分割个数,k:子图 编号 例exp2_6.m 9、semilogx:绘制以x轴为对数坐标(以10为 底),y轴为线性坐标的半对数坐标图形。 semilogy:绘制以y轴为对数坐标(以10为 底),x轴为线性坐标的半对数坐标图形。 例exp2_17.m

59

10、了解应用型绘图指令:可用于数值统计分析 或离散数据处理 bax(x,y);hist(y,x) stairs(x,y);stem(x,y) 例exp2_7.m

60

一点补充说明

? 对于图形的属性编辑同样可以通过在图 形窗口上直接进行。 ? 但图形窗口关闭之后编辑结果不会保存

61

MATLAB程序设计入门
一、MATLBA程序的基本设计原则
1、%后面的内容是程序的注解,要善于运用注解使程序更具可读性。 2、养成在主程序开头用clear指令清除变量的习惯,以消除工作空间中 其他变量对程序运行的影响。但注意在子程序中不要用clear。 3、参数值要集中放在程序的开始部分,以便维护。要充分利用 MATLAB工具箱提供的指令来执行所要进行的运算,在语句行之后 输入分号使其及中间结果不在屏幕上显示,以提高执行速度。 4、input指令可以用来输入一些临时的数据;而对于大量参数,则通过 建立一个存储参数的子程序,在主程序中用子程序的名称来调用。 5、程序尽量模块化,也就是采用主程序调用子程序的方法,将所有子 程序合并在一起来执行全部的操作。 6、充分利用Debugger来进行程序的调试(设置断点、单步执行、连续 执行),并利用其他工具箱或图形用户界面(GUI)的设计技巧,将 设计结果集成到一起。 7、设置好MATLAB的工作路径,以便程序运行。
62

8、MATLAB程序的基本组成结构 %说明 清除命令:清除workspace中的变量和图形 (clear,close) 定义变量:包括全局变量的声明及参数值的设定 逐行执行命令:指MATLAB提供的运算指令或工具箱 … … … 提供的专用命令 控制循环 :包含for,if then,switch,while等语句 逐行执行命令 … … … end 绘图命令:将运算结果绘制出来 ? 当然更复杂程序还需要调用子程序,或与simulink以 及其他应用程序结合起来。
63

二、M文件的编辑及MATLAB工作路径的设 置
? 进入MATLAB的Editor/Debugger窗口来编辑程序 ? 在编辑环境中,文字的不同颜色显示表明文字的不同 属性。 绿色:注解;黑色:程序主体;红色:属性值的设定; 蓝色:控制流程。 ? 在运行程序之前,必须设置好MATLAB的工作路径, 使得所要运行的程序及运行程序所需要的其他文件处 在当前目录之下,只有这样,才可以使程序得以正常 运行。否则可能导致无法读取某些系统文件或数据, 从而程序无法执行。 ? 通过cd指令在命令窗口中可以更改、显示当前工作路 径。 ? 通过路径浏览器(path browser)也可以进行设置

64

三、MATLAB的程序类型
MATLAB的程序类型有三种,一种是脚本M文 件;是可以存取的M文件,包括程序文件或 在命令窗口运行的命令;一种是函数 (function)文件。 1、脚本M文件 在命令窗口中输入并执行,它所用的变量都要 在工作空间中获取,不需要输入输出参数的调 用,退出MATLAB后就释放了。
65

M文本文件,以.m格式进行存取,包含一连串 的MATLAB指令和必要的注解。需要在工作空 间中创建并获取变量,也就是说处理的数据为 命令窗口中的数据,没有输入参数,也不会返 回参数。 ? 程序运行时只需在工作空间中键入其名称即可。

66

2、函数文件

例myhilb1.m

myhilb.m

? 与在命令窗口中输入命令一样,函数接受输入 参数,然后执行并输出结果。用help命令可以 显示它的注释说明。 ? 具有标准的基本结构。 (1)函数定义行(关键字function) ? function[out1,out2,..]=filename(in1,in2,..) ? 输入和输出(返回)的参数个数分别由nargin 和nargout两个MATLAB保留的变量来给出。 (2)第一行帮助行,即H1行 ? 以(%)开头,作为lookfor指令搜索的行
67

(3) 函数体说明及有关注解 以(%)开头,用以 说明函数的作用及有关内容,如果不希望显示某 段信息,可在它的前面加空行. (4)函数体语句 ? 函数体内使用的除返回和输入变量这些在 function语句中直接引用的变量以外的所有变 量都是局部变量,即在该函数返回之后,这些 变量会自动在MATLAB的工作空间中清除掉。 如果希望这些中间变量成为在整个程序中都起 作用的变量,则可以将它们设置为全局变量。

68

四、声明子程序(函数程序)变量
1、子程序与主程序之间的数据是通过参数进行传递的, 子程序应用主程序传递来的参数进行计算后,将结果 返回主程序。 例exp2_9.m 2、如果一个函数内的变量没有特别声明,那么这个变 量只在函数内部使用,即为局部变量。如果两个或多个 函数共用一个变量(或者说在子程序中也要用到主程序 中的变量,注意不是参数),那么可以用global来将它 声明为全局变量。 例exp2_10.m 全局变量的使用可以减少参数传递,合理利用全局变 量可以提高程序执行的效率。
69

五、字符串及其宏命令
? 字符串要用单引号并用括号括在里面。如:

》disp(‘text string’) %disp字符串显示命令
text string ? 在单引号里边的字符串可以作为矢量或矩阵的元素。 使用disp命令或输入变量名就可以显示它们表示的字 符串。如: 》a=['this is a';'text string'] ??? All rows in the bracketed expression must have the same number of columns.(元素1含9个字 符,包括空格;元素2含有11个字符;因此只要在元素 1中加入2个空格即可解决问题) 》aa=['this is a ';'text string'] 》disp(aa) aa = aa= this is a this is a 70 text string text string

? 宏是MATLAB语言用在常用命令部分的缩写。它可以被 存储用于建立M文件的一部分。宏命令采用字符串, 并使用eval命令去执行宏命令。下例是采用宏命令计 算阶乘的例子。 》fct=‘prod(1:n)’; %求10的阶乘 》n=10;eval(fct) ans=3628800 六、常用的编程命令(例exp2_8.m) ? pause:停止m文件的执行直至有键按下。pause(n)将 使程序暂停n秒。 ? echo on/off:控制是否在屏幕上显示程序内容。 ? keyboard:停止程序执行,把控制权交给键盘。输入 return并回车后继续程序执行。 ? x=input(‘prompt’):把输入的字符串作为提示符,等待 使用者输入一个响应,然后把它赋值到x。
71

七、关系与逻辑运算符
? MATLAB的运算符有三种类型:算术运算符、关系运算符、逻辑 运算符。它们的处理顺序依次为算术运算符、关系运算符、逻辑 运算符。在前面我们已经介绍例算术运算符,这里我们着重介绍 后两种运算符。

1、关系运算符 假设有:A=[1 2 -1 -5] < 小于 A<B > 大于 A>B <= 小于等于 A<=B >= 大于等于 A>=B == 等于 A=B ~= 不等于 A~=B

B=[0 2 3 1] ans=[0 0 1 1] A<1 ans=[0 0 1 1] ans=[1 0 0 0] A>1 ans=[0 1 0 0] ans=[0 1 1 1] ans=[1 1 0 0] ans=[0 1 0 0] ; A=1 ans=[1 0 0 0] ans=[1 0 1 1] ; A~=1 ans=[0 1 1 1]

72

2、逻辑运算符
注意:在处理逻辑运算时,运算元只有两个值 即0和1,所以如果指定的数为0,MATLAB认 为其为0,而任何数不等于0,则认为是1。 设有:A=[5 -4 0 -0.5] B=[0 1 0 9]

& 与 A&B=[0 1 0 1] A&1=[1 1 0 1] | 或 A|B=[1 1 0 1] ~ 非 ~A=[0 0 1 0]

A|1=[1 1 1 1]

~1=0
73

八、程序流程控制 1、for循环语句
基本格式 for 循环变量=起始值:步长:终止值 循环体 end 步长缺省值为1,可以在正实数或负实数范围内任意指定。 对于正数,循环变量的值大于终止值时,循环结束;对 于负数,循环变量的值小于终止值时,循环结束。循环 结构可以嵌套使用。 ? 书写格式不必太过于拘泥,在Editor编辑器中会自动进行 处理。 (例exp2_11.m)
74

2、while循环语句
基本格式 while 表达式 循环体 end 若表达式为真,则执行循环体的内容,执行后再判 断表达式是否为真,若不为真,则跳出循环体,向下 继续执行。 例exp2_12.m While循环和for循环的区别在于,while循环结构的循环 体被执行的次数不是确定的,而for结构中循环体的执 行次数是确定的。
75

3、if,else,elseif语句 当逻辑表达式的值为真时, (1)if 逻辑表达式 执行该结构中的执行语句, 执行语句 执行完之后继续向下进行; 若为假,则跳过结构中的内 end 容,向下执行。

(2)if 逻辑表达式 (3) if 逻辑表达式1 执行语句1 执行语句1 else elseif 逻辑表达式2 执行语句2 执行语句2 end … end
76

?

if-else的执行方式为:如果逻辑表达式的值 为真,则执行语句1,然后跳过语句2,向下执 行; 如果为假,则执行语句2,然后向下执行。

?

if-elseif的执行方式为:如果逻辑表达式1的 值为真,则执行语句1; 如果为假,则判断逻辑表达式2,如果为真, 则执行语句2,否则向下执行。 例exp2_13.m exp2_13_.m
77

4、switch语句
格式:switch 表达式(%可以是标量或字符 串) case 值1 执行方式:表达式的值 和哪种情况(case)的值相同, 语句1 就执行哪种情况中的语句,如 case 值2 果不同,则执行otherwise中 语句2 的语句。格式中也可以不包括 …. otherwise,这时如果表达式 otherwise 的值与列出的各种情况都不相 语句3 同,则继续向下执行。 end 例exp2_14.m
78

第5节 小结
? 本章所要掌握的是MATLAB语言的基本知识,包括 MATLAB窗口环境的使用;矩阵运算及多项式处理;基 本的绘图命令;程序设计入门。 ? MATLAB的窗口环境是基础的基础,要求熟练掌握语句 命令的输入、变量的使用、基本的数学函数及各种工作 空间与文件管理的命令。要善于运用help命令自学。 ? MATLAB具有强大的矩阵运算能力,要求熟练掌握矩阵 的输入与生成,掌握矩阵的基本运算及操作,区分带点 运算与不带点运算点的不同;掌握多项式的建立与表示 方法及多项式的基本运算。 ? MATLAB同样具有强大的图形处理能力,要求熟练掌握 基本的二维绘图命令。 ? MATLAB具有二种基本的M文件类型,要求掌握它们的 区别及基本结构,熟悉程序流程控制的使用及常用的编 程命令。

79

第一章 小结
? 仿真是对系统进行研究的一种实验方法,基本 原则是相似性原理。 ? 数字仿真具有经济、安全、快捷的特点。 ? 仿真是在模型上进行的,建立系统的模型是关 键。 ? 系统模型可以分为物理模型、数学模型及仿真 模型,也可分为物理仿真和数学仿真两大类。 ? 系统、模型、计算机是仿真的三个基本要素, 建模、仿真实验及结果分析是三项基本内容。 ? MATLAB与SIMULINK是当今广泛为人们采用 的控制系统数字仿真与CAD应用软件。

80

第二章 控制系统数学模型与 模型的相互转换、连接
? 控制系统计算机仿真与辅助设计是建立在控制系统数学 模型的基础之上的。

(1)首先就需要建立系统的数学模型,将数学模型转变为 仿真模型,即数值算法模型。 (2)通过对数学模型的求解分析实现对系统的分析与设计。

81

2.1 连续系统的数学描述
连续系统的数学模型通常可以用三种形式对系统 加以描述:微分方程、传递函数、状态空间表达式。

? 2.1.1 系统的微分方程形式模型
系统的动态特性通常可用高阶微分方程来加以描述, 因此最常用的数学模型就是微分方程,假设连续系统 为单入单出(简称SISO)系统,其输入与输出分别用 u(t)、y(t)加以表示,则描述系统的高阶微分方程为:

82

dn y
n

dt dt dt d n ?1 u d n?2 u c1 ? c2 ? ? ? ? ? cn u n ?1 n?2 dt du

? a1

d n ?1 y
n ?1

? a2

d n?2 y
n?2

? ? ? ? ? a n ?1

dy ? an y ? dt
(2-1)

其初始条件为: y(t0) = y0 ,

? ???? y= y , (t0) 0,

u(t0)=u0 ,

? ? u0) = u (t 0
83

对上式加以整理,可以得到:
n ?1

y ? u

j ?0

? a n?

i ?0 n

c n ?i p i ?
j

p

j

84

2.1.2 系统传递函数形式模型 (一)传递函数模型 对(2-1)式等号两边取拉氏变换,并假设y与u的各阶导数 的初值均为零,则存在: snY(s)+a1sn-1Y(s)+ ??? +an-1sy(s)+anY(s) =

c1sn-1U(s)+c2sn-2U(s)+???+cnU(s)
式中 Y(s) ? 输出 y(t)的拉氏变换。 U(s) ? 输入 u(t)的拉氏变换。 从而(2-1)式所描述的系统的传递函数为:
85

G(s) =

c1 s n ?1 ? c 2 s n ? 2 ? ? ? ? ? c n ?1 s ? c n Y (s) ? n U ( s) s ? a1 s n ?1 ? a 2 s n ? 2 ? ? ? ? ? a n ?1 s ? a n

(2-2)

传递函数是经典控制论描述系统的数学模型之一, 是研究线性系统动态响应和性能的重要手段与方 法。 在MATLAB语言中,可以利用分别定义的传递函 数分子、分母多项式系数向量方便地对其加以描 述。 例如对于(2-2)式所示系统可以分别定义传函的 分子分母多项式系数向量:
86

num=[c1 c2 ??? cn-1 cn]

den=[1 a1 a2 ??? an-1 an]
在这里分子、分母多项式系数向量中,系数均按s的 降幂排列。 MATLAB 5.x中可以用tf( ) 来建立传递函数的系 统模型,其基本格式为:

sys = tf (num,den)
例2-1 已知系统传递函数为: G(s) =

2s ? 9 s 4 ? 3s 3 ? 2s 2 ? 4s ? 6

87

则可以利用MATLAB将系统模型表示出来:

num = [2 9];

den = [1 3 2 4 6];
printsys(num,den,'s')
执行上述语句后得到:

num/den =
2s+9
----------------------------------

s^4 + 3 s^3 + 2 s^2 + 4 s + 6
88

对于MATLAB 5.x可以利用 tf ( )直接建立系统模型:

num = [2 9];
den = [1 3 2 4 6]; model = tf(num,den)
执行后得到:

Transfer function: 2s+9 ----------------------------------s^4 + 3 s^3 + 2 s^2 + 4 s + 6
89

MATLAB的不同版本传递函数的建立有所不同,但 结果是一致的。 对于连续时间系统可以用传递函数对其加以表示, 而对于离散时间系统则应采用脉冲传递函数对其进 行描述。 脉冲传递函数一般可表示为关于z的降幂多项式分式 形式: G(s) =

Y ( z ) c m z ? c m?1 z ? ? ? ? ? c1 z ? c0 ? n U ( z ) a n z ? a n?1 z n?1 ? ? ? ? ? a1 z ? a0
m
90

m ?1

在MATLAB 5.0 中对于离散系统同样可以建立相应的 系统模型,其基本格式为:

num=[cm,cm-1, ??? ,c1,c0];

den=[an,an-1, ??? ,a1,a0];
sys=tf(num,den,T)
其中 T为系统采样周期。 (二)系统的零极点形式模型 系统的传递函数还可表示成另一种形式,即零极 点形式。系统的零极点模型一般可表示为:
91

G(s)=

( s ? z1 )(s ? z 2 ) ? ? ? ( s ? z m ) K ( s ? p1 )(s ? p 2 ) ? ? ? ( s ? p n )

MATLAB可以使用zpk( )函数建立零极点形式的系统 模型,其基本格式为:

sys=zpk([z],[p],[k])
其中,[z],[p],[k]分别为系统的零极点和增益。

92

2.1.3 系统的状态空间表达式

状态方程是研究系统的最为有效的系统数学描述。 通常一个线性定常系统可以表示为以下形式:
(2-6)

第一式由n个一阶微分方程构成,称其为系统的状态方 程表达式,下式由l个线性代数方程组构成,称为系统的 输出方程 。

93

在MATLAB工作空间中建立系统的系数阵,形式如下: A = [a11 a12, ???,a1n;a21 ???,ann]; B = [b1; b2; ??? bn]; C = [c1, c2, ??? cn]; D = d; a22, ???,a2n; ?????,an1,an2

94

根据系统状态方程的系数阵可在MATLAB中建立系统 模型,其基本格式为: sys = ss(A,B,C,D) 上述是5.x中的格式,在MATLAB4.2中其格式为: printsys(A,B,C,D)
对于离散系统,若状态空间表达式可表示成下述形式: X(k+1) = AX(k) + BU(k) Y(k) = CX(k) + DU(k) 在MATLAB 5.x中可建立相应的系统模型,其格式为: sys = ss(A,B,C,D,T) 95 其中,T为系统采样周期

2.2 系统模型的相互转换
2.2.1 控制系统模型向传递函数或零极点增益形式 转换

(一) 状态空间方程向传递函数形式的转换 系统的状态空间方程可表示为:

等效的系统传递函数模型:

num( s ) ? C ( sI ? A) ?1 B ? D G(s) = den( s )
显然变换过程中,求取(sI-A)阵的逆比较困难, 然而 MATLAB有一系列的函数可以完成各种变换, 96

包括进行这种变换的 ss2tf ( ) 函数,其基本格式为: [num,den] = ss2tf ( A, B, C, D, iu ) 利用该函数即可实现将状态空间方程转换为传递 函数的形式, iu用于指定变换所使用的输入量.
还可以采用下述方式进行: G1 = ss ( A, B, C, D ); G2 = tf ( G1 ) 可以证明,由给定的状态方程模型转换为传递函数 形式其结果是唯一的。

97

? 例 2-3 已知连续系统 的系数矩阵是:

? A= ?2 0 0? ?0 4 1 ? ? ?
?0 0 4? ? ?

B=

?1 ? ?0 ? ? ? ?1 ? ? ?

C=[1 1 0]

D=0

求取该系统相应的传递函数模型.

98

解 应用MATLAB的ss2tf( )函数可以方便实现这种转换. A = [2 0 0;0 4 1;0 0 4]; B = [1 0 1]'; C = [1 1 0]; D = 0; [num,den] = ss2tf(A,B,C,D,1) G = tf (num,den) 利用G=tf(num,den) 将转换后的传递函数分子分母多项 式系数向量构造成传函系统模型,下面给出了运行结果:
99

num = 0 1.0000 -7.0000 14.0000 den = 1 -10 32 -32 Transfer function: s^2 - 7 s + 14 -------------------------------s^3 - 10 s^2 + 32 s – 32 考虑对于多输入系统,应用 ss2tf( )函数可以进 行指定要求的模型转换。

100

(二) 模型向零极点形式的转换
MATLAB中也提供了实现将各类系统模型转换为零 极点形式模型的函数。其基本格式为: [z,p,k] = ss2zp (A,B,C,D,IU) [z,p,k] = tf2zp (num,den) Gzp = zpk (sys) 第一式是将以状态方程形式转换为零极点模型形 式;第二式是将传递函数形式转换为零极点形式; 第三式可将非零极点形式的模型转换为零极点系统 模型。
101

2.2.2 系统模型向状态方程形式的转换
利用MATLAB函数实现系统模型向状态方程的转换。 基本格式为:
[A, B, C, D] = tf2ss ( num, den ) [A, B, C, D] = zp2ss ( z, p, k ) syss = ss(sys)

例 2-6 已知系统传递函数为:
18s ? 36 G(s) = 3 s ? 40.4s 2 ? 391s ? 150

应用MATLAB的模型转换函数将其转换为状态方 程形式的模型。
102

解: 应用tf2ss( ) 应注意由于其状态变量选取的 不同,转换结果是不唯一的。 应用tf2ss( ) 函数的转换程序与结果: num = [18,36]; den = [1 40.4 391 150]; [A,B,C,D] = tf2ss(num,den)

103

A= -40.4000 1.0000 0 -391.0000 -150.0000 0 0 1.0000 0

B=
1 0 0

104

C= D=

0 0

18

36

利用 [A,B,C,D] = zp2ss (z, p, k)可以将零极点 形式给出的模型转换为状态方程。对于 sys = ss (sys),可以将任意LTI系统模型转换为状态方程。 2.2.3 连续和离散系统之间的转换
在采样控制系统中,对于控制器的设计经常采 用模拟化设计方法,MATLAB 对于离散化转换可采 用相应的函数加以进行。

105

例 2-8 已知系统传递函数为: G(s)
s?2 = 0.5s 3 ? 8.5s 2 ? s ? 15

采样周期 T=0.2 秒,试将其进行离散化处理。

解: 可以先将其转换为相应的状态方程形式,然后 再利用ss( )函数进行离散化处理。 num = [1 2]; den = [0.5 8.5 1 15]; T = 0.2; [A,B,C,D] = tf2ss(num,den); G = ss(A,B,C,D,T); Sys = tf(G)
106

运行上述程序得到: Transfer function:
2z+4 ----------------------------z^3 + 17 z^2 + 2 z + 30 本例采用了 ss (A, B, C, D, T) 函数进行离散化处 理,还可以直接应用 tf (num, den, T) 函数直接进 行转换。 如果对离散化处理提出具体转换方式要求,则可 以采用 c2d ( )或 c2dm ( ) 函数进行,其基本格 式为: sysd = c2d ( sysc, T, method )
107

其中 sysyc ? 连续系统模型 T ? 采样周期 method - 指定转换方式。 ‘zoh’表示采用零阶保持器; ‘foh’ 采用三角形近似; ‘tustin’ 采用双线性变换;

108

例如对于上式若要求采用双线性变换方 法进行转换,则可以编写如下程序: num=[1 2];den=[0.5 8.5 1 15];T=0.2; [A,B,C,D]=tf2ss(num,den); G=ss(A,B,C,D); sysd=c2d(G,T,'tustin') Sys=tf(sysd)

运行结果为:

109

a= x1 x2 x3 b= x1 x2 x3 c= x1 -0.27273 0.07273 0.00727 x2 -0.36364 0.96364 0.19636 x3 -2.18182 -0.21818 0.97818

u1 0.36364 0.03636 0.00364
x1 0.01745 x2 0.47127 x3 0.74764
110

y1

d=
u1 y1 0.00873 Sampling time: 0.2 Discrete-time system. Transfer function: 0.008727 z^3 + 0.01164 z^2 - 0.002909 z - 0.005818 ----------------------------------------------------------------z^3 - 1.669 z^2 + 0.4982 z + 0.2582 若采用零阶保持器进行转换,其结果为: 0.01941 z^2 - 0.005024 z - 0.005344 ----------------------------------------------z^3 - 1.961 z^2 + 1.062 z - 0.03337

111

对于离散系统,如果需要对采样周期进行修改, 可以通过 MATLAB d2d ( ) 函数加以进行,其格 式为:

sys = d2d ( sysd, T )
其中 sysd 为原离散系统,T为修改后的系统采样 时间。

112

2.3 系统状态方程的变换与实现
2.3.1状态方程的相似变换
一个系统的传递函数,可以存在着众多的状 态方程实现,因此也就存在着它们之间的相互变 换,MATLAB中的 ss2ss ( ) 函数就可以实现对 系统状态方程的相似变换。其基本格式为: GT = ss2ss ( G, T ) 其中G为原系统的状态方程模型,T为非奇异变换 阵。 该变换实现相似变换:
113

z=Tx

其中

At=TAT-1, Bt=TBT-1, Ct=TCT-1, Dt=TDT-1
例 2-10 已知系统状态方程为:

?10 ? 4 2? A= ? 8 0 0? ? ? ? 0 2 0? ? ?

?2? B= ?0 ? ? ? ?0 ? ? ?

C=[0.5 -0.4375 0.4375], D=0
114

试应用ss2ss( ) 函数进行状态方程的相似变换。

解:由于变换阵T可以任意选择,只要保证其非奇异即 可,在此选择单位反对角阵作为变换阵,下面给出了 变换方法与结果。 A=[10 -4 2;8 0 0;0 2 0]; B=[2 0 0]'; C=[0.5 -0.4375 0.4375]; D=0; G1=ss(A,B,C,D) T=fliplr(eye(3)); GT=ss2ss(G1,T) 得到结果:
115

a=

x1
x1 x2 x3 b= u1 x1 x2 x3 c= 0 0 2.00000 0 0 2.00000

x2 x3 2.00000 0 0 8.00000 -4.00000 10.00000

x1
y1 d=

x2 x3 0.43750 -0.43750 0.50000

u1
y1 0
116

2.3.2 规范型状态方程的实现
MATLAB中的规范实现函数cannon( ),可以 进行LTI系统模型sys的规范状态空间表达式的实现。 其基本格式如下: G1=cannon(sys, type) 同时,还具有可以返回状态变换阵的形式。 [G1,T]=cannon(sys, type) 式中,sys表示原系统状态方程模型,字串 type 确 定规范形式的类型,它可以是模态(modal)规范 型(约当标准型),也可以是伴随矩阵(companion) 形式。T是状态变换阵返回变量,满足 z=Tx,关系, 其中要求sys为状态空间模型。
117

2.3.3 系统的均衡实现
考查系统 ,其系统的系数阵分别为:

?? 10 0 ? A= ? 0 - 25? ? ?

? 10 ? 5 -5 ,B =? 5 ? C= 10 10 ? 10 ?
?5

?

?,

D=0 在该系统中系数阵的各个元素的值相差极为悬殊, 显然对这类问题直接进行求解必然会在数值运算过 程中由于舍入的处理而带来严重的误差。 在MATLAB中提供了可以进行系统均衡变换的函 数,以消除计算中舍入对其造成的严重误差影响。 其均衡变换函数的具体格式为:
118

[Ab,Bb,Cb,G,T]=balreal(A,B,C)
其中 T为均衡变换阵,G为均衡系统的Gram阵。并 且满足下述变换关系:

Ab=T-1AT,Bb=T-1B,Cb=CT,Db=D
同时系统的状态变量也满足:Xb=TX 关系。 以前面给出的不均衡系统作为示例,利用

[Ab,Bb,Cb,G,T]=balreal(A,B,C)
得到相应的结果: Ab = -14.01374575757560 -6.64048427128458 -6.64048427128458 -20.98625424242441
119

Bb =
1.37309792179507 0.33852931507043 Cb = 1.37309792179507 0.33852931507043 G= 0.06726959142308 0.00273040857692 T= 1.0e+004 * 0.00000000085581 -0.00000000051728 5.17284303362319 8.55813618432752 由上述结果可以清楚地看出,系统的系数阵实现了预 期的系数均衡处理,注意只有稳定的系统才能进行均 120 衡变换。

2.3.4 系统的降阶实现
控制系统模型的降阶技术是简化系统分析的重要手 段,其目的是由相对低阶的模型尽可能好地近似一个高 阶原系统,使高阶模型可以按照低阶的仿真与设计的方 法加以进行。

在MATLAB中为用户提供了实现系统降阶处理的专 用函数,如modred ( )其基本格式如下:

RSYS = modred(SYS,ELIM) RSYS = modred(SYS,ELIM,‘mdc’) RSYS = modred(SYS,ELIM,‘del’)
121

其中,ELIM为待消去的状态;‘mdc’表示在降阶中保 持增益的匹配;‘del’表示在降阶中不能保证增益的匹 配。 例 2-12 已知系统的传递函数为: 180

G( s ) ?

s 4 ? 20s 3 ? 136s 2 ? 380s ? 343

应用modred( )函数进行降阶处理,保留前两个状态,降 为二阶系统。 解: 先构造modred所需要的函数,再进行降阶处理。

num=180; den=[1 20 136 380 343] [a,b,c,d]=tf2ss(num,den)

122

sys=ss(a,b,c,d) sysm=modred(sys,3:4,'del')
执行上述语句得到系统降阶后的结果:

a=

x1 x2 b=
x1 x2

x1 x2 -20.00000 -136.00000 1.00000 0 u1 1.00000 0

123

c= y1 d= y1 u1 0

x1 0

x2 0

显然在直接利用modred函数具有一定的盲目性,为 此往往将balreal函数与modred函数相结合加以使用。 由balreal函数先进行均衡变换,依据gram阵确定对 系统影响较小的状态,再应用modred函数求出降阶 后的系统。
124

例2-13 已知系统的系数阵为:
? ?3 ? ? 0 .5 A?? ? ? 1 .5 ? ? ? 1 .5 1 ?1 1 2 0 1 ?2 1 ?1 ? ?1 ? ? 0 ? ? ?4 ?

?1 ? ?0 ? B?? ? ?0 ? ? ? ?0 ?

C ? ?1 0 ? 1 0?

在尽可能保持系统基本特性的情况下进行 降阶处理。 解:对系统进行均衡变换。
125

A=[-3 1 0 -1;-0.5 -1 1 -1;-1.5 1 -2 0;-1.5 2 1 -4] B=[1 0 0 0]' C=[1 0 -1 0] [Ab,Bb,Cb,G,T]=balreal(A,B,C) 得到均衡变换后的系统模型,以及变换 阵T、gram阵。 Ab = -1.2829 0.4033 -0.1900 0.0359 0.4033 -1.7748 1.5444 -0.3144 -0.1900 1.5444 -5.9201 2.8156 -0.0359 0.3144 -2.8156 -1.0223
126

? Bb = ? 0.9849 ? -0.1577 ? 0.0730 ? 0.0138 ? Cb = ? 0.9849 -0.1577 0.0730 -0.0138 ? G= ? 0.3781 ? 0.0070 ? 0.0005 ? 0.0001
127

T= 0.6403 -1.6938 1.8373 -2.3068 -0.1591 -1.0219 -0.0109 -0.2640 -0.3446 -1.5361 1.7643 -2.2930 -0.3092 -1.3170 1.5967 -1.4304 由gram阵可以看出,变换后的模型中, 第3、4状态变量对系统的作用较小,因此 可以利用modred函数进行降阶处理,保留 第1、2状态。键入下述命令:
128

sys=ss(Ab,Bb,Cb,Db) sysr=modred(sys,3:4,'mdc') 得到降阶模型:
a= x1 -1.27804 0.36341 x2 0.36341 -1.44664

x1 x2
b= x1 x2

u1 0.98303 -0.14235
129

c=

y1
d= y1

x1 0.98303

x2 -0.14235

u1 0.00071

130

在MATLAB中还给出最小实现函数 minreal( ),它的基 本格式为: [Am,Bm,Cm,Dm] = minreal(A,B,C,D) [numm,denm] = minreal(num,den) 该函数消去了不必要的状态,从而得到系统的最小实现。

131

2.4 控制系统模型的建立与典型连接
2.4.1 基本系统模型的建立 MATLAB 提供了一些基本系统或模型建立的函数,下 面对其分别加以介绍。 (一)二阶系统的生成 在控制系统当中,二阶系统占有相当的比例, 因此 研究与讨论二阶系统具有重要的意义。 MATLAB 提供了二阶系统生成函数 ord2( ),其基本格 式为: [A,B,C,D] = ord2(Wn, z)
132

式中 Wn为自然角频率, z为阻尼因子。 自然角频率与阻尼因子的意义如下式中表示 :

G(s) ?

1 s ? 2?? n s ? ?
2 2 n

同时MATLAB也提供了生成以传递函数形式描 述的二阶系统的ord2( )函数,其格式为:

[num,den] = ord2(Wn,z)
返回变量为传递函数的分子分母多项式系数向 量。 133

例 2-14 分别生成以状态方程形式和传递函数形式描 述的二阶系统,要求=0.8, =2.2 弧度/秒。 解:直接应用ord2( ) 函数进行。

Wn=2.2; z=0.8; [A,B,C,D]=ord2(Wn,z)
得到:

A= 0 1.0000 -4.8400 -3.5200
134

B= 0 1 C= 1 0 D= 0 求取传递函数形式: Wn=2.2; z=0.8; [num,den]=ord2(Wn,z); G=tf(num,den)
135

Transfer function:
1 ------------------------s^2 + 3.52 s + 4.84

(二)随机n阶系统的模型建立
MATLAB为用户提供了建立随机n阶系统模型的函数, 其基本格式为: [num,den]=rmodel(n) [num,den]=rmodel(n,p) [A,B,C,D]=rmodel(n) [A,B,C,D]=rmodel(n,p,m) 136

[A,B,C,D] = rmodel(n,p,m) [num,den] = drmodel(n) [num,den] = drmodel(n,p) [A,B,C,D] = drmodel(n) [A,B,C,D] = drmodel(n,p,m) 其中 [num,den]=rmodel(n) 可以随机生成n阶稳定传递 函数模型。 [num,den]=rmodel(n,p) 可以随机生成单入p出的,n 阶稳定传递函数模型。 [A,B,C,D]=rmodel(n) 可以随机生成n阶稳定SISO状 态方程模型。 [A,B,C,D]=rmodel(n,p,m) 可以随机生成n阶稳定p出 m入状态空间模型。 137

drmodel(n) 函数将生成相应的离散模型。

如果欲生成一随机三阶二输入单输出状态空间表 达式,则可应用[A,B,C,D]=rmodel(n,p,m)。 运行:[A,B,C,D]=rmodel(3,1,2) 得到:

A=
-0.9688 -4.0932 0.6941 4.0638 -1.0010 -0.4993 -0.8493 0.0983 -2.0227 B=

0 0.7143 1.6236

0 0.8580 1.2540

138

C=
-1.5937 D= 0 0.6900 0 0.5711

(三) 系统模型的重构 1. 子系统选取与删除 MATLAB提供了一个从大系统中选择子系统的函数 ssselect,其基本格式为: [Ae,Be,Ce,De] = ssselect(A,B,C,D,inputs,outputs) 式中inputs指定系统输入,outputs指定系统输出。
139

[Ae,Be,Ce,De]=ssselect(A,B,C,D,inputs,outputs,state s) 利用指定的输入、输出及状态向量构建子系统。
例 2-15 已知系统的系数阵为:

A=

? 1 5 3 4? ? 7 16 7 8 ? ? ? ? 3 0 22 4? ? ? ? 42 25 7 9 ?

B

? 6 41 3 ? ?12 25 4 ? =? 3 8 26? ? ? ? ? ? 6 47 2 ?

?1 5 9 3 ? C =? ? ? 2 3 5 2?

?22 D=? 23 ?

46 12? 54 2 ? ?

利用ssselect ( ) 函数在原系统的基础上构造新系统,保 留1、3输入信号,1、2输出以及1、2、4状态变量。
140

解:根据题意要求可以编写相应的程序为: A = [1 5 3 4;7 16 7 8;3 0 22 4;42 25 7 9]; B = [6 41 3;12 25 4;3 8 26;6 47 2]; C = [1 5 9 3;2 3 5 2]; D = [22 46 12;23 54 2]; ; inputs = [1 3]; outputs = [1 2]; states = [1 2 4]; [A1,B1,C1,D1] = ssselect (A, B, C, D, inputs,outputs, states) 运行上述语句得到新构建的系统模型:
141

A1 = 1 5 4 7 16 8 42 25 9 B1 =

6 12 6 C1 =
1 2

3 4 2
5 3 12 2 3 2

D1 =
22 23
142

ssselect( )函数可应用于连续系统,与离散系统。 MATLAB还提供了与ssselect( )函数相反操作的一类 函数它们是 [Ar,Br,Cr,Dr] = ssdelete(A,B,C,D,inputs,outputs) [Ar,Br,Cr,Dr]= ssdelete(A,B,C,D,inputs,outputs,states) 该函数可从状态空间系统(A,B,C,D)中删除由 inputs、outputs和states指定的输入、输出与状态。 2.状态的增广 在对系统进行分析研究时,往往需要对状态在系统 (输出)中加以增广。
143

状态增广有时具有十分重要的实际用途,在MATLAB 中特别提供了一种状态增广函数 augstate( )。其基本 格式如下: asys = augstate(sys) augstate( ) 函数将状态增广到状态空间系统中去,构 造成新系统:
例 2-16 已知系统状态方程为:

144

A=

?10 ? 4 2? ?8 0 0? ? ? ?0 2 0? ? ?

B=

?2? ?0 ? ? ? ?0 ? ? ?

C = [0.5 -0.4375 0.4375]

D=1 将所有的状态增广到系统输出中。 解:应用 augstate( ) 函数进行输出状态增广。 A = [10 -4 2;8 0 0 ;0 2 0]; B = [2 0 0]'; C = [0.5 -0.4375 0.4375]; D = 1; SYS = SS(A,B,C,D) ASYS = augstate(SYS)
145

执行上述语句得到新构建的系统: Continuous-time system. a= x1 x2 x1 10.00000 -4.00000 x2 8.00000 0 x3 0 2.00000 b= u1 x1 2.00000 x2 0 x3 0

x3 2.00000 0 0

146

c= y1 y2 y3 y4 d= 0.50000 1.00000 0 0

x1 x2 x3 -0.43750 0.43750 0 0 1.00000 0 0 1.00000

u1 y1 1.00000 y2 0 y3 0 y4 0 显然在新系统中,状态方程的系数阵没有变化,而 147

输出方程中的C与D阵都出现了增广,满足将全部状态 添加至系统输出的要求。 2.4.2 系统组合与连接

系统组合是将两个或多个子系统按一定方式加以 连接形成新的系统。连接组合方式主要有:串联、并联、 反馈等形式。MATLAB提供了进行这类组合连接的相 关函数。 (一)模型串联 在MATLAB中提供了模型串联连接函数series( ), 两个 线性模型的串联,其基本格式为:

sys = series(sys1,sys2)

148

上述串联函数实现了sys1和sys2相串联形成新系统

sys= sys1*sys2
如果相串联的两个系统sys1(s)、sys2(s)的状态方程系 数阵分别为:(A1,B1,C1,D1),(A2,B2, C2,D2),则串联后整个系统的系数阵将变为:

? A1 A=? ? B2 C1
D=D1D2

0? ? A2 ?

? B1 ? B= ? ? ? B2 D1 ?

C=[D2C1 C2],

对于多入多出系统,上述串联函数的形式为:

sys=series(sys1,sys2,outputs1,inputs2)
149

上述函数将实现,由outputs1所指定的系统1输出端连 接到由inputs2所指定的系统2输入端,inputs2与 outputs1分别为系统2和系统1的输入、输出向量。 (二)模型并联 在MATLAB中提供了并联连接函数 parallel( ), 其基本 格式为:

sys=parallel(sys1,sys2,in1,in2,out1,out2)
该函数将两个线性系统以并连方式加以连接,

150

当sys1、sys2为SISO系统时该函数简化为: sys=parallel(sys1,sys2) Sys1与sys2在共同的输入信号作用下,将产生两 个输出信号,而并联系统的输出就是两个系统输 出之和。 若用传递函数对系统加以描述,系统总的传递函 数 G(s)=G1(s)+G2(s)。若考虑:

num1( s ) G1(s) = , den1( s )

num2( s ) G2(s) = den 2( s )

151

则系统总传函为:

num1( s)den2( s) ? num2( s)den1( s) G(s) = den1( s)den2( s)

152

例 2-17 已知两个线性系统

12s ? 4 G1(s) = 2 s ? 5s ? 2

s?6 , G2(s) = s 2 ? 7s ? 1

应用parallel( )函数进行系统并联连接。

num1 = [12 4]; den1 = [1 5 2]; num2 = [1 6]; den2 = [1 7 1]; sys1 = tf(num1,den1); sys2 = tf(num2,den2); sys = parallel(sys1,sys2)

153

运行上述语句得到连接后的总系统传函: Transfer function:
13 s^3 + 99 s^2 + 72 s + 16 ---------------------------------------s^4 + 12 s^3 + 38 s^2 + 19 s + 2

如果系统用状态方程的形式给出,则并联后的系统模 型为:

? ? x1 ? ? A1 ?x ? ? ? 0 ? ?2 ? ?

? B1 ? 0 ? ? x1 ? ? ? ? + ?B ? A2 ? ? x 2 ? ? 2?
154

? x1 ? y = [C1, C2] ? ? + (D1 + D2) u ?x?
若已知两个线性系统,应用parallel( 实现系统的连接: )函数即可方便地

A1 = [1 2 ; 3 4] ; A2 = [2 4 ; 3 5] ;

B1 = [2 6] ' ; B2 = [2 5] ' ;
C1 = [1 5] ; C2 = [2 3] ; D1 = 4 ; D2 = 1; sys1 = ss (A1,B1,C1,D1); sys2 = ss (A2,B2,C2,D2);

sys = parallel (sys1,sys2)

155

运行后得到结果: a= x1 x1 1.00000 x2 3.00000 x3 0 x4 0 b= u1 x1 2.00000 x2 6.00000 x3 2.00000 x4 5.00000

x2 x3 x4 2.00000 0 0 4.00000 0 0 0 2.00000 4.00000 0 3.00000 5.00000

156

c=

x1
x4 y1 3.00000 d= y1 1.00000

x2
5.00000

x3
2.00000

u1 5.00000

(三)反馈连接 反馈系统是控制系统中最为重要与常见的一类系 统,MATLAB提供了构造反馈系统的函数feedback( ), 其一般格式为:

sys = feedback(sys1,sys2,sign)

157

执行该语句将实现两个系统的反馈连接,sign缺省时 即为负反馈,sign=1时为正反馈。 如果sys1与sys2用传递函数描述,则反馈系统的传函 为: SYS1 ( s ) SYS(s) = 1 ? SYS1 ( s ) SYS 2 ( s ) 若sys1与sys2分别以状态方程的形式给出(A1,B1, C1,D1),(A2,B2,C2,D2),则反馈系统的状态 方程系数阵可表示为:

? B1C 2W ? ? A1 ? B1 D2 C1W ? B1W ? A =? ? ,B= ? B D W ? A2 ? B2 D1C 2W ? ? B2 C1Q ? 2 1 ?


158

C=[QC1 - D1C2W], D=D1W, 式中 W=(I+D1D2)-1 ,Q=I - D1WD2 对于MIMO系统可建立复杂的反馈系统,其格式为: sys=feedback(sys1,sys2,feedin,feedout,sign) 式中feedin为sys1的输入向量,feedout指定sys1的那 些输出端用于反馈。最终实现的反馈系统与 sys1 具有 相同的输入、输出端。 例 2-18 已知前向环节与反馈环节的状态空间表达式的 系数阵分别为:
159

?1 0? A1= ? , 0 1? ? ?

?1 1? ?1 3? B1= ?0 1? , C1= ? ? , D1 = ? ? ? 2 0?
?1 ? ?0 ? ? ?

? ? 2 0? A2= ? 1 0? , B2 = ? ?

?1 0 ? ? 2 5? ? ?

, C2=

[0 1] , D2 = 0 ,

试将前向环节的输入1和输出2,与反馈环节构成负反馈 系统。 解:显然环节1为一个多入多出系统, 可用反馈函数 feedback( )加以实现。其相应的程序给出如下:

A1 = [1 0; 0 1]; B1 = [1 1; 0 1]; C1 = [1 3; 2 0]; D1 = [1 0; 2 5];
160

A2 = [-2 0; 1 0]; B2=[1 0]'; C2 = [0 1]; D2=0; sys1 = ss (A1,B1,C1,D1); sys2 = ss (A2,B2,C2,D2); feedin = 1; feedout = 2; sign = -1; sys = feedback (sys1,sys2,feedin,feedout,sign) 运行上述程序得到相应的结果:

161

a= x1 x2 x3 x4 b= x1 x2 x3 x4 x1 x2 x3 1.00000 0 0 0 1.00000 0 2.00000 0 -2.00000 0 0 1.00000 u1 1.00000 0 2.00000 0 u2 1.00000 1.00000 5.00000 0 x4 -1.00000 0 -2.00000 0

162

c=
y1 y2 d= x1 1.00000 2.00000 u1 1.00000 2.00000 x2 3.00000 0 u2 0 5.00000 x3 0 0 x4 -1.00000 -2.00000

y1 y2

可以看出新的反馈系统,其输入输出端数 量没有改变,而状态发生了变化。
163

(三) 系统的复杂连接 在系统分析与研究的实际问题当中,往往会遇到 一些十分复杂的控制系统,需要更加有力的手段和方 法。 MATLAB为进行这种工作提供了必要的选择。如 可利用connect( )函数以及nblocks – blkbuild 所构造 的程序模块,实现由众多环节组合连接成复杂系统。

164

第三章 控制系统数字仿真
3. 1基本数值积分方法
? 连续系统的动态特性一般可由一微分方程或一组一 阶微分方程加以描述。因此对系统进行计算机仿真就 需要研究如何利用计算机对微分方程进行求解。 ? 目前所采用的方法就是应用数值积分方法对微分方 程求数值解,其中欧拉法、梯形法、龙哥库塔法等数值 积分方法是几种基本的方法。

165

3. 1. 1欧拉法(Euler)
欧拉法是一种最简单的数值积分方法,但其精 度差,实际应用很少,然而由于其算法简单且具 有明显的几何意义,因此在讨论数值积分方法时 通常先介绍欧拉法。

假设一阶微分方程为:

166

dy ? f (t , y) dt
初始条件为:

(3-1)

y(t0 ) ? y0
欧拉法的基本思想就是将式(3-1)中的积分曲线解 用直线段所组成的折线加以近似。由

dy dt


y(t ? ?t ) ? y(t ) = li m ?t ?t ?0

y(t ? ?t ) ? y(t ) li m ?t ?t ?0

= f [ t , y (t) ] ( 3-2 )
167

上式可近似表示为:

y (t ? ?t ) ? y (t ) ?t

?

f [t , y(t)]

(3-3)

若令?t = h , t = t 0 则得到, y ( t 0 + h) y ( t 0 ) + h f [ t 0 , y ( t 0 )] (3-4) 其中 h称为计算步长。

?

168

设初始给定点为 P0 , 可求出在 t = t1 = t 0 + h 上的近似值 y2 = y1 + h f ( t1 (3-6) 根据类似的思想不难推出:

,

y1 )

yk ?1 = yk + h f ( t k , yk )

(3-7)

169

170

(3-7)式就称为欧拉法递推公式。
该方法简单,计算量小,由前一点即可推出后一 点值,属于单步法。该算法精度低,适当减小计算步 长h有助于提高计算精度。 3. 1. 2 梯形法 1. 欧拉法误差原因的分析 根据欧拉法的分析可知该方法简单但较粗糙, 精度差。 考虑对于微分方程

171

dy ? f (t , y ) dt

y(t0 ) ? y0
存在精确解 y(t)为: 当 t = t1时 y ( t1) = y(t) = y0 +
t ?t 0

f (t , y ) dt

y0 + y0

t1

? f (t , y ) dt (3-8)

由前述推导可知 y ( t1 ) ?

t0

+ ( t1 - t 0 ) f (t 0 , y0 ) (3-9)

172

对照图分析
对照(3-8)与(3-9)式得到,在欧拉法中存在

t0

? f (t , y )dt ? f (t0 , y0 )(t1 ? t0 )
t1

用 f ( t 0 , y0 ) 代替变量, 由矩形近似曲边梯形, 因 此造成欧拉法误差较大。

173

2. 改进的欧拉法-梯形法 进一步改进欧拉法的方法就是利用梯形面积代替 矩形面积以逼近曲边梯形面积即
t1

? f (t , y )dt

t0

?

h ( f 0 ? f1 ) 2

推而广之,得到梯形法公式:

h yk ?1 = yk + 2 [ f (tk , yk ) ? f (tk ?1, yk ?1)]
由上式可以看出,在计算时,在等式右边存在有 yk ?1 , 因此无法求出 yk ?1 , 在这种情况下,

174

一般采用欧拉法进行一次迭代,求出 yk ?1 的预估值,再 将其代入梯形公式计算出 yk ?1 的真值,从而构成了预 估校正梯形公式如下:
(3-11)
0 y k ?1

=

yk

+ h f (t k ,

yk )

h 0 yk ?1 = yk + [ f (t k , y k ) ? f (t k ?1 , y k ?1 )] 2
(3-12)

(3-11)式称为预报式, (3-12)式称为校正式。 3. 1. 3 龙格-库塔(Runge-kutta)法 考虑如下一阶微分方程

dy ? f (t , y ) dt

175

y(t0 ) ? y0
设 t1 = t 0 + h , 这里h为计算步长,在 t1 时刻的解 为 y1 = y (t 0 + h),可在 t 0 附近展成泰勒级数,取 展开式的前三项得到:

y1 = y0 + h d y
dt

t ?t 0

= y0

1 2d y ? h 2 2 d t2

2

t ?t 0
t ?t0 y ? y0

+ h f ( t 0 , y0 )

h ?f ?f ( ? f ) + 2 ?t ?y

(3-13) 假设上式的解可表示为:
176

y1 = y0 + ( a1K1 ? a2 K 2 )h K1 ? f (t0 , y0 ) K 2 ? f (t0 ? b1h, y0 ? b2 K1h)
展开成泰勒级数,并取前三项得到:

(3-14) (3-15)

(3-16)

f ( t 0 , y0 ) + b1h ?f ? b2 K1h ?f (3-17) ?t ?y 将(3-15), (3-17) 式代入(3-14) 式中, 经整理得到:

K2 ?

?f ?f y1 ? y0 ? a1hf (t0 , y0 ) ? a2 h[ f (t0 , y0 ) ? b1h ? b2 hK1 ] ?t ?y
(3-18)
177

将 (3-18) 式与 (3-13) 式相比较得到:

? a1 ? a 2 ? 1 ? ? 2 a 2b1 ? 1 ?2 a b ? 1 ? 2 2 1 若取 b1 ? b2 ? 1 , 则得到 a1 ? a2 ? 将其代入(32 14)、(3-15)、(3-16) 式得到:
h ? ? y1 ? y0 ? 2 ( K1 ? K 2 ) ? ? K1 ? f (t0 , y0 ) ? K ? f (t ? h, y ? K h) 0 0 1 ? 2 ?

(3-19)

将 (3-19) 式写成递推形式:
178

由于(3-20)式只取了级数展开式的前三项 , 略去了 h 3 以上的高阶无穷小 , 其截断误差正比于h , 故这种方 法称之为二阶龙格-库塔法。 二阶龙格-库塔法的计算精度是不高的, 为此在实 际应用中往往采用精度更高的四阶龙格-库塔法,它的 截断误差正比于 h 5 , 四阶龙格-库塔法的递推表达式:
179

h ? ? y k ?1 ? y k ? 2 ( K 1 ? K 2 ) ? ? K 1 ? f (t k , y k ) ? K ? f (t ? h, y ? K h) k k 1 ? 2 ?

(3-20)
2

(3-20) 式可以清楚地看出,与前面导出的梯形算法 是相同的 , 因而可以知道梯形法的精度也是二阶的。 当对级数展开式只取前两项时可以得到 yk ?1 = yk + h f ( t k , yk ) , 这就是欧拉法递推式。通过对龙格一 库塔法的分析与推导,我们可以将Euler法、
180

h ? ? yk ?1 ? yk ? 6 ( K1 ? 2 K 2 ? 2 K 3 ? K 4 ) ? ? K1 ? f (t k , yk ) ? h h ? ? K 2 ? f (t k ? , yk ? K1 ) 2 2 ? h h ? ? K 3 ? f (t k ? 2 , yk ? 2 K 2 ) ? (3-21) ? K 4 ? f (t k ? h, yk ? hK 3 ) ?

梯形法与Runge-Kutta法这三种数值积分方法统一了 起来,它们都可归于龙格-库塔法只不过是阶次不同 而已。 从理论上讲,可以构造任意阶的计算方法,但截 断误差的阶数与所计算函数数值f的次数并非是等值 关系,对于四阶以上Runge-Kutta(RK)法所需计算 的f次数将高于误差阶数,这将大大增加了计算工作 量,因而选择这种高阶RK法是不适当的。 实际上,对于工程应用四阶Runge-Kutta(RK4) 法已完全能够满足仿真精度要求,所以RK4在实际应 用中具有重要的地位。

181

3. 1. 4 四阶Runge-Kutta法的向量表示 实际中系统往往是高阶的,因而描述这类系统应是 一组一阶微分方程或是状态方程,所以介绍求解状态 方程的RK4矢量式是完全必要的。考虑一n阶系统,其 状态方程的一般描述为:

? Y ? F (t ,Y ) Y (t0 ) ? Y0
四阶R-K法的向量表示为:

182

h Yk ?1 ? Yk ? ( k1 ? 2 k 2 ? 2 k3 ? k 4 ) 6 k1 ? F (t m ,Yk ) h h k 2 ? F (t m ? ,Ym ? k1 ) 2 2 h h k3 ? F (t m ? ,Ym ? k 2 ) 2 2 k 4 ? F (t m ? h ,Ym ? h k3 )

(3-22)

其中 k j ? (k1j , k2 j ,? , knj ), kij (i=1,2,…,n; j=1,2,3,4) 是微分方程组中第 i 个方程的第 j 个RK系数. 有时考虑 到计算方便, (3-22)式又可将其展开为下述形式:
183

yi ,k ?1 ? yi ,k

k i1 ? f i (t k , y1k , y 2k ,... y nk )

h ? (k i ,1 ? 2k i , 2 ? 2k i ,3 ? k i , 4 ) 6

h h h h k i 2 ? f i (t k ? , y1k ? k11 , y 2k ? k 21 ,... y nk ? k n1 ) 2 2 2 2 h h h h k i 3 ? f (t k ? , y1k ? k12 , y 2k ? k 22 ,..., y nk ? k n 2 ) 2 2 2 2

h h h k i 4 ? f (t k ? h, y1k ? k13 , y 2k ? k 23 ,..., y nk ? k n3 ) 2 2 2

其中 i = 1,2,…n ; n为系统阶数 , 下标k为递推下标 。
184

3. 1. 5 阿达母斯法(Adams) 数值积分算法的特点就是,在计算k+1时刻的值 时,只用到第 k 时刻 y k 和 f k 的值。 在逐步递推过程中,计算 yk ?1 之前已经获得了一系 列的近似值 y0 , y1 … y k 以及 f 0 , f1 … f k等,如果 能够充分利用历史上一些数据, 则可望既可加快仿真速 度又能获得较高的仿真精度,这就是构造多步法的基 本思路。 Adams 法就是利用这些信息来实现对 yk ?1 的计 算, 所以在加快仿真速度的同时提高了仿真精度。在这 里不再进行公式的推导,直接给出Adams显式表达式。
185

h = yk + ( 55 f k - 59 f k ?1 + 37 f k ? 2 24 9 f k ?3 )

yk ?1

(3-23) (3-23)即为Adams显式表达式 , 其截断误差为:

Rk

251 5 (5) = h y (? ), 720

? ? (t k -1 , t k )

186

(3-24) 给出了Adams隐式表达式。

yk ?1 = yk
+ fk ?2 ) 其截断误差为:

h + ( 9 f k ?1 +19 f k - 5 f k ?1 24
(3-24)

Rk = 19 h 5 y (5) (? ), 120

? ? (t k -1 , t k )
187

隐式Adams法的优点主要在于精度较高,然而它无 法实现自身进行递推,(3-25)式给出了由Adams 显式提供估计值,隐式加以校正的四阶Adams预估 校正法公式。 9 f k ?3 )
c yk ?1
p f k ?1=

p yk ?1

=

yk yk

h + (55 f k - 59 f k ?1 + 37 f k ? 2 24 h p + (9 f k ?1 + 19 f k- 5 f k ?1+ f k ? 2 ) 24
p yk ?1 )

=

f (t k ?1 ,

(3-25)

显然上式还需要其他方法如RK4计算出 f k 、f k ?1 、f k ? 2 等初值。
188

3. 2数值积分方法的计算稳定性 3.2.1 数值积分方法计算稳定性概念 一个稳定的系统,当用某种数值积分方法 进行仿真计算时,若该方法或应用该方法时所 选参数不适当则仿真结论出现错误,这就是由 于数值积分方法的稳定性引起的。 对于不同的数值积分方法其稳定性如何、 是否有不同的约束条件,这些都是仿真过程中 的重要问题。考察已知系统,其微分方程为:
189

dy = - 30 y dt
1 y (0) = 3

0≤t≤1.5

1 ?30t 其精确解为: y (t) = e 3
现在取步长 h = 0.1用Euler法和四阶RK4法计算t = 1.5时的 y(t) 值。 Euler 法:y (1.5) = -1.09225 ? 104 RK4 法: y (1.5) = 3.95730 ? 10 精确解: y (1.5) = 9.54173 ? 10-21
190

显然此时数值计算结果是错误的,因为数值积分 方法是一种近似积分方法,它在反复递推计算中将引入 误差,若误差积累越来越大将使计算出现不稳定。 对于系统稳定性的讨论使用微分方程或传递函数; 而后者的稳定性问题需要使用相应的差分方程来进行 讨论。考虑将下述微分方程

dy ? ?y (t ) , dt
称之为测试方程。

y (0) ? y 0

(3-26)

191

特征方程的根在S平面的左半部,原方程稳定,现在以 这种条件来讨论数值计算方法的自身稳定性问题。

3. 2. 2 欧拉法的计算稳定性分析
对测试方程根据欧拉法得到

yk ?1

=

yk

+

h?y k

(3-27)

对上式进行Z变换得到: zy (z) = (1+h) y (z) 显然其特征方程为: z - (1+h ? )=0
192

根据控制理论可知,(3-27) 式的稳定域在[Z]平面上是 单位圆。方程(3-27)稳定的条件是它的特征方程的 根处于单位圆内。即

z ? 1 ? h? 1 ? h? ? jh?
经整理得到 ? 1

< 1

设原微分方程的特征根 ? ? ? ? j? ,并代入上式得到

1 2 2 < (? ? ) ? ? h 由上式易见,在 ?- ? 坐标系下,即[ ? ]平面上的一

193

1 2 h

圆与[Z]平面上的稳定域是一一对应的,或者说[Z]平面 上的稳定域映射到原微分方程的参数平面[ ? ]上也是 1 1 一个圆,圆心在(, 0)处,半径为 , 如下图示。
h h

j [z]
1 h
-1/h

?

j

[?]

0

1

-

0

?

1 若原微分方程是稳定的,而其特征根在 [ ]平面 ? 中映射到单位圆内,则所对应的差分方程也是稳定的; 若原微分方程是稳定的,但其特征根处于
194

[ ? ]平面上映射的圆外,则差分方程是不稳定的。显 然易见[ ? ]平面中稳定域的大小与步长 h有关,h越大 稳定域越小。若考虑 ? 为负实根,即 ? = ? ? 0, ? =0 则h ? - 2 / ?。 因此当采用欧拉法仿真时h应按上式选 择。 3. 2. 3 梯形法稳定性分析 将测试方程(3-26)代入梯形法递推公式可以得到:

y k ?1

h ? y k ? ? ( y k ? y k ?1 ) 2
h ) = (1 + ? 2

对上式进行Z变换可得:
h z y (z) (1- ? 2

) y(z)
195

其特征根为: z =

h? 1? 2 h? 1? 2

将 ? ? ? ? j? 代入上式整理得到: h? 2 h? 2 (1 ? ) ?( ) 2 2 z ? h? 2 h? 2 (1 ? ) ?( ) 2 2 考察上式可知当原系统稳定。其特征根存在负实 部,即? ? 0,此时必存在 z ?1。因此只要原系统稳定 则应用梯形法计算也必稳定,故梯形法属于恒稳定数 值积分方法。 196

3. 2. 4 四阶龙哥-库塔法的计算稳定性分析 与前两种数值积分算法的稳定性分析相类似, 将 测试方程(3-26)代入RK4递推表达式后,经过整理变形 得到: 2 3 4

y k ?1 ? (1 ? ? ?

?

2

?

? ?

6 24

) yk

其中? =h? 。 再对其进行Z变换得到特征根为:

z =1 ? ? ?

?2
2

?

?3 ?4
6 24
197

根据稳定的条件 边界

z ? 1得到[ ? ]平面上稳定域的

参数方程为:

1? ? ?

?2
2

?

?3 ?4
6 24

? e j?

显然RK4法的稳定域与计算步长 h有关,如果 h过 大 ? 值有可能落在稳定域外,因此RK4法是条件稳 定的。
198

3. 3 连续控制系统仿真
3. 3. 1 数值积分方法的MATLAB实现
在控制系统仿真中 , 连续系统的数学模型一般 都可用微分方程来描述 , 对于高阶系统往往可由状态 方程来加以表示。 为了研究系统的动态过程就需要大量的递推求解 微分方程。显然必须借助计算机进行编程加以实现 , 这就是所谓系统仿真。 近年来在控制领域的分析研究中,最有影响与最 为有效的编程语言就是美国Math Works公司开发的 MATLAB,尽管最初MATLAB语言并非是为控制系统
199

的设计者们开发的.然而它的强大绘图功能、数值与矩 阵运算能力以及灵活的编程方式征服了众多的控制系 统设计与研究人员,因此借助MATLAB语言进行数值 积分方法编程,实现控制系统仿真是非常重要的。。 1. Euler法的MATLAB实现 例 3-1 已知二阶系统

? ? x1 ? ?? 0.5572 ? 0.7814? ? x1 ? ?1? ? x ? ? ? 0.7814 ? ? x ? ? ?0 ? u 0 ?? 2 ? ? ? ? ?2 ? ? ? x1 ? y ? ?1.9691 6.4493?? ? ? ?0?u ? x2 ?
200

3. 3. 4 面向结构图的系统仿真
(一) 基于典型环节的系统仿真 由于习惯上往往采用结构图形式对系统进行描述, 因此在对系统进行仿真时有必要研究如何以结构图为 基础加以进行。在此将介绍以典型环节为基本环节的 面向结构图的系统仿真方法,这种方法的主要优点是: (1)便于研究各环节参数对系统的影响。 (2)可以得到每个环节的动态响应。 (3)对于存在系统非线性环节较容易处理。 (4)不但适合单入单出系统;而且对于多入多出系统也 适用。
201

该方法只要建立起各个典型环节的以及各环节之间的 连接关系,就能方便地求取任意环节和系统的动态响 应,由于直接输入各环节的参数,因此研究参数变化 对系统的影响也十分方便。 1. 典型环节的确定
控制系统常见的环节有:比例环节K、积分环节
K s

、惯性环节

K 二阶振荡环节 2 等。为了仿真的方便,必 Ts ? 2?Ts ? 1

K Ts ? 1

T1 s ? 1 、一阶超前(滞后)环节K T2 s ? 1

须选择一种能够代表这些环节的典型环节。经过分析
可以看出只要采用下述形式的环节作为典型环节就可 以满足对系统研究的要求:
202

C i ? Di S Ai ? Bi S
显然,当典型环节中Di=0时即为惯性环节;Ai=Di=0时 典型环节即是积分环节;至于振荡环节则只要用两个 这样的典型环节串联后并加上一个负反馈即可得到。 其中典型环节中的 Ai、Bi、Ci、Di均为对角阵,其对 角线元素值由各个环节中的相关参数确定。根据典型 环节可以得到第i个环节的输出为:

C i ? Di S yi ? ui Ai ? Bi S
203

因此得到:

( Ai ? Bi S ) y i ? (Ci ? Di S )u i
其中, i=1,2,3… n 将其写为矩阵方程形式,得到: (A+BS)Y=(C+DS)U 31) 其中: 0? ? A1 ? B1 ? ? ? A2 ? ,B= ? A= ? ? ? ? ... ? ? ? An ? ?0 ?0

(3-

B2 ...

0? ? ? ? ? Bn ?
204

? C1 ? C= ? ? ? ?0

C2 ...

0? ? ? ? ? Cn ?

? D1 ? , D? = ? ? ? 0

D2 ...

0 ? ? ? ? ? Dn ?

2. 系统联接矩阵 前面给出了从微观角度对系统环节的描述,然而为 了研究一个系统,还必须从宏观角度对系统加以描述,这 就需要将各个环节之间的连接关系用一种数学的形式 加以表示。考察图3-10 所示系统,该系统由 4个典型环 节构成,根据系统结构图可以得到各环节的输入输出关 系方程组。
205

将其写作矩阵形式:

U ? WY ? W0Y0
32)

(3-

206

其中

?0 W = ?1 ?0 ?0 ?

0 ? 4 0? 0 0 ? 1? 1 0 0? 0 1 0? ?



?1 ? W ?00?= ?0 ? ?0 ? ? ?

U =[u1 u2 u3 u4]T Y =[y1 y2 y3 y4]T Y0 =y0
一般称W与W0 为系统联接矩阵. W为一方阵,元素wij 表示结构图中第j个环节对第i个环节的作用;若第j个环 节对第i个环节没有作用,则wij=0. W0 表示外部输入信 号对
207

系统的连接情况。 3.闭环系统的状态方程 将(3-32)式代入(3-31)式,消去U得到:
对上式展开、整理并取拉氏反变换得到:

( A ? BS )Y ? (C ? DS )(WY ? W0Y0 )

考虑在外输入为阶跃信号时,避免对应的初始导数无 穷大,则要求DW0 =0。如果令 Q = B – DW (333) R = CW-A (334)
208

? ? ( B ? DW )Y ? (CW ? A)Y ? CW0Y0 ? DW0Y0

可以得到 :

? QY ? RY ? CW0Y0
如果Q的逆存在,则得到: 若令: 35)

? ? Q ?1 RY ? Q ?1CW Y Y 0 0

~ A

= Q-1R ,

~ B

= Q-1C

(3-

则上式可写作:

~ ~ ? ? AY ? BW Y Y 0 0

(3-36)

显然上式即为一典型的状态方程,从而可利用前面所
209

介绍的数值积分方法对其进行求解,实现对系统的仿 真研究。下面以一个结构图形式的系统为例介绍基于 MATLAB的面向结构图的系统仿真实现。 例 3-5已知系统结构图如图 3-10 所示,采用MATLAB 实现面向结构图的系统仿真。 解:根据系统结构图首先利用MATLAB建立各个环节 参数构成的系数阵(此处的系数阵均为对角阵)与系 统连接阵,在此基础上再依(3-33)、(3-34)、 (3-35)等式建立闭环系统的相应系数阵。下面给出 在MATLAB环境下的相应程序 :

210

AA=diag([1 6.5 1 1]); BB=diag([4 1 3 2]); CC=diag([4 1 2.4 5]); DD=diag([0 0 0 1]); W=[0 0 -4 0;1 0 0 -1;0 1 0 0;0 0 1 0]; W0=[1 0 0 0]'; Q1=BB-DD*W; Q=inv(Q1); R=CC*W-AA; A=Q*R; B=Q*CC*W0;
211

一旦得到了闭环系统的系数阵,即可应用数值积分方 法实现对系统的分析仿真,应该注意的是,在该状态 方程中各个状态变量的实际含义,这里应用RK4法对 其进行仿真分析,下面给出了相应的仿真程序,图 3-11 给出了采用面向结构图法系统阶跃响应曲线。

X=[0;0;0;0];Y=0;t=0;U=25; Tf=15;h=0.01; n=Tf/h; Cb=[0 0 1 0]; for i=1:n K1=A*X+B*U; K2=A*(X+h*K1/2)+B*U;

212

K3=A*(X+h*K2/2)+B*U; K4=A*(X+h*K3)+B*U; X=X+h*(K1+2*K2+2*K3+K4)/6; Y=[Y,Cb*X]; t=[t,t(i)+h]; end plot(t,Y) grid

213

214

在这里由于采用的是面向结构图形式的系统仿真,因此 研究因各环节参数发生变化对系统的影响就十分方便, 仅对环节系数阵加以修改即可,如果将第三个环节修 改为:2.4/(s+1)时我们只需要修改B阵的相关元素即 可,其它部分无需改动,图3-12 给出了改变环节三参 数后的仿真结果,由响应曲线看出,系统的响应快速 性明显提高了,同时超调量也有所降低。另外采用面 向结构图仿真还可十分方便地研究各个环节输出的响 应情况,只要根据需要修改相应系数阵即可。图313(a)、图3-13(b)和图3-13(c)分别给出了环节1、环节 2和环节4的输出响应。

215

216

217

218

(二) 基于Connect联接函数的系统仿 真 除了前面介绍的面向结构图的系统仿真方法外,还可利 用MATLAB与控制工具箱所提供的相关函数与命令实 现以结构图形式给出的系统的仿真。下面通过一个实 例来说明如何以这种方法进行仿真工作。 例 3-6 已知系统如图 3-14所示,研究该系统的阶跃响 应。 显然该系统是一类以结构图形式给出的较复杂系统, 利用MATLAB和控制工具箱可以方便的实现系统仿真。 就基于Connect联接函数的仿真而言一般可分为四步 进行。 (1)建立系统中间状态模型
219

220

首先应对系统各环节进行编号(如图3-14中所标示编 号),然后即可据此与给定系统参数建立有关中间模型, 以便运行时在MATLAB的工作空间建立相应的系统环 节模型. 其基本语句和格式为:

nblocks=10; n1=1;d1=[0.01 1]; n2=[0.17 1];d2=[0.085 0]; n3=1;d3=[0.01 1]; n4=[0.15 1];d4=[0.051 0]; n5=70;d5=[0.0067 1]; n6=0.21;d6=[0.15 1];
221

n7=130;d7=[1 0]; n8=-0.212;d8=1; n9=-0.1;d9=[0.01 1]; n10=-0.0044;d10=[0.01 1]; blkbuild;
其中ni、di表示第i个环节传递函数的分子、分母多项 式的系数向量,如果所研究的系统是以状态方程形式给 出的各个环节,则应采用ai、bi、ci、di来表示环节状态 方程的系数阵。对于每个环节只能取其二者之一, MATLAB会自动根据所输入的阵列元素数量加以判定 并进行相应的运算。Nblocks需根据系统环节个数对 其
222

进行赋值,blkbuild 语句将对nblocks语句后的各环节 输入情况进行校验,并自动生成系统的隐含中间状态 方程模型(a,b,c,d系数阵)。 (2)建立系统连接矩阵 在MATLAB工作空间建立了系统各环节的模型之 后,还应该建立反映各环节之间相互连接关系的Q阵, 以及输入、输出向量(inputs、outputs),下面结合实例 将Q阵与inputs、 outputs 列写规则加以介绍。根据已 知系统可以得到系统连接阵Q和输入、输出向量:

223

Q=[ 1 0 0; 2 1 10; 3 2 0; 4 3 9; 5 4 0; 6 5 8; 7 6 0; 8 7 0; 9 6 0; 10 7 0]; inputs=1; outputs=7;

224

Q阵的行数与系统环节数相同,Q阵的第一列应为相应 环节编号,第i行的第二列及以后的所有元素应为进入 第 i个环节的所有环节编号。输入变量inputs为输入信 号进入的环节编号;输出变量outputs为系统信号输出 所对应的环节。如果是多入多出系统,输入变量和输 出变量应为一向量。 (3)建立整体系统模型 完成了上述工作之后就可以应用connect( ) 函数建立 整体系统的数学模型,其基本格式如下:

[A , B , C , D]=connect(a,b,c,d,Q,inputs,outputs) 或 G =connect(a,b,c,d,Q,inputs,outputs)
在connect( )函数中的a、b、c、d系数阵是由blkbuild 函
225

数自动生成的,其返回变量既可是系数阵,也可是系 统阵G,不论那种情况都将生成完整的系统模型,从 而就可利用它们对该复杂系统进行分析研究。 (4)系统仿真 由于已经求出了系统的系数阵,为了获得系统阶 跃响应,可以利用step(A,B,C,D)求取系统的阶跃相应, 图 3-15即为系统的阶跃响应曲线。当然也可应用数值 积分方法对其进行研究,可参考前面内容进行。

226

227

3.3.5 离散相似法的系统仿真
(一)离散相似的概念 假设有一连续系统如图 3-16 所示。

其输入、输出分别为u(t)与y(t),用一采样周期为T的采 样开关将输入、输出分别离散化, 要求其输出y*(t)在 采样时刻的值等于原输出y(t)在同一时刻的值, 即
228

y*(t)= y(kt), t=kT 成立,其中 k=0,1,2,… 显然,如果只是简单地在原系统的输入、输出端上人为 地加上一个采样开关,其输出y*(t)无法与原系统的输出 y(t)相同。因为输入信号u(t)经过采样开关后离散成 u*(t),此时若将u*(t)直接输入原系统其系统输出当然不 会按照原来的规律变化,因此这是我们所不希望出现的。 为了使输入信号保持连续,就需要加设相应的保持器, 如图 3-17 所示,由于保持器不能完全复现原输入信号 u(t),因此一般情况加到原系统输入端的信号u~(t)无法 等同于u(t),所以经过这样离散化得到的模型必然具有 一定的近似性.由前面的分析可以清楚地看出,在应用 离散相似法对连续系统进行仿真过程中,加设的采样开
229

关, 保持器都是虚拟的, 并不存在真实的物理装置。

(二) Z域离散相似法仿真 对于连续系统采用Z域离散相似法进行仿真其主要步 骤为: (1)画出系统结构图。 (2)在适当位置加设虚拟采样开关、保持器。
230

(3)将原系统连同引入的保持器,通过Z变换求解得到 系统的脉冲传递函数。 (4)通过求取Z反变换求得系统差分方程。 (5)依差分方程编制仿真程序。 将系统输出Z变换与输入Z变换之比称为系统的脉冲传 递函数,即 Y ( z) ? Z [G h ( s ) ? G ( s )] G(z) = U ( z) 其中 G(s)h ? 为保持器的传函 G(s) ? 为系统传函 由此可见,选择不同的保持器,将得到不同的G(z), 有关各种保持器的数学描述与特性分析请参阅有关书
231

籍,在此不再赘述。下面通过例子来说明采用Z域离 散相似法进行仿真的步骤。 例 3-7 已知连续系统的传递函数为:

Y ( s) k ? G(s) = U ( s) s ? a
采用Z域离散相似法对其进行仿真(应用零阶保持 器)。 解:首先在系统输入、输出端设置采样周期为T的虚 拟采样开关,并且在输入采样开关后面加入零阶保持 器, 如图3-18 所示。

232

据此就可以写出系统的脉冲传递函数

? 1 ? e ?Ts Y ( z) k ? ? ? z? ? G(z) = ? U ( z) s s?a? ? ? ? ? k = (1 ? z ?1 ) z ? ? ? s( s ? a) ?
k (1 ? e ? aT ) = a ( z ? e ? aT )
233

y(k+1) = e

? aT

k ? aT y(k ) ? (1 ? e )u(k ) a

利用上述方法可以容易地求出各种线性定常系统的离 散相似仿真模型,但由前面已经知道选择不同的保持器 会得到不同的模型。由所求出的上式不难看出,当在离 散相似法中选择零阶保持器时,所得到仿真递推公式与 采用Euler法得到的结果是相同的。如果选择三角形保 持器,所得到递推公式将与梯形法得到的结果一致。对 于以上述形式给出的仿真模型应用MATLAB语言很容 易编写出相应的仿真程序,实现系统仿真。 (三) 时域离散相似法仿真 对于由状态方程描述的系统进行离散化仿真就属
234

于时域离散相似方法仿真,因此也需要在连续系统中加 入虚拟采样开关和保持器,只是这里的数学模型为状态 方程形式,如图3-19 所示。

235

设系统的状态空间表达式为:

? x

= Ax + Bu

(3-37) (3-38) (3(3-

y = Cx + Du
对上式进行拉式变换并经整理后得到:
X(s)=(sI-A)-1x(0)+(sI-A)-1BU(s)

39)
Y(s)=CX(s)+DU(s) 40)

对式(3-39)进行拉示反变换并利用卷积分得到 t

x (t ) ? F (t ) x (0) ? ?0 F (t ? ? ) Bu(? )d?
(3-41)
236

其中F(t) = e At, 称为系统的状态转移矩阵。式(3-41) 就是连续系统状态方程的解。若设采样周期为T,将 t=(k+1)T 代入式(3-41)经整理得到

x[(k ? 1)T ] ? e AT x(kT) ? e AT ?kT

( k ?1)T

e

A( kT ? ? )

Bu(? )d?
(3-42)

如果采用零阶保持器,即 u(t)=u(kt) kT?t?(k+1)T 则式(3-42)经变换得到 x(k+1)=F(T)x(k)+G(T)u(k) (3-43) 其中 F(T) = e AT
237

G(T) =

?0

T

F (T ? ? ) Bd?

从而的得到采用零阶保持器的系统离散化模型 x(k+1)=F(T)x(k)+G(T)u(k) (3-44) y(k)=Cx(k)+Du(k) 45) 若采用一阶保持器则可以得到 x(k+1)=F(T)x(k)+G(T)u(k)+G1(T)(k) T (3-46) ?0 e A(T ?? ) B?d? 其中 G1(T) =

(3-

式(3-46)即为采用一阶保持器时的离散化状态方 238 程。

由前面推导可以清楚看出,在对系统数学模型进行时域 离散化过程中的计算求解状态转移阵是最关键的环节, 特别是一种适合计算机运算的数值解方法更是大家所 关心的问题。有关求解状态转移阵的数值解的方法很 多,下面不加推导地给出基于泰勒级数展开方法的状 态转移阵迭代计算式 F(T) = e
AT?

AT AT AT AT I ? AT ( I ? (I ? (I ? ? ? (I ? ) ?)) 2 3 N ?1 N

(3-

47) 在此基础上可以导出 :
G(T)=A-1(F(T)-I)B (3239

48)

利用(3-47)与(3-48)即可由计算机实现对状态转 移阵 e AT与 G (T)阵的数值递推计算。 例3-8 试求积分环节G(s)=K/s的状态方程离散化模型。 解:该环节的状态方程和输出方程分别为:

? x y

= Ku

= x

(3-49) (3-50)

因为 A=0, B=K, C=1, D=0。故离散化模型的系数阵 为 F(T) = e AT = e 0T = 1 G(T) =

?0 e 0(T ?? ) Kd?
T

= KT

G1(T) = T e 0(T ?? ) K?d? ?0

1 ? KT 2 2

240

所以采用零阶保持器时系统的离散化模型为 x(k+1) = x(k) + KTu(k) 51) y(k) = x(k) 52) 保持器为一阶时,离散化模型为:1 KT 2 u (k ) ? 2 x(k+1) = x(k) + KTu(k) + y(k) = x(k) u (k ) ? u (k ? 1) ? 其中 u?

(3(3-

T

例 3-9 试求比例积分环节
241

KI G(s) = Kp + s
的离散化模型。 解:系统的状态空间表达式为:

53)
54)

? x
y

= KI u = x + Kpu

(3(3-

由于 A=0, B=KI, C=1, D=Kp。则得到

?0 e 0(T ?? ) K I d? G(T) =

F(T) = e AT = e 0T = 1 T

= K IT
242

1 G1(T) = ? e 0(T ?? ) K I?d? ? K I T 2 0 2 所以采用零阶保持器和一阶保持器时,系统的离散化 模型分别为: x(k+1) = x(k) + KITu(k) (355) y(k) = x(k) + Kpu(k) (3-56)
T

保持器为一阶时,离散化模型为:1 ? K I T 2 u (k ) x(k+1) = x(k) + K Tu(k) +2
I

y(k) = x(k) + Kpu(k) + KpT

? u (k )

243

例3-10 程离散化模型。 解:系统的状态方程与输出方程为 :

K 试求惯性环节 G(s) = s ? a

状 态 方

? x
(3-57)

= -ax + Ku (3-58)

y = x 故 A = -a, B = K, C = 1, D = 0。 K F(T) = e – aT

G(T) = A -1 (F (T) - 1) B = a K ?aT K T ?0 e A(T ?? ) B?d? a 2 (e ? 1) ? a T G (T) = =
1
244

(1 ? e ?aT )

从而得到采用零阶保持器和一阶保持器时,系统的离 散化模型分别为: x(k+1) = e -aT x(k) +K (1 ? e ? aT ) u(k) a y(k) = x(k) 与 x(k+1) = e -aT x(k) + u(k) + K [ ] K (e ?aT ? 1) ? T u (k ) ? 2 y(k) = x(k) a a
对于比例惯性环节 G(s) = c ? ds
a ? bs 其系统的状态方程与输出方程可写做为:

59)

? x

=

a d x? u b b

(3245

60)

c a d y =( ? )x ? u d b b

(3-

从而可以得到采用零阶保持器时,系统的离散化模型 d - aT ? aT 为: (1 ? e b )u (k ) eb a x(k+1) =c ax(k) + d
(

y(k) =

d

?

b

) x(k ) ?

b

u (k )

通过上述几个例子给出了几种典型环节的离散状态方 程,显然对于复杂系统都可以转化成一些典型环节组 成,因而利用这些典型时域离散化模型就可以容易地 编制仿真程序实现对系统的仿真研究。
246

(四)离散相似模型的校正
由于在离散相似法对系统进行仿真时必须引入各种保 持器,而保持器的引入一般会引起系统幅值与相位上的 失真,为了弥补这方面的问题,可以在模型中加入校正环 节,适当调整相关参数使离散化模型尽可能地接近原系 统。对于这种失真进行校正可以采用多种方法,其中比 较有效的方法就是根匹配法。在根匹配法中主要是通 过引入一校正环节加以实现。一般校正环节取作如下 形式: ?Ts G (s)=


?e

这里 ? 主要用以实现对幅值的校正,主要用来实现对
247

相位的校正。下面通过一个例子说明根匹配校正的具 体方法。 例 3-11 已知系统其传递函数为: 1 G( s) ? s?a
其系统结构图如图3-20 所示,对其进行离散相似模型 的校正。

248

解:根据离散相似法的思想,并考虑到对失真的补偿, 其结构图如图 3-21所示。在反馈回路引入了延迟环节, 主要是考虑到在计算反馈回路之前,必须先计算前向 回路。

249

因为校正环节为指数函数,分析计算十分复杂,所以 将校正环节做一次近似为:

?e ?1 ? ?Ts
从而可以得到系统闭环脉冲传递函数:

?Ts

Y ( z) G( z ) ? GB (z) = R( z ) 1 ? G( z ) Kz ?1

?Tz(?z ? (1 ? ? )) = 2 z ? (?KT? ? 1) z ? ?KT (1 ? ? )
(3-

61)
250

显然,可以看出上式为一个二阶系统,而原系统是一 阶系统,因此首先应进行阶次匹配。在这里取 =1,则 得到一阶系统为:
GB(z) Y ( z )
?Tz ? R( z ) z ? (?KT ? 1)

(3-

62)

其极点为: z极 = 1 ? KT 在实现阶次匹配后,还应保证该极点与原系统极点相 匹配,因此应存在: e –KT = 1 - KT ?

?

1 ? e ? KT = KT
251

从而得到校正后的闭环离散相似脉冲传函:
Y ( z) z (1 ? e ? KT ) GB ( z) ? ? R ( z ) K ( z ? e ? KT )

(3-63)

这样当校正环节取 ? = 1, = 1 ? e ? KT 时,可以使 ? KT 离 散相似化系统与原系统的特性保持一致。从而得到用 于系统仿真的差分方程:

y(k ? 1) ? e

? KT

(3-64) 即可对其进行编程,实现对系统的离散相似法仿真。 值得注意的是,当系统阶次较高时,对 ? 、 的确定 ? 是较困难的,所以往往采用试探法进行。
252

1 y(k ) ? (1 ? e ? KT )r (k ? 1) K

3.3.6 含有非线性环节的系统仿 真 (一)几种典型非线性环节 一个控制系统严格地讲往往具有非线性特性,人们为了 简化问题的分析,在系统非线性不严重时,将其理想化成 线性系统,从而就可以应用线性系统理论对其进行分析 研究。然而实际系统中存在许多非线性特性明显的系 统,对这些系统的分析研究就必须按照非线性系统来 加以进行,这些非线性系统的非线相性特性一般可由 若干典型非线性特性来加以抽象和描述,首先简要讨 论和介绍几种典型非线性特性以及相关的仿真问题。 1. 饱和非线性特性
253

图3-22 给出了饱和非线性的图示描述,图3-23 给出 了可以在计算机上进行分析仿真的特性流程图,根据 图示流程在计算机上加以编程可以方便地实现饱和非 线性特性的仿真。

254

255

2.死区非线性特性 图3-24 与图3-25 分别给出了死区非线性及其相应的 仿真流程图。

256

257

3.死区继电非线性特性 图3-26 与图3-27 分别给出了死区继电非线性特性及其 对应的流程图。 常见的典型非线性环节除了上述几种之外,还有死区 分段特性、摩擦力特性、滞环以及具有滞环的继电非 线性等。在控制系统分析研究过程中,如果需要这些 环节可参考有关书籍,在此不再赘述。

258

259

260

(三)含有非线性环节的系统仿 真 因为在实际控制系统中,绝大部分非线性系统都 可以使用含有非线性环节的系统来加以描述,因此以 这种方式研究非线性系统具有重要的实际意义。在前 面曾介绍过时域离散相似法,应用这种方法可以方便 地对非线性系统进行分析研究,因为采用这种离散化 方法对系统进行分析研究时,每计算一步,各个环节 的输入与输出都需要重新计算一次,因而非线性环节 特性就能够很容易地加入到整个分析计算程序中。下 面通过一个例子说明对于含有非线性环节的系统如何 进行分析仿真。 例3-12 已知非线性系统如图3-28所示。试对其进行仿 真。
261

其中

1 G1(s) = , s

1 G2(s) = s?2
= Ax+Bu y = Cx
262

解: 首先写出线性部分的状态方程

? x

其中
A=

因此可以得到

?0 0 ? , B = 1 ? ? ?1 ? 2? ?0? ? ? ? ?

, C = [0 1]

1 0 ? ? F(T) = e AT = ? , G(T)? = T ?1 ? ?T ? ? 2 T ) e ? 2T ? 2T ) (1 ? e ?2 ? ?( 2 ? 1 ? e ? ? ? ? ? 如果采用零阶保持器,则可以得到下述形式的线性部分 的离散状态方程 x(k+1)= F(T)x(k)+ G(T)u(k) 即 x1(k+1)=x1(k)+Tu(k) (365)
263

1 T ? 2T ) x1 (k ) ? e ? 2T x 2 (k ) ? ( ? 1 ? e ? 2T )u (k ) x2(k+1) = (1 ? e 2 2 (366) y(k)=x2(k) (3-67)
在本例中非线性环节为饱和非线性,根据前面给出的饱 和非线性特性流程图可以容易地得到该环节输入输出 之间的对应关系。由系统结构图可以清楚地看到,该 环节的输入就是系统给定与反馈信号的综合,其输出 信号就是线性环节的输入信号。因而采用前面所介绍 的时域离散相似法进行含有非线性特性环节的系统仿 真是十分方便的。当然对于含有非线性特性环节的系 统应用其他方法(包括面向结构图仿真的其他方法) 进行仿真研究也是完全可行的。
264

考虑到由 (3-68)和 (3-69) 式给出的饱和非线性特性以及 信号综合的关系: u(k) = fN ( e(k) ) (3-68) e(k) = r(k) - y(k) (369) 以及(3-65)、(3-66)和(3-67)诸式即可实现该系 统的仿真。图3-29 给出了系统仿真程序框图。虽然该 系统是一较为简单的非线性系统, 然而这种对非线性 系统进行仿真研究的思路与结论,完全可以应用于更 加复杂的非线性系统。另外在本例中对于线性部分其 状态方程是按整体考虑的,如果对于线性部分按照典 型环节进行研究,则所得到仿真模型(差分方程)将会有 所不同,其原因是十分清楚的。
265

266

3.4 采样控制系统仿真
随着计算机科学技术的发展,在工业领域计算机 控制系统获得了越来越广泛的应用,因此对于这类系统 的分析与设计越显重要,这类系统的被控对象的状态 变量是连续变化,然而它的输入变量和控制变量却只 是在采样点取值的一系列脉冲序列。其数学模型为差 分方程和离散状态方程,这一类在一处和多处存在采 样脉冲序列信号控制的系统称为采样控制系统。由于 采样系统特别是数字计算机控制系统具有适应性强, 能实现各种复杂的控制等优点,其研究和应用获得了 较快的发展,相应地对采样控制系统仿真试验的要求 也越来越强烈。本节将简要地要介绍采样控制系统的 基本仿真原理与方法。
267

3. 4. 1 概述

典型的采样控系统如图3-30所示。系统的模拟信 号经过采样器与模数转换器变成计算机可以接受的数 字信号,依据选定的控制策略经计算机处理之后以数字 量输出,再经过数模变换成模拟量以控制被控对象。
268

采样控制系统一般由两部分组成:一部分是被控对象 一般属于连续部分,另一部分是速度调节器或数字控 制器属于离散部分.相对于连续系统而言,它增加了数 字控制器、采样器和保持器,但这与前面所介绍过的 连续系统离散化时所引入的采样开关或保持器是不同 的,在连续系统进行仿真时,为保证所建立的仿真模 型与原系统动静态特性一致,要对连续模型进行离散 化处理,为此人为地加入了虚拟的采样开关和保持器; 而在采样系统里,采样开关和保持器是真实存在的, 图3-31给出了一个简单的采样控制系统结构图。对这 类采样系统进行仿真分析,一般有两种方法:一种是 将连续部分的微分方程和离散部分的差分方程分别处 理,对连续离散混合模型进行仿真;另一种是将连续 部分模型和离散部分模型统一处理。 269

对于图3-31 所示的简单采样系统,在将系统的连续部 分进行离散化时一般有两种情况: 一、 如果只需要知道输出量,不需要求取计算中间状 态变量,而且Z变换又可以比较容易地求出来,则连 续部分只需要表示成一个线性环节,它的输入是经过 信号重构器后形成的连续信号,不存在信号失真和截 断误差,所以在保证能描述系统响应的前提下可以使 用较大的计算步长,但不得大于实际采样周期。为了 简化仿真计算可以采用同样的采样周期求取数字和连 续部分的差分方程,然后以采样周期作为计算步长进 行数字仿真。 二、 如果不仅需要计算系统输出量,而且还要计算系 统中的其它状态变量;或者不仅要知道各种采样时刻
270

的输出量,而且还需要求出相邻两采样点之间的输出响 应等。在以上这些情况下进行仿真时,可以将连续部 分按若干个典型环节进行计算,此时应在每个连续环 节前分别设置虚拟采样开关与信号重构器进行离散化 处理。为了提高仿真精度,连续部分的计算步长h一般 要远远小于实际采样周期Ts。为了程序实现的方便一 般取Ts=Nh(其中N为正整数)。 3. 4. 2 采样控制系统仿真的一般方法 前一部分概述中已经述及,在采样控制系统中一般 总会包括两部分:一部分是连续的被控对象,另一部分 则是离散的数字控制器或调节器。因此在进行采样控 制系统仿真时就一定会涉及到连续部分仿真时的计算 步长和离散数字部分仿真时的采样周期,自然也就涉 及采样周期与计算步长之间的关系问题, 对于采样周
271

期与计算步长之间关系一般只有两种选择方式:1 仿 真计算步长小于采样周期;2 仿真计算步长等于采样 周期。下面分别介绍这两类仿真的不同应用。 (一)以采样时刻系统输出为目标的的仿真 对于这种仿真要求只需求出系统在采样时刻的输 出变量值,所以在对采样控制系统进行仿真时可以采 用仿真步长与采样周期相等的方法进行。 1. 直接法 当所研究的采样控制系统可以得到相应的系统闭 环脉冲传递函数时,即可对其进行Z反变换得到系统 的差分方程,建立系统仿真模型对其直接加以仿真。 若考虑所研究采样控制系统的闭环脉冲传递函数为:
272

y ( z ) b0 ? b1 z ?1 ? ? ? ? ? bk z ? k ?( z ) ? ? r ( z ) 1 ? a1 z ?1 ? ? ? ? ? a m z ? m

(3-70)

则可得到系统的仿真模型
(3-71) 据此即可对系统进行仿真研究,显然这种仿真研究所得 到结果只是系统在采样时刻的输出. 2. 分别求取法 当所研究的采样控制系统是以结构图形式给出, 则可采用分别求取离散部分和连续部分的脉冲传函, 并对他们加以变换得到仿真模型,然后在系统采样时刻 分别对这些模型加以计算,根据系统信号的综合关系 再对计算值进行综合处理,为下一采样时刻的计算
i ?1 j ?0
273

y(n) ? ? ? ai y(n ? i) ? ? b j r (n ? j )

m

k

做好准备。图3-32 给出了所研究的采样控制系统结构 图。

其中 (3-72)

u( z) 1 ? az ?1 D(z) = ? k1 e( z ) 1 ? z ?1
274

1 ? e ?Ts Gh(s) = s
k0 Gp(s) = T0 s ? 1

(3-73)
(3-74)

因此可以很容易地得到控制器的仿真模型为: u(n) = u(n-1) + k1[e(n)-ae(n-1)] (3-75) 考虑到零阶保持器的存在,广义控制对象的脉冲传函 为: (z) = Z [Gh(s) Gp(s)] = k0 其中 ?

? e ?T / T0

(1 ? ? ) z ?1 1 ? ?z ?1

(3-76)

275

所以可得到相应的仿真模型为: y (n) = ? y(n-1) + k0(1-) u(n-1)

(4-77)

由系统结构图可知在信号综合点处存在: e(n)=r(n)-y(n) (4-78) 根据系统初始条件:e(-n)=u(-n)=0(n?0), 以及 y(0)=0, 考虑到(3-75)式、(3-77)式和(3-78)式之间的 关系,即可实现对采样系统的仿真研究。 (二)以受控对象状态变量为目标的仿真 如果仿真的要求不仅是要知道系统输出的情况,而 且还要求计算出受控对象内部的状态.那么就要在受 控对象的各个环节之间加设虚拟的采样开关及信号重 构器。此时的这部分工作就是采用的前面所介绍的连
276

续系统离散相似法,因此这里的虚拟采样开关的采样周 期就是仿真计算步长,由于所研究的系统是采样控制 系统,所以在系统仿真时就必须考虑控制器采样周期 与虚拟采样开关的采样周期(计算步长)之间的协调 与同步问题。考虑受控对象的传递函数为:

k0 G(s) = s (T0 s ? 1)
因为现在对所研究的采样系统,除了需要系统响应之外, 还要研究系统被控对象的中间状态,因此可构造采样系 统的仿真结构图如图3-33 所示。

277

278

图 3-33 中 Ghp1(s)=Gh1(s)Gp1(s), Ghp2(s)= h2(s)Gp2(s)。 Gh1(s)、Gh2(s)分别为零阶保持器(信号重构器) 以保持信号的连续,为了系统研究的需要将被控对象 分解成:Gp1(s)=k0/(T0s+1)、Gp2(s)=1/s 两部分相串 联的形式。由于系统存在离散与连续两部分因此这两 部分的采样周期一般是不同的,进行仿真时必须按这 种情况加以考虑,通常取T=NT’关系(N为正整数), 这里取正整倍数的关系主要是考虑计算上的简便,同 时也考虑到虚拟采样周期(计算步长)远小于系统实 际采样周期。两个采样开关每间隔时间T进行一次协 调以实现同步。在此可得到系统的脉冲传函分别为:
279

G1(z) = Z [Gh1(s) Gp1(s)]
?1 ? e ? sT k 0 ? ? =Z ? s T0 s ? 1? ?

(1 ? ? ) z ?1 = k0 ( ? ? e ?T / T0 ) 1 ? ?z ?1
G2(z) = Z [Gh2(s) Gp1(s)]
?1 ? e ? sT ' 1 ? =Z? s s? ? ?

(3-79)

T' = z ?1

(3-80)
280

由(3-79)和(3-80)式不难得到相应的仿真模型 分别为: y(n+1)=?x(n)+k0(1-?)u(n) (3-81) y(n+1)=y(n)+T’x(n) (3-82)
若在该系统中选用前面由式(3-72)给出的控制器, 其相应的仿真模型为: u(n)=u(n-1)+k1[e(n)-ae(n-1)] (3-83) 显然可以看出,在对该系统仿真时根据初始条 件对(3-82)与(3-81)、(3-83)式分别进行计 算,但是注意到他们的计算步长是不同的,因此在 进行系统仿真时,应按照两个程序循环进行,(382)式计算步长为T ’,每隔T ’计算一次作为程序内 281 环;(3-81)、

(3-83)式的计算步长为T,每隔T 间隔计算一次,即程 序内环计算N次,程序外环(3-81)、(3-83)式计 算一次。

282

第五章

连续控制系统计算机辅助设计

5.1 状态可控规范型与状态可观规范型 (一)状态可控性判定 状态可控性是控制输入u(t) 支配状态向量x(t) 的能力。 对于线性定常系统状态可控性的代数判据 是:若状态可控性矩阵Qc为

Qc ? [ B ? AB??? An?1B] (5-1) 则线性定常系统 ?(A,B,C,D)状态完全可控的 充分必要条件就是Qc的秩为n。
283

?

284

例 5-1 已知系统

判断系统状态的可控性。 解:根据可控性代数判据应用MATLAB不难得到相 应的判定程序。 A=[-1 -2 -2;0 -1 1;1 0 -1]; B=[2;0;1]; n=rank(A); Q=[B A*B A^2*B] r=rank(Q) if r==n
285

?? 1 ? 2 ? x ? ? 0 ?1 ? ?1 0 ?

? 2? ?2? 1 ? x ? ?0 ? u ? ? ? ?1 ? ? 1? ? ? ?

disp(‘System is Controllable.’) Q else r end 运行结果: r= 3 System is Controllable. Q= 2 0 1 -4 1 1 0 0 -5
286

显然可控阵的秩与系统的阶数相等,因此系统状态 是完全可控的。 例 5-2 考察多输入系统的可控性,已知系统系数 阵为: 2 3? 9? ?1 ?1 A ? ?1 4 6 ? , B ? ?0 0? ? ? ? ? ?2 ?2 1 7? 0? ? ? ? ?

解:应用MATLAB还可编写更为一般化的判定程序,
根据给定的系数阵自动确定状态可控阵。 A=[-1 -2 -2;0 -1 1;1 0 -1]; B=[2;0;1]; n=rank(A);

287

l=length(A); Q=zeros(l); Q(:,1)=B; for i=2:l Q(:,i)=A*Q(:,i-1); end r=rank(Q) if r==n disp('System is Controllable.') Q else r end
288

运行结果: r= 3 System is Controllable. Q= 2 -4 0 0 1 0 1 1 -5 与前例相同,该系统状态是完全可控的。 对于判断系统状态的可控性,MATLAB提供了 可控性判断函数ctrb( ),直接进行判断。考查 例5-2 应用 ctrb( ) 函数可直接进行系统可控性判定。
289

A=[-1 -2 -2;0 -1 1;1 0 -1];
B=[2;0;1]; C=[0 0 0]; D=[0]; sys=ss(A,B,C,D); Con=ctrb(sys) r=rank(Con) 运行结果: Con = 2 -4 0 0 1 0 1 1 -5 r = 3

290

由ctrb( )函数可直接得到系统可控阵,由此即可判定 系统的可控性。 (二)可控规范型及其MATLAB的实现 在系统进行分析与综合过程中, 可控规范型的变换 是最基本与常用的一类变换。 若系统状态是完全可 控的,则选取QC的n个线性无关的列向量构造变换阵 T1 (5-2) 对系统 ?(A,B,C)进行非奇异变换,则可以得到第 一可控规范型 ?1(Ac1,Bc1,Cc1 ),其系数阵之间满足 下述关系:
291

Qc ? [ B ? AB??? An?1B ] =

B=T1Bc1 AT1=T1Ac1 CT1=Cc1
其中

(5-3)

Ac1

? ? an ? ? ?????? ? ? ?1 ? ? a n ?1 ? ?? 1 ? ? a n?2 ? ? ? ? ? ? ? ? 1 ? ? a1 ? ? ? ?

Bc1

?1? ??? ??0? ?0? ?? ? ?0? ? ?

Cc1 ? [ CB ? CAB ? ? ? CAn?1 B ]

(5-4)
292

式 (5-4) 即构成第一可控规范型系统?1(Ac1,Bc1,Cc1 ) 的系数阵。 例 5-3 已知系统?(A,B,C)的系数阵为:
?2 A ? ?0 ? ?0 ? 0 4 0 0? ?1 ? 1 ? , B ? ?0 ? , C ? ? 1 ? ? ? ?1 ? 4? ? ? ? 1 0?

试判断系统的可控性,如果状态完全可控,则求取 系统的第一可控规范型及变换阵。 解:利用MATLAB编制程序,首先判断状态是否完 全可控,再根据(5-2)和(5-3)式确定变换阵与 第一可控规范型各系数阵。
293

A=[2 0 0;0 4 1;0 0 4]; B=[1;0;1]; C=[1 1 0]; n=rank(A); Q=[B A*B A^2*B]; r=rank(Q); if r==n T1=Q Ac1=inv(T1)*A*T1 Bc1=inv(T1)*B Cc1=C*T1 else r end
294

运行结果得到: T1 = 1 2 0 1 1 4 Ac1 = 0 0 1 0 0 1 Bc1 = 1 0 0
Cc1 = 1 3

4 8 16 32 -32 10

12
295

如果可控规范型变换阵取作下述形式 : ? 1 ? ? a1 1 ? n ?1 B ??? AB ? B ] [A ? ? ? ? ? T2= (5-5) ? ? ? a n ?1 ? a1 1? 则可导出系统的第二可控规范型系统?(A,B,C)进行 非奇异变换后,则可以得到第二可控规范型 ?2(Ac2,Bc2,Cc2 ), 其系数阵之间满足下述关系:
AT2=T2AC2 B=T2BC2 CT2=CC2 (5-6)

296

其中
? ? ? ?? ? ?? ? ? an ? ? ? ? ? ? ? ? an ?1 ? ? 1 ? ? ? ? 1? ? ?? ? a1 ? ?

Ac 2

Bc 2

?0? ??? ? ? ? ?0? ? ? ??? ?1? ? ?

例 5-4 已知系统?(A,B,C)的系数阵同例5-3,试 求取系统的第二可控规范型及变换阵。 解:由例5-3已知该系统是状态完全可控的,因 此在本例中无需再进行系统可控性判断,即可直接 求取系统的变换阵以及相应的第二可控规范型。
297

A=[2 0 0;0 4 1;0 0 4]; B=[1;0;1]; C=[1 1 0]; AS=poly(A); TA=[A^2*B A*B B]; TB=[1 0 0;AS(2)1 0;AS(3) AS(2) 1]; T2=TA*TB T=inv(T2) Ac2=T*A*T2 Bc2=T*B
Cc2=C*T2 运行结果:
298

T= 0.2500 0.5000 1.0000 Ac2 = 0 1 0 0 32 -32 Bc2 = 0 0 1 Cc2 = 14 -7

0.5000 -0.2500 2.0000 -0.5000 8.0000 0 0 1 10

1

299

5.1.2 系统状态可观性的判定与可观规范型
线性定常系统?(A,B,C,D)状态完全可观的充分 必要条件就是系统的可观性矩阵

Qo ? [ C ? CA? ? ? CA

n ?1 T

]

(5-7 )

的秩为n。 如果系统的状态完全可观,则可以构造一变 换阵P1,形成系统可观规范型。可观规范型变换阵 P1的逆阵存在下述关系: P1-1=Qo
应用可观规范型变换阵P1对系统?(A,B,C,D)进 行可观规范型变换,可以得到系统的第一可观规范型, 其系数阵形式为:
300

? 1 ? ? ? ? ? ? ? ? Ao1 ? ? ? 1? ? ? ?? ? ? ? ?? ?? an ? ? an ?1 ? ? a1 ? ? ?

? CB ? ? ? ? ? ? ? CAB ? ? ? ? Bo1 ? ? ? ? ? ? ? ? ? ? ? ? n ?1 ? ?CA B ? ? ?

Co1=

[ 1? 0 ? 0 ]

(5-8)

并满足下述变换关系: AP=PAo1 B=PBo1 CP=Co1

(5-9)
301

考察第一可观规范型,并与第一可控规范型相对照可 以看出,他们之间存在明显的对偶关系,即
T Ac1 ? Ao1 , T B c1 ? C o1 , T C c1 ? Bo1 ,

例 5-5 已知系统的系数阵为:
4 4? ? 4 A ? ?? 11 ? 12 ? 12 ? ? ? ? 13 14 13 ? ? ?

?1? B ? ?? 1? ? ? ?0? ? ?

C=[ 1 1 1 ]

试判断它的可观性,并将其变换成第一可观规范型。
302

解:根据判断系统?(A,B,C,D)状态完全可观的充 分必要条件,可以构造可观性矩阵,当判其满足状 态完全可观的充分必要条件,则可利用(5-9)式求出系 统的可观规范型。 A=[4 4 4;-11 -12 -12;13 14 13]; B=[1 -1 0]'; C=[1 1 1]; Qo=[C;C*A;C*A^2]; P=inv(Qo); n=length(A) r=rank(Qo); if r==n disp('System is Observable.')
303

r
Qo P Ao1=Qo*A*P Bo1=Qo*B Co1=C*P else r end 运行结果: System is Observable. r= 3
304

Qo = 1 1 1 6 6 5 23 22 17 P = 8.0000 -5.0000 1.0000 -13.0000 6.0000 -1.0000 6.0000 -1.0000 0.0000 Ao1 = 0 1.0000 0.0000 0.0000 0.0000 1.0000 4.0000 -8.0000 5.0000 Bo1 = 0 0 1 Co1 =1.0000 0.0000 0.0000
305

由于可观性阵的秩为3等于系统的阶次,因此可以判 断系统是完全可观的。依据第一可观规范型变换阵 的要求选取变换阵P,使系统实现按第一可观规范 型的变换。 与系统可控性判定相类似,MATLAB提供了应用 系统可观性判定函数 Obsv( ) ,使得进行系统可观 性判定更加简洁。

306

5.2系统的可控性分解与可观性分解 5.2.1 系统的可控性分解 对于系统的状态不完全可控的情况下,往往需要 对系统进行可控性分解。MATLAB提供了实现系统 进行可控性分解的函数。

若已知系统?(A,B,C),MATLAB提供的可控性分 解函数的基本格式为: [AB,BB,CB,T,K]=ctrbf(A,B,C) 其中,T为相似变换阵,K为可控子阵的阶次向量。 并存在相应变换关系满足:
307

~ =TAT’, A
~ ? Auc A?? ? A21 0? Ac ? ?

~ =TB, B
~ ?0? B?? ? ? Bc ?

~ C =CT’

~ C ? ?C uc

Cc ?

对系统进行可控性分解后得到相应的可控子系统,即 (Ac,Bc)为可控子对。 例 5-6 已知系统?(A,B,C)的系数阵分别为: ?0 ? ?1 2 ? 1? ?0 ? ?0 1 ? B?? ? C ? ?1 ? 1 1? A?? 0? ?1 ? ?1 ? 4 3 ? ? ? ? ?
308

试判断系统的可控性,若状态不完全可控则进行可 控性分解并求出可控子阵。 解:首先对系统进行可控性判断,根据前面介绍的 方法编写下述MATLAB程序。 A=[1 2 -1;0 1 0;1 -4 3]; B=[0 0 1]'; C=[1 -1 1]; n=rank(A); Q=[B A*B A^2*B]; r=rank(Q) if r==n disp('System is Controllable.') Q
309

else disp(‘System is Uncontrollable’) disp(‘rank’) r end 很容易得到系统的可控性判定结果: System is Uncontrollable rank r = 2 系统是不完全可控的,秩为2。因此就可以应用可控 性分解函数ctrbf( )对系统进行可控性分解:
310

[AB,BB,CB,T,K]=ctrbf(A,B,C) disp('controllable submatrix ') Ac=AB(2:3,2:3) Bc=BB(2:3,:) 得到: AB = 1 0 0 -2 1 1 -4 -1 3 BB = 0 0 1
311

CB = -1 -1 1 T= 0 1 0 -1 0 0 0 0 1 K= 1 1 0 controllable submatrix Ac = 1 1 -1 3 Bc = 0 1
312

5.2.2 系统的可观性分解
对于状态不完全可观测的系统,也可对其进行可观 性分解,得到可观子阵。 在MATLAB中提供了相应的函数obsvf( ),可 方便地进行可观性分解。 若已知系统?(A,B,C),则可观性分解函数obsvf( ) 的基本格式为: [AO,BO,CO,T,K]=obsvf(A,B,C) 其中,T为相似变换阵,K为可观子阵的阶次向量。 并存在相应变换关系满足:

~ A =TAT’,

~ B =TB,

~ C =CT’
313

~ ? Auo A?? ? 0

A21 ? ? Ao ?

~ ? Buo ? B?? ? ? Bo ?

~ C ? ?0

Co ?

对系统进行可观性分解后得到相应的可观子系统,即 (Ao,Co)为可观子对。 例5-7 已知系统?(A,B,C)的系数阵分别为:
1? ?1 ? ?? 2 C ? ?1 ? 1? B?? ? A?? ? ?0 ? ? 1 ? 2? 试确定系统的可观性并进行可观性分解。 解:首先进行可观性判定,编制下述程序并运行, A=[-2 1;1 -2]; B=[1 0]'; C=[1 -1];

314

Qo=[C;C*A]; n=length(A) r=rank(Qo); if r==n disp('System is Observable.') Qo else disp('System is Unobservable') disp('rank') r end 得到运行结果:
315

System is Unobservable rank r=1 显然系统是状态不完全可观的。因此可利用 MATLAB可观分解函数对系统进行可观性分解 。 [AO,BO,CO,T,K]=obsvf(A,B,C) disp('observable submatrix ') Ao=AO(2,2) Co=CO(1,2) 得到运行结果: AO = -1.0000 0.0000 0.0000 -3.0000
316

BO = 0.7071 0.7071 CO = 0.0000 1.4142 T = 0.7071 0.7071 0.7071 -0.7071 K= 1 0 observable submatrix Ao = -3.0000 Co = 1.4142
317

5.3 线性定常系统的综合 5.3.1 基于极点配置方法的系统设计 实际工程中常常提出一种要求就是被控对象是 已知的,其性能指标是预先给定的,要求设计者根据指 标要求设计控制器以使整个控制系统的性能达到设 计要求,这类问题就叫做系统的综合。 如果这一系统为线性定常系统,则这类问题就是 对线性定常系统的综合问题。综合的目的就是在系 统中引入适当的环节,使原系统的性能得到改进, 达 到一定的指标要求。
318

由于在系统的性能指标中最为重要的就是系统 的动态性能指标,而系统的动态响应一般是由闭环系 统的极点的位置所决定的,因此通过采用极点配置的 方法使系统的闭环极点配置在所期望的位置上,从而 达到一定的性能指标要求。 (一) 单输入系统的极点配置 极点配置问题就是通过对状态反馈阵的选择,使闭环 系统的极点配置在所期望的位置上. 系统通过状态反馈可以任意配置闭环极点的必要 条件是: 受控系统是完全可控的。

319

若受控系统 ?(A,B)其控制律取作线性状态反馈律: u=Rv-Kx (5-10) 其中v为1维的参考输入;x为受控系统的n维状态;R为 输入变换系数;K为n维状态反馈行阵, 即 K=[kn kn-1 ??? k1] (5-11) 其中ki为待求参数.若当已知系统的期望闭环给定极 点为: ?I (I=1,2, ???,n), 可以写出相应的闭环特征多 项式:
* * * * ? * ( s) ? ? ( s ? ?* ) = s n ? a1 s n ?1 ? a 2 s n ?2 ? ? ? ? ? a n ?1 s ? a n i
i ?1 n

(5-12)
320

受控系统 ?(A,B,)在控制律(5-10)式的作用之下,闭 环系统?c(A-BK,BR) 的状态方程可表为下述形式:

? x ? ( A ? BK )x ? BRv (5-13)
对系统进行综合的目标就是如何选取状态反馈 阵K使闭环极点配置在期望的位置上。对于一般的 ?(A,B,),很难看出如何选取K矩阵,如果先将?(A,B,) 经坐标变换T转化成第二可控规范型 ,其状态反馈 ~ 增益阵 K 的选取很直观。 设单输入系统是完全可控的,通过非奇异线性变 ~ ~ ~ 换T将?化为第二可控规范型?( A, B ) 。其中
321

? 0 ? 0 ~ ? A?? ? ? ? 0 ?? an ?

1 0 ? 0

0 1 ? 0

? ? ? ?

? an ?1 ? an ? 2 ?

0 ? 0 ? ? ? ? ? 1 ? ? a1 ? ?

?0 ? ?0 ? ~ ? ? B ? ??? ? ? ?0 ? ?1 ? ? ?

(5-14)

图5-1给出了含有状态反馈的系统结构图,其中实线部 ~ ~ ~ 分是可控规范型?( A, B )的结构,虚线部分表示线性状态 ~ ~~ 反馈部分u= Rv - Kx 的结构, ? 的每一个状态变量分别 经过一个放大环节反馈到输入端与输入信号加以综合, ~ 生成系统 ? 的控制输入信号。
322

v

R

?
- a1
~ ? k1

~ xn

?

~ xn?1

~ x2

?

~ x1

- a2
~ ? k2

- an

~ ? kn
图5-1 含有状态反馈的系统结构图
323

显然由图5-1可以看出,反馈回路中的放大环节可以 加以合并,得到: ~, 因此闭环系统 ~ ? (ai ? k i ) 仍然是可控规范型。 ~ ~~ ~ ? c ( A ? B K , BR) 其系数阵为:
0 1 0 ? 0 0 1 ~ ~~ ? ? ? ? A ? BK ? ? 0 0 0 ? ~ ~ ~ ?? a n ? k n ? a n ?1 ? k n ?1 -a n ? 2 ? k n ? 2 ? ? 0 ? ? 0 ? ? ? ? ? 1 ~? ? - a1 ? k1 ? ?

(5-15) ~ ~ ~~ ~ 则闭环系统 ? c ( A ? BK , BR) 的特征多项式为:

~ ~~ ? ( s) ? det(sI ? A ? BK ) = ~ ~ ~ n ? (a1 ? k1 ) s n ?1 ? ? ? (a n ?1 ? k n ?1 ) s ? (a n ? k n ) s (5-16) 324

令其与希望闭环特征多项式相等,则可以得到:

~ k i ? ai* ? ai

(i=1,2, n)

(5-17)

从而可以得到可控规范型的状态反馈阵:
~ * * * K ? [a n ? a n a n?1 ? a n ?1 ? a1 ? a1 ]

(5-18)

325

据此可以将单输入系统闭环极点配置的基本步骤归 纳如下: (1)判断?(A,B,)的完全可控性,确定能否完成预定的 闭环极点配置综合目标。 (2)根据给定的动态指标或闭环极点确定希望的闭环 特征多项式的各个系数 a i*. (3)确定并求取开环系统的特征多项式系数ai. ~ ~ ~ (4) 按(5-18)式确定?( A, B ) 的状态反馈阵中的诸元素. (5)确定将?化为第二可控规范型的变换阵T及其逆阵, 变换阵T可取下述形式:

326

T ? [ A n ?1 B ? AB

0 ? ? 1 ?a 1 ? B ]? 1 ? ? ?? ? ? a n ?a1 1? ?

(5-20)

(6)确定?实现所希望的极点配置的状态反馈阵K. (7)依据对系统静态指标的要求确定系统输入变换 阵R 也可以采用其他方法对系统进行极点配置,如根 据希望特征多项式和在控制律作用下的闭环特征 多项式系数相等的关系,联立求解确定状态反馈阵 等方法。 下面通过示例说明如何应用MATALB以第二可 控规范型确定状态反馈阵实现系统的极点配置
327

例5-8 已知系统?(A,B,)的状态方程为:

?3 3 ? x ? ??42? 9? x ? ?1? u ? ? ? ? ? ? ? ?
通过状态反馈使希望极点分别配置在 ? 1=–1+j2 和 ?2=-1-j2 处,试确定状态反馈阵K。 解:根据极点配置的原理与步骤,利用MATLAB 语言可以十分方便的实现极点配置。 (1)首先判断系统的可控性,利用 ctrb( )求取 系统可控阵的秩,从而即可判定系统的可控性。 下面给出了采用MATLAB求取系统可控阵秩的程 序和结果。
328

A=[-2 -3;4 -9]; B=[3;1]; C=[0 0]; D=[0]; sys=ss(A,B,C,D); Co1=ctrb(sys) Co=rank(Co1) 得到运行结果: Co1 = 3 -9 1 3 Co = 2
329

由运算结果显然可以看出系统是完全可控的,因此 可由状态反馈任意配置闭环系统的极点。 (2)确定闭环系统的希望特征多项式(系数)。 系统的希望极点分别为:?1=–1+j2 和 ?2=-1-j2, 故可以定义W阵:
?? 1 ? j 2 0 ? W ?? ? 0 ? 1 ? j 2? ?

在MATLAB工作空间建立W阵,利用poly( ) 函数即可求出对应的特征多项式(系数)。 W=[-1+j*2 0;0 -1-j*2]; WA=poly(W)
330

WA= 1

2

5

(3)确定开环系统的特征多项式系数,由于可在 poly( )求取开环系统的特征多项式系数。 AS=poly(A) 从而得到运行结果: AS = 1 11 30 ~ ~ ~ ~ (4)根据(5-18)式确定 ?( A, B ) 的状态反馈 K 阵 中的诸元素。 (5)确定将?化为第二可控规范型的变换阵T及 其逆阵。
331

(6) 确定? 实现所希望的极点配置的状态反馈K 显然上述这三步都是相关的,首先确定变换阵T并 求出其逆阵T-1,然后再根据(5-19)式求出状态反 馈阵K。利用前面求出的结果使用MATLAB可以 方便地得到极点配置的状态反馈阵。下面给出了 相应的程序和结果。 T1=[A*B B] TB=[1 0;AS(2) 1] T2=T1*TB T=inv(T2) K=[WA(3)-AS(3) WA(2)-AS(2)]*T
332

运行结果: T2 = 24 3 14 1 T= -0.0556 0.1667 0.7778 -1.3333 K= -5.6111 7.8333 其中K就是用以实现闭环极点配置所需要状态反馈 增益阵。 可以验证选取状态放馈增益阵K,可以使系统闭 环极点被配置在希望位置。 333

AT=A-B*K; BT=B; CT=[1 0]; DT=0; [z,p,gain]=ss2zp(AT,BT,CT,DT) 运行结果得到: z = -8.0000 p = -1.0000 + 2.0000i -1.0000 - 2.0000i gain = 3
334

显然,引入状态反馈后,闭环极点确实被配置在希 望位置,达到预期目标。 例 5-9 已知系统如图5-2所示,试设计该系统的状 态反馈矩阵K,使闭环系统满足下述动态指标: (1) 述出超调量 ? ? 5%; (2) 峰值时间 t ? 0.5s; (3) 静态位置误差为零

u

G1(s)

G2(s)

G3(s)

y

图 5-2 系统结构图
335

1 1 1 其中G1(s) = G2(s) = G3(s) = s s ? 12 s?6

解:(1)将综合性指标转化为系统的希望极点 由于系统无开环零点,所以闭环系统的动态 特性完全由闭环极点决定。希望的三个闭环极 点可以这样安排:选择一对主导极点?1和 ?2; 另一个极点?3远离虚轴,后者对闭环系统的动态 性能影响很小。 2 主导极点?1, 2= ? ?? n ? j 1 ? ? ? n 利用二阶系统超调量与峰值时间的计算公式可 以解出, ? ? 0.707 , ? n ? 9
336

若取 ? ? 0.707 , ? n =10, 则主导极点为: ?1, 2= -7.07 ? j7.07 选取 ?3= -100 根据特殊确定的三个极点, 定义W 阵,再求取闭 环特征多项式系数。 W=[-7.07+j*7.07 0 0;0 -7.07-j*7.07 0;0 0 -100] WA=poly(W) 从而得到运行结果 WA = 1.0e+003 * 0.0010 0.1141 1.5140 9.9970 所以可以得到希望的闭环特征多项式为 :
337

? * ( s) ? s 3 ? 114 .1s 2 ? 1514 s ? 9997
(2) 根据静态指标确定系统输入变换系数R。 考虑到原系统无开环零点,静态位置误差要求为零 ,因此存在 (3)根据系统给定的传递函数,可以容易地确定 它所对应的开环系统的状态空间表达式。下面给 出了相应的程序和运行结果。 num=[1];den=[1 18 72 0]; [A1,B1,C1,D1]=tf2ss(num,den) G1=ss(A1,B1,C1,D1);
338

R e p ? lim G L ( s ) ? * ? 1 s ?0 a3

* R ? a3 ? 9997 ? 10000

N=fliplr(eye(3)); G=ss2ss(G1,N); [A,B,C,D]=ssdata(G) 上述程序运行得到: A = 0 1 0 0 0 1 0 -72 -18 B = 0 0 1 C = 1 0 0 D = 0
339

(4)求取化系统为第二可控规范型变换阵的逆矩阵 T-1。 (5)确定实现希望系统闭环极点配置的状态反馈阵。 由(5-18)和(5-20)式可以得到化系统?为第二 可控规范型变换阵T及其逆阵T-1与实现希望系统闭 环极点配置的状态反馈增益阵K。下面给出了相应 的程序与运行结果。 W=[-7.07+j*7.07 0 0;0 -7.07-j*7.07 0;0 0 -100] WA=poly(W)

340

AS=poly(A) T1=[A^2*B A*B B] TB=[1 0 0;AS(2) 1 0;AS(3) AS(2) 1] T2=T1*TB T=inv(T2) K=[WA(4)-AS(4) WA(3)-AS(3) WA(2)AS(2)]*T 运行上述程序得到: T2 = 1 0 0 0 1 0

0

0

1
341

T2 = 1 0 0 T= 1 0 0 K= 0 1 0 0 0 1 0 1 0 0 0 1

1.0e+003 *
9.9970 1.4420 0.0961
342

(二)多输入系统的闭环极点配置 1 将完全可控的多输入系统化为对单输入完全可控 系统 对于一个完全可控的多输入系统 ?0(A, B) 通过状 态反馈的方法可以化成对某个输入分量完全可控。 对于系统的一般描述 ?0(A, B),如果欲寻找实现上 述目标的状态反馈矩阵K是十分困难的,因此可以考 ? ? ? 虑采用 ?0 的规范型 ? 0 ( A, B) ,对规范型寻找实现 相同综合目标的状态反馈阵应当是容易的。再由 化 规范型的变换阵T,就很容易确定实现原系统综合目 标的状态反馈阵K。
343

下面给出将完全可控的多输入系统化为对单输入完 全可控状态反馈阵K的求取方法和步骤。 (1)首先判断系统的可控性。 (2)确定化为第一可控规范型的变换阵T。
在可控性矩阵中挑选n个线性无关的列向量并确定 指数集,{?1,?2,??? ?r},存在?1 + ?2 + ?r = n,且r ? l, 根据此方案选出的n个列向量构造成变换阵T。
[b1 Ab1 ? A?1 ?1b1 ? ? ? br Abr ? A? r ?1br ] T=

应用变换阵T则可将?0(A, B)化成第一可控规范型。
344

? ? ? (3)确定将完全可控的多输入系统规范型 ? 0 ( A, B) 化为对第一个输入分量w1的完全可控系统的状态反 ? 馈阵 K 。 ? ? ? ? 0 ( A, B) 如果将完全可控的多输入系统规范型 化 为对第一个输入分量w1的完全可控系统,其对应的状 ? 态反馈阵K 应具有下述特点: ? 1)状态反馈阵K 为l ? n阶阵, l是输入量的维数,n是系 统的阶次。 2)其第一行为全零行。 3)非零元素“-1”出现的位置是,第2行?1列、第3 行?1+?2 列、??????,直至到第r行?1+?2+ +?r-1-列。 4)最后的?r列都是全零列。
345

则状态反馈阵可以表示成下述形式:
? K ? ?[ 0 ? 0 e 2 ? 0 ? 0 e 3 ? ? ? 0 ? 0 e r ? 0 ? 0 ] ? ? ? ? ?? ? ?? ? ? ? ?? ? ? ? ?
? 1列 ? 2列 ? -1列 ? r列

(5-21)

其中“0”均表示该列为零向量。 (4)根据原受控系统?0(A, B)的状态阵满足 ? (5-22) K = K T-1 关系确定K阵。 原受控系统?0(A, B)引入状态反馈律u= -Kx,必 可将原系统化成对输入第一分量w1完全可控的闭环 系统 ?K(A-BK, B) 。
346

2.实现多输入系统极点配置的一类方案
在此仅介绍基于化多输入系统为第一可控规范 型实现多输入变为单输入完全可控系统的方法。 此方法可分为两步进行,第一步通过特定的状 态反馈阵将多输入系统化为单个输入是完全可控的。 第二步对单输入完全可控系统,按希望的极点要求 选择它的状态反馈阵。 第一步的设计方法在前面部分已经给出了基本步 骤,第二步的设计方法在单输入系统的极点配置方 法中已经做了详细的讨论。

347

多输入系统的闭环极点配置的综合步骤如下: (1)检验?0(A, B)的可控性,如果它是状态完全可 控的,则可以任意配置闭环极点。 (2)完全可控的?0(A, B),按照将完全可控的多输 入系统化为对单输入完全可控系统的方法和步骤构 ~ ? ? K ,并令 K = K T-1,则在控制 造变换矩阵T和反馈阵 ~ ~ ~ ~ ? 律u= w - K x作用下的闭环系统 x =(A-B K)x+Bw , ~ ~ w1 1) 由第一个输入分量 K 是完全可控的。即(A-B ,b 是完全可控的。 ~ ? ? x =(A-B x )x+b K 进行极点配 (3)单输入系统 1 ~ =w - ? x,使闭环系统w =(A~ 置,选择控制律w1 1 1 k1 ~ -b ? )x+b w 的极点配置在希望的位置上。 B K 1 k1 1 1
348

~ ~ (1)其中 k1 是n维行向量,w1是l维列向量 w 的第 一个元素,w1是l维列向量w的第一个元素。 ? 如果将 k1 扩写为n? l矩阵,即
?
? ? k1 ? ? ? K ?? 0 ? ? ?0? ? ?

(5-23)

? 并且有等式BK =b1 k1 成立,那么上述控制律可以扩 ~ ~ w = w- K x,闭环方程可以写成 x =(A – BK - BK) ? 写为 x+B w 。 ~ (4)两步的状态反馈阵 K 和 K 相加就是原系统 ?0(A, B)实现闭环极点配置的状态反馈阵,即
349

K= K + K = K +

~

~

? ?k1 ? ?0? ? ? ? ?0? ? ?

(5-24)

例 5-10 已知系统的状态方程如下:
?1 ? x?? ?0 0? ?1 ? x ? ?0 1? ? 1? ?u 1?

将系统闭环极点配置在?1= -1 和 ?2= -2 处。 解:(1) 首先判断系统的可控性。 容易判断该系统可控阵的秩为2,等于系统阶次, 所以系统是完全可控的,可以任意配置闭环系统的 极点,在可控阵中选择线性无关的列向量,就是?1 =?2=1。

350

~ (2)确定状态反馈阵 K

由于变换阵T=B,所以可以求出其逆阵T-1 T-1 =
?1 ?0 ? ? 1? 1? ?

根据(5-12)式即可确定在规范型下化为由单输入 完全可控的状态反馈阵。

0? ?0 ? K = ?? 1 0? ? ? 由(5-13)是可以得到

~ ?0 K ?? ?- 1

0? ? 1?

351

所以得到闭环系统:

~ ~ ? x =(A-BK)x+B w =
从而得到单输入系统

?2 ?1 ?

? 1? x ? ?1 0? ?0 ? ?

1? w ~ 1? ?

? x ? ?2 ?1 ?

? 1? x ? ?1? w ~ 0? ?0? 1 ? ? ?

是完全可控的。由此即将多输入系统的极点配 置问题转化为单输入系统的极点配置。 (3)确定单输入系统闭环极点配置的状态反 ? k1 馈阵
352

根据已知希望闭环极点可以得到希望特征多项式:

? * ( s ) ? s2+3s+2
根据单输入系统的系数阵可以方便地计算出开环 系统的特征多项式:
? (s ) ? s2-2s+1
? k1 因此可以得出状态反馈阵 (过程可参见例5-8)

k1? =[5 1],

(4)确定原系统闭环极点配置的状 态反馈阵K ? 将 k1 扩写成2?2的多输入系统状态 反馈阵

353

? ?k1 ? ?5 K ?? ?? ? ? 0 ? ?0

1? 0? ?

从而可以求出原系统实现闭环极点配置目标的状 态反馈阵为:

~ K= K ? K ? ? 0 ?- 1 ?

0 ? ? ?5 1? ?0 ? ?

1? ? ? 5 0? ?- 1 ? ?

1? 1? ?

(5)校验 原系统闭环系数阵为:

?1 0? ? ?1 1? ? 5 1? ? ?? 3 ? 2? A-BK= 0? ?0 1? ?0 1? ?- 1 1? ? 1 ?? ? ? ? ? ? ?
354

由此可以得到其相应的特征多项式为:

? (s) ? s2+3s+2= ? * ( s )
可以看出系统极点已经被配置在给定的位置,实现 多输入系统的闭环极点配置目地。 (三)多输入系统的闭环极点配置状态反馈阵的 不 唯一性 如果设原系统的状态反馈阵为: K=

? a1 a 2 ? ?a 3 a 4 ? ? ?
355

则得到形如下述形式的闭环特征多项式,记作:

? (s ) ? s2+(a1+a3+a4-2)s+(a1a4-a2a3-a1-a3-a4+1)
? * ( s ) ? s2+3s+2

(5-25)

若设根据希望闭环极点所得到的闭环特征多项式为:

(5-26)

并令 ? ( s ) ? ? * ( s ) 从而可以得到 a1+a3+a4=5 a1a4-a2a3=6 显然其解不是唯一的, 但是所配置的闭环极点位置应 该是相同的。
356

5.3.2 基于MATLAB acker( )与place( )函数的系统闭 环极点配置的实现 在第一部分中介绍了在控制系统采用状态空间设 计方法中应用较多的极点配置法。 在应用MATLAB进行基于状态空间方法的控制 系统设计中,不但可以采用前面的方法,还可以直 接使用MATLAB自身提供的有关闭环极点配置函数。 在此通过示例介绍基于acker( )和 place( )函数的 系统闭环极点配置方法。

357

acker( )函数是基于Ackermann算法进行反馈增益 阵K的求解,该函数只适用于单输入单输出系统(SISO系 统),其基本格式为: K=acker(A,B,P) --------其中A、B是系统系数阵;P为闭环系统的 希望极点向量;K即为所求状态反馈控制律u = -Kx中 的反馈增益阵。 Place( )函数既可应用于单输入系统又适用于多输 入系统,其基本格式为: K=place(A,B,P) [K,prec,message]=place(A,B,P) 其中 A,B,P,K与acker( )函数中的意义相同, 返回变 358 量prec 为实际极点偏离希望极点位置的误差。

例5-11 已知连续控制系统,其系数阵为:
?0 1 0 ? A ? ?0 ? 12 1 ? ?0 1 6 ? ? ?
?0 ? B ? ?0 ? ?1 ? ? ?

C =[1 0 0],

D=0

试设计该系统的状态反馈阵K,使闭环系统满足下 述性能指标要求。 (1)输出超调量不大于5%; (2)峰值时间不大于0.5秒; (3)静态位置误差为零。 解: 1 对原系统动态响应进行分析
359

首先考查原系统动态响应情况,分析原系统与所要求 的动态指标存在的差异。编制MATLAB程序绘制出 原系统单位阶跃响应以及求出零极点分布情况。 A=[0 1 0;0 -12 1;0 1 6]; B=[0;0;1]; C=[1 0 0]; D=0; sys=ss(A,B,C,D); [Z,P,G]=ss2zp(A,B,C,D,1) step(sys) 程序运行得到系统零极点情况:
360

Z= Empty matrix: 0-by-1 P= 0 -12.0554 6.0554 显然原系统存在一个正实数极点(6.0554)因此系 统是不稳定的,这一点也可从图5-3 给出的系统阶跃 响应中看出。因此根本无法满足系统性能指标要求, 必须对系统加以校正设计。
361

Step Res pons e 0.18 0.16 0.14

Amplitude

0.12 0.1 0.08 0.06 0.04 0.02 0 0 0.16 0.32 0.48 0.64 0.8

Tim e (s ec.)

图 5-3

系统阶跃响应
362

2. 将综合性指标转化为系统的希望极点并确定输入 变换系数R 根据例5-9中的设计思想,将闭环系统的极点设 置在:?1, 2= -7.07 ? j7.07;?3= -100 故取系统闭环 极点向量为:
P=[ -7.07 + j7.07 -7.07 - j7.07 –100] 由于与例5-7具有相同的希望极点和静态指标, 本例 中的输入变换系数R=10000。 3. 对系统进行可控性判定,并应用acker( )函数求取 状态反馈增益阵K。 对系统进行状态反馈校正必须判定系统可控性,只有 系统状态完全可控才能利用状态反馈实现希望的极 363 点配置。

本系统虽然是不稳定的,但经系统状态的可控 性判定可知它的状态是完全可控的。因此可应用 acker( )函数求取状态增益阵K。 A=[0 1 0;0 -12 1;0 1 6]; B=[0;0;1]; C=[1 0 0]; D=0; P=[-7.07+j*7.07 -7.07-j*7.07 100]; n=length(A); sys=ss(A,B,C,D); Con=ctrb(sys); r=rank(Con);
364

if r==n disp('The system is controllable') Con disp('Feedback gain matrix K') K=acker(A,B,P) else disp('The system is uncontrollable') end 得到运行结果: The system is controllable
365

Con = 0 0 1 0 1 -6 1 6 37 The feedback gain matrix K K= 1.0e+003 * 9.9970 0.2893 0.1081

根据希望极点得到了状态反馈增益阵K,从而也 就能够达到系统性能指标要求。

366

4. 进行系统校验并考查校正后系统的阶跃响应 根据已经确定的系统反馈增益阵构建闭环系统, 同时应考虑输入变换系数R的作用。下面给出了相 关程序,图5-4给出了校正后的闭环系统的阶跃响应, csys=ss(AA,B,C,D); disp('The closed-loop poles') POLES=pole(csys) t=0:0.001:2; u=10000+0*t; lsim(csys,u,t)

367

下面给出的验证结果 The closed-loop poles POLES = 1.0e+002 * -0.0707 + 0.0707i -0.0707 - 0.0707i -1.0000 从图5-4闭环系统的阶跃响应也可看出,系统的 超调量与峰值时间均满足系统性能要求。 可以看出,对于一个不稳定的系统通过状态反馈 进行希望位置上的极点配置,可以实现预期的性能指 标。
368

系统设计与综合变得十分方便。
1.4 1.2 1 Linear Results Simulation

Amplitude

0.8 0.6 0.4 0.2 0

0

0.5

1 Time (sec.)

1.5

2

图 5-4 校正后闭环系统的阶跃响应

369

例 5-12 已知多输入系统,其系数阵为:

?0 0 5 ? A ? ?1 0 ? 4 ? ?0 1 16? ? ?
C=[1 0 0],

?2 0? B ? ?1 2 ? ?0 1 ? ? ?
D=0

试采用状态反馈方法将系统闭环极点配置在 ?1, 2= -7.07 ? j7.07;?3= -100 处。 并绘制出系统的阶跃响应,其 中u1=10,u2=5。 解:对于多输入系统进行基于状态反馈的极点配置, 可以直接应用MATLAB的 place( )函数进行。
370

1

对原系统动态响应进行分析 首先考查原系统动态响应情况,编MATLAB程 序绘制出原系统单位阶跃响应以及求出零极点 分布情况。 A=[0 0 5;1 0 -4;0 1 16]; B=[2 0;1 2;0 1]; C=[1 0 0]; D=[0 0]; sys=ss(A,B,C,D); [Z,P,G]=ss2zp(A,B,C,D,1) t=0:0.01:2; u=[10+t*0;5+t*0]; lsim(sys,u,t)
371

解出系统的零极点情况发现3个极点均位于s平面的

右半侧。 Z= 0.4171 15.5829 P= 0.1168 + 0.5509i 0.1168 - 0.5509i 15.7664 G= 2
372

8 7 6 5 4 3 2 1 0

x 10 11

12 10 8 6 4 2 0

x 10 12

0

Time a(sec.)

2

0

Time (sec.)
b

2

图5-5 原系统分别在u1(a)、u2(b)作用下的阶跃 373 响应

图5-5给出了系统分别在两个输入作用下的阶跃响应, 显然原系统是不稳定的。
2 应用place( )函数求取状态反馈增益阵K。 由于该系统为多输入系统,在进行状态反馈极 点配置时需应用place( )函数,首先应对系统可 控性进行判定,然后再求取反馈增益阵K。 分别考查对应于两个输入的系统可控性。 A=[0 0 5;1 0 -4;0 1 16]; B=[2 0;1 2;0 1]; C=[1 0 0];

374

D=[0 0]; Con1=ctrb(A,B(:,1)); Con2=ctrb(A,B(:,2)); n=rank(A); r1=rank(Con1); r2=rank(Con2); if n==r1& n==r2 disp('The system is controllable') Con1 Con2 else r1
375

r2 end 两个可控性矩阵都是非奇异的,因此可断定由 每个输入都可控制系统的全部状态分量。下面给出 程序运行结论以及相应的可控阵。 The system is controllable Con1 = 2 0 5 1 2 -4 0 1 18 Con2 = 0 5 90 2 -4 -67 1 18 284 376

由此即可对系统按期望的闭环极点进行校正, 确定状态反馈增益阵K,然后求取其阶跃响应,并 加以验证。下面给出相应的程序及执行结果。 P=[-7.07+j*7.07 -7.07-j*7.07 K=place(A,B,P) csys=ss(A-B*K,B,C,D); t=0:0.01:4; u=[10+t*0;5+t*0]; lsim(csys,u,t) disp('The closed-loop poles') poles=eig(A-B*K)
377

-100];

程序运行结果: K= 39.9708 20.2330 -11.6702 0.7133 -0.8689 31.7033 The closed-loop poles poles = -100.0000 -7.0700 + 7.0700i -7.0700 - 7.0700i 经过校正将闭环系统的极点配置在期望位置上, 图5-6给出了校正后的系统阶跃响应。
378

1.2 1 0.8 0.6 0.4 0.2 0 0 1
Time (sec.)

2

3

4
379

图 5-6 校正后系统阶跃响应

另外, acker( )函数只是用于单输入系统,而place( )函 数不但适用于多输入系统,同时也适合单输入系统。 例如在例5-11中如果选择 place( ) 来实现状态 反馈增益阵的求取,会得到与例5-11中同样的结果。 Feedback gain matrix K K= 1.0e+003 * 9.9970 0.2893 0.1081 The close loop poles POLES = -7.0700 + 7.0700i -7.0700 - 7.0700i -100.0000
380

5.3.3 具有观测器的状态反馈控制系统设计与实现 (一) 状态观测器 为了实现状态反馈需要知道所有状态信息,但是系统 的所有状态不一定都能够直接测量到,为了解决这一 困难,提出了状态重构问题。 状态重构问题的核心,就是重新构造一个系统?obv, 利用原系统可直接测量的输入和输出量作为它的输入 ? x(t ) 信号,并使它的输出信号 在一定的指标下与原系统 ?的状态向量x(t)等价。 ? 通常 x(t )叫做x(t)的估计状态或状态重构,而把实现 转状态重构的系统?obv称作观测器。 等价性指标常常满足下述渐进指标:
381

~ lim x (t) = lim [ x(t)- x(t )]=0 ? t ??
t ??

(5-27)

? = x(t)- x(t ) 叫做观测误差。如果观测器?obv 的维数与原系统?的维数相同,就称其为全维观测器; ?obv的维数少于原系统?的维数则称其为降维观测器。 根据图5-7 观测器结构图可以写出重构状态部分 的方程如下: =Ax*+Bu+H(y-y*) ? x * (A-HC)x*+Bu+Hy (5-28)
382

~ 其中 x

? x ? Ax ? Bu
y=Cx (5-29) ?=(A,B,C) H B ? x* C x* y*

u

y

A
图5-7 典型

观测器结构图
383

设(A,C)是完全可观测的,输入向量 u(t) 与输出 向量 y(t) 都可以直接测量。由(5-29)式减去(5-28)式, 得到 ? ? (5-30) x ? x * ? ( A ? HC)( x ? x * ) 只要 A-HC 是渐进稳定矩阵, x*(t)就是x(t)的重构状 ? 态 x(t )。因此通过选择 Hx? 矩阵,以达到要求的衰减速 度,其基本观测器方程为:

? ? x

=(A-HC) x +Bu+Hy ?

(5-31)

该观测器的输出就是重构状态

? x 本身。
384

从加快观测误差的衰减速率考虑,希望观测器 在复平面左半平面的极点越远离虚轴越好,这就希 望观测器应当有足够的频带宽度。另一方面观测器 的输入y(t)中不可避免地混有高频干扰,为了减小高 频分量,则希望观测器的频带应当尽可能地窄。 显然,快速性与抗干扰性的要求是一对矛盾, 只能根据具体情况加以折衷。 在应用MATLAB进行观测器的极点配置时仍可使 用place( ) 和 acker( ) 函数。其基本格式为: H=acker(A’,C’,P) H=place(A’,C’,P) 其中 P为观测器的期望极点向量。
385

例5-13 已知系统 ?(A,B,C) 的系数阵为:

?1 0 0 ? A ? ?3 - 1 1 ? ?0 2 0 ? ? ?
C=[0 0 1],

? 2? B ? ?1 ? ?1 ? ? ?
D=0

设计观测器,使其极点配置在希望极点向量 Po[-30,40,-50] 处。 解: 应用极点配置法设计观测器,首先要判断系统的可 观性,然后对系统进行极点配置,求出反馈增益阵H, 实现观测器设计。下面给出相关程序。
386

A=[1 0 0;3 -1 1;0 2 0]; B=[2;1;1]; C=[0 0 1]; D=0; Po=[-30,-40,-50]; sys=ss(A,B,C,D); Obv=obsv(sys); ro=rank(Obv); n=length(A); if ro==n disp('The system is observable') Obv disp('The gain matrix of observer is') H=acker(A',C',Po) disp('The model of observer is') Ao=A-H'*C 387 Bo=B

else disp(‘The system is unobservable’)
end

运行结果得到观测器增益阵以及所构造的观测器模 型: The system is observable Obv = 0 0 1 0 2 0 6 -2 2 The gain matrix of observer is H= 1.0e+004 * 388 1.0804 0.2352 0.0120

The model of observer is Ao = 1.0e+004 * 0.0001 0 -1.0804 0.0003 -0.0001 -0.2350 0 0.0002 -0.0120 Bo = 2 1 1 (二)降维观测器设计 设受控系统?(A,B,C)完全可观,且系数阵C存在rank (C)=m。如果某些状态分量可以由输出 y各分量
389

的线性组合构成,这样要求重构的状态分量个数将 少于n个,即观测器的维数低于n维,这就称作降维 观测器。 设系统?(A,B,C)有m个彼此独立的输出分量,应 该可以组合出m个状态分量,需要重构的只剩下n-m 个状态分量。为此可把状态向量按照输出做坐标变 换,使前m个分量等于输出,即

~ ? x1 ? ? y ? ?Cx ? ? C ? ~ x ? ? ? ? ? ? ? ? ? ? ? ? ? ??? x ~ ~ ? x 2 ? ? x 2 ? ?Gx ? ? G ? ? ? ? ? ? ? ? ?

(5-32)
390

坐标变换阵的逆阵为:

?C ? T ?1 ? ??? ?G ? ? ?

(5-33)

其中C是系统的系数阵,秩为m;G为使矩阵T-1非奇异 而任意选择的n-m个行向量组成的(n-m)?n矩阵。G 是不唯一的。为了保证C 阵的分块要求,G 阵可取作 G = ?0 ? I n?m ?形式,则存在

T ?1

?C1 ? C 2 ? ? ?? ? ? ? ? 0 ? I n?m ? ? ?

(5-34)
391

可以证明相应的变换阵是

?C1?1 ? ? C1?1C 2 ? T ??? ? ? ? ? 0 ? I n?m ? ? ?

(5-35)

显然系数阵C 经过坐标变换后得到:
~ C ? CT ? ?I m ? 0?

(5-36)

则得到量测方程为:

~~ y ? Cx ? ?I m

~ ? x1 ? ~ ? 0?? ? ? ? x1 ~ ? x2 ? ? ?

(5-37)
392

因此只要重构n-m维分状态。这时观测器的维数就 是n-m。一般地,降维观测器的最小维数是nrank(C) 维,然而在特殊情况下观测器的维数还可以 进一步降低。 设完全可观的受控系统?(A,B,C),满足 rank(C)=m,经过坐标变换系统化为: ~ ? A11 ? A12 ? ? x1 ? ? B1 ? ~ ? Ax ? Bu ? ? ? ? ? ? ? ? ? ? ? ? ? u ? ~~ ~ x ~ ? A21 ? A22 ? ? x 2 ? ? B 2 ? ? ?? ? ? ?

~ y ? x1

(5-38)
393

因为坐标变换不会改变系统的可观性,因此 ? 仍然 ~ ~ x 2 ,现将包含x 2 的子系 是完全可观测的。为了重构 统方程列写成如下形式:

~

~ ? A x ?B u? A y ~ ? x2 22 2 2 21

~ ω? A12 x 2

(5-39)

在此状态方程中可以将B2u+A21y 作为输入,子系统 的输出是n-m维向量 ω 。可以用另以子系统表示 , 即ω ? (5-40) ω? y ? A11 y ? B1u
394

将(5-40)式代入(5-39)式得到对(5-39)式子 系统设计的基本观测器:

? ? z ? ( A22 ? HA12 )z ? ( B2 u ? A21 y ) ? H ( y ? A11 y ? B1u) (5-41)
~ x 2 的重构,H 矩阵是(n-m) 式中z 是n-m维状态分量

?m 维。考虑噪声的影响对(5-41)式进行如下变 换 : z ? z ? Hy ? (5-42) 得到:
? ? z ? ( A22 ? HA12 ) z ? ( B2 ? HB1 )u ? ( A21 ? HA11 ) y (5-43)
395

? z ? z ? Hy

(5-44)

(5-43)与 (5-44)式就构成了系统?(A,B,C)的降维观测 器模型。下面通过示例具体说明降维观测器的设计 步骤。 例5-14 已知系统?(A,B,C)的系数阵为:
4 4? ? 4 ?? 11 ? 12 ? 12 ?, A?? ? ? 13 14 13 ? ? ? ? 1? ? ? 1?, C ? ?1 B?? ? ? 0? ? ? 1?

1

试设计系统?(A,B,C)的降维观测器。
396

解:1 首先对系统进行可观性判定。运行下述程序

A=[4 4 4;-11 -12 -12;13 14 13]; B=[1;-1;0];C=[1 1 1];D=0; sys=ss(A,B,C,D); Obv=obsv(sys); ro=rank(Obv); n=length(A); if ro==n disp('The system is observable') Obv else
397

disp('The system is unobservable') end 可以得到: The system is observable Obv = 1 1 1 6 6 5 3 22 17 系统是完全可观的。由于系统是3维的,输出 是1维的,故可以构造2维降维观测器。 2 确定变换阵T 因为其逆阵T-1的第一行是C阵,其余两行保证T-1 非奇异即可,故得到:
398

T ?1

? C ? ?1 ? ??? ? ?0 ? G ? ?0 ? ? ?

1 1 0

1? 0? 1? ?

对其求逆得到:
?1 ? 1 ? 1? T ? ?0 1 0? 0 1? ?0 ? ?

~ x 根据变换关系: ? Tx 得到:

~ ?x ?x ~ ~ x1 ? x1 2 3
~ x2 ? x2

~ ~ ~ ?1 ?1 又由于存在变换: ? T AT , B ? T B, C ? CT ,故 A 可以得到变换后的系数阵为: 399

~ x3 ? x3

A~= 6 0 -1 -11 -1 -1 13 1 0 B~= 0 -1 0 C~= 1 0 0

3 确定子系统 根据(5-39)式可以得到降维观测器的子系统:
400

~ ~ ? ? x 2 ? ?? 1 ? 1? ? x 2 ? ?? 1? ? ? u ? ? ? 11? y ~ ? x ? ? 1 0 ? ? x3 ? ? 0 ? ~ ? ? 13? ? ?? ? ? ? ? ? ? 3?

ω? ?0

~ ? ? x2 ? 1?? ~ ? ? x3 ?

该子系统的特征多项式的根为: roots= -0.5000 + 0.8660i -0.5000 - 0.8660i
4 确定降维观测器模型 为了满足逼近速度与抗干扰能力的要求,在此将观测 器的极点配置在-3和-4处。应用place函数容易求出实现 将观测器极点配置在期望位置上的参数H值。由
401

A22=[-1 -1;1 0]; A12=[0 -1]; p=[-3 -4]; H=place(A22',A12',p 得到参数H为: H= -5 -6 从而由(5-43)式得到所设计的降维观测器模型:

? ? ?? 6 ? 6? ? z 2 ? ? ?? 1? u ? ?19 ? y ? z ? ? 5 ? 6? ? z 3 ? ? 0 ? ?49? ? ?? ? ? ? ? ?
402

(三) 具有观测器的状态反馈控制系统设计 设状态完全可控与可观侧的受控系统?(A,B,C)的 状态方程和输出方程重写如下:

? x ? Ax ? Bu
y=Cx

(5-45)

引入状态反馈控制律 u=y-Kx 构成直接状态反馈闭环控制系统:
? x ? ( A ? BK ) x ? Bv

(5-46)

y=Cx

(5-47)
403

当受控系统 ?(A,B,C) 的状态不能直接测量时,就需 ? 要引入状态观测器,重构状态 x(t )。重构状态反馈控 制律为: (5-48) ? u=y-K x 设状态观测器为全维基本观测器?obv,即

? ? =(A-HC) x +Bu+Hy ? x

(5-49)

考虑将(5-48)与(5-49)代入(5-45)式,可以得出重构状 态反馈控制系统 ?OB F 的状态方程和输出方程:
? ? x? ? A ? BK ? ? x ? ? ? B?v ? ? ? x ? ? HC A ? HC ? BK ? ? x ? ? B ? ?? ? ?? ? ? ? ? ?

y ? ?C

0?? x ? ? x? ? ??

(5-50)
404

显然,重构状态反馈控制系统 ?OB F 的维数是受控系 统 ? 的维数加上观测器 ?obv 的维数,因此其维数应 等于2n。 考虑到带有观测器的重构状态反馈系统 的分离 性质,因此先按闭环动态要求确定(A-BK)的极点, 从而设计出反馈增益阵K;再按重构误差的衰减速 度要求确定(A-HC)的极点,从而设计出H 阵。 图 5-8 给出了带有观测器的重构状态反馈控制 系统 ?OB F 的结构图。

405

例5-15 已知受控系统同 例5-11,试进行系统闭环期 望极点为?1, 2= -0.5 ? j*0.4;?3=- 80、具有状态观 测器的反馈控制系统设计。
y

u -

? x ? Ax ? Bu
y ? Cx

y

K

? x

? ? ? x ? ( A ? HC ) x

? Bu ? Hy

图 5-8 带观测器的状态反馈 控制系统结构图
406

解:关于状态观测器的设计已在例5-11中完成,为了 保证所设计观测器响应的快速性将极点配置在30、40、 50处。

根据分离性质,观测器的极点配置与状态反馈闭环 系统的极点配置可以分别加以进行。因此直接利用例 5-13 所完成的观测器设计结果,继续进行有关设计。 首先对系统进行可控性判定,如果系统状态完全 可控即可应用acker( )函数求取状态反馈增益阵 K, 以及相应的闭环系统模型(无状态观测器)。进而求 取具有状态观测器的反馈控制系统?OB F(Ac,Bc,Cc, Dc)。下面给出了相关程序和结果。
407

A=[1 0 0;3 -1 1;0 2 0]; B=[2;1;1]; C=[0 0 1]; D=0; Po=[-30,-40,-50]; Pcs=[-0.5+j*0.4,-0.5-j*0.4,-80]; sys=ss(A,B,C,D); n=length(A); Obv=obsv(sys); ro=rank(Obv); Con=ctrb(sys); rc=rank(Con) ospole=poles(sys)' if ro==n
408

同例 5-11内容 if rc==n disp('The system is controllable') Con disp('The feedback gain of the system') K=acker(A,B,Pcs) disp('The model of closed-loop system without observer') Ac1=A-B*K Bc1=B Cc1=C Dc1=D csys=ss(A-B*K,B,C,D); 409

disp('The poles of closed-loop system') Poles=pole(csys) figure(1) step(csys,30) disp('The model of closed-loop system with observer') Ac(1:3,1:3)=A; Ac(1:3,4:6)=-B*K; Ac(4:6,1:3)=Ho*C; Ac(4:6,4:6)=A-Ho*C-B*K; Ac Bc=[B;B] Cc=[C,0 0 0] Dc=0 csyso=ss(Ac,Bc,Cc,Dc);
410

figure(2) step(csyso,30) else disp('The system is uncontrollable') end 得到系统分析与设计结果: ospoles = 1 -2 1 The system is controllable Con = 2 2 2 1 6 2 1 2 12 The feedback gain of the system
411

K= 33.1361 -3.0794 17.8072 The model of closed-loop system without observer Ac1 = -65.2722 6.1589 -35.6144 -30.1361 2.0794 -16.8072 -33.1361 5.0794 -17.8072 Bc1 = 2 1 1 Cc1 = 0 0 1 Dc1 = 0
412

The poles of closed-loop system Poles = -80.0000 -0.5000 + 0.4000i -0.5000 - 0.4000i The model of closed-loop system with observer Ac = 1.0e+004 * 0.0001 0 0 -0.0066 0.0006 -0.0036 0.0003 -0.0001 0.0001 -0.0033 0.0003 -0.0018 0 0.0002 0 -0.0033 0.0003 -0.0018 0 0 1.0804 -0.0065 0.0006 -1.0839 0 0 0.2352 -0.0030 0.0002 -0.2368 0 0 0.0120 -0.0033 0.0005 -0.0138
413

Bc = 2 1 1 2 1 1 Cc = 0 0 1 0 0 0

Dc =
0

对原系统进行的分析可以确定该系统是不稳定的,因 为系统存在一个二重正极点 1。
414

图5-9与图5-10 分别给出了无观测器的状态反馈 极点配置闭环系统的阶跃响应和具有观测器实现状 态估计的极点配置闭环系统的阶跃响应。 在前面已经看出利用acker( )和 place( )函数可以 有效地实现观测器和系统的极点配置反馈设计,对 于采用降维观测器同样也可以实现闭环反馈系统。 另外,除了连续系统,对于离散采样系统同样 可以应用这些函数进行分析与设计。

415

0.3

Amplitude

0.25 0.2 0.15 0.1 0.05 0 0 10 Time 20 30
(s)

图5-9 无观测器极点配置 反馈系统阶跃响应
416

0.3

Amplitude

0.25 0.2 0.15 0.1 0.05 0 0 10 Time 20 30 (s)

图5-10 具有观测器的极点配置 反馈系统阶跃响应
417

5.4 线性二次型最优控制
5.4.1 线性二次型最优控制的一般概念 设对象是以状态空间的形式给出的线性系统,其目 标函数是对象状态和控制输入的二次型函数。 线性二次型最优控制一般涉及两类问题:一类是具 有状态反馈的线性二次型最优控制,即LQ问题,包括 拓展目标函数二次型最优控制;另一类就是线性二次 型高斯最优控制问题,即LQG问题。 设线性定常系统的状态方程为:

? x ? Ax ? Bu

(5-51)
418

J ? ?0 [ x T (t )Qx (t ) ? u T (t ) Ru (t )]dt (5-52)
其中Q和R分别为状态和输入控制变量的加权阵,Q是半 正定对称常数阵,R为正定对称常数阵。最优控制的目 标就是求取u(t),使性能指标J 达到最小值。可以导出 最优控制律u(t)为下述形式:
u(t)= -R-1BTP(t)x(t) 令K= R-1BTP(t) 则得到: u(t)= -K x(t) (5-53) (5-54) (5-55)
419

?

其中P(t)阵满足黎卡提(Riccati)方程: ATP+PA-PBR-1BTP+Q=0

(5-56)

下面给出二次型最优控制系统设计的基本步骤: (1) 求解黎卡提方程,确定P(t)阵。如果正定阵 P(t)存在(若A-BK是稳定矩阵,则总存在正定解P 阵),则系统是稳定的。 (2)将正定阵P(t)代入(5-54)式,得到最优反馈 增益阵K。 考虑目标函数中带有交叉乘积项的拓展二次型性 能指标: ? J ? ? [ x T (t )Qx (t ) ? u T (t ) Ru (t ) ? 2 x T (t ) Nu(t )]dt
0

(5-57)
420

其中N为交叉项的权阵,是正定对称常数阵。则系 统最优控制律为:u(t)= -K x(t) (5-58) 反馈增益阵K为: K= -R-1(BTP+N) (5-59) (5-60)

其中P 阵满足黎卡提方程: ATP+PA-(PB+N)R-1(BTP+N)+Q=0

对于具有输出权重的线性二次型目标函数,只要在 (5-52)式中用输出y代替状态x即可,其目标函数为

J ? ?0 [ y T (t )Qy (t ) ? u T (t ) Ru (t )]dt

?

(5-61)
421

更一般的目标函数可表示为下述形式:

J ? ?0 [ y T (t )Qy (t ) ? u T (t ) Ru (t ) ? 2 x T (t ) Nu(t )]dt (5-62)
5.4.2 线性二次型最优控制的MATLAB实现 MATLAB提供了直接进行线性二次型最优控制设计 的函数: lqr( )与lqry( )函数 Lqr( )函数用来求解由(5-52)和(5-57)式描 述的目标函数所对应的最优反馈增益阵K; lqry( )函 数用来求解由(5-48)和(5-49)式描述的目标函数对 应的最优反馈增益阵K。 它们的基本格式分别为:
422

?

[K,P,E] =lqr(A,B,Q,R,N) [K,P,E] = lqry(SYS,Q,R,N) 其中K为最优反馈阵,P为对应的黎卡提方程的正定 解,E为闭环的特征值。 例5-16 已知系统的系数阵为
1 0? ? 0 A?? 0 0 1? ?? 14? 11? 6? ? ?

?0 ? B ? ?0 ? ?1 ? ? ?

C=[1 0 0]

其性能指标为:

J ? ?0 ( x T Qx ? u T Ru ? 2 x T Nu )dt
423

?

选取

? 200 0 0? Q ? ? 0 1 0? ? 0 0 1? ? ?

?0 ? N ? ?0 ? ?1 ? ? ?

R=1

试设计二次型最优控制器,若系统控制信号取作: u=k1r-Kx 形式,在确定最优控制律时,假设输入信 号r=0,确定K,P,E。 解:在本系统的性能指标要求有交叉项,因此在使 用lqr( )函数时应保留交叉权阵N。 A=[0 1 0;0 0 1;-14 -11 -6]; B=[0;0;1]; C=[1 0 0]; D=0;
424

[K,P,E]=lqr(A,B,Q,R,N); disp('The feedback gain matrix is') K disp('The Riccati equation solution is') P disp('The closed-loop eigenvalues is')
E

sys=ss(A,B,C,D); Q=[200 0 0;0 1 0;0 0 1]; R=1; N=[0;0;1]; AC=A-B*K;BC=B*K(1);CC=C;DC=0; sysc=ss(AC,BC,CC,DC); [y,t,x]=step(sysc,8);
425

figure(1) plot(t,y) grid figure(2) plot(t,x) grid 得到运行结果: The feedback gain matrix is K = 5.8997 4.0597 1.5577 The Riccati equation solution is P = 145.6839 52.3970 5.8997 52.3970 30.9174 4.0597 5.8997 4.0597 0.5577 The closed-loop eigenvalues is

426

E= -5.4696 -1.0441 + 1.5963i -1.0441 - 1.5963I 图(5-11)和(5-12)分别给出了系统输出与状态 的阶跃响应。

427

0.4 0.3 0.2

0.6 0.4 0.2 0

0.1 0

-0.2 -0.4

0

2

4

6

8

0

2

4

6

8

图 5-11 系 统 输 出阶跃响应

图 5-12 系 统 状 态阶跃响应 428

例5-17 已知系统的二次型最优控制性能指标要求为
J ? ?0 [ y T (t )Qy (t ) ? u T (t ) Ru (t )]dt
?

系统数学模型同例5-13,权阵Q=1,R=1,如果系 统控制信号取作:u=k1r-Kx 形式,在确定最优控制 律时,假设输入信号 r=0,试设计二次型最优控制 器,求出K,P,E。 解: 在该系统的性能指标中出现了输出权重项,因 此求取反馈增益阵需应用lqry( )函数。 A=[0 1 0;0 0 1;-14 -11 -6]; B=[0;0;1]; C=[1 0 0];D=0; sys=ss(A,B,C,D); Q=1;R=1; 429

[K,P,E]=lqry(sys,Q,R); disp('The feedback gain matrix ') K disp('The Riccati equation solution ') P disp('The closed-loop eigenvalues ') E AC=A-B*K;BC=B*K(1);CC=C;DC=0; sysc=ss(AC,BC,CC,DC); [y,t,x]=step(sysc,8); figure(1) plot(t,y) grid 430

figure(2) plot(t,x) grid

程序运行结果: The feedback gain matrix K= 0.0357 0.0247 0.0041 The Riccati equation solution P= 0.7387 0.2717 0.0357 0.2717 0.1577 0.0247 0.0357 0.0247 0.0041 The closed-loop eigenvalues
431

E= -4.1666 -0.9187 + 1.5889i -0.9187 - 1.5889I

图 5-13 和图 5-14 分别给出了系统输出阶跃响应以 及系统状态阶跃响应。

432

3 2.5 2 1.5

x 10

3

6 4 2 0

x 10

3

1 0.5 0 0 2 4 6 8 -2 -4 0

2

4

6

8

图5-13 系统输出 阶跃响应

图5-14 系统状态 433 阶跃响应

对于采样离散系统,MATLAB也提供了一类与之相 类似的函数:dlqr( )与dlqry( )函数 5.4.3 线性二次型高斯(Gauss)最优控制 当系统和测量之中存在随机扰动(噪声)时,前 面所介绍的方法就无法应用了,这就涉及到线性二 次型高斯问题(LQG),在此将对这类问题加以介 绍并给出它的MATLAB实现。 设系统状态方程和量测方程分别为 :
434

? x ? Ax ? Bu ? Gw
y=Cx+Du+Hw+v

(5-63)
(5-64)

其中,w为系统噪声,v是由传感器引起的量测噪 声,w 与 v 均为零均值高斯分布的白噪声,统计 特性为: E(w)=0 E(v)=0 E(w wT)=Qn E(v vT)=Rn E(w vT)=Nn 线性二次型高斯控制的目标函数为:
J ? ?0 [ x T (t )Qx (t ) ? u T (t ) Ru (t ) ? 2 x T (t ) Nu(t )]dt
?

(5-65)
435

通过设计确定状态反馈控制律:

u= -Kx
以使目标函数J 达到最小。 图 5-15 示出了LQG控制系统框图。

(5-66)

436

w u

受控对象

y

+ -K
Kalman滤波器

+

? x
LQG控制器

yv

v

图5-15 LQG控制系统框图
437

(一)Kalman状态估计器的设计 在MATLAB中为用户提供了进行Kalman估计器 设计的Kalman( ) 函数。它们的基本格式为: [Kest,L,P]=Kalman(sys,Qn,Rn,Nn) [Kest,L,P]=Kalman(sys,Qn,Rn,Nn, sensors,known) 在sys中所包含着的系统系数阵形式为:(A,[B G],C,[D H]),返回变量Kest为卡尔曼状态估计器模 型,L为反馈增益,P是估计误差的协方差。 所得到的估计器Kest 利用已知的输入信号u和量 ? 测信号y 进行估计,得到x、y的优化估计 x 和 y , ?

其系统模型为:
438

? ? ? ? x ? Ax ? Bu ? L( y ? Cx ? Du)
? ? y ? ?C ? ? ? D? ? x? u ? x? ? I ? ?0? ? ? ? ?? ? ?

(5-67)

(5-68)

? y 其中,x 是系统的状态估计值, 是系统实际输出y的 ? 估计值。 例 5-18 已知系统的状态方程为:
? x ? Ax ? Bu ? Gw
y=Cx+Du+Hw+v

其中 :
439

?1? ?1 0 0 ? A ? ?0 2 1 ? B ? ?0 ? ?0 0 2 ? ?1? ? ? ? ?
?0.1? G ? ?0.1? ?0.1? ? ?

C =[1 1 0]

D=0

H=0

w为系统输入噪声,v为系统量测噪声, 试设计系统 的 kalman 滤 波 器 , 其 中 Q=1 , R=1 , Qn=0.01, Rn=0.0012 解: 应用MATLAB可以十分容易地完成对卡尔曼状 态估计器模型以及反馈增益L、估计误差的协方差P 的设计。下面给出了相应的程序: A=[1 0 0;0 2 1;0 0 2]; B=[1;0;1]; 440 C=[1 1 0];

D=0; G=[0.1;0.1;0.1]; Q=1;R=1;Qn=0.01;Rn=0.0012; BG=[B,G]; DD=[D,D]; sys=ss(A,BG,C,DD); [kest,L,P]=kalman(sys,Qn,Rn); 运行结果: Model of Kalman Estimator a= x1_e x2_e x3_e x1_e -17.62522 -18.62522 0 x2_e 8.53598 10.53598 1.00000

x3_e

-49.17648

-49.17648

2.00000
441

b= x1_e x2_e x3_e c= y1_e x1_e x2_e x3_e d= y1_e x1_e x2_e x3_e

u1 u2 1.00000 18.62522 0 -8.53598 1.00000 49.17648 x1_e x2_e 1.00000 1.00000 1.00000 0 0 1.00000 0 u1 0 0 0 0 0 u2 0 0 0 0

x3_e 0 0 0

1.00000

442

the estimator gain L and the steady-state error covariance L= 18.6252 -8.5360 49.1765 P= 0.2081 -0.1857 0.3663 -0.1857 0.1755 -0.3073 0.3663 -0.3073 0.7255

(二)LQG控制系统的实现 根据给出的系统数学模型可以设计出相应的卡尔 曼状态估计器以及系统最优反馈增益阵K。从而可 以利用lqgreg( )函数构造出相应的LQG控制器。 Lqgreg( )函数的基本格式为:
443

lqg=lqgreg(Kest,K) 其中,Kest是由Kalman( )函数设计得出的卡尔 曼状态估计器模型,K是由lqr( )或lqry设计得到 的优化反馈增益阵。 例 5-19 已知系统数学模型如下: C=[0 3.125] D=0 A ? ?? 1 ? 6.25? B ? ?2? 0 ? ?16 ?0? ? ? ? ?

G ? ? 3? H=0.1, Qn=1 ?0 ? ? ? Rn=0.01,Q=10,R=1,试设计LQG控制系统。
解:根据前面的知识结合Kalman( )与Lqgreg( ) 函数,不难编制出相应的程序,给出如下:
444

A=[-1 -6.25;16 0]; B=[2;0];C=[0 3.125];D=[0]; sys=ss(A,B,C,D); clsys0=feedback(sys,1); t=0:0.01:0.8; u1=25+t*0; N=length(t); u2=idinput(N,‘rgs’); G=[3;0]; H=[0.1]; BG=[B G]; DH=[D H]; LQGsys=ss(A,BG,C,DH); [kest,L,P]=kalman(LQGsys,1,0.01); K=lqry(sys,10,1)
445

LQGF=lqgreg(kest,K) feedin=[1];feedout=[1]; ressys=feedback(LQGsys,LQGF,feedin,feedout,1 ); clsys=feedback(sys,LQGF,1); u=[u1;u2‘]; figure(1) lsim(ressys,u,t) grid figure(3) step(clsys0,10)

程序运行结果得到: K= 10.2741 7.2395
446

a= x1_e x2_e b= x1_e x2_e c= u1 20.21485 7.30753 x1_e -10.27410 u1 0 x2_e -7.23945

x1_e x2_e -21.54820 -83.90030
16.00000 -22.83604

y1
d=

y1

447

图5-16给出了无LQG控制器的闭环系统的阶跃 响应。图5-17给出了具有LQG控制器的闭环系统阶 跃响应。当考虑系统存在白噪声干扰信号时,图519示出了具有噪声干扰的系统阶跃响应,图5-18即 为系统所加入的白噪声干扰信号。

对于离散系统,MATLAB同样提供了相应的函数 实现有关最优控制的设计,基本函数为:
[K,S,E] = DLQR(A,B,Q,R,N) [K,S,E] = DLQRY(A,B,C,D,Q,R) [Ac,Bc,Cc,Dc] = DREG(A,B,C,D,K,L) [KEST,L,P,M,Z] = KALMAN(SYS,Qn,Rn,Nn) RLQG = LQGREG(KEST,K)
448

1 0.8
Amplitude
Amplitude

12 10 8 6 4 2

0.6 0.4 0.2 0 0 5 10 (sec.)

0

0

0.2

0.4

0.6

(sec.)

图5-16 原闭环系统 阶跃响应

图 5-17 LQG 闭 环 系 统阶跃响应 449

4

10
2 0 -2 -4 0 50 100
Amplitude

5 0 0 0.2 0.4 0.6 (sec.)

图5-18 系统白噪声干 扰

图 5-19 具 有 白 噪 声 干扰的系统阶跃响应 450

倒立摆系统 仿真实例
? 倒立摆是进行控制理论研究的典型实验平台
? 倒立摆系统本身所具有的高阶次、不稳定、多变 量、非线性和强耦合等特性 ? 倒立摆系统是一个典型的非线性系统 ? 对倒立摆的研究具有重要的理论价值和实际意义
451

利用极点配置设计反馈控制器进行 倒立摆运动控制仿真
? 倒立摆系统原理图:
y
l sin?

x
?

m

l cos?

?

mg

x
P

u

M

452

倒立摆系统的数学模型
? 微分方程:
? u ? ml (sin ? )? 2 ? mg cos? sin? ?? ? x M ? m ? m cos2 ?
?2 ?? ? u cos? ? ( M ? m) g sin? ? ml (cos? sin? )? ? ml cos2 ? ? ( M ? m)l

? 线性化状态方程 :
?0 1 0 0 ? ? ( M ? m) g ? ? 0 0 0? Ml ? ??Z + ?0 0 0 1 ? ? ? ? mg ? ? 0 0 0 ? ? M ? ?

d ?Z ? dt

?0 ? ? ?1? ? ? ? Ml ? ? ? ?0 ? ? 1 ? ? ? ?M ?

?u
453

状态反馈控制系统的方块图
r
?

??

?

u

?

k1
?

? x ? Ax ? Bu

y

y ? Cx

k1
k2
k3

k4

454

仿真结果
? 初始角度为+5度、目标位置为1时

455

? 初始角度为-5度、目标位置为2时

456

结论
? 当初始状态及目标位置变化时,角度的仿真曲线 均比较快速地收敛到了零,位置的仿真曲线也迅速 地达到了期望值,满足了设计预期的要求。 ? 本次设计的状态反馈倒立摆控制系统具有一定的 自适应性。

457

习 题
5.1 已知系统?(A,B,C)的系数阵分别为:
?1 0 ? A?? 0 0? ? ?

?1? B??? ?1?

C ? ?2 ? 1?

设计状态观测器将其极点配置在 ?= ?=-10处。
5.2 已知系统?(A,B,C)的系数阵分别为:
?1 0 0 ? A ? ?0 2 1 ? ? ? ?0 0 2 ? ? ?

若系统状态不能测量到,试设计状态观测器, 并将其极点配置在:-3、-4、-5。
458

?1 ? B ? ?0 ? ? ? ?1 ? ? ?

C ? ?1 1 0?

5.3 已知系统?(A,B,C)的系数阵分别为:
?? 2 1 ? A?? 0 ? 1? ? ?
?0 ? B?? ? ?1 ?

C ? ?1 0?

试判断该系统能否进行极点配置。 若指定其极点位于 -3、-3处,试求状态反馈阵K。 5.4 已知系统?(A,B,C)的系数阵分别为:
?1 1 0 ? A ? ?0 1 0 ? ? ? ?0 0 2 ? ? ?
?0 0 ? B ? ?1 0 ? ? ? ?0 ? 2 ? ? ?

(1)试判断系统是否可控, (2)对系统进行极点配置使闭环极点配置在-2和 -1? 459 j2处

5.5 试将下列系统化为第二可控规范型。
(1)
?? 1 0 ? A?? 0 ? 2? ? ?
? 1 0 0? A ? ? 2 2 3? ? ? ?? 2 0 1? ? ?

?1? B??? ?1?
?1 ? B?? 2? ? ? ? ? 2? ? ?

(2)

5.6 试将下列系统化为第一可观规范型。 (1) (2)
? 1 0? A?? ? 2 4? ? ? ?? 2 2 ? A?? 0 ? 4? ? ?

C ? ?? 1 1?

C ? ?1 1?
460

5.7 已知系统同5.2题, 试设计具有观测器的状态反 馈闭环系统,系统的期望极点自行分析设置。

461


相关文章:
控制系统计算机仿真实验报告
控制系统计算机仿真实验报告_信息与通信_工程科技_专业资料。控制系统计算机仿真实验...0.825e(k ? 2) (7) 重新编写双重循环仿真程序如下(程序 3): clear G=...
中国矿业大学《控制系统计算机仿真》实验试题答案
中国矿业大学《控制系统计算机仿真》实验试题答案_工学_高等教育_教育专区。中国...0.85s 1 0.01s+1 0.15s+1 0.051s 70 0.00167s+1 0.21 0.15s+1...
控制系统计算机仿真-实验二
控制系统计算机仿真-实验二_工学_高等教育_教育专区。实验二 面向系统结构图的连续...(1) 步长 0.01: MAX_ERR= 0.0068,MEAN_ERR=0.0013 (2) 步长 0.05:...
《控制系统仿真与CAD》学习的感想
在这一过程中我学了很多东西, 最直接的就是将控制理论和 MATLAB 软件联系起来,用计算机来仿真在《自动 控制原理》中 《控制系统仿真与 CAD》学习的感想 学习了...
《控制系统仿真与CAD》教学大纲
上机内容 八、上机内容 1. 编写 M 函数 2. 系统的 SIMULINK 仿真 3. 控制系统的时域和频域分析 九、教材及主要参考资料 1. 教材:李国勇等, 《计算机仿真技术...
控制系统计算机仿真及辅助设计实验报告
控制系统计算机仿真及辅助设计实验报告_信息与通信_工程...05 0.60 65 0.60 71 0.00 06 0.54 88 0....(见图 8) ,对校正后系统仿真,得系统的单位阶跃...
《自动控制系统计算机仿真》习题参考答案
《自动控制系统计算机仿真》习题参考答案_农学_高等教育...? ? 0.1607 -0.0536 -0.0357 ? ? ?2.0...控制系统数字仿真与CAD第... 13页 免费©...
控制系统计算机仿真-实验五
控制系统计算机仿真-实验五_工学_高等教育_教育专区。控制系统计算机仿真-实验5 东南大学-自动化学院-控制系统计算机仿真-报告实验五 采样控制系统的数字仿真实验一、...
《计算机仿真技术与CAD》习题答案
计算机仿真技术与CAD》习题答案_工学_高等教育_教育专区。计算机仿真技术与 CAD——基于 MATLAB 的控制系统(第 3 版) 第0章答: 绪论 0-1 什么是仿真?它所...
计算机仿真教案
计算机仿真教案_工学_高等教育_教育专区。研究生课程计算机仿真讲稿 计算机仿真 教材: 《控制系统计算机仿真与 CAD—MATLAB 语言应用》 ,天津大学出版社 参考书: 1...
更多相关标签:
控制系统计算机仿真 | 计算机仿真投稿系统 | 系统建模与计算机仿真 | 计算机仿真系统 | 计算机仿真技术与cad | 控制系统仿真与cad | 控制系统cad | plc控制系统设计cad |