本发明属于生产调度及运筹学领域，涉及一种生产调度问题的分布集鲁棒模型的建模及优化求解方法，考虑在工件加工时间不确定的情况下，寻求风险最小的鲁棒调度方案。
背景技术：
制造业在我国国民经济发展过程中居于非常重要的地位，其发展状况直接影响我国的综合实力。生产调度是制造系统重要的关键技术，旨在确保生产制造过程高效有序地进行。对生产过程设计合理的计划与调度策略，可以有效缩短产品的生产周期，提高准时交单率，改善设备利用率并降低库存。过去的几十年中，在确定型生产调度领域已经做了相当多的理论研究，在这种确定性的模型中，都假定工件的加工时间以及其他的参数是确定已知的。然而在现实加工过程中，由于机器或刀具条件、工人加工水平和加工环境等因素的影响，工件的加工时间往往存在着不确定性。在调度过程中忽略参数的不确定性，会导致原来的最优解并不可行。这使得在不确定环境下采用确定型模型得到的调度方案，在实际生产中难以达到决策者的事先预期。因此，针对具有不确定性的生产调度问题的研究，逐渐引起了学者们的关注。
首先出现的针对调度问题中不确定性的研究，是随机调度问题。在随机调度的模型中，不确定性参数被看作是一个分布已知的随机变量，模型的优化目标往往是系统长期的性能期望。从80年代开始，学术界和工业界在不同的方向上展开了随机调度问题的广泛研究。以单机调度问题为例，总流经时间、最大拖期时间、加权总拖期时间以及总拖期工件个数等性能指标均已出现在随机调度模型的优化目标中。虽然随机调度模型在理论上对不确定性调度问题的研究有很好的推进作用，但随机调度模型的一些固有缺点，限制了其在实际的大规模调度问题中的应用。这些缺点主要体现在以下几个方面：1)在随机模型中，不确定参数的分布是需要精确知道的。然而在实际的生产环境中，随着生产经营日益复杂，小批量个性化定制等生产模式的逐渐转化，很多情况下精确的概率分布很难获得，只能根据类似产品的加工时间来估计它的区间范围。特别是在新产品，尤其是单件或单次加工的新产品方面，这个问题尤为突出。在这种情况下，随机调度方法将不再适用。2)在随机模型中，一般采用某种系统性能的期望作为优化目标。这种目标比较适用于制定企业长期的发展计划，而并不适合解决在每次实际运营过程中的利益最大化或风险最小化的问题。3)这类随机模型通常是NP-难的，一般只能通过启发式算法或动态规划算法来进行求解，随着问题规模的逐渐扩大，求解随机调度模型的难度将以指数形式增长。
由于随机调度模型具有以上缺点，处理参数不确定性的另一种方法——鲁棒调度模型应运而生。鲁棒调度模型最早由Richard L.Daniels等人提出，其中仅知道区间信息的不确定参数通过区间数据情景进行刻画(一个情景代表不确定参数的一种可能的取值)，这种描述方法相较于随机模型中对参数分布函数的描述，更为简单且符合实际。自从Richard L.Daniels将鲁棒优化的思想引入生产调度问题中，近年来在单机调度、并行机调度以及流水车间调度问题中都有相应的发展和研究。目前的鲁棒生产调度问题都是采用基于不确定性集的鲁棒优化方法，不确定性集为有限的离散集合或者是连续的区间形式。在这种模型下，鲁棒调度的关键问题就在于如何定义最差环境，求得每个可行解在最差环境下的鲁棒费用(Robust Cost)以及如何在所有可行解的鲁棒费用中寻求最优。这种基于不确定性集的鲁棒调度模型较为符合实际生产的参数情况，可以寻求到在最差情况下也能有较好系统性能的鲁棒决策，以降低决策的风险。但由于只利用了不确定参数变化范围的边界信息，并主要考虑在最差情况下的系统性能，此种基于不确定性集的鲁棒调度模型所得到的决策可能过于保守，牺牲了在参数常态情况下的系统性能。因此，如何利用历史数据的更多信息，在保证鲁棒性的同时降低决策的保守程度是当前鲁棒调度中急需解决的问题。
上述的鲁棒调度模型均采用基于不确定性集的鲁棒优化方法，将随机变量看作属于某一不确定性集的不确定参数来进行建模，使得在不确定情况下的决策与确定情况下的决策有相似的计算复杂度，而不像随机和动态规划一样被决策变量的维度所限制。但随着鲁棒优化理论的发展，许多研究者试图将鲁棒优化与随机规划联系起来，汲取两者所长以达到相互促进的效果。比如在随机的系统里来评价鲁棒最优解的性能，以及在不确定概率分布的概念下进行鲁棒优化的研究等，后者即发展为分布集鲁棒优化。
在分布集鲁棒优化的模型中，参数的不确定性通过分布函数集来描述，不确定参数被看作一个随机变量，但其分布函数可以是某个特定分布函数集中的任意分布。因此在求解模型的鲁棒最优解时，需要考虑该分布函数集中所有可能的分布函数。当参数不确定性影响解的可行性，即约束中包含不确定参数时，分布集鲁棒优化模型需对相应的约束进行泛函映射。例如采用期望函数进行映射时，鲁棒约束集就对应于原约束在所有可能的分布函数下的期望所构成的约束集。而当不确定性在目标函数中时，需要选定一种性能测度来评价具有随机性目标函数值。这种测度可以选择期望、方差或者条件风险价值(CVaR)等。若选定了期望函数作为评价标准，则对给定的决策变量取值，其对应的鲁棒目标函数值即为随机参数的分布在分布函数集中变化时，原始目标函数期望的最大值。
技术实现要素：
本发明为了结合随机调度以及基于不确定性集的鲁棒调度问题的优点，首次将分布集鲁棒优化方法引入生产调度问题中，提出一种生产调度问题的分布集鲁棒模型的建模及优化求解方法。本发明相较于随机模型和传统的鲁棒模型，更加符合实际生产的情况，通过利用生产环境中更多的信息，可以在保证系统性能的情况下，降低决策的风险。
本发明提出的一种生产调度问题的分布集鲁棒模型的建模方法，该方法具体包括以下步骤：
1)构建生产调度问题的分布集鲁棒模型DR-SMSP
在DR-SMSP模型中，针对具有随机加工时间的单机调度问题，系统的性能指标选择为总流经时间TFT，假定所有工件均在加工开始的时刻释放，即释放时间均为零；工件的加工时间具有随机不确定性，随机加工时间的分布未知，但属于一个由支撑集、均值向量和协方差矩阵所确定的分布集中；由于工件的加工时间是随机向量，所有工件加工的总流经时间TFT是一个随机变量，系统性能TFT的随机度量选取为具有风险厌恶特性的条件风险价值CVaR；在此种设定下，DR-SMSP模型的目标为寻找一个最优的鲁棒调度序列，使得该序列的TFT在加工时间服从最差分布的情况下具有最小的CVaR；
1-1)确定模型决策变量；
该模型的决策变量为可行的调度方案，设模型中有n个工件，且工件的集合为J＝{1,2,...,n}，则一个可行的调度方案由矩阵X＝{xij,i,j＝1,...,n}表示；其中，如果工件j在加工序列的第i个位置上，则xij＝1，反之xij＝0；
1-2)加工时间的随机性表示；
该模型的加工时间为一个随机向量p，其分布Pp是未知的，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集Dp中，该分布集Dp的表达式如式(1)所示：
其中，Sup(pj)表示每个加工时间的支撑集，E(p)和Cov(p)分别表示加工时间的均值向量和协方差矩阵；
1-3)构建模型目标函数；
该模型的系统性能指标为总流经时间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的概率分布Pz属于一个由支撑集、期望以及方差信息确定的分布集Dz时，在Dz中最差情况下的CVaRα(Z)被定义为鲁棒CVaRα(Z)，即RCVaRα(Z)，其表达式如式(4)所示：
其中，sup表示取集合中的上确界；则带有CVaR风险厌恶的生产调度问题的分布集鲁棒模型的目标函数如式(5)所示：
其中，的上标p表明RCVaR所属的分布集为Dp；
1-4)约束条件；
生产调度问题的分布集鲁棒模型包含4类约束条件，其中1类约束是随机加工时间所服从的分布集约束，另外3类约束是调度方案的可行性约束，具体如下所示：
1-4-1)随机加工时间约束；
随机加工时间p的分布未知，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集中，表达式如式(6)所示：
1-4-2)工件占用位置约束；
每个工件仅可占用加工序列中的一个位置，如式(7)所示：
1-4-3)可行加工序列位置约束；
可行加工序列中的每个位置仅可被一个工件占用，如式(8)所示：
1-4-4)可行调度方案约束；
可行调度方案X中的每个元素均是0-1变量，如式(9)所示：
xij∈{0,1},i＝1,...,n,j＝1,...,n (9)
约束式(7)至式(9)均是调度方案可行性约束，将其整合到一起，形成调度方案的可行域FB，如式(10)所示：
1-5)生产调度问题的分布集鲁棒模型的数学表达式，如式(11)所示：
其中，FB为调度方案的可行域，的上标p表明RCVaR所属的分布集为Dp，min表示在可行域FB中寻找目标函数的最小值，arg表示求得最小目标函数值所对应的最优解X*；
2)对生产调度问题的分布集鲁棒模型进行转化
本模型将式(11)所示的DR-SMSP1模型分解为一个指派问题AP R1和一个整数二阶锥规划问题I-SOCP R2，给出如式(13)所示的在具有半无限支撑集的分布集中计算随机变量RCVaR的显示表达式，具体步骤如下：
2-1)计算随机变量在特定分布集下的RCVaR；
本模型给出了在具有半无限支撑集的分布集中计算随机变量的RCVaR的显示表达式，如式(13)所示；对于随机变量Z，若其分布函数属于分布集Dz：
则其RCVaRα由式(13)计算得到：
2-2)转换决策变量；
将决策变量由矩阵X转换为向量π，转换关系如式(14)所示：
其中，xi为矩阵X的第i行向量；根据此定义形式，π表示工件加工顺序的倒序，即给定一个π＝(π(1),π(2),...,π(n))，π(i)＝j表示工件j在第(n-i+1)次序加工；π相应的可行域如式(15)所示：
TFT表示为π与p的内积，如式(16)所示：
f(π,p)＝f(X,p)＝πTp (16)
2-3)分布函数集的映射关系；
由于加工时间p的随机性，对每一个确定的π来说，f(π,p)是一个随机变量，记为Fπ；基于p的均值向量和协方差矩阵，Fπ的均值μf(π)和方差分别为：
μf(π)＝πTμ；
进而令Fπ的分布集为：
对于任意一个Dp中的分布Pp，若随机向量p～Pp，则其相应投影随机变量的支撑集为[0,∞)，均值为πTμ＝μf(π)，方差为即的分布在分布集Df中；因此，在分布集Df中求得的是在Dp中求得的的一个上界，即：
进而，DR-SMSP1中的可以用其上界来替代，转化成DR-SMSP2：
根据式(13)所示的RCVaR计算表达式，按式(21)计算得到：
2-4)对生产调度问题的分布集鲁棒模型的分解；
DR-SMSP2的最优解通过求解一个指派问题AP和一个整数的二阶锥规划问题I-SOCP获得，分解后的模型DR-SMSP3如式(22)所示：
其中，
当不同工件的加工时间不相关时，DR-SMSP3中的两个子问题：指派子问题AP R1和整数二阶锥规划子问题I-SOCP R2分别转化为式(25)及式(26)的形式：
(AP R1)：
(I-SOCP R2)：
其中，μj与分别为加工时间pj的均值和方差，
基于本发明提出的一种生产调度问题的分布集鲁棒模型的优化求解方法，其特征在于，在加工时间的均值和方差排序一致时，所述模型的优化求解采用最短平均加工时间优先SAPT的求解准则；在加工时间的均值和方差不具备一致性时，针对二阶锥规划子问题I-SOCP R2，采用两种柯西松弛求解算法，即完全柯西松弛算法CCRA和部分柯西松弛算法PCRA，对模型求解；该方法具体包括以下步骤：
1)最短平均加工时间优先SAPT准则；
指派子问题AP R1通过最短平均加工时间优先SAPT准则进行求解，即按照平均加工时间由小到大排序的加工序列即为使得总流经时间TFT最小的最优序列；而对于二阶锥规划子问题I-SOCP R2，其最优序列同时受到均值和方差的影响，仅当加工时间的均值和方差排序一致时，子问题AP R1与I-SOCP R2的最优解相同；
因此，当工件加工时间的均值和方差排序一致时，通过最短平均加工时间优先SAPT准则得到的加工序列即为DR-SMSP3的最优序列；
当加工时间的均值和方差不具备一致性时，采用两种柯西松弛求解算法，即完全柯西松弛算法CCRA和部分柯西松弛算法PCRA，对模型求解；
2)采用完全柯西松弛算法CCRA求解；
2-1)松弛目标函数；
引入非负向量来松弛I-SOCP R2的均值部分，对于每个分量j，有以下柯西不等式成立：
当且仅当ηj＝μj/πj时，不等式取等；
对于I-SOCP R2的方差部分，引入一个非负变量t∈R+将根号去掉，松弛不等式如式(28)所示：
当且仅当时，不等式取等；
令fA(π)表示I-SOCP R2原来的目标函数，则通过式(27)及式(28)两个不等式，fA(π)被松弛为fB(π,η,t)，松弛目标函数如式(29)所示：
进而，I-SOCP R2被松弛为如式(30)所示：
上述松弛算法将目标函数中关于π的部分整合到了一起，使得当η和t固定时，通过排序策略来对π进行寻优；由于R2与之间在最优性方面是等价的，对原问题R2的求解可以转化为求解当π固定的时候，η和t的最优解即为松弛不等式(27)与(28)的取等条件；当η和t固定时，最优的π可以通过ξ(η,t)的非增排序获得，其中ξ(η,t)为fB(π,η,t)中π的系数向量，如式(31)所示：
2-2)采用完全柯西松弛算法求解具体流程如下：
在CCRA算法开始前，给定工件加工时间的均值向量μ以及方差向量σ，并对置信水平α、初始值个数L以及最大迭代次数Itermax设置相应的数值；在设置了初始值个数L之后，初始值的变化步长为θs＝1/(L-1)，第l个初始值由θlμ+(1-θl)σ2的非增排序得到，其中θl＝(l-1)×θs；
CCRA算法分别针对L个初始值进行迭代求解，然后在求得的L组解中选取最优的一组作为算法最终输出的最优解；在针对每个初始值的迭代过程中，对于第m次迭代，首先在固定的πm-1下，求解minη,tfB(πm-1,η,t)，最优解如式(32)所示：
然后在固定的ηm和tm下，求解minπfB(π,ηm,tm)，最优解πm由ξ(ηm,tm)的非增排序获得；如果πm＝πm-1或者m达到了最大迭代次数Itermax，则令πl＝πm，记为第l个初始值对应的解；当L个初始值所对应的解都求得之后，令作为CCRA算法的估计最优值，而π*(α)＝πl为相应的最优解；
3)采用部分柯西松弛算法PCRA求解；
3-1)松弛目标函数；
在PCRA中，原问题I-SOCP R2的目标函数fA(π)被松弛为fC(π,t)，如式(33)所示：
且fA(π)≤fC(π,t)的取等条件为
进而，原问题I-SOCP R2被松弛为如式(34)所示：
与R2在最优的情况下是等价的，即若(π*,t*)是的最优解，则π*也是的最优解，因此对原问题I-SOCP R2的求解可以转化为求解
3-2)松弛后的子问题求解；
由于xij∈{0,1}，且可以通过式(35)计算：
因此，fC(π,t)可以重新写成关于X和t的函数，如式(36)所示：
令则fC(X,t)被简化成：
当t固定时，是一个标准的指派问题，通过匈牙利Hungarian算法精确求解；
3-3)采用部分柯西松弛算法具体流程如下：
在PCRA算法开始前，给定工件加工时间的均值向量μ以及方差向量σ，并对置信水平α以及最大迭代次数Itermax设置相应的数值；算法的初始值π0通过μ的非增排序获得；在迭代过程中，对于第m次迭代，首先在固定的πm-1下，求解mintfC(πm-1,t)，最优解为：然后在固定的tm下，求解minXfC(X,tm)，最优解Xm由匈牙利算法求得，再经由式(38)求得πm：
如果πm＝πm-1或者m达到了最大迭代次数Itermax，则令为PCRA算法的估计最优值，X*(α)＝Xm,π*(α)＝πm为相应的最优解。
本发明的特点及有益效果在于：
为了能够结合随机调度以及基于不确定性集的鲁棒调度问题的优点，本发明采用基于不确定分布函数集的分布集鲁棒优化方法对生产调度问题进行建模。在分布集鲁棒优化中，不确定参数用随机变量来表示，但是该随机变量的分布函数未知，且属于某个特定的分布函数集合。在优化的过程中，需要考虑该分布函数集合中所有可能的分布函数。虽然在分布集鲁棒优化模型中，不确定参数仍被看作一个随机变量，但相较于随机优化模型，这里并不需要明确分布函数的具体形式，仅需确定一个分布函数的集合即可。而相较于基于不确定性集鲁棒优化模型，分布集鲁棒优化不仅利用了参数变化范围的信息，还将其均值和方差等更多的信息考虑进来，以降低决策的保守程度。因此，将分布集鲁棒优化方法应用于生产调度问题中，将比已有的鲁棒建模方法更加符合实际生产的情况，通过利用生产环境中更多的信息，在保证系统性能的情况下，降低决策的风险。
1)分布集鲁棒模型仅需利用随机向量的支撑集、一阶矩和二阶矩信息，并不需要精确的知道其分布信息，这一点相较于随机模型更符合实际生产的情况，实用性更强。
2)分布集鲁棒模型利用了不确定参数的一阶矩和二阶矩信息，相较于仅用区间变化范围的传统鲁棒调度模型，具有更小的保守性。通过利用更多的信息，使求得的最优鲁棒解在保证鲁棒性的前提下，有更好的系统性能。
3)分布集鲁棒生产调度模型考虑了决策者的风险厌恶特性，可以通过较小的系统平均性能的损失，极大的降低决策者所承担的风险。而且可以通过置信水平的不同设置，来平衡系统性能和鲁棒性之间的关系。使得决策者可以根据当前的需求，设置相应的参数值，以得到最合适的调度策略。
4)两种柯西松弛算法都是多项式时间的迭代下降算法，其中CCRA算法每次迭代的计算复杂度为O(nlogn)，PCRA算法每次迭代的计算复杂度为O(n3)。通过与商业求解器CPLEX进行对比，实验结果表明这两种算法均具有很高的精度以及效率。在工件数达到200个时，CCRA算法和PCRA算法分别能在0.1秒和70秒左右求得精度很高的解，非常适用于求解大规模问题。
具体实施方式
本发明提出了一种生产调度问题的分布集鲁棒模型的建模及优化求解方法，下面结合具体实施例进一步详细说明如下。
本发明提出的一种生产调度问题的分布集鲁棒模型的建模及优化求解方法，其具体实现方式包括以下步骤：
1)构建生产调度问题的分布集鲁棒模型DR-SMSP
本发明所关注的是具有随机加工时间的单机调度问题，针对该问题建立了分布集鲁棒优化模型(DR-SMSP)。在单机调度问题中，所有工件均在一台机器上加工，每个工件均有其特定的加工时间，而且在加工的过程中不能被中断。求解调度问题的目的即为寻找一个工件加工顺序，使得总流经时间、最大拖期时间或者拖期工件个数等某个系统性能指标达到最优。
在DR-SMSP模型中，系统的性能指标选择为总流经时间(Total Flow Time,TFT)，假定所有工件均在加工开始的时刻释放，即释放时间均为零。工件的加工时间具有随机不确定性，随机加工时间的分布未知，但属于一个由支撑集、均值向量和协方差矩阵所确定的分布集中。由于工件的加工时间是随机向量，所有工件加工的总流经时间TFT是一个随机变量，为了综合考虑期望和方差的影响，系统性能TFT的随机度量选取为具有风险厌恶特性的条件风险价值(CVaR，Conditional Value at Risk)。在此种设定下，DR-SMSP模型的目标为寻找一个最优的鲁棒调度序列，使得该序列的TFT在加工时间服从最差分布的情况下具有最小的CVaR。
1-1)确定模型决策变量
该模型的决策变量为可行的调度方案，设模型中有n个工件，且工件的集合为J＝{1,2,…,n}，则一个可行的调度方案可由矩阵X＝{xij,i,j＝1,…,n}表示。其中，如果工件j在加工序列的第i个位置上，则xij＝1，反之xij＝0。
1-2)加工时间的随机性表示
在考虑参数不确定性的鲁棒单机调度问题中，加工时间的不确定性大多是通过区间型的不确定性集来表示的，其中通常会采用一个预算参数来控制结果的保守程度。在本发明的DR-SMSP模型中，我们将不确定加工时间的更多信息考虑进来。把加工时间看成一个随机向量p，其分布Pp是未知的，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集Dp中，该分布集Dp的表达式如式(1)所示：
其中，Sup(pj)表示每个加工时间的支撑集，E(p)和Cov(p)分别表示加工时间的均值向量和协方差矩阵；
1-3)构建模型目标函数
模型的系统性能指标为总流经时间(Total Flow Time,TFT)，在给定一个调度方案X时，TFT由式(2)计算得到：
由于加工时间是随机向量，所有工件的TFT是一个随机变量，为了将期望和方差的影响进行综合的考虑，本发明采用具有风险厌恶特性的条件风险价值(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的概率分布属于一个由支撑集、期望以及方差信息确定的分布集Dz时，在Dz中最差情况下的CVaRα(Z)被定义为鲁棒CVaRα(Z)，即RCVaRα(Z)，其表达式如式(4)所示：
其中，sup表示取集合中的上确界。
基于以上讨论，带有CVaR风险厌恶的分布集鲁棒单机调度问题的目标函数如式(5)所示：
其中，的上标p表明RCVaR所属的分布集为Dp。
1-4)约束条件
本模型包含4类约束条件，其中1类约束是随机加工时间所服从的分布集约束，另外3类约束是调度方案的可行性约束，具体如下所示：
1-4-1)随机加工时间约束；随机加工时间p的分布未知，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集中，表达式如式(6)所示：
1-4-2)工件占用位置约束；
每个工件仅可占用加工序列中的一个位置，如式(7)所示：
1-4-3)可行加工序列位置约束；
可行加工序列中的每个位置仅可被一个工件占用，如式(8)所示：
1-4-4)可行调度方案约束；
可行调度方案X中的每个元素均是0-1变量，如式(9)所示：
xij∈{0,1},i＝1,...,n,j＝1,...,n (9)
约束式(7)至式(9)均是调度方案可行性约束，将其整合到一起，形成调度方案的可行域FB，如式(10)所示：
1-5)单机调度问题分布集鲁棒模型的数学表达式，如式(11)所示：
其中，FB为调度方案的可行域，的上标p表明RCVaR所属的分布集为Dp，min表示在可行域FB中寻找目标函数的最小值，arg表示求得最小目标函数值所对应的最优解X*。
2)对生产调度问题的分布集鲁棒模型进行转化
本模型将式(11)所示的DR-SMSP1模型分解为一个指派问题AP R1和一个整数二阶锥规划问题I-SOCP R2，给出如式(13)所示的在具有半无限支撑集的分布集中计算随机变量RCVaR的显示表达式，具体步骤如下：
2-1)计算随机变量在特定分布集下的RCVaR
本发明给出了在具有半无限支撑集的分布集中计算随机变量的RCVaR的显示表达式如式(13)所示；对于随机变量Z，若其分布函数属于分布集Dz：
则其RCVaRα由式(13)计算得到：
2-2)转换决策变量
为了方便模型的表示与计算，本发明将决策变量由矩阵X转换为向量π，转换关系如式(14)所示：
其中xi为矩阵X的第i行向量。根据此定义形式，π表示工件加工顺序的倒序，即给定一个π＝(π(1),π(2),...,π(n))，π(i)＝j表示工件j在第(n-i+1)次序加工。π相应的可行域如式(15)所示：
进而，TFT表示为π与p的内积，如式(16)所示：
(TFT)f(π,p)＝f(X,p)＝πTp (16)
2-3)分布函数集的映射关系
由于加工时间p的随机性，对每一个确定的π来说，f(π,p)是一个随机变量，记为Fπ。基于p的均值向量和协方差矩阵，Fπ的均值μf(π)和方差分别为：
μf(π)＝πTμ；
进而令Fπ的分布集为：
对于任意一个Dp中的分布Pp，若随机向量p～Pp，则其相应投影随机变量的支撑集为[0,∞)，均值为πTμ＝μf(π)，方差为即的分布在分布集Df中。因此，在分布集Df中求得的是在Dp中求得的的一个上界，即：
进而，DR-SMSP1中的可以用其上界来替代，转化成DR-SMSP2：
根据式(13)所示的RCVaR计算表达式，按式(21)计算得到：
2-4)对生产调度问题的分布集鲁棒模型的分解
由于式(20)所示的DR-SMSP2中包含分段约束，使得问题较难求解，本发明将其分解成两个独立且相对简单的锥规划问题，以避免处理分段的约束。通过严格的理论分析和证明，DR-SMSP2的最优解通过求解一个指派问题(AP)和一个整数的二阶锥规划问题(I-SOCP)获得，分解后的模型如式(22)所示：
其中，
如式(22)-式(24)所示，模型已被分解为一个指派子问题(AP R1)和一个整数二阶锥规划子问题(I-SOCP R2)。在通常情况下，DR-SMSP3模型可以用CPLEX等商业求解器来求解，然而在问题规模较大时，求解整数二阶锥规划子问题(I-SOCP R2)需要耗费巨大的计算时间。因此，本发明通过分析问题的性质，给出了在不同工件的加工时间不相关的情况下，更加高效的求解策略(最短平均加工时间优先准则，SAPT)以及算法(完全柯西松弛算法，CCRA；部分柯西松弛算法，PCRA)。在此情况下，DR-SMSP3中的两个子问题分别转化为式(25)及式(26)的形式：
(AP R1)：
(I-SOCP R2)：
其中，μj与分别为加工时间pj的均值和方差，
一种基于本发明提出的生产调度问题的分布集鲁棒模型的优化求解方法，在加工时间的均值和方差排序一致时，本发明给出了最短平均加工时间优先(SAPT)的求解准则；在加工时间的均值和方差不具备一致性时，针对二阶锥规划子问题(I-SOCP R2)，基于柯西不等式松弛以及交替优化的思路，本发明设计了两种更加通用的柯西松弛求解算法，CCRA算法与PCRA算法各有优点，可以根据实际需要选择使用。具体步骤如下：
1)最短平均加工时间优先(SAPT)准则
以总流经时间(TFT)为优化目标的确定型单机调度问题可以通过最短加工时间优先(SPT)准则精确求解，即按照加工时间由小到大排序的加工序列即为使得TFT最小的最优序列。通过与确定型问题的对比，最短加工时间优先(SPT)准则可以直接引申为最短平均加工时间优先(SAPT)准则，用于求解指派子问题(AP R1)。而对于二阶锥规划子问题(I-SOCP R2)，其最优序列同时受到均值和方差的影响。仅当加工时间的均值和方差排序一致，即较小的均值也拥有较小的方差时，子问题AP R1与I-SOCP R2的最优解相同，进而直接为DR-SMSP3的最优解。以上分析给出了在均值和方差一致的情况下，DR-SMSP3的求解策略：
当工件加工时间的均值和方差排序一致时，通过最短平均加工时间优先(SAPT)准则得到的加工序列即为DR-SMSP3的最优序列。
2)采用完全柯西松弛算法(CCRA)求解
当加工时间的均值和方差不具备一致性时，SAPT准则将不适用于I-SOCP R2。因此，本发明在这种一般的情况下，基于柯西不等式松弛以及交替优化的思路，给出了更加通用的求解算法。
2-1)松弛目标函数
首先引入非负向量来松弛I-SOCP R2的均值部分，对于每个分量j，有以下柯西不等式成立：
当且仅当ηj＝μj/πj时，不等式取等。
对于I-SOCP R2的方差部分，通过引入一个非负变量t∈R+将根号去掉，松弛不等式如式(28)所示：
当且仅当时，不等式取等。
令fA(π)表示I-SOCP R2原来的目标函数，则通过式(27)及式(28)两个不等式，fA(π)被松弛为fB(π,η,t)，松弛目标函数如式(29)所示：
进而，I-SOCP R2被松弛为如式(30)所示：
上述松弛将目标函数中关于π的部分整合到了一起，使得当η和t固定时，可通过排序策略来对π进行寻优。由于R2与之间在最优性方面是等价的，对原问题R2的求解可以转化为求解当π固定的时候，η和t的最优解即为松弛不等式(27)与(28)的取等条件；当η和t固定时，最优的π可以通过ξ(η,t)的非增排序获得，其中ξ(η,t)为fB(π,η,t)中π的系数向量，如式(31)所示：
2-2)采用完全柯西松弛算法求解具体流程如下：
根据步骤2-1)部分的目标松弛与特征分析，本发明设计了一个迭代下降算法来求解并将其命名为完全柯西松弛算法(CCRA)。
在CCRA算法开始前，需给定工件加工时间的均值向量μ以及方差向量σ，并对置信水平α、初始值个数L以及最大迭代次数Itermax设置相应的数值。其中均值向量μ以及方差向量σ根据实际加工时间的历史数据进行选取，置信水平α的默认值为0.95，可根据决策者对风险的厌恶程度进行调整，初始值个数L根据精度要求选取，L越大则算法精度越高，最大迭代次数Itermax则设置为工件个数n的30倍。在设置了初始值个数L之后，初始值的变化步长为θs＝1/(L-1)，第l个初始值由θlμ+(1-θl)σ2的非增排序得到，其中θl＝(l-1)×θs。
CCRA算法分别针对L个初始值进行迭代求解，然后在求得的L组解中选取最优的一组作为算法最终输出的最优解。在针对每个初始值的迭代过程中，对于第m次迭代，首先在固定的πm-1下，求解minη,tfB(πm-1,η,t)，最优解如式(32)所示：
然后在固定的ηm和tm下，求解minπfB(π,ηm,tm)，最优解πm由ξ(ηm,tm)的非增排序获得。如果πm＝πm-1或者m达到了最大迭代次数Itermax，则令πl＝πm，记为第l个初始值对应的解。当L个初始值所对应的解都求得之后，令作为CCRA算法的估计最优值，而π*(α)＝πl为相应的最优解。
可以证明CCRA算法是一个严格下降的算法，由于算法中的主要计算开销发生在对ξ(ηm,tm)的排序中，CCRA算法每次迭代的计算复杂度为O(nLlogn)，其中L为初始值的个数。
3)采用部分柯西松弛算法(PCRA)求解
在CCRA中，关于π的两部分均被松弛了，使得算法中的minπfB(π,ηm,tm)部分直接通过简单的排序来处理。但这种完全松弛引入了n+1个辅助变量，而且松弛后的子问题过于简单，导致算法容易陷入到局部最优中。为了解决此问题，本发明进一步提出了部分柯西松弛算法(PCRA)，其中仅引入一个变量t∈R+，且松弛后的子问题为一个指派问题。
3-1)松弛目标函数
在PCRA中，原问题I-SOCP R2的目标函数fA(π)被松弛为fC(π,t)，如式(33)所示：
且fA(π)≤fC(π,t)的取等条件为
进而，原问题I-SOCP R2被松弛为如式(34)所示：
与R2在最优的情况下是等价的，即若(π*,t*)是的最优解，则π*也是的最优解，因此对原问题R2的求解可以转化为求解
3-2)松弛后的子问题求解
由于xij∈{0,1}，且可以通过式(35)计算：
因此，fC(π,t)可以重新写成关于X和t的函数，如式(36)所示：
令则fC(X,t)可被简化成：
当t固定时，是一个标准的指派问题，可以通过匈牙利(Hungarian)算法精确求解。
3-3)采用部分柯西松弛算法具体流程如下：
当π固定的时候，t的最优解即为松弛不等式(33)的取等条件；当t固定时，子问题是一个如式(37)所示的指派问题。根据这些特性，本发明设计了部分柯西松弛算法(PCRA)求解
在PCRA算法开始前，需给定工件加工时间的均值向量μ以及方差向量σ，并对置信水平α以及最大迭代次数Itermax设置相应的数值。其中均值向量μ以及方差向量σ根据实际加工时间的历史数据进行选取，置信水平α的默认值为0.95，可根据决策者对风险的厌恶程度进行调整，最大迭代次数Itermax则设置为工件个数n的30倍。算法的初始值π0通过μ的非增排序获得。在迭代过程中，对于第m次迭代，首先在固定的πm-1下，求解mintfC(πm-1,t)，最优解为：然后在固定的tm下，求解minXfC(X,tm)，最优解Xm由匈牙利算法求得，再经由式(38)求得πm：
如果πm＝πm-1或者m达到了最大迭代次数Itermax，则令为PCRA算法的估计最优值，X*(α)＝Xm,π*(α)＝πm为相应的最优解。
PCRA算法仅引入了一个变量来松弛原问题，因此这个松弛过程相对于CCRA算法更紧。而且松弛后产生的子问题为指派问题，相较于CCRA算法中的排序策略减少了陷入局部最优的可能。这些改进使得采用PCRA算法估计的解与实际最优解之间的差距变得更小。然而，PCRA算法中的主要计算负担为匈牙利算法的实现，其计算复杂度为O(n3)，高于CCRA算法的O(nLlogn)。因此，PCRA算法相较于CCRA算法，是以更多的计算时间来换取更高的求解精度。
本发明所设计的分布集鲁棒生产调度问题的模型及优化方法的性能分析如下所示：
(1)两种柯西松弛算法的性能分析
表1为本发明所设计的两种柯西松弛算法CCRA、PCRA与商业求解器CPLEX的在求解精度与求解效率上的比较。
表1两种柯西松弛算法CCRA、PCRA与CPLEX的性能比较表
表1中数据表明，在问题规模较小(n小于15)时，CPLEX可以在3小时内求得问题的最优解，而PCRA同样可以得到最优解。当问题规模逐渐增大时，CPLEX在3小时内得不到最优解，而CCRA与PCRA算法得到的解均优于CPLEX在运行3小时后得到的结果，并且需要的计算时间非常少。因此，通过与CPLEX的比较，说明了本发明所提出的算法具有很高的精度以及运算效率，并且适用于大规模问题的求解。
(2)分布集鲁棒调度模型的鲁棒性分析
本发明将分布集鲁棒调度模型求得的鲁棒解与仅考虑均值信息得到的均值解进行对比，来说明分布集鲁棒建模方法所获得的鲁棒性。其中，mu、sig、RCVaR分别表示采用鲁棒解计算得到的TFT的均值、标准差以及鲁棒CVaR；带有“_S”标记的，为采用均值解得到的结果；带有“_Re”标记的，为两种解得到结果的相对差。
表2在5000个加工时间的均值和方差实例下，鲁棒解和均值解TFT统计对比表
表2中结果显示：
1)在相同的置信水平alpha下，鲁棒解相较于均值解使得TFT的均值变大，方差变小，RCVaR也变小。即分布集鲁棒模型以较小均值性能损失为代价，使得TFT的分散性以及RCVaR降低，以达到减少风险的目的。
2)在均值解中，置信水平alpha对TFT均值和方差没有影响，在加工时间的均值在[10,50]，标准差在[1,30]这样的范围中变化时，mu_S基本稳定在1350.29，sig_S基本稳定在351.23。alpha只影响RCVaR的求解，随着alpha越小，RCVaR越小。
3)在鲁棒解中，随着置信水平alpha的减小，TFT的均值在减小，方差在增大；不论在均值，方差，还是RCVaR方面都逐渐接近均值解的结果。即alpha越小，方差的作用越小，两种结果更加接近，当alpha为0时，鲁棒解与均值解是相同的。因此，决策者可以通过设置合适的置信水平参数来平衡期望的系统性能与承担的风险。
(3)模型针对不同分布的鲁棒性分析
表3为在不同分布下鲁棒解得到的TFT结果的统计分析表，说明了由分布集鲁棒模型求得的鲁棒解在不同分布函数设定下的鲁棒性。首先在特定的加工时间的均值及方差的设定下，得到鲁棒序列，由该序列得到的理论上的TFT均值、标准差及RCVaR如表中最后一行所示。进而在不同分布函数设定下，分别产生500,000个加工时间实例，并计算由鲁棒序列得到的实际TFT的均值、标准差及CVaR。
表3在不同分布下鲁棒解得到的TFT结果的统计分析表
由表3中的数据可以看出，在不同分布下得到的均值与标准差结果均在理论值附近，CVaR的值虽然有所不同，但均小于理论上RCVaR的值。由此可以说明，本发明所提出的分布集鲁棒建模方法所得到的最优序列，对加工时间的分布具有很好的鲁棒性。