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 。

（C）画出两种初始状态下的车速曲线。
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.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.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 的状态导数向量长度分布曲面

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

（8）运用线性化模型进行系统稳定性分析
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); % 设定仿真时的相对精度

[1] [2] 下一页