储油罐的变位识别与罐容表标定.doc

2010高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：A题 储油罐的变位识别与罐容表标定摘要通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照题意要求本文对题中所给数据进行分析，再利用Matlab曲线拟合，最后得出重新标定的罐容表。针对第一个问题，通过对倾斜油罐的罐内五种不同高度的油的容量的定性分析，得到罐内五种不同油高的油的体积的的上界值（如图a），由此可得，油罐的体积随油高度的变化函数图像从光滑的变为非光滑的。通过对已有数据进行整理分析。利用最小二乘法对其进行三次、四次、五次的曲线拟合（程序见附表1），得出来的结果和实际数据相比较可得，三次函数的拟合效果最符合实际数据（见表1）.再利用MATLAB的多项式求值即可得出每隔1厘米罐容表的值。针对第二个问题，本文在第一问的基础上对储油罐的体积进行分块。根据分出来的五个部分（如图a）建立分块求体积的基本模型.利用Mathematica对各部分分别积分求体积，最后将其依次相加，即可得到储油罐体积的表达式。在求各个部分体积与高度，偏角之间的关系时，我们采用逐步添加变量的方法，即先考虑纵向变位角度和垂直于水平面的油位高度与罐容的关系，再根据与油位高度的关系可以再次引入因变量，最终得到了罐容与关系的分段表达式。在求解和的值的时候，我们将表中的数据根据罐容表达式分为四组，然后用最小二乘法对表达式拟合，得到参量的四组估计值，取其平均值即可得的最优值。关键词：Matlab曲线拟合 最小二乘法 Mathematica曲线拟合 微积分方程方法一、问题重述通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。二、模型假设1储油罐完好无损，进、出油过程中无渗漏现象；2储油罐的厚度忽略不计；3储油罐内无油泥沉积在油罐底部和管壁上；4油位探针正常工作；5不考虑测量其间油品的挥发和热胀冷缩情况；6不考虑油浮子的宽度，即认为油浮子是一个点，油的浮力只对其中心起作用；7忽略压强对油密度的影响，即认为罐内的油不论有多少，密度适中处处相等；8油压对容器壁的弹性形变忽略不计。三、符号说明 油罐内的存油量与油位高度之间的函数关系式 油罐侧面被油覆盖的面与油位高度之间的函数关系式 油罐内油的高度 小椭圆形储油罐的总体积 实际储油罐的总体积四模型的分析与建立4.1问题一的分析与解答：4.11罐体变位后对罐容表的影响设油罐内的存油量与油位高度之间的函数关系式为，油罐侧面被油覆盖的面与油位高度之间的函数关系式为（油位高度为油表探针测得的高度）则：1）油罐变位之前，且（为油罐的长度）2）油罐变位之后，由于当油位探针测得的数据为零时其左侧仍有一部分油存在，因此不可能再与成正比。罐体变位后对罐容表的影响即为与的关系以及与的关系。下面首先求出小椭圆储油罐的体积：通过分析可得：油罐发生变位后，可以把油罐体积分成五部分（如下图a）利用微积分中的平行截面公式法逐个求出这五部分的体积，再依次把这五部分的体积相加即为油罐的总体积V。五个部分区域如下图a所示（图b为图a的一个切面图形）：油位探针图a水平线 椭圆柱体体积 ： 图b0.4m2.05m水平线V2V4油位探针V1V3以下采用数学方法计算：1体积： 底面边界曲线方程为：作截面垂直轴，则与楔形体交面是个矩形，其截面面积为：-（1） 因而所求体积为：-（2）图c2 体积：根据（1）式及（2）式可求得： 图d3体积： 同理 根据（1）式及（2）式可求得4体积：储油罐内油品液面高度的函数，某一液面高度下，罐内油品体积根据积分的概念，体积元素： 、 、 由于： 因此：V4max=3.598958()（注：以上所求积分结果均用mathematica算出。）综上所述，所建立的体积模型为：（1）（2）（3）（4）（5）由此可以看到，罐体变位后，罐内高度与储油量函数图像为非光滑的、分段、连续。用最小二乘法来进行一次模拟，结果如下：Out3= 从数据可以看出来，这次的拟合结果与本文之前拟合出来的函数接近，又和实际数据相对比可知：在第一问模拟出来的三次函数与实际的误差较小。因此，此次模拟的函数是合理的。4.12. 罐体变位后油位高度间隔为1cm的罐容表标定值的确定为了确定罐体变位后油位高度间隔为1cm的罐容表标定值，首先求出罐容量随油位高度的变化关系函数，再将每隔一厘米的数据带入即可得到罐容的标度值。（1）数据的整理：对表1中的倾斜变位进油和倾斜变位出油两组数据进行整理，得到罐内实际油量如下表：表一：倾斜变位进油倾斜变位出油实际罐内油量（/L）油位高度(/mm)实际罐内油量(/L)油位高度(/mm)962.86411.293464.741020.651012.86423.453414.741007.731062.86438.333364.74994.321112.86450.543314.74980.961162.86463.903264.74967.101212.86477.743214.74956.011262.86489.373164.74941.541312.79502.563114.74929.691362.79514.693064.74916.441412.73526.843014.74904.141462.73538.882964.74891.901512.73551.962914.74879.231562.73564.402864.74868.991612.73576.562814.74855.131662.73588.742764.74844.021712.73599.562714.74831.641762.73611.622664.74820.471812.73623.442614.74808.161862.73635.582564.74796.001912.73646.282514.74785.041962.73658.592464.74773.072012.73670.222414.74762.092062.73680.632364.74750.812112.73693.032314.74739.422162.73704.672264.74727.092212.73716.452214.74715.322262.73727.662164.74705.432312.73739.392114.74693.522362.73750.902064.74682.502412.73761.552014.74671.022462.73773.431964.74658.682512.73785.391914.74647.742562.73796.041864.74635.762612.73808.271814.74624.612662.73820.801764.74612.532712.73832.801714.74600.692762.73844.471664.74589.402812.73856.291614.74577.002862.73867.601564.74564.582912.73880.061514.74554.332962.73892.921464.74540.763012.73904.341414.74528.653062.73917.341364.74517.193112.73929.901314.74504.873162.73941.421264.74490.783212.73954.601214.74478.063262.73968.091164.74465.973312.73980.141114.74452.403362.73992.411064.74439.983412.731006.341014.74425.833462.731019.07964.74411.733512.731034.243514.741035.36（2）数据的拟合：根据实际情况和表1中的数据可的应为单调递增函数，且增长速率先快后慢，即先大于0，再等于0，最后是大于0。我们尝试用matlab拟合函数。由于三次项系数为负数的三次函数，四次，五次函数中均有一段符合此函数特征，故用多项式拟合的方法，应用MATLAB拟合出的多项式系数（程序见附录1）如下表：表二：项目x出油倾斜三次00-0.00290.61080.0563125.1129出油倾斜四次00-0.00490.8203-9.4058280.4761出油倾斜五次0-0.00020.0264-1.325662.5031-6660.31进油倾斜三次00-0.00250.53944.49142.1041进油倾斜四次00-0.00320.61910.8525102.4026进油倾斜五次0-0.00010.0132-0.525339.5775-408.598再次利用MATLAB软件中的多项式求值函数可以得到拟合函数对表中实际高度参数的函数值（程序见附录2），如下表：表三：倾斜出油时油位高度三次拟合时得到的油量四次拟合时得到的油量五次拟合时得到的油量倾斜出油时实际油量倾斜进油时油位高度三次拟合时得到的油量四次拟合时得到的油量五次拟合时得到的油量倾斜进油时实际油量1020.653.41032.6559-7.72443464.74411.290.96530.96210.9628962.861007.733.36582.6484-7.43213414.74423.451.00971.00561.00611012.86994.323.31872.6383-7.14813364.74438.331.06481.05981.05941062.86980.963.27082.626-6.88423314.74450.541.11071.10481.10341112.86967.103.22022.6109-6.62973264.74463.901.16171.15481.15161162.86956.013.17912.5971-6.43963214.74477.741.21521.20721.20161212.86941.543.12462.5769-6.2093164.74489.371.26071.25171.24371262.86929.693.07932.5587-6.03423114.74502.561.31281.30271.29131312.79916.443.02812.5364-5.85323064.74514.691.36131.351.33491362.79904.142.97992.5141-5.69823014.74526.841.41031.39781.37831412.73891.902.93142.4904-5.55592964.74538.881.45931.44541.42111462.73879.232.88072.4643-5.42062914.74551.961.51291.49751.46711512.73868.992.83942.4421-5.31992864.74564.401.56441.54731.51031562.73855.132.7832.4106-5.19512814.74576.561.6151.59631.5521612.73844.022.73742.384-5.10442764.74588.741.6661.64551.5931662.73831.642.68622.3533-5.01262714.74599.561.71151.68931.62871712.73820.472.63972.3244-4.93772664.74611.621.76261.73831.66771762.73808.162.58822.2915-4.86362614.74623.441.81281.78641.7051812.73796.002.53712.258-4.79862564.74635.581.86461.83591.74211862.73785.042.49082.2268-4.74672514.74646.281.91051.87951.77381912.73773.072.442.1919-4.69692464.74658.591.96331.92971.80891962.73762.092.39332.159-4.65712414.74670.222.01341.97711.84062012.73750.812.34512.1245-4.6222364.74680.632.05832.01951.86772062.73739.422.29642.089-4.59212314.74693.032.11192.06981.89822112.73727.092.24352.0497-4.56562264.74704.672.16222.1171.9252162.73715.322.1932.0115-4.54572214.74716.452.21322.16451.95012212.73705.432.15061.9789-4.53292164.74727.662.26172.20961.97212262.73693.522.09951.9391-4.52192114.74739.392.31252.25651.99292312.73682.502.05221.9018-4.51582064.74750.902.36232.30242.01092362.73671.022.00291.8624-4.51332014.74761.552.40822.34452.02542412.73658.681.95011.8196-4.51481964.74773.432.45952.39122.03892462.73647.741.90331.7813-4.51951914.74785.392.51092.43792.04952512.73635.761.85231.7389-4.5281864.74796.042.55662.47922.05622562.73624.611.80491.6992-4.53891814.74808.272.60892.52622.06072612.73612.531.75381.6559-4.55361764.74820.802.66232.57382.06142662.73600.691.70391.6133-4.57081714.74832.802.71322.61892.05822712.73589.401.65651.5725-4.58961664.74844.472.76242.66232.05132762.73577.001.60481.5275-4.61261614.74856.292.81212.70572.04022812.73564.581.55331.4823-4.6381564.74867.602.85932.74672.02552862.73554.331.51111.4451-4.66051514.74880.062.9112.79122.00472912.73540.761.45571.3957-4.69231464.74892.922.9642.83651.97762962.73528.651.40661.3518-4.72241414.74904.343.01062.8761.94873012.73517.191.36061.3103-4.75231364.74917.343.06332.92021.90993062.73504.871.31171.2659-4.78561314.74929.903.11372.9621.86623112.73490.781.25631.2154-4.82521264.74941.423.15962.99961.82043162.73478.061.20691.1702-4.86211214.74954.603.21143.04181.76123212.73465.971.16051.1275-4.8981164.74968.093.26393.08381.69263262.73452.401.10921.0801-4.93911114.74980.143.31023.12041.62423312.73439.981.0631.0373-4.97741064.74992.413.35673.15671.54743362.73425.831.01120.989-5.02161014.741006.343.40883.19671.45113412.73411.730.96050.9418-5.0662964.741019.073.45573.23211.35413462.731034.243.51063.27271.22713512.731035.363.51463.27561.21723514.74:（注：表中B，C，D，G，H，I列的数据皆乘以）在上表中明显可以看出四次和五次函数对两个函数的拟合程度都不好，所以选择用三次曲线去拟合这倾斜变位进油时和倾斜变位出油时油罐油量随油标高度的变化函数解析式为：其具体图形如下：图e：图f： 倾斜变位出油由于油桶的高度为1.2米，故可以设定0到120共121个分位点，利用MATLAB的多项式函数求值（程序详见附录3），带入上面的两个函数表达式可以得到在两组不同的数据下罐体变位后油位高度间隔为1cm的罐容表的罐容量，即现在罐容表的标度。由于两组数据的测定都是在同一个罐中进行，且间隔时间不长，所以为了得到精确的罐容表标度，可以将得到的两组数据平均，结果如下表：表四：油位高度（cm）进油时罐容表（L）出油时罐容表(L)罐容表（均值） (L）油位高度（cm）进油时罐容表（L）出油时罐容表(L)罐容表（均值） (L）042.10125.1083.60611755.701743.101749.40147.10125.8086.45621798.201785.401791.80253.20127.6090.40631840.801827.801834.30360.40130.7095.55641883.601870.301876.95468.50134.90101.70651926.401913.001919.70577.70140.30109.00661969.401955.701962.55687.90146.80117.35672012.501998.602005.55799.10154.40126.75682055.602041.402048.508111.30163.20137.25692098.802084.302091.559124.40173.00148.70702142.002127.302134.6510138.50183.90161.20712185.302170.202177.7511153.40195.80174.60722228.602213.102220.8512169.30208.70189.00732271.902256.002271.9013186.20222.70204.45742315.102298.902307.0014203.80237.70220.75752358.402341.602350.0015222.40253.60238.00762401.602384.302392.9516241.80270.50256.15772444.702426.902435.8017262.10288.30275.20782487.702469.402478.5518283.10307.10285.10792530.702511.802521.2519305.00326.80315.90802573.502553.902563.7020327.70347.40337.55812616.302596.002606.1521351.10368.80359.95822658.902637.802648.3522375.40391.10383.25832701.302679.402690.3523400.30414.20407.25842743.602720.802732.2024426.00438.20432.10852785.702762.002773.8525452.40463.00457.70862827.602802.902815.2526479.60488.50484.05872869.302843.502856.4027507.40514.80511.10882910.702883.802897.2528535.90541.90538.90892952.002923.902937.9529565.00569.70567.35902992.902963.602978.2530594.80598.20596.50913033.603002.903018.2531625.20627.40626.30923074.003041.903057.9532656.20657.30656.75933114.103080.503097.3033687.90687.90687.90943153.903118.703136.3034720.10719.10719.60953193.403156.503174.9535752.90751.00751.95963232.503193.903213.2036786.20783.40784.80973271.303230.803251.0537820.10816.50818.30983309.603267.303288.4538854.50850.10852.30993347.603303.303325.4539889.40884.30886.851003385.203338.703361.9540924.80919.00921.901013422.403373.703398.0541960.70954.30957.501023459.103408.103433.6042997.00990.10993.551033495.403442.003468.70431033.801026.301030.051043531.203475.303503.25441071.001063.101067.051053566.503508.003537.25451108.701100.301104.501063601.303540.103570.70461146.701137.901142.301073635.603571.603603.60471185.201175.901180.551083669.403602.403635.90481224.001214.401219.201093702.703632.603667.65491263.101253.201258.151103735.403662.103698.75501302.701292.401297.551113767.503690.903729.20511342.501332.001337.251123799.003719.003759.00521382.701371.901377.301133829.903746.403788.15531423.101412.101417.601143860.303773.003816.65541463.801452.601458.201153889.903798.903844.40551504.901493.401499.151163919.003824.003871.50561546.101534.401540.251173947.403848.303897.85571587.601575.801581.701183975.103871.703923.40581629.301617.301623.301194002.103894.403948.25591671.301659.001665.151204028.403916.203972.30601713.401701.001707.204.2问题2的分析与解答：在问题一的基础上我们把变位后的储油罐分成两个部分，一部分是中间的圆柱体部分，另一部分是左右两侧的球冠部分。由于当储油罐倾斜变位后，两边球冠部分进油量与出油量不相同，因此对其两部分分别求解。故，问题二可归结为分别对两边和中间三部分求解。体积函数的自变量均取为油位高度。最终，求完三部分后，可归为一个函数表达式，即为体积与油位高度的函数关系式。具体算法如下所述：图f1对中间的圆柱体进行求解：回归到第一问中所建立的模型，相应图f中的对应于 图b中的V1,V2,V3,V4。下面对做合理的计算：（1）的求解（如图g）可类似于第一问当中的V1的求解方法：2tan/2- 图g 注：其中的R为圆柱底面半径,（2）对求解（如图h），与的求解方法及方程相同，只是取值范围不同。 图h（3）对的求解（如图i），也与的求解方法及方程相同，只是的取值范围不同。 图i（4）对的求解（如图j），求解方法如下：如图中所示，建立空间直角坐标系，平面平行于水平面，z轴垂直与水平面。根据数学分析中对立体图形求体积的方法，用平行截面面积法求体积。截得的图形为梯形，梯形面积=（上底+下底）高/2,再根据体积公式 （注：此公式出自数学分析上册第十章第二节）可得： 图2 再对两边球状体进行求体积。（1）先对左边的球状体求体积，如图k所示，以水平面为xy平面建立的空间直角坐标系，用平行截面法来求冠状球体内油的体积，截面为图中的圆M，圆M与球冠的切面相交而得的面S即为平行截面，对其关于z轴求积分便可得到球冠内油的体积。由于的角度较小，所以可得以下关系：，由球内的几何关系可求得由此数据可以得到圆M的半径为圆心M的坐标为（0.625，0，z）；所以截面M的方程为；球冠的法向量为（，故根据点法式可以得到切面方程为由此可得截面面积积分为：由此可得球冠内油的体积为：根据空间立体几何的关系最终求得：zx球冠的切面0.625R1rOABCD MYS 图k（2）对于右边的球冠体积，根据求左边的球冠所求的体积及左右所对应的关系，可得到体积为：（3） 以上求得了各个分段区域的体积，其中的各自变量均可根据旋转角、转化为以h即油表高度来表示，因此各个部分的体积转化为以h为自变量的函数。V1除外，因为在其区域内有油量时，油表高度此时为0，即又还没浸没油表的最低刻度。求得各个区域以h为自变量的函数表达式为: 综上：可以建立体积模型为： 当h没开始取值时，当时， 当时，当时，根据实际问题，油罐中的油一般很难达到油表得得满刻度，所以V8可以不予考虑，即第四种情况一般不考虑。当时，由于h不再变化，所以不已考虑。3对于以上建立的体积模型，其中h是变量，、是变位参数。第一、 根据附表二的数据，可以分别求得、。具体算法如下：分别代入两组数据，即将h和的变化量的带入相应h取值范围的体积模型中，建立二元一次方程组，算得、的粗值。但是，当我们取值的时候，具有偶然性，并且一组数据不足以代表整个附录中的数据，这样算误差会很大，经计算，得，这种方法比较传统，由于数据量较大，计算比较麻烦，所以应该寻求最优算法。第二、由于体积函数是非线性函数，要计算出其中的两个参数，所以选择用最小二乘法，经过有限次的迭代最终估算出参数的最优值。具体做法如下：将原题表2中的数据根据上述分段函数分出四组，分别用MATLAB中的非线性最小二乘拟合函数lsqno