2007年12月 机械科学与技术 December 2007 第26卷第12期 Mechanical Science and Technology for Aerospace Engineering Vol.26 No.12 基于MCMC方法的随机加工时间研究 于晓义，褚崴，孙树栋 (西北工业大学机电学院，西安710072) 于晓义 摘要：在分析随机作业调度问题特点的基础上，建立了随机加工时间统计模型及参数估计模型， 在参数未知及参数已知的条件下，提出了基于马尔可夫链蒙特卡罗方法的随机加工时间统计技术， 并通过吉布斯抽样实现了加工时间的参数估计。通过计算机仿真实验，验证了该方法的可行性及 有效性，为随机作业调度提供更特合实际生产的数据支撑。 关键词：随机作业调度问题；随机加工时间；马尔可夫链蒙特卡罗方法；吉布斯抽样 中图分类号：TP272 文献标识码：A 文章编号：10038728(2007)12-1574-04 Study of Stochastic Processing Time Using MCMC Method Yu Xiaoyi,Chu Wei,Sun Shudong (School of Mechatronics,Northwestem Polytechnical University,Xi'an 710072) Abstract:After analyzing the characteristics of stochastic job-shop scheduling problem,we establish the stochastic processing time statistics model and the parameter estimation model.Under the two conditions that there may exist known or unknown parameters,we propose our stochastic processing time statistics technique based on the Markov chain and the Monte Carlo (MCMC)method.Then we do the Gibbs sampling to estimate the unknown parameters of stochastic processing time.Finally,we carry out the simulation of the technique,and the simulation results veri- fy its feasibility and effectiveness. Key words:stochastic job shop scheduling probelm;stochastic processing time;Markov Chain and Monte Carlo (MCMC)method;Gibbs sampling Job-shop优化调度问题已经成为提高企业生产求。随机加工时间的确定成为随机作业调度的瓶颈 效率的核心问题。工序加工时间是job-shop优化调 问题，因此迫切需要一种对随机加工时间进行准确 度的关键数据。传统的车间作业调度都是依事先确 估计的方法。由于实际可测工序加工时间通常包含 定的加工时间（常值）进行生产作业调度，而在生产 了生产过程中的各种扰动（如机器性能影响因素、 实践中，加工时间往往是随机变化的，传统的方法已 加工班组的影响因素等)，不能直接获得工序的加 不能满足现代企业车间作业调度的需求。考虑随机 工时间。本文提出了一个基于马尔可夫链蒙特卡罗 加工时间的调度问题具有重要意义，随之出现了随 方法的随机加工时间估计方法。 机作业调度的研究1，。然而，在现有的随机作业 调度研究中许多研究者都集中在对随机作业调度算 1问题建模 法的研究上，而对作为随机作业调度的主要数据支1.1随机加工时间问题建模 撵一随机加工时间的研究比较鲜见。上述研究中 以确定工件的某道工序作为描述对象。生产过 的随机加工时间都是通过从假设服从均匀分布的某 程中可测的工序生产周期包含了工件的加工时间、 一范围中取随机数来获得，而且对于其具体的取值 生产过程干扰因素（加工班组的技术因素、操作者 区间亦未作明确的阐述。笔者认为无论随机作业调 的生理因素、设备的不稳定因素等)影响产生的时 度算法效率及性能的高低，如果没有准确的加工时 间偏移。假设： 间数据，其调度方案都不能符合实际的生产调度要 (1)工序加工时间是独立的随机变量 收稿日期：2006-11-13 基金项目：国家863计划项目(2001AA412150,2003A4411110)和教育部博士点基金项目(2004699025)资助 作者简介：于晓义(1980-)，男（汉），山东，博士研究生，uxioi@mil.pu.cd血.cm 万方数据
2007年 第26卷 12月 第12期 机械科学与技术 MeclIanical Science肌d TeclIIIology for AeD∞pace Engin∞ring D∞ember 200r7 V01．26 No．12 于晓义 基于MCMC方法的随机加工时间研究 于晓义，褚崴，孙树栋 (西北工业大学机电学院，西安7100r72) 摘要：在分析随机作业调度问题特点的基础上，建立了随机加工时间统计模型及参数估计模型， 在参数未知及参数已知的条件下，提出了基于马尔可夫链蒙特卡罗方法的随机加工时间统计技术， 并通过吉布斯抽样实现了加工时间的参数估计。通过计算机仿真实验，验证了该方法的可行性及 有效性，为随机作业调度提供更符合实际生产的数据支撑。 关键词：随机作业调度问题；随机加工时间；马尔可夫链蒙特卡罗方法；吉布斯抽样 中图分类号：吼72 文献标识码：A 文章编号：1003名728(2007)12·1574讲 Study of Stochastic Processing Time Using MCMC Method Yu Xiaoyi．Chu Wei。Sun Shudong (Sch00l 0f M∞lIa_嘶IIi∞，Nortllwe8tem Pol”echrIical UIlive玛畸，)【i’肌7100r72) AbSt船ct：触er analy菌ng dle ch眦cteristics of stoch鼬tic job—Shop∞heduling problem，we establish tlle 8tcIch鹊tic pmce鼹ing dme 8tatistic8 mDdel and tlle pammeter estimati∞model．Under出e two c仰ditio璐tllat tlle陀may眈ist lcIlown or unkIlo、^rIl p盯眦eters，we prop0∞ollr stoch鹪tic processing time 8tatistics technique b船ed∞吐Ie Markov chain肌d tlIe Monte Callo(MCMC)metllod．nen we do tlle Gibbs sampling to esti眦te出e吼kIlowIl p蹦珊ete玛 0f 8toch酗tic processing time．Finally，we c邺7 out tlle siInulation of山e tecllIlique，and tlle 8imulation托叭lts Veri一 匆its fe黯ibility明d ef＆ctiveness． 15【ey words：stochastic job shop scheduling pIobehn；stochastic pI．ocessing time；MarkoV Chain肌d MoIIte Carl0 (MCMC)method；Gibbs s锄pling Job一8hop优化调度问题已经成为提高企业生产 效率的核心问题。工序加工时间是job-shop优化调 度的关键数据。传统的车间作业调度都是依事先确 定的加工时间(常值)进行生产作业调度，而在生产 实践中，加工时间往往是随机变化的，传统的方法已 不能满足现代企业车间作业调度的需求二考虑随机 加工时间的调度问题具有重要意义，随之出现了随 机作业调度的研究H．9】。然而，在现有的随机作业 调度研究中许多研究者都集中在对随机作业调度算 法的研究上，而对作为随机作业调度的主要数据支 撑——随机加工时间的研究比较鲜见。上述研究中 的随机加工时间都是通过从假设服从均匀分布的某 一范围中取随机数来获得，而且对于其具体的取值 区间亦未作明确的阐述。笔者认为无论随机作业调 度算法效率及性能的高低，如果没有准确的加工时 间数据，其调度方案都不能符合实际的生产调度要 求。随机加工时间的确定成为随机作业调度的瓶颈 问题，因此迫切需要一种对随机加工时间进行准确 估计的方法。由于实际可测工序加工时间通常包含 了生产过程中的各种扰动(如机器性能影响因素、 加工班组的影响因素等)，不能直接获得工序的加 工时间。本文提出了一个基于马尔可夫链蒙特卡罗 方法的随机加工时间估计方法。 l问题建模 1．1随机加工时间问题建模 以确定工件的某道工序作为描述对象。生产过 程中可测的工序生产周期包含了工件的加工时间、 生产过程干扰因素(加工班组的技术因素、操作者 的生理因素、设备的不稳定因素等)影响产生的时 间偏移。假设： (1)工序加工时间是独立的随机变量 收稿日期：2006—1l—13 基金项目：国家863计划项目(2001AA412150，2003从411110)和教育部博士点基金项目(2004699025)资助 作者简介：于晓义(1980一)，男(汉)，山东，博士研究生，”xiaoyi@mail．nWpu．edu．∞ 万方数据
第12期 于晓义等：基于MCMC方法的随机加工时间研究 1575 (2)生产过程中干扰因素对加工时间的影响是 次的工序生产周期样本数据，因此要求i≥2。 独立的随机变量。工序生产周期可用随机效应模 在初始阶段采用一个从三阶段分层模型简化 型2表述，即 而来的两阶段模型处理第1个加工批次的样本数 Xij=0:+Dij 据。将初始的工序加工时间O,作为一个未知参数而 0:N(u,o2) (1) 不作为随机变量考虑。第一阶段用于描述第1个加 D:N(0,c) 工批次的数据分布，第二阶段描述预先未知的参数 式中：表示工件的批次，表示批次中工件的序号； O,和o2的先验分布。两阶段分层模型(two-stage X,表示被测量工件的工序生产周期；O表示工序加 hierarchical model,.2HM)如下所示。 工时间，设其服从正态分布O:N(4,o);D表示加 第一阶段：Xy:N(O1,c) (5) 工时过程干扰产生的时间偏移，设其服从正态分布 第二阶段：O1:N(4,o2);o:IC(a1,b,) D:N(0,o)。 通常情况下参数a1,b,具有先验信息：a1= 由于过程干扰影响因素的存在使得O:D,不能通 b,=0.001。为了估计未知参数om,至少需已知该批 过直接测量得到，因此问题的关键是估计未知变量0：。 次的2个工序加工周期的实测数据，因此要求j≥2。 1.2参数估计问题建模 工序加工时间估计表达式为 在参数估计方面，运用贝叶斯理论1不仅可以 01= (6) 充分利用模型信息和样本数据信息，而且也融合了 0y=E(01x",2HM)j=2,,(7) 模型总体分布中未知参数的信息。它克服了传统的 通过以上所述可以看出单向分类随机效应模 静态模型难以处理突发事件的缺陷，具有灵活、易于 型从理论上很好地解决了参数估计问题，但是要想 适应外部变化的特点，它可以很好地解决传统统计 得到用于进行参数估计的边缘后验分布p(O:I 方法的样本不足及样本质量不佳等问题，更适合进 x),还需要求解复杂的高维积分)，其难度较大。 行模型的异常预测，反映现实问题。 由式(1)可生成一个单向分类随机效应模 2随机加工时间确定方法 型，以贝叶斯理论的观点，这个单向分类随机效 求解复杂高维函数的积分能力对贝叶斯统计显 应模型是分层模型的一个特例，即3阶段分层模型 得尤为重要。然而直接求解这些积分往往是不可能 (Three-stage Hierarchical Model,3HM) 的。为了解决计算边缘后验分布的复杂高维积分问 第一阶段：X:N(0,o) 题，本文采用Markov Chain Monte Carlo(MCMC)方 第二阶段：0：N(μ，o2);o:1C(a1,b1) (2) 法，它的基本思想是模拟一个马尔可夫链，此链的静 第三阶段：：N(,o);o2:lG(az,b2) 止分布即为要得到的后验分布（目标分布），当运行 式中：IG表示逆Gamma分布；通常情况下参数o, 这个链足够长时，从这个链得到的模拟值可视为来 o,a1,b1,a2,b2具有先验信息：41=b1=a2=b2 自目标分布的独立样本，接着可以用得到的这些模 =0.001;o=0;o后=1×100+10；X(,x2, 拟样本来总结目标分布的重要性质。这样MCMC …,x)为测量的工序生产周期样本值。 隐含地求解了积分，给出了后验分布的模拟结果。 通过计算后验分布(O:I,3HM)的数学期 MCMC方法中吉布斯(Gibbs)抽样2)是最为 望值可得到工序加工时间的贝叶斯估计O,即 行之有效的一种抽样方法，Gibbs抽样是最简单也 0,=E(0,r,3HM) 是应用最广泛的一种抽样方法，对于分层模型而言， (3) 其优势在于可以不考虑参数及阶段数量的影响，根 提到的贝叶斯分层模型的一大优势是可以从采 据以上介绍的参数模型以及观测所得的实测数据， 集的以往工序加工周期的样本信息中导出未来加工批 通过Gibbs抽样器估计分层模型中未知参数的后验 次的工序加工时间的后验预测分布，即(O:,|,3 分布，并利用参数估计得到工序加工时间O。 HM)(J表示第i个加工批次的批量)，可预测第i+1个 设Z=(21,Z2,…,23)=（01,02,…,04; 批次的加工时间的贝叶斯估计，即 σ2，σ)，其中，Z表示分层模型中的未知参数。由于 6=E(OX4,3 HM) (4) 分层模型中未知参数的分布具有共轭性)，从而可 获知相应的未知参数的全条件后验分布p(Z。Ix, 其中由于三阶段分层模型中存在未知参数σ2， Z1,…,2-1,241,…,23）(k=1,2,…,)。对于三 在求解0，以及预测ò1时至少需要2个不同加工批 阶段模型，其参数的先验联合分布为 万方数据
第12期 于晓义等：基于McMC方法的随机加工时间研究 1575 (2)生产过程中干扰因素对加工时间的影响是 独立的随机变量。工序生产周期可用随机效应模 型幢’表述，即 lxti 2 o‘+Dij {仇：Ⅳ(p，矿2) (1) o Di，：Ⅳ(o，盯：：I) 式中：f表示工件的批次，_『表示批次中工件的序号； 置，表示被测量工件的工序生产周期；D；表示工序加 工时间，设其服从正态分布D；：Ⅳ(p，盯2)；现f表示加 工时过程干扰产生的时间偏移，设其服从正态分布 D‘，：』＼，(0，盯：)。 由于过程干扰影响因素的存在使得0i以j不能通 过直接测量得到，因此问题的关键是岱汁未知变量q。 1．2参数估计问题建模 在参数估计方面，运用贝叶斯理论”’不仅可以 充分利用模型信息和样本数据信息，而且也融合了 模型总体分布中未知参数的信息。它克服了传统的 静态模型难以处理突发事件的缺陷，具有灵活、易于 适应外部变化的特点，它可以很好地解决传统统计 方法的样本不足及样本质量不佳等问题，更适合进 行模型的异常预测，反映现实问题。 由式(1)可生成一个单向分类随机效应模 型【4】，以贝叶斯理论的观点，这个单向分类随机效 应模型是分层模型的一个特例，即3阶段分层模型 (Three-stage Hiemrchical Model，3HM)o 第一阶段：五f：Ⅳ(D。，盯：) 。第二阶段：o‘：Ⅳ(肛，盯2)；盯：：佑(口l，61) (2) 第三阶段：p：Ⅳ(脚，蠢)；矿：佑(口2，62) 式中：，G表示逆Gamma分布；通常情况下参数m， 矿：，口1，6l，口2，62具有先验信息‘51：口1=61=口2=62 =o．ool；m=o；盯：=l×10W+lo；r’(菇1l，髫12， …，毛，)为测量的工序生产周期样本值。 通过计算后验分布(D；I刀J，3 HM)的数学期 望值可得到工序加工时间的贝叶斯估计伉。即 ． I Di=E(0i l刀’，3HM) (3) 提到的贝叶斯分层模型的一大优势是可以从采 集的以往工序加工周期的样本信息中导出未来加工批 次的工序加工时间的后验预测分布，即(D：+。l Ft，3 HM)(^表示第z个加工批次的批量)，可预测第i+1个 批次的加工时间的贝叶斯估计玑，即 ． oi+l=E(D：+l l∥‘，3 HM) (4) 其中由于三阶段分层模型中存在未知参数cr2， 在求解或以及预测。州时至少需要2个不同加工批 次的工序生产周期样本数据，因此要求f≥2。 在初始阶段采用一个从三阶段分层模型简化 而来的两阶段模型处理第1个加工批次的样本数 据。将初始的工序加工时间D。作为一个未知参数而 不作为随机变量考虑。第一阶段用于描述第1个加 工批次的数据分布，第二阶段描述预先未知的参数 D。和盯：的先验分布。两阶段分层模型(two—stage hierarchical model，2HM)如下所示。 第一阶段：墨，：Ⅳ(D。，盯：) (5) 第二阶段：D1：Ⅳ(肛，盯2)；仃：：，G(口1，61) 通常情况下参数口。，6。具有先验信息H】：口，= 6。=o．00l。为了估计未知参数蠢，至少需已知该批 次的2个工序加工周期的实测数据，因此要求_『≥2。 工序加工时间估计表达式为 6。。=菇。。 (6) 6v=E(Dl l茹V，2HM)．，=2，…，．，l (7) 通过以上所述可以看出单向分类随机效应模 型从理论上很好地解决了参数估计问题，但是要想 得到用于进行参数估计的边缘后验分布p(D；I ∥’)，还需要求解复杂的高维积分‘111，其难度较大。 2随机加工时间确定方法 求解复杂高维函数的积分能力对贝叶斯统计显 得尤为重要。然而直接求解这些积分往往是不可能 的。为了解决计算边缘后验分布的复杂高维积分问 题，本文采用Markov Chain Monte Carlo(MCMC)方 法，它的基本思想是模拟一个马尔可夫链，此链的静 止分布即为要得到的后验分布(目标分布)，当运行 这个链足够长时，从这个链得到的模拟值可视为来 自目标分布的独立样本，接着可以用得到的这些模 拟样本来总结目标分布的重要性质。这样MCMC 隐含地求解了积分，给出了后验分布的模拟结果。 MCMC方法中吉布斯(Gibbs)抽样[2’51是最为 行之有效的一种抽样方法，Gibb8抽样是最简单也 是应用最广泛的一种抽样方法，对于分层模型而言， 其优势在于可以不考虑参数及阶段数量的影响，根 据以上介绍的参数模型以及观测所得的实测数据， 通过Gibbs抽样器估计分层模型中未知参数的后验 分布，并利用参数估计得到工序加工时间D。。 设Z=(Zl，z2，…，zi+3)=(Dl，D2，…，0‘，肛； 矿2，矿：)，其中，z表示分层模型中的未知参数。由于 分层模型中未知参数的分布具有共轭性口】，从而可 获知相应的未知参数的全条件后验分布p(五I茗“， z。，…，乏一l，乙+l’．一，互+3)(j}=l，2，…，i)。对于三 阶段模型，其参数的先验联合分布为 万方数据
1576 机械科学与技术 第26卷 p(x,0,o4,o2) -3)(3)exp( 0+-0)+6 (8) 再根据文献[6,8]中介绍的方法可推出各参数 2鼎==p(d1,09,0,…,0,μ9o0)(11) 的全条件分布为 即Z0={09,0,…,09,0,c20,c}, p(μlx,01,…,0,o2,)c 并最终得到21)，Z2),Z3),…,Z9。计算由不同 Mg4+店0gd 的2)出发，马尔可夫链经过一段时间的迭代后，可 以认为各时刻Z的分布为平稳分布，此时它收敛。 而在收敛出现前的m次迭代中，各状态的分布还受 p(l,0,…,0,h,o)c1c(a+2,4 到初始值的影响，不能认为是平稳分布，因此应将前 +容0，-川） m个迭代值去掉，即 10aac+宫+小. 乙亡名24=123 (12) 由此可知0(k=1,2,…,)的边缘后验分布。 +含0r+含-0r川 3收敛性判断与计算机实验 p(0lx,044,o2,o)c 本文MCMC的收敛性判断方法是采用Gibbs抽 N jo' 样同时产生多个马尔可夫链，在经过一段时间后，如 +02+o40+ (9) 果这几条链稳定下来了，则其收敛。该方法需同时模 式中：0表示除0外的0：参数；元表示第k批次工 拟m条链，每条链迭代2n次，设诊断指标为R,其构 件观测数据的样本均值。两阶段模型的全条件分布 造过程如下 的计算方式与三阶段模型类似，此处不再作说明给 计算A,用以估计m条链均值间的变异 定任意初始值2={2,2，…，2,29， Z9,z9}={0o,0…,0ho,a20,c201, 4=2a-0 (13) Gibbs抽样的第一次迭代如下 2=0=p(0,l",0,0,…,0,h,o20,0) 武中0200：含a m 29=0=p(0l,0,0,,0%,μ，00,0) 计算B,它反映m条链链内变异的平均水平，用 以估计 Z0=09=p(0,1,0,0,,08,μ，09，o0) B= (14) Z=u9=pulx,00,09,…,00,d20,oo) m 2=d2n=p(d21x,00,0,…,0,,20) 式种财点玄g- 2得=m=p(01,0,0,…,00,μ"，o2) (10) 计算C=”。8+1+A,最后计算诊断 m 以上完成了由z0)={00,0,…,0, 指标R,即 u0,020,0201到2={00,00，…，0， R=C/B (15) 4,02,}的转移。经过t次选代，可得 判断标准：如果R>1,说明链不收敛，R≈1， Z9=09=p01,G0,g9,,0,h,o,2) 提示链达到了静止状态。由于~R的计算利用的是迭 29=0=p(01,0P,g…0,μ90) 代结果的后50%的数据，因此当判断为收敛时，说 明前n次可作为预迭代，再迭代n次即可用来估计。 2=0°=p(0,1,0,0,…,0唱，μ-ng2-n) Brooks和Gelman]于1998年对上述方法提出了基 2Z9=40=pul,09,0,…,09,o2-w,o-) 于该思想的另一条简便且易于实现的方法。对于每 Z=o20=p(d21,0,0,…,0,4,o-) 一条链求出区间100(1-)%的长度，这样对于m 条链可求得m个区间长度，其平均值记为l。然后， 万方数据
1576 机械科学与技术 第26卷 p(菇“，DI，盯：，弘，盯2)oc 矿撕％?一，ex小号(丝竿与塑+半+半))(8) 再根据文献[6，8]中介绍的方法可推出各参数 的全条件分布为 p(p l髫订，Dl，…，oi，盯2，盯：)oc Ⅳf，，砌+蠢荟仉仃：盯：] L ，+￡盯： ’cr2+i蠢／ p(以l髫“，01，…，oi，p，盯2)oc佑(口1+々，61 ． i +寺∑(仇一p)2) p(，l菇“，Dl，…，oi，p，以)oc，G(a2+吉(苫^+五)， 磁=群’=p(蠢I∥，，0f”，谚’，…，睇’，∥。矿’)(11) 即z‘‘)={D：”，D{”，…，D：”，p“’，盯2“’，盯警”}， 并最终得到z‘¨，z‘扪，z‘”，…，z‘”。计算由不同 的z∞’出发，马尔可夫链经过一段时间的迭代后，可 以认为各时刻Zo’的分布为平稳分布，此时它收敛。 而在收敛出现前的m次迭代中，各状态的分布还受 到初始值的影响，不能认为是平稳分布，因此应将前 m个迭代值去掉，即 磊’点薯+。z：。^=1，2，…，Ⅲ(12) 由此可知D。(后=1，2，…，f)的边缘后验分布。 ．1—1 “ ， 6：+丢蚤蚤(％一仇)2+耋(％一D；)2) o I=l p2l p2l ， 3收敛性判断与计算机实验 …一…一………一 p(DI I茗”，D—I，p，盯2，盯：)％ Ⅳ(考蛩·+荸匆吾毛)c9， 式中：0．I表示除D。外的D；参数；瓦表示第后批次工 件观测数据的样本均值。两阶段模型的全条件分布 的计算方式与三阶段模型类似，此处不再作说明。给 定任意初始值z‘0’={z：∞，乏们，…，z：∞，z船， z：2，z勰}={D：们，Dio’…，D：们，肛‘们，盯2‘∞，盯：o’}， Gibbs抽样的第一次迭代如下 z：1’=D：1’=p(DlI∥’，掣，0；0’，…，唾∞，p彻，一，硭∞) z!”=噬1’=p(D2I菇玎，D：n，硝∞，…，噬∞，p∽，∥∞，旌∞) 乏"=唾"=p(0lI∥’，D：”，Di¨，…，D盥，·p油，一∞，醒∞) Z2=肛o’=p0I∥7，D：”，D：¨，…，D：”，，‘0’，醒o’) 础=一"=p(矿I∥7，D：1)'Di”，…，D：¨，p‘1’，醒∞) 墨2=碟1’=p(蠢I∥’，硝n，噬”，…，研”，肛o’，，“’) (10) 以上完成了由z∞’={D：∞，D5∞，…，D：们， 弘‘∞，cr2‘们，叮警o’}到z‘1’ = {D：¨，D：¨，…，o：¨， 弘‘¨，，“’，盯：1’}的转移。经过t次迭代，可得 z}’=研’=p(ql∥’硝哪硝哪’．．．群由，∥…，萨埘，毋由) 霹’=谚’=p(q I∥7群’，礴哪，．··，醇由，∥‘却矿哪露哪) 零’=醇’=p(Di l∥7，钟’，谚’，…，皑，∥“’矿q’群q’) Z：：=p“’=p缸I∥，D：”，Di”，…，D：”，cr2‘“”，醒‘州) z：：：=∥‘’=p(，I∥’，DhDiI’，…，研”，p‘”，旌㈣) 本文McMC的收敛性判断方法是采用Gibbs抽 样同时产生多个马尔可夫链，在经过一段时间后，如 果这几条链稳定下来了，则其收敛。该方法需同时模 拟m条链，每条链迭代2，1次，设诊断指标为尺，其构 造过程如下 计算A，用以估计m条链均值间的变异 A=击蚤(瓦一石)2 (13) ．2^ ． m 式中：瓦。寺萎。D：；石=去善玩 o’I。n+l ¨o i=l 计算B，它反映巩条链链内变异的平均水平，用 以估计盯： (14) 计算c=生三一b+f ¨ 、 l+上1A，最后计算诊断 m， 指标R，即 R=C／B (15) 判断标准：如果R>1，说明链不收敛，尺一l， 提示链达到了静止状态。由于馏的计算利用的是迭 代结果的后50％的数据，因此当判断为收敛时，说 明前n次可作为预迭代，再迭代厅次即可用来估计。 Brooks和Gelm肌H1于1998年对上述方法提出了基 于该思想的另一条简便且易于实现的方法。对于每 一条链求出区间100(1一a)％的长度，这样对于m 条链可求得m个区间长度，其平均值记为Z。然后， 饿 《 0● 。∑㈨h∑一 上m 专^了 B = 式中 ≯ 万方数据
第12期 于晓义等：基于MCMC方法的随机加工时间研究 1577 用整个m条链的模拟值求出区间100(1-α)%的 表1MCMC收敛性判断的样本数据 长度，记为L,此时诊断指标为R=L/1,其诊断标准 批次 工序加工周期样本值(min) 与原始方法相同，下面以三阶段模型为例，应用该方 149 139 139 147 153 149 150 145 151 144 法说明MCMC的收敛性判别 154 150 155 146 151 取6个生产批次的样本数据，如表1所示，其数 139 139 154 141 149 154 158 146 158 157 据来源为某企业零件0EU73979的第15道精车工序 6 147 140 140 143 139 的实测加工周期。设α=0.2，取两组初始数据Z, mu chains 1:2 sigma2.btw chains 1:2 =10,0,00,00,0g,0g4,c20,d201= 5 15 10,0,0,0,0,0,0,0,0}及z={1,1,1,1,1,1,1, 05 0.5 S001000 500i000 1,1}用winBUGS软件)模拟生成两条链，其收敛 (a)选代次数 b)运代次数 性判断如图1（图中mu表示4，sigma2.btw表示c2, sigma2.with chains 1:2 theata[6]chains 1:2 sigma2.with表示o,theata[i]表示O)所示。从图 05 1中可以看出，在900次迭代（大概需要0.5s)后有 5001000 5001000 关9个参数的R值趋近1，说明链条收敛，相应参数 (c)运代次数 (d)迭代次数 估计值如图2所示。 图1MCMC收敛性判断图 送Node node incan d MC cmor 2.5 975% start sumple 147.5 2233 0.0094211431 1475 1520 S000 190002 sigma2.btw 23.85 4351 02224 000823 1447 1058 50 192 sigmna2.with 28.5 10.06 0.09245 15.06 26.44 53.68 5000 190002 theataf I 146.1 2.021 0.01232 142.0 1462 149.9 5000 190002 theatal2] 1477 1.935 009035143.9 147.7 1516 5000 1 theata[3] 149.9 22 001904 145.8 149.9 154.3 5000 190002 theatal4] 145.5 2.125 0.01586 1412 145.5 1494 500 1900m2 theata[5] 152.1 2782 003319 1465 1523 157.1 S000 190X002 theata6] 1438 2.506 002634139.1 143.7 1486 5000 190002 图2参数估计值 4结束语 in0p.Res.,2003,10:577-596 工序随机加工周期统计分析的核心任务是确定 [2]Gelfand A E,Smith AFM.Sampling-based approeches to caleu- 工序的随机加工时间与生产过程干扰时间偏差分布 lting marginal densities[J].Journal of the American Statisti- cal Ass0 ciation,1990.85:398-409 中未知参数的估计值及其分布密度函数。由于加工 [3】茆诗松.贝叶斯统计[M].北京：中国统计出版社，1999 过程中随机干扰的存在，使得实际的工序加工时间 [4]Mostafa S M,Ahmad R.Empirical Bayes quadratic estimators of 无法通过直接测量得到，因此工序加工时间参数估 variance components in normal linear models[J].Statistics, 计及其精度很大程度上决定了工序加工周期统计分 1986,17:337-348 析性能的优劣。考虑到MCMC方法在参数估计方 [5]Spiegelhalter D J,Thomas A,Best N G and Gilks W R.BUGS: Bayesian Inference Using Gibbs Sampling version 0.30) 面的优势且具有很强的预测能力，本文选择了MC [R].Cambridge:Medical Research Council Biostatistics MC方法，通过实例研究证明了该方法适用于解决 Uni,1994 车间零件工序随机加工时间统计问题。该问题的解 [6]Trietsch D.The harmonie rule for proceas setup adjustment with 决使得依据该方法确定的加工时间分布函数制定的 quadratic los[].Journal of Quality Technology,1998,30 车间作业调度计划将更加符合车间生产的实际情 (1),75-84 况，在随机车间作业调度领域有着很好的应用前景。 [7]Brooks S P,Gelman A.Altemative methods for monitoring conver- gence of iterative simulation[.Journal of Computational and Graphical Statistics,1998,7:434-455 [参考文献] [8]Bex G E P,Luceno A.Statistical Control by Monitoring and Feedback Adjustment[M].New York:Wiley,1997 [1]Yoehitomi Y,Yamaguchi R.A genetic agorithm and the Monte [9]Yasunari Y.A genetic algorithm approach to solving sochastic job- Carlo method for stochastic job-shop scheduling[J].Intl.Tran. shop scheduling probems[J).Intl.Tran.in Op.Res.,2002,9 万方数据
第12期 于晓义等：基于McMc方法的随机加工时间研究 1577 用整个m条链的模拟值求出区间100(1一a)％的 长度，记为￡，此时诊断指标为R=∥Z，其诊断标准 与原始方法相同，下面以三阶段模型为例，应用该方 法说明McMC的收敛性判别 取6个生产批次的样本数据，如表1所示，其数 据来源为某企业零件DEU73979的第15道精车工序 的实测加工周期。设口=o．2，取两组初始数据z{∞ ={Df0)，噬们，硝m，D：们，D；m，o?’租‘0’，cr2‘∞，醒o’}= {o，o，o，o，o，o，o，o，o}及z{o)={l，1，l，1，1，1，11， 1，l}用winBuGs软件瞪1模拟生成两条链，其收敛 性判断如图l(图中mu表示肛，sigma2．b晰表示cr2， sigma2．with表示盯：，theata[i]表示Df)所示。从图 1中可以看出，在900次迭代(大概需要O．5 s)后有 关9个参数的R值趋近l，说明链条收敛，相应参数 估计值如图2所示。 表l MCMc收敛性判断的样本数据 批次 工序加工周期样本值(IIlin) 1 149 139 139 147 153 2 149 150 145 151 144 3 154 150 155 146 151 4 139 139 154 141 149 5 154 158 146 158 157 6 147 140 14JD 143 139 1．5 髦 l 0．5 O 1．5 宅 l 0．5 0 (a)迭代次数 图2参数估计值 4结束语 工序随机加工周期统计分析的核心任务是确定 工序的随机加工时间与生产过程干扰时间偏差分布 中未知参数的估计值及其分布密度函数。由于加工 过程中随机干扰的存在，使得实际的工序加工时间 无法通过直接测量得到，因此工序加工时间参数估 计及其精度很大程度上决定了工序加工周期统计分 析性能的优劣。考虑到MCMC方法在参数估计方 面的优势且具有很强的预测能力，本文选择了MC． MC方法，通过实例研究证明了该方法适用于解决 车问零件工序随机加工时间统计问题。该问题的解 决使得依据该方法确定的加工时间分布函数制定的 车间作业调度计划将更加符合车间生产的实际情 况，在随机车间作业调度领域有着很好的应用前景。 [参考文献】 [1】Y∞hitorni Y。Y锄agIlcK R．A genetic 8l鲥tIlm衄d tl他Monte cdo mcthod细8toclIas6c job-sh叩幽edtlIing[J]．mⅡ．伽n． 1 500 lOoo ‘c)迭代次数 1．5 笔 l O．5 O l 500 looO (b)迭代次数 图1 McMc收敛性判断图 in 0p．Il馏．。2003。10：577—596 [2] Gd钿d A E，slnitlI A F M．SalnpliIlg—bascd叩pr㈣IIc8 to cal叶 l秭llg删I晒rlal d∞8id∞[J]．Jo帆址of血e Amed锄st娟sti- ∞l As8adan∞，1990，85：398—4∞ [3]茆诗松．贝叶斯统计[M]．北京：中国统计出版社，1999 [4】 Mosta缸S M。^hmBd R．Empirical Bay∞qImdIalic e硪ma姗蠢 vari蚰∞咖p哪elIls in咖哪Ial Une盯moddl8[J]．戳atisH岱， 1986。17：337—348 [5】spie嘲IIalt盯D J，n锄阳A。B鳅N G锄dcnl【8w凡双『(冷： Bayesi粕五疵咖∞Using Gibh s锄pIi豫(ve暗蛆0．30) [R]．cam晡d护： Medic8l m黔毗h c咖cd Biostali硝∞ Unit．1994 [6] 啊dt8ch D．11Ie 118珊ollic mle细p眦∞B瞅叩adj哪hnent诵tII quadl“c l∞B[J]．Jo㈣l of Qu址ty TechoIo留，1998，30 (1)。75—84 [7】 B-∞k8 sP。Gd咖A．Alt咖6vc枷od8fh咖dtoling∞nv小 gen∞蠢iter娟ve 8i删lalio鹏[J]．J伽咖叫of o珊I删伽叫 明d G哺phicm S协listi娼，1998。7：434—455 [8]B0xC E P，hl咖oA．StatistlcaIc咖砌byM呵t岫呜and F∞dbacl‘Adjm衄蜘t[M]．N删Y0rk：Wiley。199r7 [9】 Y聃uIl8ri Y．A舻n棚c a190ritIIm ap呻ch to 50lvillg sI。ch娟c job· 8hop 8ched“Ilg pfobe咖[J]．hIn．伽n．m op．R嚣．，2002，9 万方数据