18.1.1 系统平衡点和普通状态轨线图【例 8.7-7】求图 8.7-13 所示两个模型的平衡点。模型(b)输入端比模型(a)多一个输入口。2124x&212)4(xu&exm080707_1.mdl exm080707_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.89441.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.74872-1.4974ua =1.1974xb2 =0.68101.3620ub =1.6810 （2）% exm080707_1m.mclf;xx=-2,-1, 0, 1, 1, 1, 11, 1, 1, 1, 0,-1,-2;set_param(exm080707_1,InitInArrayFormatMsg, None)nxx=size(xx,2);for k=1:nxxopts=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 多初始点出发的状态轨线和平衡点38.1.2 M 码和 Simulink 模型的综合运用1 单步仿真和精良状态轨线图【例 8.7-8】绘制非线性系统 块图模型的精良状态变化轨线。2124x&（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-10DX1(jj,ii)=dx1;DX2(jj,ii)=dx2; %endendenddisp(计算结束)（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 off4hold offxlabel(x1),ylabel(x2)title(精良状态轨线斜率图)shg图 8.7-15 精良状态轨线迹斜率图2 仿真模型和优化指令的协调【例 8.7-9】题目背景：在迄今的自动控制教材中，凡讨论积分性能指标时，几乎总会提到所谓的 ITAE 传递函数标准型，并列出相应的分母多项式系数表。但值得指出的是：这些数据是 20 世纪 50 年代初期，用模拟计算机仿真得到的。因此，这些数据的准确性带有明显的时代缺陷。与 不同，ITAE 性能函数02)()(dteISEJ无法解析计算，而只能通过数值计算进行。dteITAJ0)(图 8.7-16 计算 的块图模型)(ITAEJ（1）问题的形成5（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:Krar=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 61.0000 2.1519 5.6290 6.9338 6.7925 3.7398Column 71.0000 （5）old=tf(1,1,a0,1);new=tf(1,cfmin);yold,told=step(old,50);ynew,tnew=step(new,50);6plot(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 on0 2 4 6 8 10 12 14 16 1800.20.40.60.81tITAE 6三三三三三三三三三三三三三OldNew图 8.7-17 新老标准型的阶跃响应局部放大比较图表 8.7-2 ITAE 标准型新系数（黑体）和老 “经典”系数（细体）对照阶次 ITAE值传递函数分母多项式系数2 1.991.95191 1.4 11 1.5049 13 3.1443.13831 1.75 2.15 11 1.7828 2.1715 14 4.6264.59131 2.10 3.40 2.75 11 1.9521 3.3458 2.6473 15 7.1556.32151 2.80 5.00 5.50 3.40 11 2.0667 4.4976 4.6730 3.2568 16 9.6568.33381 3.25 6.60 8.60 7.45 3.95 11 2.1519 5.6290 6.9338 6.7925 3.7398 17 15.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 18 18.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 1 研究表明：ITAE 函数搜索空间的形状非常复杂，凹凸不平，小谷很多，许多地方深谷高峰相邻。要找到真正最小值点决非易事。虽可以肯定：单点标准型的新系数比老系数具有更小的 ITAE 值；但不能断言这新系数一定指示着最小值点。78.2 数值计算方面的考虑8.2.1 微分方程解算器 Solver1 ode45 和 ode23 运作机理简要2 ode113 运作机理简要3 ode15s 和 ode23s 运作机理简要4 不同解算器解 Stiff 方程的表现【例 8.8-1】求微分方程 在 时的解。09.10xx&0)(,1x&图 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)80 1 2 3 4 5 6 7 8 9 10-0.012-0.01-0.008-0.006-0.004-0.0020txdx=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