Simulink交互式仿真6.docx
《Simulink交互式仿真6.docx》由会员分享，可在线阅读，更多相关《Simulink交互式仿真6.docx（26页珍藏版）》请在冰点文库上搜索。
Simulink交互式仿真6
.1.1系统平衡点和普通状态轨线图
【例8.7-7】求图8.7-13所示两个模型的平衡点。
模型（b）输入端比模型（a）多一个输入口。
exm080707_1.mdl
exm080707_2.mdl
（a）
（b）
图8.7-13待求平衡点的非线性系统块图模型
（1）
（2）
xa=trim（'exm080707_1',[-0.1;-0.3]）%<1>
xb=trim（'exm080707_1',[0;1]）%<2>
xa=
-0.8944
-1.7889
xb=
0.8944
1.7889
（4）
Axa=linmod2（'exm080707_1',xa）;
eig_Axa=（eig（Axa））'
Axb=linmod2（'exm080707_1',xb）;
eig_Axb=（eig（Axb））'
eig_Axa=
-1.3944-2.6457i-1.3944+2.6457i
eig_Axb=
3.4110-2.6222
（6）
[xa2,ua]=trim（'exm080707_2',[-0.1;-0.3],0）%<3>
[xb2,ub]=trim（'exm080707_2',[0;1],1）%<4>
xa2=
-0.7487
-1.4974
ua=
1.1974
xb2=
0.6810
1.3620
ub=
1.6810
（2）
%exm080707_1m.m
clf;
xx=[-2,-1,0,1,1,1,1
1,1,1,1,0,-1,-2];
set_param（'exm080707_1','InitInArrayFormatMsg','None'）
nxx=size（xx,2）;
fork=1:
nxx
opts=simset（'initialstate',xx（:
k））;%<5>
[t,x]=sim（'exm080707_1',10,opts）;%<6>
plot（x（:
1）,x（:
2））;
holdon
end
gridon,holdoff
xlabel（'x1'）;ylabel（'x2'）
title（'普通状态轨线'）
图8.7-14多初始点出发的状态轨线和平衡点
.1.2M码和Simulink模型的综合运用
101单步仿真和精良状态轨线图
【例8.7-8】绘制非线性系统
块图模型的精良状态变化轨线。
（1）
function[DX1,DX2,DP]=exm080708_zzy（x1,x2,h）
opts=simset（'solver','ode5','fixedstep',h）;%<2>
n=length（x1）;
DX1=zeros（n,n）;DX2=DX1;DP=DX1;
disp（'正在逐点计算，请稍等！
'）
forii=1:
n;
forjj=1:
n;
opts=simset（opts,'initialstate',[x1（ii）;x2（jj）]）;%<8>
[~,x]=sim（'exm080707_1',h,opts）;%<9>
dx1=x（2,1）-x1（ii）;
dx2=x（2,2）-x2（jj）;
L=sqrt（dx1^2+dx2^2）;
DP（jj,ii）=L/h;
ifL>1.e-10
DX1（jj,ii）=dx1;DX2（jj,ii）=dx2;%<15>
end
end
end
disp（'计算结束'）
（2）
%exm080708m.m
h=0.01;
x1=（-2.5:
0.25:
2.5）';x2=x1;
k=3.5;
set_param（'exm080707_1','InitInArrayFormatMsg','None'）
xs=trim（'exm080707_1',[-0.1;-0.3]）;
xus=trim（'exm080707_1',[0;1]）;
[DX1,DX2,DL]=exm080708_zzy（x1,x2,h）;
pcolor（x1,x2,DL）
shadinginterp
alpha（0.5）
colorbar
holdon
quiver（x1,x2,k*DX1,k*DX2,0）
plot（xs
（1）,xs
（2）,'bo',xs
（1）,xs
（2）,'+','MarkerSize',10）
plot（xus
（1）,xus
（2）,'bo',xus
（1）,xus
（2）,'.','MarkerSize',10）
gridoff
holdoff
xlabel（'x1'）,ylabel（'x2'）
title（'精良状态轨线斜率图'）
shg
图8.7-15精良状态轨线迹斜率图
102仿真模型和优化指令的协调
【例8.7-9】题目背景：
在迄今的自动控制教材中，凡讨论积分性能指标时，几乎总会提到所谓的ITAE传递函数标准型，并列出相应的分母多项式系数表。
但值得指出的是：
这些数据是20世纪50年代初期，用模拟计算机仿真得到的。
因此，这些数据的准确性带有明显的时代缺陷。
与
不同，ITAE性能函数
无法解析计算，而只能通过数值计算进行。
图8.7-16计算
的块图模型
（1）问题的形成
（2）
（3）
%exm080709m.m
globalaJc
amin=min（a0）;
na=length（a0）;
nd=na+2;
opts=optimset（'MaxFunEvals',300*na）;
CF=zeros（Kr,nd）;Jk=zeros（1,Kr）;
forkk=1:
Kr
ar=a0+2*amin*（rand（1,na）-0.5）;%<8>
a=fminsearch（@exm080709_itae,ar,opts）;
cf=[1,a,1];
CF（kk,:
）=cf;
Jk（kk）=Jc;
end
[Jmin,kmin]=min（Jk）;
cfmin=CF（kmin,:
）;
%exm080709_itae.m
functionJc=exm080709_itae（aa）
globalaJc
a=aa;
Tspan=[0,0.1,20];%<4>
opts=simset（'RelTol',0.0001）;
[~,~,Jt]=sim（'exm080709',Tspan,opts）;%<6>
Jc=Jt（end）;
（4）
clear
Kr=5;
a0=[3.25,6.60,8.60,7.45,3.95];
exm080709m
Jmin,cfmin
Jmin=
8.3338
cfmin=
Columns1through6
1.00002.15195.62906.93386.79253.7398
Column7
1.0000
（5）
old=tf（1,[1,a0,1]）;
new=tf（1,cfmin）;
[yold,told]=step（old,50）;
[ynew,tnew]=step（new,50）;
plot（told,yold,'b','LineWidth',1）
axis（[0,18,0,1.1]）
holdon,plot（tnew,ynew,'r','LineWidth',3）,holdoff
xlabel（'t'）
title（'ITAE6阶新老标准型的阶跃响应比较'）
legend（'Old','New',4）,gridon
图8.7-17新老标准型的阶跃响应局部放大比较图
表8.7-2ITAE标准型新系数（黑体）和老“经典”系数（细体）对照
阶次
ITAE值
传递函数分母多项式系数
2
1.99
1.9519
11.41
11.50491
3
3.144
3.1383
11.752.151
11.78282.17151
4
4.626
4.5913
12.103.402.751
11.95213.34582.64731
5
7.155
6.3215
12.805.005.503.401
12.06674.49764.67303.25681
6
9.656
8.3338
13.256.608.607.453.951
12.15195.62906.93386.79253.73981
7
15.003
10.6290
14.4810.4215.0515.5410.644.5801
12.21696.74339.346911.5778.67784.32261
8
18.680
13.2051
15.2012.8021.6025.7522.2013.305.151
12.26817.831311.847217.532516.064511.30944.80691
●研究表明：
ITAE函数搜索空间的形状非常复杂，凹凸不平，小谷很多，许多地方深谷高峰相邻。
要找到真正最小值点决非易事。
虽可以肯定：
单点标准型的新系数比老系数具有更小的ITAE值；但不能断言这新系数一定指示着最小值点。
.2数值计算方面的考虑
.2.1微分方程解算器Solver
101ode45和ode23运作机理简要
102ode113运作机理简要
103ode15s和ode23s运作机理简要
104不同解算器解Stiff方程的表现
【例8.8-1】求微分方程
在
时的解。
图8.8-1微分方程的块图模型exm080801
（1）关于exm080901.mdl的说明
（2）
symstxxd
xs=dsolve（'D2x+100*Dx+0.9999*x=0','x（0）=1,Dx（0）=0','t'）
xsd=diff（xs,'t'）
HL2=ezplot（xd-xsd,[0,10,-0.012,0]）;
set（HL2,'LineWidth',3）
title（['x''=',char（xsd）]）
xs=
9999/（9998*exp（t/100））-1/（9998*exp（（9999*t）/100））
xsd=
9999/（999800*exp（（9999*t）/100））-9999/（999800*exp（t/100））
图8.8-2微分方程的解x和它的导数dx/dt
（3）
tt=（0:
4000）/10;
xx0=subs（xsd,t,tt）;
Tspan=600;
opts=simset（'Solver','ode45'）;
[tt1,xx1,s]=sim（'exm080801',Tspan,opts）;
opts=simset（'Solver','ode15s'）;
[tt2,xx2,s]=sim（'exm080801',Tspan,opts）;
plot（tt,xx0,'k',tt1,xx1（:
2）,'b:
',tt2,xx2（:
2）,'r-.'）
axis（[246247-8.55e-4-8.35e-4]）
legend（'Symbolic','ode45','ode15s',0）
xlabel（'t'）,ylabel（'dx/dt'）
title（'Stiff方程的三种算法结果比较局部放大'）
ns1=length（xx1）
ns2=length（xx2）
ns1=
18085
ns2=
101
图8.8-3不同方法的解算结果比较
.2.2积分步长和容差
101积分步长的选择
102计算容差的选择
.2.3代数环问题
101无惯性模块和代数环
102消减代数环影响
【例8.8-2】构建由方程（8.8-2）和（8.8-3）表述系统的Simulink块图模型，讨论代数环。
（8.8-2）
（8.8-3）
图8.8-4带隐式代数方程的块图模型exm080802_1
（1）关于图8.8-4所示exm080802_1.mdl的说明
（2）
（3）
图8.8-5采用代数约束模块消减代数环影响的exm080802_2
（4）
图8.8-6采用单位延迟模块消减代数环影响的exm080802_3
（5）
%exm080802m.m
clearall
bdclose（'all'）%<2>
load_system（'simulink'）%<3>
tic;sim（'exm080802_1'）;T1=toc;
tic;sim（'exm080802_2'）;T2=toc;
t2='0.2';
open_system（'exm080802_3'）
set_param（'exm080802_3','MaxStep',t2）
tic;sim（'exm080802_3'）;T3=toc;
t002='0.002';
set_param（'exm080802_3','MaxStep',t002）
tic;sim（'exm080802_3'）;T4=toc;
disp（''）
disp（[blanks（31）,'仿真绝对耗时',blanks（5）,'仿真相对耗时']）
disp（['带代数环原模型',blanks（20）,num2str（T1）,blanks（12）,num2str
（1）]）
disp（['代数约束模块',blanks（22）,num2str（T2）,blanks（9）,num2str（T2/T1）]）
disp（['单位延迟阻断','MaxStep=',t2,blanks（9）,num2str（T3）,blanks（8）,num2str（T3/T1）]）
disp（['单位延迟阻断','MaxStep=',t002,blanks（7）,num2str（T4）,blanks（9）,num2str（T4/T1）]）
（a）含代数环的原模型
相对耗时1
（b）采用代数约束模块的修改模型
相对耗时0.25
（c）采用单位延迟模块的修改模型
最大步长取0.2
相对耗时0.18
（d）采用单位延迟模块的修改模型
最大步长取0.002
相对耗时0.34
图8.8-7采用代数约束模块的块图模型的输出曲线
.3S函数模块的创建和应用
.3.1S函数概述
.3.2S函数模块及其运作机理
（1）
（2）
图8.9-2
101开发S函数模块的一般步骤
.3.3M码S函数
101两个级别的M码S函数
102对二级M码S函数模版的注释
functionmsfuntmpl_basic（block）
%%msfuntempl_basic是二级M码S函数模版的基本型。
在大多数场合，该模版已够用。
%%用户使用该模版编写自己S函数时，绝不要沿用msfuntempl_basic名称，而应另起函数名。
%%更全面深入的模版是msfuntempl.m，它也驻留在{toolbox/simulink/blocks}文件夹上。
%%该主函数只包含如下一条指令，不得更改，不得添加。
setup（block）;
%endfunction
%%====设置Inputports、Outputports、Dialogparameters、Options等特性；必须有。
functionsetup（block）
%%
（1）设置输入输出口数目
block.NumInputPorts=1;
block.NumOutputPorts=1;
%%
（2）调用“运行对象”的SetPreCompInpPortToDynamic和SetPreCompOutPortInfoToDynamic
%%方法使模块的输入输出口继承信号的数据类型、维数、是否复数、采样模式
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;
%%（3）若模块对输入口某些属性有特别要求，则进行必要的重定义；否则，以下省略。
%%以下指令及其赋值仅是示例，用户应据需要改写。
block.InputPort
（1）.Dimensions=1;
block.InputPort
（1）.DatatypeID=0;%double
block.InputPort
（1）.Complexity='Real';
block.InputPort
（1）.DirectFeedthrough=true;%%true有直通通路；false则无。
%%（4）若模块对输出口某些属性有特别要求，则进行必要的重定义；否则，以下省略。
%%以下指令及其赋值仅是示例，用户应据需要改写。
block.OutputPort
（1）.Dimensions=1;
block.OutputPort
（1）.DatatypeID=0;%double
block.OutputPort
（1）.Complexity='Real';
%%（5）指定S函数模块的对话窗参数数目
%%以下赋值仅是示例，用户应据需要改写。
block.NumDialogPrms=0;
%%（6）指定采样时间，可取格式：
%%[0offset]，连续采样时间；[positive_numoffset]，离散采样时间；
%%[-1,0]，继承采样时间；[-2,0]，可变采样时间
%%以下赋值仅是示例，用户应据需要改写。
block.SampleTimes=[00];%%表示无偏移的连续采样
%%（7）指定仿真状态的保存和创建方法，可取选项：
%%'UnknownSimState'，先给出警告，然后采用默认设置；
%%'DefaultSimState'，采用内建模块的方法保存和重建连续状态、工作向量等
%%'HasNoSimState'，没有仿真状态要处理（如中模块不带输出口）
%%'CustomSimState'，通告Simulink有GetSimState和SetSimState方法实施
%%'DisallowSimState'，不允许保存和重建，若保存和重建则报错
block.SimStateCompliance='DefaultSimState';%%通常使用该指令及赋值。
%%（8）下面列出了块方法的最常用回调名（即单引号内的字符），它们是不可更改的。
%%函数句柄（即@及其后的字符）可以由用户自己命名，但必须与子函数名一致。
%%对于那些不需要的回调方法，用户应整行加以删除。
block.RegBlockMethod（'PostPropagationSetup',@DoPostPropSetup）;
%%设置Dwork向量的数目及其属性；仅含连续状态的S函数，不需要此回调。
block.RegBlockMethod（'InitializeConditions',@InitializeConditions）;
%%若仿真开始前及仿真过程中需要多次初始化，则使用该回调.
%%该回调对连续状态ContStates和/或Dwork向量赋初始值、配置内存等。
block.RegBlockMethod（'Start',@Start）;
%%若仅在仿真开始前需要初始化，则使用该回调。
block.RegBlockMethod（'Outputs',@Outputs）;%Required
%%任何S函数都必需该回调。
该回调计算S函数的输出，并存放于输出信号数组。
block.RegBlockMethod（'Update',@Update）;
%%若S函数有离散状态，或无直通通路，则需要该回调。
block.RegBlockMethod（'Derivatives',@Derivatives）;
%%若有连续状态，则需要该回调。
block.RegBlockMethod（'Terminate',@Terminate）;
%%二级M码S函数不必使用此回调。
%endsetup
%%====后向传递设置：
S函数含离散状态，或无直通通路时写该子函数。
functionDoPostPropSetup（block）
%%以下指令及其赋值仅是示例，用户应据需要改写。
block.NumDworks=1;
block.Dwork
（1）.Name='x1';
block.Dwork
（1）.Dimensions=1;
block.Dwork
（1）.DatatypeID=0;%double
block.Dwork
（1）.Complexity='Real';%real
block.Dwork
（1）.UsedAsDiscState=true;
%endDoPostPropSetup
%%====初始化条件：
当S函数需多次初始化时，才写该子函数
functionInitializeConditions（block）
%%（以下填写适当指令）
%endInitializeConditions
%%====启动：
当S函数仅需初始化一次，则应编写该子函数
functionStart（block）
%%以下指令行，仅是示例，用户应据需要编写。
block.Dwork
（1）.Data=0;
%endfunction
%%====输出：
任何S函数都必有该子函数
functionOutputs（block）
%%以下指令行，仅是示例，用户应据需要编写。
block.OutputPort
（1）.Data=block.Dwork
（1）.Data+block.InputPort
（1）.Data;
%endOutputs
%%====更新：
若S函数有离散状态，或无直通通路，则需要编写此子函数。
functionUpdate（block）
%%以下指令行，仅是示例，用户应据需要编写。
block.Dwork
（1）.Data=block.InputPort
（1）.Data;
%endUpda
升级会员 