CN107544251B - 一种基于分布式鲁棒模型的最小化总拖期的单机调度方法 - Google Patents

Abstract

本发明提出一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，属于生产调度和生产内外资源优化技术领域。该方法首先建立针对单机调度问题的分布式鲁棒优化模型，得到模型的目标函数表达式；然后将分布式鲁棒优化模型转化为整数二阶锥规划模型；对转化后的模型求解，将所有工件加工序列的排列组合通过枚举的方式在一个搜索树中进行表示，通过分支定界算法对搜索树进行剪枝，最终得到最小化总拖期的最优单机调度方案。本发明将生产环境中的不确定因素考虑在内，使得模型相比于假设生产环境都是确定的确定性单机模型更加符合实际生产状况，得到的调度方案能更好地应用于实际生产中。

Description

一种基于分布式鲁棒模型的最小化总拖期的单机调度方法

技术领域

本发明属于生产调度和生产内外资源优化技术领域，特别涉及一种考虑了在加工时间不确定的情况下基于分布式鲁棒模型的最小化总拖期的单机调度方法。

背景技术

生产调度问题是制造业和服务业中不可或缺的关键问题，其旨在通过合理配置和优化内外资源，从而缩短制造周期，降低生产成本，改善设备利用率，对于众多企业具有重大意义。因此生产调度自然成为了多种先进制造模式共同关注的核心内容。而随着敏捷制造、智能制造等概念的提出，人们对于产品交付的准时性的需求愈加强烈，产品的延迟交付可能会导致销售损失和顾客流失，同时也会增加企业的库存成本。基于单机模型的最小化拖期调度问题，在调度问题中最具有代表性，同时单机调度的很多理论也可以使用在其他调度模型中，因此最小化总拖期的单机调度问题具有重要的研究价值和实际意义。

过去的几十年中，大多数最小化总拖期的单机调度问题都是假设生产环境是确定的，因此建立的模型属于确定性模型。在确定性模型中，都假定生产环境中所有因素都是已知的、人为可控的，如生产时间确定，机器一直稳定运行不会存在故障等。然而实际的生产环境中往往充满了不确定因素，如机器故障，订单变动等。调度过程中参数的不确定性往往会导致原有的最优方案不可行，从而使得按照确定性模型得到的调度方案往往无法达到决策者事先的预期。因此，针对不确定环境下的生产调度问题研究逐渐引起了关注。主流研究方向包括随机规划和鲁棒优化。

在随机规划模型中，不确定的参数通常会被看作是一个分布已知的随机变量，模型的目标通常基于随机变量的分布信息优化目标函数的长期意义下的期望或者方差。从80年代开始，随机规划就陆续在生产调度领域展开使用。但是随机规划模型存在一些固有的缺点从而使得其在实际生产过程中很难得到较为广泛的应用，主要体现在：1)在随机规划模型中，需要获得参数的确切分布信息，而实际的生产过程中的参数往往是不知道精确的分布信息的，错误的分布估计很可能导致理论与实际应用的差距；2)随机规划通常优化的是系统长期的平均意义下的性能指标，但并不适合用于解决某一次生产调度过程中的优化问题。由于随机规划存在上述不足，鲁棒优化模型在生产调度领域的应用逐渐展开。

早在1950年，统计学家Wald提出悲观决策准则，其中即包含了鲁棒优化的思想，即要求决策者根据每一种方案的最坏情况进行选择。而在20世纪70年代后，鲁棒优化得到迅速发展，形成了独立完整的理论体系，广泛应用于控制、金融、供应链等领域。鲁棒优化使用区间数据对不确定参数进行描述，一个情景代表不确定参数的一种可能的取值，这种描述方法相对于随机规划模型中利用分布函数对模型参数进行的描述的方式更为简单且符合实际。相对于随机规划而言，鲁棒优化具有以下特点：1)决策关注的是不确定参数的边界情况，不需要知道参数的精确分布信息；2)鲁棒优化模型通常可以转换为等价的确定性模型进行求解，求解规模相对于随机规划方法相对较小；3)由于鲁棒优化的决策针对的是参数的最坏情况，得到的解具有一定的保守性。可以通过变化控制参数的值，弹性调节模型的鲁棒性从而控制解的保守水平。

随着鲁棒优化的发展，许多研究者试图将鲁棒优化与随机规划方法联系起来，在鲁棒优化中考虑更多的信息，减少针对参数最坏情况进行规划的保守性。一种很自然的思路即虽然无法给出参数的精确分布信息，但是关于概率分布的一些信息，如一阶矩和二阶矩信息，将满足这些概率分布信息的概率分布组成一个集合，将原有鲁棒优化中的不确定集换成概率分布集合，目标函数换成最坏情况下期望的最优值，这种模型叫做分布式鲁棒优化模型，其即考虑了分布信息，减少了传统鲁棒优化的保守性，同时相比于随机规划又更易求解。

在分布式鲁棒优化模型中，参数的不确定性通过分布集来表示，不确定参数被看作一个随机变量，其精确的分布函数未知，但是分布信息如一阶矩、二阶矩信息等已知，因此该随机变量的分布函数可以是某个特定分布函数集中的任何一个。而在目标函数中，需要选取一种性能测度来评价具有随机性的目标函数值，目前多数研究都选取了期望、方差或者条件风险价值作为测度。在生产调度问题当中，Richard L.Daniels最先将鲁棒优化的思想引入生产调度问题当中，近几年鲁棒优化在单机、并行机以及混合流水车间等模型中都有相应的发展。然而目前绝大多数的研究都是基于不确定性集的鲁棒优化，不确定性集通常为有限的离散集合或者是连续的区间形式。使用分布式鲁棒优化解决生产调度问题的研究相对较少，而其中以最小化总拖期为目标的研究更少。

目前已有技术方案中，都是采用传统的鲁棒优化方法求解实际问题，如“一种不确定性环境下发电计划的鲁棒优化方法”，其在描述不确定环境下的随机因素时，采用传统鲁棒优化方法的思路，即确定发电计划优化模型中的各不确定参数，通过参数波动变化区间对不确定参数进行描述，从而建立发电计划的不确定优化模型。然而采用参数的波动变化区间对不确定参数进行描述，进而使用传统鲁棒优化方法进行优化，得到的结果往往会过于保守，因为这种方法通常只考虑到了不确定参数的极端情况，即区间的两个端点，而在实际的生产过程中，这种情况出现的概率十分微小，因此针对极端情况得到的解往往会因为考虑了鲁棒性而损失很多性能。在实际应用领域方面，目前还没有将分布式鲁棒优化方法应用于最小化总拖期的单机调度问题上的例子。

发明内容

本发明的目的是为克服已有技术在建模方面的不足之处，提出一种基于分布式鲁棒模型的最小化总拖期的单机调度方法。本发明将生产环境中的不确定因素考虑在内，使得模型相比于假设生产环境都是确定的确定性单机模型更加符合实际生产状况，得到的调度方案能更好地应用于实际生产中。

本发明提出的一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，其特征在于，包括以下步骤：

1)建立针对单机调度问题的分布式鲁棒优化模型；

1-1)确定分布式鲁棒优化模型的不确定参数；

设不确定参数为工件的加工时间，假设有N个工件构成工件集合为N＝{1,2,…,N}，所有工件的加工时间构成一个随机向量p＝{p₁,...,p_N}，其中p₁,...,p_N分别表示工件1,...,N的加工时间；该随机向量的分布P^p未知，但是属于一个均值向量和协方差矩阵确定的分布集中，该分布集的定义如式(1)所示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (1)

其中P^p为任意一个符合均值向量为μ＝{μ₁,...,μ_N}、协方差矩阵为Σ＝diag{σ₁ ²,...,σ_N ²}的分布，μ₁,...,μ_N表示工件1,...,N的加工时间的平均值，σ₁ ²,...,σ_N ²表示表示工件1,...,N的加工时间的方差，diag{·}表示对角阵；E(.)和Cov(.)分别表示加工时间的均值向量和协方差矩阵；

1-2)确定模型目标函数与约束条件；

1-2-1)确定模型的目标函数；

在给定一个调度方案的情况下，所有工件的拖期之和可以表示为如式(2)所示：

其中y＝{y_ij,i,j＝1,...N}为模型的决策变量，一个向量y对应一个可行的调度方案，如果工件i在工件j之前加工，则y_ij＝1，否则y_ij＝0；工件j的完工时间表示为

d_j为工件j的交货期，因此每个工件的拖期表示为

考虑最坏分布情况下的总拖期之和的期望值，表达式如式(3)所示：

sup_p～(μ,Σ)E[f(p,y)] (3)

其中sup表示求取集合的上确界，p～(μ,Σ)表示所有工件加工时间向量p属于均值向量为μ和协方差矩阵为Σ的分布集，E表示期望；

模型的目标是通过求解获得一个最优的调度方案y^*，使得在该调度方案下，最坏分布情况下的总拖期之和的期望值最小，则模型的目标函数表达式如下：

y^*＝argmin_ysup_p～(μ,Σ)E[f(p,y)] (4)

1-2-2)确定模型的约束条件；

1-2-2-1)随机加工时间约束；

所有工件的加工时间p的分布未知，但属于一个均值向量、协方差矩阵已知的分布集中，表达式如式(5)所示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (5)

1-2-2-2)可行加工序列位置约束；

两个工件之间有先后顺序，任意多个工件之间不能在先后顺序上出现不合理的情况，表达式如式(6)和(7)所示：

y_ij+y_ji＝1,i,j＝1,...,N,i≠j (6)

y_ij+y_jk+y_ki≤2,i,j,k＝1,...,N,i≠j,j≠k,k≠i (7)

1-2-2-3)可行调度方案约束；

任意可行调度方案y中的每一个元素均为0-1变量，表达式如下：

y_ij∈{0,1},i≠j (8)

如式(6)-(8)表示的约束条件均为调度方案的可行性约束，构成了调度方案的可行域，表达式如式(9)所示：

1-3)建立基于分布式鲁棒优化模型的最小化总拖期的单机调度的数学表达式；表达式如下：

其中，式(10)为分布式鲁棒优化模型的目标函数，式(11)为步骤1-2)中的约束条件；

2)对分布式鲁棒优化模型进行转化；

将步骤1)建立的分布式鲁棒优化模型转化为整数二阶锥规划模型，具体步骤如下：

2-1)根据步骤1)得到的分布式鲁棒优化模型的目标函数确定分布式鲁棒优化模型的上界，表达式如下：

其中，

为工件的延迟时间，拖期的定义为工件延迟时间和0的较大值，即max{0,L_j}；

如式(12)所示的不等式将步骤1)中分布式鲁棒优化模型的目标函数转化为了模型的上界，从而将对分布式鲁棒优化模型的求解转化为对上界的求解；

2-2)根据步骤2-1)的结果，将步骤1)建立的分布式鲁棒优化模型转化为一个整数二阶锥规划模型，整数二阶锥规划模型的目标函数的表达式如下：

3)对步骤2)转化得到的整数二阶锥规划模型进行求解，使用分支定界算法得到最优单机调度方案；具体步骤如下：

3-1)构建搜索树；

将所有工件加工序列的排列组合通过枚举的方式在一个搜索树中进行表示；搜索树中，除根结点外每一个结点代表一个工件，每一个工件的所有可能的下一个工件构成了分支，从树的根结点到任意一个叶子结点所对应的一条路径即为所有工件加工顺序的一种组合，对应一个生产调度的可行解；

3-2)通过采用基于拉格朗日松弛的上界初始值估计方法得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

基于拉格朗日松弛的上界初始值估计方法采用基于2-index模型的拉格朗日的上界初始值估计方法或基于3-index模型的拉格朗日的上界初始值估计方法；具体如下：

3-2-1)基于2-index模型的拉格朗日松弛的上界初始值估计方法；

定义

w_ij＝σ_iy_ij，则式(13)所示整数二阶锥规划模型的目标函数等价为式(14)所示：

同时在约束条件式(11)基础上增加新的约束条件如下：

利用拉格朗日松弛方法对约束式(15)和式(16)进行松弛，分别引入拉格朗日乘子λ_0j,j∈N以及λ_ij，i,j∈N，记

其中

记

其中

拉格朗日函数写为：

因此，松弛约束式(15)和式(16)之后的混合整数二阶锥规划模型的目标函数写为：

式(18)所示目标函数与约束条件式(11)构成了一个新的优化模型，该优化模型称为拉格朗日对偶问题；

在给定Λ的情况下，如式(18)所示的目标函数分解为两个优化函数：等号右边加号之前的表达式为第一个优化函数，该优化函数求解采用0-1整数线性规划求解，加号之后的表达式为第二个优化函数通过下式求解：

仅考虑式(19)中的第二种情况，即λ_j≤1,j∈N；通过拉格朗日松弛方法最终得到如式(20)所示的拉格朗日优化模型：

其中，

采用约束生成的方法对式(20)进行求解，得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

3-2-2)基于3-index模型的拉格朗日松弛的上界初始值估计方法；

3-index模型的目标函数为：

约束条件包括：

其中，

表示当且仅当工件j排列在第k个位置且紧接着工件i进行加工，否则

C_k，d_k和

分别代表在第k个位置工件的完工时间，交货期以及方差；约束式(29)保证每一个工件只占据一个位置且只有一个前继工件；约束式(30)和约束式(32)表示第一个位置和最后一个位置的工件满足的条件；约束式(31)保证所有工件的前后顺序不出现环；约束式(34)-式(36)是C_k，d_k和

的定义；

记w_k＝C_k-d_k，

记X为所有可行解

i＝0,...,N,j,k＝1,...,N+1的集合；则3-index模型重写为如下形式：

x∈X (43)

松弛约束条件式(38)和式(39)，引入拉格朗日乘子θ_k,k＝1,...,N和

则拉格朗日函数如式(44)表示：

其中，

3-index的拉格朗日对偶问题写为：

对式(49)求解，求解得到的x与式(13)所示目标函数中的y是一一对应的，因此x也代表调度方案的一个可行解；将x转化为y，得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

3-3)对搜索树进行预搜索；

枚举搜索树除根结点层外前P层的所有工件序列组合，将前P层每种工件序列代入目标函数式(13)计算，保留使得目标函数值最小的一种前P层的工件序列；

3-4)确定前P层的工件序列后，对于搜索树第P+1层所有的分支进行判定并对搜索树剪枝；具体步骤如下：

3-4-1)采用问题性质作为支配准则对工件进行判定，去掉不满足支配准则的工件，从而对搜索树进行剪枝；

在加工时间为不确定参数的情况下，采用以下三条作为工作排序的支配准则，具体如下：

性质1：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k以及u_j-d_j＜u_k-d_k，则在任意最优的调度方案中，工件j在工件k之前加工；

性质2：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k，同时在某个调度方案S中满足|C_max-d_j|＜|C_max-d_k|，其中C_max＝max{C_j(S),C_k(S)}，C_j(S)和C_k(S)表示在调度方案S中工件j和工件k的完工时间，则在任意最优的调度方案中，工件j在工件k之前加工；

性质3：如果工件j和工件k属于未调度的工件集合U，同时工件j位于所有未调度工件的最后一位进行加工，则在任意一个调度方案S中，对于所有未进行调度的工件k∈U，都满足：

其中，

其中，B_k(S)表示在调度方案S中在工件k之前加工的工件集合；

采用性质1到性质3作为支配准则，剪掉搜索树中违反支配准则的分支；

3-4-2)采用基于序列的下界估计方法对搜索树中分支的目标函数值下界进行估计从而进行剪枝；

在搜索过程中，将某一个未加工的工件置于第P+1层，记为结点a，此时记仍未加工的工件集合为R，利用下界估计方法对集合R对应的目标函数值的下界进行估计，当该下界与已经排序完成的工件序列对应的目标函数值之和超出了搜索树当前最优可行解对应的目标函数值时，则该未加工工件不能位于P+1层，停止对结点a及其之后的所有分支的搜索，进而去搜索其他结点；具体如下：

式(13)的目标函数重写为如下所示：

其中

根据式(52)中RT_j的定义，RT_j(t,Σ²)是分别关于t和Σ²的增函数，因此构造如式(56)所示目标函数：

如果满足E_j(S)≤C_j(S)，

则式(56)的值是原目标函数式(13)的一个下界；定义：

其中μ_[i]和

是工件加工时间的均值和方差按从小到大排列位于第i位的取值，k₁和k₂分别是工件j在两种排列下所在的位置，k是工件j在某一个调度方案S中的位置；

3-4-3)利用搜索树的当前上界，对第P+1层所有分支进行判定：

将某一个未排序的工件置于第P+1层，记为结点b，若前P+1层工件序列对应的式(13)的目标函数值大于搜索树当前上界，则该工件不能位于当前的P+1层，停止对结点b及其之后的所有分支的搜索，进而去搜索其他结点；

3-5)对第P+1层所有分支按照步骤3-4)遍历完毕后，对于该层所有剩余的分支，对每个分支计算μ_i-d_i,i＝1,...,N；并按照差值从小到大的顺序依次从对应的分支出发，沿搜索树继续下一层工件的搜索；

3-6)重复步骤3-4)至3-5)，依次对搜索树每一层的分支进行判定并排序；

3-7)当到达搜索树的最后一层时，此时搜索得到一个完整的工件序列作为一个可行解，计算该可行解对应的式(13)的目标函数值并判定：若计算得到的该可行解的目标函数值小于当前上界，则将该可行解作为新的当前最优可行解，并将该可行解对应的目标函数值作为新的当前上界；

3-8)遍历搜索树的每一个分支；遍历完毕后，最后保留的当前最优可行解即为式(13)所示的目标函数的最优可行解，记为y^*，该最优可行解为一个具有最优目标函数值的工件序列，该工件序列为最终得到最优单机调度方案。

本发明的特点及有益效果在于：

本发明提出的一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，分布式鲁棒模型结合了传统基于不确定性集的鲁棒优化和随机规划的优点，在模型中不确定参数被看作是一个随机变量，相比于随机规划其不需要精确的分布信息，相比于传统鲁棒优化，其又考虑了分布集的均值和方差信息，从而降低了决策的保守性。这种方法更适用于实际的生产过程，保证了系统鲁棒性，降低决策的风险，同时给出的调度方案又具有一定的有效性，不过于保守。具体的优点可以总结如下：

(1)相比于随机规划，本发明方法采用分布式鲁棒模型不需要不确定参数精确的分布信息，只需要估计出一阶矩和二阶矩信息。随机规划经常假设数据服从一定的分布，如指数分布，正态分布等，然而在实际生产过程中，数据的分布并不能直接得到，通常需要大量的历史数据进行估计，估计的准确性又会影响随机规划结果的好坏，而分布式鲁棒优化方法只需要数据的一阶矩和二阶矩信息，这在实际生产过程中是十分容易获得的，因此更为实用，也更容易实现；

(2)传统基于区间的鲁棒优化方法只考虑了数据的分布区间虽然不需要知道分布信息，但同时也没有很好利用已知的信息，如一阶矩和二阶矩信息，从而使得得到的调度方案过于保守，影响生产效率，本发明中分布式鲁棒模型的优化方法减小了调度方案的保守性，使得其拥有更好的系统性能；

(3)算法方面，分析得到的问题性质可以有效的作为分支定界方法的支配准则，大幅提升分支定界算法的速度，从而使得算法能够对适当规模的问题进行快速求解，通过与商业求解器CPLEX进行对比，该算法具有更为明显的速度优势。

(4)算法中提出的基于不同模型的拉格朗日上下界估计方法保留了解的可行性，因此在得到上下界估计的同时，也能提供一个很好的初始解，可以作为分支定界算法一个很紧的上界，同时也可以作为启发式算法的初始解；

(5)目前以最小化总拖期为目标函数的单机调度问题绝大多数考虑的是确定性模型，在实际生产过程中，通常会存在很多的不确定因素，如加工时间的不确定，机器故障等。这种情况下，确定性模型得到的调度方案往往是不可行的，即使可行往往也不是一个很好的调度方案。而本发明得到的是考虑不确定因素的性能较为鲁棒的解，其对于不确定因素的任何取值都保持有可行性，同时即使是在最坏的情况下也有较好的性能。因此本发明得到的鲁棒调度方案即使在有不确定情况发生的时候依旧可以维持方案的可行性，同时保证调度方案的性能，从而减少人工调控带来的不便和性能损失。

附图说明

图1为本发明方法的整体流程框图。

图2为本发明实施例中的决策树的示意图。

具体实施方式

本发明提出的一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，下面将结合实施例进行详细的说明。

本发明针对具有不确定参数的单机调度问题，性能指标选取为总拖期。假设所有工件均在一台机器上进行加工，每个工件仅需要加工一次，加工过程不能中断，所有工件均在加工开始的时刻释放。每个工件具有加工时间和交货时间。模型的目标就是寻找一个最优的加工序列，使得总拖期最短。

在本发明提出的分布式鲁棒优化模型中，设不确定参数是互相独立的随机变量，每一个随机变量的分布未知，但属于一个均值向量和协方差矩阵所确定的分布集中；由于有不确定参数的存在，所有工件的拖期之和也是一个随机变量，因此系统性能选取为最坏情况下总拖期的期望。对于一个调度方案而言，最坏情况的定义即在所有不确定参数的分布中，使得总拖期在该调度方案下期望最大的参数分布情况；在此种设定下，分布式鲁棒优化模型的目标为寻找一个最优的鲁棒调度序列，使得该序列的最坏情况下的总拖期具有最小的期望值。

本发明方法考虑一个工件加工时间为不确定因素的单机调度问题，本发明提出的一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，包括以下步骤：

1)建立针对单机调度问题的分布式鲁棒优化模型；

1-1)确定分布式鲁棒优化模型的不确定参数；

根据不确定性环境下生产过程的总体要求以及不确定性环境下的随机因素，确定分布式鲁棒优化模型中生产调度的各不确定参数，并获取各不确定参数的一阶矩和二阶矩信息，从而构建不确定参数的分布集；

设不确定参数为工件的加工时间，各个工件的加工时间是独立的。假设有N个工件构成工件集合为N＝{1,2,...,N}，根据历史数据或者人工经验，得到工件加工时间的一阶矩和二阶矩信息，因此所有工件的加工时间构成一个随机向量p＝{p₁,...,p_N}，其中p₁,...,p_N分别表示工件1,...,N的加工时间；该随机向量的分布P^p未知，但是属于一个均值向量和协方差矩阵确定的分布集中，该分布集的定义如式(1)表示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (1)

其中P^p为任意一个符合均值向量为μ＝{μ₁,...,μ_N}、协方差矩阵为Σ＝diag{σ₁ ²,...,σ_N ²}的分布，μ₁,...,μ_N表示工件1,...,N的加工时间的平均值，σ₁ ²,...,σ_N ²表示表示工件1,...,N的加工时间的方差，diag{·}表示对角阵；E(.)和Cov(.)分别表示加工时间的均值向量和协方差矩阵；

1-2)确定模型目标函数与约束条件；

1-2-1)确定模型的目标函数：

本发明考虑的是最小化总拖期问题，采用最坏分布情况下的总拖期之和的期望值作为目标函数值。在给定一个调度方案的情况下，所有工件的拖期之和表示为如式(2)所示：

其中y＝{y_ij,i,j＝1,...N}为模型的决策变量，一个向量y对应一个可行的调度方案，如果工件i在工件j之前加工，则y_ij＝1，否则y_ij＝0；根据y的定义，工件j的完工时间可以表示为

d_j为工件j的交货期，因此每个工件的拖期即表示为

考虑最坏分布情况下的总拖期之和的期望值，表示为如式(3)所示：

sup_p～(μ,Σ)E[f(p,y)] (3)

其中sup表示求取集合的上确界，p～(μ,Σ)表示所有工件加工时间向量p属于均值向量为μ和协方差矩阵为Σ的分布集，E表示期望；

对于一个调度方案y而言，最坏的情况即如式(4)定义；该模型的目标就是通过求解获得一个最优的调度方案y^*，使得在该调度方案下，最坏分布情况下的总拖期之和的期望值最小，综上所述，模型的目标函数表达式如下：

y^*＝argmin_ysup_p～(μ,Σ)E[f(p,y)] (4)

1-2-2)确定模型的约束条件；

在构建完目标函数后，还需要考虑模型中的约束条件，本方法结合实际生产情况，具体考虑以下三种约束，具体为：

1-2-2-1)随机加工时间约束；

所有工件的加工时间p的分布未知，但属于一个均值向量、协方差矩阵已知的分布集中，表达式如式(5)所示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (5)

1-2-2-2)可行加工序列位置约束；

两个工件之间有确定的先后顺序(任意多个工件之间不能在先后顺序上出现不合理的情况，如先后顺序成环，举例来说：工件i在工件j之前，j在k之前，而k又在i之前)，表达式如式(6)和(7)所示：

y_ij+y_ji＝1,i,j＝1,...,N,i≠j (6)

y_ij+y_jk+y_ki≤2,i,j,k＝1,...,N,i≠j,j≠k,k≠i (7)

式(6)和(7)保证了任意多个工件先后顺序均不会出现成环的情况。

1-2-2-3)可行调度方案约束；

任意可行调度方案y中的每一个元素均为0-1变量，表达式如下：

y_ij∈{0,1},i≠j (8)

如式(6)-(8)表示的约束条件均为调度方案的可行性约束，共同构成了调度方案的可行域，表达式如式(9)所示：

1-3)建立基于分布式鲁棒模型的最小化总拖期的单机调度的数学表达式；

结合1-1)和1-2)所述，建立完整的基于分布式鲁棒模型的最小化总拖期的单机调度的数学表达式，如下形式：

其中，式(10)为分布式鲁棒优化模型的目标函数，式(11)即为步骤1-2)中所述的约束条件；模型目标是在所有的可行调度方案中，寻找在最坏情况下总拖期期望最小的调度方案，即确定决策变量y的值；

2)对分布式鲁棒优化模型进行转化；

将步骤1)建立的不确定的生产调度的分布式鲁棒优化模型转化为一个确定的二阶锥规划模型，具体步骤如下：

2-1)根据步骤1)得到的分布式鲁棒优化模型的目标函数确定分布式鲁棒优化模型的上界，表达式如下：

其中，根据步骤1-2-1)中工件完工时间的定义可知，

为工件的延迟时间；

如式(12)所示的不等式将步骤1)中分布式鲁棒优化模型的目标函数转化为了模型的上界，从而将对分布式鲁棒优化模型的求解转化为对上界的求解，

2-2)根据步骤2-1)的结果，将步骤1)建立的分布式鲁棒优化模型转化为一个整数二阶锥规划模型，整数二阶锥规划模型的目标函数的表达式如下：

3)使用分支定界算法对步骤2)转化得到的整数二阶锥规划模型进行求解，得到最小化总拖期的单机调度方案；

本发明在模型求解方面，提出了一种针对该生产调度问题分布式鲁棒优化模型的精确求解算法，本发明使用分支定界的方法进行搜索，采用模型分析过程得到的问题性质作为支配准则；提出了一种基于序列的上下界估计方法，提出了基于两种不同建模方式的拉格朗日松弛的上下界估计方法。最终的分支定界算法以拉格朗日松弛方法得到的初始解为上界，搜索过程使用支配准则进行剪枝，采用基于序列的下界估计方法对未进行排序的工件集合进行下界的估计，从而在较短时间内对模型进行求解，获得精确的最优解。具体步骤如下：

3-1)构建搜索树；

将所有工件加工序列的排列组合通过枚举的方式在一个搜索树中进行表示。搜索树中，除根结点外每一个结点代表一个工件，每一个工件的所有可能的下一个工件构成了分支，从树的根结点到任意一个叶子结点所对应的一条路径即为所有工件加工顺序的一种组合，对应一个生产调度的可行解；搜索树中除根结点的前l层(l＝1,...,N)代表了包含l个工件的子序列。

本发明中一个搜索树的示例如图1所示。图1表示了一个三个工件排序的搜索树的结构示意图，工件的编号分别为1，2，3，其中圆内数字表示工件号，圆外数字表示结点编号。树的根结点(Root)不具有特殊含义，只起到连接作用，除了根结点之外，每一个结点代表了一个工件。每一个结点之后可以列出的结点称为该结点的分支，如根结点0有三个分支分别为0结点到1结点，0结点到2结点，0结点到3结点。树的根结点位于第0层，之后每多一个分支就增加一层，最后一层的结点称为叶子结点，每一个从根结点到叶子结点的路径对应了一种所有工件的加工顺序的组合，即一个可行的生产调度方案。从1号结点往下搜索可以得到从根结点到10号和到11号叶子结点两个可行解，即工件序列分别为1，2，3和1，3，2的两个生产调度方案，这两个可行解是1号结点所对应的两个可行解子集，如果在搜索到1号结点时，通过上下界估计方法判断出1号结点所对应的可行解子集不如当前找到的最优解，则不再从1号结点继续往下搜索，这个过程称为剪枝；如果通过支配准则得到1号不能排在2号工件之前，则所有分支中，1号工件在2号工件之前的路径全部可以不再进行搜索，从而减小了搜索的规模。

3-2)通过采用基于拉格朗日松弛的上界初始值估计方法得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界(该上界的作用在于，在对搜索树搜索过程中，可以对每一个分支的下界进行估计，当该分支的下界超出上界时，该分支便不再进行搜索，同时搜索过程中会根据得到的最优解对上界的值进行更新)；

本发明提出了两种基于拉格朗日松弛的上界初始值估计方法，分别为基于2-index模型的拉格朗日松弛的上界初始值估计方法和基于3-index模型的拉格朗日松弛的上界初始值估计方法，在实际运用中选取其中一种即可。

3-2-1)基于2-index模型的拉格朗日松弛的上界初始值估计方法；

定义

w_ij＝σ_iy_ij，则式(13)所示整数二阶锥规划模型的目标函数等价为式(14)所示：

同时在原有约束条件式(11)基础上增加新的约束条件如下：

上述两个约束条件即进行了一次变量代换，用w_0j表示

用w_ij表示σ_iy_ij。

在求解过程中，利用拉格朗日松弛方法对约束式(15)和式(16)进行松弛，分别引入拉格朗日乘子λ_0j,j∈N以及λ_ij，i,j∈N，记

其中

记

其中

拉格朗日函数写为：

因此松弛约束式(15)和式(16)之后的混合整数二阶锥规划模型的目标函数可以写为：

式(18)所示目标函数与约束条件式(11)构成了一个新的优化模型，该模型是式(14)与约束条件式(11)、式(15)和式(16)构成的优化模型(等价于式(13))的松弛问题，即松弛了部分约束同时在目标函数中加入了违反约束的惩罚)，该优化模型称为拉格朗日对偶问题。通过优化式(18)与约束条件式(11)所构成的优化模型，得到的式(18)的最优目标函数值是小于式(14)的最优目标函数值的，因此式(18)的最优目标函数值可以作为式(14)的下界。同时式(18)的解是符合约束式(11)的，因此将该解作为步骤2)得到的整数二阶锥规划模型的一个初始可行解(代表一个可行的生产调度方案)，该初始可行解对应的式(13)的目标函数值可以作为搜索上界的初始值。

在给定Λ的情况下，如式(18)所示的目标函数可以分解为两个优化函数，等号右边加号之前的表达式为第一个优化函数，该优化函数采用0-1整数线性规划求解，可以采用CPLEX直接求解得到，加号之后的表达式为第二个优化函数通过下式求解：

由于式(18)对于任意Λ都是有界的，拉格朗日对偶问题是最大化松弛后问题，因此只需要考虑式(19)中的第二种情况，即λ_j≤1,j∈N。因此通过拉格朗日松弛方法最终可以得到如式(20)所示的拉格朗日优化模型：

其中，

如式(20)所示的LRP拉格朗日优化模型的解相对于原问题依旧是可行的，因此通过求解LRP问题得到的调度方案可行解同样可以作为原问题的可行解，将该可行解作为调度方案的初始可行解并在搜索开始时作为当前最优可行解，该可行解对应的式(13)的目标函数值可以作为搜索上界的初始值并在搜索开始时作为当前上界。由于LRP问题为一个混合整数二阶锥规划问题，因此采用约束生成的方法对式(20)进行求解。

约束生成方法是通过利用有限的约束去逼近无限的约束条件的方法，其主要过程是通过一系列的迭代过程，生成一系列的约束从而在一定程度上近似替代问题中无限个约束条件。每一个迭代过程，通过求解一个子问题生成新的约束，在将得到的新的约束加入主问题后，对主问题进行求解从而更新子问题的参数，进行重新求解子问题生成新的约束，重复这一过程直到满足一定的停止条件。

约束生成方法的具体流程如下：

重写LRP为如下形式：

(LRP2)maxγ (22)

假设约束生成方法已经进行到第k次迭代，在第k次迭代过程中，对于固定的Λ^k，求解以下子问题模型：

式(24)所示的子问题模型的最优解是式(20)所示LRP最优解的下界，同时得到的解y^k是原问题式(13)的可行解，因此可以作为原问题的上界。通过求解子问题式(24)，可以生成如式(25)所示的新的约束：

在迭代k次后，主问题中共包含k个如上所示约束，因此得到的主问题如下所示：

(MP_k)maxγ (26)

MP_k的最优解可以作为式(20)所示LRP问题最优解的近似，可以作为原问题的上界。综上所述，基于约束生成的2-index模型的通过拉格朗日松弛算法得到上界初始值的总体流程如下：

3-2-2)基于3-index模型的拉格朗日松弛的上界初始值估计方法；构建3-index模型，对该模型同样采用拉格朗日松弛方法进行上界初始值的估计。

3-index模型的目标函数为：

约束条件包括：

其中，

表示当且仅当工件j排列在第k个位置且紧接着工件i进行加工，否则

C_k，d_k和

分别代表在第k个位置工件的完工时间，交货期以及方差；约束式(29)保证每一个工件只能占据一个位置且只有一个前继工件；约束式(30)和约束式(32)表示第一个位置和最后一个位置的工件应满足的条件；约束式(31)保证所有工件的前后顺序不会出现环；约束式(34)-式(36)是C_k，d_k和

的定义。

对于3-index模型，采用类似的拉格朗日松弛和约束生成方法可以用于估计式(13)所示原问题的上界，具体的步骤如下：

记w_k＝C_k-d_k，

记X为所有可行解

i＝0,...,N,j,k＝1,…,N+1的集合；则3-index模型可以重写为如下形式：

x∈X (43)

松弛约束条件式(38)和式(39)，引入拉格朗日乘子θ_k,k＝1,...,N和

则拉格朗日函数如式(44)表示：

其中，

同样，3-index的拉格朗日对偶问题可以写为：

因此，如3-2-1)中相同的约束生成方法可以用于求解如式(49)所示优化问题，求解得到的最优解x与式(13)所示目标函数中的y是一一对应的，因此x也代表调度方案的一个可行解；将x转化为y得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

具体步骤如下：

本实施例中，采用基于3-index模型的拉格朗日松弛方法得到一个初始可行解，作为搜索过程的上界初始值，以减小搜索空间；

3-3)对搜索树进行预搜索；

枚举搜索树前P层(P的大小依据问题规模而定，该P层不包括根结点)所有可能工件序列组合，将前P层所有可能的工件序列代入目标函数式(13)计算，保留使得目标函数值最小的一种前P层的工件序列；

本实施例中，在算法开始时，先对前5层的所有可能的工件序列进行枚举，找到使得目标函数值最小的一种工件序列，并排除其他所有前5层的工件序列组合；如一共有20个工件，在所有20个工件中任意选取5个工件，如工件1,2,3,4,5，枚举这五个工件的所有排列的可能，包括(1,2,3,4,5)，(1,2,3,5,4)，(1,3,2,4,5)等，选取目标函数值最小的一个排列，则在前五个工件是1,2,3,4,5五个工件的情况下，其他的排列均不可能得到更小的结果，这样可以很快的排除掉前五层较多的可能，从而加快算法的速度。

3-4)确定前P层的工件序列后，对于搜索树第P+1层所有可能的分支进行判定并对搜索树剪枝；具体步骤如下：

3-4-1)采用问题性质作为支配准则对工件进行判定，去掉不满足支配准则的工件，从而对搜索树进行剪枝；

在加工时间为不确定参数的情况下，提出三条可以作为支配准则的问题性质，具体如下：

性质1：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k以及u_j-d_j＜u_k-d_k，则在任意最优的调度方案中，工件j一定在工件k之前加工。

性质2：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k，同时在某个调度方案S中满足|C_max-d_j|＜|C_max-d_k|，其中C_max＝max{C_j(S),C_k(S)}，C_j(S)和C_k(S)表示在调度方案S中工件j和工件k的完工时间，则在任意最优的调度方案中，工件j一定在工件k之前加工。

性质3：如果工件j和工件k属于未调度的工件集合U，同时工件j位于所有未调度工件的最后一位进行加工，则在任意一个调度方案S中，对于所有未进行调度的工件k∈U，都满足：

其中，

其中，B_k(S)表示在调度方案S中在工件k之前加工的工件集合。

采用性质1到性质3作为支配准则，剪掉搜索树中违反支配准则的分支(如在已经搜索的分支中已知工件i不能位于工件j之前，则所有工件i位于工件j之前的所有搜索树分支都不会进行访问)，从而减少搜索空间。

3-4-2)采用下界估计方法对搜索树中分支的目标函数值下界进行估计从而进行剪枝；

在搜索过程中，在确定第P+1层可以选择哪一个工件时，可以先将某一个未加工的工件置于第P+1层，记为结点a，此时记仍未加工的工件集合为R。则可以利用下界估计方法对集合R对应的目标函数值的下界进行估计，当该下界与已经排序完成的工件序列对应的目标函数值之和超出了当前在搜索树中的最优可行解对应的目标函数值时，则选择的工件不能位于当前的P+1层，停止对结点a及其之后的所有分支的搜索，进而去搜索其他结点。本发明中，提出了基于序列的下界估计方法对下界进行估计，具体如下：

式(13)的目标函数可以重写为如下所示：

其中

根据式(52)中RT_j的定义，RT_j(t,Σ²)是分别关于t和Σ²的增函数，因此构造如式(56)所示目标函数：

如果满足E_j(S)≤C_j(S)，

则式(56)的值是原目标函数式(13)的一个下界。定义：

其中μ_[i]和

是工件加工时间的均值和方差按从小到大排列位于第i位的取值，k₁和k₂分别是工件j在两种排列下所在的位置，k是工件j在某一个调度方案S中的位置。按照上式的定义，满足E_j(S)≤C_j(S)，

因此新构造目标函数最优解是原问题的一个下界。对于E_j(S)和

的计算，即根据工件加工时间均值和方差排序，确定每一个工件在某一个位置的目标值下界，之后分配工件的位置，确定工件的排序，因此新问题等价于一个指派问题，可以使用经典的求解指派问题算法(如匈牙利算法等)进行快速求解。

3-4-3)在搜索过程中，利用搜索树的当前上界，对第P+1层选择的工件的可行性进行判定：在确定第P+1层可以选择哪一个工件时，可以先将某一个未排序的工件置于第P+1层，记为结点b，若前P+1层工件序列对应的式(13)的目标函数值大于搜索树当前上界，则该工件不能位于当前的P+1层，停止对结点b及其之后的所有分支的搜索，进而去搜索其他结点；如果前P+1层工件序列对应的式(13)的目标函数值小于搜索树当前上界，则记未加工的工件集合为R'，对工件集合R'可以进行上界的估计，由于基于拉格朗日松弛的上界估计方法速度相对较慢，不适合进行大量使用，因此在实际运用中，对于每一个剩余未排序的工件集合可以按照贪婪方法，从后往前构造可行调度方案，每次选择使得式(13)目标函数值增加最小的工件，将其置于加工序列最后，重复该过程直到所有工件都被加入到序列当中。如果得到的工件的加工序列对应的的上界小于当前上界，则将该上界作为新的当前上界，该加工序列对应的可行解作为新的当前最优可行解。

3-5)对第P+1层所有可能的分支工件按照步骤3-4)遍历完毕后，对于该层所有剩余的分支，对每个分支计算μ_i-d_i,i＝1,...,N；并按照差值从小到大的顺序依次从对应的分支出发，沿搜索树继续下一层工件的搜索；该排序的方法使得本发明有更大可能得到一个更好的可行解，进而减小上界，缩小搜索空间；

3-6)重复步骤3-4)至3-5)，依次对搜索树每一层的分支进行判定并排序；

3-7)当到达搜索树的最后一层时，此时搜索得到一个完整的工件序列作为一个可行解，计算该可行解对应的式(13)的目标函数值并判定：若计算得到的该可行解的目标函数值小于当前上界，则将该可行解作为新的当前最优可行解，并将该可行解对应的目标函数值作为新的当前上界。

3-8)遍历搜索树的每一个分支；遍历完毕后，最后保留的当前最优可行解即为式(13)所示的目标函数的最优可行解，记为y^*，该最优可行解为一个具有最优目标函数值的工件序列，该工件序列为最终得到的最小化拖期的最优单机调度方案。

分支定界方法就是通过一系列的支配准则、上下界方法减少需要搜索的树的分支，从而精确的得到最优解。本发明中采用的分支定界算法的具体算法流程总结如下：

通过本实施例，生成了相关数据分别对提出的支配准则，上下界估计方法等的性能进行了对比分析，同时将提出的算法与商业求解器CPLEX在求解时间和解的质量上进行对比，最后对提出的模型的性能进行了分析。实验结果表明，算法方面，分析得到的问题性质可以大幅提升分支定界算法的速度，本发明提出的基于不同模型的拉格朗日上下界估计方法在得到上下界估计的同时，也提供一个很好的初始解，与最优解的相对误差在5％以内，因此可以作为分支定界算法一个很好的上界。

模型方面，通过分布式鲁棒优化方法得到的调度方案相比于通过确定性模型得到的调度方案具有更好的鲁棒性，即对于不同分布的数据总拖期的性能波动(方差)都远小于确定性模型的调度方法，见附表1。传统的鲁棒优化可能会损失一定的性能，但是本发明中所得到的调度方案在性能上不比确定性模型得到的方案差，多数情况反而更好，由此可以说明所提出的分布式鲁棒优化方法的优势。

表1不同分布下鲁棒调度方案与确定性调度方案的性能对比

以上实例描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (1)

1.一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，其特征在于，包括以下步骤：

1)建立针对单机调度问题的分布式鲁棒优化模型；

1-1)确定分布式鲁棒优化模型的不确定参数；

设不确定参数为工件的加工时间，假设有N个工件构成工件集合为N＝{1,2,...,N}，所有工件的加工时间构成一个随机向量p＝{p₁,...,p_N}，其中p₁,...,p_N分别表示工件1,...,N的加工时间；该随机向量的分布P^p未知，但是属于一个均值向量和协方差矩阵确定的分布集中，该分布集的定义如式(1)所示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (1)

其中P^p为任意一个符合均值向量为μ＝{μ₁,...,μ_N}、协方差矩阵为Σ＝diag{σ₁ ²,...,σ_N ²}的分布，μ₁,...,μ_N表示工件1,...,N的加工时间的平均值，σ₁ ²,...,σ_N ²表示表示工件1,...,N的加工时间的方差，diag{·}表示对角阵；E(.)和Cov(.)分别表示加工时间的均值向量和协方差矩阵；

1-2)确定模型目标函数与约束条件；

1-2-1)确定模型的目标函数；

在给定一个调度方案的情况下，所有工件的拖期之和可以表示为如式(2)所示：

其中y＝{y_ij,i,j＝1,...N}为模型的决策变量，一个向量y对应一个可行的调度方案，如果工件i在工件j之前加工，则y_ij＝1，否则y_ij＝0；工件j的完工时间表示为

d_j为工件j的交货期，因此每个工件的拖期表示为

考虑最坏分布情况下的总拖期之和的期望值，表达式如式(3)所示：

sup_p～(μ,Σ)E[f(p,y)] (3)

其中sup表示求取集合的上确界，p～(μ,Σ)表示所有工件加工时间向量p属于均值向量为μ和协方差矩阵为Σ的分布集，E表示期望；

模型的目标是通过求解获得一个最优的调度方案y^*，使得在该调度方案下，最坏分布情况下的总拖期之和的期望值最小，则模型的目标函数表达式如下：

y^*＝argmin_ysup_p～(μ,Σ)E[f(p,y)] (4)

1-2-2)确定模型的约束条件；

1-2-2-1)随机加工时间约束；

所有工件的加工时间p的分布未知，但属于一个均值向量、协方差矩阵已知的分布集中，表达式如式(5)所示：

P^p＝{P^p|E(p)＝μ,Cov(p)＝Σ} (5)

1-2-2-2)可行加工序列位置约束；

两个工件之间有先后顺序，任意多个工件之间不能在先后顺序上出现不合理的情况，表达式如式(6)和(7)所示：

y_ij+y_ji＝1,i,j＝1,...,N,i≠j (6)

y_ij+y_jk+y_ki≤2,i,j,k＝1,...,N,i≠j,j≠k,k≠i (7)

1-2-2-3)可行调度方案约束；

任意可行调度方案y中的每一个元素均为0-1变量，表达式如下：

y_ij∈{0,1},i≠j (8)

如式(6)-(8)表示的约束条件均为调度方案的可行性约束，构成了调度方案的可行域，表达式如式(9)所示：

1-3)建立基于分布式鲁棒优化模型的最小化总拖期的单机调度的数学表达式；表达式如下：

其中，式(10)为分布式鲁棒优化模型的目标函数，式(11)为步骤1-2)中的约束条件；

2)对分布式鲁棒优化模型进行转化；

将步骤1)建立的分布式鲁棒优化模型转化为整数二阶锥规划模型，具体步骤如下：

2-1)根据步骤1)得到的分布式鲁棒优化模型的目标函数确定分布式鲁棒优化模型的上界，表达式如下：

其中，

为工件的延迟时间，拖期的定义为工件延迟时间和0的较大值，即max{0,L_j}；

如式(12)所示的不等式将步骤1)中分布式鲁棒优化模型的目标函数转化为了模型的上界，从而将对分布式鲁棒优化模型的求解转化为对上界的求解；

2-2)根据步骤2-1)的结果，将步骤1)建立的分布式鲁棒优化模型转化为一个整数二阶锥规划模型；整数二阶锥规划模型的目标函数的表达式如下：

3)对步骤2)转化得到的整数二阶锥规划模型进行求解，使用分支定界算法得到最优单机调度方案；具体步骤如下：

3-1)构建搜索树；

将所有工件加工序列的排列组合通过枚举的方式在一个搜索树中进行表示；搜索树中，除根结点外每一个结点代表一个工件，每一个工件的所有可能的下一个工件构成了分支，从树的根结点到任意一个叶子结点所对应的一条路径即为所有工件加工顺序的一种组合，对应一个生产调度的可行解；

3-2)通过采用基于拉格朗日松弛的上界初始值估计方法得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

基于拉格朗日松弛的上界初始值估计方法采用基于2-index模型的拉格朗日的上界初始值估计方法或基于3-index模型的拉格朗日的上界初始值估计方法；具体如下：

3-2-1)基于2-index模型的拉格朗日松弛的上界初始值估计方法；

定义

w_ij＝σ_iy_ij，则式(13)所示整数二阶锥规划模型的目标函数等价为式(14)所示：

同时在约束条件式(11)基础上增加新的约束条件如下：

利用拉格朗日松弛方法对约束式(15)和式(16)进行松弛，分别引入拉格朗日乘子λ_0j,j∈N以及λ_ij，i,j∈N，记

其中

记

其中

拉格朗日函数写为：

因此，松弛约束式(15)和式(16)之后的混合整数二阶锥规划模型的目标函数写为：

式(18)所示目标函数与约束条件式(11)构成了一个新的优化模型，该优化模型称为拉格朗日对偶问题；

在给定Λ的情况下，如式(18)所示的目标函数分解为两个优化函数：等号右边加号之前的表达式为第一个优化函数，该优化函数求解采用0-1整数线性规划求解，加号之后的表达式为第二个优化函数通过下式求解：

仅考虑式(19)中的第二种情况，即λ_j≤1,j∈N；通过拉格朗日松弛方法最终得到如式(20)所示的拉格朗日优化模型：

其中，

采用约束生成的方法对式(20)进行求解，得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

3-2-2)基于3-index模型的拉格朗日松弛的上界初始值估计方法；

3-index模型的目标函数为：

约束条件包括：

其中，

表示当且仅当工件j排列在第k个位置且紧接着工件i进行加工，否则

C_k，d_k和

分别代表在第k个位置工件的完工时间，交货期以及方差；约束式(29)保证每一个工件只占据一个位置且只有一个前继工件；约束式(30)和约束式(32)表示第一个位置和最后一个位置的工件满足的条件；约束式(31)保证所有工件的前后顺序不出现环；约束式(34)-式(36)是C_k，d_k和

的定义；

记w_k＝C_k-d_k，

记X为所有可行解

的集合；则3-index模型重写为如下形式：

x∈X(43)

松弛约束条件式(38)和式(39)，引入拉格朗日乘子θ_k,k＝1,...,N和

则拉格朗日函数如式(44)表示：

其中，

3-index的拉格朗日对偶问题写为：

对式(49)求解，求解得到的x与式(13)所示目标函数中的y是一一对应的，因此x也代表调度方案的一个可行解；将x转化为y，得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；

3-3)对搜索树进行预搜索；

枚举搜索树除根结点层外前P层的所有工件序列组合，将前P层每种工件序列代入目标函数式(13)计算，保留使得目标函数值最小的一种前P层的工件序列；

3-4)确定前P层的工件序列后，对于搜索树第P+1层所有的分支进行判定并对搜索树剪枝；具体步骤如下：

3-4-1)采用问题性质作为支配准则对工件进行判定，去掉不满足支配准则的工件，从而对搜索树进行剪枝；

在加工时间为不确定参数的情况下，采用以下三条作为工作排序的支配准则，具体如下：

性质1：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k以及u_j-d_j＜u_k-d_k，则在任意最优的调度方案中，工件j在工件k之前加工；

性质2：如果工件j和工件k满足μ_j＜μ_k,σ_j＜σ_k,d_j＜d_k，同时在某个调度方案S中满足|C_max-d_j|＜|C_max-d_k|，其中C_max＝max{C_j(S),C_k(S)}，C_j(S)和C_k(S)表示在调度方案S中工件j和工件k的完工时间，则在任意最优的调度方案中，工件j在工件k之前加工；

性质3：如果工件j和工件k属于未调度的工件集合U，同时工件j位于所有未调度工件的最后一位进行加工，则在任意一个调度方案S中，对于所有未进行调度的工件k∈U，都满足：

其中，

其中，B_k(S)表示在调度方案S中在工件k之前加工的工件集合；

采用性质1到性质3作为支配准则，剪掉搜索树中违反支配准则的分支；

3-4-2)采用基于序列的下界估计方法对搜索树中分支的目标函数值下界进行估计从而进行剪枝；

在搜索过程中，将某一个未加工的工件置于第P+1层，记为结点a，此时记仍未加工的工件集合为R，利用下界估计方法对集合R对应的目标函数值的下界进行估计，当该下界与已经排序完成的工件序列对应的目标函数值之和超出了搜索树当前最优可行解对应的目标函数值时，则该未加工工件不能位于P+1层，停止对结点a及其之后的所有分支的搜索，进而去搜索其他结点；具体如下：

式(13)的目标函数重写为如下所示：

其中

根据式(52)中RT_j的定义，RT_j(t,Σ²)是分别关于t和Σ²的增函数，因此构造如式(56)所示目标函数：

如果满足E_j(S)≤C_j(S)，

则式(56)的值是原目标函数式(13)的一个下界；定义：

其中μ_[i]和

是工件加工时间的均值和方差按从小到大排列位于第i位的取值，k₁和k₂分别是工件j在两种排列下所在的位置，k是工件j在某一个调度方案S中的位置；

3-4-3)利用搜索树的当前上界，对第P+1层所有分支进行判定：

将某一个未排序的工件置于第P+1层，记为结点b，若前P+1层工件序列对应的式(13)的目标函数值大于搜索树当前上界，则该工件不能位于当前的P+1层，停止对结点b及其之后的所有分支的搜索，进而去搜索其他结点；

3-5)对第P+1层所有分支按照步骤3-4)遍历完毕后，对于该层所有剩余的分支，对每个分支计算μ_i-d_i,i＝1,...,N；并按照差值从小到大的顺序依次从对应的分支出发，沿搜索树继续下一层工件的搜索；

3-6)重复步骤3-4)至3-5)，依次对搜索树每一层的分支进行判定并排序；

3-7)当到达搜索树的最后一层时，此时搜索得到一个完整的工件序列作为一个可行解，计算该可行解对应的式(13)的目标函数值并判定：若计算得到的该可行解的目标函数值小于当前上界，则将该可行解作为新的当前最优可行解，并将该可行解对应的目标函数值作为新的当前上界；

3-8)遍历搜索树的每一个分支；遍历完毕后，最后保留的当前最优可行解即为式(13)所示的目标函数的最优可行解，记为y^*，该最优可行解为一个具有最优目标函数值的工件序列，该工件序列为最终得到最优单机调度方案。

Images