张惠煜,陈庆新,毛 宁
广东工业大学 广东省计算机集成制造重点实验室,广州 510006
定制型装备制造企业生产系统属于复杂的离散制造系统,具有随机不确定性(订单到达、工序工时以及返修等不确定因素)、生产周期长、拖期现象严重以及多品种大规模在制品等特点。由于生产车间物料搬运任务十分庞大,基于AGV(Automatic Guided Vehicle)的物料搬运系统(Material Handling System)成为生产车间至关重要的组成部分。AGV的数量对系统性能的优化具有显著的影响,然而这类设备的价格昂贵,在制造系统中如何以最低的成本优化配置AGV 的类型及数量,以满足系统产能及订单交货期的需求,是系统设计阶段的关键问题[1-2]。本文考虑的AGV配置问题是一个具有双重约束的随机非线性整数规划问题,然而这些约束(系统产能和订单交货期)无法用决策变量(AGV数量)的封闭形式表达。因此,建立一个有效的性能分析模型以获得系统运行指标是求解AGV配置优化问题的基础。
AGV 的数量配置问题是AGV 系统设计及控制的主要问题之一,目前求解的主要方法可以分为两类:解析方法(基于数学理论建立数学模型进行求解计算)和仿真方法(基于仿真平台建立仿真模型进行实验计算)。由于所涉及问题的复杂性,解析方法只能获得问题的近似解,基于经验统计的仿真实验可以应用于复杂系统获得较为精确的结果,然而仿真实验需要消耗大量的时间[3-4]。文献[5]根据AGV行程需求是否已知将问题划分为确定需求模型(行程需求已知)和随机需求模型(需要估计AGV 的空载行程),并综述不同的求解方法。对于随机系统主要应用排队理论进行求解,文献[6]建立排队模型计算AGV 的数量使AGV 的可用概率满足工作站的需求,但这个模型在求解稀疏车流系统时往往高估了AGV的需求数量。文献[7]根据G/G/c排队模型计算AGV 在传输和请求之间的等待时间,以等待时间为约束设计算法优化AGV 数量,该模型不考虑批量处理的情况。文献[8]基于闭排队网分析系统稳态行为,并建立线性规划模型估算AGV数量,然而闭排队网模型是在系统总在制品数不变基础上的。
本文研究的制造系统具有两个特征:(1)随机批量搬运,成批的物料运输是一个批处理的过程,为减少AGV空车等待时间缩短生产周期,在AGV到达时若不能满载也立即运往目的地,因此运输的批量具有随机性。(2)同等并行AGV,相同类型的AGV具有相同的性能参数(如额定速率、额定容量等),因此多台同类型的AGV是同等的,且并行地在各个车间之间执行物料运输任务。上述特征使得系统模型更为复杂,求解更为困难。
综上所述,目前优化AGV 配置的规划模型并未考虑随机批量处理的问题,而针对具有随机批量处理特征且缓冲区容量有限的排队网性能分析方法也尚未见报道。因此,本文针对定制型生产方式,考虑具有随机批量搬运和同等并行AGV的制造系统,基于Markov过程的排队网络理论建立系统的性能分析模型,进而求解以最小化AGV 投资成本为目标,具有系统产出率和生产周期双重约束的AGV配置优化模型。
如图1所示,工件需要先后经过粗加工车间和精加工车间进行加工,AGV 在两个车间的半成品缓存区之间并行地往返运输工件。AGV具有四种状态:装载(含“饥饿等待”)→负载行程→卸载(含“阻塞等待”)→空载返程。(1)装载时的饥饿等待:AGV 到达半成品缓冲区Buffer2 装载工件时,若 Buffer2 不为空,则 AGV 即使不能满载也立即运往半成品缓存区Buffer3;若Buffer2 为空,则AGV 需要等待直至有1 个工件在粗加工车间完工并进入Buffer2,这种状态称为“饥饿”(Starvation)。(2)卸载时的阻塞等待:AGV到达半成品缓存区Buffer3卸载工件时,若Buffer3 已满而AGV 尚不能完全卸载,则AGV 需要等待直至全部工件卸载完成才能空车返回,这种状态称为“阻塞”(Blocking)。由于阻塞、饥饿等因素的影响降低了AGV的效率。
图1 AGV物料搬运的制造系统
本文考虑以系统产出率和生产周期为约束的制造系统中AGV 配置的优化问题,其目标是优化AGV 数量,从而最小化AGV投资成本,同时满足系统产能及订单交货期的要求。
定义参数如下:
i:AGV类型,i=1,2,…,n;
ui:第i类AGV的采购单价;
xi:第i类AGV配置的数量;
X:AGV数量配置向量,X={xi};
Q:AGV投资成本;
Θ:系统平均产出率,即单位时间系统输出完成品的平均数量;
Γ:系统平均生产周期,即每个工件在系统中停留的平均时间;
E{·}:随机函数的数学期望;
ξ:随机元。
该优化问题的数学模型如下:
式(1)表示优化目标,即求解AGV优化配置的向量X*,使得AGV 投资成本最小;式(2)表示平均产能约束,实际平均产出率不小于目标平均产出率Θmin;式(3)表示平均生产周期约束,实际平均生产周期不大于目标平均生产周期Γmax;式(4)表示X为非负整数型向量。由于E{Θ(X;ξ)}和E{Γ(X;ξ)}无法用明显的数学函数式进行描述,无法应用传统的非线性整数规划方法求解。为此,针对问题的随机性和非线性,本文提出一种基于排队网求解系统性能指标的算法对该优化模型进行求解。
目前针对物料搬运系统性能分析的模型已有大量的研究成果,采用的建模方法包括整数规划模型、网络流模型、排队论模型和Markov模型等,但这些建模方法各有不足[9]。文献[10]根据文献[9]提出的扩展马尔可夫链模型建立考虑多运输AGV堵塞特性系统的排队网络模型,但其建立的均为闭环排队网模型(即系统在制品数量为常数,系统无实时输入流、输出流),而且物料运输过程中不考虑批量处理的问题。文献[11]针对物料搬运系统建立有限缓冲区开排队网模型,但未考虑批量搬运情形。针对具有有限缓冲区、随机批量搬运的问题,文献[12]考虑服务速率与批量大小相关的状态相关排队模型并提出一种Unifying Method 分析其性能指标,但均为针对单个队列的模型,并不能直接应用于具有复杂拓扑结构、队列之间存在耦合关联的网络系统。
人们所面临的系统在制品数量为随机变量且各缓冲区容量是有限的,因此采用有限缓冲区的开排队网模型描述。有限缓冲区开排队网模型在制造系统建模与分析中具有重要应用,与无限缓冲区的模型不同,堵塞和死锁(Deadlock)的现象使得其不具有乘积形式的解(Product-Form Solution)[13],因此求解难度更大。传统的精确计算方法[13]虽然可以获得精确解,但随着系统规模的增大(如系统队列数、缓存容量等)会导致状态空间的“维数灾”问题。为解决此问题,文献[14]提出一种基于Markov过程的“状态空间分解法”建立系统的性能分析模型,且仿真结果验证了该方法近似结果的精确性和有效性。文献[15]应用该分解方法以任务拒绝率为约束条件提出一种缓冲区容量配置的优化方法。
针对传统制造系统的状态空间分解法将整个系统按照节点(缓存区+工作中心)分解为若干个子系统,每个子系统是相对独立的,并通过建立连续时间马尔可夫链(Continuous-Time Markov Chain,CTMC)[16]模型建立子系统之间的耦合关联。在每个子系统中,按CTMC模型分别建立描述子系统的状态空间,并根据“生灭过程”理论分析状态转移规律,构建一系列的稳态平衡方程组,接着迭代求解获得每种状态的稳态概率,最后依据这些状态概率值计算系统的性能指标值。与传统的精确法对比,这在很大程度上降低了系统状态空间的规模,因此该方法能够求解较大规模的系统模型。基于传统状态空间分解法的思路,本文将其拓展,应用于含随机批处理、同等并行AGV的制造系统建立性能分析模型。
建立如图2所示的有限缓冲区开排队网模型,其中方框表示有限容量的缓存区,圆圈表示工作中心,弧线表示工件流,箭头表示工件流向。
图2 制造系统排队网示意模型
该模型满足以下假设条件:
工件到达系统的随机过程为泊松过程,每个工件的到达是独立的,平均到达速率为λ。
系统排队模型为“混合制(损失、阻塞)”[17]:工件到达系统时若缓存区B1已满,则工件被拒绝;否则进入B1,若工作站W1正在加工,则其加入队列排队。
每个工作站包含一台加工设备,服务时间(含设备准备时间)服从负指数分布,W1和W4的平均加工速率分别为μ1和μ4,每个工件的加工时间是独立的,且加工过程一旦开始就不可中断。
缓冲区B1、B2和B4的容量是有限的,最大容量分别为N1、N2和N4。
系统服务规则为“先到先服务(FCFS)”。
系统阻塞机制为“服务后阻塞(BAS)”。
不同类型的AGV 具有不同的行驶速率,且到达B2或B4的间隔时间是独立的,服从负指数分布,类型iAGV 的平均到达速率为Vi;AGV 行驶的平均速率为2dVi,其中d为轨道长度(即两个车间之间的距离)。
不同类型的AGV 具有不同的装载容量,第i类AGV的最大装载容量为Ci。
AGV的调度规则如下:
每台AGV 按照特定的轨道在B2或B4之间往返,即AGV之间相互独立,没有“堵车”现象发生。
AGV在B2装载时服从FCFS规则,当AGV到达B2时,若B2中存在至少一个工件,AGV 装载并立即前往B4,否则等待。
AGV在B4装载时服从FCFS规则,若B4容量达到上限,则等待直至全部卸载并返回B2。
AGV搬运批量的大小是与它的前缓存区(B2)中当前工件数有关的随机变量。建立如图3 所示的节点分解模型,节点1为不与AGV直接相关的节点,是传统的“单个到达、单个离开”的模型;在B2后添加一个加工时间为0的虚拟工作站W2构成节点2,具有“单个到达、批量离开”的特征;将AGV的批处理过程视为加工速率为V的批量处理工作站V3,即节点3,具有“批量到达、批量离开”的特征;节点4 工件转移的特征为“批量到达、单个离开”。
图3 节点分解的排队网模型
通过构造各个节点的状态空间,从而建立状态转移平衡方程,并分析V3节点的状态转移建立节点2 和节点4之间的耦合关联。
文献[18]证明对于k台同等的(相同的运行速度、容量)AGV,可当量为一台具有k倍运行速度的AGV进行分析。为简化描述,以单台AGV为例,建立节点3状态空间如下:
其中,w3表示AGV运载的工件数量。b3表示AGV的状态:b3=-1 为AGV处于“饥饿等待”状态(此时w3=0);b3=0 为 AGV 处于“空载返程”状态(此时w3=0);b3=1 为AGV 处于“负载行程”状态(此时1 ≤w3≤C);b3=2 为 AGV 处于“阻塞等待”状态(此时 1 ≤w3≤C)。该节点状态空间及状态之间的转移如图4所示。
图4 节点3的状态转移示意图
根据系统的 Markov 性,当 AGV 到达B2时,B2中具有n个工件的概率均为P2(n),0 ≤n≤N2,而P2(n)是与AGV状态无关的变量。因此,AGV从空载返程的状态S3(0,0) 转变为饥饿等待状态S3(0,-1) 或者负载行程 状态S3{(w3,1),1 ≤w3≤C} 的 概率 分 别为P3(0)load、P3(w3)load,其中:
因此AGV 状态的转移速率分别为VP3(0)load、VP3(w3)load。
当AGV 在B2处于饥饿等待状态时,AGV 将转变为负载行程状态直至W1中加工完成的一个工件到达B2(有效到达率为,将在下一节分析),因此AGV 由状态S3(0,-1)转变为状态S3(1,1)的转移速率为。
当 AGV 到达B4时,AGV 能否卸货取决于 AGV 当前装载量w3和B4的当前剩余容量k。当k≥w3时,AGV 能完全卸载,因此由负载行程状态S3{(w3,1),1 ≤w3≤C} 转变为空载返程状态S3(0,0) 的转移速率为当k <w时,尚有w-k个工件未能33卸载,AGV 由状态S3(w3,1) 转变到阻塞等待状态S3(w3-k,2)的速率为VPb4(w3)。其中Pb4(k)为B4当前剩余容量为k的概率。
当AGV 在B4处于阻塞等待状态时,AGV 上有w3个工件尚未卸载,随着在W4上完工的工件离开系统(有效加工速率为4),AGV 上未卸载工件逐渐减少直至全部卸载w3=0 并空载返回B2,即由状态S3(w3,2)转变为S3(0,0)的转移速率为。
AGV的整个运行过程是一个“生灭过程”[16],在每个循环中依次经历四种状态。系统稳态时,AGV 处于负载行程状态或空载返程状态的平均时间为T0=T1=V-1,处于饥饿等待状态或阻塞等待状态的平均等待时间为Ti=V-1×PS3(i)/PS3(0),i=-1,2,其中PS3(b3)为AGV处于各个状态S3(w3,b3)的稳态概率。
AGV 在节点2 和节点4 装载和卸载的过程可视为一个交错更新过程[16],且只在AGV 到达的更新时间点发生。对于AGV 单元前节点(节点2),AGV 由离开直至首次到达的间隔时间T2-2=T1+T2+T0内,经历负载行程S3{(w3,1),1 ≤w3≤C}、卸载(阻塞等待)S3{(w3,2),1 ≤w3≤C}、空载返程S3(0,0)三个阶段;对于AGV单元后节点(节点4),AGV由离开直至首次到达的间隔时间T4-4=T0+T-1+T2内,经历空载返程S3(0,0)、装载(饥饿等待)S3(0,-1)、负载行程S3{(w3,1),1 ≤w3≤C}三个阶段。
根据系统Markov 过程的“生灭过程”理论,对每个状态可列出对应的一个状态转移平衡方程,因此各节点可列出与其状态总数一致的状态转移平衡方程组。以状态S3(0,0)为例,输入的状态为S3(w3,1)和S3(1,2),对应的转移速率分别为输出的速率为因此,根据输入与输出状态转移平衡可以建立如下方程组:
同理,节点3共有2(C+1)个状态,针对其他状态建立的转移平衡方程组如下:
对于节点m=1,2,3,4,由其状态转移平衡方程组中未知变量(即各状态稳态概率)PSm(wm,bm)的系数可构成大型稀疏矩阵πm,因此只要求解齐次线性方程组πm xm=0就可获得各状态的稳态概率值。本文应用Gauss-Seidel迭代法[13]迭代求解各节点方程组πm xm=0,由于节点间的耦合关联,当前节点的转移速率参数可由其关联节点在上一次迭代的结果计算得出。在迭代过程收敛后,最终可获得系统各状态的稳态概率值PSm(wm,bm),据此可计算系统性能指标值:平均产出率和生产周期。
优化问题是一个随机非线性整数规划问题,而且约束函数无法用一个封闭的数学方程式表达,因此在满足预设系统产能(输出率)和产品交货期(生产周期)的情况下,采用嵌入排队网模型的禁忌搜索算法优化AGV的配置方案,其算法流程图如图5所示。
步骤1初始化
采用二进制编码,由变量松弛法获取初值。
步骤2解空间搜索
图5 禁忌搜索算法流程图
步骤2.1采用随机选择两个编码的位转变为相反值的扰动策略,产生新解。
步骤2.2调用排队网模型求解系统性能指标值,判断是否满足约束条件,若满足,则进入步骤2.3;否则返回步骤2.1。
步骤2.3通过扰动策略产生多个不同的解,并从中选择最优的解,一旦该操作所得的解被选择,禁止其执行动作的相关操作。
步骤2.4检索是否破禁:当该被禁操作处于被禁准则,但其值所得的目标函数优于历史最优解,则接受该动作执行破禁操作,并更新最优解记录。
步骤3终止判断
判断是否满足终止准则(最高扰动次数),若满足,则当前解为最终解;否则返回步骤2.1。
对于算例求解的操作系统环境为Microsoft Win10,硬件环境为CPU 3.30 GHz,4.00 GB RAM,对于解析算法的求解采用Matlab R2018a软件平台编程实现。
如图6 所示为在Tecnomatix Plant Simulation 8.2上建立的仿真模型示例。基于离散事件仿真理论,在大样本的条件下,构造统计估计量(系统平均产出率和平均生产周期)并计算采集的数据样本。仿真实验每次运行1 000天(仿真时间),并在95%的置信度下进行50次独立重复实验,最终计算系统性能指标的均值。
为验证本文提出近似解法的有效性,并分析AGV配置参数对系统性能指标的影响及规律,本文设计实验方案进行求解,并将结果与仿真实验进行对比分析。
设计自变量(影响因素)如下:小车到达速率V;小车容量C;小车数量n。
设计因变量(实验指标)如下:系统产出率Θ;系统生产周期Γ。
围绕实验目标问题,针对不同因素设计3组单因素实验和1 组双因素实验,具体参数如表1 所示。各组实验算例求解结果如表2 所示,其中“Δ/%”为近似结果相对仿真结果的偏差百分比。与仿真结果对比,本文改进的状态空间分解法对于计算系统平均产出率的偏差率均在±3%内,计算系统平均生产周期的偏差率均在±6%内,并且每次计算所需的时间远远小于仿真所需耗时。因此,可以看出本文近似方法对于系统性能的分析和计算具有良好的适应性。
表1 实验方案的系统参数
图6 系统仿真模型示例
表2 算例求解结果
如图7 所示为各组实验中因变量随自变量的变化规律曲线,由实验结果可以看出:(A)单独增加AGV 的运载速率,系统平均产出率略有提高,而系统平均生产周期则有明显的缩短;(B)单独增加AGV 的运载容量,系统平均产出率略有提高,而系统平均生产周期也有略微的缩短;(C)增加同等并行AGV 的数量,系统平均产出率略有提高,而系统平均生产周期则有明显的缩短;(D)保持AGV的当量运载能力(=AGV运载速率×AGV运载容量)不变的情况下,随着AGV运载速率的提高和AGV运载容量的降低,系统平均产出率有明显的提高,而系统平均生产周期也有明显的缩短。
图7 因变量随自变量变化规律
对于AGV 配置参数的单因素实验(A)~(C),增加AGV 的运载速率或数量,则在制品运输环节的服务速率增大,即在制品的运输时间减少,因此系统产出率增大,而生产周期减少。增加AGV的运载容量,可能增加AGV 每次运输的在制品数量,减少每个在制品的平均运输时间,但是这也可能会提高AGV 在卸载在制品时“阻塞等待”的概率,因此运载容量的变化对系统性能的影响并不显著。对于AGV运载速率和运载容量的双因素实验(D),当量运载能力不变时,提高AGV 速率并减少运载容量,增大了AGV 在运输过程中保持满载的概率,AGV 每次运输的时间缩短,从而AGV 的实际运载能力会有提高,因此系统产出率得以提高,同时缩短了生产周期。
综上所述,合理地配置AGV的数量等参数,可以有效地提高系统产出率并且缩短生产周期。
以图1所示的某装备制造企业生产系统为例,为了在两个车间之间配置AGV 执行物料运输任务,需要确定AGV配置的最优方案。两个车间的距离为d=100 m,如表3 所示,有三种类型的AGV 可供选择。系统产能及订单交货期的约束为:系统产出率不低于Θmin=0.95个/min,而生产周期不超过Γmax=18 min。根据主生产计划的任务投放速率为λ=1.0个/min,两个车间的平均加工速率分别为μ1=1.1个/min 和μ4=1.2个/min。
表3 三类AGV的性能参数
根据嵌入排队网模型的禁忌搜索算法求得算例优化的结果如表4 所示,其中标记“*”的方案为最优配置方案。由于问题规模较小,可以通过枚举搜索的方法获得全局最优解,即逐一增加每类AGV的数量,计算系统性能指标值并判断是否满足约束条件,最后从所有可行解中选择成本最低的为最优解,求解过程的系统性能参数变化曲线如图8 所示。由表4 和图8 的结果可以看出,禁忌搜索算法获得的最优解与枚举搜索获得的全局最优解一致,验证了禁忌搜索方法的有效性。基于系统产出率和生产周期的约束,以最小化AGV 投资成本为目标,配置4台类型II的AGV为该目标下的最优配置方案,以满足系统产能及订单交货期的要求。另外,表4中与仿真结果的相对偏差率进一步验证了改进状态空间分解法的有效性及精确性。由此可以看出本文所提出的求解系统性能指标值及AGV配置优化算法的合理性及有效性,这种方法可以推广到较大规模的系统分析及资源优化配置中。
表4 算例优化结果
图8 优化过程中系统性能的变化曲线
本文建立了随机非线性整数规划模型,描述了智能车间物料搬运系统的AGV数量配置问题。针对具有随机批量搬运的生产系统,建立有限缓冲区开排队网模型,并拓展状态空间分解法计算系统性能指标值,最后将其嵌入禁忌搜索算法,求解AGV配置优化问题。
通过算例求解及结果分析,与仿真结果对比,拓展的状态空间分解法计算系统性能指标值具有较高的精确度,基于该方法求解性能指标值的优化算法可以获得合理的AGV配置最优方案。因此这种方法可以应用于求解较大规模的系统资源优化配置问题。