第1章 MATLAB简介
1.1 概述
MATLAB是MATrix LABoratory(矩阵实验室)的缩写,由美国The Math Works公司于1984年推出的一种科学与工程计算语言。主要特点:
一、简单易学。 二、代码短小高效。
三、功能丰富,可扩展性强。 四、强大的图形表达功能。 五、强有力的系统仿真功能。
1.2 桌面启动
启动MATLAB桌面主要采用以下两种方法:
一、在Windows桌面上,双击MATLAB的快捷方式图标。采用这种方式打开的MATLAB桌面以matlab*\\work为当前目录。
二、双击matlab*\\bin\\win32文件夹中的MATLAB.exe。采用这种方式打开的MATLAB桌面以matlab71为当前目录。
两者区别:当前目录不同。 注:*为MATLAB的软件版本号
1.3 通用操作界面简介
一、命令窗口
缺省情况下,位于桌面右侧,是用户与MATLAB进行人机对话的主要环境。在该窗口内,可输入各种由MATLAB运行的命令、函数、表达式,显示除图形外的所有运算结果。
二、命令历史窗口
缺省情况下,位于桌面左下方的前台,该窗口记录并显示每次开启MATLAB的时间及所有MATLAB运行过的命令、函数及表达式等,允许用户对它们进行选择复制、重运行及产生M文件。
三、当前目录浏览器
缺省情况下,位于MATLAB桌面左上方的前台。在该浏览器中,可以进行当前目录的设置,展示相应目录上的.m及.mdl等文件,复制、编辑和运行M文件以及装载MAT数据文件等。
四、工作空间浏览器
缺省情况下,位于MATLAB桌面左上方的后台,该窗口列出了MATLAB工作空间中所有数据的变量信息,包括变量名、大小、字节数等。在该窗口中,可以对变量进行观察、编辑、提取及保存。
五、数组编辑器
缺省情况下,不随操作界面的出现而启动,只有在工作空间浏览器中对变量进行操作时才启动。 六、开始按钮
缺省情况下,点击按钮会出现MATLAB的现场菜单。该菜单的菜单子项列出了已安装的各类MATLAB组件和桌面工具。
七、M文件编辑器/调试器
缺省情况下,不随操作界面的出现而启动,只有当进行“打开文件”等操作时才启动。 八、帮助导航/浏览器
缺省情况下,不随操作界面的出现而启动,只有在特意选择或设置的情况下,才以独立交互界面的形式出现。该浏览器详尽展示了由超文本写成的在线帮助。
1.4 运行方式
MATLAB提供了两种运行方式,即命令行方式和M文件方式。 一、命令行运行方式
可以通过在MATLAB命令窗口中输入命令行来实现计算或绘图功能。
第 1 页/共 63页
例:已知矩阵A=??5 6??1 2?,B=,完成矩阵求和运算A+B。 ????7 8??3 4?解:在MATLAB命令窗口输入下述内容:
>>A=[5 6;7 8]; >>B=[1 2;3 4]; >>C=A+B
按下回车键后,在MATLAB命令窗口显示运行结果如下: C= 6 8 1 12
二、M文件运行方式
命令行输入方式实际上也是MATLAB语言的一种程序编制方式,但这种方式只能编写简单的程序。若程序比较复杂,就应该把程序写成一个由多行命令组成的程序文件,即程序扩展名为.m的M文件,让MATLAB语言执行这个文件。
在MATLAB命令窗口中选择菜单“File|New|M-File”,即可打开一个缺省名为Untitled.m的M文件编辑/调试器窗口。把程序输入后,选择菜单“Debug|Run”即可运行。
M文件运行方式的优点是所编写的程序是以扩展名为.m的文件形式存储的,可调试,可重复运行,特别适合于求解复杂问题。
1.5 图形窗口
在MATLAB命令窗口中选择菜单“File|New|Figure”,或在命令窗口中输入“figure”或其他绘图命令,即可打开MATLAB的图形窗口。
1.6 帮助系统
MATLAB的帮助系统包括命令行帮助、联机帮助和演示帮助。 一、命令行帮助
命令行帮助是一种“纯文本”帮助方式。利用“help”命令就可以获得命令行帮助。 格式:help 函数或工具箱名称 二、联机帮助
由帮助导航/浏览器完成。打开方法有:
①在MATLAB命令窗口中运行命令“helpbrowser”或“helpdesk”。
②在MATLAB桌面上,用鼠标左键单击工具栏帮助图标,或选择菜单“Help|MATLAB Help”。 ③在MATLAB各独立出现的交互窗口中,选择菜单“Help|MATLAB Help”。 三、演示帮助
运行演示程序的方法有两种:
①在MATLAB命令窗口中运行命令“demos”。 ②在MATLAB命令窗口中选择菜单“Help|Demos”。 四、Web帮助
官方网站:http://www.matworks.com 中文论坛:http://www.ilovematlab.cn 五、PDF帮助
官方网站下载PDF帮助文档。
1.7 工具箱
MATLAB的工具箱分为辅助功能型工具箱和专业功能性工具箱。 一、控制系统工具箱
是专门对控制系统工程设计的函数和工具的集合。该工具箱主要采用M文件形式,提供了丰富的算法
第 2 页/共 63页
程序,主要用于反馈控制系统的建模、分析与设计。
控制系统工具箱的主要作用:
1、可以创建控制系统的个各种数学模型。
2、应用控制系统工具箱能够轻松地绘制控制系统的时间响应曲线、频率特性曲线及根轨迹图。 二、Simulink
Simulink是用来进行建模、分析和仿真各种动态系统的一种交互环境,它提供了采用鼠标拖放的方法建立系统框图模型的图形交互平台。其主要功能如下:
①交互建模;②交互仿真;③扩充和定制;④与MATLAB和工具箱集成。 三、其他解决控制领域问题的工具箱 ①系统辨识工具箱 ②模糊逻辑工具箱 ③鲁棒控制工具箱 ④模型预估控制工具箱
1.8 安装和内容选择
按照安装向导即可安装,另外也可以通过网络下载绿色免安装版直接使用。
第2章 MATLAB基本使用方法及常用功能介绍
2.1 应用基础
一、最简单的计算器使用方法
MATLAB的基本特性之一就是其演草纸式的数学运算功能,用户可以在命令窗口中进行各种数学演算。 例:求算术运算?9?(10?1)?19??2的结果。
2解:在MATLAB命令窗口中输入: >>(9*(10-1)+19)/2^2
按回车键,命令被执行,显示下述结果: ans=
25
说明:①在全部输入一个命令行内容后,必须按下回车键,该命令才会被执行。无需在命令行的末尾处执行,在一个命令行的任何一处都可执行。
②运算符号均为西文字符,不能在中文状态下输入。 ③“ans”是运算答案,是MATLAB的一个默认变量。
④如果不显示计算结果,可在命令行末尾添加分号,以分号结尾的命令行语句,尽管该命令已执行,但MATLAB不会把其运算结果显示在命令窗口中。
二、矩阵
1、矩阵的生成
在MATLAB中,矩阵的生成可以以矩阵的格式输入数据,也可以用“load”命令调用已存储的矩阵数据或矩阵变量,还可以应用MATLAB提供的函数生成特殊矩阵。
在MATLAB中输入矩阵要遵循以下基本规则:
①矩阵元素之间用空格或逗号分隔,矩阵行之间用分号隔离,整个矩阵放在方括号里,且标点符号一定要在英文状态下输入。
②不必事先对矩阵维数做任何说明,存储时将自动配置。 ③MATLAB区分字母的大小写。
第 3 页/共 63页
?1 1 1???例:以矩阵格式输入数据,自定义一个三阶帕斯卡矩阵A=1 2 3。 ????1 3 6??解:在MATLAB命令窗口中输入:
>>A=[1,1,1;1,2,3;1,3,6] %或者把逗号改为空格输入 运行结果为: A= 1 1 1 1 2 3 1 3 6
A(I,j)表示矩阵A中第i行第j列元素;A(i,:)表示矩阵A中第i行全部元素;A(:,j)表示矩阵A中第j列全部元素。
2、特殊矩阵的生成 (1)空矩阵
空矩阵用“[ ]”表示。空矩阵的大小为零,但变量名却保存在工作空间中。 (2)单位矩阵
单位矩阵使用函数eye( )实现,调用格式如下: eye(n) 生成n×n维单位矩阵 eye(n,m) 生成n×m维单位矩阵 (3)零矩阵
零矩阵用函数zeros( )实现,调用格式与函数eye( )完全相同。 (4)全部是1的矩阵
元素全部为1的矩阵可用函数ones( )实现,调用格式与函数eye( )完全相同。 (5)对角矩阵的生成
对角矩阵是指对角线上的元素为任意数,其他元素为零的矩阵。用函数diag( )实现。格式为: diag(V) diag(V,K)
说明:V为某个向量,K为向量偏离主对角线的列数。K=0,V在主对角线上;K>0,V在主对角线以上;K<0,V在主对角线以下。
例:对角矩阵生成演示。
解:在MATLAB命令窗口中输入: >> v=[1 2 3 4 5]; >> diag(v)
其运行结果为: ans =
1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 0 0 0 0 5 三、MATLAB的基本要素
MATLAB的基本要素包括变量、预定义变量、数值、字符串、运算符、标点符及复数等。 1、变量
第 4 页/共 63页
MATLAB会自动依据所赋予变量的值或对变量所进行的操作来识别变量的类型。如果赋值变量已存在,将使用新值代替旧值,并以新值类型代替旧值类型。
MATLAB变量的命名遵循以下规则: (1)变量均先定义、后使用。 (2)变量名以英文字母开头。
(3)变量名可以由字母、数字和下划线混合组成。
(4)对于6.5以上版本,变量名最多可包含63个字符。 (5)变量名中不得包含空格和标点,但可以包含下划线。 (6)MATLAB区分变量大小写。 2、预定义变量
在MATLAB中存在一些固定变量(也称为常量),这就是MATLAB默认的预定义变量,也称为默认变量,每当MATLAB启动时,这些变量就被产生。
MATLAB的预定义变量 名称 ans beep bitmax eps i或j Inf或inf NaN或nan 3、数值
可以使用十进制计数法,也可以使用科学计数法,数值的有效范围为10?308变量含义 计算结果的缺省变量名 使计算机发出“嘟嘟”声 最大正整数,9.0072×10 计算机中的最小数ε,ε=2虚数单位,定义为?1 无穷大,如1 / 0 不定值,如0 / 0,∞/∞,0*∞ ?5215名称 nargin nargout pi realmin realmax varagin varagout 变量含义 函数输入变量个数 函数输出变量个数 圆周率π 最小正实数,2?1022 最大正实数,(2??)21023 可变的函数输入变量个数 可变的函数输出变量个数 ~10308。
4、字符串
创建字符串的方法:先将待建的字符串放在一个“单引号对”中,再按回车键,且该单引号对必须在英文状态下输入,但字符串内容可以为中文。
5、运算符
MATLAB的运算符包括算术运算符、关系运算符和逻辑运算符。
MATLAB的算术运算符 操作符 + - * ^ \\ 操作符 = = ~= > 6、标点符
第 5 页/共 63页
功能 算术加 算术减 算术乘 算术乘方 算术左除 操作符 / .* .^ .\\ ./ 功能 大于等于 小于 小于等于 功能 算术右除 点乘 点乘方 点左除 点右除 操作符 & | 功能 与 或 非 MATLAB的关系运算符及逻辑运算符 功能 等于 不等于 大于 操作符 >= < <= ~
所有标点符号均在西文状态下输入。
MATLAB的标点符 标点符 : ; , () [ ] { } 功能 冒号,具有多种应用功能 分号,区分行及取消运行显示灯 括号,指定运算优先级 方括号,矩阵定义的标志灯 花括号,用于构成元胞数组等 标点符 . ? ! = ‘ 续行符 百分号,注释标记 惊叹号,调用操作系统运算 等号,赋值标记 单引号,字符串的标识符,必须成对使用 功能 小数点及域访问符等 逗号,区分列及函数参数分隔符等 % (1)冒号
在MATLAB中,冒号不仅可以定义行向量,还可以截取指定矩阵中的部分元素。 例:用冒号定义增量为1的行向量。 解:在MATLAB命令窗口中输入: >> a=2:8
运行结果为: a =
2 3 4 5 6 7 8 例:用冒号定义增量为给定值的行向量。 解:在MATLAB命令窗口中输入: >> a=0:10:80 运行结果为: a =
0 10 20 30 40 50 60 70 80
例:用冒号截取指定矩阵中的部分元素。 解:在MATLAB命令窗口中输入: >> A=[1 2 3;4 5 6;7 8 9];
>> B=A(1:2,:) %取出矩阵A的第1行和第2行 运行结果为: B =
1 2 3 4 5 6 (2)分号
分号在矩阵中用来分隔行,如果不希望某些运算结果显示在屏幕中,还可以用分号作为该行结束的标志。
7、复数
复数的生成可以利用下面语句:
z=a+bi 或 z=r*exp(θ*i),其中r是复数的模,θ是复数幅角的弧度数。
例:已知复数z1?3?i4,z2?1?i2,z3?2e6,计算z?解:在MATLAB命令窗口中输入: >> z1=3+4i; >> z2=1+2i;
>> z3=2*exp((pi/6)*i); >> z=z1*z2/z3
i?z1z2。 z3第 6 页/共 63页
运行结果为: z =
0.3349 + 5.5801i
2.2 基本操作
一、命令窗口
1、命令窗口显示及设置
个性设置方法:选择菜单“File|Preference”,打开参数设置对话框。 2、命令窗口的常用控制命令 命令 功能 cd 设置当前工作目录 clf 清除图形窗口 clc 清除命令窗口中显示的内容 命令 exit quit more 功能 关闭/退出MATLAB 关闭/退出MATLAB 使其后的显示内容分页进行 clear 清除工作空间中保存的变量 type 显示指定M文件的内容 dir 列出指定目录下的文件和子目录清单 which 指出其后文件所在的目录 edit 打开M文件编辑器 3、命令窗口中命令行的编辑 键名 ↑ ↓ ← → PageUp 功能 键名 功能 使光标移到当前行的首端 使光标移到当前行的尾端 删去光标右边的字符 前寻式调回已输入过的命令行 Home 后寻式调回已输入过的命令行 End 在当前行中左移光标 在当前行中右移光标 Delete Backspace 删去光标左边的字符 前寻式翻阅当前窗口中的内容 Esc 清除当前行的全部内容 PageDown 后寻式翻阅当前窗口中的内容 二、命令历史窗口
命令历史窗口的主要应用功能及操作方法 应用功能 复制单行或多行命令 操作方法 选中单行或多行命令;单击鼠标右键打开现场菜单;选择菜单“Copy”;把选中的单行或多行命令粘贴到包括命令窗口在内的任何地方 选中单行命令;单击鼠标右键打开现场菜单;选择菜单“Evaluate Selection”;在命令窗口中运行 选中多行命令;单击鼠标右键打开现场菜单;选择菜单“Evaluate Selection”;在命令窗口中运行 选中多行命令;单击鼠标右键打开现场菜单;选择菜单“Create M-File”,打开书写这些命令的M文件编辑/调试器;进行相应操作,即建立所需的M文件 简捷操作方法 选中变量之后,按“Ctrl+C”键 运行单行命令 用鼠标左键双击单行命令 — 运行多行命令 将多行命令写成M文件 — 三、当前目录浏览器
文件详细列表区的主要应用功能及操作方法 应用功能 运行M文件 操作方法 选中文件;单击鼠标右键;选择菜单“Run” 简捷操作方法 — 第 7 页/共 63页
编辑M文件 把MAT文件的全部数据输入工作空间 把MAT文件的部分数据输入内存 选中文件,单击鼠标右键;选择菜单“Open” 双击M文件 选中数据文件,单击鼠标右键;选择菜单“Open” 双击MAT文件 选中数据文件,单击鼠标右键;选择菜单“Import — Data”,打开数据预览选择对话框“Import Wizard”;选中待装载数据变量名,单击“Finish” 四、工作空间浏览器 工作空间是指运行MATLAB的程序或命令时生成的所有变量与MATLAB提供的常量构成的空间,也称为内存空间。
工作空间浏览器的主要应用功能及操作方法
应用功能 变量的字符显示 变量的图形显示 部分变量保存为MAT文件 重命名变量名 变量复制 操作方法 选中变量;右键选择菜单“Open Selection” 选中变量;右键选择菜单“Plot all columns” — — — 选中若干变量,右键选择菜单“Save as?” 简捷操作方法 用鼠标左键双击变量 全部内存变量保存为MAT文件 右键选择菜单“Save as?” 选中欲重命名的变量;右键选择“Rename” — Ctrl+C 选中若干变量,右键选择“Copy” 五、数组编辑器
是工作空间浏览器的一个组件,用于生成数组、观察数组内容以及编辑其值。打开的三种方法: 1、选中工作空间浏览器中的任意一维或二维数组,再双击该数组。 2、单击工作空间浏览器的工具栏图标。 3、选择菜单“Open Selection”。
通常在命令窗口中输入较大规模数组时,先在命令窗口中向一个新变量赋“空”矩阵,然后打开数组编辑器逐格填写数组元素值。
六、数据文件的存取 1、数据文件的保存
save FileName 将全部变量保存为当前目录下的FileName.mat文件 save FileName v1 v2 将变量v1,v2保存为FileName.mat文件
save FileName v1 v2 –append 将变量v1,v2添加到已有的FileName.mat文件中 save FileName v1 v2 -ascii 将变量v1,v2保存为FileName 8位ASCII文件 save FileName v1 v2 –ascii -double 将变量v1,v2保存为FileName 16位ASCII文件 2、数据文件的调入
load FileName 将FileName.mat文件中的全部变量装入工作空间 load FileName v1 v2 将FileName.mat文件中的v1,v2变量装入工作空间 load FileName v1 v2 -ascii 将FileName ASCII文件中的v1,v2变量装入工作空间
2.3 数值运算
一、向量及其运算 1、向量的生成
(1)在命令窗口中直接生成向量 例:命令窗口直接生成向量演示 解:在MATLAB命令窗口中输入: >> X1=[1 2 3 4 5] 运行结果为: X1 =
1 2 3 4 5
第 8 页/共 63页
>> X2=[1;2;3;4;5]' %求列向量的转置,用右单引号,而不是一般线性代数中的上标“T” 运行结果为: X2 =
1 2 3 4 5
(2)等差元素向量的生成
①冒号生成法。基本格式为V=a:n:b。V为生成的向量,a为向量V的第一个元素,b为向量V的最后一个元素;n为步长,缺省设置为1,且n=1时可忽略。
②使用linspace( )函数。格式:X=linspace(a,b,n)
生成元素在[a,b]之间的线性等分行向量,向量元素个数为n,n的缺省值为100。 例:等差元素向量生成演示
解:在MATLAB命令窗口中输入: >> X1=1:2:9
运行结果为: X1 =
1 3 5 7 9 >> X2=linspace(10,-2,5)
运行结果为: X2 =
10 7 4 1 -2 2、向量的基本运算
(1)向量与常数的四则运算
指向量中的每个元素与常数进行的加减乘除等运算,符号分别为+-*/。当进行除法运算时,向量只能作为被除数。
(2)向量与向量之间的加减运算
指向量中的每个元素与另一个向量中相对应元素加减运算,运算符号为+-。 (3)向量的点积和叉积运算
向量的点积等于其中一个向量的模与另一个向量的模在这个向量方向上投影的乘积。向量叉积是指过两个相交向量的交点并与两向量所在平面垂直的向量,且向量维数只能为3。在MATLAB中使用函数dot()与cross()分别计算向量的点积与叉积。
例:向量的点积与叉积运算演示。 解:在MATLAB命令窗口中输入: >> A=[10 20 30]; >> B=[40 50 60];
>> C=dot(A,B) %计算点积 运行结果为: C =
3200
>> D=cross(A,B) %计算叉积 运行结果为: D =
-300 600 -300 二、数组及其运算 1、数组的概念
数组是一组实数或复数排成的长方阵列。单维数组通常指单行或单列的矩阵,即行向量或列向量。多维数组可以认为是矩阵在维数上的扩张。
第 9 页/共 63页
2、数组的基本数值运算
(1)数组与常数的四则运算
单维数组与常数的运算与向量与数的运算完全相同。 例:数组与常数的四则运算演示。 解:在MATLAB命令窗口中输入: >> A=[1 2 3;2 3 4;3 4 5]; >> B=[1 2 3;4 5 6;7 8 9]; >> s=5;
>> C=s*A-B/s+10
运行结果为: C =
14.8000 19.6000 24.4000 19.2000 24.0000 28.8000 23.6000 28.4000 33.2000
(2)数组间的四则运算
按元素与元素的方式进行。加减法运算与矩阵的加减运算完全相同;数组间的相乘、相除运算符号为“.*”、“./”或“ .\\”。
例:数组相乘运算演示。
解:在MATLAB命令窗口中输入: >> A=[1 3 5;2 4 6;3 5 7]; >> B=[2 4 6;1 3 5;3 5 7]; >> C=A.*B
运算结果为: C =
2 12 30 2 12 30 9 25 49
例:数组相除运算演示。
解:在MATLAB命令窗口中输入: >> A=[1 3 5;2 4 6;3 5 7]; >> B=[2 4 6;1 3 5;3 5 7];
>> C=A.\\B %点左除 运行结果为: C =
2.0000 1.3333 1.2000 0.5000 0.7500 0.8333 1.0000 1.0000 1.0000 >> D=A./B %点右除 运行结果为: D =
0.5000 0.7500 0.8333 2.0000 1.3333 1.2000 1.0000 1.0000 1.0000 (3)数组的乘方运算
数组的乘方运算(幂运算)符号为“.^”,按元素对元素的幂运算进行。与矩阵的幂运算完全不同。
第 10 页/共 63页
例:数组的乘方运算演示。
解:在MATLAB命令窗口中输入: >> A=[1 3 5;2 4 6;3 5 7]; >> E=A.^2
运行结果为: E =
1 9 25 4 16 36 9 25 49
说明:①数组乘、除及乘方运算符前的小黑点绝不能遗漏,否则将不按数组运算规则运算;②在执行数组与数组运算时,参与运算的数组必须维数相同,运算结果所得数组也总与原数组维数相同。
3、元胞数组(Cell Array)
元胞数组是MATLAB中一种特殊的数组,它的基本元素是元胞(Cell),每个元胞本身在数组中是平等的,只能以下标区分。元胞可以存放任何类型、任何大小的数组,包括任意维数值数组、字符串数组以及符号对象等,并且同一个元胞数组中各元胞中的内容可以不同。
在MATLAB中,元胞数组必须用花括号“{}”。生成元胞数组的方法有: (1)使用函数cell( )生成元胞数组
格式:c=cell(n) 生成n×n维空元胞数组 c=cell(m,n) 生成m×n维空元胞数组
c=cell(m,n,p,?) 生成m×n×p×?维空元胞数组
c=cell(size(A)) 生成与A维数组相同的空元胞数组,A为数值数组或元胞数组 例:使用函数cell( )生成一个2×2维元胞数组。 解:在MATLAB命令窗口中输入:
>> A=ones(3,4) %生成3×4维全部元素为1的数值矩阵 A =
1 1 1 1 1 1 1 1 1 1 1 1
>> C=cell({A,[1,2];'cell',[1;2]}) %使用cell( )函数生成元胞数组C C =
[3x4 double] [1x2 double] 'cell' [2x1 double]
>> C(:,1) %显示元胞数组C的第1列内容 ans =
[3x4 double] 'cell'
(2)使用花括号{ }生成元胞数组
例:直接用花括号生成上例中的元胞数组。 解:在MATLAB命令窗口中输入: >> A=ones(3,4);
>> C={A,[1,2];'cell',[1;2]} C =
[3x4 double] [1x2 double] 'cell' [2x1 double]
建议优先采用第二种方法生成元胞数组。
第 11 页/共 63页
三、基本数学函数运算
MATLAB的基本数学函数 三角函数和双曲函数 函数名 功能 函数名 功能 函数名 功能 指数与对数函数 函数名 功能 sin cos tan cot sec csc asin acos atan acot asec acsc sinh 正弦 余弦 正切 余切 正割 余割 反正弦 反余弦 反正切 反余切 反正割 反余割 双曲正弦 cosh tanh coth sech csch asinh acosh atanh acoth asech acsch atan2 双曲余弦 双曲正切 双曲余切 双曲正割 双曲余割 反双曲正弦 反双曲余弦 反双曲正切 反双曲余切 反双曲正割 反双曲余割 exp log log10 abs angle conj cart2sph cart2pol 指数 自然对数 常用对数 绝对值(模) 复数的相角 复数的共轭 直角坐标变为球坐标 直角坐标变为柱(或极坐标) log2 pow2 sqrt 复数函数 imag isreal 坐标变换函数 pol2cart sph2cart 以2为底的对数 以2为底的指数 平方根 复数的虚部 复数的实部 柱坐标变为直角坐标 球坐标变为直角坐标 四象限反正切 sign 符号函数 误差函数 误差补函数 其他特殊函数 expint gamma gammaln isprime 指数积分函数 Γ函数 Γ函数的对数 质数为真函数 erf erfc erfcx erfinv gammainc 不完全Γ函数 刻度误差补函数 逆误差函数 例:求tan45°的函数值。
解:在MATLAB命令窗口中输入:
>> tan(pi/4) %三角函数运算的角度单位是弧度,要先由度转化为弧度 ans =
1.0000
例:求复数z=5+i5的模和相角。 解:在MATLAB命令窗口中输入: >> z=5+5i; >> Am=abs(z) Am =
7.0711 >> Fm=angle(z) Fm =
0.7854
四、矩阵的函数运算
矩阵的基本数值运算与数组的基本数值运算基本相同,不同的是:矩阵之间进行乘、除及乘方运算的运算符没有小黑点。
实现矩阵特有运算的函数
函数名 功能 sqrtm 矩阵开方运算 expm logm cond 矩阵指数运算 矩阵对数运算 求矩阵的条件数 函数名 gsvd inv 功能 广义奇异值 矩阵求逆 norm/normest 求矩阵和向量的范数 null 矩阵的零空间 pinv 伪逆矩阵 condest 求矩阵1的范数条件估计 第 12 页/共 63页
condeig 求与矩阵特征值有关的条件数 poly det 求矩阵的行列式 polyvalm 求矩阵的特征值多项式 求矩阵特征值多项式的值 eig/eigs 求矩阵的特征值和特征向量 rank 求矩阵的秩 funm 矩阵的任意函数 trace 求矩阵的迹 1、矩阵的转置、逆运算与行列式运算 例:求矩阵的转置演示。
解:在MATLAB命令窗口中输入: >> A=[1 2 3;4 5 6;7 8 9]; >> C=A' C =
1 4 7 2 5 8 3 6 9 例:矩阵求逆演示。
解:在MATLAB命令窗口中输入: >> A=[1 2 0;2 5 -1;4 10 -1]; >> B=inv(A) B =
5 2 -2 -2 -1 1 0 -2 1
例:求矩阵的行列式演示。
解:在MATLAB命令窗口中输入: >> A=[1 2 0;2 5 -1;4 10 -1]; >> B=det(A) B = 1
2、矩阵的特征值运算 使用函数eig( )。格式:
d=eig(A) %d为矩阵A的特征值向量 D=eig(A)
[V,D]=eig(A) %V、D分别为矩阵A的特征向量矩阵与特征值矩阵
?例:求矩阵A=?1 2 0??2 5 ?1?的特征值向量、特征向量矩阵和特征值矩阵。 ??1??4 10 ??解:在MATLAB命令窗口中输入:
>> A=[1 2 0;2 5 -1;4 10 -1];
>> d=eig(A) %求矩阵A的特征值向量 d =
3.7321 0.2679 1.0000
>> [B,C]=eig(A) %求矩阵A的特征向量矩阵和特征值矩阵 B =
第 13 页/共 63页
-0.2440 -0.9107 0.4472 -0.3333 0.3333 0.0000 -0.9107 -0.2440 0.8944 C =
3.7321 0 0 0 0.2679 0 0 0 1.0000 3、矩阵的秩运算 例:矩阵求秩演示。
解:在MATLAB命令窗口中输入: >> A=[1 2 0;2 5 -1;4 10 -1]; >> B=rank(A) B = 3
五、多项式及其运算
1、多项式的表达及其构造
代数运算中,多项式一般可表示为如下形式:
P(x)?a0xn?a1xn?1?a2xn?2???an?1x?an
将上式中系数按降幂次序可存放在如下的行向量中:
P?[a0,a1,a2,?,an?1,an]
例:用MATLAB构造多项式P(x)?2x5?5x4?4x2?x?4。
解:在MATLAB命令窗口中输入:
>> P=[2 5 0 4 1 4]; %无论多项式的系数是否为零,都必须写完整 >> poly2sym(P) ans =
2*x^5+5*x^4+4*x^2+x+4 2、多项式运算函数 调用格式 p=conv(p1,p2) p=poly(AR) dp=polyder(p) dp=polyder(p1,p2) [num,den]= polyder(p1,p2) p=polyfit(x,y,n) pA=polyval(p,S) pM=polyvalm(p,S) [r,p,k]= residue(num,den) r=roots(p) 多项式求导数 功能 说明 多项式卷积(乘法) p是多项式p1和p2的乘积多项式 由根求多项式 求方阵AR的特征多项式p,或求向量AR指定根所对应的多项式 求多项式p的导数多项式dp 求多项式p1,p2乘积的导数多项式dp 对有理分式(p1/p2)求导数所得的有理分式为(num/den) 求x,y向量给定数据的n阶拟合多项式p 按数组运算规则计算多项式值,p为多项式,S为矩阵 按矩阵运算规则计算多项式值,p为多项式,S为矩阵 num、den分别是分子、分母多项式系数向量,r、p、k分别是 留数、极点、直项 r是多项式p的根向量 [q,r]=deconv(p1,p2) 多项式解卷(除法) q是p1被p2除的商多项式,r是余多项式 多项式曲线拟合 多项式求值 矩阵多项式求值 分式多项式的 部分分式展开 多项式求根 3、多项式运算举例
第 14 页/共 63页
(1)代数方程求根
例:求方程x?2x?24x?48x?25x?50?0的根。 解:在MATLAB命令窗口中输入: >> P=[1 2 24 48 -25 -50]; >> r=roots(P) r =
0.0000 + 5.0000i 0.0000 - 5.0000i 1.0000 -2.0000 -1.0000
(2)用多项式的根构造多项式
例:用多项式的根构造多项式P(x)?x5?2.5x4?2x2?0.5x?2。 解:在MATLAB命令窗口中输入: >> P=[1 2.5 0 2 0.5 2]; >> r=roots(P) r =
-2.7709 0.5611 + 0.7840i 0.5611 - 0.7840i -0.4257 + 0.7716i -0.4257 - 0.7716i >> poly(r) ans =
1.0000 2.5000 0.0000 2.0000 0.5000 2.0000
说明:当用多项式的根生成多项式时,如果某些根有虚部,则可以通过使用函数real( )抽取实部来消除。
(3)求矩阵的特征多项式
5432?1.2 3 ?0.9???例:求矩阵A=5 1.75 6的特征多项式。 ????9 0 1??解:在MATLAB命令窗口中输入:
>> A=[1.2 3 -0.9;5 1.75 6;9 0 1]; >> poly(A) ans =
1.0000 -3.9500 -1.8500 -163.2750
即矩阵A的特征多项式为f(s)?s?3.95s?1.85s?163.275 (4)多项式卷积(乘法)与多项式解卷(除法)
例:已知多项式P(x)?x?2x?3x?4和q(x)?10x?20x?30,求两个多项式的卷积p(x)*q(x),并用多项式解卷验证。
32232第 15 页/共 63页
解:①求两个多项式的卷积。 在MATLAB命令窗口中输入: >> p=[1 2 3 4]; >> q=[10 20 30]; >> c=conv(p,q) c =
10 40 100 160 170 120 ②用多项式解卷验证 >> [s,r]=deconv(c,p) s =
10 20 30 r =
0 0 0 0 0 0
说明:①s是向量c除以向量p所得的结果,r为余数。 ②函数conv( )与deconv( )互为逆运算。
③建立控制系统的传递函数模型时,会经常使用函数conv( )。 (5)分式多项式的部分分式展开
应用时域分析法分析控制系统动态性能时,常常需要求出系统在典型输入信号作用下的时间响应c(t),为此,必须首先求出c(t)的象函数C(s),并将其展开成部分分式。利用函数residue( )。
例:已知控制系统的输出象函数C(s)?解:>> num=[2 8]; >> den=[1 5 6 0];
>> [z,p,k]=residue(num,den) z =
0.6667 -2.0000 1.3333 p =
-3.0000 -2.0000 0 k = []
(6)多项式曲线拟合
进行系统设计与仿真时,常常需要用曲线拟合方法。曲线拟合就是要寻找一条光滑曲线,使其在某种准则下能最佳地拟合已知数据。使用函数polyfit( )对已知数据进行曲线拟合。拟合方法采用最小二乘法(即最小误差平方和准则)。
例:曲线拟合实例。
解:>> x=0:0.1:1; %生成用行向量表示的自变量数据 >> y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; >> p=polyfit(x,y,2) %计算二阶拟合多项式系数 p =
-9.8108 20.1293 -0.0317
2(s?4),将其展开为部分分式。 2s(s?5s?6)第 16 页/共 63页
2.4 符号运算
一、符号对象的创建和使用
进行符号运算时,首先要创建基本的符号对象,它可以是常数、变量和表达式。然后利用这些基本符号对象构成新的表达式,进而完成所需的符号运算。
符号对象的创建使用函数sym( )和syms( )来完成,他们的调用格式如下:
S=sym(A) 将数值A转换成符号对象S,A是数字(值)或数值矩阵或数值表达式 S=sym(‘x’) 将字符串x转换成符号对象S
S=sym(A,flag) 将数值A转换成flag格式的符号对象
syms arg1 arg2 ? arg1=sym(‘arg1’),arg2=sym(‘arg2’),?的简洁形式 例:创建符号变量和符号表达式演示。
解:>> y=sym('x') %定义变量y,它代表字符x y = x
>> f=sym('x^3+x^2+4*x+4') %定义变量f,它代表符号表达式x3?x2?4x?4 f =
x^3+x^2+4*x+4
例:字符表达式转换为符号变量演示。
解:>> y=sym('2*sin(x)*cos(x)') %将字符表达式转换为符号变量 y =
2*sin(x)*cos(x)
>> y=simple(y) %将已有的y符号表达式化成最简形式 y = sin(2*x)
二、符号运算中的运算符号和基本函数 例:基本运算符号与基本函数应用演示。 解:>> syms x; %定义符号变量x >> f1=x^3+x^2+4*x+4; %生成多项式f1 >> f2=x^2+4*x+10; %生成多项式f2 >> f=f1+f2 %求f1与f2之和 f =
x^3+2*x^2+8*x+14
>> y=sqrt(x^5) %求函数的平方根 y =
(x^5)^(1/2)
>> z=log10(x) %求以10为底的对数 z =
log(x)/log(10)
例:矩阵代数运算演示。求矩阵A=??a11 a12??a?的行列式值、逆和特征值。21 a
22?解:>> syms a11 a12 a21 a22; %定义符号变量a11,a12,a21,a22
>> A=[a11,a12;a21,a22] %生成矩阵A A =
[ a11, a12]
第 17 页/共 63页
[ a21, a22]
>> DA=det(A) %求矩阵A的行列式 DA =
a11*a22-a12*a21
>> IA=inv(A) %求矩阵A的逆矩阵 IA =
[ a22/(a11*a22-a12*a21), -a12/(a11*a22-a12*a21)] [ -a21/(a11*a22-a12*a21), a11/(a11*a22-a12*a21)]
>> EA=eig(A) %求矩阵A的特征值 EA =
1/2*a11+1/2*a22+1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2) 1/2*a11+1/2*a22-1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2) 三、符号表达式的操作 调用格式 collect(E,v) expand(E) factor(E) horner(E) 功能 同类项合并 表达式展开 因式分解 嵌套分局 说明 将符号表达式E中v的同幂项系数合并 对E进行多项式、三角函数、指数函数及对数函数等展开 对E(或正整数)进行因式(或因子)分解 将E分解成嵌套形式 将E通分,返回E通分后的分子N与分母D 运用多种恒等式变换对E进行综合简化 运用包括simplify在内的各种简化算法,把E转换成最简短形式 [N,D]=numden(E) 表达式通分 simplify(E) 表达式化简 simple(E) subs(E,old,new) 符号变量替换 将E中的符号变量old替换为new,new可以是符号变量、符号 常量、双精度数值与数值数组(或矩阵)等 1、符号(包括符号变量、符号常量及数值数组)替换 例:已知数学表达式y?axn?bt?c,试对其进行以下的符号替换: (1)a=sint,b=lnz,c=de的符号变量替换; (2)n=3,c=π的符号常量替换; (3)c=1:2:5的数值数组替换; (4)c=?2t
?1 2??的数值矩阵替换。 3 4??解:>> syms a b c d e t n x y z;
>> y=a*x^n+b*t+c;
>> y1=subs(y,[a,b,c],[sin(t) log(z) d*exp(2*t)]) %符号变量替换 y1 =
sin(t)*x^n+log(z)*t+d*exp(2*t)
>> y2=subs(y,[n,c],[3 pi]) %符号常量替换 y2 =
a*x^3+b*t+pi
>> y3=subs(y,c,1:2:5) %符号数组替换 y3 =
[ a*x^n+b*t+1, a*x^n+b*t+3, a*x^n+b*t+5] >> y4=subs(y,c,[1 2;3 4]) %数值矩阵替换
第 18 页/共 63页
y4 =
[ a*x^n+b*t+1, a*x^n+b*t+2] [ a*x^n+b*t+3, a*x^n+b*t+4] 2、同类项合并
例:已知数学表达式y?(x2?xe?t?1)(x?e?t),试对其同类项进行合并。 解:>> syms x t;
>> y=sym('(x^2+x*exp(-t)+1)*(x+exp(-t))');
>> y1=collect(y) %默认合并x同幂项系数 y1 =
x^3+2*exp(-t)*x^2+(1+exp(-t)^2)*x+exp(-t)
>> y2=collect(y,'exp(-t)') %合并exp(-t)同幂项系数 y2 =
x*exp(-t)^2+(2*x^2+1)*exp(-t)+(x^2+1)*x 3、因式分解
例:已知数学表达式y(x)?x4?5x3?5x2?5x?6,试对其进行因式分解。解:>> syms x;
>> y=x^4-5*x^3+5*x^2+5*x-6; >> y1=factor(y) y1 =
(x-1)*(x-2)*(x-3)*(x+1)
4、表达式展开
例:已知数学表达式y(x)=cos(3arccosx),试将其展开。 解:>> syms x;
>> y=cos(3*acos(x)); >> y1=expand(y) y1 =
4*x^3-3*x
5、表达式简化
例:已知数学表达式y(x)?2cos2x?sin2x,试对其进行简化。 解:>> syms x;
>> y=2*cos(x)^2-sin(x)^2; >> y1=simplify(y) y1 =
3*cos(x)^2-1 6、表达式通分
例:已知数学表达式y(x)?x?3x?x(x?1)?1x2(x?2),试对其进行通分。
解:>> syms x;
>> y=((x+3)/(x*(x+1)))+((x-1)/(x^2*(x+2))); >> [n,d]=numden(y) n =
第 19 页/共 63页
x^3+6*x^2+6*x-1 d =
x^2*(x+1)*(x+2) 四、符号积分变换
1、拉氏变换及其反变换
设f(t)是一个以时间t为自变量的函数,它的定义域是t>0,并设|f(t)|≤ke(a为正数),则对所有实部大于a的复数,积分
at??0f(t)edt绝对收敛。f(t)的拉氏变换F(s)定义为F(s)=
?st??0f(t)e?stdt。
F(s)称为f(t)的象函数,而f(t)称为F(s)的原函数。s=ζ+jω为复变量,称其为拉氏变换算子。从f(t)求取F(s)的过程称为拉氏正变换。
在MATLAB中,分别使用函数laplace( )与ilaplace( )进行拉氏变换与反变换。 (1)函数laplace( )
功能:求取函数的拉氏变换。格式:F=laplace(f)。 例:求单位阶跃函数f(t)=1(t)的拉氏变换。 解:>> f=sym(1); >> F=laplace(f) F = 1/s
例:求函数t,e及sinωt的拉氏变换。 解:>> syms a t w; >> F1=laplace(t^5) F1 = 120/s^6
>> F2=laplace(exp(a*t)) F2 = 1/(s-a)
>> F3=laplace(sin(w*t)) F3 =
w/(s^2+w^2)
(2)函数ilaplace( )
功能:求取函数的拉氏反变换。格式:f=ilaplace(F) 例:求象函数F(s)=1?5atab?的拉氏反变换。 2s?1(s?2)解:>> syms a b s;
>> F=1+a/(s+1)+b/(s+2)^2; >> f=ilaplace(F) f =
dirac(t)+a*exp(-t)+b*t*exp(-2*t)
说明:dirac(t)表示单位脉冲函数δ(t)。即f(t)=δ(t)+ae+bte2、Z变换及Z反变换 (1)函数ztrans( )
?1
?2t。
第 20 页/共 63页
功能:求取函数的Z变换。格式:F=ztrans(f)。 例:求单位阶跃函数f(t)=1(t)的Z变换。 解:>> n=sym(1); >> F=ztrans(n) F = z/(z-1)
(2)函数iztrans( )
功能:求取函数的Z反变换。格式:f=iztrans(F)。
z2例:求Z变换函数F(z)?的Z反变换。
(z?1)(z?0.5)解:>> syms z n;
>> F=z^2/((z-1)*(z-0.5)); >> f=iztrans(F) f =
2-(1/2)^n
2.5 图形表达功能
一、二维绘图
1、函数plot( )的调用格式 (1)plot(x,’s’)调用格式
此格式是缺省自变量绘图格式,功能如下:
①如果x为实向量,则绘制出以该向量元素下标为横坐标,元素值为纵坐标的曲线。 ②如果x为复数向量,则绘制出以该向量元素实部为横坐标,虚部为纵坐标的曲线。
③如果x为实矩阵,则按列分别绘制出以每列元素下标为横坐标,每列元素值为纵坐标的多条曲线,其曲线数等于x的列数。
④如果x为复数矩阵,则按列分别绘制出以每列元素实部为横坐标,虚部为纵坐标的多条曲线,其曲线数等于x的列数。
⑤s为选项(开关量)字符串,用于设置曲线颜色、线型、数据点型等。 (2)plot(x,y,’s’)调用格式
此格式是基本绘图格式,功能如下:
①如果x,y是相同维数向量,则绘制出以x为横坐标,以y为纵坐标的曲线。
②如果x是向量,y是矩阵,且y的行或列的维数与x的维数相同,则绘制出以x为横坐标的多条不同颜色的曲线,曲线数等于x的维数。
③如果x是矩阵,y是向量,情况与②类似,以y为横坐标。
④如果x,y是相同维数矩阵,则以x对应列元素为横坐标,以y对应列元素为纵坐标分别绘制曲线,曲线数等于矩阵的列数。
(3)plot(x1,y1,’s1’,x2,y2,’s2’,?)调用格式
此曲线是多条曲线绘图格式。其中,x1,y1,x2,y2,?为数组对,每一对(x,y)数组可绘制出一条曲线,且每一数组对的长度可以不同。
例:已知函数y(x)=sinxcosx,且x?[0,π],绘制y(x)曲线。 解:>> x=0:0.01*pi:pi; >> y=sin(x).*cos(x); >> plot(x,y)
第 21 页/共 63页
例:绘制多条不同色彩曲线的演示。 解:>> x=(0:pi/50:2*pi)'; >> k=0.4:0.1:1; >> Y=cos(x)*k; >> plot(Y)
2、多次重叠绘制曲线
为了在一张图中绘制多条曲线,即多次重叠绘制曲线,就必须使用“hold”命令。格式: hold on 使当前曲线与坐标轴具备不被刷新的功能
hold off 使当前曲线与坐标轴不再具备不被刷新的功能 hold 当前图形是否具备被刷新功能的双向切换开关
说明:当前曲线与坐标轴不被刷新是指,再次使用plot函数时,在此之前所绘制的曲线及相应坐标轴特性保持不变。
3、多窗口绘图
若需要在多个图形窗口绘制曲线,可使用“figure”命令。
第 22 页/共 63页
格式:figure(N) 创建绘图窗口,N为所创建窗口序号 例:已知函数y1(x)=sin2x,y2(x)=-15xsin2x,且x?[0,π],要求分别在两个图形窗口绘制y1(x)及 y2(x)曲线。 解:>> x=0:2*pi/90:2*pi; >> y1=sin(2*x); >> plot(x,y1,'r'); >> figure(2)
>> y2=exp(-15*x).*sin(2*x); >> plot(x,y2,'r')
4、图形窗口的分割
允许用户在同一个图形窗口里同时显示多幅独立的子图。使用函数subplot( )。格式: subplot(m,n,k) 使m×n幅子图中的第k幅成为当前图 subplot(mnk) subplot(m,n,k)的简化形式
subplot(‘position’,[left bottom width height]) 在指定位置上分割子图,并成为当前图
说明:①subplot(m,n,k)的含义是:将同一个图形窗口分割为m行×n列子窗口,k是子图的编号。 ②函数subplot( )产生的子图彼此之间独立,所有的绘图命令都可以在子图中使用。 ③使用函数subplot( )后,如果再想绘制图形窗口的单幅图,则应先使用clf命令。 ④k不能大于m与n之和。
例:在一个图形窗口中绘制函数y1=sinx,y2=sin(10x)及y12=y1y2的图形,给定x?[0,π]。 解:>> x=pi*(0:1000)/1000; >> y1=sin(x); >> y2=sin(10*x);
>> y12=sin(x).*sin(10*x);
>> subplot(2,2,1),plot(x,y1),axis([0,pi,-1,1]) %分割并绘制第一幅子图 >> subplot(2,2,2),plot(x,y2),axis([0,pi,-1,1]) %分割并绘制第二幅子图 >> subplot('position',[0.2,0.05,0.6,0.45]) %分割第三幅子图位置 >> plot(x,y12),axis([0,pi,-1,1]) %绘制第三幅子图
第 23 页/共 63页
二、图形注释
对图形进行注释,主要有三种方法:使用MATLAB图形标注函数,使用图形注释工具及使用图形窗口菜单“Insert”的注释命令。
1、使用MATLAB函数进行图形注释 函数名 title xlable ylable zlable 功能 为图形添加标题 为x轴添加标注 为y轴添加标注 为z轴添加标注 函数名 legend grid text gtext 功能 为图形添加图例 为图形坐标添加网格 在指定位置添加文本字符串 用鼠标在图形上放置文本 annontation 为图形创建特殊注释 colorbar 为图形添加颜色条 说明:①特殊注释包括:线型、箭头、文本箭头、文本框、矩形及椭圆。
②函数grid( ) 的调用格式如下:
grid on 添加坐标网络(on可以省略) grid off 去掉坐标网络
grid minor 添加更细化的坐标网络
例:为正弦函数y(x)=sinx,x?[-3π,3π]的图形添加标题、坐标轴标注及坐标网络。 解:>> x=linspace(-3*pi,3*pi,200); >> y=sin(x); >> plot(x,y)
>> title('正弦波')
>> xlabel('t'),ylabel('U') >> grid on >> grid minor
第 24 页/共 63页
2、使用属性编辑器为图形添加标题、坐标轴说明及坐标网络
属性编辑器是一种最常用的图形注释工具。打开属性编辑器之前,必须先激活图形编辑状态。 (1)激活图形编辑状态的方法
选择图形窗口菜单“Tools|Edit Plot”,或左键单击工具栏图标。 (2)打开图形编辑器的方法 ①左键双击图形窗口内区域;
②右键单击图形窗口内区域,选择“Properties”项; ③选择图形窗口菜单“View|Property Editor”。
(3)为图形添加标题、坐标轴说明及坐标网格 根据编辑器说明添加。 3、为图形添加图例 有以下三种方法:
(1)使用函数legend;(2)使用图形窗口工具栏图标;(3)选择图形窗口菜单“Insert|Legend”。 三、特殊坐标绘图 1、对数坐标绘图
常用来绘制对数频率特性曲线(即Bode图)。使用函数semilogx( )、semilogy( )和loglog( )完成。格式: semilogx(x,y) 以x轴为对数坐标,绘制(x,y)曲线(即半对数坐标) semilogy(x,y) 以y轴为对数坐标,绘制(x,y)曲线
loglog(x,y) 以x轴、y轴为对数坐标,绘制(x,y)曲线
说明:①上述函数实现的对数坐标,是指以10为底的对数坐标。 ②通常,坐标网格线函数grid( )与上述三个函数配合使用。 例:已知数组x=y=1:0.5:1000,使用函数绘制其曲线。 解:①使用函数semilogx( )绘制曲线 >> x=1:0.5:1000;y=1:0.5:1000; >> semilogx(x,y)
第 25 页/共 63页
>> xlabel('x'),ylabel('y') >> grid on
②使用函数semilogy( )绘制曲线 >> x=1:0.5:1000;y=1:0.5:1000; >> semilogy(x,y)
>> xlabel('x'),ylabel('y') >> grid on
③使用函数loglog( )绘制曲线 >> x=1:0.5:1000;y=1:0.5:1000; >> loglog(x,y)
>> xlabel('x'),ylabel('y')
第 26 页/共 63页
>> grid on
2、极坐标绘图
幅相频率特性曲线(即Nyquist图)。使用函数polar( )完成。格式:
polar(theta,rho,’s’) 绘制角度向量为theta、幅值向量为rho的极坐标图 例:使用函数polar( )绘制八叶玫瑰曲线。 解:>> theta=0:0.01:2*pi;
>> polar(theta,sin(2*theta).*cos(2*theta),'--r')
3、双纵坐标绘图
使用函数plotyy( )。格式:
plotyy(x1,y1,x2,y2) 在一个图形窗口以左右不同纵轴同时绘制(x1,y1),(x2,y2)两条曲线 plotyy(x1,y1,x2,y2,’fun’) 以左右不同纵轴将(x1,y1),(x2,y2)绘制成fun指定形式的两条曲线
第 27 页/共 63页
plotyy(x1,y1,x2,y2,’fun1’,’fun2’) 以左右不同纵轴将(x1,y1)绘制成fun1指定形式、将(x2,y2)绘制成fun2指形
式的两条曲线
说明:①左纵轴用于(x1,y1)数据对,右纵轴用于(x2,y2)数据对。
②fun、fun1可以是函数句柄或是诸如plot()、semilogx()及stem()等指定的二维绘图函数。 例:已知函数y1(x)=200e?0.05xsinx,y2(x)=0.8e?0.5xsin10x,且x?[0,20],使用函数plotyy()绘制曲线。
解:>> x=0:0.01:20;
>> y1=200*exp(-0.05*x).*sin(x); >> y2=0.8*exp(-0.5*x).*sin(10*x); >> plotyy(x,y1,x,y2)
四、三维绘图 1、三维曲线绘图
使用函数plot3(),格式:
plot3(x,y,z,’s’) 绘制三维曲线,x、y、z分别为三维坐标向量
plot3(X,Y,Z) 绘制多条三维曲线,X、Y、Z分别为三维坐标矩阵
plot3(x1,y1,z1,’s1’,x2,y2,z2,’s2’,?) 以(x,y,z,s)(四元组)结构绘制三维曲线 例:绘制三维柱面螺旋线。 解:>> t=0:pi/50:10*pi; >> plot3(sin(t),cos(t),t) >> grid on
第 28 页/共 63页
例:已知函数z(x,y)=xe?(x2?y2),且x?[-2,2],y?[-2,2],试绘制z(x,y)曲线。
解:>> [X,Y]=meshgrid(-2:.1:2); %生成x坐标与y坐标矩阵数据 >> Z=X.*exp(-X.^2-Y.^2); %生成z坐标矩阵数据 >> plot3(X,Y,Z)
2、三维曲面网线绘图
使用函数mesh(),调用格式如下:
mesh(X,Y,Z) 绘制三维曲面网线,X、Y、Z分别为三维空间坐标位置矩阵mesh(x,y,Z) 绘制三维曲面网线,以向量x,y取代矩阵X,Y 例:三维曲面网线绘图演示。 解:>> [X,Y]=meshgrid(-8:0.5:8);
第 29 页/共 63页
>> R=sqrt(X.^2+Y.^2)+eps; >> Z=sin(R)./R; >> mesh(X,Y,Z)
3、三维曲面绘图
使用函数surf()完成。格式:surf(X,Y,Z) 例:三维曲面绘图演示。
解:>> [X,Y]=meshgrid(-8:0.5:8); >> R=sqrt(X.^2+Y.^2)+eps; >> Z=sin(R)./R; >> surf(X,Y,Z)
第 30 页/共 63页
五、句柄图形简介
句柄图形(Handle Graphics)是指MATLAB使用的图形对象系统,它用于实现图形绘制和可视化函数。 1、句柄图形的概念
句柄图形是基于对象的概念,即图形的每一部分都是一个对象,每一个对象都存在与其相关的一系列句柄,而每一个对象可根据需要改变属性。
例:建立图形窗口对象演示。
解:h=figure('color','white','toolbar','none');
当用户调用一个绘图函数时,MATLAB就创建了使用各种图形对象的图形。而创建一个图形对象的同时,MATLAB也为该对象指定了一个标识符(identifier),这就是句柄。如上例中的h。
2、句柄图形的应用
用户可以使用句柄,并通过函数get()和函数set()访问对象的属性。格式: get(handle) 获得对象的全部属性
get(handle,’ ‘) 获得对象的指定属性“PropertyName”
set(handle,’PropertyName’,Property Value) 设置PropertyName的属性为Property Value 说明:①handle为图形对象句柄,PropertyName为属性名。 例:句柄图形应用演示。 解:>> x=1:10; >> y=x.^3; >> h=plot(x,y);
>> set(h,'Color','red') >> get(h,'LineWidth') ans =
0.5000
2.6 程序设计基础
一、M文件 1、M文件特点
M文件有两种形式,即M脚本文件(Script File)和M函数文件(Function File),扩展名均为.m。 2、M脚本文件
第 31 页/共 63页
脚本文件时一种简单的M文件,它没有输入、输出参数,而是包含了一系列MATLAB命令的集合,类似于DOS下的批处理文件。
例:通过M脚本文件,绘制一幅“花瓣图案”。 解:(1)编写M脚本文件。
(2)运行M脚本文件。有两种运行方式:
①选择M文件编辑/加调试器窗口菜单“Debug|Run”。
②确保M脚本文件所在路径为当前路径,然后在MATLAB命令窗口中输入命令(即文件名):e1
3、M函数文件
(1)M函数文件的概念
如果M文件的第一行包含function,那么这个文件就是M函数文件。每一个M函数文件都定义了一个函数。实际上,MATLAB提供的函数命令大部分都是由M函数文件定义的。
M函数文件比M脚本文件要复杂一些。从使用的角度来看,M函数文件是一个“黑箱”,一些数据被送进去并进行加工处理后,结果又被送出来。从形式上看,M函数文件与M脚本文件的区别在于:M函数文件的变量可以定义,但变量及其运算都仅在M函数文件内部起作用,而不在工作空间起作用,并且当M函数文件执行完后,这些内部变量将被清除。
(2)M函数文件的基本结构
第 32 页/共 63页
组成 函数定义行 H1行 函数体 格式 %注释说明 一条或若干条MATLAB命令 备注 可省略 可省略 function[输出变量列表]=函数名(输入变量列表) 函数帮助文本 %注释说明 M函数文件通常由四部分组成,具体如下:
①函数定义行,位于函数文件的首行,以MATLAB关键字function开头,函数名以及函数的输入变量和输出变量(返回变量)都在该行被定义。
②H1行,紧随函数定义行之后以%开头的第一注释行,用来概要说明该函数的功能。该行提供lookfor关键词查询和help在线帮助使用。
③函数帮助文本,位于H1行之后及函数体之前,可以有多行,每行均以%开头,用来对该函数进行注释,通常包括函数输入、输出变量的含义,函数调用格式说明,函数开发与修改的日期等。
④函数体,是函数的主要部分,由实现该M函数文件功能的MATLAB命令构成。它接收输入变量,进行程序控制,得到输出变量。
(3)M函数文件规则
①函数名必须与文件名相同,且为运行方便,脚本文件与所调用的函数文件最好放在同一文件夹内。 ②函数文件名必须以字母开头,后面可以是字母、下划线以及数字的任意组合,但不得超过31个字符。
③函数文件可以有零个或多个输入变量,也可以有零个或多个输出变量,对函数进行调用时,不能多于函数中规定的输入和输出变量个数。当函数有一个以上的输出变量时,输出变量必须包含在“[]”内,且变量之间以逗号分隔。
④函数输入和输出变量的实际个数分别 由MATLAB的两个预定义变量nargin和nargout给出,只要进入该函数,不论是否直接使用这两个变量,MATLAB都将自动生成这两个变量。
⑤函数文件中的所有变量除了事先进行特别声明外,都是局部变量。若需要使用全局变量,可使用函数global()来定义,而且在任何使用该全局变量的函数中都应该加以定义,即使在命令窗口也不例外。
(4)函数调用语句
与M脚本文件不同的是,M函数文件不能直接调用,而必须使用MATLAB的函数调用语句,该语句的基本结构为:
[输出变量列表]=函数名(输入变量列表)
例:使用M函数文件编写一个绘制任意半径和任意色彩线型的圆,并调用此函数。 解:(1)编写M函数文件。文件名为e2.m,内容如下: function sa=e2(r,s) %定义一个名为e2的函数
%circle %绘制一个以r为半径、线条属性由s定义的圆 %r %给定半径的数值
%s %给定曲线颜色的字符串 %sa %圆面积 %
a(r) %绘制半径为r的蓝色实线圆周
a(r,s) %利用字符串s给定的曲线颜色绘制半径为r的圆周 %sa=e2(r) %绘制半径为r的蓝色圆面,并计算圆面积
%sa=e2(r,s) %利用字符串s给定的曲线颜色绘制半径为r的圆面,并计算圆面积 if nargin>2
error('Too many input arguments.'); %错误信息 end
if nargin==1
第 33 页/共 63页
s='b'; end clf;
t=0:pi/100:2*pi; x=r*exp(i*t); if nargout==0 plot(x,s); else
sa=pi*r*r;
fill(real(x),imag(x),s) end
axis('square')
(2)保存M函数。将上述函数以e2.m为名保存。 (3)调用e2函数。可采用以下方法: 在MATLAB命令窗口中输入:
>> e2(2) %绘制半径为2的蓝色实线圆周 >> e2(2,'g') %绘制半径为2的绿色实线圆周
>> sa=e2(2) %绘制半径为2的蓝色圆面,并计算圆面积 sa =
12.5664
>> sa=e2(2,'g') %绘制半径为2的绿色圆面,并计算面积 sa =
12.5664
二、M文件编辑/调试器
1、启动M文件编辑/调试器的方法
(1)在MATLAB命令窗口中运行命令edit。 (2)左键单击工具栏图标。
(3)选择菜单“File|New|M=file”。 2、打开已建立的M文件的方法 (1)命令窗口中输入:edit 文件名
(2)左键单击工具栏上打开图标,再从弹出的对话框中选择所需打开的文件。 (3)选择菜单“File|Open”。
3、编写或修改后的M文件的保存方法
第 34 页/共 63页
单击工具栏保存图标,或菜单“File|Save”。 三、程序设计基础
1、全局变量和局部变量
全局变量的作用域是整个工作空间,通过命令“global”来定义,格式为:global x y z
如果想在不同的函数和MATLAB工作空间里使用同一个变量,就可以定义全局变量。如果一个变量为全局变量,则在任何一个函数里都可以对它进行赋值操作。
使用全局变量须遵循如下规则:
(1)全局变量的定义语句必须在使用该变量的语句前,为了提高程序的可读性,最好将所有全局变量的定义放在MATLAB程序的前面。
(2)一个函数如果需要调用某个全局变量,则在该函数中必须将这个变量定义为全局变量。
局部变量的作用域为函数文件所在的区域,其他函数文件无法对它实行调用。局部变量仅在其所在的函数文件运行时起作用,该函数文件运行完毕,局部变量自动消失。
2、数据类型 数据类型 数值型(Number Types) 逻辑型(Logical Types) 字符和字符串 (Characters and Strings) 日期和时间 (Dates and Times) 结构体(Structures) 元胞数组(Cell Arrays) 基本描述 包括:整型和实型、复数、不定值(NAN)、无穷大以及数据显示格式 有逻辑真(true)和逻辑假(false)两种状态 包括:字符、字符串、字符串元胞数组;字符串的比较、搜寻、替换 以及字符/数值转换 包括日期字符串,连续日期数、日期向量、日期类型转换以及输出显示格式 与C语言的结构体类似,可以存储不同类型的数据 矩阵的直接扩展,可以存储不同类型和规模的数组 函数句柄(Function Handles) 用于间接访问函数的句柄;可以很方便地调用其他函数 MATLAB类(MATLAB Classes) 使用面向对象类和方法,创建用户自己的MATLAB数据类型 Java类(Java Classes) 使用Java程序设计语言生成Java类 3、流程控制语句
主要有循环结构语句、条件(分支)结构语句、开关结构语句及试探性结构语句等。 (1)循环结构 ①for循环。格式:
for 循环控制变量=<循环次数设定> 循环体 end 说明:A.设定循环次数的数组可以是已定义数组,也可以是在for循环语句中新定义的数组,且定义的格式为:<初始值>:<步长>:<终值> 步长的缺省值为1
B.for循环可以嵌套使用。
例:使用for循环结构求
?k的值。
k?1100解:程序名为e3.m的MATLAB程序如下: sum=0; for k=1:100
sum=sum+k; end sum
第 35 页/共 63页
例:for嵌套应用实例。
解:程序名为e4.m的MATLAB程序如下: for m=1:5 for n=1:5
a(m,n)=1/(m+n-1); end end a
②while循环。格式: while <循环判断语句> 循环体 end 说明:与for循环不同的是,while循环不能指定循环的次数,循环判断语句为某种形式的逻辑判断表达式。当该表达式的逻辑值为真时,就执行循环体内的语句;当表达式的逻辑值为假时,就退出当前的循环体。
例:使用while循环结构求
?k的值。
k?1100解:程序名为e5.m的MATLAB程序如下: sum=0; k=1;
while k<=100; sum=sum+k; k=k+1; end sum
(2)分支结构
分支结构为MATLAB的条件判断语句,分为以下三种结构。 ①if-end结构。格式: if <逻辑判断语句> 执行语句 end 说明:当逻辑判断表达式为“真”时,就执行if和end之间的语句,否则不予执行。 ②if-else-end结构。格式: if <逻辑判断语句> 执行语句1 else 执行语句2 end 说明:当逻辑判断表达式为“真”时,将执行if与else内的语句,否则将执行else与end内的语句。 例:使用if-else-end结构将一数组做特殊排列。 解:程序名为e6.m的MATLAB程序如下:
第 36 页/共 63页
for k=1:9 if k<=5
b(k)=k; else
b(k)=10-k; end end b
③if-elseif-end结构。格式: if <逻辑判断语句1> 执行语句1 elseif <逻辑判断语句2> 执行语句2 else 执行语句3 end 例:使用if-elseif-end结构将矩阵a=[1 2 3;4 5 6;7 8 9]进行满足一定条件的处理。 解:程序名为e7.m的MATLAB程序如下: a=[1 2 3;4 5 6;7 8 9]; m=2;n=3; if m==n
a(m,n)=0; elseif abs(m-n)==2 a((m-1),(n-1))=-1; else
a(m,n)=-5; end a
(3)开关结构
开关结构即switch-case结构,用来解决多分支判断选择,格式为: switch <选择判断量> case 选择判断值1 选择判断语句1 case 选择判断值2 选择判断语句2 · · · otherwise 判断执行语句 end 说明:①选择判断量给出了switch-case语句的开关条件,当选择判断值与之匹配时,就执行其后的语句,如果没有选择判断值与之匹配,就执行otherwise后面的语句。在执行过程中,只有一个case命令被
第 37 页/共 63页
执行。
②执行完命令后,程序就跳出分支结构,执行end下面的语句。 例:switch-case分支结构应用实例。
解:程序名为e8.m的MATLAB程序如下: input_num=1; switch input_num case -1
disp('negative one'); case 0
disp('zero') case 1
disp('positive one'); otherwise
disp('other value'); end
(4)试探结构
试探结构try-catch为用户提供了一种错误捕获机制。格式:
try 语句段1 catch 语句段2 end 说明:试探结构首先试探性地执行语句段1,如果在此段语句执行过程中出现错误,可调用MATLAB函数lasterr()查询出错原因,并终止这段语句的执行,转而执行语句段2中的内容。如果函数lasterr()运行结果为一个空串,则表明语句段1被成功执行。若执行语句段2时又出错,则MATLAB将终止该结构。
例:试探结构应用实例。
解:程序名为e9.m的MATLAB程序如下: A=magic(3) %设置3×3维摩矩阵
B=ones(4,3) %设置4×3维元素全为1的矩阵 try
C=A*B; %取A的第N行元素 catch
C=NaN %如果C=A*B运算出错,改为C=NaN运算 end lasterr %显示出错原因
通过运行结果可以看出:在此例中,由于A和B的维数不兼容,所以在C=A*B段语句执行过程出现错误,调用lasterr()查询出错原因,并终止这段语句执行,转而执行C=NaN语句。
4、控制程序流程的其他常用命令
(1)break命令。作用是中断循环语句的执行。
(2)return命令。作用是中断函数的运行,返回到上级调用函数。
(3)pause命令。用于使程序暂时终止运行,等待用户按任意键后继续运行,适用于在调试程序时需要查看中间结果的情况。调用格式如下:
pause 暂停执行文件,等待用户按任意键继续 pause(n) 在继续执行之前暂停n秒
第 38 页/共 63页
(4)input命令。带有询问提示的输入 命令。该命令将MATLAB的控制权暂时交给用户,此后,用户通过键盘键入数值、字符串或表达式,并经回车把键入内容输入工作空间,同时把控制权交还给MATLAB。格式:
v=input(‘message’) 将用户输入的内容赋给变量v
v=input(‘message’,s) 将用户输入的内容作为字符串赋给变量v 说明:命令中的message是将要显示在屏幕中上的字符串。
(5)keyboard命令。调用键盘命令,当程序遇到keyboard时,MATLAB将控制权交给键盘,用户可以从键盘输入各种合法的MATLAB命令,只有当用户使用return命令结束输入后,控制权还给程序。
说明:keyboard与input的不同点在于:keyboard允许输入任意多个MATLAB命令,而input只能输入赋给变量的“值”,即数值、字符串或元胞数组等。
(6)警示命令。在编写M文件时,常用的警示命令有:
error(‘message’) 显示出错信息message,终止程序
lasterr 显示MATLAB自动判断的最新出错原因并终止程序 warning(‘message’) 显示警告信息message并继续运行
lastwarn 显示MATLAB自动给出的最新警告程序并继续运行
第3章 数学模型的MATLAB描述
3.1 控制系统的数学模型
一、线性定常连续系统 1、微分方程模型
设单输入单输出(SISO)线性定常连续系统的输入信号为r(t),输出信号为c(t),则其微分方程的一般形式为
dnc(t)dn?1c(t)dc(t)dmr(t)dm?1r(t)dr(t)a0?a?…?a?a?b?b?…+b?bm 1n?1n01m?1nn?1mm?1dtdtdtdtdtdt式中,系数a0,a1,?,an,b0,b1,?,bm为实常数,且m≤n。 2、传递函数(Transfer Function:TF)模型
对上式在零初始条件下求拉氏变换,并根据传递函数定义可得SISO系统传递函数的一般形式为
l[c(t)]C(s)b0sm?b1sm?1?…?bm?1s?bmM(s) G(s)????nn?1l[r(t)]R(s)a0s?a1s?…an?1s?anN(s)分子和分母多项式系数分别用向量num和den表示。
3、零极点增益(Zero-Pole-Gain:ZPK)模型
上式中分子多项式和分母多项式因式分解后,可写为如下形式:
(s?z1)(s?z2)…(s?zm)G(s)?K?K(s?p1)(s?p2)…(s?pn)?(s?z)im?(s?p)jj?1i?1n
对于SISO系统,z1,z2,?,zm为G(s)的零点,p1,p2,?,pn为G(s)的极点,K为系统的增益。 在MATLAB中,控制系统的零点和极点分别用向量Z和P表示,即 Z=[z1,z2,?,zm],P=[p1,p2,?,pn]
4、频率响应数据(Frequency Response Data:FRD)模型
设线性定常系统的频率特性为G(jω)=| G(jω)|∠ G(jω),在幅值为1,频率为ωi(i=1,2,?,n)的
第 39 页/共 63页
正弦信号r(t)=sinωit的作用下,其稳态输出为yi(t)= | G(jωi)|sin(ωit+∠ G(jωi))。频率响应数据模型就是以{ G(jωi),ωi}的形式,存储通过仿真或实验方法获得的频率响应数据值的。
5、状态空间(State-Space:SS)模型
对于多输入多输出(MIMO)系统,应用最多的是状态空间模型。一般形式为:
?(t)?Ax(t)?Bu(t)?x;u(t)为输入向量(p维);y(t)为输出向量(q维);?式中,x(t)为状态向量(n维)
y(t)?Cx(t)?Du(t)?A为系统矩阵或状态矩阵或系数矩阵(n×n维);B为控制矩阵或输入矩阵(n×p维);C为观测矩阵或输出矩阵(q×n维);D为前馈矩阵或输入/输出矩阵(q×p维)。
二、线性定常离散系统 1、差分方程模型
设SISO线性定常离散系统的输入序列为r(k),输出序列为c(k),则其差分方程的一般形式为
a0c(k?n)?a1c(k?n?1)?…?an?1c(k?1)?anc(k)?b0r(k?m)?b1r(k?m?1)?…?bm?1r(k?1)?bmr(k)
2、脉冲传递函数模型
脉冲传递函数也成为Z传递函数。SISO系统脉冲传递函数一般形式为
l[c(z)]C(z)b0zm?b1zm?1?…?bm?1z?bm G(z)???l[r(z)]R(z)a0zn?a1zn?1?…an?1z?an3、零极点增益模型
G(z)?K(z?z1)(z?z2)…(z?zm)
(z?p1)(z?p2)…(z?pn)4、状态空间模型
MIMO线性定常离散系统状态空间模型的一般形式为:
x(k?1)?Ax(k)?Bu(k)?;u(k)为输入向量序列(p维);y(k)为输?式中,x(k)为状态向量序列(n维)
y(k)?Cx(k)?Du(k)?出向量序列(q维);矩阵A、B、C、D的维数和意义与定常连续系统中相同。
3.2 数学模型的建立
线性定常系统数学模型的生成及转换函数 函数名称 功能 tf 生成(或转换)传递函数模型 ss 生成(或转换)状态空间模型 zpk frd 生成(或转换)零极点增益模型 建立频率响应数据模型 一、传递函数模型
功能:生成线性定常连续/离散系统的传递函数模型,或者将状态空间模型或零极点增益模型转换成传递函数模型。
格式:
sys=tf(num,den) 生成传递函数模型sys
sys=tf(num, den,'Property1',Value1,?, 'PropertyN',ValueN)
生成传递函数模型sys。模型sys的属性(Property)及属性值(Value)用'Property',Value指定
第 40 页/共 63页
sys=tf(num,den,Ts) 生成离散时间系统的脉冲传递函数模型sys sys=tf(num, den,Ts,'Property1',Value1,?, 'PropertyN',ValueN)
生成离散时间系统的脉冲传递函数模型sys。 sys=tf('s') 指定传递函数模型以拉氏变换算子s为自变量
sys=tf('z',Ts) 指定脉冲传递函数模型以Z变换算子z为自变量,以Ts为采样周期 tfsys=tf(sys) 将任意线性定常系统sys转换为传递函数模型tfsys
说明:①对于SISO系统,num和den分别为传递函数的分子和分母向量;对于MIMO系统,num和den为行向量的元胞数组,其行数与输出向量的维数相同,列数与输入向量的维数相同。
②Ts为采样周期,若系统的采样周期未定义,则设置Ts=-1或Ts=[]。
③缺省情况下,生成连续时间系统的传递函数模型,且以拉氏变换算子s为自变量。
s2?3s?2例:已知控制系统的传递函数为G(s)?3,用MATLAB建立其数学模型。 2s?5s?7s?3解:(1)生成连续传递函数模型。在MATLAB命令窗口中输入: >> num=[1 3 2]; >> den=[1 5 7 3]; >> sys=tf(num,den) 运行结果为: Transfer function: s^2 + 3 s + 2 --------------------- s^3 + 5 s^2 + 7 s + 3
(2)直接生成传递函数模型。在MATLAB命令窗口中输入: >>sys=tf([1 3 2],[1 5 7 3])
运行结果为: Transfer function: s^2 + 3 s + 2 --------------------- s^3 + 5 s^2 + 7 s + 3
(3)建立传递函数模型并指定输出变量名称和输入变量名称。在MATLAB命令窗口中输入: >>sys=tf(num,den,'InputName','输入端','OutputName','输出端') 运行结果为:
Transfer function from input \输入端\输出端\ s^2 + 3 s + 2 --------------------- s^3 + 5 s^2 + 7 s + 3
(4)生成离散传递函数模型(指定采样周期为0.1s)。在MATLAB命令窗口中输入: >>num=[1 3 2]; >> den=[1 5 7 3];
>> sys=tf(num,den,0.1) %指定采样周期为0.1s,缺省自变量为z 运行结果为: Transfer function: z^2 + 3 z + 2 --------------------- z^3 + 5 z^2 + 7 z + 3
第 41 页/共 63页
Sampling time: 0.1
(5)生成离散传递函数模型(未指定采样周期)。在MATLAB命令窗口中输入: >> sys=tf(num,den,-1) %生成未指定采样周期的离散系统数学模型 运行结果为: Transfer function: z^2 + 3 z + 2 --------------------- z^3 + 5 z^2 + 7 z + 3
Sampling time: unspecified
(6)生成离散传递函数模型(指定采样周期为0.1s且按照z-1排列)。在MATLAB命令窗口中输入: >> sys=tf(num,den,0.1,'variable','z^-1') 运行结果为: Transfer function:
1 + 3 z^-1 + 2 z^-2 ----------------------------
1 + 5 z^-1 + 7 z^-2 + 3 z^-3 Sampling time: 0.1
(7)生成离散传递函数模型(指定采样周期为0.1s,按照z-1排列且延迟时间为2s)。输入: >> sys=tf(num,den,0.1,'variable','z^-1','inputdelay',2) 运行结果为: Transfer function:
1 + 3 z^-1 + 2 z^-2 z^(-2) * ----------------------------
1 + 5 z^-1 + 7 z^-2 + 3 z^-3 Sampling time: 0.1
(8)生成连续时间系统传递函数模型,指定自变量为p。输入: >> num=[1 3 2]; >> den=[1 5 7 3];
>> sys=tf(num,den,'variable','p') 运行结果为: Transfer function: p^2 + 3 p + 2 --------------------- p^3 + 5 p^2 + 7 p + 3
2s2?4s?5例:系统的传递函数为G(s)?4,用MATLAB建立其数学模型。
s?7s3?2s2?6s?6解:(1)建立连续时间系统传递函数。输入: >> s=tf('s');
>> G=(2*s^2+4*s+5)/(s^4+7*s^3+2*s^2+6*s+6) 运行结果为: Transfer function:
2 s^2 + 4 s + 5 -----------------------------
s^4 + 7 s^3 + 2 s^2 + 6 s + 6
第 42 页/共 63页
(2)建立离散时间系统传递函数。输入:
>> s=tf('z',0.1); %指定使用Z变换算子生成脉冲传递函数模型(采样周期为0.1s) >> G=(2*s^2+4*s+5)/(s^4+7*s^3+2*s^2+6*s+6) 运行结果为: Transfer function:
2 z^2 + 4 z + 5 -----------------------------
z^4 + 7 z^3 + 2 z^2 + 6 z + 6 Sampling time: 0.1
?s?1??s2?2s?2?例:设多输入多输出系统的传递函数矩阵为G(s)???,应用MATLAB建立其数学 模型。
1????s??解:两种方法:(1)分别建立传递函数矩阵中每一个传递函数模型,然后按照MATLAB生成矩阵的方
式建立其模型。输入:
>> G=[tf([1 1],[1 2 2]);tf([1],[1 0])]
运行结果为:
Transfer function from input to output... s + 1 #1: ------------- s^2 + 2 s + 2
1 #2: - s
使用tf()还可以由状态空间模型或零极点增益模型得到传递函数模型。 例:系统的零极点增益模型为G(s)?(s?0.1)(s?0.2),用MATLAB建立其传递函数模型。 2(s?0.3)解:>> z=[-0.1,-0.2];p=[-0.3,-0.3];k=1;
>> sys1=zpk(z,p,k) %建立系统的零极点增益模型 运行结果为: Zero/pole/gain: (s+0.1) (s+0.2) --------------- (s+0.3)^2
>> sys2=tf(sys1) %将零极点增益模型转换为传递函数模型 运行结果为: Transfer function: s^2 + 0.3 s + 0.02 ------------------ s^2 + 0.6 s + 0.09
说明:根据零极点增益模型求取传递函数模型时,还可以使用MATLAB的求卷积函数conv(a,b)。每次调用conv(a,b)函数只能得到两个向量a和b的卷积,但conv(a,b)函数可以嵌套使用。
第 43 页/共 63页
例:线性定常系统的零极点增益模型为G(s)?s(s?6)(s?5)用MATLAB求取其
(s?3?i4)(s?3?i4)(s?1)(s?2)传递函数模型。
解:>> num=conv(conv([1 0],[1 6]),[1 5]);
>> den=conv(conv(conv([1 3+4i],[1 3-4i]),[1 1]),[1 2]); >> sys=tf(num,den) 运行结果为: Transfer function:
s^3 + 11 s^2 + 30 s --------------------------------
s^4 + 9 s^3 + 45 s^2 + 87 s + 50
二、状态空间模型
功能:生成线性定常连续/离散系统的状态空间模型,或者将传递函数模型或零极点增益模型转换为状态空间模型。
格式:
sys=ss(a,b,c,d) 生成线性定常连续系统的状态空间模型sys。a,b,c,d分别对应系统(A,B,C,D) sys=ss(a,b,c,d,’Property1’,Value1,?,’PropertyN’,ValueN)
生成连续系统的状态空间模型sys。状态空间模型sys的属性(Property)及属性值(Value)用’Property’,Value指定
sys=ss(a,b,c,d,Ts) 生成离散系统的状态空间模型sys
sys=ss(a,b,c,d,Ts, ’Property1’,Value1,?,’PropertyN’,ValueN) 生成离散系统的状态空间模型
sys_ss=ss(sys) 将任意线性定常系统sys转换为状态空间模型
??2 ?1??1 1????xx?u???例:线性定常系统的状态空间表达式为? 1 ?1??2 ?1?,应用MATLAB建立其状态空间y?[1 0]x模型。
解:(1)建立连续时间系统状态空间模型。 >> a=[-2,-1;1,-1];b=[1,1;2,-1];c=[1,0];d=0; >> sys1=ss(a,b,c,d)
运行结果为: a =
x1 x2 x1 -2 -1 x2 1 -1 b =
u1 u2 x1 1 1 x2 2 -1 c =
x1 x2 y1 1 0 d =
第 44 页/共 63页
u1 u2 y1 0 0
Continuous-time model.
(2)建立离散时间系统状态空间模型(指定采用周期为0.1s)。 >> sys1=ss(a,b,c,d,0.1) 运行结果为: a =
x1 x2 x1 -2 -1 x2 1 -1 b =
u1 u2 x1 1 1 x2 2 -1 c =
x1 x2 y1 1 0 d =
u1 u2 y1 0 0 Sampling time: 0.1 Discrete-time model.
(3)建立状态空间模型,并指定状态变量名称、输入变量名称及输出变量名称。
>> sys=ss(a,b,c,d,0.1,'statename',{'位移' '速率'},'Inputname',{'油门位移','舵偏角'},'outputname','俯仰角') 运行结果为: a =
位移 速率 位移 -2 -1 速率 1 -1 b =
油门位移 舵偏角 位移 1 1 速率 2 -1 c =
位移 速率 俯仰角 1 0 d =
油门位移 舵偏角 俯仰角 0 0 Sampling time: 0.1 Discrete-time model.
??2 ?1??1 1?x??x???2 ?1?u,用MATLAB建立其传递函数模型。例:线性定常系统状态空间表达式为 1 ?1????y?[1 0]x?[0 1]u第 45 页/共 63页
解:>> a=[-2 -1;1 -1];b=[1 1;2 -1];c=[1 0];d=[0 1]; >> sys1=ss(a,b,c,d) %建立其状态空间模型 a =
x1 x2 x1 -2 -1 x2 1 -1 b =
u1 u2 x1 1 1 x2 2 -1 c =
x1 x2 y1 1 0 d =
u1 u2 y1 0 1
Continuous-time model.
>> sys2=tf(sys1) %将状态空间模型sys1抓换为传递函数矩阵 运行结果为:
Transfer function from input 1 to output: s - 1 ------------- s^2 + 3 s + 3
Transfer function from input 2 to output: s^2 + 4 s + 5 ------------- s^2 + 3 s + 3
s?1???s3?3s2?3s?2??,例:线性定常系统的传递函数矩阵为G(s)??应用MATLAB建立其状态空间模型。 2s?3???s2?s?1???解:>> G=[tf([1 1],[1 3 3 2]);tf([1 0 3],[1 1 1])];
>> ss(G)
运行结果为: a =
x1 x2 x3 x4 x5 x1 -3 -1.5 -1 0 0 x2 2 0 0 0 0 x3 0 1 0 0 0 x4 0 0 0 -1 -1 x5 0 0 0 1 0 b =
u1 x1 1
第 46 页/共 63页
x2 0 x3 0 x4 2 x5 0 c =
x1 x2 x3 x4 x5 y1 0 0.5 0.5 0 0 y2 0 0 0 -0.5 1 d =
u1 y1 0 y2 1
Continuous-time model.
三、零极点增益模型
功能:建立线性定常连续/离散系统的零极点增益模型,或者将传递函数模型或状态空间模型转换成零极点增益模型。
格式:
sys=zpk(z,p,k) 建立连续系统的零极点增益模型sys。z,p,k分别为系统的零点向量,极点向量和增益 sys=zpk(z,p,k, ’Property1’,Value1,?,’PropertyN’,ValueN) 建立连续系统的零极点增益模型sys sys=zpk(z,p,k,Ts) 建立离散系统的零极点增益模型sys
sys=zpk(z,p,k,Ts, ’Property1’,Value1,?,’PropertyN’,ValueN) 建立离散时间系统的零极点增益模型sys
sys=zpk(‘s’) 指定零极点增益模型以拉氏变换算子s为自变量 sys=zpk(‘z’) 指定零极点增益模型以Z变换算子为自变量
zsys=zpk(sys) 将任意线性定常系统模型sys转换为零极点增益模型 说明:若系统部包含零点(或极点),则取z=[]或p=[]
例:线性定常连续系统的传递函数为G(s)?10(s?1),用MATLAB建立其零极点增益模型。
s(s?2)(s?5)解:(1)建立连续时间系统模型。 >> z=[-1];p=[0,-2,-5];k=10; >> G=zpk(z,p,k) 运行结果为: Zero/pole/gain: 10 (s+1) ------------- s (s+2) (s+5)
(2)建立离散时间系统模型(指定采样周期为0.1s)。 >> G=zpk(z,p,k,0.1) 运行结果为: Zero/pole/gain: 10 (z+1) ------------- z (z+2) (z+5)
第 47 页/共 63页
Sampling time: 0.1
(3)建立离散时间系统模型(未指定采样周期),且设定自变量为q。 >> G=zpk(z,p,k,-1,'variable','q') 运行结果为: Zero/pole/gain: 10 q^2 (1+q) ------------- (1+2q) (1+5q)
Sampling time: unspecified
(4)建立离散时间系统模型(采样周期为0.1s),且自变量按照z-1排列。 >> G=zpk(z,p,k,0.1,'variable','z^-1') 运行结果为: Zero/pole/gain: 10 z^-2 (1+z^-1) -------------------
(1+2z^-1) (1+5z^-1) Sampling time: 0.1
说明:在建立系统的零极点增益形式数学模型时,其零点向量Z和极点向量P既可以为行向量,也可以为列向量,得到的结果相同。
?1??z?0.3??,例:已知离散系统的脉冲传递函数矩阵为G(z)??应用MATLAB建立其零
2(z?0.5)????(z?0.1?j)(z?0.1?j)??极点增益模型。
解:(1)建立离散零极点增益模型(未指定采样周期)。 >> z={[ ];-0.5};p={0.3;[0.1+i,0.1-i]};k=[1;2]; >> G=zpk(z,p,k,-1) 运行结果为:
Zero/pole/gain from input to output... 1 #1: ------- (z-0.3)
2 (z+0.5) #2: --------------------
(z^2 - 0.2z + 1.01) Sampling time: unspecified
(2)建立连续零极点增益模型(指定输入变量名称及输出变量名称)。
>> G=zpk(z,p,k,'inputname','输入变量','outputname',{'输出变量1','输出变量2'}) 运行结果为:
Zero/pole/gain from input \输入变量\ 1
输出变量1: ------- (s-0.3)
2 (s+0.5)
第 48 页/共 63页
输出变量2: -------------------- (s^2 - 0.2s + 1.01)
?10s2?20s例:线性定常连续系统的传递函数为G(s)?5,应用MATLAB建立其
s?7s4?20s3?28s2?19s?5零极点增益模型。
解:>> G=tf([-10 20 0],[1 7 20 28 19 5]); %建立传递函数模型 >> sys=zpk(G)
运行结果为: Zero/pole/gain: -10 s (s-2) -----------------------
(s+1)^3 (s^2 + 4s + 5)
四、频率响应数据模型
功能:建立频率响应数据模型或将其他线性定常系统模型转换成频率响应数据模型。 格式:
sys=frd(response,frequency) 建立频率响应数据模型sys。response为存储频率响应数据的多维元胞,
frequency为频率向量,缺省单位为弧度/秒
sys=frd(response,frequency, ’Property1’,Value1,?,’PropertyN’,ValueN) 建立频率响应数据模型sys
sys=frd(response,frequency,Ts) 建立离散系统频率响应数据模型sys sysfrd=frd(sys,frequency,’Units’,units)
将其他数学模型sys转换为频率响应数据模型,并指定frequency的单 位(’Units’)为units
说明:①频率响应数据模型可以由其他三类模型转换得到,但不能将频率响应模型转换成其他模型。 ②response为复数形式。
例:设系统的频率特性为G(j?)?0.05??ej2?,计算当频率在10?1~102rad/s之间取值时的频率响
应数据模型。
解:>> freq=1:2:100; %在1与100之间产生50个等距离频率点
>> resp=0.05*(freq).*exp(i*2*freq); %计算每一个频率点freq的G(jw)值,结果为复数形式 >> sys=frd(resp,freq)
From input 1 to:
Frequency(rad/s) output 1 ---------------- --------
1 -0.020807+0.045465i 3 0.144026-0.041912i ??
97 3.452210-3.406574i 99 -4.934302-0.393914i
Continuous-time frequency response data model. 若考虑到输入变量及输出变量名称,则输入:
>> sys=frd(resp,freq,'inputname','频率','outputname','输出值')
第 49 页/共 63页
运行结果为:
From input '频率' to:
Frequency(rad/s) 输出值 ---------------- ---
1 -0.020807+0.045465i 3 0.144026-0.041912i ??
97 3.452210-3.406574i 99 -4.934302-0.393914i
Continuous-time frequency response data model.
说明:根据频率响应数据模型,可以绘制频率响应曲线。 例:设系统的传递函数为G(s)?率响应数据模型。
解:>> sys=tf([1 1],[1 4 2 6]); >> freq=0.1:100;
>> sysfrd=frd(sys,freq)
运行结果为: From input 1 to:
Frequency(rad/s) output 1 ---------------- --------
0.1 0.168158+0.011164i 1.1 1.007206+0.193739i ??
98.1 -0.000104-0.000003i 99.1 -0.000102-0.000003i Continuous-time frequency response data model. 若将频率的单位设定为赫兹(Hz),则输入语句: >> sysfrd=frd(sys,freq,'units','Hz') 运行结果为: From input 1 to:
Frequency(Hz) output 1 ------------- --------
0.1 0.245830 + 8.604151e-002i 1.1 -0.017655 - 7.168160e-003i
??
函数frd()还可以求出MIMO系统的频率响应数据模型。
s?1?12,计算当频率在10~10rad/s之间取值时的频32s?4s?2s?6s?1???s3?4s2?2s?6??12例:设线性定常系统的传递函数矩阵为G(s)??计算当频率在10~10rad/s之?,
s?1??2??s?2s?5??间取值时的频率响应数据模型。
解:>> sys=[tf([1 1],[1 2 5]);tf([1 1],[1 4 2 6])]; >> freq=0.1:100;
第 50 页/共 63页
>> sysfrd=frd(sys,freq,'units','Hz')
运行结果为: From input 1 to:
Frequency(Hz) output 1 output 2 ------------- -------- --------
0.1 0.236747+0.071835i 0.245830 + 8.604151e-002i 1.1 0.026120-0.153159i -0.017655 - 7.168160e-003i
s?1???s3?4s2?2s?6??12例:设线性定常系统的传递函数矩阵为G(s)??计算当频率在10~10rad/s之?,
s?1??2??s?2s?5??间取值时的频率响应数据模型。
解:>> sys=[tf([1 1],[1 2 5]);tf([1 1],[1 4 2 6])]; >> freq=0.1:100;
>> sysfrd=frd(sys,freq,'units','Hz')
运行结果为: From input 1 to:
Frequency(Hz) output 1 output 2 ------------- -------- --------
0.1 0.236747+0.071835i 0.245830 + 8.604151e-002i 1.1 0.026120-0.153159i -0.017655 - 7.168160e-003i
3.3 数学模型参数的获取
模型参数的获取函数
函数名称 使用方法 tfdata [num,den]=tfdata(sys) [num,den]=tfdata(sys,’v’) [num,den,Ts]=tfdata(sys) [a,b,c,d]=ssdata(sys) [a,b,c,d,Ts]=ssdata(sys) [z,p,k]=zpkdata(sys) [z,p,k]=zpkdata(sys,’v’) [z,p,kTs,Td]=zpkdata(sys) [response,freq]=frdata(sys) [response,freq,Ts]=frdata(sys) [response,freq]=frdata(sys,’v’) 功能 得到变换后的传递函数模型参数 ssdata zpkdata 得到变换后的状态空间模型参数 得到变换后的零极点增益模型参数 frddata 得到变换后的频率响应数据模型参数 3s4?2s3?5s2?4s?6例:系统的传递函数模型为G(s)?5,用MATLAB建立其传递函数模型,
s?3s4?4s3?2s2?7s?2并获取其零点向量、极点向量和增益等参数。
解:>> num=[3 2 5 4 6]; >> den=[1 3 4 2 7 2];
>> [z,p,k]=zpkdata(tf(num,den)) 运行结果为: z =
第 51 页/共 63页
[4x1 double] p =
[5x1 double] k = 3
此时只得到多维元胞数组,要显示零点向量和极点向量,必须输入: >> z1=z{1},p1=p{1} 运行结果为: z1 =
0.4019 + 1.1965i 0.4019 - 1.1965i -0.7352 + 0.8455i -0.7352 - 0.8455i p1 =
-1.7680 + 1.2673i -1.7680 - 1.2673i 0.4176 + 1.1130i 0.4176 - 1.1130i -0.2991
也可以采用一条命令直接显示零点向量和极点向量>> [z1,p1,k]=zpkdata(tf(num,den),'v')
说明:①为了方便MIMO模型或模型数组数据的获取,缺省情况下,函数tfdata()和zpkdata()以元胞数组形式返回参数。
②对于SISO模型,可在调用函数时应用第二个输入变量’v’,指定调用该函数时返回的是向量数据而不是元胞数组。
3.4 数学模型的转换
使用模型转换函数主要进行连续时间模型与离散时间模型之间的转换及离散时间模型不同采样周期之间的转换(即重新采样)。
模型转换函数 函数 功能 名称 c2d d2c d2d tf ss zpk 函数 名称 功能 由连续时间模型转换为离散时间模型 c2dm 按照指定方式将连续时间模型转换为离散时间模型 由离散时间模型转换为连续时间模型 d2cm 按照指定方式将离散时间模型转换为连续时间模型 离散时间系统重新采样 转换为传递函数模型 转换为状态空间模型 转换为零极点增益模型 tf2ss tf2zp ss2tf 将传递函数模型抓换为状态空间模型 将传递函数模型转换为零极点增益模型 将状态空间模型转换为传递函数模型 ss2ss 状态空间模型的线性变换 ss2zp 将状态空间模型转换为零极点增益模型 zp2ss 将零极点增益模型转换为状态空间模型 zp2tf 将零极点增益模型转换为传递函数模型 一、连续时间模型转换为离散时间模型 功能:将连续时间系统离散化 格式:
sysd=c2d(sys,Ts) 以采样周期Ts将线性定常连续系统sys离散化,得到离散化后的系统sysd sysd=c2d(sys,Ts,method) 以字符串“method”指定的离散化方法将线性定常系统sys离散化,包括:
(a)’zoh’——零阶保持器 (b)’foh’——一阶保持器
(c)’tustin’——图斯汀变换 (d)’matched’——零极点匹配法
第 52 页/共 63页
说明:①未指定离散化方法是,采用零阶保持器离散化方法。
②除零极点匹配法仅支持SISO系统外,其他离散化方法即支持SISO也支持MIMO。 例:连续时间系统传递函数为G(s)?s?1?0.35se,将其按Ts=0.1s进行离散化。 2s?2s?5解:(1)以零阶保持器方法。
>> sys=tf([1 1],[1 2 5],'inputdelay',0.35) %建立传递函数模型 运行结果为: Transfer function:
s + 1 exp(-0.35*s) * -------------
s^2 + 2 s + 5
>> Gd=c2d(sys,0.1) %得到离散化模型 运行结果为: Transfer function:
0.04869 z^2 + 0.002242 z - 0.04191 z^(-3) * ----------------------------------
z^3 - 1.774 z^2 + 0.8187 z Sampling time: 0.1
(2)以一阶保持器方法。 >> Gd=c2d(sys,0.1,'foh')
运行结果为: Transfer function:
0.01228 z^3 + 0.05996 z^2 - 0.05282 z - 0.0104 z^(-3) * ----------------------------------------------
z^3 - 1.774 z^2 + 0.8187 z Sampling time: 0.1
二、离散时间模型转换为连续时间模型
功能:将线性定常离散模型转换成连续时间模型 格式:
sysc=d2c(sysd) 将线性定常离散模型sysd转换成连续时间模型sysc
sysc=d2c(sysd,method) 用字符串“method”指定的方法将sysd转换成sysc 例:线性定常离散系统的脉冲传递函数为G(z)?z?1,将其转换成连续时间模型。Ts=0.1s。
z2?z?0.3解:(1)采样零阶保持器方法。
>> sysd=tf([1 -1],[1 1 0.3],0.1) %建立脉冲传递函数模型 运行结果为: Transfer function: z - 1 ------------- z^2 + z + 0.3
Sampling time: 0.1
>> sysc=d2c(sysd) %得到连续时间模型 运行结果为: Transfer function: 121.7 s - 3.215e-012
第 53 页/共 63页
---------------------
s^2 + 12.04 s + 776.7
(2)采样图斯汀方法。 >> sysc=d2c(sysd,'tustin') 运行结果为: Transfer function: -6.667 s^2 + 133.3 s --------------------
s^2 + 93.33 s + 3067
三、离散时间系统重新采样
功能:将线性定常离散时间模型重新采样或者加入输入延迟。 格式:
sys1=d2d(sys,Ts) 将离散时间模型sys按照新的采样周期Ts重新采样,得到离散时间模型sys1 例:线性定常离散系统的脉冲传递函数为G(z)?z?1,将Ts=0.1s转换成0.5s。
z2?z?0.3解:>> sysd=tf([1 -1],[1 1 0.3],0.1) %建立需采样的离散系统 运行结果为: Transfer function: z - 1 ------------- z^2 + z + 0.3
Sampling time: 0.1
>> sys_1=d2d(sysd,0.5) %对离散系统以0.5s重新采样 运行结果为: Transfer function: 0.19 z - 0.19 ----------------------
z^2 - 0.05 z + 0.00243 Sampling time: 0.5
四、传递函数模型转换为状态空间模型
格式:[A,B,C,D]=tf2ss(num,ben) 将num和den的传递函数模型转换为状态空间模型
?2s?3??2?s?2s?1??,将其转换为状态空间模型。
例:线性定常连续系统传递函数为G(s)?2s?0.4s?1解:>> num=[0 2 3;1 2 1]; %注意:分子矩阵中必须添加0,以使矩阵中两行元素个数相等 >> den=[1 0.4 1];
>> [a,b,c,d]=tf2ss(num,den) 运行结果为: a =
-0.4000 -1.0000 1.0000 0 b = 1 0
第 54 页/共 63页
c =
2.0000 3.0000 1.6000 0 d = 0 1
五、传递函数模型转换为零极点增益模型 格式:[z,p,k]=tf2zp(num,den)
2?3z?1例:线性定常离散时间系统脉冲传递函数为G(z)?,将其转换为零极点增益模型。
1?0.4z?1?z?2解:>> num=[2 3];den=[1 0.4 1]; >> [z,p,k]=tf2zp(num,den)
运行结果为: z =
-1.5000 p =
-0.2000 + 0.9798i -0.2000 - 0.9798i k = 2
六、状态空间模型转换为传递函数模型
格式:[num,den]=ss2tf(A,B,C,D,iu) 转换为num和den,得到第iu个输入向量至全部输出之间的传递
函数(矩阵)参数
??0.7524 ?0.7268??1 ?1??x???x??0 2?u0.7268 0????例:线性定常系统的状态空间模型为,将其转换为传递函数
2.8776 0???y??x?? 0 8.9463?模型。
解:>> A=[-0.7524 -0.7268;0.7268 0];B=[1 -1;0 2]; >> C=[2.8776 0;0 8.9463];D=[0 0;0 0];
>> [num,den]=ss2tf(A,B,C,D,2) %得到第2个输入至输出之间的传递函数模型 运行结果为: num =
0 -2.8776 -4.1829 0 17.8926 6.9602 den =
1.0000 0.7524 0.5282
>> [num,den]=ss2tf(A,B,C,D,1) %得到第2个输入至输出之间的传递函数模型 运行结果为: num =
0 2.8776 -0.0000 0 -0.0000 6.5022 den =
第 55 页/共 63页
1.0000 0.7524 0.5282 七、状态空间模型转换为零极点增益模型 格式:[z,p,k]=ss2zp(A,B,C,D,iu)
??0.7524 ?0.7268??1 ?1????xx???0 2?u0.7268 0例:线性定常系统的状态空间模型为????,转换为ZPK模型。 y??2.8776 8.9463?x解:>> A=[-0.7524 -0.7268;0.7268 0];B=[1 -1;0 2]; >> C=[2.8776 8.9463];D=[0 0]; >> [z,p,k]=ss2zp(A,B,C,D,1)
运行结果为: z =
-2.2596 p =
-0.3762 + 0.6219i -0.3762 - 0.6219i k =
2.8776
即第1个输入变量至输出变量间的零极点增益模型为
G1(s)?2.8776s?2.2596
(s?0.3762?i0.6219)(s?0.3762?i0.6219)八、零极点增益模型转换为传递函数模型
格式:[num,den]=zp2tf(z,p,k) z和p为列向量 例:线性定常系统的零极点增益模型为G(s)?解:>> z=[-6 -5 0]';k=1; >> p=[-3+4i -3-4i -2 -1]'; >> [num,den]=zp2tf(z,p,k)
运行结果为: num =
0 1 11 30 0 den =
1 9 45 87 50 九、零极点增益模型转换为状态空间模型 格式:[a,b,c,d]=zp2ss(z,p,k)
例:线性定常系统的零极点增益模型为G(s)?解:>> z=[-6 -5]';p=[-3+4i -3-4i -2 -1]';k=1; >> [a,b,c,d]=zp2ss(z,p,k) 运行结果为: a =
-6.0000 -5.0000 0 0
s(s?6)(s?5),转换为TF模型。
(s?3?4i)(s?3?4i)(s?1)(s?2)(s?6)(s?5),转换为SS模型。
(s?3?4i)(s?3?4i)(s?1)(s?2)第 56 页/共 63页
5.0000 0 0 0 5.0000 1.0000 -3.0000 -1.4142 0 0 1.4142 0 b = 1 0 1 0 c =
0 0 0 0.7071 d = 0
3.5 数学模型的连接
模型连接函数
函数名称 series feedback 功能 两个状态空间模型串联 函数名称 功能 parallel 两个状态空间模型并联 两个以上模型进行添加连接 两个状态空间模型按照反馈方式连接 append connect,blkbuild 将结构图转换为状态空间模型 一、优先原则 根据连接数学明显哦形式的不同,MATLAB确定的优先层级由高到低依次是FRD模型、SS模型、ZPK模型和TF模型。也就是说,如果连接的数学模型中至少有一个系统数学模型的形式为FRD模型,则无论其他系统的数学模型是上述形式的哪一种,连接后系统的数学模型为FRD模型;如果连接的数学模型中没有FRD模型,而至少有一个为SS模型,则连接得到的模型只能是SS模型。其他以此类推。
二、串联连接
两个系统sys1,sys2连接时,如果sys1的输出量作为sys2的输入量,则为串联连接。分为SISO和MIMO系统两种形式。使用函数series()实现模型的串联连接。格式:
sys=series(sys1,sys2) 相当于sys=sys1×sys2
sys=series(sys1,sys2,y1,u2) 将sys1和sys2进行广义串联连接
说明:①sys1和sys2既可以同时是连续系统模型,也可以是具有相同Ts的离散系统模型。
②y1为sys1的输出向量中与sys2输入向量串联的向量标号,u2为sys2的输入向量中与y1输入向量串联的向量标号。
③sys1和sys2为不同形式的数学模型时,sys模型的形式根据优先原则确定。 例:设两个采样周期均为Ts=0.1s的离散系统脉冲传递函数分别为
z2?3z?210,求将它们串联后的脉冲传递函数。 G1(z)?4, G(z)?2z?3z3?5z2?7z?3(z?2)(z?3)解:根据优先原则,得到的系统数学模型形式为ZPK模型。
>> G1=tf([1 3 2],[1 3 5 7 3],0.1); >> G2=zpk([ ],[-2 -3],10,0.1); >> G=series(G1,G2) 运行结果为: Zero/pole/gain:
10 (z+2) (z+1) --------------------------------------------------------
(z+2) (z+1.869) (z+3) (z+0.6245) (z^2 + 0.5063z + 2.57)
第 57 页/共 63页
Sampling time: 0.1
说明:①SISO系统模型串联连接次序不同时,所得到的状态空间模型的系数不同,但这两个系统的输出响应是相同的。
②调用函数series()时y1和u2的确定方法:设下图中sys1的输入向量u为5维,输出向量(y1+z1)为4维;sys2的输入向量(v2+u2)为2维,输出向量y为3维。串联时,sys1中的第2个和第4个输出向量分别与sys2中第1个和第2个输入向量相连接,相应的命令为:
>> y1=[2 4];
>> y1=[2 4]; %取sys1中第2个和第4个输出向量
>> u2=[1 2]; %分别与sys2第1个和第2个输入向量连接 >> sys=series(sys1,sys2,y1,u2)
sys v2 u1 sys1 y1=u2 u1
三、并联连接
两个系统sys1和sys2连接时,如果它们具有相同的输入量,且输出量是sys1输出量和sys2输出量的代数和,则系统sys1和sys2为并联连接。使用函数parallel()实现。格式:
sys=parallel(sys1,sys2) 相当于sys=sys1+sys2
sys=parallel(sys1,sys2,u1,u2,y1,y2) 将sys1和sys2进行广义并联连接。并联后得到的模型sys的输
入向量和输出向量分别为R?[v1?u?v2]T Y?[z1?y?z2]T
例:设两个采样周期均为Ts=0.1s的离散系统的脉冲传递函数分别为
sys2 y2 z2?3z?210,求并联后的脉冲传递函数。 G1(z)?4, G(z)?2z?3z3?5z2?7z?3(z?2)(z?3)解:>> G1=tf([1 3 2],[1 3 5 7 3],0.1);
>> G2=zpk([ ],[-2 -3],10,0.1); >> G=parallel(G1,G2)
运行结果为: Zero/pole/gain:
11 (z+1.869) (z+0.6673) (z^2 + 0.9178z + 3.061) --------------------------------------------------------
(z+1.869) (z+2) (z+3) (z+0.6245) (z^2 + 0.5063z + 2.57) Sampling time: 0.1
例:设系统的传递函数矩阵分别为
1.2s?1??s?1??s?2 ?(s?1)(s?3)(s?2)(s?4)??s2?2s?1 s?2??,求并联后的SS模型。 G1(s)?? G2(s)???,1s?2?s?1s?2??? ?(s?2)(s?3)(s?3)(s?4)???s2?3s?2s2?5s?6????解:>> num={[1 2] [1 1];[1] [1 2]};
第 58 页/共 63页
>> den={[1 2 1],[1 2];[1 3 2],[1 5 6]};
>> G1=tf(num,den); %建立系统G1(s)的状态空间模型 >> z={[ ] [-1];[-1] [-2]};
>> p={[-1 -3] [-2 -4];[-2 -3] [-3 -4]}; >> k=[1.2 1;1 1];
>> G2=zpk(z,p,k); %建立系统G2(s)的状态空间模型 >> G=parallel(G1,G2,2,2,1,1) 运行结果为:
Zero/pole/gain from input 1 to output... 1 #1: ----------- (s+2) (s+1) (s+2) #2: ------- (s+1)^2 #3: 0
Zero/pole/gain from input 2 to output... (s+2) #1: ----------- (s+3) (s+2)
(s+5) (s+2) (s+1) #2: ----------------- (s+2)^2 (s+4) (s+2) #3: ----------- (s+3) (s+4)
Zero/pole/gain from input 3 to output... #1: 0
1.2 #2: ----------- (s+1) (s+3) (s+1) #3: ----------- (s+2) (s+3) 四、反馈连接
使用函数feedback()实现。格式:
sys=feedback(sys1,sys2) 负反馈连接
sys=feedback(sys1,sys2,sign) 按字符串“sign”指定的反馈方式进行反馈连接 sys=feedback(sys1,sys2,feedin,feedout,sign) 构成广义反馈连接 说明:①sys的输入向量和输出向量分别于系统sys1相同。
②字符串“sign”用以指定反馈的极性,正反馈是sign=﹢1,负反馈时sign=﹣1,且负反馈时可忽略。 例:设两个线性定常系统的传递函数分别为sys1:G1(s)?馈连接后的传递函数。 解:>> G1=tf(1,[1 2 1]);
11 sys2:G(s)?,求反2s2?2s?1s?1第 59 页/共 63页
>> G2=zpk([],[-1],1); >> G3=feedback(G1,G2) 运行结果为: Zero/pole/gain: (s+1) -------------------- (s+2) (s^2 + s + 1)
>> G4=feedback(G1,G2,+1) 运行结果为: Zero/pole/gain: (s+1) ----------------- s (s^2 + 3s + 3)
五、添加连接
使用函数append()实现。格式: sys=append(G1,G2,?,GN)
?G1 0 0?设原系统模型为G1,添加系统为G2,?,GN,则增广系统为G???0 G2 0????0 0 ???例:设有4个线性定常系统,其数学模型分别为
G1(s)?10s?5 G2(s?1)2(s)?s?2 G3(S)?5x?(t)????9.0201 17.7791???1.6943 3.2138??x(t)????0.5112 0.5362???0.002 ?1.8470??u(t)
y(t)????3.2897 2.4544???0.5476 ?0.1410??13.5009 18.0745??x(t)?????0.6459 0.2958??u(t)求取将上述4个控制系统的数学模型进行添加连接后的系统数学模型。
解:4个数学模型中没有FRD模型,有1个SS模型,所以添加连接后为SS模型。首先建立4个数学模型,然后再进行模型的添加连接。 >> G1=tf(10,[1 5]); >> G2=zpk(-1,-2,2); >> G3=5;
>> A=[-9.0201 17.7791;-1.6943 3.2138];B=[-0.5112 0.5362;-0.002 -1.8470]; >> C=[-3.2897 2.4544;-13.5009 18.0745];D=[-0.5476 -0.1410;-0.6459 0.2958]; >> G4=ss(A,B,C,D);
>> sys=append(G1,G2,G3,G4) 运行结果为: a =
x1 x2 x3 x4 x1 -5 0 0 0 x2 0 -2 0 0 x3 0 0 -9.02 17.78 x4 0 0 -1.694 3.214
第 60 页/共 63页
b =
u1 u2 u3 u4 u5 x1 4 0 0 0 0 x2 0 1.414 0 0 0 x3 0 0 0 -0.5112 0.5362 x4 0 0 0 -0.002 -1.847 c =
x1 x2 x3 x4 y1 2.5 0 0 0 y2 0 -1.414 0 0 y3 0 0 0 0 y4 0 0 -3.29 2.454 y5 0 0 -13.5 18.07 d =
u1 u2 u3 u4 u5 y1 0 0 0 0 0 y2 0 2 0 0 0 y3 0 0 5 0 0 y4 0 0 0 -0.5476 -0.141 y5 0 0 0 -0.6459 0.2958 Continuous-time model. 六、复杂模型的连接
使用函数connect()实现。格式:sysc=connect(sys,q,inputs,outputs) 说明:①sys是由结构图中全部模块组成的系统;q为连接矩阵,表示结构图中模块的连接方式;inputs和outputs是复杂系统中包含输入变量和输出变量的模块编号;blkbuild为M脚本文件,用于根据传递函数或状态空间模块结构图建立对角线型状态空间结构。
②q矩阵构成如下:其行数为结构图的全部模块数;每一行第一列元素是模块的编号,该模块输入端与结构图中一些模块的输出端连接,该行q矩阵其他元素依次为与该模块相连接的其他模块编号;元素符号根据其他模块输出端是加还是减来确定;q矩阵中其他元素均为0.
③sys的求取步骤为:第一,运行脚本文件blkbuild之前,必须按照下述要求设置输入参数:(a)nblocks为结构图的总模块数;(b)若第i个模块为以传递函数模型,则分别输入该模块分子项和分母项参数ni,di;(c)若第i个模块式状态空间模型,则分别输入该模块哥哥矩阵参数ai,bi,ci,di。第二,运行blkbuild后的返回结果为系统状态空间模型(a,b,c,d),再利用函数ss()就可以建立状态空间结构。
④在③中队结构图的每一模块设置输入参数时,如果同一个模块的ni,di和ai,bi,ci,di同时存在,则会发生错误。
例:已知控制系统的结构图。其中各环节的传递函数分别为
1.2(2s?1)0.2(4s?1),G3(s)?2s4s求系统以r(t)为输入,分别以y1(t)(局部反馈
1.210.5G4(s)?,G5(s)?,G6(s)?0.5,G7(s)?4s?10.410s?1G1(s)?1,G2(s)?回路输出)和y(t)(主反馈回路输出)为输出时的传递函数(矩阵),并绘制其单位阶跃响应曲线。
第 61 页/共 63页
r G1(s) × ○- G2(s) × ○- G3(s) G5(s) G6(s) G4(s) y1 G7(s) y
解:(1)确定连接矩阵q。
该系统结构图包含7个 模块,则q矩阵有7行;G1(s)模块的输入端没有与其他模块的输出端相连接,所以q矩阵中第1行除了第1列元素为1(模块编号)外,其他列元素均为0;G2(s)模块的输入端分别于G1(s)模块和G6(s)模块的输出端连接,考虑到G6(s)模块的输出端在比较点进行减法运算,所以q矩阵第2行第2列元素为1(对应模块G1(s)),第3列元素为-6(对应模块G6(s))。据此,可以确定q矩阵如下:
q=[1 0 0;2 1 -6;3 2 -5;4 3 0;5 4 0;6 7 0;7 4 0]
(2)MATLAB求解。
>> nblocks=7; %共有7个模块
>> n1=1;d1=1; %第1个模块分子项和分母项参数(多项式形式) >> n2=1.2*[2 1];d2=[2 0]; >> n3=0.2*[4 1];d3=[4 0]; >> n4=1.2;d4=[4 1]; >> n5=1;d5=0.4; >> n6=0.5;d6=1; >> n7=0.5;d7=[10 1]; >> blkbuild;
State model [a,b,c,d] of the block diagram has 7 inputs and 7 outputs. 得到了对角线型模块,继续输入: >> sys=ss(a,b,c,d);
>> q=[1 0 0;2 1 -6;3 2 -5; 4 3 0;5 4 0;6 7 0;7 4 0]; >> inputs=1; >> outputs=[4 7];
>> sysc=connect(sys,q,inputs,outputs) 运行结果为: a =
x1 x2 x3 x4 x1 0 0 0 -0.025 x2 0.6 0 -0.75 -0.03 x3 0.12 0.05 -0.4 -0.006 x4 0 0 0.3 -0.1 b =
u1 x1 1 x2 1.2 x3 0.24 x4 0 c =
x1 x2 x3 x4
第 62 页/共 63页
y1 0 0 0.3 0 y2 0 0 0 0.05 d =
u1 y1 0 y2 0
Continuous-time model.
(3)求取系统的传递函数矩阵。 >> G=tf(sysc)
运行结果为:
Transfer function from input to output...
0.072 s^3 + 0.0612 s^2 + 0.0144 s + 0.0009 #1: ------------------------------------------------
s^4 + 0.5 s^3 + 0.0793 s^2 + 0.0051 s + 0.000225
0.0036 s^2 + 0.0027 s + 0.00045 #2: ------------------------------------------------
s^4 + 0.5 s^3 + 0.0793 s^2 + 0.0051 s + 0.000225 最后,求系统单位阶跃响应。 >> step(G(1),'*-',G(2)) 运行后得到的曲线如下图。
第 63 页/共 63页