全国大学生数学建模竞赛优秀论文

PAGEPAGE16优秀论文选编优秀论文选编按：数学建模竞赛最终的成果体现在于参赛论文，以下我们挑选我院部分获全国一等奖的优秀论文论文中难免有一些小的错误与失误。煤矿瓦斯和煤尘的监测与控制模型摘要我国煤矿每年因事故而死亡人数居世界首位！煤矿安全生产形势仍相当严峻，其大部分煤矿事故都是由瓦斯或煤尘爆炸引起的。因此，做好井下瓦斯和煤尘的监测与控制是实现煤矿安全生产的针对问题得出相应模型，并得到相应合理的结果。针对问题一，根据《煤矿安全规程》第一百三十三条的分类标准，及绝对瓦斯涌出量和相对瓦斯涌出量的计算公式。通过所给的数据，求出煤矿各监测点每天的绝对瓦斯涌出量和相对瓦斯涌出10m3/t30天的平均绝对瓦斯涌出量为9.8m3/min23.2m3/t10m3/t针对问题二，煤矿发生爆炸的可能性为相对的，而不是绝对的。假设只考虑瓦斯爆炸与煤尘爆炸，综合考虑瓦斯爆炸的可能性与煤尘爆炸的可能性，引用煤尘与在瓦斯浓度影响下煤尘的爆炸下限的偏离程度来恒量，由煤尘引起爆炸的可能性；引用瓦斯浓度与瓦斯下限的偏离程度来恒量瓦斯爆炸的可能性。综合两种发生爆炸的可能性，即为该煤矿发生爆炸的可能性。经MATLAB软件编程得出煤矿发生爆炸的不安全性（煤矿发生爆炸事故的可能性）为10.83%。并列表给出了不同瓦斯浓度与煤尘浓度对应的不安全性程度。针对问题三，根据各井巷风量的分流情况，确定最佳总通风量为进风巷、进风巷及局部通风机所在巷（包括局部通风机的风量）各井巷中风速的要求，及考虑瓦斯和煤尘等因素的影响，列出相应约束条件。经LINGO得出最佳总通风量为1415.062m3/min的风量为476.1359m3/min的风量为548.5541m3/min，局部通风机的额定风量331.8158m3/min。同时，本文还作了误差分析，对模型进行了评价及推广，并在做出相应简化假设情况下，对模型作了进一步的改进。关键字：不安全程度函数监测 瓦斯涌出量 风量（注：此文获2006年全国大学生数学建模竞赛全国一等奖）一、问题的提出基本情况生产的关键环节（见附件。瓦斯是一种无毒、无色、无味的可燃气体，其主要成分是甲烷，在矿井中它通常从煤岩裂缝中涌出。瓦斯爆炸需要三个条件：空气中瓦斯达到一定的浓度；足够的氧气；一定温度的引火源。煤尘是在煤炭开采过程中产生的可燃性粉尘。煤尘爆炸必须具备三个条件:煤尘本身具有爆炸炸浓度是302000g/m3国家《煤矿安全规程》给出了煤矿预防瓦斯爆炸的措施和操作规程，以及相应的专业标准(见附件2)。规程要求煤矿必须安装完善的通风系统和瓦斯自动监控系统，所有的采煤工作面、掘进面和回风巷都要安装甲烷传感器，每个传感器都与地面控制中心相连，当井下瓦斯浓度超标时，控制中心将自动切断电源，停止采煤作业，人员撤离采煤现场。具体内容见附件2的第二章和第三章。问题提出附图1是有两个采煤工作面和一个掘进工作面的矿井通风系统示意图，请你结合附表2的监测数据，按照煤矿开采的实际情况研究下列问题：根据《煤矿安全规程》第一百三十三条的分类标准2)矿井”还是“高瓦斯矿井根据《煤矿安全规程》第一百六十八条的规定，并参照附表1，判断该煤矿不安全的程度（即发生爆炸事故的可能性）有多大？风筒实现掘进巷的通风（见下面的注。根据附图1所示各井巷风量的分流情况、对各井巷中风速的要求（见《煤矿安全规程》第一百零一条，以及瓦斯和煤尘等因素的影响，确定该煤矿所需要的最佳（总）通风量，以及两个采煤工作面所需要的风量和局部通风机的额定风量（实际中，井巷可能会出现漏风现象。二、问题的分析背景的分析煤矿安全生产是目前社会重点关注的热点问题之一，尤其是在能源紧张，对煤碳的需求量不断根据统计资料，可知大部分煤矿事故的罪魁祸首都是瓦斯或煤尘爆炸。因此，矿井下的瓦斯和煤尘对煤矿的安全生产构成了重大威胁，做好井下瓦斯和煤尘的监测与控制是实现煤矿安全生产的关键环节。基本预备知识斯涌出量划分为：低瓦斯矿井：矿井相对瓦斯涌出量小于或等于10m3/t，且绝对瓦斯涌出量小于等于40m3/min；10m3/t，或绝对瓦斯涌出量大于40m3/min；相对瓦斯量定义：是指平均日产一吨煤所涌出的瓦斯量，单位为m3/t；绝对瓦斯量定义：是指矿井单位时间涌出的瓦斯体积，单位为m3/d或m3/min。和断电范围必须符合表3规定（具体表3见附件。2（2见附件2。问题的分析1的分析2所给监测值，可根据绝对瓦斯涌出量与相对瓦斯涌出量的计算公式，算出各监测点的绝对瓦斯涌出量与相对瓦斯涌出量。如果经考察出的监测点的相对瓦斯量有小于或等于10m3/t且绝对瓦斯量小于等于40m3/min，则鉴定该煤矿属于低瓦斯矿井。而如果经考察出的监测点的相对瓦斯量有大于10m3/t或绝对瓦斯量大于40m3/min，则鉴定该煤矿属于高瓦斯矿井。2的分析1，判断煤矿不安全的程度（即发生爆炸事故的可能性。即煤矿发生爆炸的可能性为相对的，不是绝对的。因此，假设只考虑瓦斯爆炸与煤尘爆炸，需定义不同浓度瓦斯与煤尘发生爆炸事故的可能性，并综合瓦斯爆炸的可能性与煤尘爆炸的可能性，即为该煤矿发生爆炸事故的可能性.3的分析满足各井巷中风速的要求及瓦斯和煤尘等因素的影响约束，确定煤矿所需要的最佳总通风量，以及两个采煤工作面所需要的风量和局部通风机的额定风量，这是一个有多约束条件的优化问题。三、模型的假设1、各监测站点的工作是相互独立的；2、附表中的监测值均为有效值，忽略其测量误差，且每天各班次的监测数据为该班次内的平均监测值；3、煤矿的生产是严格按照国家《煤矿安全规程》进行生产；4、煤矿爆炸只考虑由瓦斯爆炸和煤尘爆炸，不考虑其他如矿井温度，机器摩擦及一些由人为失误造成的爆炸；5、煤尘爆炸下限取其中位数40g/m3，瓦斯爆炸下限取值为5%。四、符号约定v：监测点的风速（单位：m/s；s：巷道横断面面积（单位：m2；Q：监测点的风量（单位：m3/min；A：矿井的绝对瓦斯涌出量（单位：m3/min；c：风流中的平均瓦斯浓度，即体积百分比（m：矿井中的煤尘（单位：g/m3；B：矿井的相对瓦斯涌出量（单位：m3/tR：矿井的日产量（单位：t/d；：煤尘爆炸下限（单位：g/m3；b：瓦斯爆炸下限（；k：在空气中有瓦斯时，煤尘降低系数；：在空气中有瓦斯时，煤尘发生爆炸的下限（单位：g/m3；mpg：煤尘爆炸对矿井的不安全性大小；qg：瓦斯爆炸对矿井的不安全性大小；z：煤矿的不安全性大小。五、模型的建立与求解1绝对瓦斯涌出量与相对瓦斯涌出量的计算公式由问题的分析，鉴定矿井是属于“低瓦斯矿井”还是“高瓦斯矿井与相对瓦斯涌出量值，与分类标准值进行鉴别。由绝对瓦斯涌出量与相对瓦斯涌出量的定义，结合相关的符号约定，可知风量为风速在1分钟传播的距离乘以相应巷道横断面面积，公式为：Qsv60 绝对瓦斯涌出量计算公式为：AQc/100 ……（2）一天24小时，且1小时60分钟，绝对瓦斯量的单位为m3/min，所以相对瓦斯涌出量的计算公式为：

B2460A/R ……（3）矿井的绝对瓦斯涌出量与相对瓦斯涌出量的计算2公式可得出矿井在各监测点各班次的绝对瓦斯量及相对瓦斯量。对各监测点进行编号，有i个监测点，i1,26，依次表示工作面，工作面，掘进工作面，回风巷，回风巷；对矿井日生产班次进行编号，有j个班次，j1,2,3，依次表示为早班，30nn1,230。则对应的有第n天第i个监测点第j个班次的风速、瓦斯的监测值和日产量监测值，分别记为：v(n)，cij

(n)，R(n)；第n天第i个监测点巷道记为s(n)。i30R

130R(n)30n1根据公式1，则可得第n天第i个监测点第j个班次的风速为：Q(n)s(n)v(n)60 ……（4）ij i ij根据公式2，第n天第i个监测点第j个班次的绝对瓦斯涌出量为：A(n)s(n)v

(n)

(n)60 ……（5）ij i ij ij第n天第i个监测点的平均绝对瓦斯涌出量为：A(n)13

A(n) ……（6）i 3 j根据公式3，第n天第i个监测点第j个班次的相对瓦斯涌出量为：B(n)s(n)v

(n)

(n)606024/100R ……（7）ij i ij ij第n天第i个监测点的平均相对瓦斯涌出量为：B(n)13

B(n) ……（8）i 3 j第i个监测点30天的平均绝对瓦斯涌出量为：1A1i 30

30n1

A(n) ……（9）i第i个监测点30天的平均相对瓦斯涌出量为：1B1i 30

30n1

B(n) ……（10）i则代入附表2所给的相应数据，可得各监测点30天的平均绝对瓦斯量和相均绝对瓦斯量。煤矿类型的鉴别根据附图（煤矿的通风系统示意图作面、工作面、掘进工作面的瓦斯涌出量之和与回风巷、回风巷的瓦斯涌出量之和与总回风巷的瓦斯涌出量可认为是等价的。这里取总回风巷的相对瓦斯涌出量和绝对瓦斯涌出量的值，依据矿井的分类标准来鉴别该矿井。经MATLAB软件编程得出，30天回风巷的相对瓦斯涌出量B(n)如表1所示，均大于10m3/t，6且30天的平均相对瓦斯涌出量B6

为23.2m3/t10m3/tA6

为9.8m3/min，小于40m3/min。所以，根据矿井的分类标准，该矿井属于高瓦斯矿井。130天的相对瓦斯涌出量：相对瓦斯涌出量，单位m3/t）N12345678910A23.823.3724.0822.8322.7324.2623.0822.5522.5623.14N11121314151617181920A22.7421.4424.2423.3622.2423.0824.7723.2822.6223.12N21222324252627282930A23.2523.1524.6522.8222.9722.8023.3622.4623.1823.562煤矿不安全程度（即发生爆炸事故的可能性）的定义在此假设煤矿爆炸只考虑由瓦斯浓度引起的爆炸和煤尘浓度引起的爆炸，不考虑其他如矿井温度，机器摩擦及一些由人为失误造成的爆炸。引用煤尘与在瓦斯浓度影响下煤尘的爆炸下限的偏离程度来恒量，由煤尘引起爆炸的可能性；引用瓦斯浓度与瓦斯爆炸下限的偏离程度来恒量瓦斯爆炸的可能性。若偏离值越大，煤矿的安全性越好；若偏离值越小，煤矿的安全性越差。在此采用了一个最大型心理函数计算其不安全的可能性。通过给瓦斯不安全程度函数与煤尘不安全程度函数赋予不同权系数，平衡两者的不安全程度，综合两种发生爆炸的可能性即为该煤矿发生爆炸的可能性。综合时，进行了对瓦斯的不安性与煤尘的不安全性赋权处理。煤矿安全性的计算第i监测点第j班次的瓦斯浓度、煤尘、在有瓦斯时煤矿降低系数及相应的煤尘发生爆炸的下限分别为c

，m，k，ij ij

。mij50g/m350g/m3

，取其中位值40g/m3；且瓦斯爆炸下限b，取值为5%。采用求解其偏离值的大小，即其不安全性的大小。则第i监测点第j班次煤尘对矿井的不安全性为：(1(

mij

mij))2pg1eij

mij ……（11）其中，

mij

kij

；用MATLAB软件编程中的线性最小二乘法[3]进行拟合可以将不同瓦斯浓度对应的煤尘降低系数kij

算出。30390第i监测点平均每班次煤尘爆炸对矿井的不安全性大小为：pgi

190pg90 j1

……（12）煤尘爆炸对矿井的不安全性大小先对6个监测点赋权处理综合成一个点，为：pg6

pgw

……（13）iii1则第i监测点第j班次煤尘对矿井的不安全性为：5Bqg1e(1(ij

ij))25

……（14）第i监测点平均每班次瓦斯爆炸对矿井的不安全性大小为：qgi

190qg90 j1

……（15）瓦斯爆炸对矿井的不安全性大小同样先对6个监测点赋权处理综合成一个点，为：qg6

qgw

……（16）iii1煤尘爆炸与瓦斯爆炸可以看为两个相互独立的事件，只要煤尘爆炸或瓦斯爆炸这两个事件任意有一个发生，则整个煤矿是不安全的。根据概率统计知识中任意事件概率的加法公式，得煤矿的不安全性的大小为：zpgqgpgqg ……（17）经MATLAB软件编程得出取其对六个面影响最大的作为整个矿井的不安全程度，即煤矿发生爆炸的不安全性（煤矿发生爆炸事故的可能性）大小为z0.1772。在此根据以上模型把给出的煤尘浓度与瓦斯浓度对应的矿井不安全可能性的大小如列表2如下。根据表2所得结果，人们可以从中看出煤尘浓度与瓦斯浓度对应的矿井不安全性的大小。3总通风量的定义及公式根据附图（煤矿的通风系统示意图）中各巷道的分布位置及各处风的流向（即分流情况把总通风量分为三大块，进风巷、进风巷及局部通风机所在巷（包括局部通风机的风量）的风量（分别记为QQQ。1 2 3AB2：矿井不安全性的大小（BAAB00．10．20．30．40．50．60．70．80．91．06．00.010.010.020.030.040.050.060.080.100.120.156．50.020.020.020.030.040.060.070.090.110.140.167．00.020.020.030.040.050.070.080.100.130.150.187．50.030.030.040.050.060.080.100.170.140.170.208．00.030.040.050.060.070.090.110.130.160.190.228．50.040.050.060.070.090.110.130.150.180.210.259．00.060.060.070.090.100.1230.150.180.210.240.289．50.070.080.900.110.120.150.170.200.240.280.32100.080.090.110.130.150.170.200.230.270.310.3615%（新鲜风动，否则可能会出现负压导致乏风逆流，即局部通风机将乏风吸入并送至掘进工作面。记余裕通风量为w，局部通风量为Q10

Q(1w)。3所以求最小总通风量的目标函数为：Q(QQ1 2

Q) ……（18）3其中Q

60vs；

60

s；

60

s。1 1 1 2风速的约束

2 2 3 3 3由《煤矿安全规程》第一百零一条的规定，得各巷道的风速范围约束。记风速为vi1,29，分别表示进风巷、进风巷、采煤工作面、回风巷、采煤工作i面、回风巷、总回风巷、掘进工作面的风速。处的风量应为进风巷的风量加绝对瓦斯涌出量，等于回风巷的风量。QA

QA

QA所以，v5

1 4 2s 7 s5

6，v9

10s9

9。总回风巷的风量为进风巷、进风巷的风量之和与采煤工作面、采煤工作面、局部通风机所在巷的绝对瓦斯涌出量之和。则总回风巷v8

QQ 1 Q

QA3 s8

AA6 9vi

i1,2v8，0.25v3

4（i4,5,6,7,9v8

8。还需考虑各巷道中瓦斯和煤尘等因素的影响，首先通过附表2所给的数据用MATLAB软件编程采用线性最小二乘法把风速与瓦斯及风速与煤尘的函数关系式插值拟合出来，分别记为 f(v)，F(vf4

(v)，f5

(v)，f6

(v)，f7

(v)，f8

(v)，f(v)分别表示风速采煤工作面，回风巷，采煤工作面，回风巷，总回风巷，掘进工作面对9

(v)0.11v20.79v1.91；f(v)0.3v20.36v4 4 4 5 5

1.34；f(v)3.63v215.69v17.82；f(v)6.12v226.56v

29.73；6 6 6 7 7 7f(v)0.26v22.84v8.33；f(v)0.24v20.98v

0.78；8 8 8 9 9 9由上关系式，得对应风速v，则有各巷道的瓦斯浓度，记为cc4

，c，c，c，c6 7 8

。由所给监测数据，取各监测点瓦斯浓度的最大值，作为最优瓦斯浓度的约束。即c0.78%，c4

1.11%，c6

0.33%，c7

0.83%，c8

1.18%，c9

0.71%。相应的煤尘也有一安全浓度，即煤尘应小于其对应在有瓦斯时煤尘的爆炸下限m

k。对应k值通过附表1，同样采用线性最小二乘法把瓦斯浓度与k的函数关系插值拟合出来。经MATLAB软件编程，得k0.07c20.52c0.98煤尘爆炸下限仍取中位数40g/m3，则 (0.07c20.52c0.98)40m由不同巷道的瓦斯浓度，有其对应的煤尘浓度为： 0.36v42+0.28v4+5.23；m4

=3.26v52-12.55v5+19.43 20.46v62-84.67v6+95.22；m6 m

25.55v72-107.47v7+120.23 0.57v82-5.67v8+21.21；m8

1.24v92-3.62v9+9.24。m4

，mm4 5

，mm5

，mm6

，mm7

，mm8

。m9最佳总风量的模型综上所得，求得其最佳总风量的模型如下：minQQQQ1 2 3s.t.

0.25vi

6（i1,2）v8（i3,8）i0.25vi

4（i4,5,6,7,9,9)，mica,9)，mii i

（i4,,9）,9）w15%a4

0.78%,a5

1.11%,a6

，a0.83%，a7

1.18%，a0.71%，s9

4，s2

4，s3

5，s5

4，s7

4，s8

5，s9

40.132.87模型求解经LINGO软件编程求解，求解最佳总通风量为Q1415.062m3/min，采煤工作面的风量为Q476.1359m3/minQ1

548.5541m3/min，局部通风机的额定风量为Q10

331.8158m3/min。六、误差分析误差来源：1、各监测站点在实际监测中，有观测误差，即存在监测数据与实际数据的误差。2、在模型的建立中，有模型误差。即对监测数据的处理大部分对其取平均值，由模型所得的解与实际问题的解之间存在一定的误差。3、煤矿发生爆炸，在此只考虑瓦斯爆炸和煤尘爆炸，由此所得的煤矿发生爆炸的可能性与实际煤矿生产中有一定的误差。4小，只是取其相对值。而现实中是否发生不安全事故是随机的，不确定的。5、舍入误差：在计算时取的是小数点后两位，其数据有一定的误差。七、模型的改进问题2的改进：7.2.1煤矿安全性的分析7．2．2煤矿安全性的计算第i监测点第j班次的瓦斯浓度、煤尘、在有瓦斯时煤矿降低系数及相应的煤尘发生爆炸的下限分别为c，m，k，ij ij ij

。mij50g/m350g/m3

，取其中位值40g/m3；且瓦斯爆炸下限b，取值为5%。采用最小二乘法求解其偏离值的大小，即其不安全性的大小。则第i监测点第j班次煤尘对矿井的不安全性为：ijpg(1(mij

mij))2 ……（11）ij

mij其中，

mij

kij

；用MATLABkij

可由三次样条插值法进行插值算出。30390第i监测点平均每班次煤尘爆炸对矿井的不安全性大小为：pgi

190pg90 j1

……（12）煤尘爆炸对矿井的不安全性大小取6个监测点的平均值，为：pg1pg6 ii1

……（13）第i监测点平均每班次瓦斯爆炸对矿井的不安全性大小为：qgi

190qg90 j1

……（14）瓦斯爆炸对矿井的不安全性大小取6个监测点的平均值，为：qg1qg6 ii1

……（15）为：zpgqgpgqg ……（16）经MATLAB软件编程得出z10.83%，即煤矿发生爆炸的不安全性（煤矿发生爆炸事故的可能性）为10.83%。问题3的改进：为了计算，我们进一步简化公式，即风速的约束作进一步简化，作为问题3的改进。总通风量的定义及公式根据附图（煤矿的通风系统示意图）中各巷道的分布位置及各处风的流向（即分流情况，可把总通风量分为三大块，进风巷、进风巷及局部通风机所在巷（包括局部通风机的风量）的风量（为QQ

，和各巷道漏的风量。1 2 315%15%，则对应的15%15%（新鲜风）送至掘进工作面。所以根据掘进巷道图，局部通风机所在巷的通风量为局部通风机的通风量除以（115%，记局部通风机的通风量为Q。4所以求最小总通风量的目标函数为：Q(QQ1 2

Q)(10.15) ……（17）Q3Q其中Q

60vs；

60

s；

4，

60

s。1 1 1

2 2

0.85 4 3 3总通风量的约束由《煤矿安全规程》第一百零一条的规定，各巷道的风速范围约束为：0.25v40.250.15

1v42 ……（18）v43150Q4

400还需考虑各巷道中瓦斯和煤尘等因素的影响，首先通过附表2所给的数据用MATLAB软件编程，采用线性最小二乘法把风速与瓦斯及风速与煤尘的函数关系式插值拟合出来，分别记为f(v)，F(v)。通过编f(vf1

(v)，f3

(v)分别表示风速对应进风巷、进风巷及局部通风机所在巷瓦斯的关系式，其为：f(v)0.11v20.79v1.91；f(v)3.63v215.69v

17.821 1 1 2 2 2f(v)0.24v20.98v0.783 3 3vccc1 2 3

。由《煤矿安全规程》第一百六十八条的规定，其瓦斯浓度有一安全浓度约束，得c1.5%，c1

1.5%，c3

1.0%相应的煤尘也有一安全浓度，即煤尘应小于其对应在有瓦斯时煤尘的爆炸下限m

k。对应kk软件编程，得k0.07c20.52c0.98煤尘爆炸下限仍取中位数40g/m3，则 (0.07c20.52c0.98)40m由不同巷道的瓦斯浓度，有其对应的煤尘浓度 ， ， 。m1 m2 m3m1

，mm1 2

，mm2

。m3最佳总风量的模型综上所得，得求其最佳总风量的模型如下：

minQ0w0.850.25v4(i1,2)i0.15v4s.t.c

31.5%(i1,2)ic1.0%3mi mi

(i1,2,3)八、模型的评价及推广模型的优点：1、本文建立的模型能与实际紧密联系，结合实际煤矿生产情况对所提出的问题进行求解，其模型的结果，与实际相符。这对煤矿生产管理部门具有较高的指导价值，使模型更贴近实际，通用性、推广性较强。2、模型原理简单明了，容易理解与灵活运用。3、模型的建立根据问题要求，严格按照《煤矿安全规程》的相关规定，得模型的可信度较高。模型的缺点：1、在模型的建立中，对各监测点的监测值大都取平均值处理，这在实际生产过程中，有一定的误差，使模型不能更准确的反应实际生产情况。2、在考虑煤矿不安全的程度时，忽略了导致煤矿发生爆炸的其他因素，如人为的因素，这与实际也有一定的出入。模型的推广：测站点的监测值，均可运用本模型进行鉴别工种类型。业中不同的项目，根据国际标准质量体系，判断其项目的性质。九、参考文献[J],2005,32(1):36J,2003,22(8)[M,2000十、附录第一问的程序(用MATLAB求解):clearclcloadE:\funy;loadE:\fun;v1=funy(:,1);v2=funy(:,4);v3=funy(:,7);v4=funy(:,10);v5=funy(:,13);v6=funy(:,16);C1=funy(:,2);C2=funy(:,5);C3=funy(:,8);C4=funy(:,11);C5=funy(:,14);C6=funy(:,17);s1=4;s2=5;r1=v1.*s1*60;r2=v2.*s1*60;r3=v3.*s1*60;r4=v4.*s1*60;r5=v5.*s1*60;r6=v6.*s2*60;n=length(funy);Qg1=zeros(n,1);Qg2=zeros(n,1);Qg3=zeros(n,1);Qg4=zeros(n,1);Qg5=zeros(n,1);Qg6=zeros(n,1);qg1=zeros(30,1);qg2=zeros(30,1);qg3=zeros(30,1);qg4=zeros(30,1);qg5=zeros(30,1);qg6=zeros(30,1);fori=1:nQg1(i)=r1(i).*C1(i)/100;Qg2(i)=r2(i).*C2(i)/100;Qg3(i)=r3(i).*C3(i)/100;Qg4(i)=r4(i).*C4(i)/100;Qg5(i)=r5(i).*C5(i)/100;Qg6(i)=r6(i).*C6(i)/100;endx=sum(fun)/30;forj=1:30qg2(j)=480*(Qg2(1+3*(j-1))+Qg2(2+3*(j-1))+Qg2(3+3*(j-1)))/x;qg3(j)=480*(Qg3(1+3*(j-1))+Qg3(2+3*(j-1))+Qg3(3+3*(j-1)))/x;qg4(j)=480*(Qg4(1+3*(j-1))+Qg4(2+3*(j-1))+Qg4(3+3*(j-1)))/x;end

qg5(j)=480*(Qg5(1+3*(j-1))+Qg5(2+3*(j-1))+Qg5(3+3*(j-1)))/x;qg6(j)=480*(Qg6(1+3*(j-1))+Qg6(2+3*(j-1))+Qg6(3+3*(j-1)))/x;Q=[Qg1Qg2Qg3];W1=sum(Q);W2=sum(W1);W3=sum(Qg6)/90;%绝对瓦斯涌出量的平均值R=[qg1qg2qg3];R1=sum(R);R2=sum(R1);R3=sum(qg6)/30;(1)瓦斯的不安全性加权得到偏离度(用MATLAB存函数):functiony=yuanw(x)Q=5;k=3;%加权系数y=1-exp(-k*(1-(Q-x)/Q).^2);%瓦斯的不安全性加权得到偏离度(2):煤尘的不安全性加权得到偏离度(用MATLAB存函数):functiony=yuan(g,x)k=2;u=[00.51.01.52.02.53.03.54];v=[3022.51510.56.54.532.51.5];a=polyfit(u,v,1);%用线性最小二乘法对瓦斯浓度与煤尘暴咋下限的拟合系数Q=polyval(a,g);y=1-exp(-3*(x/Q).^2.*(1-(Q-x)/Q).^2);%煤尘的不安全性加权得到偏离度(3)总考虑煤矿的不安全性(用MATLAB求解):clearclcloadE:\funy;g=zeros(90,6);c=zeros(90,6);fori=1:90forj=1:6g(i,j)=yuanw(funy(i,3*(j-1)+2));c(i,j)=yuan(funy(i,3*(j-1)+2),funy(i,3*(j-1)+3));%在加权得到将瓦斯的偏离度与煤尘的偏离度的结合endendfori=1:90G(i,1)=max(g(i,:));C(i,1)=max(c(i,:));endGC=mean(G+C-G.*C);%将瓦斯与煤尘的偏离度对煤矿不安全性x=[0:0.1:1];y=[6:0.5:10];%是煤尘与瓦斯对应的不安全度Z=zeros(length(x),length(y));fori=1:length(x)forj=1:length(y)Z(i,j)=1-(1-yuanw(x(i)))*(1-yuan(x(i),y(j)));endend第三问的程序:MATLABclearclcx1=[00.500.751.01.502.03.04.0];k=[10.750.600.500.350.250.10.05];A=polyfit(x1,k,2);%用最小二乘法得出瓦斯浓度对煤尘爆炸下限影响函数的系数poly2str(A,'x');%用最小二乘法得出瓦斯浓度对煤尘爆炸下限影响函数loadE:\funy;v4=funy(:,1);v5=funy(:,10);v6=funy(:,4);v7=funy(:,13);v8=funy(:,16);v9=funy(:,7);c4=funy(:,2);c5=funy(:,11);c6=funy(:,5);c7=funy(:,14);c8=funy(:,17);c9=funy(:,8);M4=funy(:,3);M5=funy(:,12);M6=funy(:,6);M7=funy(:,15);M8=funy(:,18);M9=funy(:,9);f4=polyfit(v4,c4,2);f5=polyfit(v5,c5,2);f6=polyfit(v6,c6,2);%用最小二乘法得出风速的瓦斯浓度影响函数系数f7=polyfit(v7,c7,2);f8=polyfit(v8,c8,2);f9=polyfit(v9,c9,2);fg4=polyfit(v4,M4,2);fg5=polyfit(v5,M5,2);fg6=polyfit(v6,M6,2);%用最小二乘法得出风速的煤尘浓度影响函数系数fg7=polyfit(v7,M7,2);fg8=polyfit(v8,M8,2);fg9=polyfit(v9,M9,2);l4=poly2str(f4,'v');l5=poly2str(f5,'v');l6=poly2str(f6,'v');%用最小二乘法得出风速的瓦斯浓度影响函数l7=poly2str(f7,'v');l8=poly2str(f8,'v');l9=poly2str(f9,'v');lg4=poly2str(fg4,'v');lg5=poly2str(fg5,'v');lg6=poly2str(fg6,'v');%用最小二乘法得出风速的煤尘浓度影响函数lg7=poly2str(fg7,'v');lg8=poly2str(fg8,'v');lg9=poly2str(fg9,'v');LINGOmodel:min=Q1+Q2+Q3;!目标使总风量最少;s1=4;s2=4;s3=4;s4=4;s6=4;s7=4;s8=5;s9=4-0.1256;A4=3.7945;A6=4.4021;A9=1.2339;r=40;w=0.15;Q1=60*v1*s1;!各个井道需要的风量;Q2=60*v2*s2;Q3=60*v3*s3;Q10=Q3*(1-w);v4=(Q1+A4)/(60*s4);v6=(Q2+A6)/(60*s6);v7=(Q2+A6)*0.85/(60*s7);v8=(Q3*(1-w)+Q10+A9+Q1+A4+Q2+A6)*(1-w)/(60*s8);v9=(Q10+A9)/(60*s9);c4=0.11013*v4^2-0.78536*v4+1.9101;!有MATLAB的的关系函数;c5=0.029394*v5^2-0.35587*v5+1.3427;c6=3.6264*v6^2-15.6863*v6+17.8211;c7=6.1184*v7^2-26.562*v7+29.7329;c8=0.26158*v8^2-2.8424*v8+8.3301;c9=-0.23669*v9^2+0.98382*v9-0.78487;m4=0.35661*v4^2+0.28382*v4+5.2326;m5=3.26*v5^2-12.5501*v5+19.4258;m6=20.4598*v6^2-84.6745*v6+95.2166;m7=25.5506*v7^2-107.4708*v7+120.2322;m8=0.56646*v8^2-5.6671*v8+21.2051;m9=1.2367*v9^2-3.6172*v9+9.2355;r4=(0.0738*c4^2-0.5234*c4+0.9791)*r;r5=(0.0738*c5^2-0.5234*c5+0.9791)*r;r6=(0.0738*c6^2-0.5234*c6+0.9791)*r;r7=(0.0738*c7^2-0.5234*c7+0.9791)*r;r8=(0.0738*c8^2-0.5234*c8+0.9791)*r;r9=(0.0738*c9^2-0.5234*c9+0.9791)*r;Q10>150;!局部通风机额定功率风速;Q10<400;v1>0.25;v1<6;!各个井道风速的约束条件;v2>0.25;v2<6;v3>0.25;v3<6;v4>0.25;v4<4;v5>0.25;v5<4;v6>0.25;v6<4;v7>0.25;v7<6;v8<8;v9>0.25;v9<4;c4<0.78;!瓦斯浓度关系函数小于给出数据中最大的瓦斯浓度;c5<0.83;c6<1.11;c7<1.18;c8<0.71;c9<0.33;m4<r4;!煤尘浓度关系函数小于煤尘爆炸下限函数;m5<r5;m6<r6;m7<r7;m8<r8;m9<r9;end附： 2006高教社杯全国大学生数学建模竞赛题目（请先阅读“对论文格式的统一要求”）D题:煤矿瓦斯和煤尘的监测与控制关键环节（见附件。瓦斯是一种无毒、无色、无味的可燃气体，其主要成分是甲烷，在矿井中它通常从煤岩裂缝中涌出。瓦斯爆炸需要三个条件：空气中瓦斯达到一定的浓度；足够的氧气；一定温度的引火源。2000g/m3，而当矿井空气中瓦斯浓度增加时，会使煤尘爆炸下限降低，结果如附表1(2)规程要求煤矿必须安装完善的通风系统和瓦斯自动监控系统，所有的采煤工作面、掘进面和回风巷都要安装甲烷传感器，每个传感器都与地面控制中心相连，当井下瓦斯浓度超标时，控制中心将自动切断电源，停止采煤作业，人员撤离采煤现场。具体内容见附件2的第二章和第三章。1是有两个采煤工作面和一个掘进工作面的矿井通风系统示意图，请你结合附表2按照煤矿开采的实际情况研究下列问题：根据《煤矿安全规程》第一百三十三条的分类标准见附件224页斯矿井”还是“高瓦斯矿井1，判断该煤矿不安全的程度（生爆炸事故的可能性）有多大？为了保障安全生产，利用两个可控风门调节各采煤工作面的风量，通过一个局部通风机和风筒实现掘进巷的通风（见下面的注。根据附图1所示各井巷风量的分流情况、对各井巷中风速的要求（《煤矿安全规程》第一百零一条页)，以及瓦斯和煤尘等因素的影响，确定该煤矿所需要的最佳（总（。注 掘进巷需要安装局部通风机，其额定风量一般为150~400。局部通风机所在的巷道中至少需要有15%的余裕风量（新鲜风）才能保证风在巷道中的正常流动，否则可能会出现负压导致乏风逆流，即局部通风机将乏风吸入并送至掘进工作面。名词解释()。掘进巷：用爆破或机械等方法开凿出的地下巷道，用以准备新的采煤区和采煤工作面。掘进工作面：掘进巷尽头的开掘现场。新鲜风：不含瓦斯和煤尘等有害物质的风流。乏风：含有一定浓度的瓦斯和煤尘等有害物质的风流。附表1:瓦斯浓度与煤尘爆炸下限浓度关系采煤工作面Ⅱ风筒采煤工作面Ⅱ风筒局部通风机掘进巷道放大示意图待采煤层掘进工作面进风巷主巷道II回风巷II掘进工作面进风巷I回风巷I总进风巷进风口总回风巷出风口待采煤层待采煤层图例局部通风机17风量调节风门（进风）采煤工作面Ⅰ风）空气中瓦斯浓度(%)00.51.01.52.02.53.03~53.54煤尘爆炸下限浓度(g/m3) 30~50 22.5~37.5 15~25 10.5~17.5 6.5~12.5 4.5~7.52.5~3.5 1.5~2.5PAGEPAGE22附图1某煤矿的通风系统示意图图注:主巷道断面大约为5m2，其他各采煤区的进风巷、回风巷和掘进巷的断面大约为4m2，掘进巷道中的风筒直径为400mm。一类手机资费方案的优选模型摘要手机现已成为人们用于通信的重要工具之一，手机资费问题一直是人们关注的焦点，研究手机资费方Maltab软件进行求解得出结论，为用户提供根据自己的情况选择最优资费方案的依据。1，2得出了同种套餐中的不同服务类型之间的关系，从而得出了不同的用户选择不同资费方案。就问题2，评价准则为在全球通、神州行和动感地带三种服务中用户在各服务时间相同的情况下，用户所花费的资金为越少越好，即用户的满意度越高越好，运用综合评价的方法对各种资费方案进行评价，得出了各种类型的用户所占权重系数的大小。根据附表1所给的数据总的时间进行了处理，运用求平均值法，求出了不同类型的用户对不同资费方案的满意度。3，22就问题4，基于用户所处地的经济水平、用户的稳定性、用户的不同需求、公司的利益的因素，根据北京、上海现推行的套餐，对以上因素都考虑的较全面，所以直接在原有基础上对利润进行调整，得出了更使用户满意的方案。2关键词：资费方案；综合评价；数学模型（注：此文获2007年全国大学生数学建模竞赛全国一等奖）一、问题重述随着我国经济的飞速发展，手机现已成为人们日常生活、工作等社会活动中必备的工具之一，同时，通信业务量也随之飞速增长。手机自费问题一直是人们关心的热点话题，经过长时间的研究资费方案始终没有实质性变化。但是20071各种品牌得套餐，手机套餐得花样琳琅满目，让人眼花缭乱。人们不禁要问：手机套餐究竟优惠几何？9968案（见附表2，建立数学模型分析研究一下问题：合的套餐方案并说明理由。资费标准作分析、比较，并给出评价。对全球通“被叫全免计划”方案进行评价并给出理由。根据你们可能会考虑到的因素，设计一个合理的套餐方案。二、问题分析问题一]据移动电话资费标准和基本费，其他服务不考虑。分别对同种“套餐”中的不同类型服务进行两两比较，以后者的月基本费为相比较的临界值，可以通过找出它们两两之间的关系，得出各种类型的用户所适用的“套餐”方案。问题二对各种资费方案提出评价准则和方法，将北京、上海“套餐”与现行的资费标准进行分析、比较且给出评价。资费评价准则是指在全球通、神州行和动感地带三种服务中用户在各服务时间相同的情况下，用户所花费的资金为越少越好，即用户的满意度越高越好。根据题目中所给附表1根据不同类型的用户对手机所使用的功能要求不同，即有的用户手机通话量大；有的用户手机短信量较大等。根据用户对手机不同功能的需求，不同类型的用户对通话量和短信发送量的需求不一样。因此，引入偏大型柯西函数确定它们所占有的权系数值，从而来确定不同类型的用户每月通话量和短信量。用户自己的资费方案（全球通、神州行和动感地带。问题三对新推出的全球通“被叫全免费计划”套餐进行评价，用户必需在网一年，即交一年的月租，考虑该型套餐的月租为第一个月交清，即第一个月交600元月租，接打电话按单项收费。“被叫全免费”是第一个月交的费用多，以后每月费用为是资费标准方案的一半。比较两种类型收费相等时月份是否大于12个月，即整年内用户的使用时间若大于两者收费相等时的月数，则选择被叫全免费，若使用时间小于两者收费相等时的使用月数，则选用资费标准方案。问题四要为移动公司设计一个全球通的资费方案，就应该综合考虑移动公司和用户两方面的利益，移动公司在在资费方案达到用户的满意情况下获取利益，所以应从多种相关因素考虑制定套餐方案。计时都考虑了能注意到的因素，所以认为合理的套餐方案时在9968成本的不变性设计套餐。三、符号说明T——各套餐中第i种形式的本地主叫的通话时间；i,4minia——第i种套餐中用户在本地的免费主叫通话时间；i1,i

,4minb——超出套餐的主叫时间；i,4元miniQ——所使用不同套餐所需缴纳的费用；j1,2,3;i1,ji四、模型假设

,4元1、用户若使用某种资费方式则使用期至少一个月，未满一个月按一个月计费；2、套餐的免费时间不可积累至下个月，未用完的免费时间不返还给用户；3、只考虑用户使用某一种资费方式的卡，不可同时使用两种或两种以上的卡。五、模型建立与求解问题一假设各套餐中第i种形式的本地主叫通话量为T，单位为min；国内IP长途通话时间为T

，单位i 2i为min；移动的月基本费为A，单位为元min；超出套餐的主叫时间为b，单位为元min；套i i餐中免费主叫通话时间为a，单位为min；则资费的计算方法为i本地主叫通话资费A为ATi

0Ta

（1）i iAi ATab长途通话资费f为T

i i iiTif0.1TT

Tai i（2）模型的建立北京移动公司全球通“畅听99由移动电话资费的费用标准和附表中所给的数据，得出同种套餐中不同类型的收费方式分别为99Q990.T

为0

288Q

i

（3）i 99

t2881i

i

t2881i当用户选择该套餐中139元的套餐时，所需缴纳的费用Q为2i1390.T

0t

560Q

i

（4）2i 139

t5601i

i

t5601i当用户选择该套餐中199元的套餐时，所需缴纳的费用Q为3i1990.T

0t

1000Q3i 199

it10001i

i

1it10001i

（5）当用户选择该套餐中299元的套餐时，所需缴纳的费用Q为2i2990.T

0t

2000Q2990.25i20004i20004i

2000t

（6）当用户拨打本地电话的费用小于Q

i 中每月规定所要收取的费用时，就选择该套餐服务。若当用户拨打本地电话的费用Q服务，

恰好等于或大于Qj1i

中每月规定所要收取的费用时，就在下个月选择Qj1i

的套餐当用户在本月所选择的套餐服务是Q即

服务时，用户所缴纳的费用恰好为下一套餐的最少费用Q时，2iQQ1i 2i

（7）当用户在本月所选择的套餐服务是Q2i即

服务时，用户所缴纳的费用恰好为下一套餐的最少费用Q时，3iQQ2i 3i

（8）当用户在本月所选择的套餐服务是Q即

服务时，用户所缴纳的费用恰好为下一套餐的最少费用Q时，4iQQ3i 4i

（9）68由移动电话资费的费用标准和附表中所给的数据，得出同种套餐中不同类型的收费方式分别为当用户选择该套餐中68元的套餐时，所需缴纳的费用Q为j168Q

0t

（10）j1 680.18

t360j1

t 360j1当用户选择该套餐中139元的套餐时，所需缴纳的费用Q为2i128Q

0t j

（11）j1 128

t360j1

t 800j2当用户选择该套餐中199元的套餐时，所需缴纳的费用Q为3i188Q

0t j

（12）j3 188

t j3

t 1200j3当用户拨打本地电话的费用小于Q

中每月规定所要收取的费用时，就选择该套餐服务。若当用户拨打本地电话的费用Qji恰好等于或大于Q

Q

的套餐服务，

ji1

ji1当用户在本月所选择的套餐服务是Q服务时，用户所缴纳的费用恰好为下一套餐的最少费用Q时，j1 j2即QQ （13）j1 j2当用户在本月所选择的套餐服务是Q服务时，用户所缴纳的费用恰好为下一套餐的最少费用Q时，j2 j3即Q Q （14）j2 j3模型的求解北京移动公司全球通“畅听99Matlabt394min；t1i

800min；t3i

1500min；Q

与通话时间所要缴纳的费用Q2i

间的临界值t394min；通话时间所要缴纳的费用Q1i

与通话时间所要缴纳的费用Q

间的临界值t2i

800min；通话时间所要缴纳的费用Q

与通话时间所要缴纳的费用Q4i

间的临界值t3i

1500min。每个用户每月的本地主叫时间t不相同，则有i当0ti

394min时，用户选用的99套餐为99元的套餐服务；当394ti当800t

800min时，用户选用的99套餐为139元的套餐服务；1500min时，用户选用的99套餐为199元的套餐服务；i当1500ti

min时，用户选用的99套餐为299元的套餐服务。北京99套餐费用与通话时间的关系如图1。图1北京99套餐费用与通话时间建议：当用户的通话时间恰好为它们的临界值时，选用后者的套餐服务。2.68Matlabt 694min；t 1175min；j1 j2套餐中不同的类型两两相比较，每月通话时间所要缴纳的费用Q与通话时间所要缴纳的费用Q间j1 j2的临界值t 394min；每月通话时间所要缴纳的费用Q与通话时间所要缴纳的费用Q间的临界值j1 j2 j3t 1175min；每个用户每月的本地主叫时间t不相同，则有j2 j当0tj

694min时，用户选用的68套餐为68元的套餐服务；当694tj

1175min时，用户选用的68套餐为128元的套餐服务；当1175tj

min时，用户选用的68套餐为188元的套餐服务；建议：当用户的通话时间恰好为它们的临界值时，选用后者的套餐服务。上海68套餐费用与通话时间的关系如图2。23图2上海68套餐费用与通话时间的关系问题二综合评价模型建立由评价准则在全球通、神州行和动感地带三种服务中用户在各服务时间相同的情况下，用户所花费的资金为越少越好，即用户的满意度越高越好。对全球通、神州行、动感地带三种服务，每个月所需交纳的费用的大小，来判断用户的满意度，满意度可以相应地评分集为{很满意，满意，较满意，不满意}对应的数值分别为5，4，3，2，1。根据实际情况取偏大型柯西分布属函数 1x2

1x3fxalnxb

3x

（15）其中；ab为待定常数。1，f51；0.8f30.8当评价“不满意”时，隶属度取0.01，即f10.01；fff代入可以推导出1.10860.8942a0.3915，b0.3699；分别代入（15）f20.524524PAGEPAGE27较满意，不满意}的量化值为1,0.5245,0.01；对属函数进行归一化处理分别得出 1 ；1 10.52450.010.01

0.5245 ；10.52450.010.01 3

0.01 ；10.52450.010.01不同的用户对手机的功能不相同，因此可以把不同的用户大致地分为一下四种类型：①发短信为主，偶尔打电话型；②打电话为主，偶尔发短信型；③偶尔打电话，偶尔发短信型；④经常打电话，经常发短信型。则可以从两个方面考虑进行加权处理,即10.010.5245；用户以发短信为主，四种用户类型所占的权重分别为类型①用户所占的权重为1；类型②用户所占的权重为0.01；类型③用户所占的权重为0.01；类型④用户所占的权重为0.5245。12006R和各通话量时长t11

之间的比值，求出平均每个月的通话时长t，即

tt11 （16）tR把2006年移动总用户的人数平均分成4等分，即RR （17）4由权重系数与平均人数R之间的乘积进行求和，得iQ4R （1）ii1偏重打电话型的用户每月的通话时间w为iw T （19）12Q12006Rx11求出每个用户平均每个月的发送短信量x，即

之间的比值，xx11 （20）xR由权重系数与平均人数R之间的乘积进行求和，得iQ4ii1偏重打电话型的用户每月的发短信量w为j

x （21）w X （22）12Q问题二的模型求解通话费的求解从附件1中可得到移动电话方面200646108.215882.7977.84.24.0IP1465.6亿分钟，国外IP12.3亿分钟，港澳台IP10.7亿分钟。按照模型中用户对打电话和发短信的不同偏重将总用户分为四类，据分布函数求得的满意指标，以及移动各方面的总通话时间，用Matlab求得各种类型的用户各方面的每月总的平均通话时间如表1。表1 各类用户各方面的每月平均通话时间（单位：min）本地 国

国外主叫长途主叫长途长途长途长途台偏重打电话743.424744.22340.1966 0.187268.60060.57570.5008偏重发短信7.434200.44220.0020 0.00190.686000.00580.0050都不偏重7.434200.44220.0020 0.00190.686000.00580.0050

IP

IP

IP港澳（续表1）都偏重 都偏重 389.9263 23.19520.10310.098235.9810.30200.2627从附件中可得到5种卡型对各方面得通话收费如表2。表2 五种卡型各方面通话收费（单位：/分）本地 国内 国外 IP国内IP国外IP港澳主叫长途长途港澳台长途长途台全球通0.417.52.00.73.91.9神舟行0.61.38.62.60.92.14.1动感地带0.490.897.50.590.294.121.7299套餐————7.52.00.12.11.968套餐——1.07.52.00.73.91.999套餐和68Matlab求得各类用户使用各种卡时月总费用如表3。表3各用户使用各种卡的月总通话费（单位：元/月）全球通神舟行动感地带99套餐68套餐偏重打电话839.31921141.5410.6442152.7684225.2897偏重发短信57.893211.41454.106499.119168.9729都不偏重57.893211.41454.1064105.279168.9729都偏重463.9979598.6914215.3829145.2446135.9484短信费的求解考虑实际情况可规定经常发短信每个用户每个月发短信500300条，两方面都不偏重的每一用户每月发短信6010条。同样可以2599680.1元/0.15元/条，动感地带按四种用户发短信条数多少选择短信套餐。据Matlab可求得各种用户使用各种卡的每个月总的短信费如表4。表4各用户使用各种卡的月总短信费（单位：元/月）偏重打电话

全球通1.0

神舟行1.5

动感地带15

991.0

68套餐1.0偏重发短信偏重发短信都不偏重都偏重506.030759.045301520506.030506.030对各种服务进行评价月总话费即为当月通话费加短信费，各用户使用各种卡的月总话费如表5。表5 各用户使用各种卡的月总话费（单位：月）全球通神舟行动感地带99套餐68套餐偏重打电话840.31921143.0425.6442153.7684226.2897偏重发短信107.893286.414534.1064149.1191118.9729都不偏重63.893220.414519.1064111.279174.9729都偏重493.9979643.6914235.3829175.2446165.9484比较结果数据得出评价结论，当用户经常打电话偶尔发短信时，若在北京、上海则对应的选用套餐，若在其它地方则选用全球通；当用户经常发短信偶尔打电话时，选用动感地带；当用户打电话和发短信都很少时，选用神州行，当用户常打电话也常发短信时，若在北京、上海则对应的选用套餐，若在其它地方则选用全球通。问题三5012＝600i种用户用该卡一年前n个月共缴纳的费用y为ni

ynqni

600 （23）每个月另行的缴纳的费用q为iq(wi

50)/2

(n1, ,12;i1, ,4) （24）其它资费方案前n个月总共缴纳的费用y(k)nwni

(k1, ,5;i1, ,4;n1, ,12) （25）wik

为第i种用户每个月用第k种卡的费用。仍据第二问的评价准则进行评价，与北京9968较，以月为单位作出各种卡的增长曲线如图。图3被叫全面与其它比较观察图象与第二问的评价结果进行比较，得出结论，从第二个月开始新推出的“被叫全免费”就比全球通的更省钱。问题四模型建立若为移动公司设计全球通的资费方案，我们会考虑的因素有公司的利益、当地的经济水平、顾客的需求、顾客源的稳定型等。将超出套餐的主叫费看成两部分分别为利润和成本，每个层次的套餐都是确定的，成本与利润的关系是国定的为25％，通过这些来确定每种套餐的成本费x0.25xyi i

（26）其中，x为成本；y为超出部分每分钟主叫资费；i i求出每种套餐时间内的每分钟的单价

Q50P ii ai

（27）其中，Q为月基本费；a为套餐免费主叫时间；P各种套餐的单价；i i i求出每种套餐时间内的每分钟的利润w为iwPx

（28）i i i确定制定套餐的现行利润w(现利润为原来的90％)为i28PAGEPAGE36w90%w

（29）其中，现制定资费标准的利润；确定制定套餐的单价P为i

i iPxw

（30）其中，现制定的套餐时间内的单价；确定套餐的月基本费Q为i

i i iQPt

（31）超出套餐部分主叫资费y为i

i i iyx

0.250.9

（32）i i i其中，y为现制定的套餐超出部分主叫资费。i模型求解根据Matlab求出北京的新套餐资费如表6，上海的新套餐资费如表7。原套餐单价原套餐利润新套餐利润新套餐单价原套餐单价原套餐利润新套餐利润新套餐单价新套餐月基本 新套餐超出部费分主叫费0.53210.25210.22690.506991.940.14700.33750.13750.12380.3238131.30.19600.24900.08900.08010.2401190.10.24500.17450.5450.04900.1690288.10.3430表7 新设计的上海套餐资费表新套餐月基本新套餐超出部套餐单价 原套餐利润 新套餐利润 新套餐单价费 分主叫费0.32780.09430.16540.188961.38400.12740.22250.09450.08510.1231120.440.15680.19830.18380.08490.3094176.680.1764原六、结果分析原问题一中由北京套餐的各种形式的分别求出两两之间的临界值分别为394min800min1500min这里只作为用户选择套餐的一种参考的估计，因为在计算是忽略了GPRS12590*，12586*等通话费，求得的结果与实际临界值有差别。根据实际中用户每月的打电话时间并相同，即用户对自己大电话时间也是估计值，所以此方法虽做不到完全准确，仍具有可行性。问题二，得到各种资费方式各类用户每月的总话费比较可作为用户选择资费方式，同问一中的数据结果与实际相符合。结果显示，对打电话和发短信都不偏重的用户神州行月费用为18.42元，动感地带的月总费用为19.104元。即动感地带的月总费用只比神州行的大了0.7元，导致评价后果并不是很准确。七、模型改进针对问题二中求得的月总费用差别不大，评价不能保证准确，所以改为用资费来计算，根据问题一种给出的资费计算方法分别计算出各种资费方式本地通话部分、长途部分、漫游部分、短信部分的资费。p,3),3)

（i9968）p(q

xt)/ti i

(i1,(50s)fp i i

50s(tii

wi

(1f)

(i4,5)1i 2ti

2t ii其中q为第i种卡的月租费，x为第i种卡的每分钟资费，s为第i种卡的月基本费，i i iw为第ifi

0为使用第i种卡超出套餐时间，fi

1为使用第i卡未超出套餐时间。长途部分资费为p2i

p (1.1t50)/t21p 1.322p (151.42t)/t23p (500.7t)/t24p,5)3i,5)

qpiq

mti1i

nti2i

(i1,3i tt2i其中m、n分别为使用第i种卡漫游地当地资费、漫游地以外资费。i ip4ip 0.1，p41

0.15，p43

15/n，p41

0.1，p42

0.15八、模型评价与推广模型的优点为通话费和短信费两部分，先分开后总体计算费用对各种费用方案进行评价；为单种形式资费，具有实用性；对问二、问三结合现实将用户分为四类，具有普遍适用性。模型的缺点1资费，会造成最总结果费用的不准确，但不会影响最终的评价结果；2过实际调查对四种不同用户给出短信量更合理；3．GPRS12590，12586结果费用。模型的推广1动、电信和联通等各种通信业务；2．3．推广到公路、铁路运输费用的选择。九、参考文献[1]苏龙.移动电话合理资费研究[J].工商管理硕士学位论文，2001，附录问题一的Matlabclearclcsymst1=[0,280,2500];q1=[99,99,876];t2=[0,560,2500];q2=[139,139,624];t3=[0,1000,2500];q3=[199,199,499];t4=[0,2000,2500];q4=[299,299,376];plot(t1,q1,'r',t2,q2,'--',t3,q3,'-.',t4,q4,':')title('北京套餐时间与月总费关系')xlabel('时间'),ylabel('月总费')legend('99','139','199','299clearclcsymst1=[0,360,1500];q1=[68,68,273];t2=[0,800,1500];q2=[128,128,240];t3=[0,1200,1500];q3=[188,188,227];plot(t1,q1,'r',t2,q2,':',t3,q3,'-.')title('上海套餐时间与月总费关系')xlabel('时间'),ylabel('月总费')legend('68','128','188clearclcsymst1=[0,280,394.3];q1=[99,99,139];t2=[394,560,800];q2=[139,139,199];t3=[800,1000,1500];q3=[199,199,299];t4=[1500,2000,2500];q4=[299,299,376];plot(t1,q1,'r',t2,q2,'--',t3,q3,'-.',t4,q4,':')title('北京套餐时间与资费关系')xlabel('时间'),ylabel('资费')legend('99','139','199','299clearclct1=[0,360,693];q1=[68,68,128];t2=[693,800,1175];q2=[128,128,188];t3=[1175,1200,1500];q3=[188,188,227];plot(t1,q1,'r',t2,q2,':',t3,q3,'-.')title('上海套餐时间与资费关系')xlabel('时间'),ylabel('资费')legend('68','128','188问题二的Matlabclearclcb1=15882.7*4/4.61082/12/(0.01+0.5245+0.01+1);n1=944.8*4/4.61082/12/(0.01+0.5245+0.01+1);w1=4.2*4/4.61082/12/(0.01+0.5245+0.01+1);g1=4.0*4/4.61082/12/(0.01+0.5245+0.01+1);pn1=1465.6*4/4.61082/12/(0.01+0.5245+0.01+1);pw1=12.3*4/4.61082/12/(0.01+0.5245+0.01+1);pg1=10.7*4/4.61082/12/(0.01+0.5245+0.01+1);x1=[b1,n1,w1,g1,pn1,pw1,pg1]%偏重打电话的用户各方面通话时间jiahe=b1+n1+w1+g1+pn1+pw1+pg1b2=15882.7*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;n2=944.8*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;w2=4.2*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;g2=4.0*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;pn2=1465.6*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;pw2=12.3*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;pg2=10.7*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;x2=[b2,n2,w2,g2,pn2,pw2,pg2];jiahe2=b2+n2+w2+g2+pn2+pw2+pg2b3=15882.7*4/4.61082/12/(0.01+0.5245+0.01+1)*0.01;n