本發(fā)明屬于生產(chǎn)調度及運籌學領域，涉及一種生產(chǎn)調度問題的分布集魯棒模型的建模及優(yōu)化求解方法，考慮在工件加工時間不確定的情況下，尋求風險最小的魯棒調度方案。

背景技術：

制造業(yè)在我國國民經(jīng)濟發(fā)展過程中居于非常重要的地位，其發(fā)展狀況直接影響我國的綜合實力。生產(chǎn)調度是制造系統(tǒng)重要的關鍵技術，旨在確保生產(chǎn)制造過程高效有序地進行。對生產(chǎn)過程設計合理的計劃與調度策略，可以有效縮短產(chǎn)品的生產(chǎn)周期，提高準時交單率，改善設備利用率并降低庫存。過去的幾十年中，在確定型生產(chǎn)調度領域已經(jīng)做了相當多的理論研究，在這種確定性的模型中，都假定工件的加工時間以及其他的參數(shù)是確定已知的。然而在現(xiàn)實加工過程中，由于機器或刀具條件、工人加工水平和加工環(huán)境等因素的影響，工件的加工時間往往存在著不確定性。在調度過程中忽略參數(shù)的不確定性，會導致原來的最優(yōu)解并不可行。這使得在不確定環(huán)境下采用確定型模型得到的調度方案，在實際生產(chǎn)中難以達到?jīng)Q策者的事先預期。因此，針對具有不確定性的生產(chǎn)調度問題的研究，逐漸引起了學者們的關注。

首先出現(xiàn)的針對調度問題中不確定性的研究，是隨機調度問題。在隨機調度的模型中，不確定性參數(shù)被看作是一個分布已知的隨機變量，模型的優(yōu)化目標往往是系統(tǒng)長期的性能期望。從80年代開始，學術界和工業(yè)界在不同的方向上展開了隨機調度問題的廣泛研究。以單機調度問題為例，總流經(jīng)時間、最大拖期時間、加權總拖期時間以及總拖期工件個數(shù)等性能指標均已出現(xiàn)在隨機調度模型的優(yōu)化目標中。雖然隨機調度模型在理論上對不確定性調度問題的研究有很好的推進作用，但隨機調度模型的一些固有缺點，限制了其在實際的大規(guī)模調度問題中的應用。這些缺點主要體現(xiàn)在以下幾個方面：1)在隨機模型中，不確定參數(shù)的分布是需要精確知道的。然而在實際的生產(chǎn)環(huán)境中，隨著生產(chǎn)經(jīng)營日益復雜，小批量個性化定制等生產(chǎn)模式的逐漸轉化，很多情況下精確的概率分布很難獲得，只能根據(jù)類似產(chǎn)品的加工時間來估計它的區(qū)間范圍。特別是在新產(chǎn)品，尤其是單件或單次加工的新產(chǎn)品方面，這個問題尤為突出。在這種情況下，隨機調度方法將不再適用。2)在隨機模型中，一般采用某種系統(tǒng)性能的期望作為優(yōu)化目標。這種目標比較適用于制定企業(yè)長期的發(fā)展計劃，而并不適合解決在每次實際運營過程中的利益最大化或風險最小化的問題。3)這類隨機模型通常是NP-難的，一般只能通過啟發(fā)式算法或動態(tài)規(guī)劃算法來進行求解，隨著問題規(guī)模的逐漸擴大，求解隨機調度模型的難度將以指數(shù)形式增長。

由于隨機調度模型具有以上缺點，處理參數(shù)不確定性的另一種方法——魯棒調度模型應運而生。魯棒調度模型最早由Richard L.Daniels等人提出，其中僅知道區(qū)間信息的不確定參數(shù)通過區(qū)間數(shù)據(jù)情景進行刻畫(一個情景代表不確定參數(shù)的一種可能的取值)，這種描述方法相較于隨機模型中對參數(shù)分布函數(shù)的描述，更為簡單且符合實際。自從Richard L.Daniels將魯棒優(yōu)化的思想引入生產(chǎn)調度問題中，近年來在單機調度、并行機調度以及流水車間調度問題中都有相應的發(fā)展和研究。目前的魯棒生產(chǎn)調度問題都是采用基于不確定性集的魯棒優(yōu)化方法，不確定性集為有限的離散集合或者是連續(xù)的區(qū)間形式。在這種模型下，魯棒調度的關鍵問題就在于如何定義最差環(huán)境，求得每個可行解在最差環(huán)境下的魯棒費用(Robust Cost)以及如何在所有可行解的魯棒費用中尋求最優(yōu)。這種基于不確定性集的魯棒調度模型較為符合實際生產(chǎn)的參數(shù)情況，可以尋求到在最差情況下也能有較好系統(tǒng)性能的魯棒決策，以降低決策的風險。但由于只利用了不確定參數(shù)變化范圍的邊界信息，并主要考慮在最差情況下的系統(tǒng)性能，此種基于不確定性集的魯棒調度模型所得到的決策可能過于保守，犧牲了在參數(shù)常態(tài)情況下的系統(tǒng)性能。因此，如何利用歷史數(shù)據(jù)的更多信息，在保證魯棒性的同時降低決策的保守程度是當前魯棒調度中急需解決的問題。

上述的魯棒調度模型均采用基于不確定性集的魯棒優(yōu)化方法，將隨機變量看作屬于某一不確定性集的不確定參數(shù)來進行建模，使得在不確定情況下的決策與確定情況下的決策有相似的計算復雜度，而不像隨機和動態(tài)規(guī)劃一樣被決策變量的維度所限制。但隨著魯棒優(yōu)化理論的發(fā)展，許多研究者試圖將魯棒優(yōu)化與隨機規(guī)劃聯(lián)系起來，汲取兩者所長以達到相互促進的效果。比如在隨機的系統(tǒng)里來評價魯棒最優(yōu)解的性能，以及在不確定概率分布的概念下進行魯棒優(yōu)化的研究等，后者即發(fā)展為分布集魯棒優(yōu)化。

在分布集魯棒優(yōu)化的模型中，參數(shù)的不確定性通過分布函數(shù)集來描述，不確定參數(shù)被看作一個隨機變量，但其分布函數(shù)可以是某個特定分布函數(shù)集中的任意分布。因此在求解模型的魯棒最優(yōu)解時，需要考慮該分布函數(shù)集中所有可能的分布函數(shù)。當參數(shù)不確定性影響解的可行性，即約束中包含不確定參數(shù)時，分布集魯棒優(yōu)化模型需對相應的約束進行泛函映射。例如采用期望函數(shù)進行映射時，魯棒約束集就對應于原約束在所有可能的分布函數(shù)下的期望所構成的約束集。而當不確定性在目標函數(shù)中時，需要選定一種性能測度來評價具有隨機性目標函數(shù)值。這種測度可以選擇期望、方差或者條件風險價值(CVaR)等。若選定了期望函數(shù)作為評價標準，則對給定的決策變量取值，其對應的魯棒目標函數(shù)值即為隨機參數(shù)的分布在分布函數(shù)集中變化時，原始目標函數(shù)期望的最大值。

技術實現(xiàn)要素：

本發(fā)明為了結合隨機調度以及基于不確定性集的魯棒調度問題的優(yōu)點，首次將分布集魯棒優(yōu)化方法引入生產(chǎn)調度問題中，提出一種生產(chǎn)調度問題的分布集魯棒模型的建模及優(yōu)化求解方法。本發(fā)明相較于隨機模型和傳統(tǒng)的魯棒模型，更加符合實際生產(chǎn)的情況，通過利用生產(chǎn)環(huán)境中更多的信息，可以在保證系統(tǒng)性能的情況下，降低決策的風險。

本發(fā)明提出的一種生產(chǎn)調度問題的分布集魯棒模型的建模方法，該方法具體包括以下步驟：

1)構建生產(chǎn)調度問題的分布集魯棒模型DR-SMSP

在DR-SMSP模型中，針對具有隨機加工時間的單機調度問題，系統(tǒng)的性能指標選擇為總流經(jīng)時間TFT，假定所有工件均在加工開始的時刻釋放，即釋放時間均為零；工件的加工時間具有隨機不確定性，隨機加工時間的分布未知，但屬于一個由支撐集、均值向量和協(xié)方差矩陣所確定的分布集中；由于工件的加工時間是隨機向量，所有工件加工的總流經(jīng)時間TFT是一個隨機變量，系統(tǒng)性能TFT的隨機度量選取為具有風險厭惡特性的條件風險價值CVaR；在此種設定下，DR-SMSP模型的目標為尋找一個最優(yōu)的魯棒調度序列，使得該序列的TFT在加工時間服從最差分布的情況下具有最小的CVaR；

1-1)確定模型決策變量；

該模型的決策變量為可行的調度方案，設模型中有n個工件，且工件的集合為J＝{1,2,...,n}，則一個可行的調度方案由矩陣X＝{x_ij,i,j＝1,...,n}表示；其中，如果工件j在加工序列的第i個位置上，則x_ij＝1，反之x_ij＝0；

1-2)加工時間的隨機性表示；

該模型的加工時間為一個隨機向量p，其分布P^p是未知的，但屬于一個由支撐集、均值向量和協(xié)方差矩陣確定的分布集D^p中，該分布集D^p的表達式如式(1)所示：

其中，Sup(p_j)表示每個加工時間的支撐集，E(p)和Cov(p)分別表示加工時間的均值向量和協(xié)方差矩陣；

1-3)構建模型目標函數(shù)；

該模型的系統(tǒng)性能指標為總流經(jīng)時間TFT，在給定一個調度方案X時，TFT由式(2)計算得到：

由于加工時間是隨機向量，所有工件的TFT是一個隨機變量，本模型采用具有風險厭惡特性的條件風險價值CVaR作為隨機TFT的度量；隨機損失Z的CVaR表示其在最差1-α概率下的期望，由式(3)計算得到：

CVaR_α(Z)＝E[Z|Z≥inf{z:Prob(Z＞z)≤1-α}] (3)

其中，α∈(0,1)，表示CVaR的置信水平，Prob表示概率取值，inf表示求取集合中的下確界；

當隨機損失Z的概率分布P^z屬于一個由支撐集、期望以及方差信息確定的分布集D^z時，在D^z中最差情況下的CVaR_α(Z)被定義為魯棒CVaR_α(Z)，即RCVaR_α(Z)，其表達式如式(4)所示：

其中，sup表示取集合中的上確界；則帶有CVaR風險厭惡的生產(chǎn)調度問題的分布集魯棒模型的目標函數(shù)如式(5)所示：

其中，的上標p表明RCVaR所屬的分布集為D^p；

1-4)約束條件；

生產(chǎn)調度問題的分布集魯棒模型包含4類約束條件，其中1類約束是隨機加工時間所服從的分布集約束，另外3類約束是調度方案的可行性約束，具體如下所示：

1-4-1)隨機加工時間約束；

隨機加工時間p的分布未知，但屬于一個由支撐集、均值向量和協(xié)方差矩陣確定的分布集中，表達式如式(6)所示：

1-4-2)工件占用位置約束；

每個工件僅可占用加工序列中的一個位置，如式(7)所示：

1-4-3)可行加工序列位置約束；

可行加工序列中的每個位置僅可被一個工件占用，如式(8)所示：

1-4-4)可行調度方案約束；

可行調度方案X中的每個元素均是0-1變量，如式(9)所示：

x_ij∈{0,1},i＝1,...,n,j＝1,...,n (9)

約束式(7)至式(9)均是調度方案可行性約束，將其整合到一起，形成調度方案的可行域FB，如式(10)所示：

1-5)生產(chǎn)調度問題的分布集魯棒模型的數(shù)學表達式，如式(11)所示：

其中，F(xiàn)B為調度方案的可行域，的上標p表明RCVaR所屬的分布集為D^p，min表示在可行域FB中尋找目標函數(shù)的最小值，arg表示求得最小目標函數(shù)值所對應的最優(yōu)解X^*；

2)對生產(chǎn)調度問題的分布集魯棒模型進行轉化

本模型將式(11)所示的DR-SMSP1模型分解為一個指派問題AP R₁和一個整數(shù)二階錐規(guī)劃問題I-SOCP R₂，給出如式(13)所示的在具有半無限支撐集的分布集中計算隨機變量RCVaR的顯示表達式，具體步驟如下：

2-1)計算隨機變量在特定分布集下的RCVaR；

本模型給出了在具有半無限支撐集的分布集中計算隨機變量的RCVaR的顯示表達式，如式(13)所示；對于隨機變量Z，若其分布函數(shù)屬于分布集D^z：

則其RCVaR_α由式(13)計算得到：

2-2)轉換決策變量；

將決策變量由矩陣X轉換為向量π，轉換關系如式(14)所示：

其中，x_i為矩陣X的第i行向量；根據(jù)此定義形式，π表示工件加工順序的倒序，即給定一個π＝(π(1),π(2),...,π(n))，π(i)＝j表示工件j在第(n-i+1)次序加工；π相應的可行域如式(15)所示：

TFT表示為π與p的內積，如式(16)所示：

f(π,p)＝f(X,p)＝π^Tp (16)

2-3)分布函數(shù)集的映射關系；

由于加工時間p的隨機性，對每一個確定的π來說，f(π,p)是一個隨機變量，記為F_π；基于p的均值向量和協(xié)方差矩陣，F(xiàn)_π的均值μ_f(π)和方差分別為：

μ_f(π)＝π^Tμ；

進而令F_π的分布集為：

對于任意一個D^p中的分布P^p，若隨機向量p～P^p，則其相應投影隨機變量的支撐集為[0,∞)，均值為π^Tμ＝μ_f(π)，方差為即的分布在分布集D^f中；因此，在分布集D^f中求得的是在D^p中求得的的一個上界，即：

進而，DR-SMSP1中的可以用其上界來替代，轉化成DR-SMSP2：

根據(jù)式(13)所示的RCVaR計算表達式，按式(21)計算得到：

2-4)對生產(chǎn)調度問題的分布集魯棒模型的分解；

DR-SMSP2的最優(yōu)解通過求解一個指派問題AP和一個整數(shù)的二階錐規(guī)劃問題I-SOCP獲得，分解后的模型DR-SMSP3如式(22)所示：

其中，

當不同工件的加工時間不相關時，DR-SMSP3中的兩個子問題：指派子問題AP R₁和整數(shù)二階錐規(guī)劃子問題I-SOCP R₂分別轉化為式(25)及式(26)的形式：

(AP R₁)：

(I-SOCP R₂)：

其中，μ_j與分別為加工時間p_j的均值和方差，

基于本發(fā)明提出的一種生產(chǎn)調度問題的分布集魯棒模型的優(yōu)化求解方法，其特征在于，在加工時間的均值和方差排序一致時，所述模型的優(yōu)化求解采用最短平均加工時間優(yōu)先SAPT的求解準則；在加工時間的均值和方差不具備一致性時，針對二階錐規(guī)劃子問題I-SOCP R₂，采用兩種柯西松弛求解算法，即完全柯西松弛算法CCRA和部分柯西松弛算法PCRA，對模型求解；該方法具體包括以下步驟：

1)最短平均加工時間優(yōu)先SAPT準則；

指派子問題AP R₁通過最短平均加工時間優(yōu)先SAPT準則進行求解，即按照平均加工時間由小到大排序的加工序列即為使得總流經(jīng)時間TFT最小的最優(yōu)序列；而對于二階錐規(guī)劃子問題I-SOCP R₂，其最優(yōu)序列同時受到均值和方差的影響，僅當加工時間的均值和方差排序一致時，子問題AP R₁與I-SOCP R₂的最優(yōu)解相同；

因此，當工件加工時間的均值和方差排序一致時，通過最短平均加工時間優(yōu)先SAPT準則得到的加工序列即為DR-SMSP3的最優(yōu)序列；

當加工時間的均值和方差不具備一致性時，采用兩種柯西松弛求解算法，即完全柯西松弛算法CCRA和部分柯西松弛算法PCRA，對模型求解；

2)采用完全柯西松弛算法CCRA求解；

2-1)松弛目標函數(shù)；

引入非負向量來松弛I-SOCP R₂的均值部分，對于每個分量j，有以下柯西不等式成立：

當且僅當η_j＝μ_j/π_j時，不等式取等；

對于I-SOCP R₂的方差部分，引入一個非負變量t∈R₊將根號去掉，松弛不等式如式(28)所示：

當且僅當時，不等式取等；

令f_A(π)表示I-SOCP R₂原來的目標函數(shù)，則通過式(27)及式(28)兩個不等式，f_A(π)被松弛為f_B(π,η,t)，松弛目標函數(shù)如式(29)所示：

進而，I-SOCP R₂被松弛為如式(30)所示：

上述松弛算法將目標函數(shù)中關于π的部分整合到了一起，使得當η和t固定時，通過排序策略來對π進行尋優(yōu)；由于R₂與之間在最優(yōu)性方面是等價的，對原問題R₂的求解可以轉化為求解當π固定的時候，η和t的最優(yōu)解即為松弛不等式(27)與(28)的取等條件；當η和t固定時，最優(yōu)的π可以通過ξ(η,t)的非增排序獲得，其中ξ(η,t)為f_B(π,η,t)中π的系數(shù)向量，如式(31)所示：

2-2)采用完全柯西松弛算法求解具體流程如下：

在CCRA算法開始前，給定工件加工時間的均值向量μ以及方差向量σ，并對置信水平α、初始值個數(shù)L以及最大迭代次數(shù)Iter_max設置相應的數(shù)值；在設置了初始值個數(shù)L之后，初始值的變化步長為θ_s＝1/(L-1)，第l個初始值由θ_lμ+(1-θ_l)σ²的非增排序得到，其中θ_l＝(l-1)×θ_s；

CCRA算法分別針對L個初始值進行迭代求解，然后在求得的L組解中選取最優(yōu)的一組作為算法最終輸出的最優(yōu)解；在針對每個初始值的迭代過程中，對于第m次迭代，首先在固定的π_m-1下，求解min_η,tf_B(π_m-1,η,t)，最優(yōu)解如式(32)所示：

然后在固定的η_m和t_m下，求解min_πf_B(π,η_m,t_m)，最優(yōu)解π_m由ξ(η_m,t_m)的非增排序獲得；如果π_m＝π_m-1或者m達到了最大迭代次數(shù)Iter_max，則令π^l＝π_m，記為第l個初始值對應的解；當L個初始值所對應的解都求得之后，令作為CCRA算法的估計最優(yōu)值，而π^*(α)＝π^l為相應的最優(yōu)解；

3)采用部分柯西松弛算法PCRA求解；

3-1)松弛目標函數(shù)；

在PCRA中，原問題I-SOCP R₂的目標函數(shù)f_A(π)被松弛為f_C(π,t)，如式(33)所示：

且f_A(π)≤f_C(π,t)的取等條件為

進而，原問題I-SOCP R₂被松弛為如式(34)所示：

與R₂在最優(yōu)的情況下是等價的，即若(π^*,t^*)是的最優(yōu)解，則π^*也是的最優(yōu)解，因此對原問題I-SOCP R₂的求解可以轉化為求解

3-2)松弛后的子問題求解；

由于x_ij∈{0,1}，且可以通過式(35)計算：

因此，f_C(π,t)可以重新寫成關于X和t的函數(shù)，如式(36)所示：

令則f_C(X,t)被簡化成：

當t固定時，是一個標準的指派問題，通過匈牙利Hungarian算法精確求解；

3-3)采用部分柯西松弛算法具體流程如下：

在PCRA算法開始前，給定工件加工時間的均值向量μ以及方差向量σ，并對置信水平α以及最大迭代次數(shù)Iter_max設置相應的數(shù)值；算法的初始值π₀通過μ的非增排序獲得；在迭代過程中，對于第m次迭代，首先在固定的π_m-1下，求解min_tf_C(π_m-1,t)，最優(yōu)解為：然后在固定的t_m下，求解min_Xf_C(X,t_m)，最優(yōu)解X_m由匈牙利算法求得，再經(jīng)由式(38)求得π_m：

如果π_m＝π_m-1或者m達到了最大迭代次數(shù)Iter_max，則令為PCRA算法的估計最優(yōu)值，X^*(α)＝X_m,π^*(α)＝π_m為相應的最優(yōu)解。

本發(fā)明的特點及有益效果在于：

為了能夠結合隨機調度以及基于不確定性集的魯棒調度問題的優(yōu)點，本發(fā)明采用基于不確定分布函數(shù)集的分布集魯棒優(yōu)化方法對生產(chǎn)調度問題進行建模。在分布集魯棒優(yōu)化中，不確定參數(shù)用隨機變量來表示，但是該隨機變量的分布函數(shù)未知，且屬于某個特定的分布函數(shù)集合。在優(yōu)化的過程中，需要考慮該分布函數(shù)集合中所有可能的分布函數(shù)。雖然在分布集魯棒優(yōu)化模型中，不確定參數(shù)仍被看作一個隨機變量，但相較于隨機優(yōu)化模型，這里并不需要明確分布函數(shù)的具體形式，僅需確定一個分布函數(shù)的集合即可。而相較于基于不確定性集魯棒優(yōu)化模型，分布集魯棒優(yōu)化不僅利用了參數(shù)變化范圍的信息，還將其均值和方差等更多的信息考慮進來，以降低決策的保守程度。因此，將分布集魯棒優(yōu)化方法應用于生產(chǎn)調度問題中，將比已有的魯棒建模方法更加符合實際生產(chǎn)的情況，通過利用生產(chǎn)環(huán)境中更多的信息，在保證系統(tǒng)性能的情況下，降低決策的風險。

1)分布集魯棒模型僅需利用隨機向量的支撐集、一階矩和二階矩信息，并不需要精確的知道其分布信息，這一點相較于隨機模型更符合實際生產(chǎn)的情況，實用性更強。

2)分布集魯棒模型利用了不確定參數(shù)的一階矩和二階矩信息，相較于僅用區(qū)間變化范圍的傳統(tǒng)魯棒調度模型，具有更小的保守性。通過利用更多的信息，使求得的最優(yōu)魯棒解在保證魯棒性的前提下，有更好的系統(tǒng)性能。

3)分布集魯棒生產(chǎn)調度模型考慮了決策者的風險厭惡特性，可以通過較小的系統(tǒng)平均性能的損失，極大的降低決策者所承擔的風險。而且可以通過置信水平的不同設置，來平衡系統(tǒng)性能和魯棒性之間的關系。使得決策者可以根據(jù)當前的需求，設置相應的參數(shù)值，以得到最合適的調度策略。

4)兩種柯西松弛算法都是多項式時間的迭代下降算法，其中CCRA算法每次迭代的計算復雜度為O(nlogn)，PCRA算法每次迭代的計算復雜度為O(n³)。通過與商業(yè)求解器CPLEX進行對比，實驗結果表明這兩種算法均具有很高的精度以及效率。在工件數(shù)達到200個時，CCRA算法和PCRA算法分別能在0.1秒和70秒左右求得精度很高的解，非常適用于求解大規(guī)模問題。

具體實施方式

本發(fā)明提出了一種生產(chǎn)調度問題的分布集魯棒模型的建模及優(yōu)化求解方法，下面結合具體實施例進一步詳細說明如下。

本發(fā)明提出的一種生產(chǎn)調度問題的分布集魯棒模型的建模及優(yōu)化求解方法，其具體實現(xiàn)方式包括以下步驟：

1)構建生產(chǎn)調度問題的分布集魯棒模型DR-SMSP

本發(fā)明所關注的是具有隨機加工時間的單機調度問題，針對該問題建立了分布集魯棒優(yōu)化模型(DR-SMSP)。在單機調度問題中，所有工件均在一臺機器上加工，每個工件均有其特定的加工時間，而且在加工的過程中不能被中斷。求解調度問題的目的即為尋找一個工件加工順序，使得總流經(jīng)時間、最大拖期時間或者拖期工件個數(shù)等某個系統(tǒng)性能指標達到最優(yōu)。

在DR-SMSP模型中，系統(tǒng)的性能指標選擇為總流經(jīng)時間(Total Flow Time,TFT)，假定所有工件均在加工開始的時刻釋放，即釋放時間均為零。工件的加工時間具有隨機不確定性，隨機加工時間的分布未知，但屬于一個由支撐集、均值向量和協(xié)方差矩陣所確定的分布集中。由于工件的加工時間是隨機向量，所有工件加工的總流經(jīng)時間TFT是一個隨機變量，為了綜合考慮期望和方差的影響，系統(tǒng)性能TFT的隨機度量選取為具有風險厭惡特性的條件風險價值(CVaR，Conditional Value at Risk)。在此種設定下，DR-SMSP模型的目標為尋找一個最優(yōu)的魯棒調度序列，使得該序列的TFT在加工時間服從最差分布的情況下具有最小的CVaR。

1-1)確定模型決策變量

該模型的決策變量為可行的調度方案，設模型中有n個工件，且工件的集合為J＝{1,2,…,n}，則一個可行的調度方案可由矩陣X＝{x_ij,i,j＝1,…,n}表示。其中，如果工件j在加工序列的第i個位置上，則x_ij＝1，反之x_ij＝0。

1-2)加工時間的隨機性表示

在考慮參數(shù)不確定性的魯棒單機調度問題中，加工時間的不確定性大多是通過區(qū)間型的不確定性集來表示的，其中通常會采用一個預算參數(shù)來控制結果的保守程度。在本發(fā)明的DR-SMSP模型中，我們將不確定加工時間的更多信息考慮進來。把加工時間看成一個隨機向量p，其分布P^p是未知的，但屬于一個由支撐集、均值向量和協(xié)方差矩陣確定的分布集D^p中，該分布集D^p的表達式如式(1)所示：

其中，Sup(p_j)表示每個加工時間的支撐集，E(p)和Cov(p)分別表示加工時間的均值向量和協(xié)方差矩陣；

1-3)構建模型目標函數(shù)

模型的系統(tǒng)性能指標為總流經(jīng)時間(Total Flow Time,TFT)，在給定一個調度方案X時，TFT由式(2)計算得到：

由于加工時間是隨機向量，所有工件的TFT是一個隨機變量，為了將期望和方差的影響進行綜合的考慮，本發(fā)明采用具有風險厭惡特性的條件風險價值(CVaR，Conditional Value at Risk)作為隨機TFT的度量。隨機損失Z的CVaR_α表示其在最差1-α概率下的期望，由式(3)計算得到：

CVaR_α(Z)＝E[Z|Z≥inf{z:Prob(Z＞z)≤1-α}] (3)

其中，α∈(0,1)，表示CVaR的置信水平，Prob表示概率取值，inf表示求取集合中的下確界。

當隨機損失Z的概率分布屬于一個由支撐集、期望以及方差信息確定的分布集D^z時，在D^z中最差情況下的CVaR_α(Z)被定義為魯棒CVaR_α(Z)，即RCVaR_α(Z)，其表達式如式(4)所示：

其中，sup表示取集合中的上確界。

基于以上討論，帶有CVaR風險厭惡的分布集魯棒單機調度問題的目標函數(shù)如式(5)所示：

其中，的上標p表明RCVaR所屬的分布集為D^p。

1-4)約束條件

本模型包含4類約束條件，其中1類約束是隨機加工時間所服從的分布集約束，另外3類約束是調度方案的可行性約束，具體如下所示：

1-4-1)隨機加工時間約束；隨機加工時間p的分布未知，但屬于一個由支撐集、均值向量和協(xié)方差矩陣確定的分布集中，表達式如式(6)所示：

1-4-2)工件占用位置約束；

每個工件僅可占用加工序列中的一個位置，如式(7)所示：

1-4-3)可行加工序列位置約束；

可行加工序列中的每個位置僅可被一個工件占用，如式(8)所示：

1-4-4)可行調度方案約束；

可行調度方案X中的每個元素均是0-1變量，如式(9)所示：

x_ij∈{0,1},i＝1,...,n,j＝1,...,n (9)

約束式(7)至式(9)均是調度方案可行性約束，將其整合到一起，形成調度方案的可行域FB，如式(10)所示：

1-5)單機調度問題分布集魯棒模型的數(shù)學表達式，如式(11)所示：

其中，F(xiàn)B為調度方案的可行域，的上標p表明RCVaR所屬的分布集為D^p，min表示在可行域FB中尋找目標函數(shù)的最小值，arg表示求得最小目標函數(shù)值所對應的最優(yōu)解X^*。

2)對生產(chǎn)調度問題的分布集魯棒模型進行轉化

本模型將式(11)所示的DR-SMSP1模型分解為一個指派問題AP R₁和一個整數(shù)二階錐規(guī)劃問題I-SOCP R₂，給出如式(13)所示的在具有半無限支撐集的分布集中計算隨機變量RCVaR的顯示表達式，具體步驟如下：

2-1)計算隨機變量在特定分布集下的RCVaR

本發(fā)明給出了在具有半無限支撐集的分布集中計算隨機變量的RCVaR的顯示表達式如式(13)所示；對于隨機變量Z，若其分布函數(shù)屬于分布集D^z：

則其RCVaR_α由式(13)計算得到：

2-2)轉換決策變量

為了方便模型的表示與計算，本發(fā)明將決策變量由矩陣X轉換為向量π，轉換關系如式(14)所示：

其中x_i為矩陣X的第i行向量。根據(jù)此定義形式，π表示工件加工順序的倒序，即給定一個π＝(π(1),π(2),...,π(n))，π(i)＝j表示工件j在第(n-i+1)次序加工。π相應的可行域如式(15)所示：

進而，TFT表示為π與p的內積，如式(16)所示：

(TFT)f(π,p)＝f(X,p)＝π^Tp (16)

2-3)分布函數(shù)集的映射關系

由于加工時間p的隨機性，對每一個確定的π來說，f(π,p)是一個隨機變量，記為F_π?；趐的均值向量和協(xié)方差矩陣，F(xiàn)_π的均值μ_f(π)和方差分別為：

μ_f(π)＝π^Tμ；

進而令F_π的分布集為：

對于任意一個D^p中的分布P^p，若隨機向量p～P^p，則其相應投影隨機變量的支撐集為[0,∞)，均值為π^Tμ＝μ_f(π)，方差為即的分布在分布集D^f中。因此，在分布集D^f中求得的是在D^p中求得的的一個上界，即：

進而，DR-SMSP1中的可以用其上界來替代，轉化成DR-SMSP2：

根據(jù)式(13)所示的RCVaR計算表達式，按式(21)計算得到：

2-4)對生產(chǎn)調度問題的分布集魯棒模型的分解

由于式(20)所示的DR-SMSP2中包含分段約束，使得問題較難求解，本發(fā)明將其分解成兩個獨立且相對簡單的錐規(guī)劃問題，以避免處理分段的約束。通過嚴格的理論分析和證明，DR-SMSP2的最優(yōu)解通過求解一個指派問題(AP)和一個整數(shù)的二階錐規(guī)劃問題(I-SOCP)獲得，分解后的模型如式(22)所示：

其中，

如式(22)-式(24)所示，模型已被分解為一個指派子問題(AP R₁)和一個整數(shù)二階錐規(guī)劃子問題(I-SOCP R₂)。在通常情況下，DR-SMSP3模型可以用CPLEX等商業(yè)求解器來求解，然而在問題規(guī)模較大時，求解整數(shù)二階錐規(guī)劃子問題(I-SOCP R₂)需要耗費巨大的計算時間。因此，本發(fā)明通過分析問題的性質，給出了在不同工件的加工時間不相關的情況下，更加高效的求解策略(最短平均加工時間優(yōu)先準則，SAPT)以及算法(完全柯西松弛算法，CCRA；部分柯西松弛算法，PCRA)。在此情況下，DR-SMSP3中的兩個子問題分別轉化為式(25)及式(26)的形式：

(AP R₁)：

(I-SOCP R₂)：

其中，μ_j與分別為加工時間p_j的均值和方差，

一種基于本發(fā)明提出的生產(chǎn)調度問題的分布集魯棒模型的優(yōu)化求解方法，在加工時間的均值和方差排序一致時，本發(fā)明給出了最短平均加工時間優(yōu)先(SAPT)的求解準則；在加工時間的均值和方差不具備一致性時，針對二階錐規(guī)劃子問題(I-SOCP R₂)，基于柯西不等式松弛以及交替優(yōu)化的思路，本發(fā)明設計了兩種更加通用的柯西松弛求解算法，CCRA算法與PCRA算法各有優(yōu)點，可以根據(jù)實際需要選擇使用。具體步驟如下：

1)最短平均加工時間優(yōu)先(SAPT)準則

以總流經(jīng)時間(TFT)為優(yōu)化目標的確定型單機調度問題可以通過最短加工時間優(yōu)先(SPT)準則精確求解，即按照加工時間由小到大排序的加工序列即為使得TFT最小的最優(yōu)序列。通過與確定型問題的對比，最短加工時間優(yōu)先(SPT)準則可以直接引申為最短平均加工時間優(yōu)先(SAPT)準則，用于求解指派子問題(AP R₁)。而對于二階錐規(guī)劃子問題(I-SOCP R₂)，其最優(yōu)序列同時受到均值和方差的影響。僅當加工時間的均值和方差排序一致，即較小的均值也擁有較小的方差時，子問題AP R₁與I-SOCP R₂的最優(yōu)解相同，進而直接為DR-SMSP3的最優(yōu)解。以上分析給出了在均值和方差一致的情況下，DR-SMSP3的求解策略：

當工件加工時間的均值和方差排序一致時，通過最短平均加工時間優(yōu)先(SAPT)準則得到的加工序列即為DR-SMSP3的最優(yōu)序列。

2)采用完全柯西松弛算法(CCRA)求解

當加工時間的均值和方差不具備一致性時，SAPT準則將不適用于I-SOCP R₂。因此，本發(fā)明在這種一般的情況下，基于柯西不等式松弛以及交替優(yōu)化的思路，給出了更加通用的求解算法。

2-1)松弛目標函數(shù)

首先引入非負向量來松弛I-SOCP R₂的均值部分，對于每個分量j，有以下柯西不等式成立：

當且僅當η_j＝μ_j/π_j時，不等式取等。

對于I-SOCP R₂的方差部分，通過引入一個非負變量t∈R₊將根號去掉，松弛不等式如式(28)所示：

當且僅當時，不等式取等。

令f_A(π)表示I-SOCP R₂原來的目標函數(shù)，則通過式(27)及式(28)兩個不等式，f_A(π)被松弛為f_B(π,η,t)，松弛目標函數(shù)如式(29)所示：

進而，I-SOCP R₂被松弛為如式(30)所示：

上述松弛將目標函數(shù)中關于π的部分整合到了一起，使得當η和t固定時，可通過排序策略來對π進行尋優(yōu)。由于R₂與之間在最優(yōu)性方面是等價的，對原問題R₂的求解可以轉化為求解當π固定的時候，η和t的最優(yōu)解即為松弛不等式(27)與(28)的取等條件；當η和t固定時，最優(yōu)的π可以通過ξ(η,t)的非增排序獲得，其中ξ(η,t)為f_B(π,η,t)中π的系數(shù)向量，如式(31)所示：

2-2)采用完全柯西松弛算法求解具體流程如下：

根據(jù)步驟2-1)部分的目標松弛與特征分析，本發(fā)明設計了一個迭代下降算法來求解并將其命名為完全柯西松弛算法(CCRA)。

在CCRA算法開始前，需給定工件加工時間的均值向量μ以及方差向量σ，并對置信水平α、初始值個數(shù)L以及最大迭代次數(shù)Iter_max設置相應的數(shù)值。其中均值向量μ以及方差向量σ根據(jù)實際加工時間的歷史數(shù)據(jù)進行選取，置信水平α的默認值為0.95，可根據(jù)決策者對風險的厭惡程度進行調整，初始值個數(shù)L根據(jù)精度要求選取，L越大則算法精度越高，最大迭代次數(shù)Iter_max則設置為工件個數(shù)n的30倍。在設置了初始值個數(shù)L之后，初始值的變化步長為θ_s＝1/(L-1)，第l個初始值由θ_lμ+(1-θ_l)σ²的非增排序得到，其中θ_l＝(l-1)×θ_s。

CCRA算法分別針對L個初始值進行迭代求解，然后在求得的L組解中選取最優(yōu)的一組作為算法最終輸出的最優(yōu)解。在針對每個初始值的迭代過程中，對于第m次迭代，首先在固定的π_m-1下，求解min_η,tf_B(π_m-1,η,t)，最優(yōu)解如式(32)所示：

然后在固定的η_m和t_m下，求解min_πf_B(π,η_m,t_m)，最優(yōu)解π_m由ξ(η_m,t_m)的非增排序獲得。如果π_m＝π_m-1或者m達到了最大迭代次數(shù)Iter_max，則令π^l＝π_m，記為第l個初始值對應的解。當L個初始值所對應的解都求得之后，令作為CCRA算法的估計最優(yōu)值，而π^*(α)＝π^l為相應的最優(yōu)解。

可以證明CCRA算法是一個嚴格下降的算法，由于算法中的主要計算開銷發(fā)生在對ξ(η_m,t_m)的排序中，CCRA算法每次迭代的計算復雜度為O(nLlogn)，其中L為初始值的個數(shù)。

3)采用部分柯西松弛算法(PCRA)求解

在CCRA中，關于π的兩部分均被松弛了，使得算法中的min_πf_B(π,η_m,t_m)部分直接通過簡單的排序來處理。但這種完全松弛引入了n+1個輔助變量，而且松弛后的子問題過于簡單，導致算法容易陷入到局部最優(yōu)中。為了解決此問題，本發(fā)明進一步提出了部分柯西松弛算法(PCRA)，其中僅引入一個變量t∈R₊，且松弛后的子問題為一個指派問題。

3-1)松弛目標函數(shù)

在PCRA中，原問題I-SOCP R₂的目標函數(shù)f_A(π)被松弛為f_C(π,t)，如式(33)所示：

且f_A(π)≤f_C(π,t)的取等條件為

進而，原問題I-SOCP R₂被松弛為如式(34)所示：

與R₂在最優(yōu)的情況下是等價的，即若(π^*,t^*)是的最優(yōu)解，則π^*也是的最優(yōu)解，因此對原問題R₂的求解可以轉化為求解

3-2)松弛后的子問題求解

由于x_ij∈{0,1}，且可以通過式(35)計算：

因此，f_C(π,t)可以重新寫成關于X和t的函數(shù)，如式(36)所示：

令則f_C(X,t)可被簡化成：

當t固定時，是一個標準的指派問題，可以通過匈牙利(Hungarian)算法精確求解。

3-3)采用部分柯西松弛算法具體流程如下：

當π固定的時候，t的最優(yōu)解即為松弛不等式(33)的取等條件；當t固定時，子問題是一個如式(37)所示的指派問題。根據(jù)這些特性，本發(fā)明設計了部分柯西松弛算法(PCRA)求解

在PCRA算法開始前，需給定工件加工時間的均值向量μ以及方差向量σ，并對置信水平α以及最大迭代次數(shù)Iter_max設置相應的數(shù)值。其中均值向量μ以及方差向量σ根據(jù)實際加工時間的歷史數(shù)據(jù)進行選取，置信水平α的默認值為0.95，可根據(jù)決策者對風險的厭惡程度進行調整，最大迭代次數(shù)Iter_max則設置為工件個數(shù)n的30倍。算法的初始值π₀通過μ的非增排序獲得。在迭代過程中，對于第m次迭代，首先在固定的π_m-1下，求解min_tf_C(π_m-1,t)，最優(yōu)解為：然后在固定的t_m下，求解min_Xf_C(X,t_m)，最優(yōu)解X_m由匈牙利算法求得，再經(jīng)由式(38)求得π_m：

如果π_m＝π_m-1或者m達到了最大迭代次數(shù)Iter_max，則令為PCRA算法的估計最優(yōu)值，X^*(α)＝X_m,π^*(α)＝π_m為相應的最優(yōu)解。

PCRA算法僅引入了一個變量來松弛原問題，因此這個松弛過程相對于CCRA算法更緊。而且松弛后產(chǎn)生的子問題為指派問題，相較于CCRA算法中的排序策略減少了陷入局部最優(yōu)的可能。這些改進使得采用PCRA算法估計的解與實際最優(yōu)解之間的差距變得更小。然而，PCRA算法中的主要計算負擔為匈牙利算法的實現(xiàn)，其計算復雜度為O(n³)，高于CCRA算法的O(nLlogn)。因此，PCRA算法相較于CCRA算法，是以更多的計算時間來換取更高的求解精度。

本發(fā)明所設計的分布集魯棒生產(chǎn)調度問題的模型及優(yōu)化方法的性能分析如下所示：

(1)兩種柯西松弛算法的性能分析

表1為本發(fā)明所設計的兩種柯西松弛算法CCRA、PCRA與商業(yè)求解器CPLEX的在求解精度與求解效率上的比較。

表1兩種柯西松弛算法CCRA、PCRA與CPLEX的性能比較表

表1中數(shù)據(jù)表明，在問題規(guī)模較小(n小于15)時，CPLEX可以在3小時內求得問題的最優(yōu)解，而PCRA同樣可以得到最優(yōu)解。當問題規(guī)模逐漸增大時，CPLEX在3小時內得不到最優(yōu)解，而CCRA與PCRA算法得到的解均優(yōu)于CPLEX在運行3小時后得到的結果，并且需要的計算時間非常少。因此，通過與CPLEX的比較，說明了本發(fā)明所提出的算法具有很高的精度以及運算效率，并且適用于大規(guī)模問題的求解。

(2)分布集魯棒調度模型的魯棒性分析

本發(fā)明將分布集魯棒調度模型求得的魯棒解與僅考慮均值信息得到的均值解進行對比，來說明分布集魯棒建模方法所獲得的魯棒性。其中，mu、sig、RCVaR分別表示采用魯棒解計算得到的TFT的均值、標準差以及魯棒CVaR；帶有“_S”標記的，為采用均值解得到的結果；帶有“_Re”標記的，為兩種解得到結果的相對差。

表2在5000個加工時間的均值和方差實例下，魯棒解和均值解TFT統(tǒng)計對比表

表2中結果顯示：

1)在相同的置信水平alpha下，魯棒解相較于均值解使得TFT的均值變大，方差變小，RCVaR也變小。即分布集魯棒模型以較小均值性能損失為代價，使得TFT的分散性以及RCVaR降低，以達到減少風險的目的。

2)在均值解中，置信水平alpha對TFT均值和方差沒有影響，在加工時間的均值在[10,50]，標準差在[1,30]這樣的范圍中變化時，mu_S基本穩(wěn)定在1350.29，sig_S基本穩(wěn)定在351.23。alpha只影響RCVaR的求解，隨著alpha越小，RCVaR越小。

3)在魯棒解中，隨著置信水平alpha的減小，TFT的均值在減小，方差在增大；不論在均值，方差，還是RCVaR方面都逐漸接近均值解的結果。即alpha越小，方差的作用越小，兩種結果更加接近，當alpha為0時，魯棒解與均值解是相同的。因此，決策者可以通過設置合適的置信水平參數(shù)來平衡期望的系統(tǒng)性能與承擔的風險。

(3)模型針對不同分布的魯棒性分析

表3為在不同分布下魯棒解得到的TFT結果的統(tǒng)計分析表，說明了由分布集魯棒模型求得的魯棒解在不同分布函數(shù)設定下的魯棒性。首先在特定的加工時間的均值及方差的設定下，得到魯棒序列，由該序列得到的理論上的TFT均值、標準差及RCVaR如表中最后一行所示。進而在不同分布函數(shù)設定下，分別產(chǎn)生500,000個加工時間實例，并計算由魯棒序列得到的實際TFT的均值、標準差及CVaR。

表3在不同分布下魯棒解得到的TFT結果的統(tǒng)計分析表

由表3中的數(shù)據(jù)可以看出，在不同分布下得到的均值與標準差結果均在理論值附近，CVaR的值雖然有所不同，但均小于理論上RCVaR的值。由此可以說明，本發(fā)明所提出的分布集魯棒建模方法所得到的最優(yōu)序列，對加工時間的分布具有很好的魯棒性。