东北大学薛定宇 控制系统及计算机仿真chap 9
控制系统仿真_薛定宇第四章 线性控制系统的数学模型
控制系统仿真与CAD 国家级精品课程
2014-12-31
44/138
4.3.3 控制系统的状态方程实现
由传递函数到状态方程的转换 不同状态变量选择，结果不唯一 默认变换方式，采用MATLAB函数
G可以是传递函数、状态方程和零极点模型 适用于有延迟的、离散的或多变量模型 可以将延迟传递函数模型转成内部延迟
2014-12-31
5/138
4.1.1 线性连续系统数学模型及 MATLAB 表示
线性系统的常系数线性常微分方程模型
为阶次， 为常数， 物理可实现 线性定常系统 LTI （linear time invariant）
控制系统仿真与CAD 国家级精品课程
2014-12-31
6/138
传递函数的理论基础 ——Laplace变换
2014-12-31
47/138
均衡实现 （banlanced realization）
由一般状态方程输入输出关系显著程度不 明显，需要进一步变换 均衡实现是一种很有用的方式 用MATLAB直接求解
得出均衡实现的模型 得出排序的 Gram 矩阵
控制系统仿真与CAD 国家级精品课程
2014-12-31
离散系统也有内部延迟模型
控制系统仿真与CAD 国家级精品课程
2014-12-31
33/138
4.3 系统模型的相互转换
前面介绍的各种模型之间的相互等效变换
主要内容
连续模型和离散模型的相互转换 系统传递函数的获取 控制系统的状态方程实现 状态方程的最小实现 传递函数与符号表达式的相互转换
控制系统计算机仿真课后答案
控制系统计算机仿真课后答案参考答案说明:1( 对于可以用文字或数字给出的情况，直接给出参考答案。
2( 对于难以用文字或数字给出的情况，将提供MATLAB程序或Simulink模型。
第 1 章1.1 系统是被研究的对象，模型是对系统的描述，仿真是通过模型研究系统的一种工具或手段。
1.2 数学仿真的基本工具是数字计算机，因此也称为计算机仿真或数字仿真。
将数学模型通过一定的方式转变成能在计算机上实现和运行的数学模型，称之为仿真模型。
1.3 因为仿真是在模型上做试验，是一种广义的试验。
因此，仿真基本上是一种通过试验来研究系统的综合试验技术，具有一般试验的性质。
而进行试验研究通常是需要进行试验设计。
1.4 解析法又称为分析法，它是应用数学推导、演绎去求解数学模型的方法。
仿真法是通过在模型上进行一系列试验来研究问题的方法。
利用解析法求解模型可以得出对问题的一般性答案，而仿真法的每一次运行则只能给出在特定条件下的数值解。
，解析法常常是围绕着使问题易于求解，而不是使研究方法更适合于问题，常常因为存在诸多困难而不能适用。
从原则上讲，仿真法对系统数学模型的形式及复杂程度没有限制，是广泛适用的，但当模型的复杂程度增大时，试验次数就会迅速增加，从而影响使用效率。
1.5 仿真可以应用于系统分析、系统设计、理论验证和训练仿真器等方面。
1.6,8,20,71，，，，,,,,，x,100x，0u,,,, ,,,,0100，，，，y,,,002x注:本题答案是用MATLAB中tf2ss()函数给出的，是所谓“第二能控标准型”(下同)。
11.7,3,3,11，，，，,,,,，x,100x，0u,,,, ,,,,0100，，，，y,,,013x1.82s，3s，3G(s), 32s，4s，5s，21.91.368,0.36801，，，，,,,,x(k，1),100x(k)，0u(k),,,, ,,,,0100，，，，y(k),,,00.3680.264x(k)1.10 仿真模型见praxis1_10_1.mdl;MATLAB程序见praxis1_10_2.m。
控制系统计算机辅助设计-MATLAB语言与应用
20
国际上出版了关于 MATLAB及 CACSD 的专著和教材，但它们大都是MATLAB的
入门教材，并没有真正深入、系统地探讨 CACSD 技术及 MATLAB实现，将MATLAB
的强大功能与控制领域成果有机结合是本 书力图解决的主要问题。
2020/4/12
控制系统计算机辅助设计－MATLAB语言与应用
2020/4/12
控制系统计算机辅助设计－MATLAB语言与应用
17
除了经典的多变量频域方法之外，还出现了一些 基于最优化技术的控制方法，其中比较著名的是 英国学者 John Edmunds 提出的多变量参数最优化 控制方法和英国学者 Zakian 提出的不等式控制方 法等。
与此同时，美国学者似乎更习惯于状态空间的表 示与设计方法。此方法往往又称为时域方法 (timedomain)，首先在线性二次型指标下引入了最优控 制的概念，并在用户的干预下(如人工选择加权矩 阵)得出某种最优控制的效果，这样的控制又往往 需要引入状态反馈或状态观测器新的控制概念。
辨识工具箱、鲁棒控制工具箱、多变量频域设计工 具箱、µ分析与综合工具箱、神经网络工具箱、最 优化工具箱、信号处理工具箱、以及仿真环境
Simulink。
2020/4/12
控制系统计算机辅助设计－MATLAB语言与应用
16
1.5 控制系统计算机辅助设计 领域的新方法
早期的 CACSD 研究侧重于对控制系统的计算机辅 助分析上，开始时人们利用计算机的强大功能把 系统的频率响应曲线绘制出来，并根据频率响应 的曲线及自己的控制系统设计经验用试凑的方法 设计一个控制器，然后利用仿真的方法去观察设 计的效果。
2020/4/12
控制系统计算机辅助设计－MATLAB语言与应用
控制系统仿真 薛定宇第九章 分数阶系统的分析与设计课件
n 启动optimfopid界面
控制器控制器
nn
19
分数阶控制系统小结n 给出了分数阶微积分的定义与计算方法n 微积分计算：glfdiff()、fode_sol()n 线性分数阶系统分析与设计n 类的建立：@fotf；fotf.m、display.mn 重载函数 (FOTF互连) 支持 *、+、feedback()n 时域与频域分析：重载函数尽量控制系统工具箱 函数同名，且调用方式尽可能保持一致nstep(), lsim(), bode(), nyquist(), nichols(), isstable(), norm()n 其他内容：解析解法、成比例阶系统、微分算子近似，基于框图的分数阶非线性系统 制 统 真 AD国家级精品课程
控制系统仿真与CAD 国家级精品课程
6
分数阶微积分计算举例n 函数调用n 例11-1 常数的微积分是什么？
n 例11-2 正弦函数的分数阶微分
n 信息量比整数阶微分丰富
控制系统仿真与CAD 国家级精品课程
7
11.2 分数阶线性微分方程的求解分数阶线性微分方程
控制系统仿真与CAD 国家级精品课程
nMATLAB实现
n Gottfried Leibniz的记号 n微分、积分n 扩展：统一微积分
控制系统仿真与CAD 国家级精品课程
3
分数阶微积分的历史n 法国数学家Guillaume François AntoineL’Hôpitaln 1695年，询问Leibniz， n=1/2？
n Leibniz给出了t 的1/2导数n 彼此关于此问题的通信标志分数阶微积分研究的开始n 早期研究是纯数学研究n 控制方面的研究nManabe，1960开始
n display.m 模型显示 (略) n 例11-6
控制系统仿真_薛定宇第二章_MATLAB语言程序设计基础
控制系统仿真与CAD 国家级精品课程
2014-12-31
24/70 24
2.4.4 匿名函数与inline函数
不用编写*.m，定义函数的方法
可以动态地定义函数，在科学运算中有用 二者功能重复，后者远差于前者，不再使用
匿名函数： 举例
控制系统仿真与CAD 国家级精品课程
2014-12-31
控制系统仿真与CAD 国家级精品课程
2014-12-31
3/70 3
2.1.1 常量与变量
常量（保留字符串）
eps： 机器的浮点运算误差限，默认值 i和j：虚数单位，-1的平方根，i=sqrt(-1) Inf： 无穷大量，-Inf为负无穷大 NaN：不定式（Not a Number） pi：圆周率p的双精度浮点表示
25/70 25
2.5 二维图形绘制
可以绘制出各种二维曲线 可以由数据绘图，也可以由格式绘图 本节主要内容
基本二维图形绘制语句 带有其他选项的绘图函数 二维曲线的标注方法 在MATLAB图形上添加文字标注 特殊图形绘制函数及举例 隐函数的曲线绘制
控制系统仿真与CAD 国家级精品课程
符号型数据结构
syms 命令申明符号变量
控制系统仿真与CAD 国家级精品课程
2014-12-31
5/70 5
vpa(A) 显示符号变量
例2-1 显示e的前300位 其他数据结构
int8() 、int16() 、int32() 、uint16() 和 uint32() 字符串型数据：MATLAB字符串是用单引号括起来的 多维数组：三维数组是一般矩阵的直接拓展。 单元数组：多变量系统的频域响应数据存储 类与对象：传递函数对象、状态方程对象
MATLAB 历年真题薛定宇
东 北 大 学 考 试 试 卷 2009－2010学年第一学期 课程名称：MATLAB 语言与应用（70分） 一、（11分）给出下面问题的MATLAB 语句实现。
1.(3分) 求解能转换成多项式方程的联立方程，并检验得出的高精度数值解的精度。
2222222323224--4-23+4-5+3=x y zxy x yz xz xy yz x z xz xy xy xz x zy ⎧=⎪=⎨⎪⎩ [x,y,z]=solve('x^2*y^2-z*x*y-4*x^2*y*z^2=x*z^2','x*y^2-2*y*z^2=3*x^3*z^2+4*x*z','y^3*x-5*x*y^2+3*x*z^2=x^4*z*y','x,y,z') (2分)norm(double([x.^2.*y.^2-z.*x.*y-4*x.^2.*y.*z.^2-x.*z.^2,x.*y.^2-2*y.*z.^2-3*x.^3.*z.^2-4*x.*z, y.^3.*x-5*x.*y.^2+3*x.*z.^2-x.^4.*z.*y])) (1分)2.绘制下列图形（1）（3分）请在一幅图中绘制出22/2/3sin()0xy x x e e xy -+=和222cos()0x yx x y y e -++=的曲线，其中x ∈[-10,10], y ∈[-9,9]。
ezplot('x^2*exp(-x*y^2/2)+exp(-x/3)*sin(x*y)',[-10,10,-9,9])；(1分)hold on; (1分)ezplot('y^2 *cos(y+x^2) +x^2*exp(x+y)',[-10,10,-9,9])；(1分)（2）（5分）已知函数---=--222(,)(2)n m nm g n m n m e n 。
要求使用subplot 命令生成2行2列的子窗口，并分别在第1、2、3、4个子窗口中绘制该函数的三维曲面图、俯视图、侧视图和正视图。
《控制系统计算机辅助设计MATLAB语言与应用第2版》薛定宇_课后习题答案
【17】
（1)z＝ｘy
＞>[x,y]=mesｈgrid(-3:0.01:3，-3：0.０1:３);
z=x.*y；
mesh(x,ｙ，z);
>> contoｕr3(x,ｙ,z,５0)；
(１)ｚ=sin（xｙ)
>> [x,y］=mesｈgrid(-3:０.01:3，-3:０.０1：3);
【２】
相应的ＭATＬAB命令:B=A(２：2:end,：）
>>Ａ=ｍagic(8)
A=
64２361６06７57
95５5412１3 51５01６
1７47４6 20２１43 42２4
40262737３6 30313３
3２３435 29２838 392５
41２32244４5１91８48
4９１514 5253 11 10 5６
【１0】
fｕｎcｔion y＝fｉb(ｋ)
if nargin~=1,ｅrror（＇出错：输入变量个数过多，输入变量个数只允许为１!')；ｅｎdﻭif narｇout＞１,ｅｒrｏr(＇出错：输出变量个数过多！'）;end
if k<=０,error('出错:输入序列应为正整数！＇）;endﻭiｆk==1|k＝=2,y＝１;ﻭｅｌsｅy=ｆib（k-１)＋fib(k-２);ｅnｄﻭeｎｄ
8５8 5954 62 6３1
＞＞Ｂ=Ａ(2:2:ｅnd，:)
B =
9５5 54１21３515０16
402６27３7 3６30３1３3
4１２3２24４4519１8 4８
858 5954 6２631
控制系统仿真_薛定宇第一章_控制系统计算机辅助设计概述
PID控制器的自整定方法 智能控制
自适应控制：在线辨识和调整控制器参数 模糊逻辑控制器、仿人控制、神经网络控制等
分数阶控制：非整数阶微积分在控制中应用
控制系统仿真与CAD 国家级精品课程
2014-12-31
16/17
1.5 MATLAB版本选择及学习方法
2008a的符号运算工具箱是目前的最高水平
12/17
1999 年MATLAB 5.3，全新最优化工具箱和 Simulink 3 2000 年，MATLAB 6.0，在操作界面上有了很大改观，抛 弃了LINPACK、EISPACK，改用更具优势的 LAPACK 软 件包和 FFTW 系统，速度变得更快，数值性能也更好 2004 年，MATLAB 7.0 版：多领域物理建模仿真策略 MathWorks 公司每年在 3 月和 9 月分别推出 a 版和 b 版， 2008a的符号运算功能达到顶峰 2012年9月，MATLAB 2012b暨8.0版，全新的界面，强大 的仿真功能、更多的工具箱
利用底层语言开发专用 Melsa & Jones：McGraw-Hill，1973 Karl Astrom： INTRAC（IDPAC、MODPAC、SYNPAC、 POLPAC 等，及仿真语言 SIMNON） CLADP（Cambridge linear analysis and design programs） NASA Armstrong：ORACLS（optimal regulator algorithms for the control of linear systems） 孙增圻、袁曾任：清华大学出版社，1988及早期讲义 韩京清等：自然科学基金重大项目：CACSDC
matlab概论.
• 2、按计算机类型分 • （一）、模拟仿真：采用数学模型在模拟计算机上进行的 实验研究。
• 特点：
• • • • • 1）描述连续的物理系统的动态过程比较自然逼真。 2）仿真速度极快，失真小，结果可靠。 3）受元器件性能影响，仿真精度较低。 4）对计算机控制系统（采样控制系统）的仿真较难。 5）实验过程的自动化程度低。
2018/9/24 16
• 二、控制系统CAD的主要内容 • 1、频域法：属经典控制理论范畴，适用SISO系统。借助传递函数、
劳斯判据、伯德图、乃氏图及根轨迹来分析系统动态和稳态性能，设 计系统校正装置的结构，确定参数。
• 2、时域法：属现代控制理论，适用多变量系统。包括：最优控制、
极点配置、状态反馈及观测器设计、能控性能观性分析等。
2018/9/24 5
第二节 仿真实验的分类与性能比较
• 1、按模型分类：物理仿真、数学仿真 • 物理仿真特点：总有实物介入，具有实 时性和在线性。构造复杂，造价高。 • 数学仿真特点：在（模拟或数字）计算 机上进行，具有非实时性和离线性。经 济、快捷。 • 本课主要讨论基于数学模型的数字仿真
2018/9/24 6
的创造性结合起来，用以缩短设计周期、提高设计质量的技术。并不 是简单地用计算机来代替人工计算，而是通过CAD系统与人之间的信 息交互作用，从本质上提高设计者的能力和设计结果的水平。如汽车 设计，服装设计。
• 2、CAD系统的组成：应用软件、计算机、外围设备、设计者
• 核心：应用软件，如控制系统设计的MATLAB,机械设计的AutoCAD • 基础：计算机，各种PC机，小型机、大型机等 • 外围设备：包括输入输出等设备，是人机信息交换的手段
2018/9/24
控制系统计算机辅助设计_MATLAB语言与应用(第2版)薛定宇_课后复习题答案
第1章控制系统计算机辅助设计概述第2章 MATLAB语言程序设计基础第3章线性控制系统的数学模型第4章线性控制系统的计算机辅助分析第5章 Simulink在系统仿真中的应用第6章控制系统计算机辅助设计第1章控制系统计算机辅助设计概述【1】/已阅，略【2】已阅，略【3】已经掌握help命令和Help菜单的使用方法【4】区别：MATLAB语言实现矩阵的运算非常简单迅速，且效率很高，而用其他通用语言则不然，很多通用语言所实现的矩阵运算都是对矩阵维数具有一点限制的，即使限制稍小的，但凡维数过大，就会造成运算上的溢出出错或者运算出错，甚至无法处理数据的负面结果【5】【8】(1)输入激励为正弦信号(2)输入激励为脉冲模拟信号(3)输入激励为时钟信号(4) 输入激励为随机信号(5) 输入激励为阶跃信号δ=0.3δ=0.05δ=0.7结论：随着非线性环节的死区增大，阶跃响应曲线的围逐渐被压缩，可以想象当死区δ足够大时，将不再会有任何响应产生。
所以可以得到结论，在该非线性系统中，死区的大小可以改变阶跃响应的幅值和超调量。
死区越大，幅值、超调量将越小，而调整时间几乎不受其影响第2章 MATLAB语言程序设计基础【1】>> A=[1 2 3 4;4 3 2 1;2 3 4 1;3 2 4 1]A =1 2 3 44 3 2 12 3 4 13 24 1>>B=[1+4i,2+3i,3+2i,4+i;4+i,3+2i,2+3i,1+4i;2+3i,3+2i,4+i,1+4i;3+2i,2+3i,4+i,1+4i]B =1.0000 + 4.0000i2.0000 +3.0000i 3.0000 + 2.0000i4.0000 + 1.0000i4.0000 + 1.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 1.0000 + 4.0000i2.0000 +3.0000i 3.0000 + 2.0000i4.0000 + 1.0000i 1.0000 + 4.0000i3.0000 + 2.0000i 2.0000 + 3.0000i4.0000 + 1.0000i 1.0000 + 4.0000i >> A(5,6)=5A =1 2 3 4 0 04 3 2 1 0 02 3 4 1 0 03 24 1 0 00 0 0 0 0 5∴若给出命令A(5,6)=5则矩阵A的第5行6列将会赋值为5，且其余空出部分均补上0作为新的矩阵A，此时其阶数为5×6【2】相应的MATLAB命令：B=A(2:2:end,:)>> A=magic(8)A =64 2 3 61 60 6 7 579 55 54 12 13 51 50 1617 47 46 20 21 43 42 2440 26 27 37 36 30 31 3332 34 35 29 28 38 39 2541 23 22 44 45 19 18 4849 15 14 52 53 11 10 568 58 59 5 4 62 63 1>> B=A(2:2:end,:)B =9 55 54 12 13 51 50 1640 26 27 37 36 30 31 3341 23 22 44 45 19 18 488 58 59 5 4 62 63 1∴从上面的运行结果可以看出，该命令的结果是正确的【3】>> syms x s; f=x^5+3*x^4+4*x^3+2*x^2+3*x+6f =x^5 + 3*x^4 + 4*x^3 + 2*x^2 + 3*x + 6>> [f1,m]=simple(subs(f,x,(s-1)/(s+1)))f1 =19 - (72*s^4 + 120*s^3 + 136*s^2 + 72*s + 16)/(s + 1)^5m =simplify(100)【4】>> i=0:63; s=sum(2.^sym(i))s =615【5】>> for i=1:120if(i==1|i==2) a(i)=1;else a(i)=a(i-1)+a(i-2);endif(i==120) a=sym(a); disp(a); endend[ 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049, , , , , , 5, 7, 2, 9, 1, 20, 61, 81, 42, 723, 565, 288, 853, 141, 0994, 9135, 0129, 9264, 9393, 28657, 78050, 06707, 84757, 91464, , , , , , 8, 5, 3, 8, 31, 89, 20, 09, 29, 738, 167, 905, 072, 977, 6049, 9026, 5075, 4101, 9176, 83277, 82453, 65730, 48183, 413913, 662096, 076009, 738105, 814114, 0552219, 6366333, 6918552, 3284885, 0203437, 93488322, 23691759, 17180081, 40871840]【6】>>k=1;for i=2:1000for j=2:iif rem(i,j)==0if j<i, break;endif j==i, A(k)=i; k=k+1; break; endendendenddisp(A);Columns 1 through 132 3 5 7 11 13 17 19 23 29 31 37 41 Columns 14 through 2643 47 53 59 61 67 71 73 79 83 89 97 101 Columns 27 through 39103 107 109 113 127 131 137 139 149 151 157 163 167 Columns 40 through 52173 179 181 191 193 197 199 211 223 227 229 233 239 Columns 53 through 65241 251 257 263 269 271 277 281 283 293 307 311 313 Columns 66 through 78317 331 337 347 349 353 359 367 373 379 383 389 397 Columns 79 through 91401 409 419 421 431 433 439 443 449 457 461 463 467 Columns 92 through 104479 487 491 499 503 509 521 523 541 547 557 563 569 Columns 105 through 117571 577 587 593 599 601 607 613 617 619 631 641 643 Columns 118 through 130647 653 659 661 673 677 683 691 701 709 719 727 733 Columns 131 through 143739 743 751 757 761 769 773 787 797 809 811 821 823 Columns 144 through 156827 829 839 853 857 859 863 877 881 883 887 907 911 Columns 157 through 168919 929 937 941 947 953 967 971 977 983 991 997【7】说明：h和D在MATLAB中均应赋值，否则将无法实现相应的分段函数功能syms x; h=input(‘h=’); D=input(‘D=’);y=h.*(x>D)+(h.*x/D).*(abs(x)<=D)-h.*(x<-D)【10】function y=fib(k)if nargin~=1,error('出错：输入变量个数过多，输入变量个数只允许为1！');endif nargout>1,error('出错：输出变量个数过多！');endif k<=0,error('出错：输入序列应为正整数！');endif k==1|k==2,y=1;else y=fib(k-1)+fib(k-2);endend【13】-1-0.500.51-1-0.8-0.6-0.4-0.200.20.40.60.81【14】>> t=[-1:0.001:-0.2,-0.1999:0.0001:0.1999,0.2:0.001:1]; y=sin(1./t); plot(t,y);grid on;-1-0.8-0.6-0.4-0.20.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81【15】(1) >> t=-2*pi:0.01:2*pi; r=1.0013*t.^2;polar(t,r);axis('square')90270180(2) >> t=-2*pi:0.001:2*pi; r=cos(7*t/2);polar(t,r);axis('square')2700902701800(3) >> t=-2*pi:0.001:2*pi;r=sin(t)./t;polar(t,r);axis('square')90180【17】(1)z=xy>> [x,y]=meshgrid(-3:0.01:3,-3:0.01:3);z=x.*y;mesh(x,y,z);>> contour3(x,y,z,50);-2-112-22-10-5510(1)z =sin(xy )>> [x,y]=meshgrid(-3:0.01:3,-3:0.01:3); z=sin(x.*y);mesh(x,y,z);>> contour3(x,y,z,50);-2-112-22第3章 线性控制系统的数学模型【1】(1) >> s=tf('s');G=(s^2+5*s+6)/(((s+1)^2+1)*(s+2)*(s+4)) Transfer function:s^2 + 5 s + 6--------------------------------s^4 + 8 s^3 + 22 s^2 + 28 s + 16(2) >> z=tf('z',0.1);H=5*(z-0.2)^2/(z*(z-0.4)*(z-1)*(z-0.9)+0.6) Transfer function:5 z^2 - 2 z + 0.2---------------------------------------z^4 - 2.3 z^3 + 1.66 z^2 - 0.36 z + 0.6 Sampling time (seconds): 0.1【2】(1)该方程的数学模型>> num=[6 4 2 2];den=[1 10 32 32];G=tf(num,den)Transfer function:6 s^3 + 4 s^2 + 2 s + 2------------------------s^3 + 10 s^2 + 32 s + 32(2)该模型的零极点模型>> G=zpk(G)Zero/pole/gain:6 (s+0.7839) (s^2 - 0.1172s + 0.4252)-------------------------------------(s+4)^2 (s+2)(3)由微分方程模型可以直接写出系统的传递函数模型【5】(1) >> P=[0;0;-5;-6;-i;i];Z=[-1+i;-1-i];G=zpk(Z,P,8)Zero/pole/gain:8 (s^2 + 2s + 2)-------------------------s^2 (s+5) (s+6) (s^2 + 1)(2) P=[0;0;0;0;0;8.2];Z=[-3.2;-2.6];H=zpk(Z,P,1,'Ts',0.05,'Variable','q')Zero/pole/gain:(q+3.2) (q+2.6)---------------q^5 (q-8.2)Sampling time (seconds): 0.05【8】(1)闭环系统的传递函数模型>> s=tf('s');G=10/(s+1)^3;Gpid=0.48*(1+1/(1.814*s)+0.4353*s/(1+0.4353*s));G1=feedback(Gpid*G,1)Transfer function:7.58 s^2 + 10.8 s + 4.8-------------------------------------------------------------- 0.7896 s^5 + 4.183 s^4 + 7.811 s^3 + 13.81 s^2 + 12.61 s + 4.8(2)状态方程的标准型实现>> G1=ss(G1)a =x1 x2 x3 x4 x5x1 -5.297 -2.473 -2.186 -0.9981 -0.7598x2 4 0 0 0 0x3 0 2 0 0 0x4 0 0 2 0 0x5 0 0 0 0.5 0b =u1x1 2x2 0x3 0x4 0x5 0c =x1 x2 x3 x4 x5y1 0 0 0.6 0.4273 0.3799d =u1y1 0Continuous-time state-space model.(3)零极点模型>> G1=zpk(G1)Zero/pole/gain:9.6 (s^2 + 1.424s + 0.6332)--------------------------------------------------------(s+3.591) (s^2 + 1.398s + 0.6254) (s^2 + 0.309s + 2.707)【11】>> Ga=feedback(s/(s^2+2)*1/(s+1),(4*s+2)/(s+1)^2);Gb=feedback(1/s^2,50);G=3*feedback(Gb*Ga,(s^2+2)/(s^3+14))Transfer function:3 s^6 + 6 s^5 + 3 s^4 + 42 s^3 + 84 s^2 + 42 s---------------------------------------------------------------------------s^10 + 3 s^9 + 55 s^8 + 175 s^7 + 300 s^6 + 1323 s^5 + 2656 s^4 + 3715 s^3+ 7732 s^2 + 5602 s + 1400【13】c1=feedback(G5*G4,H3)=G5*G4/(1+G5*G4*H3)c2=feedback(G3,H4*G4)=G3/(1+G3*H4*G4)c3=feedback(c2*G2,H2)=c2*G2/(1+c2*G2*H2)=G3*G2/(1+G3*H4*G4+G3*G2*H1)G=feedback(G6*c1*c3*G1,H1)=G6*c1*c3*G1/(1+ G6*c1*c3*G1*H1)=G6*G5*G4*G3*G2*G1/(1+G3*H4*G4+G3*G2*H1+G5*G4*H3+G5*G4*H3*G3*H4*G4+G5*G4*H3*G3* G2*H1+G6*G5*G4*G3*G2*G1*H1)【14】>> s=tf('s');c1=feedback(0.21/(1+0.15*s),0.212*130/s);c2=feedback(c1*70/(1+0.0067*s)*(1+0.15*s)/(0.051*s),0.1/(1+0.01*s));G=(1/(1+0.01*s))*feedback(130/s*c2*1/(1+0.01*s)*(1+0.17*s)/(0.085*s),0.0044/(1+ 0.01*s))Transfer function:0.004873 s^5 + 1.036 s^4 + 61.15 s^3 + 649.7 s^2 + 1911 s--------------------------------------------------------------------------- 4.357e-014 s^10 + 2.422e-011 s^9 + 5.376e-009 s^8 + 6.188e-007 s^7+ 4.008e-005 s^6 + 0.001496 s^5 + 0.03256 s^4 + 0.4191 s^3+ 2.859 s^2 + 8.408 s 第4章线性控制系统的计算机辅助分析【1】(1) >> num=[1];den=[3 2 1 2];G=tf(num,den);eig(G)ans =-1.00000.1667 + 0.7993i0.1667 - 0.7993i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(2) >> num=[1];den=[6 3 2 1 1];G=tf(num,den);eig(G)ans =-0.4949 + 0.4356i-0.4949 - 0.4356i0.2449 + 0.5688i0.2449 - 0.5688i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(3) >> num=[1];den=[1 1 -3 -1 2];G=tf(num,den);eig(G)ans =-2.0000-1.00001.00001.0000分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(4) >> num=[3 1];den=[300 600 50 3 1];G=tf(num,den);eig(G)ans =-1.9152-0.14140.0283 + 0.1073i0.0283 - 0.1073i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(5) >> s=tf('s');G=0.2*(s+2)/(s*(s+0.5)*(s+0.8)*(s+3)+0.2*(s+2));eig(G)ans =-3.0121-1.0000-0.1440 + 0.3348i-0.1440 - 0.3348i分析：由以上信息可知，系统的所有极点都在s域的左半平面，因此系统是稳定的【2】(1) >> num=[-3 2];den=[1 -0.2 -0.25 0.05];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =0.5000 0.5000 0.2000分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(2) >> num=[3 -0.39 -0.09];den=[1 -1.7 1.04 0.268 0.024];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =1.1939 1.1939 0.1298 0.1298分析：由以上信息可知，由于前两个特征根的模均大于1，因此该系统是不稳定的(3) >> num=[1 3 -0.13];den=[1 1.352 0.4481 0.0153 -0.01109 -0.001043];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =0.8743 0.1520 0.2723 0.2344 0.1230分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(4) >> num=[2.12 11.76 15.91];den=[1 -7.368 -20.15 102.4 80.39 -340];H=tf(num,den,'Ts',0.5,'Variable','q');abs((eig(H))')ans =8.2349 3.2115 2.3415 2.3432 2.3432分析：由以上信息可知，所有特征根的模均大于1，因此该系统是不稳定的【3】(1) >>-4-3.5-3-2.5-2-1.5-1-0.50x 10-6P ole-Zero Map Real Axis (seconds -1)I m a g i n a r y A x i s (s e c o n d s -1)A=[-0.2,0.5,0,0,0;0,-0.5,1.6,0,0;0,0,-14.3,85.8,0;0,0,0,-33.3,100;0,0,0,0,-10]; eig(A) ans =-0.2000 -0.5000 -14.3000 -33.3000 -10.0000分析：由以上信息可知，该连续线性系统的A 矩阵的所有特征根的实部均为负数，因此该系统是稳定的(2)>>F=[17,24.54,1,8,15;23.54,5,7,14,16;4,6,13.75,20,22.5589;10.8689,1.2900,19.099,…21.896,3;11,18.0898,25,2.356,9];abs(eig(F)') ans =63.7207 23.5393 12.4366 13.3231 19.7275分析：由以上信息可知，该离散系统的F 矩阵的所有特征根的模均大于1，因此该系统是不稳定的【4】>> A=[-3 1 2 1;0 -4 -2 -1;1 2 -1 1;-1 -1 1 -2]; B=[1 0;0 2;0 3;1 1];C=[1 2 2 -1;2 1 -1 2];D=[0 0;0 0];G=ss(A,B,C,D); tzero(G)pzmap(G)ans =-3.6124-1.2765结论：∴可以得到该系统的 零点为-3.6124、-1.2765分析：由以上信息可知，系统的特征根的实部均位于s 域的左半平面，因此该系统是稳定的>> s=tf('s');G=0.2*(s+2)/(s*(s+0.5)*(s+0.8)*(s+3)+0.2*(s+2)); Gc=sscanform(G,'ctrl')Go=sscanform(G,'obsv')a =x1 x2 x3 x4x1 0 1 0 0x2 0 0 1 0x3 0 0 0 1x4 -0.4 -1.4 -4.3 -4.3b =u1x1 0x2 0x3 0x4 1c =x1 x2 x3 x4y1 0.4 0.2 0 0d =u1y1 0Continuous-time state-space model.a =x1 x2 x3 x4x1 0 0 0 -0.4x2 1 0 0 -1.4x3 0 1 0 -4.3x4 0 0 1 -4.3b =u1x1 0.4x2 0.2x3 0x4 0c =x1 x2 x3 x4y1 0 0 0 1d =u1y1 0Continuous-time state-space model.(1)>> num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320];[R1,P1,K1]=residue(num,[den 0]);[R1,P1]ans =-1.2032 -8.0000-1.0472 -7.00000.2000 -6.00000.7361 -5.0000-2.8889 -4.00002.2250 -3.0000-2.0222 -2.00003.0004 -1.00001.0000 0>> [n,d]=rat(R1);sym([n./d]')ans =[ -379/315, -377/360, 1/5, 53/72, -26/9, 89/40, -91/45, 7561/2520, 1][阶跃响应的解析解]y(t)=(-379/315)*e-8t+(-377/360)*e-7t+(1/5)*e-6t+(53/72)*e-5t+(-26/9)*e-4t+(89/40)*e-3t +(-90/45)*e-2t+(7561/2520)*e-t+1(2) >> num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320];[R2,P2,K2]=residue(num,den);[R2,P2]ans =9.6254 -8.00007.3306 -7.0000-1.2000 -6.0000-3.6806 -5.000011.5556 -4.0000-6.6750 -3.00004.0444 -2.0000-3.0004 -1.0000>> [n,d]=rat(R2);sym([n./d]')ans =[ 3032/315, 887/121, -6/5, -265/72, 104/9, -267/40, 182/45, -7561/2520][脉冲响应的解析解]y(t)=(3032/315)*e-8t+(887/121)*e-7t+(-6/5)*e-6t+(-265/72)*e-5t+(104/9)*e-4t+(-267/40) *e-3t+Linear Simulation ResultsA m p l i t u d e(182/45)*e -2t +(-7561/2520)*e -t(3) >> syms t;u=sin(3*t+5); Us=laplace(u) Us =(3*cos(5) + s*sin(5))/(s^2 + 9) >> s=tf('s');Us=(3*cos(5)+s*sin(5))/(s^2+9);num=[18 514 5982 36380 122664 222088 185760 40320]; den=[1 36 546 4536 22449 67284 118124 109584 40320]; G=tf(num,den); Y=Us*G; num=Y.num{1}; den=Y.den{1};[R3,P3,K3]=residue(num,den); [R3,P3] ans =1.1237 -8.0000 0.9559 -7.0000 -0.1761 -6.0000 -0.6111 -5.00002.1663 -4.0000 -1.1973 - 0.0010i 0.0000 +3.0000i -1.1973 + 0.0010i 0.0000 - 3.0000i -1.3824 -3.0000 0.8614 -2.0000 -0.5430 -1.0000 >> [n,d]=rat(R3); sym([n./d]') ans =[109/97, 282/295, -59/335, -965/1579, 951/439, - 449/375 + (18*i)/17981, - 449/375 - (18*i)/17981, -1663/1203, 317/368, -82/151] [正弦信号时域响应的解析解]y(t)=(109/97)*e -8t+(282/295)*e -7t+(-59/335)*e -6t+(-965/1579)*e -5t+(-449/375)*e -4t+(-1663/1203)*e -3t +(317/368)*e -2t +(-82/151)*e -t-2.3947sin(3t) [输出波形]>> num=[18 514 5982 36380 122664 222088 185760 40320]; den=[1 36 546 4536 22449 67284 118124 109584 40320];G=tf(num,den); t=[1:.1:20]';u=sin(3*t+5); lsim(G,u,t);分析：由解析解可知，输出信号的稳态 部分是振荡的，并且其幅值与相位始终 在到达稳态的时候保持不变，因此 右图所示的输出波形与解析解所得的结论是一致的【10】(1)因为PI 或PID 控制器均含有Ki/s 项，这是一个对误差信号的积分环节，假设去掉这一环节，则当Kp →∞，即|e(t)|很小也会存在较大扰动，这会影响到系统的动态特性；当加入这一环节后，如果要求|e(t)|→0，则控制器输出u(t)会由Ki/s 环节得到一个常值，此时系统可以获得较好的动态特性，因此这两个控制器可以消除闭环系统的阶跃响应的稳态误差(2)不稳定系统能用PI 或PID 控制器消除稳态误差。
计算机仿真技术薛定宇教材
[U,S,V]=svd(A) 例：a =
98
68
[u,s,v]=svd(a)
u=
0.7705 -0.6375
0.6375 0.7705
s=
15.5765 0
0
1.5408
2020/12/7
求矩阵A的奇异值及 分解矩阵，满足 U*S*V’=A，其中U、 V矩阵为正交矩阵 （U*U’=I），S矩阵 为对角矩阵，它的对 角元素即A矩阵的奇 异值。
>>a=[1 2 3;4 5 6]'
>>b=[1+2i 2-7i].'
a=
b=
14
1.0000 + 2.0000i
25
2.0000 - 7.0000i
36
>>b=[1+2i 2-7i]'
b = 1.0000 - 2.0000i
2020/12/7
2.0000 + 7.0000i23
2、四则运算与幂运算
单元数组：矩阵的每个元素不是数值，而是 可以存储任意类型的信息
类与对象：自行编写包含各种复杂详细信息 的变量，即类变量，变量可包含各种下级的 信息，也可重新定义运算
符号变量：区别与常规数值变量，便于公式 推导和解析求解
2020/12/7
7
补充：数值型数据的显示格式
任何MATLAB的语句的执行结果都可以 在屏幕上显示，同时赋值给指定的变量，没 有指定变量时，赋值给一个特殊的变量ans， 数据的显示格式由format命令控制。
2020/12/7
30
消除子块：如果将矩阵的子块赋值为空矩阵[ ]， 则相当于消除了相应的矩阵子块。
2、矩阵的大小
《控制系统仿真与CAD》课程教学大纲
《控制系统仿真与CAD》课程教学大纲一、课程基本信息二、理论教学内容及基本要求1 控制系统仿真与计算机辅助设计概述（1.5学时）1.1 为什么采用MATLAB语言？1.2 课程的主要结构1.3 控制系统计算机辅助设计方法概述2 MATLAB 必备基础知识（2.5学时）2.1 MATLAB 语言程序设计基础了解MATLAB 数据结构并掌握基本语句结构，包括直接赋值语句和函数调用语句。
了解并掌握循环语句结构、条件语句结构、开关语句结构以及与众不同的试探语句。
2.2 MATLAB 语言图形绘制方法掌握二维图形和三维图形绘制命令与函数与各种图形的绘制方法，可以将实验结果和仿真结果用可视的形式显示出来。
2.3 图形用户界面设计技术3 科学运算问题的MATLAB求解（4学时）掌握MATLAB 语言在线性代数问题、微分方程问题及最优化问题中的应用，系统地了解数学运算问题的MATLAB 求解方法。
4 线性控制系统的数学模型（7学时）4.1 线性连续系统的数学模型输入了解基本的数学模型，并掌握这些模型在MATLAB 环境下的表示方法及多变量系统的表示方法。
4.2 线性离散系统的数学模型输入了解基本的数学模型，并掌握这些模型在MATLAB 环境下的表示方法。
4.3 系统模型的相互转换掌握基于MATLAB 的系统模型转换方法，如连续与离散系统之间的相互转换和状态方程转换成传递函数模型的方法，以及转换成状态方程模型的各种实现方法。
4.4 方框图模型的化简了解三类典型的连接结构：串联、并联和反馈连接，并掌握模块输入、输出从一个节点移动到另一个节点所必需的等效变换及复杂系统的等效变换和化简。
4.5 线性系统的模型降阶了解线性系统的模型降阶的原理与实现。
4.6 线性系统的模型辨识了解线性系统的模型辨识的原理与实现。
5 线性控制系统的计算机辅助分析（5学时）5.1 线性系统定性分析了解系统稳定性的判定方法和系统的可控性和可观测性等系统性质的分析，并掌握其他的各种标准型实现。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
微积分计算：glfdiff()、fode_sol()
线性分数阶系统分析与设计
类的建立：@fotf；fotf.m、display.m 重载函数（FOTF互连）支持 *、+、feedback()
时域与频域分析：重载函数尽量控制系统工具箱 函数同名，且调用方式尽可能保持一致
step(), lsim(), bode(), nyquist(), nichols(), isstable(), norm()
国家级精品课程
控制系统仿真与CAD
第十一章 分数阶系统的分析与设计
东北大学信息学院 薛定宇
第十一章 分数阶系统的分析与设计
基于传统微积分理论的控制都是整数阶控制 误用的词：分数阶更确切的名称是非整数阶 本章主要内容
分数阶微积分的定义与计算 分数阶线性微分方程的求解 分数阶传递函数模型 —— FOTF 类
16/20
例11-13 分数阶微分方程 传递函数
输入信号
检验：用步长0.0005
控制系统仿真与CAD 国家级精品课程
2015/12/10
17/20
11.5 分数阶系统的设计
分数阶PID控制器
Igor Podlubny
分数阶对象最优分数阶PID控制器设计界面
控制系统仿真与CAD 国家级精品课程
微分、积分
扩展：统一微积分
控制系统仿真与CAD 国家级精品课程
2015/12/10
3/20
分数阶微积分的历史
法国数学家Guillaume Franç ois Antoine L’Hô pital
1695年，询问Leibniz， n=1/2？ Leibniz给出了t 的1/2导数 彼此关于此问题的通信标志 分数阶微积分研究的开始
2015/12/10
11/20
fotf.m 文件——创建类，@fotf 目录中
FOTF对象的输入
控制系统仿真与CAD 国家级精品课程
2015/12/10
12/20
display.m 模型显示 （略） 例11-6
重载函数 （置于@fotf目录）
系统连接：mtimes.m、plus.m、feedback.m等 时域分析：step.m、lsim.m 频域分析：bode.m、nyquist.m、nichols.m 稳定性、范数：isstable.m、norm.m
其他内容：解析解法、成比例阶系统、微分算子 近似，基于框图的分数阶非线性系统仿真方法
控制系统仿真与CAD 国家级精品课程
2015/12/10
20/20
控制系统仿真与CAD 国家级精品课程
2015/12/10
5/20
分数阶微积分的计算
给定函数采样值 y、时间 t、阶次 g MATLAB函数的直接调用
g 可以为正数也可以为负数
时间等间距
控制系统仿真与CAD 国家级精品课程
2015/12/10
6/20
分数阶微积分计算举例
函数调用 例11-1 常数的微积分是什么？
例11-10
控制系统仿真与CAD 国家级精品课程
2015/12/10
15/20
系统的范数
MATLAB 实现 例11-11 上例的系统范数
系统的频域分析：bode()、nyquist()、nichols() 例11-12 上例频域分析
控制系统仿真与CAD 国家级精品课程
2015/12/10
控制系统仿真与CAD 国家级精品课程
2015/12/10
13/20
例11-7 分数阶PID控制器
例11-8
例11-9 反馈系统
控制系统仿真与CAD 国家级精品课程
2015/12/10
14/20
11.3.2 分数阶线性系统分析
稳定性分析
基阶 稳定区域 重载函数编写 Nhomakorabea
isstable.m
例11-2 正弦函数的分数阶微分
信息量比整数阶微分丰富
控制系统仿真与CAD 国家级精品课程
2015/12/10
7/20
11.2 分数阶线性微分方程的求解
分数阶线性微分方程
闭式算法
MATLAB实现
控制系统仿真与CAD 国家级精品课程
2015/12/10
8/20
例11-4 分数阶微分方程求解
2015/12/10
18/20
例11-23 分数阶受控对象
输入模型
启动optimfopid界面
控制器 控制器 整数阶PID控制器 局限性：受控对象不能含有延迟
控制系统仿真与CAD 国家级精品课程
2015/12/10
19/20
分数阶控制系统小结
给出了分数阶微积分的定义与计算方法
输入信号单位阶跃 MATLAB求解
解的验证：选择步长为0.001
控制系统仿真与CAD 国家级精品课程
2015/12/10
9/20
11.3 分数阶传递函数模型与分析
Laplace变换在分数阶系统中应用
Laplace变换的重要性质
零初值：
分数阶传递函数
控制系统仿真与CAD 国家级精品课程
如何在MATLAB下定义并编程新的类、对象
线性分数阶系统的分析 分数阶系统控制器设计
控制系统仿真与CAD 国家级精品课程
2015/12/10
2/20
11.1 分数阶微积分的定义与计算
Newton 和 Leibniz 分别创立了微积分学
Sir Issac Newton的记号
微分
Gottfried Leibniz的记号
早期研究是纯数学研究 控制方面的研究
Manabe，1960开始 Igor Podlubny 1999
控制器、专著
控制系统仿真与CAD 国家级精品课程
2015/12/10
4/20
分数阶微积分的定义
Grü nwald-Letnikov定义
Riemann-Liouville定义
Caputo定义——非零初值
2015/12/10
10/20
11.3.1 分数阶传递函数类的编程
分数阶传递函数模型
定义一个类：FOTF类
建立 @fotf 目录 在此目录下至少建立两个文件
fotf.m：创建FOTF对象 display.m：显示该对象的内容
其他该对象的重载函数，置于@fotf 中
控制系统仿真与CAD 国家级精品课程