2013全国大学生数学建模竞赛 清华大学获奖论文A,3728

1 卧式储油罐的变位识别与罐容表标定模型 摘 要 本文针对卧式储油罐变位后罐容重新标定问题建立相应的数学模型 问题一 首先对罐内储油分五种情况进行讨论 利用定积分分段求出不同油位高 度时的实际油量 即罐内储油体积 并依此建立储油罐变位后罐容表重新标定的数学 模型 其次 用附件 1 的实验数据计算模型的理论值 从而拟合出变位后的罐容标定 偏差 发现理论值和实验值拟合得较好 最后给出变位后的罐容表标定值 问题二 罐体的变位分为纵向和横向两个方向 首先将油浮子探测高度转换为等 效的纵向倾角为的油位高度 再分五种情况对纵向变位做类似问题一的讨论 并依 此建立变位后罐容表标定的数学模型 然后根据附件 2 中进油前的实际数据做非线性 最小二乘拟合 并确定变位参数 进而利用附件 2 中进油后的实际数4 57 2 56 据对模型做单样本 Kolmogorov Smirnov 检验 显著性水平达到 0 861 以上 从而较好 的验证了模型的合理性和可靠性 最后给出变位后的罐容标定表 关键词 罐容标定 罐体变位 非线性最小二乘拟合 单样本 K S 检验 2 1 问题的提出 通常加油站都有若干个储存燃油的地下储油罐 并且一般都有与之配套的 油位 计量管理系统 采用流量计和油位计来测量进 出油量与罐内油位高度等数据 通过 预先标定的罐容表 即罐内油位高度与储油量的对应关系 进行实时计算 以得到罐 内油位高度和储油量的变化情况 许多储油罐在使用一段时间后 由于地基变形等原因 使罐体的位置会发生纵向倾 斜和横向偏转等变化 以下称为变位 从而导致罐容表发生改变 按照有关规定 需要定期对罐容表进行重新标定 请用数学建模方法研究解决储油罐的变位识别与罐 容表标定的问题 1 为了掌握罐体变位后对罐容表的影响 利用小椭圆型储油罐 两端平头的椭 圆柱体 分别对罐体无变位和倾斜角为 4 10的纵向变位两种情况做了实验 实验 数据如附件1所示 请建立数学模型研究罐体变位后对罐容表的影响 并给出罐体变位 后油位高度间隔为1cm的罐容表标定值 2 对于实际的储油罐 试建立罐体变位后标定罐容表的数学模型 即罐内储油 量与油位高度及变位参数 纵向倾斜角度 和横向偏转角度 之间的一般关系 请利 用罐体变位后在进 出油过程中的实际检测数据 根据所建立的数学模型确定变位参数 并给出罐体变位后油位高度间隔为10cm的罐容表标定值 进一步利用附件中的实际检 测数据来分析检验模型的正确性与方法的可靠性 2 模型假设及符号约定 2 1 模型的假设 1 储油罐是刚体 2 忽略油的物理变化 3 忽略注油管 出油管和油位探测装置的体积 2 2 符号的约定 椭圆长半轴 a 椭圆短半轴 b 油位探针到储油罐主体较近侧面的距离 m 油位探针到储油罐主体较远侧面的距离 n 储油罐主体的半径 同时也是球冠底面的半径 r 球冠厚度 d 球冠所在球体的半径 R 油罐左下角到左端油面点的距离 h 油位高度 H 储油罐主体的长度 l 油面所在平面 储油罐的油体的空间 3 储油罐横向偏转时的等效油位高度 1 H 无变位储油罐中油的体积 简记 HVV 变位储油罐中油的体积 简记 HVV 无变位储油罐左边的球冠在下面油的体积 简记 HV冠Hz 冠 V 无变位储油罐的圆柱体在下面油的体积 简记 HV柱Hz 柱 V 变位储油罐左边的球冠中油的体积 简记 1 HV冠 1冠 V 变位储油罐右边的球冠中油的体积 简记 2 HV冠 2冠 V 变位储油罐的圆柱体中油的体积 简记 HV圆柱 圆柱 V 4 3 模型的准备 3 1 命题 1 由和构成的图形的面积1 2 2 2 b by a x 20 bhhy 2 1 1 1 1arcsin 2b h b h b h abS 证明 h dy b by aS 0 2 12 h bydbyb b a 0 22 2 db b a bh b 22 2 bh b b bb b a arcsin 2 12 222 2 1 1 1 1arcsin 2b h b h b h ab 3 2 对附件 1 数据的处理 无变位进油数据表 将油罐初值跟累加进油量相加得出罐内实际油量 将油位高 度的单位转化为米 对于倾斜变位进油的数据也作同样的处理 部分数据见下表 完 整数据见附件 1 表 无变位进油数处理后部分数据 油位高度 m 0 1590 1760 1930 2090 2240 2390 2540 2680 282 罐内油量 L 312362412462512562612662712 表 倾斜变位进油数处理后部分数据 油位高度 m 0 4110 4230 4380 4500 4630 4770 4890 5020 514 罐内油量 L 962 861012 861062 861112 861162 861212 861262 861312 791362 79 无变位出油数据表 由于无法确切知道罐内的实际油量 但根据出油量跟相应油 高的变化 得出罐容对油高的差商 对于倾斜变位出油的数据也作同样的处理 表 无变位出油数处理后部分数据 油位高度 m 1 151 121 101 081 061 041 021 01 罐容对油高的 差商 L m 1870 562189 142422 482610 972767 022905 293022 973130 87 图 1 5 表 倾斜变位出油数处理后部分数据 油位高度 m 1 021 010 990 980 970 960 940 930 93 罐容对油高的 差商 L m 3869 973728 563742 513607 504508 573455 434219 413773 584065 04 3 3 Kolmogorov Smirnov 检验 单样本 K S 检验是一种拟合优度的非参数检验方法 单样本 K S 检验是利用样本 数据推断总体是否服从某一理论分布的方法 适用于探索连续型随机变量的分布形态 单样本 K S 检验可以将一个变量的实际频数分布与正态分布 Normal 均匀分布 Uniform 泊松分布 Poisson 指数 Exponential 分布进行比较 其零假设 H0 为样本来自的总体与指定的理论分布无显著差异 SPSS 实现 K S 检验的过程如下 1 根据样本数据和用户的制定构造出理论分布 差分布表得到相应的理论累计概 率分布函数 0 F x 2 利用样本数据计算各样本数据点的累计概率得到检验累计概率分布函数 0 xS 3 计算和在相应的变量值点上的差 得到差值序列 单样本 0 xF 0 xS x D xD K S 检验主要对差值序列进行研究 SPSS 在统计中将计算 K S 的 Z 统计量 并依据 K S 分布表 小样本 或正态分布 表 大样本 给出对应的相伴概率值 如果相伴概率小于或等于用户的显著性水平 则应拒绝零假设 H0 认为样本来自的总体与指定的分布有显著差异 如果相伴概 率值大于显著性水平 则不能拒绝零假设 H0 认为样本来自的总体与指定的分布无显 著差异 4 问题一的分析 4 1 罐体无变位 储油的示意图如图 2 所示 用油体的纵向横 截面积与储油罐的长度 相乘计算罐容 Sl 4 2 罐体发生纵向变位 以油罐左下角为原点建立笛卡尔坐标系 储油罐中的储油如图 3 所示 图 2 6 储油罐由方程由 ly b bz a x 0 1 2 2 2 2 确定 油面所在的平面方程 coscossinhzy 罐容 罐体里面的油量 其中为由 dvV 1 2 2 2 2 b bz a x ly 0 coscossinhzy 所围成的空间 bh20 油体在平面上的截面如图 4 所示 yz 容易确定 tanmHh 与轴交于点 y 0 cot 0 hQ 与轴交于点 z 0 0 hP 与直线相交于点 Rz x 2 0 2 cot 2 0 RRhR 与直线相交于点 Ly x0 tan 0 lhlT 于任意一点处作垂直于轴的平面 它截油体的截面面积 cot 0 hy y 2 1 1 1 1arcsin 2 b z b z b z abzS 图 3 图 4 7 当时罐内并非无油 当油较高时 0 H 油面与储油罐的位置关系如图 5 所示 当 时油罐内也并非全满 纵向倾斜的储油RH2 罐容油量的计算可分为 5 种情况 1 当时 罐内有油 小于一临界值 0 H 但不确定 2 当时 油面截储油罐成三角形 tan0nH 罐容dyysHV h cot 0 dz b z b z b z ab h 0 2 1 1 1 1arcsin 2 cot 1 3 1 1 arcsin 2 cot 2 3 2 2 1 22 kVkkkkkab 其中 1 tan 1 b mH b h k 3 当时 tan2tanmbHn 沿轴方向 延长椭圆柱体 可分为两部分 和的差 其中 可归y 1 2 1 2 为 2 中的情况 对 1 1 tan 1 1 b mH b h k 对 2 1 tan 1 tan 2 b nH b lh k 2121 kVkVVVV 4 当时 bHmb2tan2 平面将分成两部分 其中在左边的是椭圆圆柱体 右边部分可归 cot 2 bhy 结为 4 2 3 中的情况 其中 11 2 1 b b k tan1 2 b l k cot 2 21121 kVbhab kVkVVVVHV 5 当时 罐内油未满 大于一临界值 但不确定 bH2 5 问题一的模型 5 1 罐体无变位时的数学模型 图 5 图 5 8 由 4 1 的分析 罐体无变位时的数学模型 2 1 1 1 1arcsin 2b H b H b H ablSlV 将以上模型得出的理论数据与附件 1 无变位进油处理后的数据作比较 结果如图 6 所 示 图 6a 图 6b 由图 6 可得到理论值稍偏大于处理后的实验值 这可能是由于我们忽略了测量仪 器等实验器材对油位高度所造成的影响 从无变位出油数据表我们无法知道罐内的实际油量 但根据出油量跟相应油高的 变化 可以得出罐容对油高的差商 为了进一步对检验模型的灵敏度 求关于 HV 的导数 得H 2 112 b H al dH HdV 由上式得出的理论数据与附录 1 无变位出油处理后的差商数据作比较 结果如图 7 所 示 9 图 7a 图 7b 5 2 罐体发生纵向变位的数学模型 由 4 2 的分析 建立小椭圆储油罐发生纵向变位时的数学模型 其中 2 3 2 2 1 22 1 3 1 1 arcsin 2 cot kkkkkabkV 2 1arcsin 2 cot kkkab dH HdV tan 1 nH tan2 2 mbH 1 tan b mH k 1 tan 1 b mH k 1 tan 2 b nH k tan1 3 b nm k 4 0 m05 2 n 2 78 1 a 2 2 1 b 360 21 4 1 4 利用以上模型给出罐体变位后油位高度间隔为 1cm 的罐容表标定值 bHHkVbmHab HHHkVkV HHkV HV 2 cot tan 0 23 2121 1 10 油位高度 cm 12 363738 118119 油量 L 3 506 30 812 20849 97888 15 3981 104026 22 详细罐容表见附录 表 1 和处理小椭圆储油罐无变位时相同的方式 对模型进行检验 利用 MATLAB 得到如 下图 8 和图 9 图 8a 图 8b 由图 8 得到理论值偏大于实验值 误差值服从偏大型正态分布 发生该现象的原因可 能是我们忽略了油位探针等仪器以及储油罐的形变所引起的体积影响 引进纵向倾斜角变量以研究罐体变位后对罐容表的影响 由 图 9a图 9b 11 22 V HV H V H H V HVHHV 原先制定的罐容表 在发生纵向形变后 表中某一高度对应的罐容数值与实际同 一高度对应的罐容数值的误差主项 H V 6 问题二的分析 如图 10 在储油罐上建立笛卡尔坐标系 图 10 不难确定 d rd R 2 22 球冠 1 方程 0 2222 yd RrzdRyx 球冠 2 方程 2 222 dlyl RrzdRLyx 圆柱方程 ly rrzx 0 222 油面方程 coscossinhzy 从而 0 0 hP 0 cot 0 hQ 2 cot 2 0 RRhR tan 0 LhLT 12 球冠在间的部分的体积为Hz 0 Hz yd RrzdRyx 0 0 2222 0 0 222 d H dzrzdRyRdyHV冠 arccos1 4 arccos 1 3 1 arctan2sin arctan 12 1 arctan 23 arctan 23 1 2r 2 222 2 22 22 3 kkk rd dr k K K MM K Kk kK M K k M K 其中 1 r H k 22 2 dr rd K 2 22 1 1 k Kk M 圆柱在间的部分的体积为ly 0 ly rrzx 0 222 22 1 1 1 1arcsin 2r H r H r H lrHV 柱 6 1 当罐体只发生纵向倾斜 在实际情况中罐体的纵向倾斜角一般都是小角度 与 4 2 的分析相似 分 5 种情 况讨论罐容和油位高度的关系 1 当时 罐内有油 油量小于某一临界值 0 H 2 当时 tan0nH 圆柱冠 VVV 1 图 6 1 2 图 6 1 13 其中 tan 2 0tan 0 222 1 mHVdzrzdRyRdyV d yh 冠冠 tan 2 tan 0 22 cot 0 mHVdzrzrdyV yhh 柱圆柱 3 当时 tan2tanmbHn 图 6 1 3 21冠圆柱冠 VVVV tan 1 mHVHV 冠冠 其中 tan 2 mHVHV 冠冠 tantanmHVmHVV 柱柱圆柱 4 当时bHmb2tan2 图 6 1 4 21冠圆柱冠 VVVV 3 3 2 1 3 1 3 1 RdRdRV 冠 其中 tan tan cottan2 2 nHVmHVmRHrV 柱柱圆柱 tan 2 mHVV 冠冠 5 当时 罐内油未满 大于一临界值 但不确定 bH2 6 2 当罐体只发生横向偏转 14 只考虑横向变位 小偏角 取过探针且与油面垂直的平面为截面 储油罐发生横 向偏转倾斜 不影响油位探针与储油罐的相对位置 因此 发生横向偏转后 油位探 针仍过正截面圆的圆心 如图 6 2 1 2 所示 AC 是液面 AB 是油位高度 CD 是实际高 度 AB CD 过圆心 O ACCD 下面推导油位高度与实际高度 半径 横向倾斜角的关系H 1 Hr 1 油面高于圆心 OB ODDC OB OC OBAOAB coscos 即 cos 1 rH rH 2 油面低于圆心 coscos CDOD OB OC OBOAOBAB 即 coscos 11 rH r Hr rH 当油面过圆心时 关系式显然成立 r rH rH cos 1 综上所述 或 cos 1 rH rH cos 1 rHrH 7 问题二的模型 7 1 模型的建立 当储油罐发生纵向倾斜和横向偏转时 将储油罐横向逆转使得储油罐相当于只发 生纵向倾斜 油位高度与实际高度的关系 H 1 H cos 1 rHrH 将实际高度当作等效油位高度 现在储油罐相当于只发生纵向倾斜 且油位高度为 1 H 根据 6 1 的讨论 建立罐体变位后标定罐容表的数学模型 1 H 图 6 2 1 图 6 2 2 15 2tan2 tan tan tan cottan2 3 1 3 1 tan2tan tan tan tan tan tan0 tan tan cos 111 11 2232 111 11 111 1 rHmrmHVmHV mHVmRHrRdRdR mrHnmHVmHV mHVmHV nHmHVmHV rHrH HV 柱柱 冠 柱柱 冠冠 柱冠 其中 HV冠 arccos1 4 arccos 1 3 1 arctan2sin arctan 12 1 arctan 23 arctan 23 1 2r 2 222 2 22 22 3 kkk rd dr k K K MM K Kk kK M K k M K 1 r H k 22 2 dr rd K 2 22 1 1 k Kk M 22 1 1 1 1arcsin 2r H r H r H lrHV 柱 1 5 1 6 2 drnm 7 2 模型参数的确定 由上述油位高度与储油量的数学模型知 罐内油量是关于油位高度的非线 V HH 性函数 传统的最小二乘法无法准确估计模型的参数 故采用非线性最小二乘拟 合来估计模型的未知参数 其 MATLAB 实现如下 非线性最小二乘拟合是已知输入向量 xdata 和输出向量 ydata 并且知道输入与输 出的函数关系为 ydata F x xdata 但不知道系数向量 x 今进行曲线拟合 求 x 使得 下式成立 i ii ydataxdataxFydataxdataxF 2 2 2 x 2 1 2 1 min 16 函数 lsqcurvefit 格式 x resnorm residual exitflag output lambda lsqcurvefit fun x0 xdata ydata lb ub options 其具体程序见附录 最终求得问题二模型的变位参数为 4 57 56 2 7 3 模型的检验 现利用附件 2 中进油后的实际数据通过 SPSS 软件对模型做单样本 K S 检验 其 结果为 表一 Descriptive Statistics Percentiles NMeanStd DeviationMinimumMaximum25th50th Median 75th 实际观测值30131272 988216056 131825036 2657704 79 16800 230031282 550045538 9650 模型计算值30132917 078912004 3479613148 8752706 11 22014 335033403 360043268 7500 表二 One Sample Kolmogorov Smirnov Test 实际观测值模型计算值 N301301 Mean31272 988232917 0789Normal Parametersa b Std Deviation16056 1318212004 34796 Absolute 075 074 Positive 075 074 Most Extreme Differences Negative 070 062 Kolmogorov Smirnov Z1 2931 292 Asymp Sig 2 tailed 071 071 a Test distribution is Normal b Calculated from data 表三 One Sample Kolmogorov Smirnov Test 2 实际观测值模型计算值 17 N301301 Minimum5036 2613148 87Uniform Parametersa b Maximum57704 7952706 11 Absolute 033 035 Positive 033 035 Most Extreme Differences Negative 027 024 Kolmogorov Smirnov Z 569 602 Asymp Sig 2 tailed 903 861 a Test distribution is Uniform b Calculated from data 由表二和表三知 双尾渐进显著性水平都大于 0 故认为由非线性最小二乘拟合估 计出的参数是合理的 即模型二是合理的 7 4 间隔为 10cm 的灌容表标定值 罐体变位后油位高度间隔为 10cm 的罐容表标定值 油位高度 m 0 10 2 1 51 61 7 2 82 9 油量 L3392 15930 1 35553 438304 740798 3 60879 860916 9 详细罐容表见附录 表 2 8 模型的优缺点和改进 问题一的模型 优点 模型可以调节的可控参数较多 具有一定的普遍性和灵活性 由图 6 9 可以 看出模型可以较好的拟合实验数据 缺点 未考虑储油罐的变形及油的物理变化 如温度影响 并且理论值总是稍偏大于 实验值 问题二的模型 优点 模型具有一定的普遍性和灵活性 容易推广 容易利用计算机实现 缺点 变位参数被限定为小角度 18 改进 可以考虑变位参数被限定为大角度 也可以考虑主体为圆柱体 两端 为球冠体的罐体变位后标定罐容表的数学模型 参考文献 1 华东师范大学数学系 数学分析 M 北京 高等教育出版社 2008 2 吕林根 许子道 解析几何 M 北京 高等教育出版社 2008 3 陈佰军 左振滨 编制卧式贮油罐容量表的计算机模型 J 黑龙江北一农垦大学学 报 1993 No 2 56 60 4 薛薇 统计分析与 SPSS 的应用 M 北京 中国人民大学出版社 2008 19 附附 录 录 表 1 油位 高度 cm 油量 L 油位 高度 cm 油量 L 油位 高度 cm 油量 L 油位 高度 cm 油量 L 油位 高度 cm 油量 L 13 5025429 66 491330 11732362 27973341 69 26 3026461 49 501371 88742405 37983378 51 310 0027494 00 511413 85752448 37993414 86 414 8028527 14 521456 02762491 261003450 72 520 7029560 90 531498 35772534 021013486 06 627 9030595 25 541540 85782576 641023520 87 736 3031630 15 551583 50792619 121033555 11 846 1032665 58 561626 28802661 421043588 77 957 4033701 53 571669 19812703 551053621 81 1070 1034737 96 581712 21822745 491063654 20 1184 4035774 86 591755 32832787 221073685 91 12100 3036812 20 601798 52842828 741083716 92 13117 7037849 97 611841 80852870 021093747 17 14136 9038888 15 621885 13862911 061103776 64 15157 8239926 72 631928 51872951 831113805 27 16180 2640965 66 641971 93882992 331123833 01 17204 00411004 95652015 37893032 531133859 82 18228 91421044 58662058 82903072 431143885 62 19254 88431084 53672102 28913112 001153910 33 20281 86441124 79682145 71923151 231163933 86 21309 76451165 34692189 13933190 111173956 06 22338 54461206 16702232 50943228 611183981 10 23368 14471247 23712275 82953266 721194026 22 24398 53481288 56722319 09963304 42 20 表 2 油位高度 m油量 L油位高度 m油量 L油位高度 m油量 L 0 13392 11 125279 92 149310 0 0 25930 11 227730 42 251487 9 0 38425 51 330262 02 353643 3 0 41042 841 432871 02 455739 6 0 512356 71 535553 42 557718 9 0 614314 01 638304 72 659476 9 0 716337 71 740798 32 760772 4 0 818442 61 842865 22 860879 8 0 920634 21 944981 72 960916 9 1 022913 72 047134 93 0 模型二算法程序 建立数据 M 文件 Data 2010A bigchange Data 2304 14 53747 27 2490 06 57782 53 1476 98 31680 00 1218 44 24411 40 2290 87 53439 37 994 32 18338 22 842 18 14435 27 2454 51 57054 65 2348 13 54750 40 653 11 9941 94 1338 49 27766 55 2404 05 55984 22 1004 54 18607 78 2610 29 60065 11 2280 46 53196 16 2067 14 47928 87 1038 89 19519 72 637 24 9586 64 1583 65 34700 98 2190 91 51047 20 524 27 7172 62 1958 83 45083 88 2382 50 55514 45 1231 17 24764 66 1457 25 31121 10 1291 34 26443 03 21 1642 93 36374 74 2116 53 49191 94 1878 40 42914 84 1696 61 37883 79 795 32 13280 49 707 24 11181 08 1830 91 41615 55 910 97 16172 90 2127 37 49465 96 899 88 15889 56 1067 66 20290 19 726 89 11640 87 821 69 13927 29 1269 19 25823 64 2474 40 57464 67 1913 71 43872 31 2548 47 58927 69 613 59 9064 22 2510 23 58185 31 1350 95 28117 30 1114 10 21545 60 1095 36 21037 34 2209 13 51492 32 1170 06 23075 57 1432 35 30416 13 2528 01 58534 01 1777 51 40140 74 1735 97 38984 85 2158 40 50243 62 1977 87 45590 80 751 99 12235 51 1989 59 45901 47 771 35 12699 59 1594 13 34997 30 1932 64 44382 33 2024 06 46808 88 2362 44 55070 68 1856 54 42318 32 1818 96 41286 70 2587 60 59657 02 1000 建立函数 M 文件 function y fun2010A 3 x xdata y 169 64 pi 169 32 atan 25 169 4 3 x 1 2 3 x 2 xdata 22 3 2 2 1 4 3 x 1 2 3 x 2 xdata 3 2 2 1 2 15 32 acos 2 3 x 2 xdata 3 2 4 3 x 1 15 32 a