CN107622338A - 一种基于分布集鲁棒并行机调度模型的生产调度方法 - Google Patents

Abstract

本发明体提出一种基于分布集鲁棒并行机调度模型的生产调度方法，属于生产调度及运筹学领域。该方法首先构建具有风险厌恶特性的分布集鲁棒优化模型DR‑PMSP‑RA，根据模型的目标函数和约束条件，得到初始模型DR‑PMSP‑RA1的表达式；对DR‑PMSP‑RA模型的目标函数进行转化，得到该目标函数的估计上界，并将初始模型转化为估计模型DR‑PMSP‑RA2，该估计模型可分解为两个独立的子模型，对子模型分别求解，更优的子模型解即为整个模型的最优解，得到最优的生产调度方案。通过本方法建立的模型更加符合实际生产的情况，通过利用生产环境中更多的信息，可以在保证系统性能的情况下，降低决策的风险。

Description

一种基于分布集鲁棒并行机调度模型的生产调度方法

技术领域

本发明属于生产调度及运筹学领域，特别涉及一种基于分布集鲁棒并行机调度模型的生产调度方法，考虑在工件加工时间具有随机不确定性的情况下，寻求风险最小的鲁棒调度方案。

背景技术

制造业在我国国民经济发展过程中居于非常重要的地位，其发展状况直接影响我国的综合实力。生产调度是制造系统重要的关键技术，旨在确保生产制造过程高效有序地进行。对生产过程设计合理的计划与调度策略，可以有效缩短产品的生产周期，提高准时交单率，改善设备利用率并降低库存。

在针对同型并行机问题的生产调度中，以往一般采用的是基于确定性模型的生产调度方法。在这种模型中，都假定工件的加工时间以及其他的参数是确定已知的。然而在现实加工过程中，由于机器或刀具条件、工人加工水平和加工环境等因素的影响，工件的加工时间往往存在着不确定性。在实际调度过程中忽略参数的不确定性，会导致模型求得的最优解并不可行。这使得在不确定环境下采用确定型模型得到的调度方案，在实际生产中难以达到决策者的事先预期。因此，针对具有不确定性的并行机调度问题的研究，逐渐引起了学者们的关注。

首先出现的针对调度问题中不确定性的研究，是随机调度问题。在随机调度模型中，不确定性参数被看作是一个分布已知的随机变量，模型的优化目标往往是系统长期的性能期望。从80年代开始，学术界和工业界在不同的方向上展开了随机调度问题的广泛研究。总流经时间、最大拖期时间、加权总拖期时间以及总拖期工件个数等性能指标均已出现在随机调度模型的优化目标中。虽然随机调度模型在理论上对不确定性并行机调度问题的研究有很好的推进作用，但随机调度模型的一些固有缺点，限制了其在实际的大规模生产调度中的应用。这些缺点主要体现在以下几个方面：1)在随机调度模型中，不确定参数的分布是需要精确知道的。然而在实际的生产环境中，随着生产经营日益复杂，小批量个性化定制等生产模式的逐渐转化，很多情况下精确的概率分布很难获得，只能根据类似产品的加工时间来估计它的区间范围。特别是在新产品，尤其是单件或单次加工的新产品方面，这个问题尤为突出。在这种情况下，基于随机调度模型的生产调度方法将不再适用。2)在随机调度模型中，一般采用某种系统性能的期望作为优化目标。这种目标比较适用于制定企业长期的发展计划，而并不适合解决在每次实际运营过程中的利益最大化或风险最小化的问题。3)这类随机调度模型的求解通常是NP-难的，一般只能通过启发式算法或动态规划算法来进行求解，随着问题规模的逐渐扩大，求解随机调度模型的难度将以指数形式增长。

由于随机调度模型具有以上缺点，处理参数不确定性的另一种方法——鲁棒调度模型应运而生。鲁棒调度模型最早由Richard L.Daniels等人提出，其中仅知道区间信息的不确定参数通过区间数据情景进行刻画(一个情景代表不确定参数的一种可能的取值)，这种描述方法相较于随机调度模型中对参数分布函数的描述，更为简单且符合实际。自从Richard L.Daniels将鲁棒优化的思想引入生产调度问题中，近年来在单机调度、并行机调度以及流水车间调度问题中都有相应的发展和研究。目前的鲁棒生产调度问题都是采用基于不确定性集的鲁棒优化方法，不确定性集为有限的离散集合或者是连续的区间形式。在这种模型下，鲁棒调度的关键问题就在于如何定义最差环境，求得每个可行解在最差环境下的鲁棒费用(Robust Cost)以及如何在所有可行解的鲁棒费用中寻求最优。这种基于不确定性集的鲁棒调度模型较为符合实际生产的参数情况，可以寻求到在最差情况下也能有较好系统性能的鲁棒决策，以降低决策的风险。但由于只利用了不确定参数变化范围的边界信息，并主要考虑在最差情况下的系统性能，此种基于不确定性集的鲁棒调度模型所得到的决策可能过于保守，牺牲了在参数常态情况下的系统性能。因此，如何利用历史数据的更多信息，在保证鲁棒性的同时降低决策的保守程度是当前鲁棒调度中急需解决的问题。

发明内容

本发明的目的是为克服已有技术的不足之处，提出一种基于分布集鲁棒并行机调度模型的生产调度方法。本发明建立的模型更加符合实际生产的情况，通过利用生产环境中更多的信息，可以在保证系统性能的情况下，降低决策的风险，得到最优的生产调度方案。

本发明提出的一种基于分布集鲁棒并行机调度模型的生产调度方法，其特征在于，包括以下步骤：

1)构建具有风险厌恶特性的分布集鲁棒优化模型DR-PMSP-RA，得到初始模型DR-PMSP-RA1的表达式；

在DR-PMSP-RA模型中，系统的性能指标选择为总流经时间TFT；假定所有工件均在加工开始的时刻释放，即释放时间均为0，工件的加工时间具有随机不确定性，随机加工时间的分布未知，但属于一个由支撑集、均值向量和协方差矩阵所确定的分布集中；系统性能指标TFT的随机度量选取为条件风险价值CVaR；DR-PMSP-RA模型的目标为寻找一个最优的鲁棒调度方案，该调度方案的TFT在工件加工时间服从最差分布的情况下具有最小的CVaR；

1-1)确定模型决策变量；

DR-PMSP-RA模型的决策变量为可行的调度方案，设该模型中有J个工件和M个机器，工件和机器的集合分别为J＝{1,2,...,J}和M＝{1,2,...,M}，则一个可行的调度方案由一个三维矩阵X∈{0,1}^J×M×J＝{x_jml∈{0,1}|j∈J,m∈M,l∈L＝J}表示；其中，如果工件j被指派到第m个机器上，并以倒数第l的次序加工，则x_jml＝1，反之x_jml＝0；

1-2)加工时间的随机向量表示；

所有工件的加工时间为一个随机向量p，该向量所服从的分布为未知，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集中，该分布集的表达式如式(2)所示：

式中，表示每个工件的加工时间均为非负的，E[p]和Cov[p]分别表示所有工件加工时间向量的均值向量和协方差矩阵；

1-3)构建DR-PMSP-RA模型目标函数；

在给定一个可行调度方案X和所有工件加工时间向量p时，TFT由式(3)计算得到：

式中，p_j表示工件j的加工时间。

所有工件的TFT是一个随机变量，采用具有风险厌恶特性的条件风险价值CVaR作为随机TFT的度量；随机损失Z的CVaR_α表示其在最差1-α概率下的期望，由式(4)计算得到：

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

式中，α∈(0,1)表示CVaR的置信水平，P表示概率取值，inf表示求取集合中的下确界；将在某一分布集上的最大CVaR_α(Z)定义为鲁棒CVaR_α(Z)，即RCVaR_α(Z)；

DR-PMSP-RA模型的目标函数表达式如式(5)所示：

式中，的上标p表明RCVaR所属的所有工件加工时间向量的分布集为D^p，sup表示取集合中的上确界；

1-4)确定DR-PMSP-RA模型的约束条件；

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

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

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

可行调度方案X中的每个元素均是0-1变量，表达式如式(7)所示：

1-4-3)工件占用位置约束；

每个工件仅可占用一台机器上的一个位置，表达式如式(8)所示：

1-4-4)位置被工件占用约束；

每台机器上的每个位置最多可被一个工件占用，表达式如式(9)所示：

1-4-5)排序紧凑约束；

每台机器上被占用的位置是连续的，且从1开始，表达式如式(10)所示：

如式(7)-式(10)所示的后四类约束均是约束调度方案可行性的，将其整合到一起，形成调度方案的可行域X，如式(11)所示：

1-5)建立具有风险厌恶特性的分布集鲁棒并行机初始模型DR-PMSP-RA1的表达式，如式(12)所示：

式中，X为调度方案的可行域，的上标p表明RCVaR所属的分布集为D^p，min表示在可行域X中寻找目标函数的最小值，arg表示求得最小目标函数值所对应的最优解X^*；

2)对DR-PMSP-RA模型的目标函数进行转化与估计；具体步骤如下：

2-1)转换决策变量；

DR-PMSP-RA模型的决策变量由三维矩阵X等价转化为二维矩阵Y，转换关系如式(13)所示：

Y的可行域表达式如式(14)所示：

将二维矩阵Y表示为向量π，表达式如式(15)所示：

π表示忽略机器序号后工件加工顺序的倒序；

π的可行域表达式如式(16)所示：

TFT表示为π与p的内积，表达式如式(17)所示：

f(π,p)＝f(X,p)＝π^Tp。 (17)

2-2)将等价转化为协正定规划；

的值等于如式(18)所示的协正定规划问题RCVaR-COP的最优值：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (18)

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊。

式中，‘·’表示两个矩阵的内积，s.t.代表约束条件，‘±_co 0’表示该符号左侧的矩阵是一个协正定矩阵；

2-3)利用半定松弛得到的估计上界；

通过将式(18)中的两个协正定矩阵约束松弛为正定矩阵约束，RCVaR-COP问题被松弛为一个半正定规划问题RCVaR-SDP，表达式如式(19)所示：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (19)

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊，

式中，‘±0’表示该符号左侧的矩阵是一个半正定矩阵；

令半正定规划问题RCVaR-SDP的最优值为则根据半定松弛的关系，为的一个上界；

2-4)利用分布集映射关系得到的估计上界；

由于所有工件加工时间向量p的随机性，对每一个确定的π来说，f(π,p)是一个随机变量，记为f_π；基于p的均值向量和协方差矩阵，f_π的均值μ_f(π)和方差表达式如式(20)所示：

进而令f_π的分布集为：

对于一维非负随机变量f_π，其RCVaR通过式(22)计算得到：

对于任意一个D^p中的分布，若随机向量p服从于分布集D^p中的分布F，则其相应投影随机变量f_π＝π^Tp的支撑集为[0,∞)，均值为π^Tμ＝μ_f(π)，方差为则f_π的分布在分布集D^f中；因此，在分布集D^f中求得的是在D^p中求得的的一个上界，即：

与之间的关系表达式如式(24)所示：

3)对步骤1)建立的DR-PMSP-RA1模型进行转化；

3-1)替换

利用步骤2)得到的上界替代DR-PMSP-RA1模型转化为估计模型DR-PMSP-RA2，表达式如下：

3-2)分解估计模型DR-PMSP-RA2；

将DR-PMSP-RA2模型分解成两个子模型，分解后的模型表达式如式(26)所示：

式中，DR-PMSP-RA2模型分解后成为DR-PMSP-RA3模型，DR-PMSP-RA3模型包含R₁和R₂两个子模型，子模型R₁的最优解为子模型R₂的最优解

4)对DR-PMSP-RA模型进行求解，得到最优的生产调度方案；

对步骤3-2)分解得到的子模型R₁和子模型R₂求解，分别得到两个子模型的最优解；其中，更小的一个最优解即为DR-PMSP-RA模型的最优解；DR-PMSP-RA模型的最优解为一个最优的向量π值，其对应的忽略机器序号后所有工件加工顺序的倒序即为最优的生产调度方案。本发明所的特点及有益效果如下：

本发明提出的一种基于分布集鲁棒并行机调度模型的生产调度方法，采用基于不确定分布函数集的分布集鲁棒优化方法对同型并行机调度问题进行建模。在分布集鲁棒优化中，不确定参数用随机变量来表示，但是该随机变量的分布函数未知，且属于某个特定的分布函数集合。在优化的过程中，需要考虑该分布函数集合中所有可能的分布函数。虽然在分布集鲁棒优化模型中，不确定参数仍被看作一个随机变量，但相较于随机优化模型，这里并不需要明确分布函数的具体形式，仅需确定一个分布函数的集合即可。而相较于基于不确定性集鲁棒优化模型，分布集鲁棒优化不仅利用了参数变化范围的信息，还将其均值和方差等更多的信息考虑进来，以降低决策的保守程度。因此，将分布集鲁棒优化方法应用于生产调度问题中，比已有的鲁棒建模方法更加符合实际生产的情况，通过利用生产环境中更多的信息，在保证系统性能的情况下，降低决策的风险，得到最优的生产调度方案。

1)本发明将不确定的工件加工时间看作分布未知的随机变量，其分布被限制在由支撑集，均值和协方差矩阵所确定的分布函数集中。系统的性能指标选做工件总流经时间，并采用随机变量的CVaR度量来考虑决策者的风险厌恶特性。

2)本发明所建立的分布集鲁棒并行机调度模型仅需利用随机向量的支撑集、一阶矩和二阶矩信息，并不需要精确的知道其分布信息，这一点相较于随机模型更符合实际生产的情况，实用性更强。

3)本发明所建立的分布集鲁棒并行机调度模型利用了不确定参数的一阶矩和二阶矩信息，相较于仅用区间变化范围的传统鲁棒调度模型，具有更小的保守性。通过利用更多的信息，使求得的最优鲁棒解在保证鲁棒性的前提下，有更好的系统性能。

4)本发明所建立的分布集鲁棒并行机调度模型考虑了决策者的风险厌恶特性，可以通过较小的系统平均性能的损失，极大的降低决策者所承担的风险。而且可以通过置信水平的不同设置，来平衡系统性能和鲁棒性之间的关系。使得决策者可以根据当前的需求，设置相应的参数值，以得到最合适的调度策略。

5)在实际应用中，本发明充分考虑了工件加工时间的不确定性，可以在仅知道工件加工时间的均值和方差的情况下，求得使系统性能最优的调度方案。本发明所建立的模型既不需要精确知道加工时间分布信息，又可以充分利用历史数据，可使得系统即使在分布最差的情况下也有较好的性能，提升了整个生产系统的鲁棒性，降低了在极端情况下损失利益的风险。

具体实施方式

本发明提出的一种具有风险厌恶特性的分布集鲁棒并行机调度建模方法，下面结合具体实施例进一步详细说明如下。

本发明提出的一种基于分布集鲁棒并行机调度模型的生产调度方法，包括以下步骤：

1)针对同型并行机调度问题，构建具有风险厌恶特性的分布集鲁棒优化模型DR-PMSP-RA，得到初始模型DR-PMSP-RA1的表达式；

本发明所关注的是具有随机加工时间的同型并行机调度问题，针对该问题建立了具有风险厌恶特性的分布集鲁棒优化模型(DR-PMSP-RA)。在同型并行机调度问题中，所有工件均可在任一台机器上加工，但每个机器在同一时间仅可加工一个工件，而一个工件也仅可分配给一个机器。每个工件均有其特定的加工时间，而且在加工的过程中不能被中断。求解同型并行机调度问题的目的即为寻找一个调度方案，确定所有工件的机器分配情况以及在每个机器上的加工顺序，使得总流经时间、最大拖期时间或者拖期工件个数等系统性能指标中的某个系统性能指标达到最优。

在DR-PMSP-RA模型中，系统的性能指标选择为总流经时间(Total Flow Time,TFT)，假定所有工件均在加工开始的时刻释放，即释放时间均为0。工件的加工时间具有随机不确定性，随机加工时间的分布未知，但属于一个由支撑集、均值向量和协方差矩阵所确定的分布集中。由于工件的加工时间是随机向量，所有工件加工的总流经时间(TFT)是一个随机变量，为了考虑决策者的风险厌恶特性以及综合考虑期望和方差的影响，系统性能指标TFT的随机度量选取为条件风险价值(CVaR，Conditional Value at Risk)。在此种设定下，DR-PMSP-RA模型的目标为寻找一个最优的鲁棒调度方案，使得该调度方案的TFT在工件加工时间服从最差分布的情况下具有最小的CVaR。

1-1)确定模型决策变量；

DR-PMSP-RA模型的决策变量为可行的调度方案，设该模型中有J个工件和M个机器，工件和机器的集合分别为J＝{1,2,...,J}和M＝{1,2,...,M}，则一个可行的调度方案由一个三维矩阵X∈{0,1}^J×M×J＝{x_jml∈{0,1}|j∈J,m∈M,l∈L＝J}表示；其中，如果工件j被指派到第m个机器上，并以倒数第l的次序加工，则x_jml＝1，反之x_jml＝0。

例如，对于一个具有5个工件，2台机器的并行机调度问题，一个可行的调度方案S_k可写为：

机器1：工件5—工件2—工件4；机器2：工件3—工件1。

根据调度方案与X的对应关系，相应的调度方案S_k的三维矩阵X_k表达式为：

式中，X_k(m＝1)与X_k(m＝2)分别表示机器1对应的二维矩阵和机器2对应的二维矩阵；这两个二维矩阵共同构成了代表可行的调度方案S_k的三维矩阵X_k。

1-2)加工时间的随机向量表示；

在考虑参数不确定性的鲁棒并行机调度问题中，加工时间的不确定性大多是通过区间型的不确定性集来表示的，其中通常会采用一个预算参数来控制结果的保守程度。在本发明的DR-PMSP-RA模型中，我们将不确定加工时间的更多信息考虑进来，将所有工件的加工时间看作一个随机向量p。它所服从的分布是未知的，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集中，该分布集的表达式如式(2)所示：

其中，表示每个工件的加工时间均为非负的，E[p]和Cov[p]分别表示工件加工时间向量的均值向量和协方差矩阵。

1-3)构建DR-PMSP-RA模型目标函数；

DR-PMSP-RA模型的系统性能指标为总流经时间(Total Flow Time,TFT)，在给定一个可行调度方案X和所有工件加工时间向量p时，TFT由式(3)计算得到：

式中，p_j表示工件j的加工时间。

由于所有工件的加工时间p是一个随机向量，所有工件的TFT是一个随机变量，为了将期望和方差的影响进行综合的考虑，本发明采用具有风险厌恶特性的条件风险价值(CVaR，Conditional Valueat Risk)作为随机TFT的度量。随机损失Z的CVaR_α表示其在最差1-α概率下的期望，由式(4)计算得到：

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

式中，α∈(0,1)表示CVaR的置信水平，P表示概率取值，inf表示求取集合中的下确界。将在某一分布集上的最大CVaR_α(Z)定义为鲁棒CVaR_α(Z)，即RCVaR_α(Z)；

带有CVaR风险厌恶特性的分布集鲁棒并行机调度模型(即DR-PMSP-RA模型)的目标函数表达式如式(5)所示：

式中，的上标p表明RCVaR所属的所有工件加工时间向量的分布集为D^p，sup表示取集合中的上确界。

1-4)确定DR-PMSP-RA模型的约束条件；

DR-PMSP-RA模型包含5个约束条件，其中1个是约束随机加工时间所服从的分布集，另外4个是约束调度方案的可行性，具体如下所示：

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

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

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

可行调度方案X中的每个元素均是0-1变量，表达式如式(7)所示：

1-4-3)工件占用位置约束；

每个工件仅可占用一台机器上的一个位置，表达式如式(8)所示：

1-4-4)位置被工件占用约束；

每台机器上的每个位置最多可被一个工件占用，表达式如式(9)所示：

1-4-5)排序紧凑约束；

每台机器上被占用的位置必须是连续的，且从1开始，表达式如式(10)所示：

如式(7)-式(10)所示的后四类约束均是约束调度方案可行性的，将其整合到一起，形成调度方案的可行域X，如式(11)所示：

1-5)建立具有风险厌恶特性的分布集鲁棒并行机初始模型DR-PMSP-RA1的表达式，如式(12)所示：

式中，X为调度方案的可行域，的上标p表明RCVaR所属的分布集为D^p，min表示在可行域X中寻找目标函数的最小值，arg表示求得最小目标函数值所对应的最优解X^*(X^*即最优的可行调度方案)。

式(12)所建立的表达式为本发明建立的初始模型，下面通过步骤(2)和(3)对该初始模型进行转化，将步骤1)建立的模型转化为更容易求解的形式。

2)对DR-PMSP-RA模型的目标函数进行转化与估计；

本发明通过对的对偶问题进行分析，将其等价转化为一个协正定规划(COP)问题，并分别通过半定松弛(SDR)以及分布函数集的多维到一维映射关系，给出了两个可以求解的估计模型，具体步骤如下：

2-1)转换决策变量；

由于为工件分配的机器序号对目标函数值没有影响，DR-PMSP-RA模型的决策变量可以由三维矩阵X等价转化为二维矩阵Y，转换关系如式(13)所示：

关于Y的可行域可直接根据X写出，表达式如式(14)所示：

为了进一步方便模型的表示与计算，二维矩阵Y表示为向量π，表达式如式(15)所示：

π表示忽略机器序号后工件加工顺序的倒序，例如π_j＝l表示工件j在某一台机器上的倒数第l次序加工。

π的可行域表达式如式(16)所示：

TFT可表示为π与p的内积，表达式如式(17)所示：

f(π,p)＝f(X,p)＝π^Tp。 (17)

2-2)将等价转化为协正定规划；

根据CVaR的定义以及优化的对偶原理，的值等于如式(18)所示的协正定规划问题RCVaR-COP的最优值：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (18)

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊。

式中，‘·’表示两个矩阵的内积，s.t.代表约束条件，‘±_co 0’表示该符号左侧的矩阵是一个协正定矩阵。

2-3)利用半定松弛得到的估计上界；

虽然求解的值可以转化为协正定规划问题RCVaR-COP，由于判断一个矩阵是否为协正定矩阵是NP-难的，其精确值仍然很难获得。本发明通过对RCVaR-COP问题进行半定松弛得到一个估计上界。

通过将式(18)中的两个协正定矩阵约束(即式(18)的约束条件中带有±_co 0的两个约束)松弛为正定矩阵约束，RCVaR-COP问题被松弛为一个半正定规划问题RCVaR-SDP，表达式如式(19)所示：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (19)

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊，

式中，‘±0’表示该符号左侧的矩阵是一个半正定矩阵。

令半正定规划问题RCVaR-SDP的最优值为则根据半定松弛的关系，为的一个上界。此上界相当于将分布集约束式(6)中的非负约束松弛掉了，在实际应用中可以为提供较好的估计。

2-4)利用分布集映射关系得到的估计上界；

由于所有工件加工时间向量p的随机性，对每一个确定的π来说，f(π,p)是一个随机变量，记为f_π。基于p的均值向量和协方差矩阵，f_π的均值μ_f(π)和方差表达式如式(20)所示：

进而令f_π的分布集为：

对于一维非负随机变量f_π，其RCVaR可以通过式(22)计算得到：

对于任意一个D^p中的分布，若随机向量p服从于分布集D^p中的分布F，则其相应投影随机变量f_π＝π^Tp的支撑集为[0,∞)，均值为π^Tμ＝μ_f(π)，方差为这说明f_π的分布在分布集D^f中。因此，在分布集D^f中求得的是在D^p中求得的的一个上界，即：

由于保留了分布集D^p中的非负约束，它相较于半定松弛得到的上界有更好的估计效果，与之间的关系表达式如式(24)所示：

3)对步骤1)建立的DR-PMSP-RA1模型进行转化，将DR-PMSP-RA1模型转化为可求解的估计模型；

3-1)替换

在实际应用中，对的估计较为精确，且在某些特殊的均值与协方差情况下，两者等价。因此，通过将用其上界来替代，DR-PMSP-RA1模型转化为估计模型DR-PMSP-RA2，表达式如下：

3-2)分解估计模型DR-PMSP-RA2；

由于式(22)所示的表达式中包含分段约束，使得问题较难求解，本发明通过对其性质的分析，将DR-PMSP-RA2模型分解成两个独立且相对简单的子模型，以避免处理分段的约束。通过严格的理论分析和证明，DR-PMSP-RA2模型的最优解可通过分别求解一个混合整数线性规划子模型和一个整数二阶锥规划(I-SOCP)子模型获得，分解后的模型表达式如式(26)所示：

式中，DR-PMSP-RA2模型分解后成为DR-PMSP-RA3模型，DR-PMSP-RA3模型包含R₁和R₂两个子模型，子模型R₁的最优解为子模型R₂的最优解

4)对DR-PMSP-RA模型进行求解，得到最优的生产调度方案；

步骤3-2)已经将步骤1)建立的模型分解为两个可独立求解的子模型，因此对式(26)所示的两个子模型分别求解，分别得到两个向量π值，再选取其中向量π值更小的子模型的最优解的作为DR-PMSP-RA3模型的最优解即可，该最优解即为DR-PMSP-RA模型的最优解。其中，子模型R₁为混合整数线性规划模型，可通过CPLEX等求解器直接求解；子模型R₂为整数的二阶锥规划模型，既可通过CPLEX求解，也可以根据问题的性质设计更高效的迭代下降算法进行求解。

DR-PMSP-RA模型最终求解出一个最优的向量π值。向量π表示忽略机器序号后所有工件加工顺序的倒序，通过最优的向量π即可给出一个完整的最优生产调度方案。

例如，对于一个具有5个工件，2台机器的同型并行机调度问题，若DR-PMSP-RA模型求得的最优向量π值为π＝(1，2，1，3，2)，则对应的一个最优调度方案为：

机器1：工件4—工件2—工件1；机器2：工件5—工件3。

由于在同型并行机中，所有的机器均相同且加工效果等价，因此在确定最优调度方案时，仅需保证工件的加工次序符合最优解π的要求即可，而对于机器序号来说，则可以在保证每个机器在同一时间仅加工一个工件的前提下任意选择。

本发明的一个具体实施例将通过本发明建立的分布集鲁棒并行机调度模型求得的鲁棒解与仅考虑均值信息得到的均值解进行对比，以说明本发明方法可提升系统的鲁棒性，具体分析如下所示：

本实施例中，随机加工时间的均值在区间[10,60]中任意选取，方差通过一个相关性水平ρ来控制，选取为对于每一对均值和方差实例，分布集鲁棒并行机调度模型可以得到一个鲁棒解，其相应的理论上的TFT可由式(20)直接算出。表1展示了在ρ＝1情况下进行5000次实例所得到的结果，给出了不同α取值下鲁棒解所获得的鲁棒性提升以及均值损失。

表1在不同置信水平下，鲁棒解和均值解TFT统计对比表

表1中结果显示：

1)在相同的置信水平下，通过本发明方法得到的鲁棒解相较于均值解使得TFT的均值变大，方差变小。即分布集鲁棒模型以较小均值性能损失为代价，使得TFT的分散性降低，以达到减少风险的目的。

2)在均值解中，置信水平对TFT均值和方差没有影响。在鲁棒解中，随着置信水平的减小，TFT的均值在减小，方差在增大；不论在均值还是方差方面都逐渐接近均值解的结果。即置信水平越小，方差的作用越小，两种结果更加接近，当置信水平为0时，鲁棒解与均值解是相同的。因此，决策者可以通过设置合适的置信水平参数来平衡期望的系统性能与承担的风险。

表2展示了在α＝0.95情况下进行5000次实例所得到的结果，给出了不同相关性水平设置下鲁棒解所获得的鲁棒性提升以及均值损失。

表2在不同均值和方差相关性水平下，鲁棒解和均值解TFT统计对比表

ρ	鲁棒解	均值解	平均损失	最大损失	鲁棒解	均值解	平均提升	最大提升
									0	13133.34	12071.39	8.09％	14.45％	1855.26	2423.83	30.65％	61.08％
1	12352.84	12066.22	2.32％	4.80％	910.39	1086.97	19.40％	42.42％
									2	12322.11	12060.49	2.12％	4.22％	1084.06	1252.07	15.50％	31.81％
3	12290.87	12055.97	1.91％	4.10％	1261.48	1419.03	12.49％	26.06％
									4	12268.99	12057.72	1.72％	3.51％	1441.42	1589.38	10.26％	21.12％
5	12266.65	12075.87	1.56％	3.03％	1626.21	1765.72	8.58％	16.46％
									6	12222.78	12051.93	1.40％	3.48％	1801.83	1931.55	7.20％	15.09％
7	12226.07	12070.39	1.27％	2.98％	1987.96	2110.45	6.16％	12.29％
									8	12210.65	12068.85	1.16％	2.40％	2170.65	2286.48	5.34％	11.44％
9	12078.21	12078.21	0％	0％	1617.57	1617.57	0％	0％

表2中结果显示：

在所有相关性水平的设置下，鲁棒解均可以达到显著降低TFT方差的效果，而在均值方面的损失一直保持在较低的水平。另外，随着相关性水平越来越高，鲁棒解与均值解的效果越来越接近，说明了在方差越不可控的情况下，本发明所设计的分布集鲁棒模型的效果更加显著。

Claims (1)

1.一种基于分布集鲁棒并行机调度模型的生产调度方法，其特征在于，包括以下步骤：

1)构建具有风险厌恶特性的分布集鲁棒优化模型DR-PMSP-RA，得到初始模型DR-PMSP-RA1的表达式；

在DR-PMSP-RA模型中，系统的性能指标选择为总流经时间TFT；假定所有工件均在加工开始的时刻释放，即释放时间均为0，工件的加工时间具有随机不确定性，随机加工时间的分布未知，但属于一个由支撑集、均值向量和协方差矩阵所确定的分布集中；系统性能指标TFT的随机度量选取为条件风险价值CVaR；DR-PMSP-RA模型的目标为寻找一个最优的鲁棒调度方案，该调度方案的TFT在工件加工时间服从最差分布的情况下具有最小的CVaR；

1-1)确定模型决策变量；

DR-PMSP-RA模型的决策变量为可行的调度方案，设该模型中有J个工件和M个机器，工件和机器的集合分别为J＝{1,2,...,J}和M＝{1,2,...,M}，则一个可行的调度方案由一个三维矩阵X∈{0,1}^J×M×J＝{x_jml∈{0,1}|j∈J,m∈M,l∈L＝J}表示；其中，如果工件j被指派到第m个机器上，并以倒数第l的次序加工，则x_jml＝1，反之x_jml＝0；

1-2)加工时间的随机向量表示；

所有工件的加工时间为一个随机向量p，该向量所服从的分布为未知，但属于一个由支撑集、均值向量和协方差矩阵确定的分布集中，该分布集的表达式如式(2)所示：

<mrow> <msup> <mi>D</mi> <mi>p</mi> </msup> <mo>=</mo> <mo>{</mo> <mi>F</mi> <mo>|</mo> <mi>P</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>&amp;Element;</mo> <msubsup> <mi>R</mi> <mo>+</mo> <mi>J</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <mi>p</mi> <mo>&amp;rsqb;</mo> <mo>=</mo> <mi>&amp;mu;</mi> <mo>,</mo> <mi>C</mi> <mi>o</mi> <mi>v</mi> <mo>&amp;lsqb;</mo> <mi>p</mi> <mo>&amp;rsqb;</mo> <mo>=</mo> <mo>&amp;Sigma;</mo> <mo>}</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

式中，表示每个工件的加工时间均为非负的，E[p]和Cov[p]分别表示所有工件加工时间向量的均值向量和协方差矩阵；

1-3)构建DR-PMSP-RA模型目标函数；

在给定一个可行调度方案X和所有工件加工时间向量p时，TFT由式(3)计算得到：

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>,</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>lp</mi> <mi>j</mi> </msub> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

式中，p_j表示工件j的加工时间。

所有工件的TFT是一个随机变量，采用具有风险厌恶特性的条件风险价值CVaR作为随机TFT的度量；随机损失Z的CVaR_α表示其在最差1-α概率下的期望，由式(4)计算得到：

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

式中，α∈(0,1)表示CVaR的置信水平，P表示概率取值，inf表示求取集合中的下确界；将在某一分布集上的最大CVaR_α(Z)定义为鲁棒CVaR_α(Z)，即RCVaR_α(Z)；

DR-PMSP-RA模型的目标函数表达式如式(5)所示：

<mrow> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>X</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mi>s</mi> <mi>u</mi> <mi>p</mi> </mrow> <mrow> <mi>F</mi> <mo>&amp;Element;</mo> <msup> <mi>D</mi> <mi>p</mi> </msup> </mrow> </munder> <msub> <mi>CVaR</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>X</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>

式中，的上标p表明RCVaR所属的所有工件加工时间向量的分布集为D^p，sup表示取集合中的上确界；

1-4)确定DR-PMSP-RA模型的约束条件；

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

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

<mrow> <msup> <mi>D</mi> <mi>p</mi> </msup> <mo>=</mo> <mo>{</mo> <mi>F</mi> <mo>|</mo> <mi>P</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>&amp;Element;</mo> <msubsup> <mi>R</mi> <mo>+</mo> <mi>J</mi> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <mi>p</mi> <mo>&amp;rsqb;</mo> <mo>=</mo> <mi>&amp;mu;</mi> <mo>,</mo> <mi>C</mi> <mi>o</mi> <mi>v</mi> <mo>&amp;lsqb;</mo> <mi>p</mi> <mo>&amp;rsqb;</mo> <mo>=</mo> <mo>&amp;Sigma;</mo> <mo>}</mo> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

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

可行调度方案X中的每个元素均是0-1变量，表达式如式(7)所示：

<mrow> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;Element;</mo> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>}</mo> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>m</mi> <mo>&amp;Element;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>&amp;Element;</mo> <mi>L</mi> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

1-4-3)工件占用位置约束；

每个工件仅可占用一台机器上的一个位置，表达式如式(8)所示：

<mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

1-4-4)位置被工件占用约束；

每台机器上的每个位置最多可被一个工件占用，表达式如式(9)所示：

<mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;le;</mo> <mn>1</mn> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>m</mi> <mo>&amp;Element;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>&amp;Element;</mo> <mi>L</mi> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>

1-4-5)排序紧凑约束；

每台机器上被占用的位置是连续的，且从1开始，表达式如式(10)所示：

<mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;GreaterEqual;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>m</mi> <mo>&amp;Element;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>J</mi> <mo>-</mo> <mn>1</mn> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

如式(7)-式(10)所示的后四类约束均是约束调度方案可行性的，将其整合到一起，形成调度方案的可行域X，如式(11)所示：

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>X</mi> <mo>=</mo> <mo>{</mo> <mi>X</mi> <mo>&amp;Element;</mo> <msup> <mrow> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>}</mo> </mrow> <mrow> <mi>J</mi> <mo>&amp;times;</mo> <mi>M</mi> <mo>&amp;times;</mo> <mi>J</mi> </mrow> </msup> <mo>|</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;le;</mo> <mn>1</mn> <mo>;</mo> <mo>&amp;ForAll;</mo> <mi>m</mi> <mo>&amp;Element;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>&amp;Element;</mo> <mi>L</mi> <mo>;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;GreaterEqual;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>m</mi> <mo>&amp;Element;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>J</mi> <mo>-</mo> <mn>1</mn> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>

1-5)建立具有风险厌恶特性的分布集鲁棒并行机初始模型DR-PMSP-RA1的表达式，如式(12)所示：

<mrow> <msup> <mi>X</mi> <mo>*</mo> </msup> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>X</mi> <mo>&amp;Element;</mo> <mi>X</mi> </mrow> </munder> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>X</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>

式中，X为调度方案的可行域，的上标p表明RCVaR所属的分布集为D^p，min表示在可行域X中寻找目标函数的最小值，arg表示求得最小目标函数值所对应的最优解X^*；

2)对DR-PMSP-RA模型的目标函数进行转化与估计；具体步骤如下：

2-1)转换决策变量；

DR-PMSP-RA模型的决策变量由三维矩阵X等价转化为二维矩阵Y，转换关系如式(13)所示：

<mrow> <msub> <mi>y</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>j</mi> <mi>m</mi> <mi>l</mi> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>&amp;Element;</mo> <mi>L</mi> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>

Y的可行域表达式如式(14)所示：

<mrow> <mi>Y</mi> <mo>=</mo> <mo>{</mo> <mi>Y</mi> <mo>&amp;Element;</mo> <msup> <mrow> <mo>{</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>}</mo> </mrow> <mrow> <mi>J</mi> <mo>&amp;times;</mo> <mi>J</mi> </mrow> </msup> <mo>|</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>y</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>y</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;le;</mo> <mi>M</mi> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>&amp;Element;</mo> <mi>L</mi> <mo>;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>y</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>&amp;GreaterEqual;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>y</mi> <mrow> <mi>j</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>J</mi> <mo>-</mo> <mn>1</mn> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

将二维矩阵Y表示为向量π，表达式如式(15)所示：

<mrow> <msub> <mi>&amp;pi;</mi> <mi>j</mi> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>ly</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>

π表示忽略机器序号后工件加工顺序的倒序；

π的可行域表达式如式(16)所示：

<mrow> <mo>&amp;Pi;</mo> <mo>=</mo> <mo>{</mo> <mi>&amp;pi;</mi> <mo>&amp;Element;</mo> <msup> <mi>R</mi> <mi>J</mi> </msup> <mo>|</mo> <msub> <mi>&amp;pi;</mi> <mi>j</mi> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>ly</mi> <mrow> <mi>j</mi> <mi>l</mi> </mrow> </msub> <mo>,</mo> <mo>&amp;ForAll;</mo> <mi>j</mi> <mo>&amp;Element;</mo> <mi>J</mi> <mo>;</mo> <mi>Y</mi> <mo>&amp;Element;</mo> <mi>Y</mi> <mo>}</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>

TFT表示为π与p的内积，表达式如式(17)所示：

f(π,p)＝f(X,p)＝π^Tp。 (17)

2-2)将等价转化为协正定规划；

的值等于如式(18)所示的协正定规划问题RCVaR-COP的最优值：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (18)

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> </mrow> </mtd> <mtd> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>k</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>-</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>-</mo> <mi>&amp;pi;</mi> <mo>)</mo> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mi>Z</mi> </mtd> </mtr> </mtable> </mfenced> </mtd> </mtr> </mtable> <msub> <mo>&amp;PlusMinus;</mo> <mrow> <mi>c</mi> <mi>o</mi> </mrow> </msub> <mn>0</mn> <mo>,</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>r</mi> <mn>0</mn> </msub> </mtd> <mtd> <mrow> <msup> <msub> <mi>r</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mi>Z</mi> </mtd> </mtr> </mtable> </mfenced> <msub> <mo>&amp;PlusMinus;</mo> <mrow> <mi>c</mi> <mi>o</mi> </mrow> </msub> <mn>0</mn> <mo>,</mo> </mrow>

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊。

式中，‘·’表示两个矩阵的内积，s.t.代表约束条件，‘±_co 0’表示该符号左侧的矩阵是一个协正定矩阵；

2-3)利用半定松弛得到的估计上界；

通过将式(18)中的两个协正定矩阵约束松弛为正定矩阵约束，RCVaR-COP问题被松弛为一个半正定规划问题RCVaR-SDP，表达式如式(19)所示：

min k+(1-α)^-1[r₀+μ^Tr₁+(∑+μμ^T)·Z] (19)

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> </mrow> </mtd> <mtd> <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>k</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>-</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>-</mo> <mi>&amp;pi;</mi> <mo>)</mo> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mi>Z</mi> </mtd> </mtr> </mtable> </mfenced> <mo>&amp;PlusMinus;</mo> <mn>0</mn> <mo>,</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>r</mi> <mn>0</mn> </msub> </mtd> <mtd> <mrow> <msup> <msub> <mi>r</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mn>1</mn> </msub> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mi>Z</mi> </mtd> </mtr> </mtable> </mfenced> <mo>&amp;PlusMinus;</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> <mo>,</mo> </mrow>

r₀∈R,r₁∈R^J,Z∈R^J×J,k∈R₊，

式中，‘±0’表示该符号左侧的矩阵是一个半正定矩阵；

令半正定规划问题RCVaR-SDP的最优值为则根据半定松弛的关系，为的一个上界；

2-4)利用分布集映射关系得到的估计上界；

由于所有工件加工时间向量p的随机性，对每一个确定的π来说，f(π,p)是一个随机变量，记为f_π；基于p的均值向量和协方差矩阵，f_π的均值μ_f(π)和方差表达式如式(20)所示：

<mrow> <msub> <mi>&amp;mu;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&amp;pi;</mi> <mi>T</mi> </msup> <mi>&amp;mu;</mi> <mo>;</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&amp;pi;</mi> <mi>T</mi> </msup> <mi>&amp;Sigma;&amp;pi;</mi> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>

进而令f_π的分布集为：

<mrow> <msup> <mi>D</mi> <mi>f</mi> </msup> <mo>=</mo> <mo>{</mo> <msup> <mi>F</mi> <mi>f</mi> </msup> <mo>|</mo> <mi>P</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>&amp;Element;</mo> <msub> <mi>R</mi> <mo>+</mo> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>&amp;rsqb;</mo> <mo>=</mo> <msub> <mi>&amp;mu;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>V</mi> <mi>a</mi> <mi>r</mi> <mo>&amp;lsqb;</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>&amp;rsqb;</mo> <mo>=</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow>

对于一维非负随机变量f_π，其RCVaR通过式(22)计算得到：

<mrow> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>f</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <mi>&amp;alpha;</mi> </mrow> </mfrac> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> </mrow> </mtd> <mtd> <mrow> <mn>0</mn> <mo>&amp;le;</mo> <mi>&amp;alpha;</mi> <mo>&amp;le;</mo> <mfrac> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;mu;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;mu;</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>+</mo> <msqrt> <mfrac> <mi>&amp;alpha;</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>&amp;alpha;</mi> </mrow> </mfrac> </msqrt> <mo>&amp;CenterDot;</mo> <msqrt> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> </mrow> </mtd> <mtd> <mrow> <mfrac> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;mu;</mi> <mi>f</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&amp;le;</mo> <mi>&amp;alpha;</mi> <mo>&amp;le;</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>

对于任意一个D^p中的分布，若随机向量p服从于分布集D^p中的分布F，则其相应投影随机变量f_π＝π^Tp的支撑集为[0,∞)，均值为π^Tμ＝μ_f(π)，方差为则f_π的分布在分布集D^f中；因此，在分布集D^f中求得的是在D^p中求得的的一个上界，即：

<mrow> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>f</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mi>s</mi> <mi>u</mi> <mi>p</mi> </mrow> <mrow> <msup> <mi>F</mi> <mi>f</mi> </msup> <mo>&amp;Element;</mo> <msup> <mi>D</mi> <mi>f</mi> </msup> </mrow> </munder> <msub> <mi>CVaR</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;GreaterEqual;</mo> <munder> <mrow> <mi>s</mi> <mi>u</mi> <mi>p</mi> </mrow> <mrow> <mi>F</mi> <mo>&amp;Element;</mo> <msup> <mi>D</mi> <mi>p</mi> </msup> </mrow> </munder> <msub> <mi>CVaR</mi> <mi>&amp;alpha;</mi> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>&amp;pi;</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>&amp;pi;</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow>

与之间的关系表达式如式(24)所示：

<mrow> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>&amp;pi;</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>&amp;le;</mo> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>f</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;le;</mo> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>p</mi> </msubsup> <msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>(</mo> <mrow> <mi>&amp;pi;</mi> <mo>,</mo> <mi>p</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mrow> <mi>S</mi> <mi>D</mi> <mi>P</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow>

3)对步骤1)建立的DR-PMSP-RA1模型进行转化；

3-1)替换

利用步骤2)得到的上界替代DR-PMSP-RA1模型转化为估计模型DR-PMSP-RA2，表达式如下：

<mrow> <msup> <mi>&amp;pi;</mi> <mo>*</mo> </msup> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>&amp;pi;</mi> <mo>&amp;Element;</mo> <mo>&amp;Pi;</mo> </mrow> </munder> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>f</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>25</mn> <mo>)</mo> </mrow> </mrow>

3-2)分解估计模型DR-PMSP-RA2；

将DR-PMSP-RA2模型分解成两个子模型，分解后的模型表达式如式(26)所示：

<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mrow> <mi>&amp;pi;</mi> <mo>&amp;Element;</mo> <mo>&amp;Pi;</mo> </mrow> </munder> <msubsup> <mi>RCVaR</mi> <mi>&amp;alpha;</mi> <mi>f</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>&amp;pi;</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mo>{</mo> <msubsup> <mi>R</mi> <mn>1</mn> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>R</mi> <mn>2</mn> <mo>*</mo> </msubsup> <mo>}</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>26</mn> <mo>)</mo> </mrow> </mrow>

式中，DR-PMSP-RA2模型分解后成为DR-PMSP-RA3模型，DR-PMSP-RA3模型包含R₁和R₂两个子模型，子模型R₁的最优解为子模型R₂的最优解

4)对DR-PMSP-RA模型进行求解，得到最优的生产调度方案；

对步骤3-2)分解得到的子模型R₁和子模型R₂求解，分别得到两个子模型的最优解；其中，更小的一个最优解即为DR-PMSP-RA模型的最优解；DR-PMSP-RA模型的最优解为一个最优的向量π值，其对应的忽略机器序号后所有工件加工顺序的倒序即为最优的生产调度方案。