Matlab参考教程SIMULINK交互式仿真集成环境1.doc

9 SIMULINK交互式仿真集成环境
9.1 引导
SIMULINK是一个进行动态系统建模、仿真和综合分析的集成软件包。它可以处理的系统包括：线性、非线性系统；离散、连续及混合系统；单任务、多任务离散事件系统。
在SIMULINK 提供的图形用户界面GUI上，只要进行鼠标的简单拖拉操作就可构造出复杂的仿真模型。它外表以方块图形式呈现，且采用分层结构。从建模角度讲，这既适于自上而下（Top-down）的设计流程（概念、功能、系统、子系统、直至器件），又适于自下而上（Bottum-up） 逆程设计。从分析研究角度讲，这种SIMULINK模型不仅能让用户知道具体环节的动态细节，而且能让用户清晰地了解各器件、各子系统、各系统间的信息交换，掌握各部分之间的交互影响。
在SIMULINK环境中，用户将摆脱理论演绎时需做理想化假设的无奈，观察到现实世界中摩擦、风阻、齿隙、饱和、死区等非线性因素和各种随机因素对系统行为的影响。在SIMULINK环境中，用户可以在仿真进程中改变感兴趣的参数，实时地观察系统行为的变化。由于SIMULINK环境使用户摆脱了深奥数学推演的压力和烦琐编程的困扰，因此用户在此环境中会产生浓厚的探索兴趣，引发活跃的思维，感悟出新的真谛。
在MATLAB5.3版中，可直接在SIMULINK环境中运作的工具包很多，已覆盖通信、控制、信号处理、DSP、电力系统等诸多领域，所涉内容专业性极强。本书无意论述涉及工具包的专业内容，而只是集中阐述：SIMULINK 3.0的基本使用技法和相关的数值考虑。
节9.1虽是专为SIMULINK初学者写的，但即便是熟悉SIMULINK以前版本的读者也值得快速浏览这部分内容，因为新版的界面、菜单、工具条、模块库都有较大的变化。第9.2节比较详细地阐述建模的基本操作：通用模块的具体化设置、信号线勾画、标识、模型窗参数设置。这部分内容是进一步深入的前提。从第9.3节起，由浅入深地讲述SIMULINK对各种数学、工程问题的建模、仿真和分析的基本方法。
本章采用“算例”作为主体，配以适量的归纳性表述。本章包含了31个“尽量简单”又“独立完整”的“典型”算例，而这正是SIMULINK在线PDF文件之所缺。读者通过“手、眼、脑”并用地练习算例，掌握SIMULINK的一般使用规则和操作技法。
鉴于SIMULINK的本质，本节算例必定涉及数学、物理、和若干工程考虑。本书已采取“无量纲记述”、“注释”等措施使算例尽可能易读易懂，读者只要稍微耐心，就可以从这些有背景的内容体验到SIMULINK仿真之细腻和切实，从这些带背景性的算例品出SIMULINK的精妙之处。
9.1.1 SIMULINK的安装
9.1.2 SIMULINK入门
9.1.3 SIMULINK库浏览器界面
9.2 模型的创建
9.2.1 模型概念和文件操作
9.2.1.1 SIMULINK模型是什么
9.2.1.2 模型文件的操作
9.2.2 模块操作
9.2.2.1 模块的基本操作
9.2.2.2 向量化模块和标量扩展
【*例9.2.2.2-1】演示“示波”模块的向量显示能力。

图 9.2.2.2-1-1 示波器显示向量波形

【*例9.2.2.2-2】演示“求和”模块的向量处理能力：输入扩展。

图 9.2.2.2-2-1 输入的标量扩展

【*例9.2.2.2-3】演示“增益”模块的向量处理能力：参数扩展。

图 9.2.2.2-3-1 模块参数的标量扩展

9.2.2.3 参数设置
9.2.3 信号线操作
9.2.3.1 产生连线
9.2.3.2 信号线的分支和折曲
9.2.3.3 信号线宽度显示
9.2.3.4 彩色显示信号线
9.2.3.5 插入模块
9.2.3.6 信号线标识（label）
【例9.2.3.6-1】演示：信号线标识的传播

9.2.4 对模型的注释
9.2.5 常用的Sourse库信源
【*例9.2.5-1】如何调用MATLAB工作空间中的信号矩阵作为模型输入。本例所需的输入为 。
（1）编写一个产生信号矩阵的M函数文件
function TU=source925_1(T0,N0K)
t=linspace(0,K*T0,K*N0+1);
N=length(t);
u1=t(1:(N0+1)).^2;
u2=(t((N0+2):(2*N0+1))-2*T0).^2;
u3(1:(N-(2*N0+2)+1))=0;
u=[u1,u2,u3];
TU=[t',u'];

（2）构造简单的接受信号用的实验模型（如图9.2.5-1的左图）

图9.2.5-1 接受信号用的实验模型

（3）模块的参数设置
（4）在指令窗中，运行以下指令，在MATLAB工作空间中产生TU信号矩阵。
TU=source925_1(1,100,4);

（5）选中exm925_1模型窗菜单【Simulation：Start】，示波器呈现图9.2.5-1右图信号。

9.2.6 常用的Sink库信宿
9.2.6.1 库信宿一览表
9.2.6.2 示波器
9.2.7 仿真的配置
9.2.7.1 解算器参数的设置（Solver）
9.2.7.2 仿真数据的输入输出设置（Workspace I/O）
9.2.7.3 仿真中异常情况的诊断（Diagnostics）
9.3 连续系统建模
9.3.1 线性系统
9.3.1.1 积分模块的功用
【*例9.3.1.1-1】复位积分器的功用示例。

图9.3.1.1-1

9.3.1.2 积分模块直接构造微分方程求解模型
【*例9.3.1.2-1】假设从实际自然界（力学、电学、生态等）或社会中，抽象出有初始状态为0的二阶微分方程，是单位阶跃函数。本例演示如何用积分器直接构搭求解该微分方程的模型。
（1）改写微分方程
（2）利用SIMULINK库中的标准模块构作模型

图9.3.1.2-1-1 求解微分方程的SIMULINK模型exm9312_1.mdl

（3）仿真操作
（4）保存在MATLAB工作空间中的数据
clf
tt=ScopeData.time; %为书写简单，把构架域的时间数据另赋给tt
xx=ScopeData.signals.values; %目的同上。
[xm,km]=max(xx);
plot(tt,xx,'r','LineWidth',4),hold on
plot(tt(km),xm,'b.','MarkerSize',36),hold off
strmax=char('最大值',['t = ',num2str(tt(km))],['x = ',num2str(xm)]);
text(6.5,xm,strmax),xlabel('t'),ylabel('x')

图9.3.1.2-1-2 利用存放在MATLAB工作空间中的仿真数据所绘制的曲线

93.1.3 传递函数模块
【例9.3.1.3-1】直接利用传递函数模块求解方程(9.3.1.3-1)。
（1）根据式(9.3.1.3-3)构造如图9.3.1.3-1所示的模型exm9313_1.mdl

图9.3.1.3-1 由传递函数模块构成的仿真模型exm9313_1.mdl

（2）仿真操作

9.3.1.4 状态方程模块和单位脉冲输入的生成
【*例9.3.1.4-1】假设式(9.3.1.4-1)中的输入函数是单位脉冲函数，研究该系统的位移变化。本例演示：（A）状态方程模块的使用；（B）脉冲函数的生成方法。
（1）单位脉冲函数的数学含义及近似实现

（2）利用库模块构造如图9.3.1.4-1-1所示的仿真模型exm9314_1.mdl

图9.3.1.4-1-1 带近似单位脉冲的状态方程模块构成的仿真模型exm9314_1.mdl

（3）仿真结果

图9.3.1.4-1-2 仿真结果

9.3.2 非线性系统
9.3.2.1 建立非线性仿真模型的基本考虑
【*例9.3.2.1-1】物理背景：如图9.3.2.1-1-1所示喷射动力车的定位控制问题。

图9.3.2.1-1-1 装置左右喷射发动机的车辆示意图

（1）根据车辆的动态方程，构作基本仿真模型

图9.3.2.1-1-2 基本仿真模型

（2）为观察仿真时间进程，引入仿真时钟显示。
（3）为模仿“车辆速度与位移小于某阈值时被认为控制目标达到”，引入仿真终止环节。

图9.3.2.1-1-3 完整仿真模型exm9321_1.mdl

（4）为模仿“符号切换的非瞬时性”，本例采用“定步长”Solver解算器
（5）为便于分析接近原点时的颤振现象，通过模型窗I/O设置保存仿真数据
（6）仿真操作和仿真结果观察
（7）保存数据的利用和分析
subplot(1,2,1),plot(xout(:,2),xout(:,1))
grid on,axis([-0.2,1,-1,0.2]),axis square
xlabel('\fontsize{14}位移'),ylabel('\fontsize{14}速度'),
subplot(1,2,2),plot(xout(:,2),xout(:,1))
grid on,axis([-0.1,0.05,-0.05,0.1]),axis square

图9.3.2.1-1-4 完整相轨迹及局部放大图

9.3.2.2 任意非线性函数模块及其应用
【例9.3.2.2-1】轿车沿直线山坡路向前行驶。要求设计一个简单的比例放大器，使轿车能以指定的速度运动。本例演示：（A）仿真系统的创建。（B）非线性模块的使用。（C）任意函数模块的应用D)体现“自下而上”的建模方式。（E）本例将作为下面章节多个算例的基础，读者切莫跳略此题。
（1）轿车的运动方程
（2）轿车动力学的仿真模型

图9.3.2.2-1-2 行驶轿车动力学仿真模型

（2）比例控制及其仿真模型

图9.3.2.2-1-3 比例控制器仿真模型

（3）构成完整的仿真模型

图9.3.2.2-1-4 受控轿车的完整仿真模型

（4）仿真操作和说明

9.4 子系统的创建、装帧及受控执行
9.4.1 简装子系统及其应用
9.4.1.1 创建简装子系统的“先有内容后套包装”法
【例9.4.1.1-1】题目的背景和参数与例9.3.2.2-1完全相同，要求创建利用比例控制器使轿车的运动速度稳定在期望车速的分层仿真模型。本例演示：如何从非分层模型获得分层模型；创建简装子系统的“先有内容后套包装”法。
（1）将待“包装”模型另存为一个新模型
（2）生成第一个简装子系统
（3）生成第二个简装子系统
（4）更改子系统的标识名
（5）重命名输入输出口的名称

图9.4.1.1-1 采用简装子系统的完整仿真模型

9.4.1.2 创建简装子系统的“先有包装后置内容”法
【例9.4.1.2-1】本例演示：如何自上而下构造分层模型；产生简装子系统的“先有包装后置内容”法。
（1）在新建模型窗中复制包括子系统模块在内的所有标准模块（见图9.4.1.2-1）

图 9.4.1.2-1 带标准子系统模块的新建模型窗

（2）在Subsystem标准模块中“装入”轿车动态模型
（3）在Subsystem 1标准模块中“装入”比例控制器模型
（4）对Untitled模型窗中滑键模块进行适当设置，进行各模块间的连接，就形成与图9.4.1.1-1相同的完整仿真模型。
（5）对完整仿真模型窗进行仿真参数设置，为进行仿真作最后准备

9.4.2 精装子系统
9.4.2.1 精装子系统的制作过程
9.4.2.2 装帧示例
【*例9.4.2.2-1】目标：把图94.1.1-1所示轿车速度控制模型中的轿车动态模型简装子系统变成精装子系统。

（2）打开装帧编辑器

图9.4.2.2-1-1 显示初始化页面的装帧编辑器
（3）填写初始化页面中“块类别”
（4）使用初始化页面中的“对话窗参数设置区”为未来对话框设计提示及设置变量
（5）填写初始化页中“初始指令区”
（6）利用Icon页为精装子系统模块制作图标

图 9.4.2.2-1-2 装帧编辑器的图标制作页

（7）利用文档页为精装子系统模块编写模块性质描述和在线帮助说明

图 9.4.2.2-1-3 装帧编辑器的文档制作页
（8）点击装帧编辑器左下角的【OK】键整个装帧便告成功。

图9.4.2.2-1-4 制作完成的简装子系统

9.4.2.3 精装子系统的使用特点
【例9.4.2.3-1】本例演示：精装子系统参数对话窗的来源和外形特点；如何打开精装子系统自身的“下层”结构模型；精装子系统如何从外界获得参数。
（1）精装子系统参数对话窗的来源和工作特点

图9.4.2.3-1-1 精装子系统的参数设置对话窗

（2）利用精装子系统所在窗口的【Look under Mask】菜单项，可打开精装子系统的结构模型窗。
（3）参数值的传递路径

9.4.3 条件执行子系统
9.4.3.1 使能子系统
【*例94.3.1】利用使能原理构成一个半波整流器。本例演示使能子系统的创建及工作机理。
（1）打开SIMULINK的新建模型窗口
（2）从SIMULINK库中提取三个模块Sine wave、 Subsystem 、 Scope 到新建窗。然后进行文件保存操作，并起文件名为exm9431_1（这保存操作只为以后调用方便，并非必要）。
（3）双击空子系统模块Subsystem ，打开其结构模型窗。
（4）从SIMULINK库中拷贝In输入口模块、Out输出口模块、Enable使能模块到子系统的结构模型窗；把In 模块的输出直接送到Out模块的输入端；Enable模块无须进行任何连接，且本例采用它的缺省设置；便实现了题目所需使能子系统（图9.4.3.1-1 b）。
（5）按图9.4.3.1-1-1 a 那样，完成exm9431_1窗口中各模块间的连接。
（6）双击示波器模块，打开显示窗。然后选择exm9431_1窗口菜单项【Simulation:Start】，就可看到如图9.4.3.1-1-2的半波整流后的波形。

图9.4.3.1-1-1 利用使能子系统实现半波整流的仿真模型

图9.4.3.1-1-2 所得仿真波形

9.4.3.2 触发子系统
【*例9.4.3.2-1】利用触发子系统获得零阶保持的采样信号。本例演示：触发子系统工作原理；在MATLAB指令窗中运行SIMULINK模型。
（1）构造如图9.4.3.2-1-1所示的仿真模型exm9432_1.mdl

图9.4.3.2-1-1 利用触发子系统实现零保持采样的仿真模型exm9432_1

（2）在MATLAB指令窗中运行SIMULINK模型
[t,x,y]=sim('exm9432_1',10);
clf,hold on
plot(t,y(:,1),'b')
stairs(t,y(:,2),'r')
stairs(t,y(:,3),'c:'),hold off
axis([0 10 -1.1 1.1]),box on
legend('sinewave','output','trigger',4)

图 9.4.3.2-1-2 exm9432_1.mdl模型仿真运行结果

9.4.3.3 触发使能子系统
9.4.3.4 使能子系统和出发子系统综合运用示例
【*例9.4.3.4-1】本例是前面例9.3.2.2-1 , 9.4.1.1-1 , 9.4.2.2-1的继续，使得汽车速度受两种不同的控制器操纵。具体要求是：（A）当汽车实际速度与期望速度的误差绝对值，且时，将切换为PI比例-积分控制器；（B）一旦PI控制器被使用，只要仍满足，那么PI将继续起控制作用；（C）除以上情况外，则都使用简单的P比例控制器
（1）建立P Controller使能子系统和PI Controller使能子系统

图 9.4.3.4-1-1 P和PI控制器使能子系统的结构图

（2）建立Model Selector子系统模块

图 9.4.3.4-1-2 Model Selector子系统的结构图

（3）为保存“控制器切换时间”创建触发子系统

图 9.4.3.4-1-3 SaveSwitchTime触发子系统的结构图

（4）建立完整的仿真模型

图 9.4.3.4-1-4 轿车速度双模式控制系统的仿真模型

（5）仿真结果

图 9.4.3.4-1-5 两个仿真示波器上的曲线记录

9.4.3.5 交替执行子系统
【*例9.4.3.5-1】在例9.4.3.4-1中，比例控制器和比例-积分控制器的工作切换是借助Model Selector子系统产生的两个输出切换信号Choose PI和Choose P实现的。本例将演示：如何依靠一个Choose PI信号和merge汇合模块的配合使用，实现同样的控制器切换。
（1）把exm9434_1.mdl另存为exm9435_1.mdl 。
（2）把图9.4.3.4-1-4中虚线框所围的那部分改成图9.4.3.5-1所示那样。

图9.4.3.5-1 等效于图9.4.3.4-1-5虚线框所围部分的模型

9.5 离散时间系统和混合系统
9.5.1 若干基本模块
【*例9.5.1-1】用组合逻辑模块产生的“逻辑和”结果及“逻辑或”结果。
（1）建立输入输出关系（见表9.5.1-2。真值表用黑框所围）
（2）建立模型exm951_1.mdl（见图9.5.1-1）

图9.5.1-1 逻辑运算模块和组合逻辑模块运行对照模型

9.5.2 多速率离散时间系统
【*例9.5.2-1】在离散控制系统中，控制器的更新频率一般低于对象本身的工作频率。而显示系统的更新频率总比显示器的可读速度低得多。假设有某过程的离散状态方程

式中是输入。该过程的采样周期为秒。控制器应用采样周期为秒的比例控制器；显示系统的更新周期为秒。
（1） 建立模型exm952_1.mdl（见图9.5.2-1-1）

图9.5.2-1-1 经着色的多速率离散系统

（2）在exm952_1模型窗中，选取菜单项【Format:Sample time colors】后，模型中不同采样周期的模块和连线就会以不同颜色表示。在本例中，采样速度最快的被控过程部分显示为红色；速度次之的控制器部分显示为绿色；显示x1历史记录的那部分显示为兰色。

（3）仿真结果
tt=TX.time;
x1=TX.signals.values;
plot(tt,x1)grid on,
xlabel('kT'),ylabel('x1(kT)')

图9.5.2-1-1 的历史记录

9.5.3 离散-连续混合系统
【例9.5.3-1】本例是在例9.4.2.2-1的基础上进行的。目标是：设计一个离散PID控制器子系统对轿车速度进行控制。本例演示：（A）离散PID的构成；（B）展示仿真模型在研究控制器各参数影响上的能力。
（1）打开例9.4.2.2-1中模型exm9422_1.mdl ，选用菜单项【File:Save as】产生“另存”模型exm953_1.mdl 。它将用做以下改建的“原始”模型。
（2）把exm953_1窗口中的P Controller改名为PID Controller。（参见图9.5.3-1-1）

图9.5.3-1-1 采用离散PID的轿车速度控制系统

（3）双击PID Controller子系统“原始”模块，打开该子系统的exm953_1/PID Controller结构模型窗。
（4）把exm953_1/PID Controller模型窗中模型改建成图9.5.3-1所示。

图9.5.3-1-2 exm953_1/PID Controller子系统的结构模型

（5）在exm953_1模型窗中，选取菜单项【Format:Sample time colors】，使该模型及其所有子系统模型按采样速率着色：连续时间部分用黑色；离散时间部分用红色；离散、连续混合的子系统被着黄色。

（6）仿真实验

图9.5.3-1-3 exm953_1/PID Controller子系统的结构模型

9.6 SIMULINK的分析工具
9.6.1 确定模型的特征
【例9.6.1-1】观察例9.5.3-1中所建模型exm953_1.mdl中状态向量的结构。
在指令窗或Notebook中运行以下指令。
[sizes,x0,StateCell]=exm953_1;
SIZES=sizes',X0=x0',StateCell
SIZES =

2 2 0 0 0 0 3
X0 =
0 0 0 0
StateCell =
'exm953_1/Automobile Model/Int1'
'exm953_1/Automobile Model/Int2'
'exm953_1/PID Controller/DD'
'exm953_1/PID Controller/DT-I'

9.6.2 用MATLAB指令运行SIMULINK模型
9.6.2.1 运行SIMULINK模型的sim指令
9.6.2.2 设置编辑仿真参数的simset 指令
9.6.2.3 获取模型仿真参数的simget指令
9.6.2.4 MATLAB指令运行SIMULINK模型的示例
【*例9.6.2.4-1】以例9.5.3-1中所建模型exm953_1.mdl为基础进行本题解算。演示：（A）显示模型窗中的初始状态设置。（B）把初始车速重置为120 ，而其他初始值仍为0 。（B）画出两种初始状态下的车速曲线。
InInit=simget('exm953_1','InitialState') %获取模型窗对初始值的设置
[t,x,y]=sim('exm953_1',100); %在模型内设置参数下进行仿真
opts=simset('InitialState',[120,0,0,0]); %初始值的重置
[tt,xx,yy]=sim('exm953_1',100,opts); %在重置初值下仿真
plot(t,x(:,1),':b',tt,xx(:,1),'r')
legend('\fontname{隶书}\fontsize{16}内初值','外初值',4)
InInit =
[]

图9.6.2.4-1 两种初值下的仿真曲线

9.6.3 模型的线性化问题
9.6.3.1 线性化的数学描述
9.6.3.2 连续系统的线性化模型
9.6.3.3 离散系统的线性化模型
9.6.3.4 模型线性化的算例
【*例9.6.3.4-1】求非线性系统在坐标原点处的线性化模型。
（1）根据式（9.6.3.1-2）的可“手工”求得。
（2）创建图9.6.3.4-1所示的SIMULINK模型exm9634_1.mdl 。

图9.6.3.4-1 SIMULINK模型exm9634_1.mdl

（3）用指令求坐标原点处的线性化模型（结果与手工计算一致）
[A,B,C,D]=linmod('exm9634_1');A
A =
0.0000 0.0000
2.0000 -1.0000

（4）用指令求[ 1 , 0.5]坐标处的线性化模型（注意模型与前一个不同）
[A1,B1,C1,D1]=linmod('exm9634_1',[1,0.5]);A1
A1 =
2.0000 1.0000
2.0000 -1.0000

（5）由特征根看模型稳定性
eA=eig(A)',eA1=eig(A1)'
eA =
0.0000 -1.0000
eA1 =
2.5616 -1.5616

9.6.4 系统平衡点的求取
9.6.5 综合算例
9.6.5.1 “一步仿真”和精良状态轨迹斜率图
【*例9.6.5.1-1】求非线性系统的相平面轨迹、平衡点，并进行稳定性分析。本例综合演示：（A）SIMULINK模型和MATLAB指令的配合使用。（B）sim , simset , trim 指令的应用。（C）“一步仿真”计算方法。（D）二阶系统相轨迹的精良图形。
（1）非线性系统由SIMULINK模型exm9634_1.mdl表达（见例9.6.3.4-1）
（2）编写绘制传统状态轨迹（State trajectory）的M文件exm9651_1.m
[exm9651_1.m]
% exm9651_1.m 画普通状态轨迹
clf;hold on
xx=[2,1;-1,1;0,1;1,1;1,0;1,-1;1,-2]; %轨线起点
nxx=size(xx,1); %起点数
for k=1:nxx
opts=simset('initialstate',[xx(k,1),xx(k,2)]); %设置仿真初值
[t,x,y]=sim('exm9634_1',10,opts);
plot(x(:,1),x(:,2)); %画状态轨线
end
xlabel('x1');ylabel('x2'),grid,hold off

（3）运行exm9651_1.m 得如图9.6.5.1-1-1所示的传统状态轨迹图

图 9.6.5.1-1-1 系统exm9634_1.mdl的状态轨迹图

（4）编写采用“一步仿真”计算状态轨线斜率的函数M文件portraitzzy.m
[portraitzzy.m]
function [DX1,DX2,DP]=portraitzzy(x1,x2,h)
% PORTRAITZZY 采用“一步仿真”计算状态变量斜率和状态导数的二次方根
% x1,x2 分别给定“状态平面”上的格点坐标
% h 给定积分计算采用的时间步长
% DX1,DX2 轨线斜率在状态坐标轴上的投影长度
% DP 状态导数向量的长度
opts=simset('solver','ode5','fixedstep',h); %采用ode5定步长积分算法 <7>
n=length(x1);
X1=zeros(n,n);X2=X1; %预置空间
for ii=1:n;
for jj=1:n;
opts=simset(opts,'initialstate',[x1(ii),x2(jj)]);%设置状态初值 <12>
[t,x,y]=sim('exm9634_1',h,opts); %步长为 h 的“一步仿真” <13>
dx1=x(2,1)-x1(ii); %计算x1的变化率
dx2=x(2,2)-x2(jj); %计算x2的变化率
L=sqrt(dx1^2+dx2^2); %计算状态轨线长度的变化率。
Z(jj,ii)=L;
if L>1.e-10 %若状态轨线变化率大于“零”阈值
DX1(jj,ii)=dx1/L;DX2(jj,ii)=dx2/L; %计算各状态变量的近似斜率 <19>
%注意下标次序。这是绘图指令格式要求
end
end
end
DP=Z/h; %状态导数向量的长度

（5）绘制精良的状态轨迹斜率图（Quiver plot）
h=0.01; %设置仿真步长
x1=-2.5:0.25:2.5;x2=x1; %轨线起始点
k=0.15;
[X1,X2,Z]=portraitzzy(x1,x2,h);
quiver(x1,x2,k*X1,k*X2,0) %调用quiver指令绘制平面上各点处的变化率图
xlabel('x1'),ylabel('x2')

图 9.6.5.1-1-2 系统exm9634_1.mdl的精良状态轨迹斜率图

6）运行以下指令绘制状态导数向量的长度分布曲面（见图9.6.5.1-1-3）
surfc(x1,x2,Z),view([18,32]),xlabel('x1'),ylabel('x2')

图 9.6.5.1-1-3 系统exm9634_1.mdl的状态导数向量长度分布曲面

（6）利用trim指令求系统平衡点
xa=trim('exm9634_1',[-1,-2]')
xb=trim('exm9634_1',[1,2]')
xa =
-0.8944
-1.7889
xb =
0.8944
1.7889

（7）运用线性化模型进行系统稳定性分析
Axa=linmod('exm9634_1',xa);eig_Axa=(eig(Axa))'
Axb=linmod('exm9634_1',xb);eig_Axb=(eig(Axb))'
eig_Axa =
-1.3944 - 2.6457i -1.3944 + 2.6457i
eig_Axb =
3.4110 -2.6222

9.6.5.2 仿真模型和优化指令的协调
【*例9.6.5.2-1】题目背景：在迄今的自动控制教材中，凡讨论积分性能指标时，几乎总会提到所谓的ITAE传递函数标准型，并列出相应的分母多项式系数表。但值得指出的是：这些数据是20 世纪50年代初期，用模拟计算机仿真得到的。因此，这些数据的准确性带有明显的时代缺陷。与不同，ITAE性能函数无法解析计算，而只能通过数值计算进行。本例将综合运用SIMULINK模型和优化指令求取新的ITAE传递函数标准型。
本例演示：（A）如何用SIMULINK模块计算性能函数。（B）SIMULINK方块模型、目标函数和优化程序之间的协调和参数传递。（C）跨空间交换数据、跨空间计算表达式。
（1）问题的形成
（2）构作计算的SIMULINK模型（如图9.6.5.2-1-1所示）

图9.6.5.2-1-1 计算的SIMULINK模型smodel.mdl

（3）编写实现寻优的M文件exm9652_1.m和itae.m
[itae.m]
function ss=itae(aa)
%ITAE 计算ITAE积分性能值
global a ss
a=aa;
Tspan=evalin('base','Tspan'); %从主空间获得采样时间向量 <5>
opts=simset('RelTol',0.0001); %设定仿真时的相对精度
[tt,x,s]=sim('smodel',Tspan,opts); %对SIMULINK模型smodel仿真 <7>
ss=s(end); %优化过程中产生的ITAE值。 <8>

[exm9652_1.m]
%exm9652_1.m 用于ITAE标准型系数的寻优
clear
global a ss
a0=[3.3 6.6 8.6 7.5 3.9]; %被优化系数的初值 <4>
Tspan=(0:500)/10; %计算ITAE值用的时间采样向量 <5>
options.TolFun=0.0001; %设定优化目标函数的精度
options.TolX=0.0001; %设定被优化参数的精度
a=fminunc('itae',a0,options); %用拟牛顿法寻优 <8>
coeff=[1,a,1] %输出完整的标准型系数 <9>
ss %输出相应的ITAE性能值

（4）在指令窗中运行以下指令并得到结果
exm9652_1
coeff =
1.0000 2.1516 5.6284 6.9323 6.7917 3.7394 1.0000
ss =
8.3338

（5）ITAE标准传递函数“经典”系数和修正系数的阶跃响应比较
old=tf(1,[1 3.25 6.60 8.60 7.45 3.95 1]);
new=tf(1,coeff);
[yold,told]=step(old,50);
[ynew,tnew]=step(new,50);
plot(told,yold,'b','LineWidth',1)
axis([3,18,0.95,1.05])
hold on,plot(tnew,ynew,'r','LineWidth',3),hold off
legend('Old','New',4),grid on

图9.6.5.2-1-1 新老标准型的阶跃响应局部放大比较图

9.7 数值计算方面的考虑
9.7.1 微分方程解算器Solver
9.7.1.1 ODE45和ODE23运作机理简要
9.7.1.2 ODE113运作机理简要
9.7.1.3 ODE15S和ODE23S运作机理简要
9.7.1.4 不同解算器处理Stiff系统时表现
【*例9.7.1.4-1】求微分方程在时的解。本例演示：对于Stiff方程，如果解算方法选择不当将产生严重后果。
（1）用符号计算求此方程解析解和它的导函数
xsym=dsolve('D2x+100*Dx+0.9999*x=0','x(0)=1,Dx(0)=0','t')
dxsym=diff(xsym,'t')
xsym =
9999/9998*exp(-1/100*t)-1/9998*exp(-9999/100*t)
dxsym =
-9999/999800*exp(-1/100*t)+9999/999800*exp(-9999/100*t)

（2）构作SIMULINK模型exm9714_1.mdl （如图9.7.1.4-1-1所示）

图9.7.1.4-1-1 Stiff系统仿真模型

（3）运行以下指令画出三种计算结果的局部放大图
t=(0:5000)/10;x=eval(char(dxsym)); %产生绘解析解图形的数据
Tspan=500;
opts=simset('Solver','ode45');
[tt1,xx1,s]=sim