《自控系统数字仿真实验.doc》由会员分享，可在线阅读，更多相关《自控系统数字仿真实验.doc（33页珍藏版）》请在金锄头文库上搜索。
1、实验一面向微分方程的数值积分法仿真一、实验目的1掌握数值积分法的基本概念、原理及应用；2用龙格-库塔法解算微分方程，增加编写仿真程序的能力；3分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；4.对数值算法中的“病态问题”进行研究。二、实验内容1、已知系统微分方程及初值条件取步长h=0.1，试分别用欧拉方程法和RK4法求t=2h时的y值，并将求得的值与解析解比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原因。（编程完成；选用MATLAB ode函数完成。） 程序如下：t0=0;tf=2;h=0.1;y1=1;y2=1;y3=1;t1=0;t2=0;t3=0;n=rou
2、nd(tf-t0)/h;for i=1:n %欧拉法 y1(i+1)=y1(i)+h*(2*h+y1(i); t1=t1,t1(i)+h;endfor i=1:n %四阶龙格-库塔法 k1=y2(i)+t2(i); k2=y2(i)+h*k1/2+t2(i)+h/2; k3=y2(i)+h*k2/2+t2(i)+h/2; k4=y2(i)+h*k3+t2(i)+h; y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6; t2=t2,t2(i)+h;endfor i=1:n y3(i+1)=2*exp(t3(i)-t3(i)-1; t3=t3,t3(i)+h;endplot(
3、t1,y1,r,t2,y2,g,t3,y3,k)分析:红线为用欧拉法得到的结果，绿线为用四阶龙格库塔法得到的结果，蓝线为根据解析方程得到的结果。其差异原因主要有两个：（1）二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格库塔法是根据四阶微分方程得到的；（2）由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。 2、已知系统的传递函数为在单位阶跃输入下，系统响应的解析解为试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：1）比较两种数值积分解与解析解的逼近程度；（绘图）2）改变步长，分析步长对数值解精度的影响；3）不断加大步长，分析计算稳
4、定性的变化。程序如下： a=0 1 0; 0 0 1; -22.06 -27 -10;b=0;0;1; c=40.6 0 0;X1=0;0;0; t=0;Y1=0;X=0; u=1; Y2=0; Y3=0;X2=0;0;0; x=0; h=0.1;t0=0; tf=2; t1=0; t2=0; t3=0;N=(tf-t0)/h;for i=1:N k1=a*X1+b; k2=b+a*(h*k1/2+X1); k3=b+a*(h*k2/2+X1); k4=b+a*(h*k3+X1); X1=X1+h*(k1+2*k2+2*k3+k4)/6; Y1=Y1,c*X1; t1=t1,t1(i)+h;e
5、ndfor i=1:N x=X2(:,i)+h*(a*X2(:,i)+b*u); y=c*x; X2=X2,x; Y2=Y2,y; t2=t2,t2(i)+h;endfor i=1:N y=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.24*i); Y3=Y3,y; t3=t3,t3(i)+h;endplot(t1,Y1,r,t2,Y2,g,t3,Y3,b) 当h=0.01时，仿真图如下： 可见当步长很小时，数值积分解逼近解析解，其数值解得的精度也很高。 当h=0.1时，仿真图如下： 可见随着步长的增加，其数值积分解与解析解的逼近程
6、度减小，精度降低。 当h=0.3时，仿真图如下： 可见其数值积分解已变得不稳定了，即步长不能取得太大3、求下图所示系统的阶跃响应y(t)的数值解。（v=1,k=1,开始时间t0=0,结束时间tf=10,h=0.25）分析k、v对系统响应的影响。（编程用RK4求解；Simulink）(修正：G(s)表达式的分母中，0.25s+1改为（0.25s+1）的二次方) 程序如下： k=1; a=conv(1 0 0,conv(0.25 1,0.25 1) b=2*k k; X0=0 0 0 0 ;v=1; n0=4; tf=10; h0=0.25; r=1; V=v; n=n0; T0=0; Tf=tf
7、; h=h0; R=r; b=b/a(1); a=a/a(1); A=a(2: n+1); %首一化处理 A=rot90(rot90(eye(n-1,n); -fliplr(A); %形成能控标准型A阵 B=zeros(1,n-1),1; %形成输入阵B(n维列向量)m1=length(b); %分子系数向量维数m+1 C=fliplr(b),zeros(1,n-m1); %形成输出阵C(n维行向量) Ab=A-B*C*V; %形成闭环系数阵AbX=X0;y=0;t=T0; %赋初值，准备开始递推运算 N=round(Tf-T0)/h; %确定输出点数for i=1:N %四阶龙格-库塔法 k
8、1=Ab*X+B*R; k2=Ab*(X+h*k1/2)+B*R; k3=Ab*(X+h*k2/2)+B*R; k4=Ab*(X+h*k3)+B*R; %求各阶次斜率 X=X+h*(k1+2*k2+2*k3+k4)/6; %求状态 y=y,C*X; %求输出并以向量形式保存 t=t,t(i)+h; %输出对应时刻以向量形式保存endt,y %输出数据形式结果plot(t,y) %输出曲线形式结果 k=0.5 K=1，v=1 K=1，V=0.5 K=1,v=5时， K=2，v=1时分析：当k取值增大，v值不变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k增加到一定程度时系统便发散了（即不
9、稳定了）。当v值增大，k值不变时，波头也是变多变陡，当v值增大到一定程度时系统便不稳定了。4、已知系统状态状态方程为（修正：x的系数矩阵中最后一个元素为-40，而非40）采用RK4法，步长分别取h=0.01、0.04、0.06，求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。（提示：病态系统）程序如下：a=-21 19 -20;19 -21 20;40 -40 -40;x=1;0;-1;X=x;t=0;t0=0;tf=2;h=0.01;n=round(tf-t0)/h;for i=1:n x=X(:,i)+h*a*X(:,i); X=X,x; t=t,t(i)+h;endb=X
10、(1,:);c=X(2,:);d=X(3,:);plot(t,b,r,t,c,g,t,d,b) h=0.01时 H=0.02时 H=0.04时分析：如图，当h=0.01时，在t=0.2s以后系统输出便趋于平稳，当取h=0.02时，系统输出振荡剧烈，趋于稳定的时间也变长，当取h=0.04后，系统输出呈发散振荡形式。当h=0.06后系统仍然是发散的，即当h的取值改变时，原先稳定的系统变得不稳定了，这便是病态系统。但是这个结果与书上的不同。三、思考题1、在进行仿真计算时，是否选用的数值积分法的阶次越高越好？ 不是，因为阶次越高，虽然精度高了，但其计算的复杂程度也越高，有时甚至无法完成计算。2、选用数
11、值积分法进行仿真的原则。(1) 精度 仿真结果的精度主要受三项误差影响： 截断误差：由算法本身的精度阶次所决定； 舍入误差：由计算机字长决定； 累积误差：由以上两项误差随计算时间长短累积情况决定。(2) 计算速度 取决于所用数值方法和计算步长。(3) 稳定性 数值稳定性主要与计算步长h有关，不同的数值方法对h都有不同的稳定性限制范围，且与被仿真对象的时间常数有关。一般所选步长与系统最小时间常数有以下数量级关系： 而多步法、隐式算法有较好的数值稳定性，在对稳定性较注重时，应予以优先选择。 实验二第一部分面向结构图的数值积分法仿真一、实验目的加深理解连续系统面向结构图仿真的原理及特点，进一步掌握数字积分法解算微分方程的方法，增加编写仿真程序的能力。二、实验内容1、用面向方框图的数字仿真方法对下列系统进行仿真。2、求解下图所示系统在f=-1(t)阶跃扰动作用下第、第环节的动态过程。分别用面向框图的数值积分法（RK4法）、MATLAB中有关系统建模的命令和Simulink三种方法求解。三实验程序及结果1.P=0,0.07,1,0.14;1,0.012,1,0;0,0.05,1,0.15;10,1,1,0;1,0.01,0,0.0008;WIJ=1,0,1;1,4,-1;2,1,1;3,2,1;3,5,-1;4,3,1;5,4,1n=5;Y0=1;Yt0=0
