清华大学论文模板

来源：个人技术集锦

题目：储油罐的变位识别与罐容表标定

【摘要】

本文根据题目要求，针对具体问题建立相应的数学模型，并对问题进行精确求解，逐层深入比较题给数据与理论求解值同时分析它们之间的关系，据此得出合理的结论。

对于问题一,首先，采用截面法建立数学模型对无变位及有变位的理论油容量与油位高度之间的关系表达式。然后，根据附件1所给的无变位进油时油高求出对应理论油量并与所给实际油量值相比较，发现两者的相对偏差基本成一稳定值3.371%，并根据两者之间的关系求出理论油量的修正函数，同时将偏差函数运用于无变位出油中，通过比较相邻两时刻绝对出油量的相对误差（均值为3.21104）检验修正函数的正确性。接着，根据对应倾斜变位进油油高的校正理论油量与实际油量相比较，分析得出纵向变位时液位较小时罐容表读数偏小，液位较高时读数偏大的结论。最后运用MATLAB编写程序对变位后罐容按油位高度间隔为1cm的进行标定，结果见附表1。

对于问题二，由于横向倾斜不影响体积的计算，只影响油面高度实际值h与罐容表测量值h*之间的关系。于是，先根据三重积分与截面法计算只有纵向倾斜的情况下罐容与倾斜角度和h的关系Vyf(,h)。然后，根据油面高度与横向偏转角度及h*关系hf(,h*)得罐内储油量与油位高度及纵向倾斜角度和横向偏转角度之间的一般关系Vyf(,,h*)。取0,0（无变位时），根据附件2中油位高度计算出相应油量，发现与附件所给油量容积几乎完全一致。接着利用附件2中实际出油量与理论出油量的差值的平方和最小的条件，运用遍历法和NDEPSO算法两种方法分别求得结果如下： α（度） β（度）误差遍历法 2.1 4.4 4.9e-5 NDEPSO算法 2.129 4.244 3.7e-6 根据误差最小原则选取(2.129,4.244)，求得罐容表标定值部分数据如下（其余数据见附表4）： h(m) v(L) 0 45.88 0.1 355.42 0.2 1068.68 0.3 2228.22 0.4 0.5 0.6 3709.18 5440.16 7379.69 利用第二次注油后的数据，带入已求出的α、β，和表中的h求出实际出油量与理论出油量近似相等，验证模型的正确性和可靠性。并进行了模型推广，进行了基于MATLAB的GUI设计。

关键词：截面法，三重积分，遍历法，NDEPSO算法，GUI设计

一问题的重述

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

许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。

请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。

（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。

（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度 ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。

二模型的假设

1、假设储油罐纵向偏斜时偏斜角较小，符合地基变化的实际情况。 2、假设忽略储油罐偏斜时浮子测量所带来的误差。 3、假设忽略出油与进油速度对油位高度的影响 4、假设图1中的储油罐尺寸上不存在工艺误差

三符号约定

St -----截椭圆面面积 Vt -----截椭圆柱体体积

------纵向倾斜角度

Sy------截圆面面积

Hg-----球冠高度

Vy1-----左侧球冠处油体积 Vy2-----中间圆柱体处油体积

Vy3-----右侧球冠处油体积

------横向偏转角度 Vy1z-----整个球冠体积

四模型的分析与建立

4.1对问题一的分析与模型建立 4.1.1问题一的分析

为了掌握罐体变位后对罐容表的影响，首先对无变位的情况进行分析，通过建立无变位体积求解模型，求得理论值。与附件一中实验测得值进行比较，分析得出两者之间的关系，然后据此对理论值进行修正，同利用无变位出油数据对其进行检验。

截椭圆柱体体积与油位高度之间的关系：采用截面法建立数学模型对无变位及有变位的理论油容量进行计算，得出了两端平头的椭圆柱体油罐的罐容与油位高度之间的关系。

另外，以所得的偏差为依据对所计算有变位时的理论油容量进行校正，再将其与试验测得数据进行比较，分析罐体变位后对罐容表的影响并对给出罐体变位后油位高度间隔为1cm的罐容表标定值。最后，运用曲线拟合的方法来形象直观地说明理论值、校正值以及实验测试值之间的关系。 4.1.2问题一的模型建立

4.1.2.1椭圆截面面积公式的推导：

通过对问题一的分析，我们建立求解截椭圆面面积的模型：

设椭圆的长短轴长分别为2a、2b；截面的油位高度为H。

图1 截面椭圆示意图

2za22xzbz 设椭圆方程为：221；xa12bbab22被截椭圆面的液面高为H,其中0H2b，图中带阴影部分为储油横截面，先用微元法求解微元面积ds：

HbHbHbds2xdz；则：Stbds2xdzbba222bzdz b采用换元法，令zsin，则dzcosd

z与积分区间变换关系如下表：

表1 积分区间变化关系 z1 b Hb 221 2 arcsin((Hb)/b) z2 22 则：St2abcosdab(cos(2)1)d

11积分得：Stab2121absin(2) 21带上下限可得：StHba(Hb)H(2bH) ababarcsin()2bb既得截面面积St与截面油高H之间的关系。

4.1.2.2储油罐无变位时体积公式推导：

储油罐无变位时，当油位高度一定时，储油罐的横截面是一个面积恒定的截椭圆，故此时的体积即为等截面柱体体积的求解。

设柱体高为L，其中Ll1l2 此时的横截面图形如图1所示：

SHba(Hb)H(2bH) ababarcsin()2bb可得无变位时油罐储油量与油位高度之间的关系，如下式所示：

a(Hb)H(2bH)HbVtSLL(ababarcsin())

2bb4.1.2.3储油罐纵向变位时体积公式推导：

当储油罐纵向变位时，其中储油横截面为变面积截面，故考虑使用截面法采用积分来对体积进行求解，具体推导过程如下图所示：

图2 平头储油罐纵切面示意图

其中h为油位高度，H为对应每个y时所截截面的油面高度；

根据上述图形所建立的坐标系，沿y轴对截面S进行积分既得体积，即：

Hba(Hb)H(2bH)VtSdy(ababarcsin())dy

2bby1y1其中HH(y)

根据y的积分限将积分区间分成三段进行积分，如下所示：

当油位高度0hl2tan（即，hAB）时，罐内油可看成一个截面面积一端为0的截椭圆柱体；

积分下限y10，积分上限y2l1htan

y2y2Hh(l1y)tan

l1htanVt10(ababarcsin(2Hba(Hb)H(2bH)))dy bb当油位高度l2tanh2bl1tana（即，ABhAC）时。积分下限y10，积分上限y2l1l2

Hh(l1y)tan

l1l2Vt20(ababarcsin(2Hba(Hb)H(2bH)))dy bb当油位高度2bl1tanah2b（即，ACh）时，储油罐内的油体积可以看做整个椭圆柱体的体积减去上面一部分一个一端为0的小截椭圆柱体所得

的部分；

积分下限y1l1(2bh)/tana，积分上限y2l1l2

VtzabL

H2bh(l1x)tan

Vt3VtzHba(Hb)H(2bH)(ababarcsin())dy2bbl1(2bh)/tanal1l2 综上所述,建立椭圆油罐体油容积与油位高度h 的关系式为：Vtf(h)

4.2问题二的分析与模型建立 4.2.1问题二的分析

对于问题二，它既有纵向的倾斜又有横向偏转，首先，考虑只有纵向倾斜的情况下得出罐容与倾斜角度和油位高度的关系。然后，由于在纵向倾斜角一定的情况下，横向偏斜不影响油体积的求取，其实际油位高度与横向偏斜与否无关，根据实际油位高度与油浮子测得的高度及其横向偏转角度之间的关系，既可得，罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度 ）之间的一般关系。通过分析附件2的数据我们看出给出的显示油高和显示容积应该为无偏向时的罐容表，我们就可以利用先前得到的一般关系式，令两偏斜角度为零时，即无偏斜的情况下，计算相应油高时的理论容积，再与附件2中的显示容积比较，得出两者偏差，并可达到间接检验设立的模型的正确性的目的。再根据模型，将一定的偏斜角度带入，根据对应的油高，即可得对应的容积，于是我们得到两对应高度下之差的理论容积之差，再与附件2的对应的高度下的实际出油量（容积之差）相比，若两偏斜角度选取合适，则相比之差应很小，于是我们可以采用某种算法，寻找这样的两偏斜角度，使得所有的理论容积之差与实际出油量偏差很小，故我们可以计算多个这样的偏差平方和，将问题转化为求解使得该偏差平方和最小的两偏斜角度。再通过这两偏斜角度带入附件2中数据，如上所示的求解所谓的平方和，即这两偏斜角度对应下的总偏差，检验其值，即可检验模型的正确性和可靠性。 4.2.2问题二的模型建立

4.2.2.1圆截面面积公式的推导：

通过对问题二的分析，我们建立求解截圆面面积的模型：设圆半径为R；截面的油位高度为R。

图3 截圆面示意图

设圆方程为：x2z2R2

根据椭圆与圆之间的对称关系根据椭圆截面面积St与油位高度H之间的关系：

Hba(Hb)H(2bH) Stababarcsin()2bb将椭圆的长短轴a、b均用圆的半径R进行替换，可得圆截面面积Sy与油位高度H之间的关系： HRSyR2R2arcsin()(HR)H(2RH)

2R4.2.2.2体积公式的推导：

现要对一个主体为圆柱体，两端为球冠体的储油罐储油体积进行计算。主体圆柱体的半径为R；圆柱体高为L；两端球冠高为Hg；球冠对应球体半径为r；油浮子距圆柱体两端的距离分别为l1、l2；罐体纵向倾斜角度为。

此时油容积由三部分组成，分别为左侧球冠处油体积Vy1，中间圆柱体处油体积Vy2，右侧球冠处油体积Vy3。现分别对这三部分体积公式进行推导。 4.2.2.3罐体纵向变位时体积公式推导

先只考虑罐体只有纵向变位时的体积的推导：（1）、左侧球冠处油体积公式推导

对左侧球冠处油体积采用三重积分进行计算，由于左侧球冠的球心与与其相交的圆柱体圆的圆心不共心，现对球冠所在球的半径进行求解：

图4 球冠半径求取示意图

如图所示可得r满足方程：OC2AC2(OBAB)2 即：rR(rHg)，解得：r222(R22Hg)2Hg

现建立空间直角坐标系如下图所示：以球冠球心为坐标原点，垂直于圆柱高向上方向为Z轴，平行于圆柱高方向指向球冠为Y轴，垂直圆柱高向里方向为X轴。

图5 纵向倾斜储油罐坐标系建立示意图

由以球心为坐标原点，得：球面方程为：x2y2z2r2

球冠与圆柱体的交面平行于XOZ平面（x0，z0），距坐标原点的距离

drHg，故交面方程为：yrHg

现将过球心的纵切面绘制如下：

图6 过球心纵切面示意图

据图分析可知，油平面平行与X轴，油平面法线为图中n所示，n垂直于平面且与Z轴成角，则平面的法向量为（0，-sin,cos）,从图可知平面

x00过点N(x0,y0,z0)，其中：y0rHg

z0(R(hl1tan))据平面的点法式方程可知，油平面方程为：

sin(yy0)cos(zz0)0

积分次序为：先Z后X最后Y；根据Z的积分限将积分分成两部分，如下所述：

当油平面与球冠的交线在球冠左端点以下(即B点以下)时，油平面以下，交面以左与球面所围成的球冠面处体积即为所求，积分区域为：

x2y2z2r2yrHg

sin(yy0)cos(zz0)0y2x2z2Vy11dydxdzdydxdz

y1x1z1其中Z轴方向的积分限下限为球面，上限为油平面；

由球面方程x2y2z2r2解出可得Z轴方向的积分下限z1，通过油平面方程解出Z的表达式既得Z轴方向的积分上限z2：

zr2x2y21 sin(yy0)z0z2cosX轴方向的积分限上、下限为球面与油平面的交线在XOY平面的投影，由

x2y2z2r2球面方程与油平面联立如下所示：

sin(yy0)cos(zz0)0由上述联立方程消去Z即得在XOY平面的投影方程求解可得X轴方向的积

xr2y2z221分下限x1，上限x2：

222x2ryz2Y轴方向的积分限下限为球冠与圆柱体的交面所对应的Y坐标,即

y1rHg，积分限上限为过球心的纵切面、球面、油平面的交点对应的Y坐标：

图7 Y轴方向积分限求取示意图

由图可知直线方程斜率为tan，且过点N(y0,z0)，y2为直线方程与圆方

y2z2r2程的交点：，联立方程求解可得：

ztan(yy)z0022(1tan2)y2(2z0tan2x0tan2)y(x0tan2z02z0x0tanr2)0BB24AC 即，AyByC0，解得：y2，又有y1rHg

2A2综上将三重积分化为一重积分可得：

y2x2z2y2Vy11dydxdz(z2(x2x1)y1x1z1y1121(ry2)[(21)(sin22sin21)]dy22221arcsin(x1/ry) 其中：

222arcsin(x2/ry)当油平面与球冠的交线在球冠左端点以上(即B点以上)时，油平面以下，交面以左与球面所围成的球冠面处体积即为整个球冠的体积减去积分区域为部分的体积：

12(rHg) 整个球冠的体积Vy1zHg3x2y2z2r2积分区域yrHg

sin(yy0)cos(zz0)0y2x2z3Vy12Vy1zdydxdz（注：与第一种情况相比被积函数区别只在于Z轴

y1x1z2的积分上、下限分别为z2、z3），化三重积分为一重积分得：

Vy122111H(rHg)(z2(x2x1)(r2y2)[(21)(sin22sin21)]dy322y1y2g（注：与第一种情况相比被积函数区别只在于z2前的系数为负）。

其中：z3r2x2y2，其他积分上、下限同第一种情况。

由于积分限y1、y2为油位高度h跟纵向偏角的函数，故得到了左侧球冠油容积Vy1f1(,h)。综上两种情况可得：当油位高度h在（0，2R）之间变化

Vy11时，左球冠容积表达式为：Vy1f1(,h)Vy12;;h(0,R(Hgl1)tan)h(R(Hgl1)tan,2R)

（2）、中间圆柱体处油体积公式推导

由已求得圆截面面积Sy与油位高度H之间的关系：

Sy2R2R2arcsin(HR)(HR)H(2RH) R

图8 积分限分段示意图

同样与截椭圆柱体体积公式计算雷同（在此不做详细推导），跟据其对称关

系，可得出：Vy2Vy21;h(0,l2tan)f2(,h)Vy22;h(l2tan,2Rl1tan)

Vy23;h(2Rl1tan,2R)其中：

H1Vy21Sydy;H1l1h/tan;Hh(l1y)tan0l1l2Vy22Sydy;l1l2L;Hh(l1y)tan0L12Vy23RLSydy;H2l1(2Rh)tan;H2Rh(l1y)H2（3）、右侧球冠处油体积公式推导：

与左侧球冠体积求法类似，将坐标原点建立于球冠对应球的球心。垂直于圆柱高向上方向为Z轴，平行于圆柱高方向指向球冠为Y轴，垂直圆柱高向外方向为X轴。

得球面方程、球冠与圆柱面交面方程、油平面方程如下:

x2y2z2r2; yrHg;sin(yy0)cos(zz0)0其中：y0rHg ；z0(R(hl2tan)) 。

Vy3Vy31;h(0,l2tan)f3(,h)Vy32;h(l2tan,R(l2Hg)tan)

Vy33;h(R(l2Hg)tan,2R)其中：

Vy310y21212Vy32(z2(x2x1)(ry)[(21)(sin22sin21)]22y1y21VH2(rH)(z(xx)1(r2y2)[()1(sin2sin2)]dygg21212212y3332y1 其中：z2积表达式。

sin(yy0)z0，其余积分上、下限表达式同左侧球冠体

cos由于纵向偏斜时，左侧球冠处油体积Vy1，中间圆柱体处油体积Vy2，右侧球冠处油体积Vy3。三部分体积表达式均为f(,h)以 、h为变量且按照h分段的表达式，而与坐标系的建立无关。故可得主体为圆柱体，两端为球冠体的储油罐纵向倾斜时，油容积Vy与油位高度h之间的关系为三部分体积表达式的和，即：Vy=f(,h)=Vy1+Vy2+Vy3

4.2.2.4罐体纵向变位加上横向偏斜时体积公式推导

由于当罐体纵向偏斜角时，在此基础上横向偏斜，此时圆截面所对应的液面高度保持不变，故横向偏斜以后油实际容积不发生改变，改变的只是油浮子与测的的油位高度发生变化，油浮子从原来跟油面垂直变换到与油面成()角，从而导致油浮子测得的油位高度与实际油位高度发生了偏差，而在2 11

上步罐体纵向偏斜时体积求解是油浮子所测的油位高度的函数，当时油浮子测量高度与实际高度相等。因此，此时，只要将油面实际高度h用罐体横向偏斜时油浮子所测得的油位高度h*来表示，即，hf(,h*)，而Vy=f(,h)

*hf(,h)联立  ，可得Vyf(,,h*)

Vyf(,h)油面实际高度h用罐体横向偏斜时油浮子所测得的油位高度h两者之间的关系用下图表示：

图9 横向偏转时液面高度变化关系示意图

由图可知：h与h*均过圆心，油面高度高于圆心或者低于圆心，均可得：

h(h*R)cosR 带入表达式

VyVyf(,,h*)f(,h)=，便可得，既得罐内储油量与油位

高度以变位参数（纵向倾斜角及横向偏转角之间的）一般关系。

五模型的求解

5.1模型一的求解

根据问题一可知，两端平头的圆柱体具体参数为：小椭圆油罐横截面椭圆长轴2a1.78，短轴2b1.2，l10.4,l22.05,油罐总长Ll1l22.45。 5.1.1确定储油罐实际储油量与理论储油量之间的偏差

根据对问题一所建立的模型，罐体无变位时油罐储油量与油位高度之间的关系为：

Hba(Hb)H(2bH)VtSLL(ababarcsin())

2bb其中H 即为油位高度，将附件1中所给的无变位进油累加进油量与罐内油

量初值相加，得到各个油位高度对应的油量。对应附件1中无变位进油的油位

高度根据上式运用Excell，求出对应进油高度所计算的理论油量，将理论油量

VT与附件中所给的实际油量VS进行比较分析，得出两者之间的相对偏差

VTVS，从数据（见附表2）可以看出相对偏差很稳定，偏差的方差很小，VT故通过取平均得出平均相对偏差3.71%。

偏差分析：从实测数据可知理论油量大于实际油量，据分析可知，造成这一现象主要是因为题给的是一个简化模型，忽略了壁厚将外径等同于内径，导致带来较大的偏差。同时由于储油罐内的油浮子、管道等也占据一定的空间（这部分只占很少的一部分）造成理论值大于实际值。

偏差修正：根据理论油量VT与实际油量VS之间的偏差关系可得油量的修正函数为：VS(1)VT。

偏差校正检验（见程序2），对应附件附件1中无变位出油的油位高度求出与之对应的理论油容量，并按照上述修正函数进行修正后，将两两相邻的油容量相减即可得出两个时间点之间的绝对出油量VT与实际所给绝对出油量VS（相邻两累计出油量之间的差值）取差值，求其相对误差VTVS，VTVTVSVT，再

对相对误差求和取平均得3.21104，故验证了所求得的修正函数关系是正确的。

5.1.2罐体变位后对罐容表影响的研究

根据对问题一所建立的模型，罐体无变位时油罐储油量与油位高度之间的关系：

Vt1Vtf(h)Vt2 (详见4.1.2问题一模型建立)

Vt3根据上述罐体变位后油容积与油位高度h之间的函数关系，对应附件1中倾斜变位进油中各个油位高度所对应的倾斜变位时的理论油量进行求解，同样运用5.1.1中的偏差校正函数进行校正得出各油位高度对应的油容量实际值进行比较分析，计算出两者之间的现对偏差（见附表3）。从数据可以看出随着油位高度从411.29mm到1035.36mm逐渐增加，两者之间的偏差满足从一个正值不断减少到零，再继续递减到一个负值的趋势。

变位对罐容表影响分析：从数据分析可知在液位低于一定值时理论值大于实际显示值，并满足随着液位降低两者之间的偏差逐渐增大；当液位高于一定值时理论值小于实际显示值，且满足随着液位高度的增加两者之间的偏差呈增大趋势。这是由于储油罐发生了纵向倾斜而罐容表的读数没有进行相应的修正，

这就导致当油位高度h=0时，实际上是存在油容量的而罐容表标注误以为油容量为0；当h达到最大时，实际油并没有满罐，而罐容表误以为油满，故出现h较小时理论值大于显示值，当h较大时理论值小于显示值，符合两者偏差的变化规律。故可得纵向变位时液位较小时读数偏小，液位较高时读数偏大。 5.1.3罐体变位后罐容表标定

罐体按照倾斜角为=4.10的纵向变位后，当h从0到1.2之间按1cm间隔递增时，根据模型一所建立的变位后容积与油位高度的函数关系，运用MATLAB运算出其对应的油容积，并对其进行标注。标注结果见（附表1）。 5.2问题二的求解

问题要去对一个主体为圆柱体，两端为球冠体的储油罐储油体积进行计算。主体圆柱体的半径为R=1.5；圆柱体高为L=8；两端球冠高为Hg1；球冠对应球体半径为r；油浮子距圆柱体两端的距离分别为l12、l26；罐体纵向倾斜角度为，横向偏转角为。

根据球冠对应球体半径满足r2R2(rHg)2关系式，解得r=1.625。根据4.2.2建立的模型，所得主体为圆柱体，两端为球冠体的储油罐储油体积与纵向倾斜角度为，横向偏转角为以及油面高度h之间的关系满足：。 Vyf(,,h*)（详见问题二模型建立）

对应于附件2中显示各个油高，运用MATLAB（见程序4）求出其对应的油量容积，与附件2中的数据进行比较发现与显示油容量容积几乎完全一致，可以准确到小数点后两位，据此可以验证我们所建立的罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度 ）之间的一般关系模型是正确的。

得到

*Vyf(,,h*)后，我们利用附件二两次注油间的数据，采用用两种

方法分别计算出了、。利用第二次注油后的数据检验了α、β的准确性。

计算、：

方法一：采用遍历法穷举、。根据JJG 266-1996卧式金属罐容积检定规程，使、分布在在可能出现的极大范围内，故认为其分别在(0°，10°）间变化,步长均为0.1°。对于每组确定的、β再由附件2中两个相邻的高h﹡,可求得相邻的Vy之差，和表中给出的相邻V之差比较，越使两者接近的一组、约可靠，故可转化为求fmin(VVy)最小值问题。附件二中，两次注油间的数据约

i1n2三百组，从而找出最可靠的、，计算结果见下表。计算程序见（程序6）。

方法二：采用改进差分进化算法的基础上与粒子群算法进行融合，在方法一的基础上求解最小值。其融合寻优策略如下：

PSO算法与DE算法都是基于群体的启发式算法，它们的最主要区别是新个体的产生方式。PSO算法在优化求解过程中前期收敛速度较快，后期易陷入局部最优。差分进化算法在选择操作中采用的是一种“贪婪”搜索策略，即经过变异和交叉操作后的个体与父代个体进行竞争，只有当其适应度好于父代时才被选作为子代，否则直接进入下一代，该机制可以增加算法的收敛速度，但也容易陷入局部最优点而使算法停滞[5]。

针对DE和PSO单个算法均异陷入局部最优的缺点，本文考虑将两种算法进行融合。由于两种算法存在差异，故一种算法的局部最优可能对另一种算法不适应，即两者均有能力跳出对方的局部最优解，通过两种算法的结合便可以很好的跳出局部最优去搜索全局最优解。同样两者之间信息的交流，给双方一个扰动，有利于跳出局部最优。下面通过图示来说明其中的关系：

全局最优当前最优局部最优全局最优当前最优局部最优来自粒子群的个体

图10 DE模型中个体运动方式图11 NDEPSO模型中个体运动方式图10给出了在基本NDE模型中个体运动方式,随着迭代次数的增加种群的平均适应度方差不断减小，即所有个体都趋同于群体中最优个体（gebest）,若此时的最优个体为非全局最优，同时所给扰动不够，则会陷入局部最优而无法逃逸。若将 NDE跟PSO进行融合，则个体运动方式如图11所示。若陷入局部最优，但是存在来自于PSO的信息交流，此时则会引导个体逃逸原来局部最优点，不断去搜寻全局最优解，从而实现精确求解。

NDE和PSO算法融合策略：将随机产生种群按适应度升序排列后分成大小相等的较小规模种群，一半适应度排序较劣称为NDE子种群，另一半适应度排序较优称为PSO种群；NDE子种群按照NDE进行进化，PSO子种群按照PSO进行进化，两个子种群各自进化形成的新子种群相组合形成我们的下一代种群，然后继续按照上述步骤反复操作，直至寻找到最优解。算法流程如图12所示：

第i个个体第i个个体 15

初始化各参数形成初始种群令T=1计算适应度进化能力排序选择一半当前种群形成DE子种群选择另一半种群形成PSO子种群计算当前代数缩放因子进行变异操作对PSO子种群进行速度更新进入选择池并计算适应度进行串行选择计算当前代数交叉概率因子进行交叉操作对PSO子种群进行位置更新进入选择池并计算适应度进行串行选择混合形成新种群判断t是否大于Tmax得出最优解算法结束

图12 NDEPSO算法流程

由于两种算法结合带来计算的时间复杂度的增加，为了克服这一问题，本文在NDEPSO算法上继续改进。通过改进选择策略，提高收敛速度。改进选择策略为：将父代与试验个体一起入“选择池”按适应度排序，适应度优的留下，称之为串行选择。串行选择的优点在于若父代与对应试验个体在适应度排序中都较优则都被留下来，而在基本算法中父代与试验个体只能选其一，这样加速了进化速度，加快了收敛速率。

参数设置：函数维数为2，种群规模设为20，最大迭代次数100，DE的缩放因子

F[0.4,0.8]，交叉概率因子CR[0.3,0.9]，学习因子均为1.4962，惯性权重为

0.7298。种群中每个个体都含有变量（、），经过优化寻优，可得结果如下表所示。计算程序见(程序7、8)，响应曲线(见附图8) α β 误差方法一 2.1 4.4 4.9e-5 方法二 2.129 4.244 3.7e-6 由于方法一采用步长较大，显然方法二较有优势。故采用方法二数据，认定α=2.219、β=4.244。将所得的（α、β）带入Vyf(,,h*)中便较容易得出罐容比，计算程序见（程序5）。结果见（附表4）。

模型二可靠性分析：利用附件二第二次注油后的数据对α、β值的可靠性进行分析。由计算出的α、β以及标中每个油标高h带入Vyf(,,h*)，相邻

Vy差值与表中对应相邻v的差值比较，VTVS。求解程序见（程序12）VT=1.5*e-4。

六模型的评价

6.1 模型优点

本文分别采用了截面法、三重积分法等基本原理，使得计算值更为精确，较好的阐明了罐内储油量、油位高度以及变位参数这一动态平衡数学模型，同时研究了一种混合智能计算方法NDEPSO，使得计算机能够精确地寻找变位参数，并控制误差在极低的范围，以实现储油罐的准确变位识别，并能根据所得到的变位参数，高精度地进行罐容表标定。

同时根据本命题本文进行了基于MATLAB储油罐罐容表标定GUI设计，图形开发界面充分发挥了MATLAB的强大计算功能，使得用户可以在直观简洁的操作界面上，只需输入相关数据，按按鼠标就可以完成复杂的计算，并可以得到可视化的计算结果。所设计的模型系统具有：

1）适用性强，涉及范围广，针对不同规格的储油罐，可灵活根据参数进行罐容表，且操作方便、简单、明确。

2）能针对不同的进/出油及对应油位高度之差，进行储油罐变位识别，且计算精确、可靠性高。 6.2 模型缺点及改进

任何模型、系统都受到实际生活中的各种限制，本模型也不例外，为了简化模型，基本假设很多都是理想状态。实际上由于受到外界温度，以及进\\出油所产生的压力变化和流速等因素的影响，油的密度等自身属性将发生轻微变化，造成系统与实际值产生偏差。同时计算机计算的精度同样会影响到最后的结果值。此外，本模型在运行及GUI设计方面还存在缺陷，算法复杂度较高，运行所需时间较长，这些都有待于以后再作进一步研究。

七模型的推广

本文中所建立的模型不仅对于解决题目中的问题具有很强的实用性，而且具有很好的推广性。例如，本文提出了一种NDEPSO算法，该算法结合了改进的DE和PSO算法的优良性能，在解决局部最优问题上得到了很好的应用，寻优能力较强，并且选择操上采用了串行选择策略，加速了进化速度，同遍历算法进行比较，在寻找最优解的过程中具有一定的优越性。

由于本文是在基于大量的MATLAB实验获得的，程序较多，运行操作较复杂。本文创造性的进行了基于MATLAB的GUI设计，使得用户可以进行直观的操作，并得到准确的罐容表。经设计后的GUI设计，只需在MATLAB命令窗口中输入Question即可进入选择系统界面，经选择可进入相应的的系统界面，这时可在文本框中输入储油罐的各参数，按下“求解罐容表”按钮，即可输出相应的结果值，下一步可选择进入另一系统，亦可选择退出。

八参考文献

【1】刘卫国 MATLAB程序设计教程，中国水利水电出版社，2006年3月第

三版

【2】韩中庚数学建模方法及其应用，高等教育出版社，2005年6月第一版【3】Frank R.Giordano等，数学建模机械工程出版社，2005年1月第三版【4】朱德通，优化模型与实验，统计大学出版社，2003年4月第一版

【5】Price K，Storn R，Lampinen J．Differential Evolution-A Practical Approach

to Global Optimization [M]．Berlin：Springer-Verlag Press，2005．

【6】栾丽君、谭立静、牛奔。一种基于粒子群优化算法和差分进化算法的新

型混合全局优化算法[J]。信息与控制,2007,36(6):708-714。

【7】Zhang M, Luo W J, Wang X F. Differential evolution with dynamic stochastic

selection for constrained optimization[J]. In- formation Sciences: an International Journal, 2008, 178(15):3043-3074.

附件：

附表一：第一问罐容标定表 h(mm) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 v(L) 1.617955976 3.412163248 6.052647302 9.638936008 14.25954153 19.99427942 26.91523492 35.09207753 44.5869397 55.45876512 67.76301883 81.55178387 96.87453429 113.7782318 132.3074223 152.4983417 174.1828556 197.1226769 221.1901585 246.29273 272.3562769 299.3187634 327.1265604 355.732513 385.0941043 415.1729726 445.9337519 477.3440712 509.3737825 541.9943801 575.1794843 h(mm) 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.6 0.61 v(L) 608.9039716 643.1440712 677.877172 713.0818223 748.7372469 784.8236369 821.3220529 858.2140388 895.4819115 933.1083744 971.0768074 1009.370977 1047.975229 1086.874006 1126.052427 1165.495805 1205.189742 1245.120033 1285.273055 1325.635085 1366.192788 1406.932928 1447.842651 1488.908913 1530.119249 1571.461 1612.921702 1654.489178 1696.151255 1737.895563 1779.710313 h(mm) 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 v(L) 1821.583427 1863.502923 1905.456916 1947.433617 1989.421043 2031.4076 2073.381208 2115.330079 2157.242425 2199.106263 2240.909708 2282.640777 2324.287296 2365.837187 2407.278273 2448.598282 2489.784558 2530.82473 2571.706044 2612.415649 2652.940498 2693.267452 2733.382981 2773.273462 2812.925172 2852.323811 2891.45498 2930.30409 2968.855872 3007.094963 3045.005515 h(mm) 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18 1.19 1.2 v(L) 3082.571295 3119.775393 3156.600801 3193.029355 3229.04279 3264.621781 3299.746519 3334.396036 3368.548783 3402.181859 3435.271397 3467.791887 3499.716659 3531.017111 3561.662612 3591.620211 3620.854155 3649.325407 3676.99087 3703.802132 3729.704502 3754.63401 3778.514225 3801.247743 3822.697062 3842.602346 3860.847738 3877.475269 附表二：无变位进油时理论容积值与实际值的偏差

x轴 h 159.02 176.14 192.59 208.50 223.93 238.97 253.66 268.04 y轴 v V 理论值影响 312.00 322.8809 362.00 374.6313 412.00 426.3632 462.00 478.13 512.00 529.8501 562.00 581.604 612.00 633.3503 662.00 685.0793 19

3.370% 3.372% 3.369% 3.374% 3.369% 3.371% 3.371% 3.369%

282.16 296.03 309.69 323.15 336.44 349.57 362.56 375.42 388.16 400.79 413.32 425.76 438.12 450.40 462.62 474.78 486.89 498.95 510.97 522.95 534.90 546.82 558.72 570.61 582.48 594.35 606.22 618.09 629.96 641.85 653.75 665.67 677.63 678.54 690.53 690.82 702.85 714.91 727.03 739.19 751.42 763.70 764.16 712.00 762.00 812.00 862.00 912.00 962.00 1012.00 1062.00 1112.00 1162.00 1212.00 1262.00 1312.00 1362.00 1412.00 1462.00 1512.00 1562.00 1612.00 1662.00 1712.00 1762.00 1812.00 1862.00 1912.00 1962.00 2012.00 2062.00 2112.00 2162.00 2212.00 2262.00 2312.00 2315.83 2365.83 2367.06 2417.06 2467.06 2517.06 2567.06 2617.06 2666.98 2668.83 736.845 788.5759 840.3271 892.0544 943.8006 995.5405 1047.295 1099.053 1150.806 1202.552 1254.292 1306.03 1357.773 1409.49 1461.234 1512.977 1564.736 1616.484 1668.241 1719.982 1771.728 1823.457 1875.191 1926.953 1978.677 2030.432 2082.196 2133.95 2185.672 2237.431 2289.161 2340.885 2392.669 2396.603 2448.369 2449.619 2501.394 2553.113 2604.882 2656.589 2708.336 2760.008 2761.938 20

3.372% 3.370% 3.371% 3.369% 3.369% 3.369% 3.370% 3.371% 3.372% 3.372% 3.372% 3.371% 3.371% 3.369% 3.369% 3.369% 3.370% 3.371% 3.371% 3.371% 3.371% 3.370% 3.370% 3.371% 3.370% 3.370% 3.371% 3.372% 3.371% 3.371% 3.371% 3.370% 3.371% 3.370% 3.371% 3.370% 3.371% 3.371% 3.371% 3.370% 3.370% 3.371% 3.371%

776.53 788.99 801.54 814.19 826.95 839.83 852.84 866.00 879.32 892.82 892.84 906.53 920.45 934.61 949.05 963.80 978.91 994.43 1010.43 1026.99 1044.25 1062.37 1081.59 1102.33 1125.32 1152.36 1193.49

2718.83 2768.83 2818.83 2868.83 2918.83 2968.83 3018.83 3068.83 3118.83 3168.83 3168.91 3218.91 3268.91 3318.91 3368.91 3418.91 3468.91 3518.91 3568.91 3618.91 3668.91 3718.91 3768.91 3818.91 3868.91 3918.91 3968.91 2813.663 2865.418 2917.168 2968.917 3020.665 3072.41 3124.143 3175.89 3227.633 3279.382 3279.458 3331.178 3382.936 3434.673 3486.424 3538.165 3589.919 3641.668 3693.417 3745.136 3796.889 3848.648 3900.384 3952.137 4003.861 4055.612 4107.36 3.370% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.370% 3.371% 3.370% 3.371% 3.371% 3.371% 3.371% 3.371% 3.370% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 3.371% 附表三：倾斜变位时标高对容积的影响

x 轴h 411.29 423.45 438.33 450.54 463.90 477.74 489.37 502.56 514.69 526.84 538.88 551.96 564.40 Y轴 v（实验）校正V（理论）影响 962.86 975.999 1.3462% 1012.86 1022.655 0.9578% 1062.86 1080.358 1.6197% 1112.86 1128.176 1.3576% 1162.86 1180.947 1.5316% 1212.86 1236.076 1.8782% 1262.86 1282.737 1.5496% 1312.79 1336.000 1.7373% 1362.79 1385.278 1.6234% 1412.73 1434.898 1.5449% 1462.73 1484.302 1.4534% 1512.73 1538.212 1.6566% 1562.73 1589.690 1.6959% 21

576.56 588.74 599.56 611.62 623.44 635.58 646.28 658.59 670.22 680.63 693.03 704.67 716.45 727.66 739.39 750.90 761.55 773.43 785.39 796.04 808.27 820.80 832.80 844.47 856.29 867.60 880.06 892.92 904.34 917.34 929.90 941.42 954.60 968.09 980.14 992.41 1006.34 1019.07 1034.24 1035.36

1612.73 1662.73 1712.73 1762.73 1812.73 1862.73 1912.73 1962.73 2012.73 2062.73 2112.73 2162.73 2212.73 2262.73 2312.73 2362.73 2412.73 2462.73 2512.73 2562.73 2612.73 2662.73 2712.73 2762.73 2812.73 2862.73 2912.73 2962.73 3012.73 3062.73 3112.73 3162.73 3212.73 3262.73 3312.73 3362.73 3412.73 3462.73 3512.73 3514.74 1640.179 1690.897 1736.057 1786.490 1835.999 1886.910 1931.816 1983.501 2032.331 2076.025 2128.034 2176.800 2226.077 2272.883 2321.750 2369.572 2413.691 2462.741 2511.924 2555.537 2605.386 2656.174 2704.522 2751.243 2798.243 2842.892 2891.689 2941.593 2985.491 3034.954 3082.197 3125.028 3173.408 3222.197 3265.117 3308.141 3356.108 3399.077 3449.131 3452.775 1.6735% 1.6658% 1.3437% 1.3300% 1.2674% 1.2814% 0.9880% 1.0472% 0.9645% 0.6404% 0.7192% 0.6463% 0.5996% 0.4467% 0.3885% 0.2887% 0.0398% 0.0004% -0.0321% -0.2815% -0.2819% -0.2468% -0.3035% -0.4175% -0.5177% -0.6978% -0.7276% -0.7186% -0.9124% -0.9152% -0.9906% -1.2064% -1.2391% -1.2579% -1.4582% -1.6501% -1.6871% -1.8727% -1.8439% -1.7947% 附表四：第二问罐容标定表 h(m) v(L) h(m) v(L) h(m) 22

v(L) h(m) v(L)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 45.88 355.42 1068.68 2228.22 3709.18 5440.16 7379.69 9496.79 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 11765.98 14165.19 16674.52 19275.63 21951.26 24684.93 27460.70 30263.01 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 33076.54 35886.04 38676.23 41431.69 44136.67 46774.98 49329.76 51783.21 2.4 2.5 2.6 2.7 2.8 2.9 3 2.4 54116.25 56307.93 58334.55 60168.13 61773.08 63097.32 64029.43 54116.25 附图一：附表一游标高度与容积的直观参考

附图一：问题一、二罐容表标定曲线图

附图二：进入GUI选择系统页面

附图三：GUI问题一系统主页面

附图四：GUI问题一系统求解页面

附图五：GUI问题二系统主页面

附图六：GUI问题一系统求解主页面

附图七：NDEPSO求解结果截图

附图八：NDEPSO求解结果截图

%程序1

%拟合附件一四张表的油位于容量的关系曲线,仅作参考用 X1=[]; X2=[]; X3=[]; X4=[];

[P,S]=polyfit(X1,Y1,3); [Q,S]=polyfit(X2,Y2,3); [M,S]=polyfit(X3,Y3,3); [N,S]=polyfit(X4,Y4,3);

x1=100:0.01:1200;

y1=P(1).*x1.*x1.*x1+P(2).*x1.*x1+P(3).*x1+P(4); y2=Q(1).*x1.*x1.*x1+Q(2).*x1.*x1+Q(3).*x1+Q(4); y3=M(1).*x1.*x1.*x1+M(2).*x1.*x1+M(3).*x1+M(4); y4=N(1).*x1.*x1.*x1+N(2).*x1.*x1+N(3).*x1+N(4);

plot(x1,y1,'r'); hold on

plot(x1,y2,'b'); hold on

plot(x1,y3,'g'); hold on

plot(x1,y4,'y'); hold on

程序2 无偏位时，通过出油数据，检验进油时h与V的模型可靠性 global h c=[]; d=[]; for i=1:73

a=1.78/2;b=0.6;ai=0*pi/180; h=d(i)/1000; if h<2.05*tan(ai)

[S(i),n]=quadl('fx1',0,0.4+h/tan(ai)); elseif 2.05*tan(ai)<=h&h<2*b-0.4*tan(ai) [S(i),n]=quadl('fx1',0,2.45); elseif 2*b-0.4*tan(ai)<=h<=2*b

[S(i),n]=quadl('fx2',0.4-(2*b-h)/tan(ai),2.45);S(i)=pi*a*b*2.45-S(i); end

S(i)=S(i)*1000*(1-0.0337);

h=d(i+1)/1000; if h<2.05*tan(ai)

[T(i),n]=quadl('fx1',0,0.4+h/tan(ai)); elseif 2.05*tan(ai)<=h&h<2*b-0.4*tan(ai)

[T(i),n]=quadl('fx1',0,2.45); elseif 2*b-0.4*tan(ai)<=h<=2*b

[T(i),n]=quadl('fx2',0.4-(2*b-h)/tan(ai),2.45);T(i)=pi*a*b*2.45-T(i); end

T(i)=T(i)*1000*(1-0.0337); P(i)=(S(i)-T(i))-(c(i+1)-c(i)); end e=0; for i=1:73

e=e+abs(P(i)/(c(i+1)-c(i))); end e=e/73 程序3

% 解决第一问罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。 global h

a=1.78/2;b=0.6;ai=4.1*pi/180;i=1; for h=0:0.01:1.2 if h<2.05*tan(ai)

[S(i),n]=quadl('fx1',0,0.4+h/tan(ai));i=i+1; elseif 2.05*tan(ai)<=h&h<2*b-0.4*tan(ai) [S(i),n]=quadl('fx1',0,2.45);i=i+1; elseif 2*b-0.4*tan(ai)<=h<=2*b

[S(i),n]=quadl('fx2',0.4-(2*b-h)/tan(ai),2.45);S(i)=pi*a*b*2.45-S(i);i=i+1; end

S(i-1)=S(i-1)*1000; End %程序4

% 根据所读的油位高度h、变位系数α、β计算总油量V global h r ai R bt H

r=1.625;ai=2.1*pi/180;R=1.5;bt=4.2*pi/180; i=1;

for H=0:0.01:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1(i),n]=quadl('vzuojifen1',0.625,xm);i=i+1; elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1(i),n]=quadl('vzuojifen2',0.625,xm);S1(i)=pi*(r-1/3)-S1(i);i=i+1; elseif 3-2*tan(ai)<=h&h<=3 S1(i)=pi*(r-1/3);i=i+1; end

S1(i-1)=S1(i-1)*1000; end i=1;

for H=0:0.01:3

h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,2+h/tan(ai));i=i+1; elseif 6*tan(ai)<=h&h<3-2*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,8);i=i+1; elseif 3-2*tan(ai)<=h

[S2(i),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2(i)=pi*R*R*8-S2(i);i=i+1; end

S2(i-1)=S2(i-1)*1000; end i=1;

for H=0:0.01:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3(i)=0;i=i+1;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[S3(i),n]=quadl('vyoujifen1',0.625,xm);i=i+1; elseif 1.5+6*tan(ai)<=h&h<=3

[S3(i),n]=quadl('vyoujifen2',0.625,xm);S3(i)=pi*(r-1/3)-S3(i);i=i+1; end

S3(i-1)=S3(i-1)*1000; end

S4=S1+S2+S3; %程序5

%解决第二问罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。 global h r ai R bt H

r=1.625;ai=2.1*pi/180;R=1.5;bt=4.4*pi/180; i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1(i),n]=quadl('vzuojifen1',0.625,xm);i=i+1; elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1(i),n]=quadl('vzuojifen2',0.625,xm);S1(i)=pi*(r-1/3)-S1(i);i=i+1; elseif 3-2*tan(ai)<=h&h<=3 S1(i)=pi*(r-1/3);i=i+1; end

S1(i-1)=S1(i-1)*1000; end i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,2+h/tan(ai));i=i+1; elseif 6*tan(ai)<=h&h<3-2*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,8);i=i+1; elseif 3-2*tan(ai)<=h

[S2(i),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2(i)=pi*R*R*8-S2(i);i=i+1; end

S2(i-1)=S2(i-1)*1000; end i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3(i)=0;i=i+1;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[S3(i),n]=quadl('vyoujifen1',0.625,xm);i=i+1; elseif 1.5+6*tan(ai)<=h&h<=3

[S3(i),n]=quadl('vyoujifen2',0.625,xm);S3(i)=pi*(r-1/3)-S3(i);i=i+1;

end

S3(i-1)=S3(i-1)*1000; end

S4=S1+S2+S3; 程序6

%根据附件2两次加油之间的数据求α、β的第一种方法，算法不够优化 tic

global h r ai R bt H D=[ ]; E=[ ];

P=zeros(50,50); for i=1:301 r=1.625;R=1.5; for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1(s,t),n]=quadl('vzuojifen1',0.625,xm); elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1(s,t),n]=quadl('vzuojifen2',0.625,xm);S1(s,t)=pi*(r-1/3)-S1(s,t); elseif 3-2*tan(ai)<=h&h<=3 S1(s,t)=pi*(r-1/3); end

S1(s,t)=S1(s,t)*1000; end end

for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i)/1000; h=(H-R)*cos(bt)+R;

if h<6*tan(ai)

[S2(s,t),n]=quadl('vzhongjifen1',0,2+h/tan(ai)); elseif 6*tan(ai)<=h&h<3-2*tan(ai)

[S2(s,t),n]=quadl('vzhongjifen1',0,8); elseif 3-2*tan(ai)<=h

[S2(s,t),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2(s,t)=pi*R*R*8-S2(s,t); end

S2(s,t)=S2(s,t)*1000; end end

for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3(s,t)=0;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[S3(s,t),n]=quadl('vyoujifen1',0.625,xm); elseif 1.5+6*tan(ai)<=h&h<=3

[S3(s,t),n]=quadl('vyoujifen2',0.625,xm);S3(s,t)=pi*(r-1/3)-S3(s,t); end

S3(s,t)=S3(s,t)*1000; end end S4=S1+S2+S3;

for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i+1)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[T1(s,t),n]=quadl('vzuojifen1',0.625,xm); elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[T1(s,t),n]=quadl('vzuojifen2',0.625,xm);T1(s,t)=pi*(r-1/3)-T1(s,t); elseif 3-2*tan(ai)<=h&h<=3 T1(s,t)=pi*(r-1/3); end

T1(s,t)=T1(s,t)*1000; end end

for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i+1)/1000; h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[T2(s,t),n]=quadl('vzhongjifen1',0,2+h/tan(ai)); elseif 6*tan(ai)<=h&h<3-2*tan(ai)

[T2(s,t),n]=quadl('vzhongjifen1',0,8); elseif 3-2*tan(ai)<=h

[T2(s,t),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);T2(s,t)=pi*R*R*8-S2(s,t); end

T2(s,t)=T2(s,t)*1000; end end

for s=1:50 for t=1:50

ai=s*pi/1800;bt=t*pi/1800; H=E(i+1)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) T3(s,t)=0;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[T3(s,t),n]=quadl('vyoujifen1',0.625,xm); elseif 1.5+6*tan(ai)<=h&h<=3

[T3(s,t),n]=quadl('vyoujifen2',0.625,xm);T3(s,t)=pi*(r-1/3)-T3(s,t); end

T3(s,t)=T3(s,t)*1000; end end T4=T1+T2+T3;

P=P+abs(S4-T4-D(i)); i=i+1; end

m(i)=1;o(i)=1;min=P(1,1); for p=1:50 for q=1:50 if min>P(p,q)

m(i)=p;o(i)=q;min=P(p,q); end end end toc %程序7

%根据附件2求α、β的第二种方法，即 % NDEPSO混合优化算法 tic %格式标准化 clear all; clc; format long; tic

%初始化各个因子

Fmin=0.3;Fmax=0.9;CRmin=0.3;CRmax=0.9;%DE初始化参数 c1=1.4962; %学习因子c1 c2=1.4962; %学习因子c2 w=0.7298; %惯性权重w N=20; %粒子群规模 D=2; %搜索空间维数 %eps=10^(-6); %满足条件限制的误差 MaxDT=100; %粒子群繁殖的代数

%初始化粒子的速度和位置，数据结构用矩阵A表示 for i=1:N for j=1:2*D

A(i,j)=rand*5*pi/180;

end end for i=1:N

for j=2*D+1:3*D A(i,j)=A(i,j-2*D); end end

%计算各个粒子的适应度 for i=1:N

A(i,3*D+1)=fitness(A(i,1:D)); end

%对粒子的适应度进行排序 B=sortrows(A,3*D+1);

%排序后适应度低的前面一半粒子直接进入下一代 %下面进入主要循环，循环到最大次数得到最优解和最小值 %DDTime=1; for time=1:MaxDT for i=1:N/2 for j=1:3*D+1

NextGeneration(i,j)=B(i,j); end end

for i=1:N/2 for j=1:3*D+1

Cross(i,j)=B(i+N/2,j); end end %进行DE操作

CR(time)=0.27/(CRmax+(CRmin-CRmax)*exp(-60*(1-time/MaxDT)^3)); F(time)=0.27/(Fmax+(Fmin-Fmax)*exp(-10*(1-time/MaxDT)^4)); rot = (0:1:N/2-1);

ind = randperm(4); % index pointer array

a1 = randperm(N/2); % shuffle locations of vectors rt = rem(rot+ind(1),N/2); % rotate indices by ind(1) positions a2 = a1(rt+1); % rotate vector locations rt = rem(rot+ind(2),N/2); a3 = a2(rt+1);

pm1 = Cross(a1,1:D); % shuffled population 1 pm2 = Cross(a2,1:D); % shuffled population 2 pm3 = Cross(a3,1:D); % shuffled population 3

mui = rand(N/2,D) < CR(time); % all random numb1ers < CR are 1, 0 otherwise mpo = mui < 0.5; % inverse mask to mui

Cross(1:N/2,1:D) = pm3 + F(time)*(pm1 - pm2); % differential variation

Cross(1:N/2,1:D) = Cross(1:N/2,2*D+1:3*D).*mpo + Cross(1:N/2,1:D).*mui; % binomial crossover

%变异交叉结束，进行更新 for i=1:N/2

Cross(i,3*D+1)=fitness(Cross(i,1:D)); if Cross(i,3*D+1)Cross(i,j)=Cross(i,j-2*D); end else

for j=2*D+1:3*D Cross(i,j)=B(i,j); end end end

%下面选择最好的粒子N/2个进入下一代 Pool=zeros(N,3*D+1); for i=1:N/2 for j=1:3*D+1

Pool(i,j)=B(i+N/2,j); end end

for i=1+N/2:N for j=1:3*D+1

Pool(i,j)=Cross(i-N/2,j); end end

PoolX=sortrows(Pool,3*D+1); for i=1+N/2:N for j=1:3*D+1

NextGeneration(i,j)=PoolX(i-N/2,j); end end

Pbest=NextGeneration(1,2*D+1:3*D); for i=2:N/2

if NextGeneration(i,3*D+1)%根据粒子群公式进行迭代 for i=1:N/2 for j=D+1:2*D

A(i,j)=w*NextGeneration(i,j)+c1*rand*(NextGeneration(i,j+D)-NextGeneration(i,j-D))+c2*rand*(

Pbest(j-D)-NextGeneration(i,j-D)); end end for i=1:N/2 for j=1:D

A(i,j)=NextGeneration(i,j)+A(i,j+D); end

A(i,3*D+1)=fitness(A(i,1:D));

if A(i,3*D+1)for j=2*D+1:3*D

A(i,j)=NextGeneration(i,j-2*D); end end end

for i=N/2+1:N

A(i,:)=NextGeneration(i,:); end

B=sortrows(A,3*D+1); Pg(time)=fitness(Pbest); end toc for i=2:100 if Pg(i)>Pg(i-1) Pg(i)=Pg(i-1); end end

%算法结束，得到的结果显示如下：

disp('****************************************************') disp('最后得到的最优位置为：') X=Pbest'*180/pi

disp('得到的函数最小值为：') Minimize= Pg(time)

disp('****************************************************') %绘制进化代数和适应度关系曲线图 xx=linspace(1,MaxDT,MaxDT); yy=Pg(xx); plot(xx,yy,'b-') hold on grid on

title('f1目标函数收敛曲线图')

legend('NDEPSO算法适应度曲线走势'); xlabel('迭代代数') ylabel('目标函数对数值') toc

%程序8 NDEPSO混合优化算法的主调用函数 function result=fitness(xxx) global h r ai R bt H

if xxx(1)>5*pi/180|xxx(2)>5*pi/180|xxx(1)<0|xxx(2)<0 P=inf; result=P; else D=[ ]; E=[ ]; P=0; for i=1:301 r=1.625;R=1.5;

ai=xxx(1);bt=xxx(2); H=E(i)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1,n]=quadl('vzuojifen1',0.625,xm); elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1,n]=quadl('vzuojifen2',0.625,xm);S1=pi*(r-1/3)-S1; elseif 3-2*tan(ai)<=h&h<=3 S1=pi*(r-1/3); end

S1=S1*1000; ai=xxx(1);bt=xxx(2); H=E(i)/1000; h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[S2,n]=quadl('vzhongjifen1',0,2+h/tan(ai)); elseif 6*tan(ai)<=h&h<3-2*tan(ai) [S2,n]=quadl('vzhongjifen1',0,8);

elseif 3-2*tan(ai)<=h

[S2,n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2=pi*R*R*8-S2; end S2=S2*1000;

ai=xxx(1);bt=xxx(2); H=E(i)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3=0;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai) [S3,n]=quadl('vyoujifen1',0.625,xm); elseif 1.5+6*tan(ai)<=h&h<=3

[S3,n]=quadl('vyoujifen2',0.625,xm);S3=pi*(r-1/3)-S3; end

S3=S3*1000; S4=S1+S2+S3;

ai=xxx(1);bt=xxx(2); H=E(i+1)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[T1,n]=quadl('vzuojifen1',0.625,xm); elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[T1,n]=quadl('vzuojifen2',0.625,xm);T1=pi*(r-1/3)-T1; elseif 3-2*tan(ai)<=h&h<=3 T1=pi*(r-1/3); end

T1=T1*1000; ai=xxx(1);bt=xxx(2); H=E(i+1)/1000; h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[T2,n]=quadl('vzhongjifen1',0,2+h/tan(ai));

elseif 6*tan(ai)<=h&h<3-2*tan(ai) [T2,n]=quadl('vzhongjifen1',0,8); elseif 3-2*tan(ai)<=h

[T2,n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);T2=pi*R*R*8-S2; end T2=T2*1000;

ai=xxx(1);bt=xxx(2); H=E(i+1)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) T3=0;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai) [T3,n]=quadl('vyoujifen1',0.625,xm); elseif 1.5+6*tan(ai)<=h&h<=3

[T3,n]=quadl('vyoujifen2',0.625,xm);T3=pi*(r-1/3)-T3; end

T3=T3*1000; T4=T1+T2+T3;

P=P+((S4-T4)-D(i))^2; i=i+1; end result=P; end end

%程序9 第一种算法以及NDEPSO混合优化算法都用到的调用函数（共6个） function f = vuojifen1( x )

%VUOJIFEN Summary of this function goes here % Detailed explanation goes here global h r ai

z2=(x+1.375).*tan(ai)+h-1.5; y1=-(r*r-x.*x-z2.*z2).^0.5; y2=(r*r-x.*x-z2.*z2).^0.5; st1=asin(y1./((r*r-x.*x).^0.5)); st2=asin(y2./((r*r-x.*x).^0.5));

f=z2.*(y2-y1)+0.5.*(r*r-x.*x).*(st2-st1)+0.25.*(r*r-x.*x).*(sin(2.*st2)-sin(2.*st1)); end

function f = vzuojifen2( x )

%VZUOJIFEN2 Summary of this function goes here % Detailed explanation goes here global h r ai

z2=(x+1.375).*tan(ai)+h-1.5; y1=-(r*r-x.*x-z2.*z2).^0.5; y2=(r*r-x.*x-z2.*z2).^0.5; st1=asin(y1./((r*r-x.*x).^0.5)); st2=asin(y2./((r*r-x.*x).^0.5));

f=-z2.*(y2-y1)+0.5.*(r*r-x.*x).*(st2-st1)+0.25.*(r*r-x.*x).*(sin(2.*st2)-sin(2.*st1)); end

function f = vzhongjifen1( x )

%VZHONGJIFEN1 Summary of this function goes here % Detailed explanation goes here global h ai R H=h+(2-x).*tan(ai);

f=pi/2.*R.*R+R.*R.*asin((H-R)/R)+(H-R).*sqrt((2.*R.*H-H.*H)); end

function f = vzhongjifen2( x )

%VZHONGJIFEN2 Summary of this function goes here % Detailed explanation goes here global h R ai

H=3-h-(2-x).*tan(ai);

f=pi/2.*R.*R+R.*R.*asin((H-R)/R)+(H-R).*sqrt((2.*R.*H-H.*H)); end

function f = vyoujifen1( x )

%VYOUJIFEN1 Summary of this function goes here % Detailed explanation goes here global h r ai

z2=(-x-5.375).*tan(ai)+h-1.5; y1=-(r*r-x.*x-z2.*z2).^0.5; y2=(r*r-x.*x-z2.*z2).^0.5; st1=asin(y1./((r*r-x.*x).^0.5)); st2=asin(y2./((r*r-x.*x).^0.5));

f=z2.*(y2-y1)+0.5.*(r*r-x.*x).*(st2-st1)+0.25.*(r*r-x.*x).*(sin(2.*st2)-sin(2.*st1)); end

function f = vyouujifen2( x )

%VYOUUJIFEN2 Summary of this function goes here

% Detailed explanation goes here global h r ai

z2=(-x-5.375).*tan(ai)+h-1.5; y1=-(r*r-x.*x-z2.*z2).^0.5; y2=(r*r-x.*x-z2.*z2).^0.5; st1=asin(y1./((r*r-x.*x).^0.5)); st2=asin(y2./((r*r-x.*x).^0.5));

f=-z2.*(y2-y1)+0.5.*(r*r-x.*x).*(st2-st1)+0.25.*(r*r-x.*x).*(sin(2.*st2)-sin(2.*st1)); end %程序10

%根据求的α、β 解决第二问罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为10cm的 %罐容表标定值。 global h r ai R bt H

r=1.625;ai=2.1*pi/180;R=1.5;bt=4.4*pi/180; i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1(i),n]=quadl('vzuojifen1',0.625,xm);i=i+1; elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1(i),n]=quadl('vzuojifen2',0.625,xm);S1(i)=pi*(r-1/3)-S1(i);i=i+1; elseif 3-2*tan(ai)<=h&h<=3 S1(i)=pi*(r-1/3);i=i+1; end

S1(i-1)=S1(i-1)*1000; end i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,2+h/tan(ai));i=i+1; elseif 6*tan(ai)<=h&h<3-2*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,8);i=i+1; elseif 3-2*tan(ai)<=h

[S2(i),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2(i)=pi*R*R*8-S2(i);i=i+1; end

S2(i-1)=S2(i-1)*1000;

end i=1;

for H=0:0.1:3

h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3(i)=0;i=i+1;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[S3(i),n]=quadl('vyoujifen1',0.625,xm);i=i+1; elseif 1.5+6*tan(ai)<=h&h<=3

[S3(i),n]=quadl('vyoujifen2',0.625,xm);S3(i)=pi*(r-1/3)-S3(i);i=i+1; end

S3(i-1)=S3(i-1)*1000; end

S4=S1+S2+S3; %程序11

%画两问的罐容表标定图 x轴为h ,y为V x1=[ ] y1=[ ]; X2=[]; Y2=[];

subplot(2,1,1), plot(x1,y1,'b-.') grid on;

subplot(2,1,2), plot(x2,y2,'r-.') grid on; box on; %程序12

%根据附表2第二次加油后的数据对求的ai=2.1 bt=4.2 进行验算 global h r ai R bt H D=[]; E=[];

r=1.625;ai=2.1*pi/180;R=1.5;bt=4.2*pi/180;e=0; for i=1:301 H=D(i)/1000; h=(H-R)*cos(bt)+R;

x0=0.625;

z0=h+2*tan(ai)-1.5; A=1+tan(ai)^2;

B=2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0-2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<1.5-3*tan(ai)

[S1(i),n]=quadl('vzuojifen1',0.625,xm); elseif 1.5-3*tan(ai)<=h&h<3-2*tan(ai)

[S1(i),n]=quadl('vzuojifen2',0.625,xm);S1(i)=pi*(r-1/3)-S1(i); elseif 3-2*tan(ai)<=h&h<=3 S1(i)=pi*(r-1/3); end

S1(i)=S1(i)*1000; end

for i=1:301 H=D(i)/1000; h=(H-R)*cos(bt)+R; if h<6*tan(ai)

[S2(i),n]=quadl('vzhongjifen1',0,2+h/tan(ai)); elseif 6*tan(ai)<=h&h<3-2*tan(ai) [S2(i),n]=quadl('vzhongjifen1',0,8); elseif 3-2*tan(ai)<=h

[S2(i),n]=quadl('vzhongjifen2',2-(3-h)/tan(ai),8);S2(i)=pi*R*R*8-S2(i); end

S2(i)=S2(i)*1000; end i=1; for i=1:301 H=D(i)/1000; h=(H-R)*cos(bt)+R; x0=0.625;

z0=h-6*tan(ai)-1.5; A=1+tan(ai)^2;

B=-2*z0*tan(ai)-2*x0*tan(ai)^2;

C=x0^2*tan(ai)^2+z0*z0+2*z0*x0*tan(ai)-r*r; xm=(-B+(B*B-4*A*C)^0.5)/2/A; if h<6*tan(ai) S3(i)=0;

elseif 6*tan(ai)<=h&h<1.5+6*tan(ai)

[S3(i),n]=quadl('vyoujifen1',0.625,xm); elseif 1.5+6*tan(ai)<=h&h<=3

[S3(i),n]=quadl('vyoujifen2',0.625,xm);S3(i)=pi*(r-1/3)-S3(i);

end

S3(i)=S3(i)*1000; end

S4=S1+S2+S3; for i=1:300

F(i)=E(i)-(S4(i)-S4(i+1)); e=e+F(i)/E(i); end

e=real(e/300)

附程序中的数据 x1=[ 0

0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000

1.10000000000000 1.20000000000000 1.30000000000000

2.10000000000000 2.20000000000000

0.500000000000000 0.600000000000000 0.700000000000000 0.8000000000000000.900000000000000 1

1.40000000000000 1.50000000000000 1.60000000000000 1.700000000000001.80000000000000 1.90000000000000 2

2.30000000000000 2.40000000000000 2.50000000000000 2.600000000000002.70000000000000 2.80000000000000 2.90000000000000 3]

3709.18000000000 5440.16000000000 7379.69000000000 9496.7900000000011765.9800000000 14165.1900000000 16674.5200000000 19275.630000000021951.2600000000 24684.9300000000 27460.7000000000 30263.010000000033076.5400000000 35886.0400000000 38676.2300000000 41431.690000000044136.6700000000 46774.9800000000 49329.7600000000 51783.210000000054116.2500000000 56307.9300000000 58334.5500000000 60168.130000000061773.0800000000 63097.3200000000 64029.4300000000]

3.41216324800000 6.05264730200000 9.63893600800000

14.2595415300000 19.9942794200000 26.9152349200000 35.092077530000044.5869397000000 55.4587651200000 67.7630188300000 81.551783870000096.8745342900000 113.778231800000 132.307422300000 152.498341700000174.182855600000 197.122676900000 221.190158500000 246.292730000000272.356276900000 299.318763400000 327.126560400000 355.732513000000385.094104300000 415.172972600000 445.933751900000 477.344071200000509.373782500000 541.994380100000 575.179484300000 608.903971600000643.144071200000 677.877172000000 713.081822300000 748.737246900000

y1=[ 45.8800000000000 355.420000000000 1068.68000000000 2228.22000000000

x2=[ 45.8800000000000 355.420000000000 1068.68000000000 2228.22000000000

y2=[1.61795597600000

784.823636900000 821.322052900000 858.214038800000 895.481911500000933.108374400000 971.076807400000 1009.37097700000 1047.975229000001086.87400600000 1126.05242700000 1165.49580500000 1205.189742000001245.12003300000 1285.27305500000 1325.63508500000 1366.192788000001406.93292800000 1447.84265100000 1488.90891300000 1530.119249000001571.46100000000 1612.92170200000 1654.48917800000 1696.151255000001737.89556300000 1779.71031300000 1821.58342700000 1863.502923000001905.45691600000 1947.43361700000 1989.42104300000 2031.407600000002073.38120800000 2115.33007900000 2157.24242500000 2199.106263000002240.90970800000 2282.64077700000 2324.28729600000 2365.837187000002407.27827300000 2448.59828200000 2489.78455800000 2530.824730000002571.70604400000 2612.41564900000 2652.94049800000 2693.267452000002733.38298100000 2773.27346200000 2812.92517200000 2852.323811000002891.45498000000 2930.30409000000 2968.85587200000 3007.094963000003045.00551500000 3082.57129500000 3119.77539300000 3156.600801000003193.02935500000 3229.04279000000 3264.62178100000 3299.746519000003334.39603600000 3368.54878300000 3402.18185900000 3435.271397000003467.79188700000 3499.71665900000 3531.01711100000 3561.662612000003591.62021100000 3620.85415500000 3649.32540700000 3676.990870000003703.80213200000 3729.70450200000 3754.63401000000 3778.514225000003801.24774300000 3822.69706200000 3842.60234600000 3860.847738000003877.47526900000]

152.72

202.72

252.72

302.72

352.72

402.72

452.72

502.72

552.72

602.72

652.72

702.72

752.72

802.72

852.72

902.72

952.72

c=[ 52.72 102.72

1002.72 1052.72 1102.72 1152.72 1202.72 1252.72 1302.72 1352.72 1402.721452.72 1502.72 1552.72 1602.72 1652.72 1702.72 1752.72 1802.72 1852.721902.72 1952.72 2002.72 2052.72 2102.72 2152.72 2202.72 2252.72 2302.722352.72 2402.72 2452.72 2502.72 2552.72 2602.72 2652.72 2702.72 2752.722802.72 2852.72 2902.72 2952.72 3002.72 3052.72 3102.72 3152.72 3202.723252.72 3302.72 3352.72 3402.72 3452.72 3502.72 3552.72 3602.72 3652.723702.72] 978.08 852.15 738.98 630.96 523.95 414.36 297.18 160.48

962.99 839.14 726.81 619.08 511.97 401.84 283.33 142.62]

2483.61000000000 2474.79000000000 2472.37000000000948.26 826.27 714.7 607.21 499.96 389.22 269.24

933.84 813.52 702.64 595.35 487.9 376.49 254.88

919.69 800.87 690.61 583.48 475.8 363.64 240.21

905.78 788.33 678.63 571.61 463.65 350.67 225.21

892.1 775.88 666.68 559.72 451.43 337.55 209.81

878.61 763.51 654.75 547.82 439.15 324.27 193.94

865.3751.21642.84535.9426.8310.82177.54

d=[ 1150.72 1123.99 1101.15 1080.51 1061.36 1043.29 1026.08 1009.54 993.57

D=[ 2486.21000000000

2469.66000000000 2462.44000000000 2459.47000000000 2449.290000000002439.11000000000 2429.89000000000 2426.39000000000 2417.800000000002410.61000000000 2398.88000000000 2390.39000000000 2380.39000000000

2373.85000000000 2367.52000000000 2357.27000000000 2354.440000000002351.11000000000 2347.37000000000 2341.46000000000 2331.150000000002321.11000000000 2318.51000000000 2312.52000000000 2305.250000000002299.08000000000 2290.51000000000 2282.23000000000 2277.31000000000 22712268.84000000000 2257 2253.33000000000 2250.27000000000 2244.540000000002240.56000000000 2233.66000000000 2228.27000000000 2216.750000000002205.55000000000 2203.02000000000 2193.64000000000 2188.950000000002182.72000000000 2175.25000000000 2163.82000000000 2157.640000000002145.81000000000 2140.80000000000 2131.79000000000 2123.120000000002115.73000000000 2091.02000000000 2065.38000000000 2036.28000000000 2003.74000000000 1979.33000000000 1957.15000000000 1934.05000000000 1899.16000000000 1873.50000000000 1836.85000000000 1815.31000000000 1798.73000000000 1766.08000000000 1747.47000000000 1715.34000000000 1681.91000000000 1654.72000000000 1638.80000000000 1607.80000000000 1569.23000000000 1542.45000000000 1522.03000000000 1487.66000000000 1462.95000000000 1436.91000000000 1405.53000000000 1383.48000000000 1359.73000000000 1326.14000000000 1304.81000000000 1274.49000000000 1240.42000000000 1211.43000000000 1187.07000000000 2106.75000000000 2079.03000000000 2054.56000000000 2029.67000000000 1995.29000000000 1972.51000000000 1951.30000000000 1925.88000000000 1889.86000000000 1862.44000000000 1828.91000000000 1811.69000000000 1790.75000000000 1757.69000000000 1740.16000000000 1706.16000000000 1678.85000000000 1651.81000000000 1632.40000000000 1600.62000000000 1564.82000000000 1533.50000000000 1513.36000000000 1478.91000000000 1451.79000000000 1430.30000000000 1398.82000000000 1376.75000000000 1355.82000000000 1316.48000000000 1297.05000000000 1266.05000000000 1228.97000000000 1208.23000000000 1178.45000000000 2098.08000000000 2075.32000000000 2045.87000000000 2017.86000000000 1989.53000000000 1969.31000000000 1943.47000000000 1921.23000000000 1884.42000000000 1851.64000000000 1826.68000000000 1807.90000000000 1784.04000000000 1755.36000000000 1731.62000000000 1694.47000000000 1670.74000000000 1647.15000000000 1625.13000000000 1589.18000000000 1556.06000000000 1530.82000000000 1502.91000000000 1476.84000000000 1449.78000000000 1420.60000000000 1396.46000000000 1373.22000000000 1346.43000000000 1312.60000000000 1288.22000000000 1257.57000000000 1224.88000000000 1200.16000000000 1168.75000000000 48

2094.300000000002072.990000000002041.970000000002014.290000000001985.620000000001961.410000000001938.960000000001910.980000000001876.580000000001841.460000000001820.430000000001801.670000000001775.080000000001752.670000000001725.540000000001687.160000000001660.960000000001643.610000000001618.550000000001580.810000000001551.170000000001526.270000000001497.470000000001468.820000000001443.160000000001415.380000000001392.700000000001367.810000000001342 1330.830000000001307.720000000001280.750000000001248.780000000001215.790000000001193.660000000001163.25000000000

1154.63000000000 1148.47000000000 1138.05000000000 1127.720000000001123.15000000000 1115.02000000000 1107.20000000000 1099.790000000001089.09000000000 1084.44000000000 1079.26000000000 1076.070000000001064.67000000000 1056.21000000000 1049.42000000000 1041.030000000001033.58000000000 1025.11000000000 1017.67000000000 1008.460000000001001.23000000000 989.300000000000 985.110000000000 982.050000000000978.950000000000 976.320000000000 970.270000000000 963.790000000000958.130000000000 948.500000000000 940.220000000000 930.500000000000919.170000000000 907.440000000000 903.520000000000 900.130000000000891.170000000000 888.230000000000 880.980000000000 873.670000000000863.060000000000 856.210000000000 850.280000000000 841.560000000000832.150000000000 824.950000000000 819.470000000000 815.970000000000808.110000000000 803.490000000000 801.050000000000 791.500000000000787.070000000000 780.650000000000 771.770000000000 766.180000000000756.810000000000 750.870000000000 742.030000000000 732.990000000000726.570000000000 724.370000000000 719.060000000000 712.820000000000708.120000000000 704.150000000000 693.930000000000 687.630000000000676.750000000000 670.840000000000 661.150000000000 655.180000000000645.100000000000 635.550000000000 629.770000000000 625.610000000000615.710000000000 604.210000000000 598.940000000000 590.230000000000583.840000000000 573.500000000000 563.820000000000 560.140000000000549.520000000000 537.630000000000 530.480000000000 519.640000000000511.760000000000 508.210000000000 504.210000000000 498.140000000000488.660000000000 478.400000000000 468.500000000000 463.310000000000455.970000000000 453.070000000000 449.960000000000 446.590000000000437.810000000000 430.860000000000 426.960000000000 420.010000000000416.530000000000 413.980000000000]

189.530000000000 52.9600000000000 59.9600000000000 154.21000000000065.2700000000000 222.810000000000 223.490000000000 205.20000000000076.9600000000000 192.390000000000 162.330000000000 263.860000000000193.370000000000 228.180000000000 151.830000000000 146.070000000000236.810000000000 65.6000000000000 77.5200000000000 87.9500000000000138.590000000000 241.940000000000 236.550000000000 62.4700000000000143.050000000000 173.970000000000 147.010000000000 206.670000000000199.370000000000 118.740000000000 153.100000000000 52.7900000000000288.440000000000 89.5700000000000 77.6700000000000 140.75000000000097.4700000000000 171.120000000000 132.290000000000 286.230000000000279.240000000000 63.3100000000000 236.080000000000 118.130000000000156.010000000000 189.740000000000 290.800000000000 156.390000000000302.070000000000 128.830000000000 230.700000000000 223.700000000000191.250000000000 232.520000000000 224.940000000000 99.020000000000085.9100000000000 312.610000000000 96.9600000000000 60.8200000000000200.150000000000 285.120000000000 228.640000000000 103.620000000000150.010000000000 175.950000000000 313.960000000000 95.1700000000000

E=[ 55

282.160000000000 226.080000000000 154.420000000000 102.990000000000170.380000000000 182.700000000000 86.6000000000000 212.040000000000115.700000000000 159.690000000000 211.380000000000 121.440000000000133.700000000000 220.640000000000 128.180000000000 280.390000000000321.380000000000 255.440000000000 148.380000000000 214.32000000000085.1200000000000 304.040000000000 297.740000000000 280.740000000000127.280000000000 221.330000000000 59.4500000000000 172.550000000000143.400000000000 100.470000000000 103.570000000000 174.33000000000081.7300000000000 219.970000000000 187.550000000000 248.020000000000252.460000000000 145.220000000000 257.720000000000 83.0200000000000 83.8800000000000 179.470000000000 203.480000000000 122.320000000000 252.170000000000 244.880000000000 246.250000000000 312.540000000000 184.050000000000 188.110000000000 187.060000000000 106.680000000000 130.530000000000 79.2600000000000 171.540000000000 227.090000000000 119.200000000000 177.880000000000 230.850000000000 121.640000000000 280.730000000000 298.420000000000 191.840000000000 185.720000000000 78.4400000000000 142.360000000000 281.260000000000 217.700000000000 256.150000000000 225.360000000000 183.130000000000 233.450000000000 204.360000000000 327.080000000000 229.380000000000 129.820000000000 204.590000000000 321.650000000000 248.050000000000 73.9400000000000 292.620000000000 58.6400000000000 55.6600000000000 273.680000000000 63.6300000000000 98.5500000000000 261.040000000000 266.270000000000 214.490000000000 230.760000000000 313.020000000000 87.3500000000000 230.610000000000 165.420000000000 216.120000000000 123.170000000000 219.880000000000 219.720000000000 305.060000000000 65.3400000000000 241.740000000000 291.060000000000 72.7700000000000 164.020000000000 169.900000000000 107.360000000000 64.4100000000000 74.8200000000000239 170.270000000000 284.730000000000204.550000000000 149.650000000000273.310000000000 175.17000000000097.2400000000000 137.820000000000184.450000000000 302.620000000000235.050000000000 326.240000000000137.650000000000 246.300000000000128.710000000000 118.680000000000154.800000000000 273.280000000000226.020000000000 164.190000000000187.100000000000 173.450000000000145.310000000000 276.050000000000105.470000000000 257.940000000000151.440000000000 225.700000000000122.990000000000 309.810000000000106.760000000000 135.630000000000243.090000000000 206.410000000000233.210000000000 239.810000000000111.180000000000 246.200000000000217.580000000000 177.850000000000262.130000000000 148.120000000000278.400000000000 275.320000000000206.560000000000 196.220000000000135.540000000000 84.5100000000000176.670000000000 218.930000000000191.820000000000 235.360000000000107.030000000000 78.4400000000000155.920000000000 162.750000000000208.690000000000 240.84000000000096.0600000000000 84.0300000000000177.340000000000 178.200000000000143.510000000000 208.740000000000129.180000000000 84.110000000000059.9300000000000 223.060000000000

101.510000000000 148.070000000000 205.310000000000 128.350000000000214.020000000000 135.840000000000 200.570000000000 201.280000000000144.710000000000 50.1900000000000 116.750000000000 140.180000000000103.130000000000 86.9700000000000 225.390000000000 137.100000000000235.680000000000 126.010000000000 208.710000000000 127.780000000000213.280000000000 200.680000000000 119.580000000000 87.6300000000000205.050000000000 234.920000000000 107.450000000000 175.480000000000128.200000000000 206 191.450000000000 72.2000000000000 206.040000000000231.250000000000 135.860000000000 205.870000000000 145.64000000000066.6500000000000 74.3100000000000 185.460000000000 177.260000000000 50.1300000000000 54.9100000000000 118.460000000000 65.8100000000000 43.1300000000000]

112.800000000000 91.5200000000000 57.1700000000000 115.300000000000 51

171.720000000000128.340000000000149.83000000000057.0900000000000

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文

全部栏目

清华大学论文模板