基于批量分流生产的供应链网络调度研究

赫 超 唐 亮，2

1(沈阳航空航天大学机电工程学院 辽宁 沈阳 110000)2(大连海事大学交通运输工程学院 辽宁 大连 116026)

0 引 言

现如今市场经济模式的快速转变对企业的制造环境产生了翻天覆地的影响。从制造业漫长的发展历史来看，总结制造业的发展阶段，其中主要包括劳动密集型、设备密集型、信息密集型、知识密集型等几个重要的发展阶段。目前制造业正历经着从信息密集型过渡向知识密集型的更科学的发展阶段。在21世纪的新型市场环境中，优势竞争企业应该具备的最核心能力是对市场的快速响应速度，因此，以协作分工和高度自适应性为特征的协同制造模式通过数字化、网络化信息与其他优势企业在生产经营业务活动中各个环节进行合作，以实现供应链中各个企业间的生产组合利用和资源共享。协同制造能根据用户的不同需要迅速地与其他企业组成联盟，进行产品的协同设计和制造，从而快速响应用户的需求。这使得拥有协同生产能力的制造商们共同组成了关于产品生产的供应链。

目前在生产调度、智能调度优化方面，国内外许多学者已经从多个角度做出了大量有意义的工作并获取了丰硕的研究成果[1-6]。近年来，研究者不再局限于仅对生产企业内部生产环节进行优化，而是越来越多地从供应链的角度综合考虑包括原料采购、生产、库存以及物流配送等环节在内的集成优化。这些研究[7-9]从不同的角度对集成调度模型进行了构建，即制造企业的不同分布结构、交货时间约定、交货方式、车辆运输方式等，并考虑了总成本和客户服务水平之间的权衡。总体而言，该领域大部分研究可以分为两阶段：供应链调度[10-14]和三阶段供应链调度[15-22]。

归纳现有针对供应链调度问题的文献，我们发现针对协同生产方式的局部特点所开展的研究较少。而结合实际生产过程，我们发现当供应链中某一制造企业完成产品中间工序的加工任务后，可能会将手中的半成品以成批量的形式委派到具有相应生产能力的多个下游企业处进行生产。这种协同制造方式与上游制造企业将全部数量的半成品委任给单个下游制造企业的传统协同制造模式相比，更加符合以快速响应客户需求为核心的市场经济模式，并使各个制造企业的优势生产资源得到更充分的利用。

基于上述内容，本文将针对基于批量分流生产的供应链网络调度问题进行研究，建立了以总最短完工时间为优化目标的数学模型。此外，对于求解算法的选择上，由于本文研究问题的数学模型中0-1变量较多，其中不仅包括协同企业选择变量，还包括生产路径选择变量。而Benders算法是一种专门求解问题中含有不同类型变量的分解算法，故本文选用Benders分解算法对问题模型中0-1变量和纯整数变量进行分离后,通过反复求解对偶子问题和Benders主问题不断修正上下界值最终得到问题的最优解。此外，为提高算法求解效率，通过应用pareto最优割代替Benders最优割的方法对传统Benders算法进行改进。

1 数学模型及求解算法

1.1 问题描述

在协同制造模式下，企业通常以标准化和模块化生产运作，具备一种或几种优势生产能力。为了充分利用这些协同企业的某种优势资源，对产品关键部件的协同工序进行分解，不同的协同工序可以交由具备该工序生产能力的协同企业完成。同时，由于具备某种协同工序生产能力的协同企业可能有多个，因而构成协同供应链有向网络，网络中包含多条可行加工路径。一旦出现多种不同类型订单需求的情况，由于其工序的不同，协同企业会发生相应变化，导致部分协同企业处于不同协同供应链网络的现象，因此有必要扩展至多个交互的协同供应链网络进行研究。针对本文供应链网络调度问题给出如下假设：

1) 每个协同企业可完成某一项协同工序或几项协同工序的加工；

2) 各个订单的每个协同工序只能由一个协同企业完成，且同一个协同企业不能并行处理不同订单的生产任务；

3) 各个订单按照先到达先加工规则，加工结束立即运输；

4) 各个订单产品的中间工序在加工过程不允许打断,也不存在返回再加工的情况；

5) 整个生产系统运作开始时间为0时刻；

6) 订单必须成批生产，不接受某几个产品的单独生产。

1.2 模型构建

假定多个协同制造企业组成的供应链现已接受多个订单的生产，其中包含多个种类的产品。为了构建分批制造的供应链调度问题模型，首先分别对模型中参数与变量进行描述。

在相关参数中，K={1,2,…,k}为全部订单集合，k为订单总数；Ni={1,2,…,ni}为第i个订单全部可行加工路径集合，ni为第i个订单可行加工路径总数；M={1,2,…,m}为全部协同企业集合，m为协同企业总数；Tij(p,q)为运输时间矩阵，表示i订单第j条可行路径的加工矩阵中第p行q列元素，若i订单在第j条可行路径上企业p与企业q之间存在先后相邻加工关系则大于0，否则为0；b为每批产品所包含的产品数量；di为订单i包含的产品数量；elij和slij分别代表订单i在可行加工路径j处首工序和尾工序制造企业序号；wip为订单i在协同企业p处单位产品的加工时间。

在相关变量中，yij为0-1变量，i订单选择在第j条可行路径上加工为1，否则为0；xij表示i订单在第j条可行路径上的加工批量；cijp表示i订单在第j条路径上协同企业p的完工时间；lii′jj′p为0-1变量，若订单i、i′分别在第j、j′条可行路径上加工且在这两个加工路径上的交叉企业p处，订单i先于订单i′完工则为1，否则为0；ei表示订单i的完工时间。综上所述，基于批量分流生产的供应链网络调度问题模型构建如下：

(1)

(2)

b×xij≤di×yij∀i∈K,j∈Ni

(3)

b×xij≥b-di×(1-yij) ∀i∈K,j∈Ni

(4)

yij≤cijp∀i∈K,j∈Ni,p=elij

(5)

cijp≤G×yij∀i∈K,j∈Ni,p=elij

(6)

cijp≥b×xij×wip∀i∈K,j∈Ni,p=slij

(7)

cijp+t(p,q)×yij≤cijq-b×xij×wiq

∀Tij(p,q)>0

(8)

cijp≤ci′j′p-b×xi′j′×wi′p+G×(1-lii′jj′p)+

G×(2-yij-yi′j′)

∀i、i′∈K,j∈Ni,j′∈Ni′,p∈Hij∩Hi′j′,i≥i′,i≠i′‖j≠j′

(9)

ci′j′p≤cijp-b×xij×wip+G×(2-yij-yi′j′+lii′jj′p) ∀i、i′∈K,j∈Ni,j′∈Ni′,p∈Hij∩Hi′j′,

i≥i′,i≠i′‖j≠j′

(10)

ei≥cijp∀i∈K,j∈Ni,p=elij

(11)

模型中目标函数式(1)表示所有订单的完工时间总和的最小值；式(2)表示同一订单在不同的路径上加工数量的总和必须达到既定数量；式(3)和式(4)表示上游企业将半成品后续工序的工作托付给下游企业时，若选择的下游企业不止一个，则分流的批量至少为1;式(5)和式(6)表示产品在所选择的可行加工路径中的各个协同企业的完工时间均大于0，在未被选择的路径上均等于0，其中G代表一个很大的正整数；式(7)限制了订单首加工工序的最小完工时间，保证了整个生产系统运作开始时间为0时刻；式(8)结合了运输时间限制了同一路径中相邻两个企业的完工时间关系；式(9)和式(10)表示企业不能同时加工两个订单，这里要额外说明一点，为了减少变量lii′jj′p以及相关约束的数目，我们约束变量lii′jj′p在满足i≤i′的前提下，必须额外满足i≠i′或者j≠j′;式(11)表示产品最终的完工时间大于部分数量产品的完工时间。

2 改进Benders分解算法及其应用

2.1 Benders分解算法及其改进

Benders分解算法是J.F.Benders在1962年提出的，它是一种专门求解混合整数规划的算法。其主要思想是通过固定一部分变量的值将原始复杂问题转换成相对简单的子问题，固定变量的取值不同相应产生的子问题也不同。每次求解子问题都会根据不同求解情况得到一个Benders最优割或Benders可行割加入到新构建的Benders主问题中，通过一定次数的迭代不断缩小最优解的搜索范围后得到最优解。以如下形式的0-1混合整数规划为例说明Benders算法流程：

mincTx+dTy

(12)

s.t.Ay+Bx≤b

y∈{0,1}x∈Z+

式中：c、d、b分别为m、n、l维列向量，A、B分别为l×n、l×m维系数矩阵，当变量y取值为y*时，式(12)可以转变为以下形式：

min{dTy*+min{cTx|Bx≤b-Ay*,x∈Z+}}

(13)

式(13)中内层问题min{cTx|Bx≤b-Ay*,x∈Z+}被称为原始子问题，其对偶问题Q为：

max(b-Ay*)Tλ

(14)

s.t.BTλ≥c

λ∈Z+

记S={λ|BTλ≥c,λ∈Z+}是由对偶子问题Q的可行域构成的多面体，它不随y取值的不同而发生改变，up和ud分别为多面体S的极点集与极方向集。通过引用一个辅助变量θ构造Benders主问题，其最优值与原始问题相同,具体形式如下：

minθ

(15)

s.t. (b-Ay)α≤θα∈up

(b-Ay)β≤0β∈ud

y∈{0,1}θ∈Z+

α和β分别多面体S的极点和极方向，每个极点或极方向在Benders主问题中都代表一条约束，由极点产生的约束叫做Benders最优割平面，由极方向产生的约束叫做Benders可行割平面。当多面体S全部极点与极方向被添加进up和ud中，求解问题Benders主问题得到的目标函数值与原始问题的最优值相等。但在实际求解过程中想要得到全部极点和极方向是很难做到的，尤其当问题规模较大时，多面体S中极点和极方向的总数量也随着约束数大幅增长，而且在这些约束中只有一部分是积极约束，故没有必要求出全部极点与极方向。

以下给出一种基于Benders表示的约束生成算法：在求解前需要设置初始上UB、下界LB和程序终止阈值ε，并且给定一组列变量y的初始可行解y*，初始up和ud为空集。在求解过程中，首先求解y值初始固定后的对偶子问题Q，并且根据对偶子问题Q的求解情况判断添加Benders最优割或可行割，若此时对偶子问题存在最优解，则将λ添加进up中，并且根据对偶定理，其目标函数与原问题目标函数值相等，故可以将该值作为临时上界；若对偶子问题存在无界解，则将λ添加进up中。在更新up与ud中元素后下一步求解Benders主问题式(15)，将求解得到的θ最优值作为临时下界，并用最优解y**替换给定的初始值y*产参与下一次流程运算。经过反复的迭代，上下界逐渐逼近，当两者逼近到一定程度满足终止条件时，问题求解结束并得到原始问题最优值与最优解。

上述过程为传统Benders算法计算流程，此流程中不需要求得多面体S中全部极点和极方向。但传统Benders算法在每次迭代过程中无法保证每次得到的Benders最优割的质量，从而影响上下界的收敛速度。故本文利用pareto最优割替换Benders最优割提高算法求解效率，并重新安排算法流程。pareto最优割也是在求解对偶子问题的过程中得到的，具体方式如下：

传统Benders算法中，当对偶问题Q有唯一解时，其解对应一个极点被添加到Benders主问题中，在改进Benders算法中此部分有所不同，具体流程如下：

首先我们假设在对偶子问题中y取值为y1。按照算法流程接下来计算对应的对偶子问题：

max(b-Ay1)Tλ

(16)

s.t.BTλ≥c

λ∈Z+

当问题式(16)存在最优解时,记最优解集为G(y1)。下面通过引入一个新的y值y2构造新的问题：

max(b-Ay2)Tλ

(17)

s.t.BTλ≥c

λ∈Z+,λ∈G(y1)

问题式(17)必定存在最优解，此时求解得到的λ值成为约束添加到Benders主问题中称为pareto最优割。

2.2 基于改进算法求解供应链调度模型

前文阐述了研究问题的数学模型以及其中各部分约束含义，下面结合研究问题阐述计算流程。

在利用改进Benders算法求解之前我们要将原问题模型转化成形如式(12)较为规范又有明确划分的矩阵形式。首先我们要得到系数矩阵A、B和常数列向量。本文将约束式(2)-式(10)中各项系数转换为统一形式的系数矩阵MA，其中不同变量对应着相应的列，不同约束对应着相应的行。若某条约束中不包含个别变量，则在MA中对应行列的系数值为0。比如在式(3)中并不包含变量yij、lii′jj′p、cijp、ei,那么在MA中代表式(3)所有展开的行中，以上变量对应的列中系数均为0。对于式(2)来说，我们可以将其转换为双向关系约束，即将“等于”转换为“大于等于”和“小于等于”。常数列向量则是各个约束中的常数移到不等式右边得到的，通过上述的方法很容易得到，这样我们就可以得到一个系数矩阵MA和常数列向量。

下一步我们将矩阵根据变量的不同类型将其划分成两个部分，即MA=(A，B)。A为0-1变量列代表的MA中部分矩阵，B为纯整数变量列代表的MA中部分矩阵。在本文建立的数学模型中，0-1变量包括yij、lii′jj′p，纯整数变量包括xij、ei、cijp。可以根据实际问题的规模将这两种变量对应的列分离成矩阵A和矩阵B。目标函数中系数行向量也以此方法得到，就本文模型而言，目标函数中系数行向量中除变量ei所对应的列的系数为1，其余列系数均为0。最终我们得到形如问题式(12)的0-1混合整数规划：

(18)

s.t.AY+BX≤r3

Y∈{0,1}X∈Z+

式中:Y=(y,l),X=(x,e,c),且均为列变量。

综上所述，改进Benders算法具体流程如下：

Step1设置初始UB、LB和程序终止阈值ε，并赋予两组0-1列变量(y*、l*)和(y**、l**)初始值，up和ud均为空集。转Step2。

Step2求解利用(y*，l*)产生的问题的对偶问题Q。若问题Q存在最优解，将(y*，l*)值赋予(y**，l**)。并相应得到一个pareto最优割平面，同时更新up。此时子问题最优值若小于当前UB则其值代替UB，并通过引入列变量(optimy，optiml)记录当前最优解(y*，l*)；若子问题存在无界解，则相应得到一个Benders可行割平面，并更新ud。转Step3；

Step3根据当前UB、LB与ε值判断程序是否终止。若满足终止条件，则UB为原始问题最优值。转Step5；否则转Step4。

Step4求解更新up和ud后的Benders主问题，根据θ最优值修正当前LB，并将最优解赋给(y*,l*)。转Step2；

Step5(optimy,optiml)为原始问题中列变量(y,l)最优值，分别代入原始问题约束式(2)-式(10)中，此时求解原始问题最优解，问题求解结束。

算法具体流程如图1所示。

图1 算法流程图

3 仿真算例与分析

3.1 参数设计

(1) 协同制造网络设计 一般来说，不同类型产品加工网络差距较大。结合实际生产情况，本文设计3种不同类型的加工网络。每类加工网络中每一个节点都代表不同的协同制造企业，每个企业都可以根据自身生产资源参与不同订单产品的加工。图2为平衡型加工网络。这种网络是生产比较平衡的一种情况，每个网络节点都有多个均匀的上下游节点可供选择，一般此类网络的物流比较平均(规则网络)，网络压力较小。图3为瓶颈型加工网络。瓶颈型加工网络表明由于某些制造工序的特殊性，网络中可行加工路径必须经过某些节点，因此生产网络的总生产能力会受到这些节点的制约。图4为跳跃型加工网络。对于跳跃型加工网络，主要是考虑到不同的企业的加工能力不同，针对相同的产品有的企业可以完成连续多个工序的加工任务，由此一旦加工路线经过此类网络节点，可以较快到达网络终点。

图2 平衡型加工网络

图3 瓶颈型加工网络

图4 跳跃型加工网络

(2) 相关加工参数设计 在本文的仿真算例中假设接受生产三种不同类型的产品，平衡型、瓶颈型、跳跃型三种协同制造网络各代表一种产品，可行路径数目为8、8、4。这三种产品的需求数量分别为30、40、30，每10个产品为1批，共11家协同制造企业参与提供生产资源。每家企业生产不同订单单件产品所需时间与企业间运输时间分别如表1、表2所示。

表1 不同企业对不同订单加工时间 h

表2 企业间运输时间 h

鉴于篇幅，本文将全部订单的每条可行加工路径中各个协同企业间的运输时间总结在一起。为了避免混淆，下面本文通过例子详细说明，比如本文将生产路径m2→m5→m3→m7→m11作为订单1的第1条生产路径，则在矩阵T11(p,q)中只有2行5列、5行3列、3行7列为非0值，其值等于表2中相应位置值(表头除外)。

3.2 仿真结果描述

表3 具体调度方案

以上具体求解结果中，订单1-订单3的最终完工时间分别360、400、235 h，总最小完工时间为995 h。从路径y15、y18中在企业6和后续企业的加工时段可以看出，企业6在连续生产20件订单1产品后，分别向企业3、8分配了10件产品的生产任务，这将减少产品在协同制造网络中的流通时间。为了验证分流协同加工方式可以有效缩短产品在生产中流通时间，我们取消分流生产方式后求解该仿真算例，得到最小总加工时间为1 840 h。在算法性能上，我们也使用传统Benders算法对算例进行求解，并与基于pareto最优割改进的Benders算法进行比较。通过仿真，传统Benders算法耗时22 337 s，迭代次数为5 728次；改进Benders算法耗时14 152 s，迭代次数3 819次。仿真结果表明改进算法具有明显的优势，从而验证了改进算法的有效性。

4 结 语

本文以多种产品协同制造网络对供应链调度问题进行了研究。为了使研究更具普适性，设计三种不同类型的协同制造网络。同时，本文考虑产品可以成批量分流加工，以总最小完工期为目标建立了0-1混合整数规划模型。最后针对本文问题特点采用Benders分解算法进行精确求解。为了提高求解效率,本文利用pareto最优割替换Benders最优割改进传统Benders算法。仿真结果表明,利用改进Benders算法可有效提高问题求解效率。