Simulink交互式仿真(6).doc

8.1.1 系统平衡点和普通状态轨线图【例8.7-7】求图8.7-13所示两个模型的平衡点。模型(b)输入端比模型(a)多一个输入口。exm080707_1.mdlexm080707_2.mdl(a)(b)图8.7-13 待求平衡点的非线性系统块图模型（1）（2）xa=trim(exm080707_1,-0.1;-0.3)%xb=trim(exm080707_1,0;1) %xa = -0.8944 -1.7889xb = 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.6457ieig_Axb = 3.4110 -2.6222 （6）xa2,ua=trim(exm080707_2,-0.1;-0.3,0)%xb2,ub=trim(exm080707_2,0;1,1)%xa2 = -0.7487 -1.4974ua = 1.1974xb2 = 0.6810 1.3620ub = 1.6810 （2）% exm080707_1m.mclf;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);for k=1:nxx opts=simset(initialstate,xx(:,k);% t,x=sim(exm080707_1,10,opts);% plot(x(:,1),x(:,2); hold onendgrid on,hold offxlabel(x1);ylabel(x2)title(普通状态轨线)图8.7-14 多初始点出发的状态轨线和平衡点8.1.2 M码和Simulink模型的综合运用 1 单步仿真和精良状态轨线图【例8.7-8】绘制非线性系统块图模型的精良状态变化轨线。（1）function DX1,DX2,DP=exm080708_zzy(x1,x2,h)opts=simset(solver,ode5,fixedstep,h);%n=length(x1);DX1=zeros(n,n);DX2=DX1;DP=DX1;disp(正在逐点计算，请稍等！)for ii=1:n; for jj=1:n; opts=simset(opts,initialstate,x1(ii); x2(jj);% ,x=sim(exm080707_1,h,opts);% dx1=x(2,1)-x1(ii); dx2=x(2,2)-x2(jj); L=sqrt(dx12+dx22); DP(jj,ii)=L/h; if L1.e-10 DX1(jj,ii)=dx1;DX2(jj,ii)=dx2;% end endenddisp(计算结束)（2）%exm080708m.mh=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)shading interpalpha(0.5)colorbarhold onquiver(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)grid offhold offxlabel(x1),ylabel(x2)title(精良状态轨线斜率图)shg图 8.7-15 精良状态轨线迹斜率图 2 仿真模型和优化指令的协调【例8.7-9】题目背景：在迄今的自动控制教材中，凡讨论积分性能指标时，几乎总会提到所谓的ITAE传递函数标准型，并列出相应的分母多项式系数表。但值得指出的是：这些数据是20 世纪50年代初期，用模拟计算机仿真得到的。因此，这些数据的准确性带有明显的时代缺陷。与不同，ITAE性能函数无法解析计算，而只能通过数值计算进行。图8.7-16 计算的块图模型（1）问题的形成（2）（3）%exm080709m.mglobal a Jcamin=min(a0);na=length(a0);nd=na+2;opts=optimset(MaxFunEvals,300*na);CF=zeros(Kr,nd);Jk=zeros(1,Kr);for kk=1:Kr ar=a0+2*amin*（rand(1,na)-0.5);% a=fminsearch(exm080709_itae,ar, opts); cf=1,a,1; CF(kk,:)=cf; Jk(kk)=Jc;endJmin,kmin=min(Jk);cfmin=CF(kmin,:);% exm080709_itae.mfunction Jc=exm080709_itae(aa)global a Jca=aa;Tspan=0,0.1,20;%opts=simset(RelTol,0.0001);,Jt=sim(exm080709,Tspan,opts);%Jc=Jt(end);（4）clearKr=5;a0=3.25,6.60,8.60,7.45,3.95;exm080709mJmin,cfmin Jmin = 8.3338cfmin = Columns 1 through 6 1.0000 2.1519 5.6290 6.9338 6.7925 3.7398 Column 7 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)hold on,plot(tnew,ynew,r,LineWidth,3),hold offxlabel(t)title(ITAE 6阶新老标准型的阶跃响应比较)legend(Old,New,4),grid on图8.7-17 新老标准型的阶跃响应局部放大比较图表8.7-2 ITAE标准型新系数（黑体）和老“经典”系数（细体）对照阶次ITAE值传递函数分母多项式系数21.991.95191 1.4 11 1.5049 133.1443.13831 1.75 2.15 11 1.7828 2.1715 144.6264.59131 2.10 3.40 2.75 11 1.9521 3.3458 2.6473 157.1556.32151 2.80 5.00 5.50 3.40 11 2.0667 4.4976 4.6730 3.2568 169.6568.33381 3.25 6.60 8.60 7.45 3.95 11 2.1519 5.6290 6.9338 6.7925 3.7398 1715.00310.62901 4.48 10.42 15.05 15.54 10.64 4.580 11 2.2169 6.7433 9.3469 11.577 8.6778 4.3226 1818.68013.20511 5.20 12.80 21.60 25.75 22.20 13.30 5.15 11 2.2681 7.8313 11.8472 17.5325 16.0645 11.3094 4.8069 1l 研究表明：ITAE函数搜索空间的形状非常复杂，凹凸不平，小谷很多，许多地方深谷高峰相邻。要找到真正最小值点决非易事。虽可以肯定：单点标准型的新系数比老系数具有更小的ITAE值；但不能断言这新系数一定指示着最小值点。8.2 数值计算方面的考虑8.2.1 微分方程解算器Solver 1 ode45和ode23运作机理简要 2 ode113运作机理简要 3 ode15s和ode23s运作机理简要 4 不同解算器解Stiff方程的表现【例8.8-1】求微分方程在时的解。图8.8-1 微分方程的块图模型exm080801（1）关于exm080901.mdl的说明（2）syms t x xdxs=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(246 247 -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 = 18085ns2 = 101 图8.8-3 不同方法的解算结果比较8.2.2 积分步长和容差 1 积分步长的选择 2 计算容差的选择8.2.3 代数环问题 1 无惯性模块和代数环 2 消减代数环影响【例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.mclear allbdclose(all)%load_system(simulink)%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 采用代数约束模块的块图模型的输出曲线8.3 S函数模块的创建和应用8.3.1 S函数概述8.3.2 S函数模块及其运作机理（1）（2）图8.9-2 1 开发S函数模块的一般步骤8.3.3 M码S函数 1 两个级别的M码S函数 2 对二级M码S函数模版的注释function msfuntmpl_basic(block)% msfuntempl_basic 是二级M码S函数模版的基本型。在大多数场合，该模版已够用。% 用户使用该模版编写自己S函数时，绝不要沿用msfuntempl_basic名称，而应另起函数名。% 更全面深入的模版是msfuntempl.m ，它也驻留在toolbox/simulink/blocks文件夹上。% 该主函数只包含如下一条指令，不得更改，不得添加。setup(block);%endfunction% =设置Input ports、Output ports、Dialog parameters、Options等特性；必须有。function setup(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; % doubleblock.InputPort(1).Complexity = Real;block.InputPort(1).DirectFeedthrough = true;% true有直通通路；false则无。%(4)若模块对输出口某些属性有特别要求，则进行必要的重定义；否则，以下省略。% 以下指令及其赋值仅是示例，用户应据需要改写。block.OutputPort(1).Dimensions = 1;block.OutputPort(1).DatatypeID = 0; % doubleblock.OutputPort(1).Complexity = Real;%(5)指定S函数模块的对话窗参数数目% 以下赋值仅是示例，用户应据需要改写。block.NumDialogPrms = 0;%(6)指定采样时间，可取格式：%0 offset，连续采样时间；positive_num offset ，离散采样时间；%-1, 0 ，继承采样时间；-2, 0，可变采样时间% 以下赋值仅是示例，用户应据需要改写。block.SampleTimes = 0 0;%表示无偏移的连续采样%(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函数不必使用此回调。 %end setup%= 后向传递设置：S函数含离散状态，或无直通通路时写该子函数。function DoPostPropSetup(block)% 以下指令及其赋值仅是示例，用户应据需要改写。block.NumDworks = 1;block.Dwork(1).Name = x1;block.Dwork(1).Dimensions= 1;block.Dwork(1).DatatypeID = 0; % doubleblock.Dwork(1).Complexity = Real; % realblock.Dwork(1).UsedAsDiscState = true;%end DoPostPropSetup%= 初始化条件：当S函数需多次初始化时，才写该子函数function InitializeConditions(block)%（以下填写适当指令）%end InitializeConditions%= 启动：当S函数仅需初始化一次，则应编写该子函数function Start(block)% 以下指令行，仅是示例，用户应据需要编写。 block.Dwork(1).Data = 0;%endfunction%= 输出：任何S函数都必有该子函数function Outputs(block)% 以下指令行，仅是示例，用户应据需要编写。 block.OutputPort(1).Data = block.Dwork(1).Data + block.InputPort(1).Data;%end Outputs%= 更新：若S函数有离散状态，或无直通通路，则需要编写此子函数。function Update(block)% 以下指令行，仅是示例，用户应据需要编写。 block.Dwork(1).Data = block.InputPort(1).Data;%end Update%= 导数计算：连续时间状态更新function Derivatives(block)% 以下填写适当指令%end Derivatives%= 终止：对C MEX S函数必需，但对二级M码S函数则不必。function Terminate(block)%end Terminate 3 二级M码S函数模块设计示例【例8.9-1】为图8.9-3所示单摆设计一个M码S函数模块，并进而利用该模块构建一个在“周期方波力”作用下，单摆的摆动块图模型。该模块以外力为输入，摆角为输出，等效摩擦系数、等效重力系数以及摆的初始条件都以S函数模块的对话窗参数出现。图 8.9-3 单摆示意图（1）（3）function exm080901_simpend(block)% exm080901_simpend is a M-file S-function rewritten from template msfuntmpl_basic.m.% Produced by zzy.% The setup method is used to set up the basic attributes. setup(block);% end simpendzzy2function setup(block)% (1)Register number of input and output portsblock.NumInputPorts = 1;block.NumOutputPorts = 1;% (2)Setup port properties to dynamically inheritedblock.SetPreCompInpPortInfoToDynamic;block.SetPreCompOutPortInfoToDynamic;% (5