煤矿瓦斯和煤尘的检测与控制_大学生数学建模竞赛论文.doc

高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛的题目是： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：1高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：煤矿瓦斯和煤尘的检测与控制摘 要 随着社会的发展，煤矿安全事故也在不断的增加。本文就某一个煤矿的具体情况，对煤矿中瓦斯和煤尘爆炸进行分析，并建立模型解决几个典型问题。 问题一：按照国家煤矿安全规程中给出的评价标准，先对两个关键的指标绝对瓦斯涌出量和相对瓦斯涌出量进行合理的定义。运用附表2中给出的数据求出该煤矿的绝对瓦斯涌出量，由于相对瓦斯涌出量与每天工作时间有关，所以根据临界的相对瓦斯涌出量求出临界时间。由临界时间来判断该煤矿是属于“高瓦斯煤矿”还是“低瓦斯煤矿”。用excel求解得到：当该煤矿的工作时间大于10.105时，该煤矿属于“高瓦斯煤矿”。当工作时间小于10.105 时，该煤矿属于“低瓦斯煤矿”。问题二：附表2中给出的是安全的瓦斯浓度和煤尘量，所以对附表中的数据按各个工作面进行拟合预测，采用线性拟合和灰色预测的方法预测出该煤矿中各个工作面的速度与瓦斯浓度及速度与煤尘的关系函数。再求出该煤矿的爆炸下限区间，以及资料中给出的条件求出的爆炸下限区间。再用几何概率的方法求出该煤矿发生爆炸的可能性。结果为：采煤工作面i和掘进工作面发生爆炸的概率为0.当时，采煤工作面可能发生爆炸。当值取时对应得值0.250.350.450.550.650.750.850.9510.40678.28306.62375.31044.25803.40482.70522.12521.051.151.251.351.451.551.651.63921.22760.87550.57100.30510.07040问题三：在问题二的基础上，以该煤矿所需要的总通风量最佳为目标函数构建非线性规化模型，在相应的约束条件下运用matlab软件求解，采煤工作面的风量，采煤工作面的风量，局部通风机的额定风量，最佳总通风量。本文巧妙的运用线性规划和灰色理论知识对数据进行拟合并得到了很好的效果。关键词：拟合；灰色预测；非线性规划；煤矿安全规程一、问题的重述煤矿安全生产是我国目前亟待解决的问题之一，做好井下瓦斯和煤尘的监测与控制是实现完全生产的关键环节（见附件1）。瓦斯是无毒、无色、无味的可燃气体，其主要成分是甲烷，在矿井中它通常从煤岩中涌出。瓦斯爆炸需要三个条件：空气中瓦斯达到一定浓度；足够的氧气；一定温度的引火源。煤尘是在煤炭开采过程中产生的可燃性粉尘。煤尘爆炸必须具备三个条件：煤尘本身具有爆炸性；煤尘悬浮于空气中并达到一定的浓度；存在引爆的高温热源。实验证明，一般情况下煤尘的爆炸浓度是30-2000,而当矿井中瓦斯浓度增加时，会使煤尘爆炸下限降低，结果如附表1所示。国家煤矿安全规程给出了煤矿预防瓦斯爆炸的措施和操作规程，以及相应的专业标准（见附录2）。规程要求煤矿必须安装完善的通风系统和瓦斯自动监控系统，所有采煤工作面、掘进面、回风巷都要安装甲烷传感器，每个传感器都与地面传感器相连，当井下瓦斯浓度超标时，控制中心将自动切断电源，停止采煤作业，人员撤离采煤现场。具体内容见附录2的第二章和第三章。附图1是有两个采煤工作面和一个掘进工作面的矿井通风系统示意图，请结合附表2的监测数据，按照煤矿开采的实际情况研究下列问题：（1）根据煤矿安全规程第一百三十三条的分类标准（见附录2），鉴别该矿是属于“低瓦斯矿井”还是“高瓦斯矿井”。（2）根据煤矿安全规程第一百六十八条的规定，并参照附表一，判断该煤矿不安全程度（即发生爆炸事故的可能性）有多大？（3）为了保障安全生产，利用两个可控风门调节各采煤工作面的风量，通过一个局部通风机和风筒实现掘进巷的通风（见下面的注）。根据附图1所示各井巷中风速的要求（见煤矿安全规程第一百零一条），以及瓦斯和煤尘等因素的影响，确定给煤矿所需要的最佳（总）通风量，以及各个采煤工作面所需要的风量和局部通风机的额定风量（实际中，井巷可能会出现漏风现象）。注 掘进巷需要安装局部通风机，其额定风量一般为150-400。局部通风机所在巷道中至少需要有15%的余裕风量（新鲜风）才能保证风在巷道中中正常流动，否则可能会出现负压导致乏风逆流，及局部通风机将乏风吸入并送至掘进工作面。 附表1：瓦斯浓度和煤尘爆炸下限浓度关系瓦斯浓度（%）00.511.522.53爆炸下限浓度305022.537.510.57.56.512.56.412.554.57.5 35二、模型假设和符号说明模型假设：1、假设在工作时间内单位时间内煤的产量相同，即用可以用平均产量来表示单位时间的产量；2、假设给出的30天的数据具有普遍性，可以代替整体；3、假设该矿不存在其他的安全隐患，爆炸只考虑瓦斯和煤层的爆炸；4、假设只要瓦斯或煤尘达到一定浓度就会爆炸，不考虑其他爆炸条件；5、假设瓦斯浓度只考虑,不考虑瓦斯中包含的其他物质。6、假设井巷中不漏风；7、假设采煤量影响瓦斯浓度和煤尘的量忽略不计。符号说明：第天第班的绝对瓦斯涌出量（单位：）； ：第天第班的相对瓦斯涌出量（单位：）； ：第天第班总回风巷中得风速（单位：；）； ：第天第班的总回风巷中的瓦斯体积比（%）； ：工作面或巷道的断面的横截面积（单位：）,分别是采煤工作面，采煤工作面，主巷道，掘进工作面;：第天的工作时间（单位：）; ：第天的采煤量（单位：）; ：通风量（单位:）; 分别是采煤工作面、采煤工作面、局部通风机所在的巷道的通风量，局部通风机的额定风量；：风速（单位：；）， 分别是采煤工作面、采煤工作面、局部通风机所在的巷道、局部通风机的风速;：空气中有瓦斯时的煤尘爆炸下限（单位：）;: 降低系数；：当瓦斯浓度为0时，煤尘的爆炸范围；：发生爆炸的概率，时分别为采煤工作面、掘进工作面；:空气中瓦斯的浓度（%），时分别为采煤工作面、掘进工作面；：空气中煤尘的浓度（），时分别为采煤工作面、掘进工作面。三、问题分析问题一是要求我们根据煤矿安全规程第一百三十三条的分类标准鉴别该矿是属于“低瓦斯矿井”还是“高瓦斯矿井”。在该标准中提到了相对瓦斯涌出量和绝对瓦斯涌出量，并且没有给出明确的定义。而我们又要根据这两个指标来判断该矿井的级别。所以我们先从这个标准中的信息对这两个评价指标进行了合理的定义。瓦斯是从煤岩裂缝中涌出的，然后再涌入巷道，最后通过风口排到外界。所谓“涌出量”,从字面上理解即是从一个空间中涌出到外面的量。于是我们将瓦斯涌出量类比于物理中体积流量，将绝对瓦斯涌出量定义为单位时间内通过某一横截面积的瓦斯，单位是。相对瓦斯涌出量定义为每生产1t的煤从煤矿中向外涌出的瓦斯量，单位为。当相对瓦斯涌出量和绝对瓦斯涌出量达到或超过一定值时，则该煤矿属于“高瓦斯煤矿”，于是我们的目的开始转为求解相对瓦斯涌出量和绝对瓦斯涌出量。在计算瓦斯量时，我们是直接算从回风巷中涌出的，还是分别算出各个工作面和回风巷的瓦斯量之和，抑或是都可以。从附图1中分析可以发现，我们只能算从回风巷中涌出的瓦斯量。因为工作面和工作面中的瓦斯都会涌到回风巷中，如果直接加起来的话会出现重复，这一点我们也可以从数据中得到验证。在理清这些之后，就可以对附表2中的数据进行处理，用excel表格进行相应的计算。先计算出绝对瓦斯涌出量，然后再根据绝对瓦斯涌出量算出对应的相对瓦斯涌出量，再运用煤矿安全规程中的分类标准进行分类，从而解决第一个问题。问题二要求我们判断该煤矿的不安全程度（即事故发生的可能性），在煤矿事故中，主要是瓦斯爆炸和煤尘爆炸。而爆炸条件是有一定浓度的瓦斯或煤尘、充足的氧气和一定的温度，从附件1中的资料中可以知道，一般情况下矿井中的氧气以及爆炸的温度都可以达到。所以假设发生爆炸的条件只需考虑瓦斯和煤尘浓度。从附表2中的数据分析可以发现，每天监测的瓦斯浓度和煤尘都远远不足以达到瓦斯爆炸的瓦斯浓度，而有瓦斯存在时煤尘爆炸下限会降低。所以考虑有瓦斯存在的煤尘爆炸，先根据附件1背景资料中的瓦斯浓度和煤尘下限的关系。求出煤尘爆炸下限的区间。然后对附表二中的数据进行处理，附表2中的数据都是安全的数据从表中的数据和我们的生活经验都可以知道，风速越大，巷道或工作面中的瓦斯浓度就越低，而煤尘则会越多。另外采煤量越多，瓦斯和煤尘也越多。从附表2中的数据可以知道，采煤量对瓦斯浓度和煤尘的影响不大，所以我们可以假设采煤量不影响瓦斯浓度和煤尘，只考虑速度的影响。因此可以通过拟合的方法求出个工作面上瓦斯浓度与风速以及煤尘与风速之间的函数关系。然后预测出当不同风速时爆炸的瓦斯浓度和煤尘量。根据拟合的曲线求出该煤矿的煤尘爆炸下限。然后用几何概率的方法，则求得的该煤矿的爆炸下限的区间长度与根据资料中的下限公式求出的下限区间的长度的比值即是该煤矿的爆炸可能性。问题三是根据附图1中各井巷风量的分流情况、对各井巷的要求，以及瓦斯和煤尘等因素，确定最佳总通风量，以及两个采煤工作面所需要的风量和局部通风机的额定风量。实际中，井巷可能会出现漏风现象，由于漏风是很难估计的，与很多因素有关。而且相对于总风量，漏风只是少数，可以忽略不计。所以假设井巷不漏风。在煤矿中，降低瓦斯浓度的唯一办法是通过通风系统将瓦斯排除到井外，从降低瓦斯浓度的角度考虑的话，巷道中风速越大越好。但是风速越大，煤尘又会增加。而煤尘增加同样会爆炸。所以要通过控制一定量的通风量来使得瓦斯浓度和煤尘达到一定的平衡使其不发生爆炸。在达到这个条件的情况下，总风量越小越好。所以以总风量最小为目标函数，在根据题目中的条件确定约束条件，再用lingo软件进行求解从而解决第三个问题。四、模型建立与求解4.1 模型一的建立与求解对于问题的求解，我们直接根据煤矿安全规程第一百三十三条的分类标准进行判断和求解。前面我们对绝对瓦斯涌出量和相对瓦斯涌出量已经进行了定性定义，先用表达式表示如下：绝对瓦斯涌出量： 相对瓦斯涌出量： 是指对应工作面或巷道的横截面积，分别是工作面，回风巷，工作面，回风巷，掘进工作面，总回风巷。是指第天第班的工作时间。是指第天第班的采煤量。经过分析，我们只考虑从总回风巷中的绝对瓦斯涌出量和相对瓦斯涌出量，用表示。分类标准是：当时，该矿井属于“低瓦斯矿井”；当时，该矿井属于“高瓦斯矿井”。根据前面的定义式，可以得出第天第班的绝对瓦斯涌出量表达式：根据附表2中的数据运用excel可以求出，由于每天早，中晚班的值不同。又因为我们考虑的是安全问题，所以我们选择每天三个班次中最大的。三十天的绝对瓦斯量最大的值如下表：表4.1.1 三十天的绝对瓦斯涌出量日期123456710.351510.348810.587610.54689.808510.342810.1568日期8910111213149.69579.979210.1019.75369.10810.526410.8054日期151617181920219.78369.98410.48569.876610.60810.454410.0584日期2223242526272810.041610.975210.291210.43259.88210.09929.6579日期293010.86310.1985由上表可知，该矿井的绝对瓦斯涌出量全都小于。=： 表示第天的最大绝对瓦斯涌出量；：第天早中晚班绝对瓦斯涌出量的最大值。由于与每天的工作时间有关，所以不能直接判断，因此要鉴别该矿属于“低瓦斯矿井”还是“高瓦斯矿井”，即令=10。可以求出临界时间，将原先的评价标准转化为用时间来评价该矿的属性。根据附表2中的数据，用excel求出每天的临界时间如下：对这些数据取平均值，得到。所以当该煤矿每天的工作时间时，该煤矿属于“高瓦斯煤矿”； 当时，该煤矿属于“低瓦斯煤矿”。4.2.1模型建立与求解针对问题，根据煤矿安全规程第一百六十八条的规定，并参照附表1中的数据，经分析知该煤矿的瓦斯浓度达到最大时也只在1%左右，而根据附件1中在新鲜空气中瓦斯爆炸界限一般为，虽然在实际矿井中，空气中的含氧量略低，但该矿井中的瓦斯浓度仍无法 达到其爆炸界限。设为空气中有瓦斯时的煤尘爆炸下限，为煤尘的爆炸下限，一般为且降低系数为。由于瓦斯参与煤尘爆炸下限降低，当瓦斯浓度低于4%时，煤尘的爆炸下限.根据表1：瓦斯浓度对煤尘爆炸下限的影响系数 空气中的瓦斯深度% 0 0.50 0.75 1.0 1.50 2.0 3.0 4.0 k 1 0.75 0.60 0.50 0.35 0.25 0.1 0.05 测试结果表明，随着空气中的瓦斯浓度的增加，其降低系数k减小，为了确定它们之间的回归函数的类型，我们将这8组数据作为坐标在平面直角坐标系中描出它们相应的点，即散点图（见下图1） （ 图1 ）从图中可以看出，所有散点大体上散布在一条指数形式的曲线上，因此估计大致是指数的形式，即认为与具有如下关系：其中是常数.由于与之间的关系不是线性关系，我们通过变量代换把非线性回归化成线性回归.作变量代换： .对空气中瓦斯浓度（%）进行独立观测，得到n个空气中的瓦斯浓度及对应的，把看作随机变量表示成：相互独立.对未知参数的估计，利用最小二乘法：求使为最小。利用求极值的方法来求，令整理得：求得的估计值为：其中 ，.用这两个估计值得回归方程：利用matlab中的统计工具箱（程序见附录1，所得图形如上图1）可以得到回归系数及其置信区间（置信水平），检验统计量值的结果如下表：参数参数估计值置信区间 0.0668-0.0089 0.1424-0.7655-0.8027 -0.7283由检验统计量可知有明显的线性关系。通过以上分析得下降系数与瓦斯浓度关系为：由于该煤矿每个部分都可能发生爆炸事故，现分别对采煤工作面i，采煤工作面ii，掘进工作面进行分析。为了分析空气中瓦斯浓度，空气中煤尘浓度m与风速v之间的关系，需要对附件1的数据进行拟合。经过分析发现，附件1中的风速v有许多重复值，为了提高准确度，对应的取平均值，整理的数据（见附录2）.利用整理后的数据分别作出对v的散点图如下：（图2：采煤工作面1，w与v拟合图） （图3：采煤工作面1，m与v拟合图） （图4：采煤工作面2，w与v散点图） （图5：采煤工作面2，m与v散点图） （图6：掘进工作面，w与v散点图） （图7：掘进工作面，m与v拟合图）从图2，3，7可以发现，随着v的增加，w,m的值有比较明显的线性关系，用线性模型 直接利用matlab统计工具箱中的命令regress求解（程序见附录3），得到的回归系数估计值及其置信区间（置信水平），检验统计量值的结果如下表：对于采煤工作面1，w与v参数参数估计值置信区间 1.24041.0969 1.3840-0.2416-0.3019 -0.1814则回归方程为：对于采煤工作面1，m与v参数参数估计值置信区间 3.23722.6671 3.80721.97991.7405 2.2193回归方程为：对于掘进工作面1，m与v参数参数估计值置信区间 2.99552.3162 3.67481.93871.7405 2.2409回归方程为：4.2.2 灰色预测拟合从图4,5,6可以发现，随着v的增加，w,m的值没有明显的线性关系，所以不能用线性拟合的方法拟合。并且数据量也不多，正好符合灰色预测的条件，所以用灰色预测的方法进行拟合。一、工作面的瓦斯与速度的拟合关系工作面瓦斯与速度对应1.010.8830.8630.8780.8690.8670.8720.890.792.062.072.082.092102.112.122.132.14由于速度是以0.01的值增加的，我们可以运用灰色理论的预测方法，记时为第一次，即，则。第一步：级比检验建立煤尘数据时间序列如下： ；求级比：判断级比：由于所有的都在即内，所以可以用做满意的模型。第二步：建模 对原始数据进行一次累加，可得： 其中： 构造数据矩阵和数据向量: 计算的值： 可得： 建立模型： 求解得时间响应式： 求生成数列值及模型还原值：令由上面的时间响应式可得：由得：第四步：模型检验模型各种检验指标的计算结果如下表： 表1 模型检验表序号速度原始值模型值残差相对误差级比偏差12.061.011.010022.070.8830.8862-0.02990.0036-0.135532.080.8630.8762-0.00370.0194-0.015742.090.8780.8734-0.00250.00530.024252.100.8690.867-0.00130.0023-0.003062.110.8670.86070.00690.00720.005072.120.8720.85450.00110.02010.013082.130.890.84830.00230.04690.027392.140.790.8421-0.02660.0660-0.1184对表中数据进行分析，当残差小于0.2时，则达到一般要求。当残差小于0.1时则达到较高要求。而通过这种方法拟合的曲线残差远小于0.1，另外相对误差也非常低，级比偏差小于0.1时达到较高要求。上表中的级比偏差也远小于0.1.所以我们验证得到该拟合曲线的符合程度非常高。计算程序见附录4。将代入时间响应式，可以求出以前的所有瓦斯浓度。我们就可以得到速度为时的瓦斯浓度，即二、工作面速度与煤尘的关系：用同样的方法可以拟合出工作面的速度与煤尘之间的关系，具体如下：级比检验：原始数据： 级比 满足级比检验。可以直接建立模型，运用matlab软件进行求解（程序见附录4）可以得到：检验表如下：表2 模型检验表序号速度原始值模型值残差相对误差级比偏差12.067.667.660022.077.597.58340.00660.00360.000932.087.617.60140.00860.01940.001142.097.627.61950.00050.00530.000152.107.6597.63770.02130.0023-0.002862.117.617.6559-0.04590.00720.006072.127.6337.6741-0.04110.02010.005482.137.7287.69230.03570.04690.004692.147.7257.71070.01430.06600.0019可以得到工作面2中速度与煤尘的关系式：三、掘进巷速度与瓦斯浓度关系拟合掘进巷道中的速度与瓦斯浓度，由于掘进巷道中的速度不是完全由0.01的值每次增加的，但是大部分是按这个规律。所以我们补充一些数据，使得其符合灰色预测的条件，。原始数据：x0=（0.235 0.215 0.24 0.24 0.24 0.247 0.24 0.24 0.21 0.226 0.29 0.27 0.25 0.23 0.22 0.18 0.23 0.24 0.25 0.23 0.22 0.22 0.22 0.22 0.22 0.2 0.2）级比检验：其中，则；把上面得到的级比进行分析发现，这组数据不满足级比检验，所以我们对该组数据进行处理，在原数据的基础上使每一个数据加1,。即：所以，x0=（1.235 1.215 1.24 1.24 1.24 1.247 1.24 1.24 1.21 1.226 1.29 1.27 1.25 1.23 1.22 1.18 1.23 1.24 1.25 1.23 1.22 1.22 1.22 1.22 1.22 1.2 1.2）对应的级比如下：级比验证得： 所以可以用建模。用matlab软件进行求解，程序见附录4。求得结果如下： 部分检验表如下：表三 模型检验表序号速度模型值残差相对误差级比偏差12.091.2351.2350022.101.2151.245-0.030.0247-0.015532.111.241.2438-0.00380.00310.021142.121.241.2427-0.00270.00210.00152.131.241.2403-0.00150.0012-0.00162.141.2471.23910.00670.00540.006672.151.241.23790.0090.0007-0.004782.161.241.23670.00210.00170.00192.171.211.2356-0.02670.00221-0.0238说明：完整表格见附录3。经过检验，精度较高。拟合度较高，可以用拟合的式子来预测。4.2.3 问题求解 设为空气中瓦斯浓度，为空气中煤尘浓度，其中分别对应采煤工作面i，采煤工作面ii，掘进工作面。由公式其中一般为。利用几何概率求解，即发生爆炸事故的概率：其中的关系式如上由煤矿安全规程第一百零一条表2井巷中的允许风速，即。利用matlab软件求解得：，即煤工作面i，掘进工作面发生的爆炸的概率为0.当时，即可能发生爆炸。我们给出值取时对应得值列表如下：0.250.350.450.550.650.750.850.9510.40678.28306.62375.31044.25803.40482.70522.12521.051.151.251.351.451.551.651.63921.22760.87550.57100.30510.07040根据附件2所给数据，采煤工作面的风速都大于1.65，所以爆炸的概率为0。4.3 问题三的求解针对问题，设分别为采煤工作面i，采煤工作面ii，局部通风机所在巷道的通风道,为局部通风机的额定风量。定义：通风量为单位时间通过某一横截面的体积，即在安全生产和其他约束条件的前提下，使该煤矿所需要总通风量达到最佳，即最小。故在此以总通风量最少为目标，以其他的要求作为约束建立目标函数， 即：由于通风量所以将该目标函数转化为：根据附图1所示各井巷风量的分流情况及煤矿安全规程第一百零一条中对各井巷中风速的要求；采煤工作面，掘进中的煤巷和半煤岩巷的允许风流速度最大为最小为，即风速的约束：由题目中局部通风机的额定分量一般为即：又局部通风机所在的巷道中至少需要有的余裕分量（新鲜风）来保证风在巷道中的正常流动，避免出现负压导致乏风逆流，则：由于煤矿中的瓦斯和煤尘等因素的影响，风速过小虽煤尘数比较小然瓦斯浓度会太大而爆炸，风速过大虽瓦斯浓度比较小但煤尘数会太大而达到煤尘爆炸下限，所以每个采煤工作面及掘进巷的瓦斯浓度、煤尘数都要收到相应约束。根据 煤矿安全规程第一百六十八条中对各采煤工作面、掘进巷瓦斯的警告浓度要求：瓦斯浓度应不超过1%，再根据附表1：瓦斯浓度与煤尘爆炸下限浓度关系空气中瓦斯浓度(%)00.51.01.52.02.53.03.54煤尘爆炸下限浓度(g/m3)305022.537.5152510.517.56.512.54.57.5352.53.51.52.5为了安全起见，我们把断电情况下的瓦斯浓度1.5%对应的煤尘爆炸下限作为煤尘约束。利用问题中通过分析得到的瓦斯浓度、煤尘数和风速之间的函数关系，那么：对于采煤工作面i：瓦斯约束： 煤尘约束：对于采煤工作面ii： 瓦斯约束： 煤尘约束：对于掘进巷： 瓦斯约束： （注：上式左边减1是因为在数据处理的时候进行了加1处理）煤尘约束：综上所述并整理得非线性优化模型： 用lingo软件进行求解（程序见附录5），结果如下： 最佳总通风量:五、模型评价1、本文从理论上定量的讨论了风速与瓦斯浓度和煤尘浓度之间的关系，并得到了较好的结果；2、模型中通过分析各个量之间的关系，找到了在不同风速下爆炸事故发生的概率；3、模型讨论了该煤矿的日工作时间，能客观的判断煤矿是高瓦斯矿井还是低瓦斯矿井；4、模型中所得的关系式均通过了相关检验，保证了结果的可靠性；5、对于回风巷考虑欠佳，可能对第二问中爆炸事故的概率有一定的影响。参考文献：1刘思峰 灰色系统理论及其应用 上海：科学出版社 2010.5.12姜启源 谢金星 叶俊 数学模型（第三版） 北京：高等教育出版社 2006.53胡剑良 孙晓军 matlab数学实验 北京：高等教育出版社 2006.6.1附 录附录1：clear;clc r=0.0001,0.5,0.75,1.0,1.5,2.0,3.0,4.0; k=1,0.75,0.6,0.5,0.35,0.25,0.1,0.05; plot(r,log(k),ko); title(rk); xlabel(r); hold on rr=ones(8,1),r; b,bint,ri,rint,stats=regress(log(k),rr); b,bint,stats plot(r,b(1)+b(2)*r,k:);hold on text(1.5,0.53,r-k); plot(r,k,k*);hold on plot(r,exp(b(1)+b(2)*r),k); text(1.7,-1.1,r-ln(k); 附录2：采煤工作面1采煤工作面2掘进工作面 风速 瓦斯浓度煤尘风速瓦斯浓度煤尘风速瓦斯浓度煤尘2.210.7337.582.061.017.662.090.2357.062.220.6957.6452.070.8837.592.10.2157.192.240.677.652.080.8637.612.110.247.122.250.77.822.090.8787.622.140.2477.1672.270.7017.6982.10.8697.6592.150.247.192.290.6957.8032.110.8677.612.170.217.082.310.6757.8092.120.8727.6332.180.2267.1842.340.6787.8732.130.897.7282.190.297.162.360.6697.8532.140.797.7252.210.257.32.370.6747.9662.220.237.242.390.6437.8452.230.227.282.40.6687.9662.240.187.412.410.6638.0072.250.237.232.430.638.022.270.257.412.450.648.142.290.227.482.460.678.172.310.227.522.490.648.192.340.227.492.520.648.1432.350.27.512.560.628.332.360.27.622.60.618.422.370.247.72.40.297.622.440.187.73附录3：clear;clcv=2.21 2.22 2.24 2.25 2.27 2.29 2.31 2.34 2.36 2.37 2.39 2.4 2.41 2.43 2.45 2.46 2.49 2.52 2.56 2.6;m=7.58 7.645 7.65 7.82 7.698 7.803 7.809 7.873 7.853 7.966 7.845 7.966 8.007 8.02 8.14 8.17 8.19 8.143 8.33 8.42;plot(v,m,ko);hold onvv=ones(length(v),1),v;b,bint,r,rint,st=regress(m,vv);b,bint,stplot(v,b(1)+b(2)*v,k)title(1,v-m)clear;clcv=2.21 2.22 2.24 2.25 2.27 2.29 2.31 2.34 2.36 2.37 2.39 2.4 2.41 2.43 2.45 2.46 2.49 2.52 2.56 2.6;w=0.733 0.695 0.67 0.7 0.701 0.695 0.675 0.678 0.669 0.674 0.643 0.668 0.663 0.63 0.64 0.67 0.64 0.64 0.62 0.61;plot(v,w,ko);hold onvv=ones(length(v),1),v;b,bint,r,rint,st=regress(w,vv);b,bint,stplot(v,b(1)+b(2)*v,k)title(1,v-w)v=2.09 2.1 2.11 2.14 2.15 2.17 2.18 2.19 2.21 2.22 2.23 2.24 2.25 2.27 2.29 2.31 2.34 2.35 2.36 2.37 2.4 2.44;m=7.06 7.19 7.12 7.167 7.19 7.08 7.184 7.16 7.3 7.24 7.28 7.41 7.23 7.41 7.48 7.52 7.49 7.51 7.62 7.7 7.62 7.73;w=235 215 240 247 240 210 226 290 250 230 220 180 230 250 220 220 220 200 200 240 290 180/1000;%w=cumsum(w);plot(v,w,ko);hold onvv=ones(length(v),1),v;b,bint,r,rint,st=regress(w,vv);b,bint,stplot(v,m,ko);hold onvv=ones(length(v),1),v;b,bint,r,rint,st=regress(m),vv);b,bint,stplot(v,b(1)+b(2)*v,k)title(,v-m)附录4：工作面瓦斯与速度：clc,clearx0=1.01 0.883 0.863 0.878 0.869 0.867 0.872 0.89 0.79n=length(x0);lamda=x0(1:n-1)./x0(2:n)range=sort(lamda)x1=cumsum(x0)for i=2:n z(i)=0.5*(x1(i)+x1(i-1)endb=-z(2:n),ones(n-1,1);y=x0(2:n);u=byx=dsolve(dx+a*x=b,x(0)=x0);x=subs(x,a,b,x0,u(1),u(2),x1(1);yuce1=subs(x,t,0:n-1);digits(6),y=vpa(x)yuce=x0(1),diff(yuce1)epsilon=x0-yuce %计算残差delta=abs(epsilon./x0) %计算相对误差rho=1-(1-0.5*u(1)/(1+0.5*u(1)*lamda %计算级比偏差工作面煤尘与速度clc,clearx0= 7.66 7.59 7.61 7.62 7.659 7.61 7.633 7.728 7.725n=length(x0);lamda=x0(1:n-1)./x0(2:n)range=sort(lamda)x1=cumsum(x0)for i=2:n z(i)=0.5*(x1(i)+x1(i-1)endb=-z(2:n),ones(n-1,1);y=x0(2:n);u=byx=dsolve(dx+a*x=b,x(0)=x0);x=subs(x,a,b,x0,u(1),u(2),x1(1);yuce1=subs(x,t,0:n-1);digits(6),y=vpa(x)yuce=x0(1),diff(yuce1)epsilon=x0-yuce %计算残差delta=abs(epsilon./x0) %计算相对误差rho=1-(1-0.5*u(1)/(1+0.5*u(1)*lamda %计算级比偏差掘进工作面速度与瓦斯：clc,clearx0=1+0.235 0.215 0.24 0.24 0.24 0.247 0.24 0.24 0.21 0.226 0.29 0.27 0.25 0.23 0.22 0.18 0.23 0.24 0.25 0.23 0.22 0.22 0.22 n=length(x0);lamda=x0(1:n-1)./x0(2:n)range=sort(lamda)x1=cumsum(x0)for i=2:n z(i)=0.5*(x1(i)+x1(i-1)endb=-z(2:n),ones(n-1,1);y=x0(2:n);u=byx=dsolve(dx+a*x=b,x(0)=x0);x=subs(x,