数模全国一等奖储油罐的变位识别与罐容表标定

储油罐的变位识别与罐容表的标定
摘要
本文研究储油罐的变位识别与罐容表的标定。

分别以小椭圆型油罐和实际卧式储油罐为研究对象，运用高等数学的积分的知识，分别建立罐体变位前后罐内油体积与油高读数之间的积分模型，使用Matlab 软件得出结论。

对于问题一，以小椭圆型储油罐为研究对象，在无变位时，小椭圆型储油罐为规则的椭球柱体，可利用解析几何与高等数学的知识建立油罐内体积与油高读数之间的积分模型，得出罐体无变位时的理论值。

当罐体发生纵向变位时，小椭圆型储油罐的截面不再是规则的几何形体，但根据倾角α及所给小椭圆型罐体的尺寸，可得其截面面积的表达式，利用高等数学中积分的方法，根据不同油高，建立了模型一，得到了储油量和油高的关系公式。

最后，根据实验数据的处理，用拟合的方法，修正了某些系统误差的影响，计算出罐体变位后油位高度间隔1cm 的罐容表的标定值。

对于问题二，由于实际储油罐内没油的高度不同，我们将其分为五种情况分别讨论，并对每种情况建立积分公式，得出罐内油体积与油位高度及变位参数（纵向倾斜角α和横向偏转角β）之间的函数关系式，利用所给的实验数据，运用最小二乘法，建立非线性规划模型
2
1
2
arg ,(((,,)(,,)))min (,,)n
i
i i i V H V H
OilData error OilData αβ
αβαβαβ-==--∑用Matlab 非线性规划求解得出使得总体误差最小的α与β值：α=2.12°，β=4.06°。

通过α与β的数值计算出出油量理论值与实测值的平均相对误差小于0.5% 。

对模型进行了较为充分的正确性验证和稳定性验证：在α与β的值为0时，其计算出来的罐容值与理论值完全吻合，说明模型在体积计算上是正确的；当对油高进行0.1%的扰动时，α的值变化也在0.1%左右，说明α的稳定性很好，但是β的值从4.06°变成了3.75°，变化了大约8%，所以我们详细分析了β的数学表达式，从理论上分析了影响其稳定性的因素。

根据得到的变位参数计算出实际罐体变位后油位高度间隔为10cm 的罐容表的标定值。

最后，本文对模型的优缺点进行了评价，并讨论模型的推广。

关键字：储油罐；变位识别；罐容表标定；非线性规划
一．问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。

许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。

按照有关规定，需要定期对罐容表进行重新标定。

根据上述所述，求解下列问题：
（1）为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为=4.1°的纵向变位两种情况做了实验。

请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。

（2）对于实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度）之间的一般关系。

请利用罐体变位后在进/出油过程中的实际检测数据，根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。

进一步利用实际检测数据来分析检验你们模型的正确性与方法的可靠性。

二．问题分析
本文研究罐容表的读数与储油罐的变位的关系。

借助高等数学积分的方法，求出储油量与油高读数的函数关系式，并对倾斜的储油罐进行容量标定。

1.对问题一的分析
问题一中用小椭圆储油罐分别对罐体无变位和纵向倾斜进行实验，研究变位对罐容表的影响，因此我们分别建立变位前和变位后的罐容表读数与罐内油体积的函数关系式，通过函数关系式计算出理论值，再与所给的实际值相比较，得出其相对误差，然后通过分析系统误差进行修正，出罐体变位后油位高度间隔为1cm的罐容表的标定值。

2.对问题二的分析
问题二中是以实际储油罐为研究对象，不仅考虑了储油罐的纵向倾斜，而且还考虑了横向偏转，为了使问题简化，我们先只考虑纵向倾斜，由于储油罐的形体不规则，所以我们将它分成如图1所示的三部分，分别算出每部分的体积与罐容表读数的函数关系式，然后对其求和。

再考虑横向偏转，建立它与所给的油高的函数关系式。

然后将二者进行综合考虑得出变位后罐容表读数与储油罐内油体积的函数关系式，通过关系式和所给数据，运用最小二乘法，通过MATLAB程序，搜索出α和β的最小误差解，再对模型的稳定性和正确性进行评定，最后给出高度间隔10cm的罐容表的标定值。

图1 油罐分区域积分示意图
三．模型假设
假设一：数据是储油罐的内壁参数。

假设二：忽略温度、压力对汽油的密度的影响。

假设三：储油罐在偏移的过程中，油位探针始终与油罐底面垂直。

假设四：对卧式储油罐来说，不考虑其长期埋在地下所发生的蠕变。

假设五：累加进出油量数据是准确可靠的。

四．符号说明
H: 对应于罐容表读数的液面实际高度。

H: 球冠中与油罐圆柱左侧底面距离为x处的油高。

1
R: 球冠中与油罐左侧底面相距为x处的小圆半径。

2
H：球冠中与油罐圆柱右侧底面距离为x处的油高。

2
R：球冠中与油罐右侧底面相距为x处的小圆半径。

3
R: 储油罐圆柱部分的底面半径。

1
R: 球冠所在球体的大圆半径。

H:第i条数据所对应的罐容表读数。

i
OilData：用于分析的油量进出数据。

a: 椭圆长半轴长。

b: 椭圆短半轴长。

n: 用于分析的进出油测量数据个数。

h：罐容表读数。

五． 模型的建立与求解
5.1 模型一的建立与求解
问题一要求研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm 的罐容表标定值。

5.1.1 计算未变位和变位的理论罐内油位高度与储油量的关系
利用高等数学中微元法求体积的方法建立罐容表读数与罐内油体积的函数关系式的模型。

（1） 在无变位的情况下，储油罐内的油所占空间为柱体，其体积为 V S L = (1) 其中S 为柱体底面面积，L 为柱体的长度。

2h
b
S x dy -=⎰
(2)
底面椭圆方程为 22
22
1x y a b += (3)
22a x b y b
=
- (4)
将（4）代入（2），得到
2
22
h
b
a S
b y dy b
-=-⎰ (5)
其积分解析表达式为 22221
(arcsin )2
a h S h
b h b b b b π=
-++ (6) 其中， h H b =-
(7)
如图
图2微元法求椭圆切面面积
221
[()(2)arcsin(1)]2a H S H b H b H b b b b π=--+-+
(8) 221
[()(2)arcsin(1)]2
a H V L H
b H b H b b b b π=--+-+
(9)
图3 油罐无倾斜时示意图
（2）当油罐发生纵向偏转时，油罐中油所占空间为一倾斜柱体，如图4所示：
图4 油罐偏移示意图 如图4所示，根据几何关系可知，
'(0.4)tan h H x α=--
(10)
又根据油面的高度不同，可分为以下三种情况：
图5 情况1：低油位
若油面位于图5所示位置，则：
22 10
22
[(
0.4tan tan)(0.4tan)
0.4tan tan1
arcsin]
2
a
V H x b b H b
b
H x b
b b dx
b
αα
ααα
αα
π
=+---+-+
+--
+
⎰（H+0.4tan）/tan
(11)
图6 情况2：正常油位
若油面位于图6所示位置，则：
2.45
22
20
22
[(0.4tan tan)(0.4tan)
0.4tan tan1
arcsin]
2
a
V H x b b H b
b
H x b
b b dx
b
ααα
αα
π
=+---+-+
+--
+
⎰
(12)
图7 情况3高油位
若油面位于图7位置，则：
2.45
30.4(1.2H)/tan 22
[(tan tan 1arcsin ]2
a
V abL x b b
x b b b dx
b απααπ--=---+⎰
(13)
由上述公式知，油罐的变位会对罐内油高与储油量的对应关系（罐容表），产生
较大的影响。

综合式（11）-（13），可以得到模型1如下：
(
)0
222.45
02[(0.4tan tan 0.4tan tan 1arcsin ] H<2.05*tan 2[(0.4tan tan ()arcsin a H x b b
H x b b b dx b a H x b b
V H H b ααααααπααα+--+--+++--=++⎰⎰（H+0.4tan ）/tan ，当
0<22.450.4(1.2H)/tan 22
0.4tan tan 1] (H 1.2-0.4tan )2[(tan tan 1arcsin
] (H 1.2)2x b b dx b a abL x b b x b b b dx b αααπααπααπα--⎧⎪⎪
⎪⎪⎨
--+≤≤---++≤⎰，当2.05*tan ，当1.2-0.4tan <⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎪⎩
(14)
5.1.2 应用试验数据对理论关系式进行修正
当无变位进油时，我们可以根据式（9）
221
[(arcsin(1)]2
a H V L H
b b b b b π=--+
对每一个油位高度求出其理论储油量；另根据累加进油量和罐内油量初值，可求得实际储油量。

由于理论储油量和实测数据之间存在一定的系统误差，所以我们用线性回归方式得到修正系数 m = 1.035。

因此，无变位实际体积的修正计算公式为：
221
[(arcsin(1)]/2f a H V L H b b b m b b π=--+
(15) 对不同高度用式（14）计算对应的体积f V 和实测值进行对比验证，平均误差为0.01%，达到较好的计算精度（图8）。

参考数据见附表1
图8
当罐体有α=4.1°倾斜角的纵向变位时，利用模型1我们对每一实验数据给出的油高计算其理论储油量。

系统误差校正：
所谓系统误差，是由于原始读数不准确造成的，其原因可能是仪表不准确、罐体变形或者进油出油管道和仪表占据一定的容积。

虽然我们不知道具体的原因，但是我们通过统计分析可以一定程度上消除系统误差。

方法如下：根据实验数据中累加进油量和罐内油量初始值求出实际储油量，与模型计算值进行比较，用二阶多项式拟合储油量差值和油高。

20.00040.5834124.24,
0.9773
dv H H Relative =-+-=
(16)
这两列数据的相关系数达到0.977，有理由采用此多项式对模型的计算值进行系统误差修正。

图9 系统误差和油高的拟合
对于试验中变位时，数据中的油高均处在5.1中正常油位情况，模型一中其他两种情况没有涉及。

所以，模型简化为：实际储油量（f
V ）=模型储油量（
2
V ）
-系统误差（dv ），即
2.45
22
022
2[(0.4tan tan )(0.4tan )0.4tan tan 1arcsin ]2
(0.00040.5834124.24)
f a
V H x b b H b b
H x b b b dx
b H H αααααπ=+---+-+--++--+-⎰ (17)
用实验数据验证，拟合效果良好，平均误差为0.059%。

图10
根据模型一，对系统误差进行修正后，我们可以计算求得模型所需的罐容表，详见下表。

小椭圆型储油罐罐容表标定值
油位高度（mm ） 储油量（L ） 油位高度（mm ） 储油量（L ） 油位高度（mm ） 储油量（L ） 油位高度（mm ） 储油量（L ）
0.00 1.67 300.00 580.47 600.00 1716.72 900.00 2995.61 10.00 3.53 310.00 611.97 610.00 1759.00 910.00 3036.59 20.00 6.26 320.00 644.09 620.00 1801.42 920.00 3077.31 30.00 9.97 330.00 676.80 630.00 1843.97 930.00 3117.75 40.00 14.76 340.00 710.08 640.00 1886.63 940.00 3157.90 50.00 20.69 350.00 743.91 650.00 1929.40 950.00 3197.73 60.00 27.85 360.00 778.26 660.00 1972.26 960.00 3237.24 70.00 36.32 370.00 813.12 670.00 2015.20 970.00 3276.39 80.00 46.14 380.00 848.46 680.00 2058.20 980.00 3315.18 90.00
57.39
390.00
884.28
690.00
2101.26
990.00
3353.58
5.2 模型二的建立与求解
5.2.1 建立实际储油罐储油量和油位高度的模型
首先只考虑纵向倾角α。

由于实际储油罐相当于圆柱体与球冠体组成，故用垂直于油罐的平面切割油罐，与罐中的油相交，所截的平面为弓形。

劣弧弓形的面积公式为：
21cos (r h
S r r h r
--=--
(18) 优弧弓形的面积公式为：
21
(cos )(h r
S r h r r
π--=-+-(19) 其中r 为弓形所在圆的半径；h 为弓形的高。

所以罐中油的体积微元为
21[cos (r h
dv r r h dx r
--=--
或者：
21[(cos )(h r
dv r h r dx r
π--=-+-
由于实际储油罐是不规则的几何形体，故我们在计算罐内油的体积时，将卧
式储油罐分为三部分，如图11所示：
图11
由于罐内油的高度不同，可分为以下五种情况：
图12
图13
图14
图15
图16
为了保证罐内油体积的一般性，我们先对图13所示情况进行求解。

则：
V V V V
=++
(20)
ⅠⅡⅢ
对于第二部分的体积的求解，可类比模型一中的方法的第二部分体积的微元为
121
1111(2tan tan )
[cos (2tan tan )
(2tan tan )(22tan tan )]R H x dv R R H x R H x R H x dx
αααααααα--+-=---++---+ (21)
故其体积为：
18
21
1110
1(2tan tan )
[cos (2tan tan )
(2tan tan )(22tan tan )]R H x V R R H x R H x R H x dx
αααααααα--+-=---++---+⎰Ⅱ (22)
对于第一部分体积可用图17求解油高：1H GT =。

图17
图中:
R OA = 1R AI = 2R MN = 1CI =
2221(1)R R R --= ⇒ 1.625R =
2222(1)R R R x =--+ ⇒ 222(1)R R R x =--+1H MT MG =- 1MG R BQ PQ =--
2tan BQ H α=+
tan PQ x α=
所以，
121(2tan tan )H R R H x αα=----
(23)
则： ()21
21
2211212
[cos
()2]R H dv R R H H R H dx R --=---Ⅰ 然后再利用下图确定积分区间[0，1x ]：
图18
图中 CM 即为所要求解的1x 。

CI=1 OI=R-1 IQ=()12tan R H α-+ KI=(R-1)tan α
∴ KN=KQ cos α=()()1[1tan 2tan ]cos R R H ααα-+-+=OP
22R OP -EQ=PE-PN+NQ NQ=KQ sin α
OK=PN=1
cos R α
-
联立以上各式可得：
()()111
cos 1tan 2tan sin }cos R x R R H ααααα
-=+-+-+⎡⎤⎣⎦
(24)
1
21
21
2210
2
[cos
(x R H V R R H dx R --=--⎰Ⅰ (25)
对于第三部分的体积，方法与第一部分的体积求解相似。

其体积为：
(26)
式中：
()()211
cos 1tan 6tan sin }cos R x R R H ααααα
-=-----⎡⎤⎣⎦
(27)
3R =(28) ()2316tan tan H R R H x αα=--++
(29)
所以：
(1
2
21
21
2210
2
8
21
11101
x 21323
320
3[cos ((2tan tan )
[cos (2tan tan )
cos x R H V R R H dx R R H x R R H x R dx R H R R H dx R αααα----=---+-+---++
⎡---⎢⎣⎰⎰⎰
(30)
其他几种情况均可用类似的方法进行求解。

对于图13所示情况：
1
21
21
221022tan 21
1tan 10
1
1[cos ((2tan tan )
[cos (2tan tan x H R H V R R H dx R R H x R R R H x dx
αααααα-+--=---+-+-
--+⎰⎰
(31)
对于图14所示的情况：
()()()133
21122120221111018
121
1
11
cos (2tan tan [cos 2tan tan 2tan tan [cos 2tan tan
x x x H R V R H
R dx R H x R R H x R R R H x R R H x R πααπαααααα---⎡⎛⎫
-=-+-⎢ ⎪⎝⎭⎣⎛⎫+--+-++-- ⎪⎝⎭-+-+---+⎰⎰⎰(
21
x 1
21232
332
20
3cos x dx
R H R R H dx R dx R π-⎡-+--+⎢⎣⎰
⎰ (32)
积分限中出现的表达式
1
32tan tan H R x αα
+-=
(33)
对于图15所示情况，
()(1212112212028
2111101121223
323203cos (2tan tan [cos 2tan tan (cos )x x x H R V R H R dx R H x R R H x R R
H
R R H R dx R dx R πααπααππ---⎡⎛⎫-=-+-⎢ ⎪⎝⎭⎣⎛⎫+--+-++-- ⎪⎝⎭
⎡-+-+-+⎢⎣⎰⎰⎰2
123x R dx π+⎰⎰ (34)
对于图16所示情况：
1
2
1
221
21
212210
2
2tan 21
1tan 10
1
128{[cos ((2tan tan )
[cos (2tan tan }
x H R H V R dx R R R H dx R R H x R R R H x dx ααππαααα-+--=+----+-+-
--+⎰⎰⎰
(35)
因此,我们总结出7个分段被积函数如下:
()21
21
12211212
cos ()2R H fvol R R H H R H R -----(x)= (36)
121
1211
1(2tan tan )
()cos (2tan tan )(2tan tan )(22tan tan )
R H x fvol x R R R H x H x R H x αααααααα--+-=---++---+ (37)
()()21
32
33322323
()cos 2R H fvol x R R H H R H R --=--- (38)
()211242121212()cos ()2H R fvol x R H R H R H R π-⎛⎫
-=-+-- ⎪⎝⎭
(39)
()211511112tan tan ()cos 2tan tan (2tan tan )(22tan tan )
H x R fvol x R H x R R H x R H x ααπαααααα-⎛⎫
+--=-++-- ⎪⎝⎭+---+
(40) ()()21
23
63232323
()(cos )2H R fvol x R H R H R H R π--=-+-- (41)
272()fvol x R π=
(42)
下面再对横向偏移角β进行分析研究，如图19所示：
图19
图中 0h 为罐容表的读数，所以真实液面高度为:
10111cos ()cos cos H R h R R R βββ=+-+-
即：
011()cos H h R R β=-+
(43)
所以最终所得的体积关系式只需将上述体积关系式中的H 换为式(30)，即可得出。

经以上分析,我们得到模型二如下:
()()()()()()()()1
1
2
32132tan tan 1200
81230
001847523000()(), 06tan (), 6tan 1.53tan ()V(,,)H x x x x x x x fvol x dx fvol x dx H fvol x dx fvol x dx fvol x dx H fvol x dx fvol x dx fvol x dx fvol x dx fvol x dx h α
α
ααααβ++≤≤++≤≤-++++=⎰⎰⎰⎰⎰⎰⎰⎰当当()()()()()11
212
1812475630
00,
1.53tan 1.57tan ,
x x
x x x H fvol x dx fvol x dx fvol x dx fvol x dx R dx ααπ-≤≤+++++⎰⎰⎰⎰⎰⎰⎰当()()()12tan 12tan 7112000
1.57tan 32tan 28()(),32tan 3H x H fvol x dx R fvol x dx fvol x dx H α
α
ααπα+⎧⎪⎪⎪⎪⎪⎪⎪
⎪⎨⎪⎪⎪+≤≤-⎪
⎡⎤⎪+-+⎢⎥⎪⎣⎦-≤≤⎩
⎰⎰⎰当 当⎪
⎪
(44)
2
12
arg ,(((,,)(,,)))min (,,)n
i i i i V H V H OilData error OilData αβ
αβαβαβ-==--∑
(45)
式(44)中H 由式（43）确定，其他变量的定义参见符号说明：
使得式（45）取得最小值的
,αβ即是待定的变位参数。

5.2.2 模型的求解
我们利用MATLAB 软件，将模型二中的式（44）和式（45）编程（相关程序见附录），利用非线性规划的方法求解α和β值，α=2.12°，β=4.06°。

根据得到的变位参数值，得到罐体变位后油位间隔为10cm 的罐容表标定值。

011
()cos H h R R β=-+
罐容表标定值
H(mm) V （L ） H （mm ） V(L) 0 46.40662 1500 30232.31 100 354.5192 1600 33046.07 200 1061.754 1700 35855.97 300 2214.65 1800 38646.71 400 3691.607 1900 41402.87 500 5419.525 2000 44108.69 600 7356.6 2100 46747.98 700 9471.698 2200 49303.87 800 11739.27 2300 51758.56 900 14137.18 2400 54092.97 1000 16645.5 2500 56286.15 1100 19245.85 2600 58314.42 1200 21920.95 2700 60149.79 1300 24654.29 2800 61756.69 1400
27429.93
2900
63083
六． 模型的分析
6.1模型的正确性
如图所示，横坐标为题目所作实验的出油量，竖坐标为经过所建模型求解的两次体积差，即模型的出油量，二者相对误差较小，线性拟合较好。

6.2模型的灵敏度分析
当对油高进行0.1%的随机扰动，即让()()'10.001*1,1H H rand =+-时，模型二求得α=2.119°，比较无扰动的α=2.12°，变化幅度也在0.1%左右，说明α的稳定性较好，但是β的值从4.06°变成了3.75°，变化了大约8%。

所以模型在β方向上抗干扰能力较低，因此我们对β进行理论分析。

由011()cos H h R R β=-+可知：
1
1
01
cos H R h R β--=- (46)
所以当01h R 在附近变化时对β的影响较大。

用Matlab 对其进行数据模拟，结果如下图所示：表明β在1R 附近比较敏感。

七． 模型的评价与推广
7.1 模型的优缺点
7.1.1 模型的优点
（1）本文借助高等数学微积分的思想，建立罐体在变位前后标定罐容表的数学模型，得出罐内储油量与油位高度及变位参数的函数关系式，理论基础成熟，可信度较高。

（2）该模型以微积分为基础，简单易懂，又有相应的软件（Matlab 软件）支持，算法简单，容易推广。

（3）运用多种拟合方法使结果更加精确，通过灵敏度和误差分析使模型更具有实际意义，增加了应用价值。

（4）模型二的误差分析运用最小二乘法，使得模型得出的结果更加准确。

7.1.2 模型的缺点
（1）模型二中的β值的稳定性不是很好
（2）在用拟合法处理数据的时候，由于模型假设具有一定的主观性，导致拟合的曲线不是十分精确。

7.2 模型的推广
本文所建立的微分模型不仅适用于储油罐的顶部为球冠的情况，还可以推广到顶部为弧形顶、平顶、椭球顶、锥顶的情况，易于计算，在实际应用中具有延伸和推广的价值。

八． 模型的改进
本文所建立的积分模型中α与β是独立无相关的，而在模型二的求解中用β
的值调整了油高读数,然后再去求解α,因此使得这两个方向的变化独立了,影响了模型的精度。

因此我们必须考虑α方向上对油高读数的影响。

不妨假设
011(cos )cos H h R R αβ=-+，则模型二的形式如下：
()()()()()()()()1
1
2
12tan tan 1200
81230
0014752300arg min ,()(), 06tan (), 6tan 1.53tan ()V(,,)H x x x x x fvol x dx fvol x dx H fvol x dx fvol x dx fvol x dx H fvol x dx fvol x dx fvol x dx fvol x dx fvol x dx h α
α
αβααααβ++≤≤++≤≤-++++=⎰⎰⎰⎰⎰⎰当当()()()()()13231
212
80
01812475630
00,
1.53tan 1.57tan ,
x x x x
x x x H fvol x dx fvol x dx fvol x dx fvol x dx R dx ααπ-≤≤+++++⎰⎰⎰⎰⎰⎰⎰⎰⎰当()()()12tan 12tan 7112000
1.57tan 32tan 28()(),32tan 3H x H fvol x dx R fvol x dx fvol x dx H α
α
ααπα+⎧⎪⎪⎪⎪⎪⎪⎨+≤≤-⎡⎤+-+⎢⎥⎣⎦-≤≤⎰⎰⎰当 当⎪
⎪⎪⎪⎪⎪
⎪⎪⎪
⎪⎩
(47)
式中：
011(cos )cos H h R R αβ=-+
21
2
arg ,(((,,)(,,)))min (,,)n
i
i i i V H V H
OilData error OilData αβ
αβαβαβ-==--∑
用matlab 软件搜索求得：
α=2.0925° β=3.5055°
误差的平方和为2.4574，较前个模型二的误差平方和2.4570大了一点点，但这微小的差别几乎可以忽略不计，但是我们得到的α与β更小，从实际意义来说，更加符合实际，我们认为这个两个值更加准确。

当然还有很多的假设，他们都满足一定的数学条件，有自己的物理意义。

如可以假设011[(1tan )]cos H h R R αβ=--+，或者011()cos cos H h R R βα=-+等等。

但是α与β真正意义上的关系式还需要很复杂的数学推导和耐心的去解决，在此由于时间有限，我们不予解答出来。

等以后有时间再仔细推敲。

模型二数据的误差检校，由于模型二以实际卧式储油罐为研究对象，则不能忽略注油管和出油罐等管道对罐内油体积的影响。

因此我们必须对模型二的数据进行检校，排除系统误差的影响。

罐内油体积的相对误差为17%左右，用油高读数和实际油的体积与理论算出来的体积之差拟合，如下图所示：
其相对误差大约为0.6%，较前面的相对误差减小了100多倍，数据更加精确。

使模型更加优化。

九．参考文献
【1】姜启源，数学模型（第三版）.高等教育出版社，2003
【2】石辛民，基于MATLAB的实用数值计算.清华大学出版社，2006
【3】黄永建，MATLAB语言在运筹学中的应用，湖南大学出版社，2005
附录
附表一
无变位进油无变位进油
实际油量修正后油量相对误差实际油量修正后油量相对误差
312311.96390.00011617621761.7960.000116
362361.96439.87E-0518121811.7810.000121
412411.94680.00012918621861.7920.000112
462461.96318E-0519121911.7670.000122
512511.93420.00012919621961.7720.000116
562561.93790.0001120122011.7850.000107
612611.93430.00010720622061.7890.000102
662661.91410.0001321122111.7620.000113
712711.92929.94E-0521622161.7710.000106
MATLAB 程序
function [volumes]= totalEval(x)
tic;
oilTank = getOil;
global H alpha
alpha = x(1);
beta = x(2);
n = size(oilTank, 1);
volumes = zeros(n, 1);
for i = 1:31
H0 = (i - 1) / 10;
H = (H0 - 1.5) * cos(beta) + 1.5;
[x1 x2] = calcL1(H, alpha);
if H < 6 * tan(alpha)
volumes(i) = real(quad(@intV2, 0, x1) + quad(@intV1, 0, (H + 2 * tan(alpha))/ tan(alpha)));
elseif H < 1.5 - 3 * tan(alpha)
volumes(i) = real(quad(@intV1, 0,8) + quad(@intV2, 0, x1) + quad(@intV3, 0, x2));
elseif H < 1.5 + 7 * tan(alpha)
lmax = (H + 2 * tan(alpha) - 1.5) / tan(alpha);
volumes(i) = real(quad(@intV4, 0, 1) + quad(@intV5, 0, lmax) + quad(@intV1, lmax, 8) + quad(@intV3, 0, x2));
elseif H < 3 - 2 * tan(alpha)
volumes(i) = real(quad(@intV4, 0, 1) + quad(@intV1, 0,8) + quad(@intV6, 0,1));
else
volumes(i) = real(quad(@intV7, 0,1)) + 8 * pi * 2.25 - real(quad(@intV2, 0, x1) + quad(@intV1, 0, (H + 2 * tan(alpha))/ tan(alpha)));
end
end
toc;
function [xmax1 xmax2] = calcL1(H, alpha)
R = 1.625;
R1 = 1.5;
talpha = tan(alpha);
Rc = (R - 1) ./ cos(alpha);
xmax1 = ((R ^ 2 - (( R - 1) * talpha + R1 - (H + 2 * talpha)) .^ 2 * cos(alpha) ^2) .^ 0.5 - Rc + ...
((R -1) * talpha + R1 - (H + 2 * talpha)) * sin(alpha)) * cos(alpha);
xmax2 = ((R ^ 2 - ( -( R - 1) * talpha + R1 - (H - 6 * talpha)) .^ 2 * cos(alpha) ^2) .^ 0.5 - Rc - ...
( -(R -1) * talpha + R1 - (H - 6 * talpha)) * sin(alpha)) * cos(alpha);
%function V1 = intV1(H, x, alpha)
function V1 = intV1(x)
R1 = 1.5;
global H alpha
%H = 1.29;
%alpha = 0.07;
HH = H + 2 * tan(alpha);
V1 = R1 ^ 2 .* acos( (R1 - ( HH - x .* tan(alpha))) / R1) - ...
(R1 - HH + x .* tan(alpha)) .* ((HH - x .* tan(alpha)) .* (2 * R1 - HH + x .* tan(alpha))) .^ 0.5;
function V = intV2(x)
global H alpha
R = 1.625;
R1 = 1.5;
R2 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
%R3 = R2;
talpha = tan(alpha);
H1 = R2 - ( R1 - H - 2 * talpha - x * talpha);
%H2 = R3 - ( R1 - H + 6 * talpha + x .* talpha);
V = R2 .^ 2 .* acos((R2 - H1) ./ R2) - (R2 - H1) .* ( H1 .* ( 2 .* R2 - H1)) .^ 0.5;
function V = intV3(x)
global H alpha
R = 1.625;
R1 = 1.5;
R3 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
talpha = tan(alpha);
H2 = R3 - ( R1 - H + 6 * talpha + x .* talpha);
V = R3 .^ 2 .* acos((R3 - H2) ./ R3) - (R3 - H2) .* ( H2 .* ( 2 .* R3 - H2)) .^ 0.5;
function V = intV4(x)
global H alpha
R = 1.625;
R1 = 1.5;
R2 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
%R3 = R2;
H1 = R2 + H + 2 * tan(alpha) + x * tan(alpha) - R1;
%H2 = R3 - ( R1 - H + 6 * talpha + x .* talpha);
V = R2 .^ 2 .* ( pi - acos(-(R2 - H1) ./ R2)) - (R2 - H1) .* ( H1 .* ( 2 .* R2 - H1)) .^ 0.5;
function V = intV5(x)
global H alpha
R = 1.625;
R1 = 1.5;
R2 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
%R3 = R2;
%H1 = R2 + H + 2 * tan(alpha) + x * tan(alpha) - R1;
%H2 = R3 - ( R1 - H + 6 * talpha + x .* talpha);
V = R1 .^ 2 .* ( pi - acos( (H + 2 .* tan(alpha) - x .* tan(alpha) - R1) / R1)) + ...
(H + 2 * tan(alpha) - x .* tan(alpha) - R1) .* ...
((H + 2 * tan(alpha) - x .* tan(alpha)) .* (2 * R1 - H - 2 * tan(alpha) + x .* tan(alpha))) .^ 0.5;
function V = intV6(x)
global H alpha
R = 1.625;
R1 = 1.5;
R3 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
%R3 = R2;
%H1 = R2 + H + 2 * tan(alpha) + x * tan(alpha) - R1;
talpha = tan(alpha);
H2 = R3 - ( R1 - H + 6 * talpha + x .* talpha);
V = R3 .^ 2 .* (pi - acos((H2 - R3) ./ R3)) + (H2 - R3) .* (H2 .* ( 2 * R3 - H2)) .^ 0.5;
function V = intV7(x)
global H alpha
R = 1.625;
R1 = 1.5;
R2 = (R^2 - (R - 1 + x) .^ 2) .^ 0.5;
%R3 = R2;
%H1 = R2 + H + 2 * tan(alpha) + x * tan(alpha) - R1;
V = pi * R2 .^ 2;。