《连续系统仿真》第五次课后讨论记要

《连续系统仿真》第五次课后讨论记要讨论主题：水箱实例建模分析时间：2015.4.11地点：A区动力楼508讨论人：廖立杨鹏志沈融分工：廖立和沈融负责资料收集和记要整理，杨鹏志负责ppt制作讨论内容：水箱问题的模型建立过程，编写程序进行分析1.仿真对象图1.1水箱实例进口流量：由阀门1和2调节。水位控制器：接受来自传感器的水位信号，并与设定值进行比较，然后发出控制信号调节阀门3的开度。控制器的作用是保证容器的液位。出口流量：由阀门3调节。2.对象模型的建立2.1质量守衡方程质量平衡方程：容器内质量的变化等于进口质量流量减去出口质量流量，其数学式表示为：321(2.1)式中：m容器内液体质量（kg）,21WW和是进口质量流量，3W是出口质量流量（kg/s）。2.2流动方程vpyCWv333(2.2)式中：阀门3的前后压差Δp(Pa)，容器中液体的比容υ(m3/kg)，阀门开度y3，定义为全开时的流动面积与工作时的流动面积之比。vplg（2.3）如果容器中液体上部的压力和出口压力均为大气压的话，压差只取决于液体的高度l。g是重力加速度（9.81m/s2）。假定阀门1和2的压差的改变没有很大影响，密度也不变化。那么流量21WW和只是决定于阀门的开度，y1，y2：1'11yCWv（2.4）2'22yCWv（2.5）每一个开度yi依赖于两个因素：阀程、流动面积11xy（2.6）22xy（2.7）233xy（2.8）2.3阀运动方程对于电动阀门而言，每一个阀门都由一个伺服器来驱动的，或称为调节器、执行机构。以使阀程从xi到要求值xdi。总需要一定的时间才能开启/关闭阀门，采用一个最简单的动力学模型，即是一个一阶指数滞后：1111xxdtdxd（2.9）2222xxdtdxd（2.10）3333xxdtdxd（2.11）i是与第i个阀程调节器相关的时间常数，通常是几秒钟。2.4控制器阀门3的阀程受水位控制器的支配。假设水位控制器完全是比例作用，以致能使要求的阀程sdllkx3(2.12)式中：l是测量的水位（m）；ls是设定的水位（m）；k是比例常数（m-1）。水位与以下因素有关：容器的横截面积A，液体的密度υ，容器中液体质量m，表示为：Amvl（2.13）2.5初始条件0.330.220.110xxxxxxmm（2.14）2.6计算（2.15）3.仿真模型包含的内容3.1数学模型中的变量常量AgvkCCCvvv,,,,,,,,,3213'2'1初始条件0.30.20.10,,,xxxm代数变量lxyyyp微分dtdxdtdxdtdxdtdm321,,,状态变量321,,,xxxm3.2方程转换mvAgxCxCxCdtdmmAkvkldtdxxxdtdxxxdtdxvvvsdd2332'21'133322222111113.3矩阵形式的方程mxxxx321sddlxxu21uxfuxfuxfuxfuxf,,,,,4321将微分方程组用矩阵形式的方程表述，选择数值积分算法，求解常微分方程组：ttutxf,,dtdx————ttdtttutxftxtx0,,0实质：将常微分方程（组）变换成差分方程（组）。4.四阶龙格库塔法求解4.1所给数据onstantsC常量'1vC'2vC3vCkvgA123)/(skg)/(skg)(2m)(1m)/(3kgm)/(2skg)(2m)(s)(s)(s1010004517.00.1001.08.933610Forcingiablesvar1dx2dxsl)(m5.05.02Stateiablevar状态变量1x2x3xm5.05.0707.08121ebraicsAlg1W2W3Wpl1y2y3y3dx)(kg)(kg)(kg)(Pa)(m551026529707.25.05.05.0707.04.2四阶龙格库塔算法原理四阶龙格库塔算法的具体公式为：112341213243/622t,t/2,/2*,t/2,/2*,t,nnnnnnnnnnXXhkkkkkfXkfhXhkkfhXhkkfhThk计算过程：先分别求出k1,k2,k3,k4，,然后再将其代入到1nX中即可。4.3具体求解过程将所给常量代入微分方程组，选择四阶龙格库塔数值积分算法，求解常微分方程组：mxxxdtdmmdtdxxdtdxxdtdx3001.08.9004517.01010103001.011021665.0335.0232132211(1)状态变量1x051015202530354045500.10.20.30.40.5迭代次数状态变量x1控制器比例系数k=1051015202530354045500.10.20.30.40.5迭代次数状态变量x1控制器比例系数k=0.5图4.1状态变量1xA:当控制器比例系数k=1时：选取初值x10=0.1，迭代步长为1，迭代不到20次，状态变量x1达到稳定,x1稳定值为0.5，与所给数据一致。B:当控制器比例系数k=0.5时，选取初值x10=0.1，迭代步长为1，状态变量x1达到稳定的迭代次数和稳定值与控制器比例系数k=1的情况相同。(2)状态变量2x051015202530354045500.10.20.30.40.5迭代次数状态变量x2控制器比例系数k=1051015202530354045500.10.20.30.40.5迭代次数状态变量x2控制器比例系数k=0.5图4.2状态变量2xA：当控制器比例系数k=1时：选取初值x20=0.1，迭代步长为1，迭代不到30次，状态变量x2达到稳定,x2稳定值为0.5，与所给数据一致。B:当控制器比例系数k=0.5时，状态变量x2达到稳定的迭代次数和稳定值与比例系数k=1的情况相同。(3)状态变量3x050001000015000-0.500.511.5迭代次数状态变量x3控制器比例系数k=105000100001500000.511.5迭代次数状态变量x3控制器比例系数k=0.5图4.3状态变量3xA：当控制器比例系数k=1，选取初值30x=0.5，迭代步长为1，大概迭代15000次左右，状态变量3x达到稳定，3x稳定值为0.707，与所给数据一致。B：当控制器比例系数k=0.5，，选取初值30x=0.5，迭代步长为1，大概迭代12000次左右，状态变量3x达到稳定，3x稳定值为0.707。(4)状态变量mC:控制器比例系数k=1，状态变量3x，选取初值30x=0.5，迭代步长为1，大概迭代15000次左右，状态变量3x达到稳定，3x稳定值为0.707，与所给数据一致0500010000150000.40.50.60.70.80.9迭代次数状态变量x3图4.3状态变量3xD:控制器比例系数k=1，状态变量m，选取初值m=6800，迭代步长为1，大概迭代15000次左右，状态变量m达到稳定，m稳定值为6002，与所给数据有误差。050001000015000520054005600580060006200640066006800状态变量m图4.3状态变量m（2）将控制器比例系数改为k=0.5，求解结果：A:状态变量1x051015202530354045500.10.20.30.40.5迭代次数状态变量x1控制器比例系数k=1051015202530354045500.10.20.30.40.5迭代次数状态变量x1控制器比例系数k=0.5B:状态变量2x051015202530354045500.10.20.30.40.5迭代次数状态变量x2控制器比例系数k=1051015202530354045500.10.20.30.40.5迭代次数状态变量x2控制器比例系数k=0.5C:状态变量3x050001000015000-0.500.511.5迭代次数状态变量x3控制器比例系数k=105000100001500000.511.5迭代次数状态变量x3控制器比例系数k=0.5D:状态变量m05000100001500050005500600065007000迭代次数状态变量m控制器比例系数k=105000100001500050005500600065007000迭代次数状态变量m控制器比例系数k=0.5程序：%状态变量x1，k=1%分别用x，y，z，w表示x1,x2,x3,m（x1=x，x2=y，x3=z，m=w）x=0.1;y=0.1;z=0.5;w=6800;h=1;n=50;t=0;A=[t,x,y,z,w];fori=1:n-1%i=1:n定义龙格库塔程序k1=0.5/3-x/3;l1=0.5/6-y/6;m1=-0.2+1*0.001/30*w;n1=10*x+10*y-0.004517*z^2*sqrt(9.8*w/0.003);k2=0.5/3-(x+k1*h/2)/3;l2=0.5/6-(y+l1*h/2)/6;m2=-0.2+1*0.001/30*(w+n1*h/2);n2=10*(x+k1*h/2)+10*(y+l1*h/2)-0.004517*(z+m1*h/2)^2*sqrt(9.8*(w+n1*h/2)/0.003);k3=0.5/3-(x+k2*h/2)/3;l3=0.5/6-(y+l2*h/2)/6;m3=-0.2+1*0.001/30*(w+n2*h/2);n3=10*(x+k2*h/2)+10*(y+l2*h/2)-0.004517*(z+m2*h/2)^2*sqrt(9.8*(w+n2*h/2)/0.003);k4=0.5/3-(x+k3*h)/3;l4=0.5/6-(y+l3*h)/6;m4=-0.2+1*0.001/30*(w+n3*h);n4=10*(x+k3*h)+10*(y+l3*h)-0.004517*(z+m3*h)^2*sqrt(9.8*(w+n3*h)/0.003);x=x+h*(k1+2*k2+2*k3+k4)/6;y=y+h*(l1+2*l2+2*l3+l4)/6;z=z+h*(m1+2*m2+2*m3+m4)/6;w=w+h*(n1+2*n2+2*n3+n4)/6;t=t+h;A=[A;t,x,y,z,w];endsubplot(2,1,1)plot(A(:,1),A(:,2),'b-+','MarkerSize',3)%xx1legend('控制器比例系数k=1')xlabel('迭代次数');ylabel('状态变量x1');holdon%状态变量x1，k=0.5%分别用x，y，z，w表示x1,x2,x3,m（x1=x，x2=y，x3=z，m=w）x=0.1;y=0.1;z=0.5;w=6800;h=1;n=50;t=0;A=[t,x,y,z,w];fori=1:n-1%i=1:n定义龙格库塔程序k1=0.5/3-x/3;l1=0.5/6-y/6;m1=0.5*(-0.2+1*0.001/30*w);n1=10*x+10*y-0.004517*z^2*sqrt(9.8*w/0.003);k2=0.5/3-(x+k1*h/2)/3;l2=0.5/6-(y+l1*h/2)/6;m2=0.5*(-0.2+1*0.001/30*(w+n1*h/2));n2=10*(x+k1*h/2)+10*(y+l1*h/2)-0.004517*(z+m1*h/2)^2*sqrt(9.8*(w+n1*h/2)/0.003);k3=0.5/3-(x+k2*h/2)/3;l3=0.5/6-(y+l2*h/2)/6;m3=0.5*(-0.2+1*0.001/30*(w+n2*h/2));n3=10