2022年认证杯SPSSPRO杯数学建模A题(第一阶段)人员的紧急疏散求解全过程文档及程序

2022年认证杯SPSSPRO杯数学建模

A题 人员的紧急疏散求解

原题再现：

  在过去的几十年里，由于大规模集会活动的数量和规模的增加，紧急疏散的问题变得越来越重要。通过有限宽度的门或狭窄通道进行疏散是最值得关注的情况之一。为了更好地理解各种情况下的人群行为，已经进行了许多实验和数值模拟。门的宽度、人群组成、出口的位置都会影响疏散的效果。组织真人进行实验固然是一种有效的方法，但由于环境和背景的真实性、压力水平、人员密度、样本大小等无法和真实情况一致，所以许多不同的研究结果都是有争议甚至是互相矛盾的。
  基于模型的研究不需要面对实验遇到的困难。但如果模型的假设是不可靠的，它们可能会产生与现实不符的结果。我们需要针对不同的情形和参数进行研究，来了解人群在疏散逃生时可能呈现出的不同动力学行为。
  第一阶段问题：
  1. 人们的经验表明，情况的紧急程度，即人群的焦虑程度对疏散的过程会有明显的影响。有的实验结果显示紧急程度的增加会加快疏散速度，有的实验结果则发现过高的焦虑程度会减慢疏散的过程。我们假设一个内部无障碍的方形房间（或大厅）中有若干人，只有一个门可用于疏散。请你建立合理的数学模型，确定关键因素，研究人群的焦虑程度对疏散效果的影响。请说明哪些关键参数造成了实验结果的差异。
  2. 房间情况同上，由于疏散的过程显然不是匀速的，请研究疏散速率随时间的变化曲线，以及各关键因素将对曲线产生怎样的影响。
  3. 如果在足够大的大厅内有两个有一定间距的出口，经验表明，在人们选择疏散出口的时候，不一定只会选择距离更近的那个，而是预期时间最短的那个。请建立合理的模型，来模拟疏散路线的选择过程。并研究如果门的宽度不同，将会产生怎样的动态过程。如果大厅内的能见度降低，又将对这个过程产生怎样的影响？

整体求解过程概述(摘要)

  人员疏散问题是国内外研究热点与难点，而找出关键因素对人员疏散的影响是技术难题。人员疏散模型是研究人员疏散的重要手段，如何更加准确模拟疏散过程是当前研究的关键。为使人员疏散离散模型更好描述行人微观特性，提出一种改进元胞自动机的人员疏散模型。首先，需要对数据进行预处理，为了简化后续运算量，需要通过量化分析行人间的作用力或建立基于行人速度的移动方法来改进现有离散模型。将社会力模型中力的概念引入到元胞自动机模型和格子气模型中，建立了一种基于行人间、行人与障碍物间相互作用力的规则。
  其次，根据数据分析结果提出一种用于结合社会力与离散模型优点的元胞自动机模型，该模型可以更精确地模拟疏散过程。模拟当行人移动速度大于１ｍ／ｓ时的疏散情景。讨论模型中不同行人最大运动速度在疏散过程中对行人密度分布以及人群形状的影响。尽管上述研究量化了行人间、行人与障碍物间的摩擦力与排斥力，但是模型中行人间产生摩擦力时行人并未真正接触，而且没有给出作用力大小和行人移动距离之间的关系，行人受到的作用力大小仅与移动方向有关，与移动距离无关。同时对疏散模型进行讨论速度因素对疏散结果的影响。
  再者，对建立的模型进行定量评估，为更好地解决上述问题，本文在引入行人社会力模型的基础上对元胞自动机进行修改，提出一种改进的元胞自动机模型。同时，该模型改进了基于速度的移动规则并引入出口选择机制。根据模拟结果，发现该模型不但具有离散模型的特点，同时具备连续模型更真实模拟个体行为的优点。
  最后，通过案例分析检验模型的有效性和可行性分析，发现模型具有一定的可行性。另外针对模型的不足提出了一些优化意见。
  本文主要是在传统模型基础上，将空间划分为更细小网格，令单个行人占据多个网格以便展现行人间微观行为。引入社会力模型中自身驱动力、排斥力、摩擦力的运算规则改进模型，以模拟行人间碰撞现象。通过建立行人出口选择机制并改进移动规则，真实模拟行人动态出口选择过程。疏散过程中考虑了速度因素，能更加具体展现行人不同的移动方式。模拟结果显示，改进的模型能够真实地反映人员疏散的具体过程。

问题分析：

  通过对题目背景的分析，本题要求我们进行合理的数据处理之后建立一个基于改进元胞自动机的建筑室内人员快速疏散模型，本文针对建筑内多障碍物房间楼梯和电梯等典型区域内的人员疏散过程开展研究，重点研究多障碍物房间内人员疏散特性、紧急疏散时楼梯上的行人流流动规律以及楼梯和电梯用于紧急疏散时疏散效率的比较。出口选择行为的相关性、室内疏散路线进行分析，并基于元胞自动机理论建立考虑多障碍物房间内过道因素的疏散模型，并与实验结果对比，利用提出的模型研究不同出口位置和不同内部布局对多障碍物房间疏散效率的影响。
  对高层建筑内的竖直疏散，利用经典的格子气理论建立可以模拟紧急疏散时楼梯上行人流流动的模型，重点考虑楼梯转换平台上人员的转弯行为，以及楼梯入口处的汇流过程。此外，基于改进的楼梯和电梯疏散疏散时间的计算方法，比较不同待疏散规模和不同楼层高度时分别利用电梯、楼梯和混合使用时疏散单层人员和同时疏散整多层人员所需的疏散时间，并探讨混合疏散时利用电梯疏散的最佳人员比例问题。最后，基于贝叶斯分析提出三阶段的高层建筑竖直疏散的人群动态管理框架，并讨论以疏散指示信息为媒介的疏散人员与管理人员之间的信息交互。

模型假设：

  1. 假设出口 1 和 2 风险评估结果相似
  2. 假设路径的不同和人员焦虑程度对疏散速度影响小
  3. 假设大厅能见度对疏散过程无影响
  4. 假设门宽度的大小对疏散时间没有影响
  5. 假设所选情况符合题目要求
  6.假设在数据分析过程均在误差允许的范围内

模型的建立与求解

实验数据的预处理

  人们的经验表明，情况的紧急程度，即人群的焦虑程度对疏散的过程会有明显的影响。有的实验结果显示紧急程度的增加会加快疏散速度，有的实验结果则发现过高的焦虑程度会减慢疏散的过程。我们假设一个内部无障碍的方形房间（或大厅）中有若干人，只有一个门可用于疏散。请你建立合理的数学模型，确定关键因素，研究人群的焦虑程度对疏散效果的影响。请说明哪些关键参数造成了实验结果的差异。
  通过分析题意可知，要研究在内部无障碍的方形大厅里人群的疏散效果，需要明确影响人群疏散的关键因素。对于人群而言，在发生紧急情况下，当门口宽度过小，大厅内部人群数量过多易形成拥挤，进而严重影响人群疏散速度。此外，当大厅内部人群中老年人和小孩数量较多时，同样影响紧急情况下人群疏散速度。综上所述，疏散速度的快慢主要取决于当前状态下的人群数量、当前状态下门口一次性通过的人群数、人群平均年龄，另外，发生灾难的等级会影响到人群焦虑程度，进而影响疏散效果。

考虑过到因素的疏散模型

  元胞自动机(Cellular Automaton)被广泛应用于模拟运动、扩散、蔓延和演化等自然现象，在交通流与行人疏散中也得到了广泛应用，它是一种空间和时间均离散的可扩展的多维系统，模拟空间由一系列的元胞格点(Cell)构成，每个元胞具有有限个状态，利用局部规则对系统内所有的元胞状态进行更新。
  元胞自动机模型并没有严格定义的物理方程或数学函数，整个系统的状态和相互作用均通过局部规则来实现。因此，元胞自动机具有算法简单，计算效率高，容易理解等优点，可用于模拟大规模的人员疏散过程。最早的元胞自动机模型是由 VonNeumann 提出来的，到现在已经有了近 70 年的历史。
  基于元胞自动机理论的人员疏散模型建模过程大致相同，如图所示大体分为四个步骤:即模拟区域离散化，进行网格划分;初始人员分布，并对人员进行编号标识:第三步为状态更新计算:最后为输出结果和疏散过程演示。

  模拟区域离散化元胞自动机将模拟区域离散成一系列的格点，人员在格点间移动以实现运动过程，每一个格点最多可被一人占据。格点的形状常见的有三角形、正方形和六边形，如图 5-2 所示。不同形状的格点各有优缺点，三角形格点邻居数目少但在表达和显示方面比较困难;方形格点直观简单，易于表达，但是在模拟各向同性的现象时无能为力;而六边形格点可以解决各向同性的问题，但是在计算机上的表达较为困难和复杂。最为常见的为正方形格点，本文模型亦选取方形格点对区域进行离散化。

梯度提升树回归模型

  GBDT 又叫 MART，是一种迭代的决策树算法，该算法由多棵决策树组成，所有树的结论累加起来做最终答案。它在被提出之初就和 SVM 一起被认为是泛化能力较强的算法。GBDT 中的树是回归树（不是分类树），GBDT 用来做回归预测，调整后也可以用于分类。GBDT 的思想使其具有天然优势可以发现多种有区分性的特征以及特征组合。业界中，Facebook 使用其来自动发现有效的特征、特征组合作为 LR 模型中的特征，以提高 CTR 预估的准确性

回归树的模型描述

  回归树的总体流程类似于分类树，区别在于回归树的每一个节点都会得一个预测值，以年龄为例，该预测值等于属于这个节点的所有人年龄的平均值。分枝时穷举每个阈值找最好的分割点，但衡量最好的标准不再是最大熵，而是最小化平方误差。也就是被预测出错的人数越多，错的越离谱，平方误差就越大，通过最小化平方误差能够找到最可靠的分枝依据。分枝直到每个叶子节点上人的年龄达到预设的终止条件（如叶子个数上限），若最终叶子节点上人的年龄未达到预设的终止条件，则以该节点上所有人的平均年龄作为该叶子节点的预测年龄。

梯度提升树回归模型结果

  （1）模型参数输出结果

  （2）关键参数对实验结果的影响

  根据图 5-3 研究报告可知，发生灾难的等级对人群疏散效果的影响高达 44%，说明在发生紧急情况时，发生灾难的等级越高，就越能影响人群疏散效果。同样，门的宽度对人群疏散效果的影响高达 44.1%，说明在发生紧急情况时，门的宽度越小，也会严重影响人群疏散效果。此外，在发生紧急情况时，人群的紧张度和平均年龄对疏散效果的影响分别为 2%和 2.6%，而人群的数量对疏散效果的影响为 7.3%，说明相比于人群在灾难面前的恐惧，人群数量越多就越会影响疏散效果。
  （3）模型评估结果

  上表中展示了交叉验证集、训练集和测试集的预测评价指标，通过量化指标来衡量梯度提升树回归模型的预测效果。其中，均方误差是预测值与实际值之差的平方，其取值越小，模型准确度越高；均方根误差是均方误差的平方根，其取值越小，模型准确度越高；平均绝对差是绝对差的平均值，能反映预测值误差的实际情况，其取值越小，模型准确度越高；平均绝对百分比误差是平均绝对误差的变形，它是一个百分比值，其取值越小，模型准确度越高；拟合优度将预测值跟只使用均值的情况下相比，结果越靠近1，模型准确度越高。从这些量化指标树脂大小可以看出梯度提升模型的准确度较高，实验结果比较可靠。

人员疏散速率的影响因素分析

  在人员紧急疏散过程中，疏散初期，出口处比较通畅，出口附近人员迅速向出口靠拢，极易造成小范围堵塞。与锥形瓶实验相似，如图 5-4 所示

  紧急情况下，房间内所有人员均以最快速度向安全出口疏散，当同时抵达出口人员数量大于该出口单位时间最大疏散人员数量时，将形成小范围堵塞，后续人员迅速赶到，使出口处出现滞留现象。此时，因人员冲撞，在出口处形成不规则且连续的弧形堵塞，即形成僵持局面。人员滞留表征线性疏散阶段开始，出口将以最大流动系数疏散建筑内人员，线性疏散示意图如图 5-5 所示

疏散宽度

  假如一个内部无障碍的方形房间(或大厅）中有若干人，只有一个门可用于疏散。当建筑中发生火灾或者晃动时等突发事件时，能否快速疏散建筑内的人员变得非常关键，接下来通过建模分析随着时间的变化其疏散速率的影响因素。
  通过研究表明，人员的疏散行走速度和人员密度，如式（5）所示
              S＝Ｋ（1 - 0.266ρ） （5）
式中：S 为人员的行走速率，ｍ/ｓ；ρ为人员密度，人／ｍ２，0.54＜ρ＜3.8；Ｋ为常
数（根据不同类型疏散通道取值，1.00～1.40）。
  人员到达安全出口所需的步行时间为： T1＝L／S （6）式中：Ｌ为人员从防火分区到安全出口的最大行走距离，ｍ；S 为人员的行走速率，ｍ/s。
  在大厅中人员疏散过程中，人员疏散速率是由周围环境的拥挤程度决定的，数据拟合表明，人员前后距离的影响占支配地位。通过分析研究人员的数据得出了人员疏散移动速率的计算方法，见式（7）：Uj（ρ）＝Ｕｍ（αA＋βB＋γ）A＝1.32－0.82ln（ρ）B＝3.0－0.76ρ（7）:
  当疏散人员较少时，疏散者可以自由奔跑，此时疏散速度上界记为Ｕｍ ，Uj 为拥挤状态下人员疏散速率，α取值范围为 0.25～0.44；β取值范围为 0.014～0.088；γ取值范围为 0.15～0.26，ρ为人员密度。式（7）考虑了疏散时人员左右及前后拥挤的综合效应，对于大厅疏散比较贴合。
  下面研究安全疏散通道的宽度对人员疏散效率的影响，不同安全疏散通道宽度对人员疏散效率有影响。本次模拟试验中，将疏散通道宽度设置为 3～7ｍ，疏散通道宽度每隔 1ｍ递增变化，每次的模拟人数分别为 150 人、200 人、250 人、300 人、350 人、400人。为了获得准确的疏散结果数据，带有相同参数的疏散模拟每组进行 6 次，并取疏散时间的平均值。

  由图 5-6 可见，当疏散通道宽度为 3～5ｍ时，人员疏散时间随疏散通道宽度的增加呈线性递减，而当疏散通道宽度为 5～7ｍ时，其疏散时间大致相同，变化不明显，表明当疏散通道宽度达到 5ｍ这个临界值时，人员疏散时间随疏散通道宽度的变化曲线趋势变缓，说明疏散通道宽度变化的增加对于人员疏散效率的影响降低；在疏散通道宽度一定的情况下，随着疏散人数的增加，人员疏散时间也呈近线性逐渐递增；当人员疏散人数一定时，人员疏散时间与疏散通道宽度呈反比关系。

  对上述典型的现象进行剖析，只有一个安全出口，出口宽度小，疏散人员较少时，疏散区域的人群密度小，行人间的运动受到他人的影响较小，在有限的时间行人可以自主的疏散出去，且不产生拥挤行为:疏散人员较多时，疏散区域的人员密度大，在行人运动一小段时间时，整个人群就陷入了拥堵状态，大量的人员聚集在出口处，而由于宽度小，人流量有限，当宽度相差不大时，人流量也相差不大，最后模拟所得的时间基本相等，如图所示出口宽度为 1m、1.5m 的情况。对于只有一个安全出口，且出口宽度较大的情况下，疏散人员较少时，疏散区域的人流密度小，人员自主运动不受外界的影响;而疏散人员较多时，疏散区域的人流密度大，人员自主运动受外界的影响，人员拥挤严重，出口宽度较大时，人流量大大，人员疏散时间较短;另外一方面，对于同一个人在相同大小的疏散区域内，当出口宽度不同时，人员所行走的距离也不同，如图 5-7 所示，当出口宽度小时，人员行走的距离为 L2，当出口宽度大时，人员行走的距离为 LI，根据三角形三边定律，可得 L2>L1，即出口变宽度，疏散距离相对的减小。

人员密度

  人员疏散速率与人员密度紧密相关，人员密度体现的是一定范围空间内人员的稠密程度，是人员疏散过程的重要空间指标，通过观察人员疏散速率可以发现，在疏散初期阶段，人员分散在展厅内各个位置，人员相互间的间距较大，人群密度较小，不易出现拥挤现象，靠近安全出口的人员首先完成疏散，疏散速率较高；随着时间的增加，当疏散进行到一定阶段时，展厅内不同位置的人员分别通过疏散主要通道、疏散次要通道向安全出口移动，此时汇集人数增加，由于安全出口宽度有限，易形成“瓶颈”，造成人员排队等待、拥挤，甚至堵塞的现象，节点处人员密度较高，疏散速度降低；到疏散最后阶段时，大部分人员已顺利通过安全出口，仅剩少部分人员未撤离，此时人员拥堵程度降低，疏散速率逐渐提高。

  图 5-8 所示,每种工况只有一个安全出口，安全出口的宽度依次为 1m1.5m、2m、2.5m、3m、3.5m、4m、45m、5m，每种工况下均进行 7 次模拟，模拟人数 100 人递增至 400 人，每次增加 50 人。对比出口宽度为 1m、1.5m 和 45m、5m，不管疏散人数任何的变化，两组情况下所用的疏散时间基本相等。当疏散人数为 400 人时，出口宽度为 1m 或 1.5m，疏散时间大致相等，大约为 180s，而出口宽度为 45m 或 5m 时，所用的时间约为 55s，两者的时间差为 125s;而当疏散人数 100 人时，疏散宽度为 1m 或 1.5m，所用的时间大致相等，大约为 48s，而出口宽度为 45m 或 5m 时，所用的时间为 20s，两者的时间差28s。可以得出，人数越少时，疏散时间随宽度的变化不明显，人数越多时，疏散时间随宽度的变化很明显。由图 5-5 可以看出，当疏散宽度相同时，疏散时间随疏散人数的递增而递增，并且接近线性增长。在人数相同的情况下，安全疏散所需要的时间与出口宽度呈反比例关系。图 5-9 给出了 7 种类型特征人群的行走速率如下：

焦虑恐慌心理

  突发事件下群体的恐慌传播是一个互相直接接触的过程，群体性恐慌的传播和疾病范畴中传染病的传播相似，群体的恐慌能够在人群中扩散，传染病也能够在人群中扩散。说明其都具有“传染性”，如“SARS”病毒具有非常强烈的传染性，同时“SARS”造成的社会恐慌也具有传染性，甚至在病毒未侵袭地区也形成了一定程度的社会恐慌。本文借鉴传染病动力学中的相关概念，将传染病传播过程中的疾病传染者、疾病易感者、疾病潜伏者、疾病免疫者引入突发地铁火灾的乘客群体中，将地铁突发火灾中的乘客群体分为恐慌感染者、恐慌易感者、恐慌潜伏者、恐慌免疫者，以此对群体性恐慌心理的传播机理进行分析。
  恐慌感染者是突发地铁火灾时已经受到恐慌心理影响的乘客群体，并具备将恐慌心理传播给其他乘客的能力。恐慌易感者是对恐慌心理没有免疫力，还未受到恐慌情绪传染但很容易受到其影响的乘客群体。恐慌潜伏者是已了解恐慌信息，即情绪表达、肢体语言、言论等，短时间内对信息的真实性无法做出有效处理，等待官方信息发布或权威人士发言来帮助自身对信息正确判断的乘客群体。
  恐慌免疫者是火灾知识储备、火灾经验充足、对突发火灾持积极态度或对官方正面信息发布完全相信的乘客群体，能够有效控制自身情绪，排除恐慌心理的影响，从而对恐慌产生“免疫力”。根据以上对地铁乘客群体的分类，本文提出如下突发地铁火灾时群体性恐慌产生的机理模型。

  为验证本文所提模型的实用性，通过模拟实验进行模型应用分析。选取 14mx9m 的房旧作为疏散环境,设房间有 2 个出口，左上角出口为 A,右下角出口为 B,每个出口大小均为 1m,即每次只允许通过 1 个行人。设疏散人员为 100 人行人质量 65至75kg,半径为0.25m,每个行人实际占据面积为 0.5mx0.5m,行人随机分布在房间中。由以往社会力模型研究结论知:初始行人速度为 0.9~1.5m/s,方向任意,在应急疏散情况下期望速度大小v=2.5m/s。社会力模型的初始参数设置。

模拟疏散过程与分析

  具体模拟过程 5-12 如图所示。疏散开始后，行人逐渐向出口移动。最初，行人根据出口选择机制分别向两出口汇聚，随着出口附近行人密度越来越大,位于后方的行人推挤前方的行人，此时可以观察到行人闻发生碰撞、排斥与摩擦现象，且在出口附近出现了堵塞效应“弓形”结构开始形成。随着后方行人到达出口,出口聚集区边缘受到碰撞挤用的行人会在合力作用下向着周边绕行。人群最终为均匀地分布在出口附近,从而使出口附近行人“弓形”结构更加明显。在疏散后期,出口附近行人逐渐变少，“弓形”结构逐渐减小直至最终消失。

行人出口选择分析

  在整个疏散过程中,各个出口疏散的人数以及疏散的总人数随时间的变化趋势如图5-12 所示结合 5-11 图可以看出，出口附近人群密度和距离效应共同影响行人对出口的选择。两出口附近人员然摩较小,行人主要依据距离效应选择出口(图 5-11(a))。因距离A 出口近的行人较多,在疏散第 13 s 时(图 5-11(b)),行人在 A 出口汇集且密度大于 B出口。尽管考虑出口选择的惯性作用，图 5-11(a)中部分最初选择 A 出口的行人重新选择 B 出口。
  随着 B 出口处汇集人数增加(图 5-11©)两出口队近人群密度近似，距离又成为影响出口选择的主要因素,部分选择 B 出口的行人最终又选择 A 出口作为目的地。在此期间，行人会依据出口附近的人员密度与行人到出口的距离效应不断重新做出选择，从而使两出口附近的人数得到动态调整，待疏散人数大致相等。虽然 A 出口疏散人数较多，但在疏散结束时两出口几乎同时完成人员疏散工作，没有出现一个出口闲置而另一出口仍处于拥堵的情况,说明实验符合真实模拟场景。

行人碰撞现象分析

  在疏散过程中，必定会出现行人碰撞挤压现象。行人间碰撞比例越高,行人受伤概率越大因此研究行人碰撞比例对研究人群整体安全有着重要作用。经典元胞自动机模型无法描述该过程,但在改进模型中可以分析行人碰撞现象。
  图 5-12 是碰撞人数占总疏散人数、剩余疏散人数的比例与时间的关系对比图。结合模拟过程发现,随着部分行人开始汇聚到出口附近，行人间确撞次数明显增加。在离开出口人数小于新到来出口人数时，两出口附诉行人感度不断增大行人间碰撞比例达到整个蔬散过程最高。之后行人间碰撞次数逐渐减少，但因剩余疏散行人数量也在减少, 碰撞人数占剩余疏散人数的比例反而升高，在 45s 达到峰值 30%。若某一时刻，图中两种比例都处于较大值时，意味着行人碰撞总人数较多且待疏散行人碰撞比例较高，此时行人受伤概率最大。综合两条曲线可以发现,在疏散中期(23~40 s),行人间碰撞比例较高,有可能发生安全事故,需要疏散引导。

论文缩略图：

程序代码：

clc
%%%定义参数
renshu=200;
buxingshudu=3;
renweizhi=rand(1,200)*50;
shijian=renweizhi./buxingshudu;
zuidachurenliang=3;
H=0;
for i=1:1:1000
c=find(shijian<=i);
[n,m]=size(c);
H=m;
w(i)=H;
end
for i=1:1:999
meiyimiaodaodarenshu(i)=w(i+1)-w(i);
end
meiyimiaodaodarenshu=[w(1) meiyimiaodaodarenshu]
xianyourenshu=0
for i=1:1:1000
xianyourenshu=xianyourenshu+meiyimiaodaodarenshu(i)
if xianyourenshu>zuidachurenliang
chumenrenshu(i)=zuidachurenliang
else
chumenrenshu(i)=xianyourenshu
end
xianyourenshu=xianyourenshu-chumenrenshu(i) ;
paiduirenshu(i)= xianyourenshu;
end
clc
%%%定义参数
renshu=200;
buxingshudu=3;
renweizhi=rand(1,200)*50;
shijian=renweizhi./buxingshudu;
zuidachurenliang=3;
H=0;
for i=1:1:1000
c=find(shijian<=i);
[n,m]=size(c);
H=m;
w(i)=H;
end
for i=1:1:999
meiyimiaodaodarenshu(i)=w(i+1)-w(i);
end
meiyimiaodaodarenshu=[w(1) meiyimiaodaodarenshu]
xianyourenshu=0
for i=1:1:1000
xianyourenshu=xianyourenshu+meiyimiaodaodarenshu(i)
if xianyourenshu>zuidachurenliang
chumenrenshu(i)=zuidachurenliang
else
chumenrenshu(i)=xianyourenshu
end
xianyourenshu=xianyourenshu-chumenrenshu(i) ;
paiduirenshu(i)= xianyourenshu;
end

clc
%%%定义参数
renshu=200;
buxingshudu=3;
renweizhi=rand(1,200)*50;
shijian=renweizhi./buxingshudu;
zuidachurenliang=3;
H=0;
for i=1:1:1000
c=find(shijian<=i);
[n,m]=size(c);
H=m;
w(i)=H;
end
for i=1:1:999
meiyimiaodaodarenshu(i)=w(i+1)-w(i);
end
meiyimiaodaodarenshu=[w(1) meiyimiaodaodarenshu]
xianyourenshu=0
for i=1:1:1000
xianyourenshu=xianyourenshu+meiyimiaodaodarenshu(i)
if xianyourenshu>zuidachurenliang
chumenrenshu(i)=zuidachurenliang
else
chumenrenshu(i)=xianyourenshu
end
xianyourenshu=xianyourenshu-chumenrenshu(i) ;
paiduirenshu(i)= xianyourenshu;
end

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可