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

摘要

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

主权项

一种基于分布式鲁棒模型的最小化总拖期的单机调度方法，其特征在于，包括以下步骤：1)建立针对单机调度问题的分布式鲁棒优化模型；1‑1)确定分布式鲁棒优化模型的不确定参数；设不确定参数为工件的加工时间，假设有N个工件构成工件集合为N＝{1,2,...,N}，所有工件的加工时间构成一个随机向量p＝{p1,...,pN}，其中p1,...,pN分别表示工件1,...,N的加工时间；该随机向量的分布Pp未知，但是属于一个均值向量和协方差矩阵确定的分布集中，该分布集的定义如式(1)所示：Πp＝{Pp|E(p)＝μ,Cov(p)＝Σ}  (1)其中Pp为任意一个符合均值向量为μ＝{μ1,...,μN}、协方差矩阵为Σ＝diag{σ12,...,σN2}的分布，μ1,...,μN表示工件1,...,N的加工时间的平均值，σ12,...,σN2表示表示工件1,...,N的加工时间的方差，diag{·}表示对角阵；E(.)和Cov(.)分别表示加工时间的均值向量和协方差矩阵；1‑2)确定模型目标函数与约束条件；1‑2‑1)确定模型的目标函数；在给定一个调度方案的情况下，所有工件的拖期之和可以表示为如式(2)所示：f(p,y)=Σj=1Nmax{0,Σi=1,i≠jNyijpi+pj-dj}---(2)]]>其中y＝{yij,i,j＝1,...N}为模型的决策变量，一个向量y对应一个可行的调度方案，如果工件i在工件j之前加工，则yij＝1，否则yij＝0；工件j的完工时间表示为dj为工件j的交货期，因此每个工件的拖期表示为考虑最坏分布情况下的总拖期之和的期望值，表达式如式(3)所示：supp～(μ,Σ)E[f(p,y)]  (3)其中sup表示求取集合的上确界，p～(μ,Σ)表示所有工件加工时间向量p属于均值向量为μ和协方差矩阵为Σ的分布集，E表示期望；模型的目标是通过求解获得一个最优的调度方案y*，使得在该调度方案下，最坏分布情况下的总拖期之和的期望值最小，则模型的目标函数表达式如下：y*＝argminysupp～(μ,Σ)E[f(p,y)]  (4)1‑2‑2)确定模型的约束条件；1‑2‑2‑1)随机加工时间约束；所有工件的加工时间p的分布未知，但属于一个均值向量、协方差矩阵已知的分布集中，表达式如式(5)所示：Πp＝{Pp|E(p)＝μ,Cov(p)＝Σ}  (5)1‑2‑2‑2)可行加工序列位置约束；两个工件之间有先后顺序，任意多个工件之间不能在先后顺序上出现不合理的情况，表达式如式(6)和(7)所示：yij+yji＝1,i,j＝1,...,N,i≠j  (6)yij+yjk+yki≤2,i,j,k＝1,...,N,i≠j,j≠k,k≠i  (7)1‑2‑2‑3)可行调度方案约束；任意可行调度方案y中的每一个元素均为0‑1变量，表达式如下：yij∈{0,1},i≠j  (8)如式(6)‑(8)表示的约束条件均为调度方案的可行性约束，构成了调度方案的可行域，表达式如式(9)所示：Ψ={y|yij+yji=1,i,j=1,...,N,i≠j;yij+yjk+yki≤2,i,j,k=1,...,N,i≠j,j≠k,k≠i;yij∈{0,1},i≠j}---(9)]]>1‑3)建立基于分布式鲁棒优化模型的最小化总拖期的单机调度的数学表达式；表达式如下：minysuppi~(μi,σi2),i=1,...,NE(Σi=1Nmax{0,Σi=1,i≠jNyijpi+pj-dj})---(10)]]>Πp={Pp|E(p)=μ,Cov(p)=Σ}y∈Ψ---(11)]]>其中，式(10)为分布式鲁棒优化模型的目标函数，式(11)为步骤1‑2)中的约束条件；2)对分布式鲁棒优化模型进行转化；将步骤1)建立的分布式鲁棒优化模型转化为整数二阶锥规划模型，具体步骤如下：2‑1)根据步骤1)得到的分布式鲁棒优化模型的目标函数确定分布式鲁棒优化模型的上界，表达式如下：suppi∈(μi,σi2),i=1,...,NE(Σi=1Nmax{0,Σi=1,i≠jNyijpi+pj-dj})=suppi∈(μi,σi2),i=1,...,NΣi=1NE(max{0,Σi=1,i≠jNyijpi+pj-dj})≤Σi=1Nsuppi∈(μi,σi2),i=1,...,NE(max{0,Σi=1,i≠jNyijpi+pj-dj})=Σj=1N12[μ(Lj)+σ2(Lj)+(μ(Lj))2]---(12)]]>其中，为工件的延迟时间，拖期的定义为工件延迟时间和0的较大值，即max{0,Lj}；如式(12)所示的不等式将步骤1)中分布式鲁棒优化模型的目标函数转化为了模型的上界，从而将对分布式鲁棒优化模型的求解转化为对上界的求解；2‑2)根据步骤2‑1)的结果，将步骤1)建立的分布式鲁棒优化模型转化为一个整数二阶锥规划模型；整数二阶锥规划模型的目标函数的表达式如下：minyΣj=1N12[μ(Lj)+σ2(Lj)+(μ(Lj))2]y∈Ψ---(13)]]>3)对步骤2)转化得到的整数二阶锥规划模型进行求解，使用分支定界算法得到最优单机调度方案；具体步骤如下：3‑1)构建搜索树；将所有工件加工序列的排列组合通过枚举的方式在一个搜索树中进行表示；搜索树中，除根结点外每一个结点代表一个工件，每一个工件的所有可能的下一个工件构成了分支，从树的根结点到任意一个叶子结点所对应的一条路径即为所有工件加工顺序的一种组合，对应一个生产调度的可行解；3‑2)通过采用基于拉格朗日松弛的上界初始值估计方法得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；基于拉格朗日松弛的上界初始值估计方法采用基于2‑index模型的拉格朗日的上界初始值估计方法或基于3‑index模型的拉格朗日的上界初始值估计方法；具体如下：3‑2‑1)基于2‑index模型的拉格朗日松弛的上界初始值估计方法；定义wij＝σiyij，则式(13)所示整数二阶锥规划模型的目标函数等价为式(14)所示：miny∈ΨΣj=1N12(w0j+w0j2+Σi=1Nwij2+σj2)---(14)]]>同时在约束条件式(11)基础上增加新的约束条件如下：w0j=Σi=1Nyijpi+pj-dj,∀j∈N---(15)]]>wij=σiyij,∀i,j∈N---(16)]]>利用拉格朗日松弛方法对约束式(15)和式(16)进行松弛，分别引入拉格朗日乘子λ0j,j∈N以及λij，i,j∈N，记其中记其中拉格朗日函数写为：L(y,W,Λ)=Σj=1N[(1+λ0j)(Σi=1Nyijμi+μj-dj)+Σi=1Nλijσiyij+w0j2+Σi=1Nwij+σj2-λ0jw0j-Σi=1Nλijwij]=Σj=1N[(1+λ0j)(Σi=1Nyijμi+μj-dj)+Σi=1Nλijσiyij+wjTwj+σj2-λjwj---(17)]]>因此，松弛约束式(15)和式(16)之后的混合整数二阶锥规划模型的目标函数写为：miny∈Ψ,wj∈RN+1L(y,W,Λ)=miny∈ΨΣj=1N[(1+λ0j)(Σi=1Nyijμi+μj-dj)+Σi=1Nλijσiyij]+minwj∈RN+1,j=1,...,NΣj=1N(wjTwj+σj2-λjwj)---(18)]]>式(18)所示目标函数与约束条件式(11)构成了一个新的优化模型，该优化模型称为拉格朗日对偶问题；在给定Λ的情况下，如式(18)所示的目标函数分解为两个优化函数：等号右边加号之前的表达式为第一个优化函数，该优化函数求解采用0‑1整数线性规划求解，加号之后的表达式为第二个优化函数通过下式求解：仅考虑式(19)中的第二种情况，即λj≤1,j∈N；通过拉格朗日松弛方法最终得到如式(20)所示的拉格朗日优化模型：(LRP)max||λj||≤1,j∈Nminy∈ΨΣj=1N(h(λj,y)+σj1-||λj||2)---(20)]]>其中，h(λj,y)=[(1+λ0j)(Σi=1Nyijμi+μj-dj)+Σi=1Nλijσiyij]---(21)]]>采用约束生成的方法对式(20)进行求解，得到调度方案的初始可行解并在搜索开始时作为当前最优可行解，将该初始可行解对应的式(13)的目标函数值作为搜索上界的初始值并在搜索开始时作为当前上界；3‑2‑2)基于3‑index模型的拉格朗日松弛的上界初始值估计方法；3‑index模型的目标函数为：(3-index)min12Σk=1N(Ck-dk)+(Ck-dk)2+σk2---(28)]]>约束条件包括：Σi=0NΣk=1N+1xi,jk=1,∀j=1,...,N+1---(29)]]>Σj=1Nx0,j1=1---(30)]]>Σi=0Nxi,jk-Σm=1N+1xj,mk+1=0,∀j,k=1,...,N---(31)]]>Σi=1Nxi,N+1N+1=1---(32)]]>xi,jk∈{0,1},∀i=0,...,N,j,k=1,...,N+1---(33)]]>Ck=Σl=1kΣi=0NΣj=1Npjxi,jl,∀k=1,...,N---(34)]]>dk=Σi=0NΣj=1Ndjxi,jk,∀k=1,...,N---(35)]]>σk2=Σl=1kΣi=0NΣj=1Nσj2xi,jl,∀k=1,...,N---(36)]]>其中，表示当且仅当工件j排列在第k个位置且紧接着工件i进行加工，否则Ck，dk和分别代表在第k个位置工件的完工时间，交货期以及方差；约束式(29)保证每一个工件只占据一个位置且只有一个前继工件；约束式(30)和约束式(32)表示第一个位置和最后一个位置的工件满足的条件；约束式(31)保证所有工件的前后顺序不出现环；约束式(34)‑式(36)是Ck，dk和的定义；记wk＝Ck‑dk，记Ξ为所有可行解的集合；则3‑index模型重写为如下形式：min12Σk=1N(Ck-dk)+wk2+Σl=1kΣi=0NΣi=1N(wijl)2---(37)]]>s.t.wk=Ck-dk---(38)]]>wijl=σjxijl---(39)]]>Ck=Σl=1kΣi=0NΣj=1Npjxi,jl,∀k=1,...,N---(40)]]>dk=Σi=0NΣj=1Ndjxi,jk,∀k=1,...,N---(41)]]>σk2=Σl=1kΣi=0NΣj=1Nσj2xi,jl,∀k=1,...,N---(42)]]>x∈Ξ  (43)松弛约束条件式(38)和式(39)，引入拉格朗日乘子θk,k＝1,...,N和则拉格朗日函数如式(44)表示：L(x,W,Θ)=Σk=1N[(1+θk)(Ck-dk)+Σl=1kΣi=0NΣj=1Nθijlσjxijl+wk2+Σl=1kΣi=0NΣj=1Nwijl-θkwk-Σl=1kΣi=0NΣj=1Nθijlwijl]=Σk=1N[(1+θk)(Ck-dk)+Σl=1kΣi=0NΣj=1Nθijlσjxijl+||wk||-θkwk---(44)]]>其中，W={w1,...,wN,w011,...,wN,N+1N+1}---(45)]]>Θ=[θ1,...,θN,θ011,...,θN,N+1N+1]---(46)]]>wk=[wk,wijl],i=0,...,N,j=1,...,N+1,l=1,...,k---(47)]]>θk=[θk,θijl],i=0,...,N,j=1,...,N+1,l=1,...,k---(48)]]>3‑index的拉格朗日对偶问题写为：(LRP3-index)max||θk≤1||,k=1,...,Nminx∈ΞΣk=1N[(1+θk)(Ck-dk)+Σl=1kΣi=0NΣj=1Nθijlσjxijl]---(49)]]>对式(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,dj＜dk以及uj‑dj＜uk‑dk，则在任意最优的调度方案中，工件j在工件k之前加工；性质2：如果工件j和工件k满足μj＜μk,σj＜σk,dj＜dk，同时在某个调度方案S中满足|Cmax‑dj|＜|Cmax‑dk|，其中Cmax＝max{Cj(S),Ck(S)}，Cj(S)和Ck(S)表示在调度方案S中工件j和工件k的完工时间，则在任意最优的调度方案中，工件j在工件k之前加工；性质3：如果工件j和工件k属于未调度的工件集合U，同时工件j位于所有未调度工件的最后一位进行加工，则在任意一个调度方案S中，对于所有未进行调度的工件k∈U，都满足：RTk(S)≤RTk(PU,ΣU2)-ΔRTj---(50)]]>其中，RTk(S)=Ck(S)-dk+(Ck(S)-dk)2+ΣBk(S)2+σk2---(51)]]>RTk(PU,ΣU2)=PU-dk+(PU-dk)2+ΣU2---(52)]]>ΔRTj=RTj(PU,ΣU2)-RTj(PU-μk,ΣU2-σk2)---(53)]]>PU=Σi∈Uμi,ΣU2=Σi∈Uσi2,ΣBk(S)=Σi∈Bk(S)σi2---(54)]]>其中，Bk(S)表示在调度方案S中在工件k之前加工的工件集合；采用性质1到性质3作为支配准则，剪掉搜索树中违反支配准则的分支；3‑4‑2)采用基于序列的下界估计方法对搜索树中分支的目标函数值下界进行估计从而进行剪枝；在搜索过程中，将某一个未加工的工件置于第P+1层，记为结点a，此时记仍未加工的工件集合为R，利用下界估计方法对集合R对应的目标函数值的下界进行估计，当该下界与已经排序完成的工件序列对应的目标函数值之和超出了搜索树当前最优可行解对应的目标函数值时，则该未加工工件不能位于P+1层，停止对结点a及其之后的所有分支的搜索，进而去搜索其他结点；具体如下：式(13)的目标函数重写为如下所示：minS∈Φ12Σi=1NRTj(Cj(S),Σj2(S))---(55)]]>其中根据式(52)中RTj的定义，RTj(t,Σ2)是分别关于t和Σ2的增函数，因此构造如式(56)所示目标函数：minS∈Φ12Σi=1NRTj(Ej(S),Vj2(S))---(56)]]>如果满足Ej(S)≤Cj(S)，则式(56)的值是原目标函数式(13)的一个下界；定义：Ej(S)=Σi=1kμ[i],k1≤kΣi=1k-1μ[i]+μjk1>k---(57)]]>Vj2(S)=Σi=1kσ(i)2,k2≤kΣi=1k-1σ(i)2+σj2k2>k---(58)]]>其中μ[i]和是工件加工时间的均值和方差按从小到大排列位于第i位的取值，k1和k2分别是工件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‑di,i＝1,...,N；并按照差值从小到大的顺序依次从对应的分支出发，沿搜索树继续下一层工件的搜索；3‑6)重复步骤3‑4)至3‑5)，依次对搜索树每一层的分支进行判定并排序；3‑7)当到达搜索树的最后一层时，此时搜索得到一个完整的工件序列作为一个可行解，计算该可行解对应的式(13)的目标函数值并判定：若计算得到的该可行解的目标函数值小于当前上界，则将该可行解作为新的当前最优可行解，并将该可行解对应的目标函数值作为新的当前上界；3‑8)遍历搜索树的每一个分支；遍历完毕后，最后保留的当前最优可行解即为式(13)所示的目标函数的最优可行解，记为y*，该最优可行解为一个具有最优目标函数值的工件序列，该工件序列为最终得到最优单机调度方案。