B纯净水安全监控问题的研究与应用
摘要：本文我们根据某地区提供的近年来关于各公司的纯净水检测报告，考虑电导率、菌落总数、大肠困群、霉菌四个危害因素对纯净水的安全风险评价的影响建立了数学模型，并利用此模型对纯净水安全监控问题做进一步探讨。
针对问题一，提出简单的线性综合评价模型。将数据按各批次进行不同因素的加权均值处理，再无量纲化后，利用变异系法求得权向量
,再根据本文已知建立特征矩阵，得到权向量
，采用综合法得到四个危险指标的权向量
。从而确定纯净水安全风险度模型.
利用此模型对所有批次的纯净水进行安全风险评估。
针对问题二，将数据按各公司进行不同因素的加权均值处理，利用模型一对该城区的九家水公司按风险度进行排序评价，
；然后再按各公司的不同因素的加权均值，分类排序评价得表六，并分析表六中的排序名次（排在前面的危害因素危害程度更大）得出各公司相对其他公司而言的主要可能危害因素。依据公司主要危害因素、各公司的不同因素的加权均值数据、公司风险度、饮用水安全方面有关知识，指出公司的实际特点。
针对问题三，我们先对纯净水生产流通环节即仓库和销售网点两类的不合格数进行统计计算出不合格危害因素在仓库的比例占76.19%，销售网点占3.15%，作直方图，可以直观的看出仓库比销售点存在更大的隐患。利用excel对各批次不合格危害因素进行统计。再用Excel将表一的无量纲化后的数据绘制成图表，得到四个危害指标随着检验批次的变化图。根据图表总结分析，我们得出A公司在电导率方面存在严重问题，B在不断完善中，D公司在各方面危害指标的管理存在严重问题。
针对问题四，我们运用了规划模型，运用分层求解法求得了公司的规模划分问题，再将将A公司的各点的危害因素进行加权平均，然后进行无量纲化的处理，运用模型一求出各存放点的安全风险度值进而得出A公司各存放点批次的分布。
针对问题五，我们根据以上建立的模型分析出了此城区饮用水的在各项危害指标中都存在隐患，并提出了5条相关性的监控意见。
【关键词】纯净水 安全风险评价 安全监督 线性综合评价模型 变异系数法 层次分析法 综合权向量法 规划 安全形势 监控措施
一、问题的重述
日趋加剧的水污染，已对人类的生存安全构成重大威胁，成为人类健康、经济和社会可持续发展的重大障碍。据世界权威机构调查，在发展中国家，各类疾病有8%是由于饮用了不卫生的水而传播的，每年因饮用不卫生水至少造成全球2000万人死亡，因此，水污染被称作"世界头号杀手"。
我国政府对纯净水安全问题十分重视，已将纯净水安全作为一项重要的公共管理目标，采取了一系列措施，强化纯净水安全的监管，并取得了初步成效。但纯净水安全问题的总体形势仍不容乐观，依然存在一系列隐忧，近年来食品安全方面的恶性、突发性事件屡屡发生。20##年07月12日，南通一纯净水厂发生造假事件。20##年3月底，贵阳市发生数百人感染甲肝事件，经卫生部中国疾控中心专家组核查，确认“竹源牌”桶装水是造成疫情爆发的主因。20##年03月25日，某大学B区学生饮用了“清清”牌桶装纯净水后，百余学生先后出现集体腹泻事件。20##年2月26日，湖南师范某寝室在长沙爱高普纯净水有限公司订购的桶装纯净水中出现了黑色虫子事件。生物性和化学性污染对纯净水安全的影响愈来愈严重。
本问题主要考虑纯净水的以下危害因素: （按照危害的严重性依次给出）
“电导率”: 是纯净水的特征性指标，反映的是纯净水的纯净程度，以及生产工艺控制的好坏，“电导率”根本达不到国家卫生标准要求，与自来水无异，根本不能算做纯净水。
菌落总数: 是指纯净水检样经过处理，在一定条件下培养后所取1ml（g）检样中所含菌落的总数。它可以作为判定纯净水被污染程度的指标之一。
大肠菌群：反映纯净水加工过程中对大便污染程度的一个指标。数值越高证明污染越严重。
霉菌：食物霉变后产生，直接引起中毒，或产生致癌物质，毒害人体。
纯净水的安全危机的爆发，往往是日常的监控机制和管理长期存在漏洞的反映。完整、有效的纯净水安全风险分析监测预控，为政府及有关部门实施控制措施提供决策依据和技术支持，可以有效提高纯净水安全监管效率和管理水平，及时化解可能出现的安全危机。近年来，我国在从国家宏观层面探讨建立纯净水安全预警机制的研究方面，已取得了不少理论成果但由于我国地域辽阔，经济社会发展水平很不平衡，如何构建有效的预警机制并应用到饮用水安全监控过程还处于起步阶段。
某城区共有九家生产并销售纯净水的公司，其中A公司和B公司规模较大，其余均为小公司。针对该城区提供的近年的关于各公司的纯净水检测报告（见附件），请你利用数学建模的方法回答以下问题：
1、结合本问题所给数据，给出纯净水安全风险分析的科学评价方法，确定评价的标准和评价的规则，对该城区所有批次的纯净水进行评判排序。
2、对该城区范围内的监控对象（各公司）按风险度进行排序评价，并对它们分类综合评价，指出各公司产品的主要可能的危害因素，并指出同类公司的实际特点。
3. 对检测出的不合格的样品成因分析：评价纯净水生产流通环节（归为仓库和销售网点两类）的危害因素以及各个危害指标的分布规律，并通过四类危害指标的分析，讨论A、B、D公司的管理状况。
4.国家相关部门每年要面对各种专项检验，对于纯净水专项检验的投入经费有限，已知该城区下一年度投在纯净水方面的检验总批次为100个批次，在现有历史数据的基础上，并考虑各公司的实际运行状况，如何设置各公司检验批次的分布，使得抽检方案的针对性最优（即检出的风险性为最大）。
5.结合你的工作，请你给该城区食品安全委员会写一篇短文，阐述你的观点，评价该城区的饮用水安全形势并给出监控对策。
二、模型假设
(1)各个样品有一个危害因素不合格则视为该样品不合格；
(2)采样的地点、样品和数量均是随机的；
(3)各公司的规模与被检验的总批次数量成正比；
(4)假设题目所给数据真实有效；
(5)本文只考虑所检测的四个危害因素，其他因素忽略；
(6)各个公司在一年内的生产模式基本保持不变；
(7)各个公司放在仓库和销售点的产品相同；
(8)本批次检验的产品均是待出售的现货或存货；
(9)在一年内没有被检查，则视为该批次在这一年为合格。
三、符号说明
：在第k时间时第j批次的样品经项目i检验的结果；
：在第k时间时第j批次的样品经项目i检验的样品数量；
：第j批次的样品经项目i检验的结果；
：数据经过无量纲化处理后对应的数值；
：i因素的权重；
：权重向量；
: i地点抽检了n批次；
: 纯净水的安全风险度；
: 市场占有率
:各公司的风险度
:各公司的批次数
四、模型的分析与准备
文中已给出纯净水检测的数据，对此数据进行分析。考虑到以产品的批次来分类，所以需对原始数据进行简单处理。其中，对霉菌和酵母一项“没有检验出”时看作0来处理。按批次进行平均加权：
由于各组数据的单位不同，需对其进行无量纲化处理 ：
得到yij。相关数据见下表：
Tab.1各批次不同因素的加权平均
五、模型的建立与求解
5.1模型一的建立与求解（问题一）
本模型我们考虑最简单的线性综合模型：
我们在原题所给数据的基础上，利用变异系数法计算各项指标的权重:
其中：
是第
项指标的变异系数、也称为标准差系数；
是第项指标的标准差；
是第项指标的平均数。
各项指标的权重为：
由题目所给数据，利用以上公式和无量纲化后的数据计算各指标的权重系数如下表：
Tab.2各指标的权重系数
，此权向量完全是由数据计算而得到，客观性较强。但是，由于确定各项指标的权重需要综合考虑评价对象的属性机理等因素，因而不能仅仅依赖数学计算处理后的数据，必须结合各项因素的社会属性等，给出另一个主观性较强的权重向量。这一点题目中给我们提供了各项指标的重要性程度，因此我们可以利用层次分析法，根据“1-9尺度”原则，建立成对比较矩阵 ：
经一致性检验，CR=CI/RI=0.0433<0.1（程序见附录 road1.m）,所以 A 的特征向量可用作权向量。
从而得到另一个权重向量：
(程序见附录road2.m)。此权向量由生活中的经验常理等得出，主观性较强，与 ?1可以互补。
最后综合考虑这两方面我们可以得到权重向量 ? 。利用综合法，计算 ? ：
? =??1+（1??）?2
其中0≤a≤1，为调整主客观因素比例最终决定各目标权重系数的因子，它可以根据不同系统和实际需要而调整。
此处，我们取a=0.6。得到 ? =
(程序见附录road3.m)
。将 ? 代入（3）式可得模型一：
Q表示安全风险度，Q值越大，此纯净水的质量越不安全。由此，我们可以对各个批次的纯净水进行评判排序。数据如下表：
Tab.3所有批次的纯净水安全风险评判排序
从上表可以看出，批次20070512的纯净水的安全风险最小，最安全；批次20090111的纯净水安全风险最大，最不安全。
进一步分析：结合表一可以看出，批次20080902的水在电导率、霉菌和酵母菌两项检验上严重超标，水质量严重不合格。批次20090110的水在菌落总数一项上超标及其严重，水质量也严重不合格。考虑到每一批次中，有一项检验不合格，则此水质量即为不合格。以上被检验的所有批次中合格的水有11批20070415，20070425，20070512，20070609，20090109，20070810，20071216，20071218，20071220，20080826，20070809。其他9个批次均不合格。
5.2模型二的建立与求解（问题二）
5.2.1公司的安全风险综合评价
Tab.4各公司不同因素的加权平均
将量纲化后的数值代入（5）式可得各个公司的值(程序见附录road4.m)，以此来排名：
Tab.5所有公司的纯净水安全风险评判排序
由上表可以看出，公司G的纯净水安全风险最小；公司D的纯净水安全风险最大，最不安全。统计各个公司各批次在抽查中有不合格的项目次数，得出：
Tab.6各个公司各批次不合格项目次数统计
分析此表可得：C、E、F、G公司纯净水质量最好，不合格项目为0记录。公司A，D不合格记录最多，其次是公司H、I，B较少。
5.2.2公司分类的评价
同样再对各个公司分类进行评价。根据表三的数据，我们分别按这四个危害因素将这九个公司进行排序评价：
Tab.7各公司分类分析评价表
上表中，排名靠前的公司在此项目上的安全风险度越大。从表中可以看出，各个公司在各项的评价排名上有很大区别。其中显而易见，A、D两家公司多次排名靠前，产品安全风险较大，不安全。而E、F、G三家公司排名较后，产品较安全。
对表格七的数据进一步分析：可以找出各公司产品的主要可能的危害因素，如下表：
Tab.8各个公司产品主要可能的危害因素
由上可以看出，有些公司有多种危害因素，但只按照其危害程度相对其他公司较大的来计。有些公司虽然没有检验出不合格的产品，但是它的某项指标的数值较大，比较接近标准值，是潜在的危害。
5.2.3公司的实际特点分析
结合表八，做进一步的分析可以找出同类公司的实际特点。利用文中给出的对电导率菌落总数大肠杆菌及霉菌和酵母菌的定义和解释，结合我们查找的有关饮用水安全等方面的知识，对这九家公司的实际特点分别作出说明。结论见下表
Tab.9公司的实际特点
5.3模型三的建立与求解（问题三）
5.3.1生产环节的评价
以样品为单位，进行不合格产品的成因分析。数据中共有9个公司，20个批次和70个样本被检验。按照纯净水生产流通环节来分可分为仓库和销售网点两类。对其不合格样品数量进行统计分析，可得以下数据：
Tab.10不合格率分类统计表
由上表可以清楚地看出：位于仓库的水不合格率高达76.19%；相比之下，销售点的产品质量要好些，不合格率只有3.57%；可见销售点水的质量问题比较严重，可能是时间过长。
由上表数据将两者进行对比做出下表：
Fig.1仓库与销售店的不合格率的比较
由上图可以看出：销售网点的水少部门在电导率方面存在问题，其他方面还好；而仓库的纯净水在电导率和菌落总数方面存在着很严重的问题，不合格率很高，而且在霉菌和酵母这一项上有较大问题；不过在大肠杆菌上两个都还好，没有不合格的问题。所以可知，导电率偏高是各个地方普遍存在的问题。
5.3.2危害指标的分析
我们把各个不合格因素进行分类，按照检验的批次顺序，将各个危害指标中不合格的公司进行汇总。得出的结果如下表：
Tab.11各批次不合格危害因素统计表
由上表可以的到各个危害因素在各个公司的不合格批次产品中的数量分布，但由不合格次数不足以判断出完整的各个危害因素的分布规律，还需结合不合格的程度进行综合衡量才行。
对于表一中的数据，将无量纲化后数据绘制成图标，得到四个危害指标随着检验批次及公司的变化图。图中的横坐标是以批次的顺序即以时间的变化的。所以此图可以看出四个危害指标随时间的变化的大体走向。
Tab.12各公司检验次数
Fig.2检验指标变化趋势图
在上图中，有四种不同颜色的实折线和虚直线，当实折线超过相同颜色的虚直线代表超标，检验不合格。各点即为检测数据。结合表九可以得出以下规律，针对四个危害因素分别给出：
（1）电导率的分布规律
结合表九，从图中的蓝线可以看出，从不合格样本数量来看，电导率偏高的纯净水样本主要分布在A公司，其次是H和I，D公司最少。其他公司在此项没有检验出不合格。从具体数值来看，D厂只抽检了一批电导率不合格的纯净水，但数值高达84.4，若认为抽到电脑率极高的纯净水是小概率事件，在不考虑特殊原因的情况下，恰恰在D公司发生了，所以可以认为D公司此批次大多数的纯净水电导率都没控制好，都有严重的问题。而D公司其他批次的水没有检测出问题，说明D公司在总体上来看，在此项上问题不是很大。
结论：电导率这项危害主要分布在A公司，其次是H、I。20##年此危害较少存在，但自20##年此项危害出现率较高。
（2）菌落总数的分布规律
结合表九，从图中的红线可以看出菌落总数不合格的情况在A、B、D、H、I公司都存在，而A、B、D、H公司都有严重超标现象。其中，同（1）中的D公司类似，B公司只检验出一批菌落总数超标的纯净水，但该值高达800，说明B公司在此批次的纯净水在菌落总数上没有控制好，导致大都出现非常严重的问题，而B公司总体上在此项没太大问题。
结论：菌落总数这项危害因素主要分布在A、H、D公司，其次是I。此危害在20##年有一个极大值，但自20##年来得到较好的控制。
（3）大肠杆菌的分布规律
结合表九，从图中可以看出，检验中并没有发现此项超标的情况，说明总体上纯净水在此方面没有太大问题。但是除了E、C、G非常干净外，其他的所有是都存在不同程度的轻微污染。此危害指标程度不定，随时间的波动性较大。
（4）霉菌和酵母的分布规律
结合表九，从图中的紫线可以看出，霉菌和酵母超标只出现在A、D两公司，由由上表数据可以看出该危害因素主要分布在D公司，其次是A。此危害指标随时间出现间歇性超标现象
5.3.3A、B、D公司的管理状况
由上面的分析A、B、D厂生产的纯净水分别在电导率、菌落总数及霉菌和酵母这三项危害因素上有严重问题。
A公司生产的纯净水在三年的纯净程度均不是很理想，纯净水受到不同程度的污染，该企业的管理状况不容乐观，但在20##年生产的纯净水受污染程度明显减少，说明管理得到了一定的改善。但在导电率方面非常差，生产工艺控制方面的缺陷依旧严重，急需改善。
B公司生产的纯净水于2007在生产流程中受到了一定的污染，但20##年以来菌落总数已得到很好的控制，其他方面该公司做得很好，所以企业的管理在不断的完善之中。
D公司生产的纯净水在电导率、菌落总数、霉菌和酵母这三项指标上都有严重的问题，而且3年以来这些问题一直存在，说明该企业在管理方面有很大的缺陷，生产工艺控制方面做得很糟，其产品对人的健康是一较大的安全隐患。
5.4模型四的建立与求解（问题四）
5.4.1按公司规模划分
此问是典型的规划问题，目的就是使各批次检出的总风险最大，即抽检出的各个公司的水总风险性最大，我们这里运用分层求解的办法求解。
首先把100批次分配到各个公司。题目说9个公司中A公司和B公司规模较大，其余均是小公司，所以在分配抽检批次时，必须考虑此项因素，这里我们做一下粗略的划分，参照前面对各公司的抽检次数将各个公司规模划分如下：
在第二问中我们依旧求解出各个公司的风险度（见第二问求解），分配时综合考虑市场占有率和风险度两因素：
Tab.13公司市场占有率统计表
设市场占有率为，风险度为(数据见表四)，则
=*,
因此当批次总数为100时，分配到各公司的批次数为:=100
结果见下表：
Tab.14各公司检验批次的分布
5.4.2对公司批次进行划分
首先我们根据数据表发现各公司的不同时期风险度是动态变化的，如图所示：
Fig.3 A公司各批次仓库纯净水的风险度
风险度随着时间的变化先增加后减少，说明该公司在风险度低时会疏忽对生产工艺等方面的监督，而使风险度变大，当风险度逐步增大时因为城区监管部门的检验有逐步减少到极小值，如此周期循环下去，符合我们的一般经验。这里取它的极大值作为它的风险度最大值，极小值为它的风险度最小值。
我们根据题目所给数据，将A公司的各点的危害因素进行加权平均，然后进行无量纲化的处理，如下：
再根据安全风险度公式 求得各点的安全风险度值，具体数据如下表：
将上表中计算得到的Q值，进行归一化处理，：
从而得到A公司各存放点的市场占有率，进而得到A公司各点的批次分布情况，如下：
同理，我们根据以上方法，分别求得B、C、D、E、F、G、H、I公司各存放点的批次的分布，如下表：
5.5 给食品安全委员会的一封信
尊敬的食品委员会的一封信：
您好！
最近我们对贵城区的纯净水的安全形势做了研究分析，对纯净水的安全风险评判建立了数学模型，利用此模型对纯净水安全监控问题提出一些建议，并对城区九家纯净水公司的产品卫生安全情况作了分析评价。
针对该城区近三年来九个纯净水公司的产品质量检测报告，我们经过建模分析发现纯净水中危害因素电导率超标现象越来越普遍，问题很严重；菌落总数出现间歇性超标，但其总体有下降趋势；菌落和霉菌控制方面也一直控制得不好，问题一直存在。同时我们对城区的九家纯净水公司做了安全风险评判，其结果为，G公司产品最安全，其排名依次为
。按照目前趋势如果不采取措施，贵城区的饮用水安全环境会逐渐降低，居民的用水安全将受到严重威胁。
为了解决此问题我们小组提出几条纯净水安全监控对策：
1、缩短对于葛根素生产纯净水的检测周期，检验各公司生产纯净水的水源的各项指标，然后找到各公司的相对薄弱环节并进行科学的指导与改进，并加强薄弱指标的检查力度；
2、在生物性和化学性污染日趋严重的大环境下，提倡各公司强化传统生产水的工艺、替换传统的消毒剂、吸附、膜过滤和生物预处理等深度净水技术
3、为了提高居民对自身饮用水的安全意识，对每期纯净水抽样检验结果进行公布，并将各公司进行排名，促进各公司之间的竞争，带动纯净水行业的进步，加大排名靠后的公司的检查力度；
4、加强各公司对人员培训和生产过程的监督与管理。安全优质的引用水不仅取决于优质的水源、优质的净水工艺以及安全的供水网络，也取决于制水人员的卫生意识以及对水质变化的应急处理能力。因此需加强水厂制水人员、管理人员及监督人员的专业训练；
5、对于多次检查不合格的且不接受科学技术的公司依照法律进行停产整治，保障人民饮水安全的保障
上面是我们小组在数学建模中获得的结果与想法，希望能给您们带来工作上的帮助。
六、模型的评价与改进
6.1 模型的评价
优点：
（1） 用合理的方式处理不同类型的数据，使其尽量符合实际情况
（2） 在权重系数的确定中，运用变异系数和层次分析法相结合的形式，变异系数法仅仅用数据求权重，客观性较大；而层次分析法利用了各个指标的重要性程度，主观性强，两者结合使问题的考虑更全面更符合实际
（3） 在问题求解中运用了大量的表格和图，使结果明了、清晰
（4） 在求解最有分配方案时，我们利用了层次分析，巧妙的把复杂的问题转化为单个的线性的规划问题，简化了计算过程。
缺点：
（1） 对每个批号的检验结果取平均值，精确度不高
（2） 在确定建模一的权重时，虽然运用了综合权重法，但尚有不足之处不是很精确
6.1.2模型的改进
对于我们给出的权重模型，在考虑权重时，我们可以再利用几种权重分析法，例如主成分分析法，负相关系数法等，最后再综合各权重向量得到一个比较贴近实际的权重分配，这样可以增加模型的可信度。
对于纯净水安全问题，它会随着人的主观意志而改变，不会是静态的，也不是很确定，单一模型很难满足实际情况，故我们可以建立一个动态的或多种机制综合来衡量。
七、参考文献
[1]姜启源,谢金星,叶俊,数学模型（第三版）.高等教育出版社2003.8
[2]赵静,但琦,数学建模与数学实验（第三版）. 高等教育出版社2008.1
[3]陈杰,MATLAB宝典.电子工业出版社2007.1
八、附录
程序road1.m
%层次分析法求权重
clear;
clc
a=[1 3 5 7
1/3 1 3 5
0.2 1/3 1 3
1/7 0.2 1/3 1];
[x,y]=eig(a)
eigenvalue=diag(y);
lamda=eigenvalue(1)
ci1=(lamda-4)/3;cr1=ci1/0.9
w1=x(:,1)/sum(x(:,1))
程序road2.m
%迭代求权重
clc;
clear;
A=[1 3 5 7
1/3 1 3 5
0.5 1/3 1 3
1/7 0.5 1/3 1];
[V,B]=eig(A)
e0=[0.25 0.25 0.25 0.25]';
b=A*e0;
x=b(1,:)+b(2,:)+b(3,:)+b(4,:);
e1=b/x
m=A*e1;
x=0;
for n=1:4
x=x+m(n,:);
end
e2=m/x
m=A*e2;
x=0;
for n=1:4
x=x+m(n,:);
end
e3=m/x
m=A*e3;
x=0;
for n=1:4
x=x+m(n,:);
end
e4=m/x
m=A*e4;
x=0;
for n=1:4
x=x+m(n,:);
end
e5=m/x
m=A*e5;
x=0;
for n=1:4
x=x+m(n,:);
end
e6=m/x
m=A*e6;
x=0;
for n=1:4
x=x+m(n,:);
end
e7=m/x
m=A*e7;
x=0;
for n=1:4
x=x+m(n,:);
end
e8=m/x
程序road3.m
%安全风险度计算
clear;
clc;
W1=[0.1548 0.3314 0.1703 0.3435]';
W2=[0.5289 0.2558 0.1477 0.0676]';
w=0.6*W1+0.4*W2
y=[
0.058098 0 0.5 0
0.038633 0.003125 0 0
0.014859 0.009375 0 0
0.025557 0.015625 1 0
0.107429 0.109375 0.5 0
0.080684 1 0.25 0.25
0.05951 0.0015625 0.5 0
0.070282 0.00625 0.25 0
0.054978 0.00521875 0.1665 0
0.061516 0.0015625 0 0
0.40416 0.003125 0 0
0.050074 0 0 0
0.067756 0 0 0
0.122511 0.05 0.125 0.25
0.055721 0.10625 0 0
0.316493 0.140625 0 0
0.027043 0 0 0
0.413076 0 0 0
1 0.0234375 0 1
0.381575 0.01875 0.2 0
];
Q=y*w
程序road4.m
%公司安全风险值
clear;
clc;
W1=[0.1548 0.3314 0.1703 0.3435]';
W2=[0.5289 0.2558 0.1477 0.0676]';
w=0.6*W1+0.4*W2
y=[
0.7145 0.2087 0.1902 0.1305
0.1440 1.0000 0.3281 0.0000
0.1003 0.0217 1.0000 0.0000
0.8676 0.6472 0.4375 1.0000
0.1244 0.0066 0.0000 0.0000
0.1871 0.0000 0.8750 0.0000
0.0478 0.0198 0.0000 0.0000
0.8034 0.5920 0.1750 0.0000
1.0000 0.2022 0.2916 0.0000 ];
Q=y*w