CN105929690B - 一种基于分解多目标进化算法的柔性车间鲁棒调度方法 - Google Patents

Abstract

本发明公开了一种基于分解多目标进化算法的柔性车间鲁棒调度方法,(1)读取柔性作业车间的作业和机器属性等输入信息；定义优化目标，设定约束条件；(2)初始化算法的参数；(3)确定每个子问题的邻域，产生初始父代群体，从初始群体中确定出所有的Pareto非支配解构成外部存储器；(4)生成子代群体。进行交配选择，采用自适应变异算子和基于修复的交叉算子繁殖子代个体，并更新外部存储器；(5)利用生成的子代群体对各子问题的当前最优个体进行更新，构成新的父代群体；(6)判断个体目标评价次数若达到最大，则输出外部存储器，即一组Pareto非支配的柔性作业车间调度解；未达到则跳转至(4)。本发明快速高效地实现柔性作业车间中的调度任务。

Description

一种基于分解多目标进化算法的柔性车间鲁棒调度方法

技术领域

本发明涉及柔性作业车间调度控制技术领域，特别是一种基于分解多目标进化算法的柔性车间鲁棒调度方法。

背景技术

柔性作业车间调度问题是对经典车间调度问题的一般化，它是一类NP难问题。在柔性作业车间调度问题中，允许一道工序被多台机器加工，因此，需为每项作业的每道工序分配适当的机器，并确定各机器上工序的加工顺序，以在满足各种约束条件的前提下，实现作业的完工时间最短、各机器的负载均衡等优化目标。

实际柔性作业车间的生产环境中存在着很多不确定因素，如由于对作业规格要求的改变导致部分工序加工时间发生改变；由于转包商的运输延迟，导致作业的投放时间推迟等。当面临这些扰动时，根据初始数据产生的最优调度方案的性能可能大大降低，因此亟需研究一种能够处理不确定因素的新型柔性作业车间鲁棒调度方法，以在保证完工时间较短、各机器的负载较为均衡等性能的同时，增强车间调度方案对不确定性的抗干扰能力，降低调度性能对扰动的敏感程度。柔性作业车间的鲁棒调度技术对实际生产***的成功实施具有重要意义。

进化算法是模拟生物在自然环境中的进化过程而形成的一类自适应全局优化概率搜索算法。进化算法可以处理传统优化方法难以解决的复杂优化问题，例如非连续、多模态等问题，它对整个群体实施选择、交叉、变异等操作，可以在算法的一次运行中并行搜索到多个解，加之它具有较强的环境自适应能力，因此，进化算法特别适用于求解柔性作业车间调度这类同时存在多个Pareto非支配解的多目标鲁棒优化问题。基于分解的多目标进化算法(MOEA/D)是由张青富教授和李辉博士在2007年首次提出的。该算法给出了一种新颖且通用的多目标处理框架。它将多目标优化问题分解为一组单目标子问题，并采用基于群体的进化算法协同求解这些子问题。由于它的高效性，近年来，MOEA/D已经成为国际学术界的研究热点之一，并且已有学者将MOEA/D成功应用于实际优化问题的求解中。

目前已有的柔性作业车间调度方法存在以下不足：

(1)大多仅考虑了静态的生产环境。它们假设柔性作业车间中的所有信息都是预先可知且确定不变的，显然，当实际生产环境存在不确定因素时，依据静态方法产生的调度方案不再适用。

(2)对多个优化目标的处理方式比较单一。大多已有方法采用加权求和法将多个目标转换为一个目标，这种方法将会引入较多的参数，并且需要事先对各个目标进行归一化处理。由于柔性作业车间调度问题的多个目标之间往往是相互矛盾的，因此更好的方式是采用多目标进化算法对多个目标并行处理，从而为生产管理者提供一组反映目标间不同折中程度的调度方案，为其做出最终的决策提供参考。

(3)局部搜索能力较弱，易于陷入局部最优，调度效率低下。

发明内容

本发明所要解决的技术问题是克服现有技术的不足而提供一种基于分解多目标进化算法的柔性车间鲁棒调度方法，以在实际的车间生产环境中，在优化作业的完工时间、促使各机器负载均衡的同时，尽可能地降低不确定性对调度方案性能的负面影响，并提高调度方法的局部搜索能力，避免陷入局部最优，实现快速高效的车间调度。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，包括以下步骤：

步骤(1)、读取柔性作业车间的输入信息，定义优化目标，设定约束条件：

车间的输入信息包括每项作业的工序数、每项作业中每道工序可分配的机器集合、各工序在相应机器上的加工时间、正常工作的机器数目；优化目标包括作业完工时间、机器的最大负载、调度性能的鲁棒性；约束条件包括工序优先级约束和禁止抢先占有的约束；

步骤(2)、初始化基于分解的多目标进化算法的参数：

设置迭代次数NmbEvl、群体规模N即为子问题的个数、产生N个均匀分布的目标权向量λ¹,…,λ^N、子问题的邻域规模T，繁殖算子中的父代个体从邻域里选取的概率δ、邻域个体允许被每个子代个体替换的最大数目n_r，评价鲁棒性时采样的不确定场景个数Q，以及交叉概率CR；

步骤(3)、确定每个子问题的邻域，产生初始父代群体：

(3.1)、对第i个子问题，i＝1,…,N，确定邻域B(i)＝{i₁,…,i_T}；分别计算第i个子问题的权向量λⁱ和其它(N-1)个子问题权向量间的欧拉距离，将T个距离λⁱ最近的权向量所对应的子群体序号{i₁,…,i_T}构成B(i)；

(3.2)、随机产生初始父代群体{x¹,…,x^N}，其中，个体xⁱ为第i个子群体的当前解；群体中的每个个体包括工序序列向量和机器分配向量；随机采样一组不确定性场景θ_q，q＝1,2,...,Q；计算初始父代群体中每个个体的目标值f₁＝makespan_I、f₂＝workload_I和f₃＝robustness；makespan_I和workload_I分别表示初始景象下的作业完工时间和机器最大负载，robustness表示调度性能的鲁棒性；从初始父代群体中确定出所有的Pareto非支配解构成外部存储器Arc；设置目标评价次数计数变量ct＝N；归一化各目标值，将每个个体xⁱ的第j个目标值f_j(xⁱ)映射到区间[0,1]之间，产生xⁱ在第j个目标上的归一化值fno_j(xⁱ)，即：

其中，和分别表示当前父代群体中的所有个体在第j个目标上的最大值和最小值，j＝1,2,3；

(3.3)、设置进化代数t＝0；

步骤(4)、生成子代群体：

随机采样一组不确定性场景θ_q；设置子代群体令i＝1；

(4.1)、交配选择；

产生一个均匀分布的随机数rand₁∈[0,1]；设置第i个子问题的更新邻域P(i)：

基于以下规则产生3个不同的父代个体和

(4.2)、繁殖；

采用基于微分进化算法的繁殖算子，由和xⁱ生成子代个体uⁱ；

(4.3)、更新外部存储器；

评价uⁱ的目标值，并归一化各目标值；令ct＝ct+1；将uⁱ加入外部存储器Arc；并删除当前Arc中所有重复解和Pareto支配解；设Chipop＝Chipop∪uⁱ；如果i＜N，则令i＝i+1，转至(4.1)，否则执行步骤(5)；

步骤(5)、解的更新：

设混合群体Mixpop＝{x¹,…,x^N}∪Chipop，令i＝1，

(5.1)、令计数器c＝0；

(5.2)、从第i个子问题的更新邻域P(i)中随机选取一个子群体的序号k，k∈P(i)；

(5.3)、对于第k个子问题，从Mixpop中确定最佳解；

设标记变量mark＝0；给定第k个子问题的权向量 为第j个目标的权值， 和参考向量 为第j个目标的参考点，j＝1,2,3，采用切比雪夫法求得任意一个个体x在第k个子问题的合成目标函数为：

由于采用归一化目标值，设置参考向量z^*＝(0,0,0)；对于Mixpop中的每个个体xy^l和第k个子问题的当前解x^k，如果下述3个条件中有一个满足：(i)g^te(xy^l|λ^k,z^*)＜g^te(x^k|λ^k,z^*)；(ii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l Pareto支配x^k；(iii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l在鲁棒性能f₃＝robustness上的目标值小于x^k，则令x^k＝xy^l，FV^k＝F(xy^l)，Fno^k＝Fno(xy^l)，且mark＝1；其中，FV^k表示x^k的目标向量，即FV^k＝F(x^k)＝[f₁(x^k),f₂(x^k),f₃(x^k)]，F(xy^l)表示xy^l的目标向量，即F(xy^l)＝[f₁(xy^l),f₂(xy^l),f₃(xy^l)]，Fno^k表示x^k经归一化后的目标向量，即Fno^k＝[fno₁(x^k),fno₂(x^k),fno₃(x^k)]，Fno(xy^l)表示xy^l经归一化后的目标向量，即Fno(xy^l)＝[fno₁(xy^l),fno₂(xy^l),fno₃(xy^l)]；

(5.4)、如果mark＝1，则令c＝c+1；

(5.5)、从P(i)中删除序号k；如果c＜n_r且P(i)非空集，则转至(5.2)；否则，如果i＜N，则令i＝i+1，转至(5.1)，否则转至步骤(6)；

步骤(6)、终止准则判断：

如果ct＞NmbEvl，则终止迭代，输出外部存储器Arc，即一组Pareto非支配的柔性作业车间调度解；否则，令t＝t+1，转至步骤(4)。

作为本发明所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法进一步优化方案，步骤(1)中所述优化目标中的作业完工时间表示完成柔性作业车间中的所有作业所花费的时间开销，定义为：

其中，C_v和S_v分别表示在第v项作业中，最后一道工序的加工完成时间和第一道工序的加工开始时间，v＝1,2,…,n，n为柔性作业车间中所有作业的总数；I表示初始景象，即针对不确定属性，设它的值为某一预估值确定不变，在此情况下计算作业完工时间；

步骤(1)所述优化目标中的机器最大负载表示各台机器加工时间的最大值，定义为：

其中，m是柔性作业车间内机器的总台数；是根据调度方案，在第s台机器上按第r个顺序进行加工的工序O^sr的加工时间；w_s是分配到第s台机器上的工序个数；

步骤(1)所述优化目标中，调度性能的鲁棒性是衡量调度方案的作业完工时间和机器的最大负载对不确定因素的敏感程度，定义为：

鲁棒性采用基于景象的方法定义，将一个调度方案在不确定属性的多种采样值{θ_q|q＝1,2,…,Q}下进行仿真，以比较作业完工时间和机器最大负载的实际值与初始预估值之间的差值；其中，makespan_q和workload_q分别是θ_q下相应的作业完工时间和机器最大负载目标值；

步骤(1)所述的工序优先级约束是指每项作业的各道工序是按事先确定的顺序进行加工的；在柔性作业车间问题中，每道工序可以在它允许的机器集合中任一台上进行加工；

步骤(1)所述的禁止抢先占有的约束包括：(i)每道工序的加工，只能在同一作业中排在它之前的所有工序都完成之后才能开始进行；(ii)如果将一道工序分配给了某台机器，只有在这台机器完成之前调度的所有工序之后，才能开始该道工序的加工。

作为本发明所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法进一步优化方案，步骤(3.2)所述的随机产生初始群体是指，每个个体中的工序序列向量通过随机排列所有作业中的所有工序生成；对于机器分配向量，将每道作业工序随机分配到它机器集合中的任一台上进行加工。

作为本发明所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法进一步优化方案，步骤(4.3)和步骤(5.3)所述的Pareto支配是指：设x₁和x₂为多目标优化问题的两个解，问题的目标个数为m₁，假设所有目标均需最小化，称x₁Pareto支配x₂当且仅当

且

其中，g和h分别表示目标的某一个序号，f_g(x₁)和f_g(x₂)分别表示x₁和x₂在第g个目标f_g上的目标值，f_h(x₁)和f_h(x₂)分别表示x₁和x₂在第h个目标f_h上的目标值。

作为本发明所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法进一步优化方案，步骤(3)和步骤(6)所述的Pareto非支配解是指，如果在某一集合中不存在任何其它解x'Pareto支配解x，则称x为该集合中的Pareto非支配解。

作为本发明所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法进一步优化方案，步骤(4.2)所述的繁殖算子是指，基于微分进化算法中的变异算子和交叉算子，设计一种改进的自适应变异算子和基于修复的交叉算子，由父代个体和第i个子问题的当前个体xⁱ，生成子代个体uⁱ；其中，改进的自适应变异算子的实施方法如下：

β＝e^-0.015t,

其中，F₁，F₂和F₃是变异因子，它们的值分别从[0.5,1]中均匀随机生成；是当前外部存储器Arc中，距离xⁱ最近的解，参数β的值随着进化代数t的取值而变化；yⁱ是由改进的自适应变异算子生成的个体；

对于个体的工序序列向量和机器分配向量，分别实施上述所述改进的自适应变异算子；对于生成的yⁱ中的工序序列向量ysⁱ，将它的各个元素按从小到大的顺序进行排列，得到ysⁱ的排序向量tysⁱ，并将排序后的每个元素在原向量ysⁱ中的位置索引记录在索引向量Iⁱ中；将xⁱ的工序序列向量xsⁱ中各个元素按从小到大的顺序进行排列，得到xsⁱ的排序向量txsⁱ，并令ysⁱ(Iⁱ)＝txsⁱ，ysⁱ(Iⁱ)表示ysⁱ中所有在Iⁱ记录的位置上的元素；对于生成的yⁱ中的机器分配向量ymⁱ，针对每个元素，搜索其对应工序的机器集合，从中确定出与该元素值最接近的机器，并将该机器替代此元素；如果有两台以上机器同时满足此条件，则从中随机选择一台机器替换对应元素；更新后的工序序列向量ysⁱ和机器分配向量ymⁱ构成新个体yⁱ；

所述基于修复的交叉算子的实施方法如下：

(a)、将由改进的自适应变异算子生成的个体yⁱ与xⁱ组合，构成子代个体uⁱ：

其中，分别是uⁱ,yⁱ,xⁱ的第ll个元素，L是个体xⁱ的长度，CR∈[0,1]是交叉概率，是从[0,1]中均匀随机产生的一个数，Irandⁱ是从集合{1,2,…,L}中随机选择的一个整数；

在机器分配向量的表示方法中，yⁱ的机器分配向量ymⁱ与xⁱ的机器分配向量xmⁱ在相同位置上的机器对应于同一道工序；上述交叉算子直接作用于ymⁱ和xmⁱ,交叉得到的结果即为子代个体uⁱ的机器分配向量umⁱ；

针对yⁱ的工序序列向量ysⁱ与xⁱ的工序序列向量xsⁱ，进一步实施以下步骤：

(b)、将步骤(a)求得的uⁱ的工序序列向量usⁱ中，将来自于xsⁱ的元素位置索引记录在1号索引向量index₁中，将来自于ysⁱ的元素位置索引记录在2号索引向量index₂中；令usⁱ＝xsⁱ，临时索引向量tempindex₂＝index₂,测试向量test＝usⁱ(index₁)，usⁱ(index₁)表示usⁱ中所有在index₁记录的位置上的元素,索引删除向量donor_delete＝[]，[]表示空集合；

(c)、对test中的每个元素ue，确定在ysⁱ(index₁)中是否存在一些元素等于ue；如果存在，则从中均匀随机地选择一个元素，将选中的元素在ysⁱ中的位置索引加入donor_delete，并将该索引从index₁中删除；否则，找出在ysⁱ(index₂)中等于ue的元素，从中均匀随机地选取一个，将它在ysⁱ中的位置索引加入donor_delete，并将该索引从index₂中删除；其中，ysⁱ(index₁)表示ysⁱ中所有在index₁记录的位置上的元素，ysⁱ(index₂)表示ysⁱ中所有在index₂记录的位置上的元素；

(d)、将donor_delete中的元素从集合{1,2,…,L}中移除，即令索引保留向量donar_reserve＝{1,2,…,L}\donar_delete，且usⁱ(tempindex₂)＝ysⁱ(donar_reserve)；其中，\表示去除，usⁱ(tempindex₂)表示usⁱ中所有在tempindex₂记录的位置上的元素，ysⁱ(donar_reserve)表示ysⁱ中所有在donar_reserve记录的位置上的元素。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1)本发明能够处理柔性作业车间生产环境中存在的不确定因素，在保证完工时间较短、机器负载较为均衡等性能的同时，增强车间调度方案对不确定性的抗干扰能力，降低调度性能对扰动的敏感程度；它能够同时处理工序排序和机器分配这两种调度策略，因此，与现有技术相比，本发明更适合处理现实世界中存在不确定因素的柔性作业车间调度问题。本发明提出的鲁棒调度技术对实际生产***的成功实施具有重要意义；

2)本发明同时优化了柔性作业车间的效率指标(完工时间、机器最大负载)和鲁棒性，并采用基于分解的多目标进化算法对多个目标并行处理，从而能够为生产管理者提供一组反映目标间不同折中程度的调度方案，为其做出最终的决策提供有力的参考；

3)本发明设计了一种新型的，能够有效利用全局信息的子问题更新策略；设置外部存储器保存Pareto非支配个体，并让这些精英个体有机会参与到子代个体的生成中；采用一种改进的自适应变异算子和基于修复的交叉算子进行繁殖操作，以助于更好地维护算法在“探索”和“利用”之间的平衡；这些机制能够有效地提高本发明的局部搜索能力，避免陷入局部最优，快速高效地实现柔性作业车间中的调度任务。

附图说明

图1是本发明的主体流程图。

图2为个体的表示方法示例图。

图3为调度方案解Solution_a在初始景象下的甘特图。

图4为调度方案解Solution_a在扰动情况下的甘特图。

图5为调度方案解Solution_b在初始景象下的甘特图。

图6为调度方案解Solution_b在扰动情况下的甘特图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

一个柔性作业车间中，有5台机器，4项待加工的作业。每项作业的每道工序均可以在5台机器上分别加工。这些工序中，有部分工序的加工时间存在不确定性。4项作业包含的工序个数、每道工序在各台允许加工机器上的初始预估加工时间如表1所示。

表1

	O11	O12	O13	O21
					预估加工时间	2,5,4,1,2	5,4,5,7,5	4,5,5,4,5	2,5,4,7,8
	O22	O23	O31	O32
					预估加工时间	5,6,9,8,5	4,5,4,54,5	9,8,6,7,9	6,1,2,5,4
	O33	O34	O41	O42
					预估加工时间	2,5,4,2,4	4,5,2,1,5	1,5,2,4,12	5,1,2,1,2

使用本发明提出的基于分解多目标进化算法的柔性作业车间鲁棒调度方法求解该柔性作业车间实施例的调度方案，主体流程图如图1所示，具体步骤如下：

(1)初始化。读取柔性作业车间的输入信息，包括每项作业的工序数、每项作业中每道工序可分配的机器集合、各工序在相应机器上的初始预估加工时间(见表1)；定义优化目标，设定约束条件：

优化目标中的“作业完工时间”表示完成柔性作业车间中的所有作业所花费的时间开销，它定义为：

其中，C_v和S_v分别表示第v项作业J_v中，最后一道工序的加工完成时间和第一道工序的加工开始时间，v＝1,2,…,n；n为柔性作业车间中所有作业的总数，本实例中，n＝4；I表示初始景象，本实例中，针对加工时间这一不确定属性，假设它的值为初始预估值确定不变，在此情况下计算初始景象下的作业完工时间；

优化目标中的“机器最大负载”表示各台机器总的加工时间的最大值，它定义为：

其中，m是柔性作业车间内机器的总台数，本实例中，m＝5；是根据当前调度方案，在第s台机器上按第r个顺序进行加工的工序O^sr的加工时间；w_s是分配到第s台机器上的工序个数；I表示初始景象；

优化目标中“调度性能的鲁棒性”衡量调度方案的作业完工时间和机器的最大负载对不确定因素的敏感程度，它定义为：

鲁棒性采用基于景象的方法定义，将一个调度方案在不确定工序加工时间的多种采样值{θ_q|q＝1,2,…,Q}下进行仿真，以比较作业完工时间和机器最大负载的实际值与初始预估值之间的差值；其中，θ_q是不确定工序加工时间的第q个采样值，Q是样本个数；makespan_q和workload_q分别是采样值θ_q下相应的作业完工时间和机器最大负载目标值；

工序优先级约束是指每项作业的各道工序是按事先确定的顺序进行加工的；在柔性作业车间问题中，每道工序可以在它允许的机器集合中任一台上进行加工；

禁止抢先占有的约束包括：(i)每道工序的加工，只能在同一作业中排在它之前的所有工序都完成之后才能开始进行；(ii)如果将一道工序分配给了某台机器，只有在这台机器完成之前调度的所有工序之后，才能开始该道工序的加工；

(2)初始化基于分解的多目标进化算法的参数：

设置目标评价次数NmbEvl为15000、群体规模(即子问题的个数)N为500、产生N个均匀分布的目标权向量λ¹,…,λ^N、子问题的邻域规模T为50，繁殖算子中的父代个体从邻域里选取的概率δ为0.7、邻域个体允许被每个子代个体替换的最大数目n_r为5，评价鲁棒性时采样的不确定场景个数Q为30，交叉概率CR为0.3；

(3)确定每个子问题的邻域，产生初始父代群体：

(3.1)对第i个子问题，i＝1,…,N，确定邻域B(i)＝{i₁,…,i_T}。分别计算第i个子问题的权向量λⁱ和其它(N-1)个子问题权向量间的欧拉距离，将T个距离λⁱ最近的权向量所对应的子群体序号{i₁,…,i_T}构成B(i)；

(3.2)随机产生初始父代群体{x¹,…,x^N}，其中个体xⁱ为第i个子群体的当前解；群体中的每个个体包括工序序列向量和机器分配向量；随机采样一组不确定性场景θ_q，q＝1,2,...,Q；计算初始父代群体中每个个体的目标值f₁＝makespan_I、f₂＝workload_I和f₃＝robustness；从初始父代群体中确定出所有的Pareto非支配解构成外部存储器Arc；设置目标评价次数计数变量ct＝N；归一化各目标值，将每个个体xⁱ的第j个目标值f_j(xⁱ)映射到区间[0,1]之间，产生xⁱ在第j个目标上的归一化值fno_j(xⁱ)，即：

其中，和分别表示当前父代群体中的所有个体在第j个目标上的最大值和最小值，j＝1,2,3。

(3.3)设置进化代数t＝0；

初始父代群体由N个随机生成的个体组成。本发明中，一个个体包括两个向量：(i)工序序列向量；(ii)机器分配向量。图2给出了个体表示方法的示例。对于工序序列向量，采用基于作业的表示方法，同一作业的工序均用该作业号表示。例如，图2中，工序O₁₁、O₁₂、O₁₃均用作业号1表示。每道工序依据它在工序序列向量中出现的顺序进行转换。例如，工序序列向量中第一个出现的数字3代表工序O₃₁,第二个出现的3代表O₃₂，依此类推。因此，可将图2中的工序序列向量解释为：

其中，表示将工序a首先加入它所分配机器的等待队列，然后再调度工序b。

机器分配向量表示给每道工序分配的机器。它的分配顺序是：从当前序号最小作业的第一道工序到序号最大作业的最后一道工序。例如，在图2的机器分配向量中，第一个元素2表示将作业1的第一道工序O₁₁分配给机器2，第二个元素3表示将作业1的第二道工序O₁₂分配给机器3，依此类推，机器分配向量可解释为：

O₁₁→机器2,O₁₂→机器3,O₁₃→机器5,O₂₁→机器4,

O₃₁→机器5,O₃₂→机器3,O₄₁→机器1,O₄₂→机器2

其中，→表示将工序分配给相应的机器。

随机产生初始群体是指，每个个体中的工序序列向量通过随机排列所有作业中的所有工序生成；对于机器分配向量，将每道作业工序随机分配到它机器集合中的任一台上进行加工；

Pareto支配是指：设x₁和x₂为多目标优化问题的两个解，问题的目标个数为m₁，假设所有目标均需最小化，称x₁Pareto支配x₂当且仅当

其中，g和h分别表示目标的某一个序号，f_g(x₁)和f_g(x₂)分别表示x₁和x₂在第g个目标f_g上的目标值，f_h(x₁)和f_h(x₂)分别表示x₁和x₂在第h个目标f_h上的目标值。

Pareto非支配解是指，如果在某一集合中不存在任何其它解x'Pareto支配x，则称x为该集合中的Pareto非支配解(简称非支配解)。提高Pareto非支配解在任何一个目标上的性能，都必然会导致它在剩余的至少一个目标上的性能降低；

(4)生成子代群体：

随机采样一组不确定性场景θ_q；设置子代群体令i＝1；

(4.1)交配选择；

产生一个均匀分布的随机数rand₁∈[0,1]；设置第i个子问题的更新邻域P(i)：

基于以下规则产生3个不同的父代个体和

(4.2)繁殖；

采用基于微分进化算法的繁殖算子，由和xⁱ生成子代个体uⁱ；

(4.3)更新外部存储器；

评价uⁱ的目标值，并归一化各目标值；令ct＝ct+1；将uⁱ加入外部存储器Arc；并删除当前Arc中所有重复解和Pareto支配解；设Chipop＝Chipop∪uⁱ；如果i＜N，则令i＝i+1，转至(4.1)，否则执行步骤(5)；

步骤(4.2)所述的繁殖算子是指，基于微分进化算法中的变异算子和交叉算子，设计一种改进的自适应变异算子和基于修复的交叉算子，由父代个体和第i个子问题的当前个体xⁱ，生成子代个体uⁱ；其中，改进的自适应变异算子的实施方法如下：

β＝e^-0.015t,

其中，F₁，F₂和F₃是变异因子，它们的值分别从[0.5,1]中均匀随机生成；是当前外部存储器Arc中，距离xⁱ最近的解，参数β的值随着进化代数t的取值而变化；yⁱ是由改进的自适应变异算子生成的个体；

对于个体的工序序列向量和机器分配向量，分别实施上述所述改进的自适应变异算子；在实施完成后，生成向量中的部分元素可能为非整数，即不可行值。为了解决该问题,对于生成的yⁱ中的工序序列向量ysⁱ，将它的各个元素按从小到大的顺序进行排列，得到ysⁱ的排序向量tysⁱ，并将排序后的每个元素在原向量ysⁱ中的位置索引记录在索引向量Iⁱ中；将xⁱ的工序序列向量xsⁱ中各个元素按从小到大的顺序进行排列，得到xsⁱ的排序向量txsⁱ，并令ysⁱ(Iⁱ)＝txsⁱ(ysⁱ(Iⁱ)表示ysⁱ中所有在Iⁱ记录的位置上的元素)；对于生成的yⁱ中的机器分配向量ymⁱ，针对每个元素，搜索其对应工序的机器集合，从中确定出与该元素值最接近的机器，并将该机器替代此元素；如果有两台以上机器同时满足此条件，则从中随机选择一台机器替换对应元素；更新后的工序序列向量ysⁱ和机器分配向量ymⁱ构成个体yⁱ；

所述基于修复的交叉算子的实施方法如下：

(a)将由改进的自适应变异算子生成的个体yⁱ与xⁱ组合，构成子代个体uⁱ：

其中，分别是uⁱ,yⁱ,xⁱ的第ll个元素，L是个体xⁱ的长度，CR∈[0,1]是交叉概率，是从[0,1]中均匀随机产生的一个数，Irandⁱ是从集合{1,2,…,L}中随机选择的一个整数；

在本发明机器分配向量的表示方法中，yⁱ的机器分配向量ymⁱ与xⁱ的机器分配向量xmⁱ在相同位置上的机器对应于同一道工序。因此，上述交叉算子可以直接作用于ymⁱ和xmⁱ,交叉得到的结果即为子代个体uⁱ的机器分配向量umⁱ。

对于工序序列向量，上述交叉算子的实施可能产生不可行结果，即在交叉后的工序序列向量中，出现冗余的作业工序，而丢失了另外一些工序。为了解决该问题，本发明针对yⁱ的工序序列向量ysⁱ与xⁱ的工序序列向量xsⁱ，需进一步实施以下步骤：

(b)将步骤(a)求得的uⁱ的工序序列向量usⁱ中，来自于xsⁱ的元素位置索引记录在1号索引向量index₁中，将来自于ysⁱ的元素位置索引记录在2号索引向量index₂中；令usⁱ＝xsⁱ，临时索引向量tempindex₂＝index₂,测试向量test＝usⁱ(index₁)(usⁱ(index₁)表示usⁱ中所有在index₁记录的位置上的元素),索引删除向量donor_delete＝[]，[]表示空集合；

(c)对test中的每个元素ue，确定在ysⁱ(index₁)(ysⁱ(index₁)表示ysⁱ中所有在index₁记录的位置上的元素)中是否存在一些元素等于ue；如果存在，则从中均匀随机地选择一个元素，将选中的元素在ysⁱ中的位置索引加入donor_delete，并将该索引从index₁中删除；否则，找出在ysⁱ(index₂)(ysⁱ(index₂)表示ysⁱ中所有在index₂记录的位置上的元素)中等于ue的元素，从中均匀随机地选取一个，将它在ysⁱ中的位置索引加入donor_delete，并将该索引从index₂中删除；

(d)将donor_delete中的元素从集合{1,2,…,L}中移除，即令索引保留向量donar_reserve＝{1,2,…,L}\donar_delete(\表示去除)，且usⁱ(tempindex₂)＝ysⁱ(donar_reserve)(usⁱ(tempindex₂)表示usⁱ中所有在tempindex₂记录的位置上的元素，ysⁱ(donar_reserve)表示ysⁱ中所有在donar_reserve记录的位置上的元素)。

(5)解的更新：

设混合群体Mixpop＝{x¹,…,x^N}∪Chipop，令i＝1，

(5.1)令计数器c＝0；

(5.2)从第i个子问题的更新邻域P(i)中随机选取一个子群体的序号k，k∈P(i)；

(5.3)对于第k个子问题，从Mixpop中确定最佳解；

设标记变量mark＝0；给定第k个子问题的权向量(为第j个目标的权值，)和参考向量 为第j个目标的参考点，采用切比雪夫法求得任意一个个体x在第k个子问题的合成目标函数为：

由于采用归一化目标值，设置参考向量z^*＝(0,0,0)；对于Mixpop中的每个个体xy^l和第k个子问题的当前解x^k，如果下述3个条件中有一个满足：(i)g^te(xy^l|λ^k,z^*)＜g^te(x^k|λ^k,z^*)；(ii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l Pareto支配x^k；(iii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l在鲁棒性能f₃＝robustness上的目标值小于x^k，则令x^k＝xy^l，FV^k＝F(xy^l)，Fno^k＝Fno(xy^l)，且mark＝1；其中，FV^k表示x^k的目标向量，即FV^k＝F(x^k)＝[f₁(x^k),f₂(x^k),f₃(x^k)]，F(xy^l)表示xy^l的目标向量，即F(xy^l)＝[f₁(xy^l),f₂(xy^l),f₃(xy^l)]，Fno^k表示x^k经归一化后的目标向量，即Fno^k＝[fno₁(x^k),fno₂(x^k),fno₃(x^k)]，Fno(xy^l)表示xy^l经归一化后的目标向量，即Fno(xy^l)＝[fno₁(xy^l),fno₂(xy^l),fno₃(xy^l)]；

(5.4)如果mark＝1，则令c＝c+1；

(5.5)从P(i)中删除序号k；如果c＜n_r且P(i)非空集，则转至(5.2)；否则，如果i＜N，则令i＝i+1，转至(5.1)，否则转至步骤(6)；

步骤(6)终止准则判断：

如果ct＞NmbEvl，则终止迭代，输出外部存储器Arc，即一组Pareto非支配的柔性作业车间调度解，将这些解提供给生产管理者进行参考，并由他从中挑选出一个符合生产要求的满意解作为调度方案；否则，令t＝t+1，转至步骤(4)。

本发明的效果可以通过以下仿真实验进一步说明：

1.实验条件：

在CPU为Intel core Duo 2.2GHz、内存4GB、WINDOWS 7***上使用Matlab 2010进行仿真。

2.实验内容：

本发明针对上述具有5台机器，4项作业的柔性作业车间实施例求解生产调度方案。本实施例中，每项作业的每道工序均可以在5台机器上分别加工。这些工序中，有部分工序的加工时间存在不确定性。4项作业包含的工序个数、每道工序在各台允许加工机器上的初始预估加工时间、以及受扰动后的实际加工时间如表1所示。

3.实验结果

本发明运行一次，可以求得一组Pareto非支配的柔性作业车间调度解。从这组解中挑选2个调度方案解Solution_a和Solution_b，当面临相同的加工时间不确定性时，给出它们在效率性能(作业完工时间、机器最大负载)和鲁棒性能上的比较。每道工序在各台允许加工机器上的初始预估加工时间、以及受扰动后的实际加工时间如表2所示。两个解在初始预估景象下的作业完工时间makespan_I、机器最大负载workload_I、受扰动情况下的作业完工时间makespan_q、机器最大负载workload_q，以及它们的鲁棒目标值如表3所示。从表3可以看出，为了获取更好的鲁棒性，Solution_b的初始作业完工时间makespan_I和机器最大负载workload_I劣于Solution_a。然而，当面临相同的加工时间扰动时，Solution_a的作业完工时间makespan_q和机器最大负载workload_q却不仅劣于Solution_b在扰动情况下的相应值，也劣于Solution_b在初始景象中的makespan_I和workload_I。由此说明，Solution_b所具有的优良鲁棒性，补偿了它在初始景象中效率性能(makespan_I和workload_I)的弱势，使得它具有更强的抗干扰能力，降低了它的效率性能对扰动的敏感程度。图3为调度方案解Solution_a在初始景象中的甘特图，图4为调度方案解Solution_a在扰动情况下的甘特图，图5为调度方案解Solution_b在初始景象中的甘特图，图6为调度方案解Solution_b在扰动情况下的甘特图。从甘特图中可以获取为每道工序分配的加工机器，每道工序的开始加工时间和结束时间、每项作业的开始加工时间和结束时间，以及每台机器上各工序的先后加工顺序。

表2

每个单元格中每行的5个值分别为各工序在各台机器上的预估或扰动加工时间

表3

分别使用本发明方法，以及两个现有的经典多目标进化算法MOEA/D-DE和NSGA-II求解本实施例，将求得的Pareto非支配解集在收敛性能和分布性能上进行比较。收敛性测度采用超体积比率(hypervolume ratio，简称HVR)度量。HVR的值越大，说明算法求得的Pareto非支配解集在目标空间上的收敛性越好，分布越广。分布性能用测度Spread度量，Spread越小，说明算法求得的Pareto非支配解集的分布越宽广且越均匀。将本发明方法和其它两种已有算法分别在本实施例中运行30次，采用显著性水平为0.5的Wilcoxon秩和检验法对3种算法在30次运行中的性能进行统计测试，结果如表4所示。表4中，p值是Wilcoxon秩和检验法的返回值，假设两组样本的分布具有相同的中值，则p值表示样本数据支持该假设的证据,p值越大，证据越强。符号‘+/-/＝’分别表示，在算法A vs.算法B中，依据所采用的显著性水平为0.5的Wilcoxon秩和检验法，算法A显著优于B，或算法A显著劣于B，或算法A和B之间没有显著差别。由表4可见，在本实施例中，对于收敛性测度HVR和分布性能用测度Spread，本发明均显著优于MOEA/D-DE和NSGA-II，说明与两种已有的方法相比，本发明能够为生产管理者提供一组分布更宽广，且收敛性更好的Pareto非支配解。

表4

综上，本发明提出的基于分解多目标进化算法的柔性作业车间鲁棒调度方法，由于采用了新型的子问题更新策略，利用外部存储器中的精英个体参与子代个体生成，以及使用改进的自适应变异算子和基于修复的交叉算子进行繁殖，从而可以更好地利用全局信息，并维护算法在“探索”和“利用”之间的平衡，提高了算法的收敛特性和分布性能。这些机制有效地提高了本发明的搜索能力，克服了传统的生产调度方法局部搜索能力较弱、易于陷入局部最优、调度效率低下的缺点，能够快速高效地实现柔性作业车间中的调度任务。此外，通过引入鲁棒性能指标，本发明能够处理柔性作业车间生产环境中存在的不确定因素，在保证完工时间较短、机器负载较为均衡等性能的同时，增强车间调度方案对不确定性的抗干扰能力。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (6)

1.一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，包括以下步骤：

步骤(1)、读取柔性作业车间的输入信息，定义优化目标，设定约束条件：

车间的输入信息包括每项作业的工序数、每项作业中每道工序可分配的机器集合、各工序在相应机器上的加工时间、正常工作的机器数目；优化目标包括作业完工时间、机器的最大负载、调度性能的鲁棒性；约束条件包括工序优先级约束和禁止抢先占有的约束；

步骤(2)、初始化基于分解的多目标进化算法的参数：

设置迭代次数NmbEvl、群体规模N即为子问题的个数、产生N个均匀分布的目标权向量λ¹,…,λ^N、子问题的邻域规模T，繁殖算子中的父代个体从邻域里选取的概率δ、邻域个体允许被每个子代个体替换的最大数目n_r，评价鲁棒性时采样的不确定场景个数Q，以及交叉概率CR；

步骤(3)、确定每个子问题的邻域，产生初始父代群体：

(3.1)、对第i个子问题，i＝1,…,N，确定邻域B(i)＝{i₁,…,i_T}；分别计算第i个子问题的权向量λⁱ和其它(N-1)个子问题权向量间的欧拉距离，将T个距离λⁱ最近的权向量所对应的子群体序号{i₁,…,i_T}构成B(i)；

(3.2)、随机产生初始父代群体{x¹,…,x^N}，其中，个体xⁱ为第i个子群体的当前解；群体中的每个个体包括工序序列向量和机器分配向量；随机采样一组不确定性场景θ_q，q＝1,2,...,Q；计算初始父代群体中每个个体的目标值f₁＝makespan_I、f₂＝workload_I和f₃＝robustness；makespan_I和workload_I分别表示初始景象下的作业完工时间和机器最大负载，robustness表示调度性能的鲁棒性；从初始父代群体中确定出所有的Pareto非支配解构成外部存储器Arc；设置目标评价次数计数变量ct＝N；归一化各目标值，将每个个体xⁱ的第j个目标值f_j(xⁱ)映射到区间[0,1]之间，产生xⁱ在第j个目标上的归一化值fno_j(xⁱ)，即：

其中，和分别表示当前父代群体中的所有个体在第j个目标上的最大值和最小值，j＝1,2,3；

(3.3)、设置进化代数t＝0；

步骤(4)、生成子代群体：

随机采样一组不确定性场景θ_q；设置子代群体令i＝1；

(4.1)、交配选择；

产生一个均匀分布的随机数rand₁∈[0,1]；设置第i个子问题的更新邻域P(i)：

基于以下规则产生3个不同的父代个体和

(4.2)、繁殖；

采用基于微分进化算法的繁殖算子，由和xⁱ生成子代个体uⁱ；

(4.3)、更新外部存储器；

评价uⁱ的目标值，并归一化各目标值；令ct＝ct+1；将uⁱ加入外部存储器Arc；并删除当前Arc中所有重复解和Pareto支配解；设Chipop＝Chipop∪uⁱ；如果i＜N，则令i＝i+1，转至(4.1)，否则执行步骤(5)；

步骤(5)、解的更新：

设混合群体Mixpop＝{x¹,…,x^N}∪Chipop，令i＝1，

(5.1)、令计数器c＝0；

(5.2)、从第i个子问题的更新邻域P(i)中随机选取一个子群体的序号k，k∈P(i)；

(5.3)、对于第k个子问题，从Mixpop中确定最佳解；

设标记变量mark＝0；给定第k个子问题的权向量 为第j个目标的权值，和参考向量 为第j个目标的参考点，j＝1,2,3，采用切比雪夫法求得任意一个个体x在第k个子问题的合成目标函数为：

由于采用归一化目标值，设置参考向量z^*＝(0,0,0)；对于Mixpop中的每个个体xy^l和第k个子问题的当前解x^k，如果下述3个条件中有一个满足：(i)g^te(xy^l|λ^k,z^*)＜g^te(x^k|λ^k,z^*)；(ii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l Pareto支配x^k；(iii)g^te(xy^l|λ^k,z^*)＝g^te(x^k|λ^k,z^*)且xy^l在鲁棒性能f₃＝robustness上的目标值小于x^k，则令x^k＝xy^l，FV^k＝F(xy^l)，Fno^k＝Fno(xy^l)，且mark＝1；其中，FV^k表示x^k的目标向量，即FV^k＝F(x^k)＝[f₁(x^k),f₂(x^k),f₃(x^k)]，F(xy^l)表示xy^l的目标向量，即F(xy^l)＝[f₁(xy^l),f₂(xy^l),f₃(xy^l)]，Fno^k表示x^k经归一化后的目标向量，即Fno^k＝[fno₁(x^k),fno₂(x^k),fno₃(x^k)]，Fno(xy^l)表示xy^l经归一化后的目标向量，即Fno(xy^l)＝[fno₁(xy^l),fno₂(xy^l),fno₃(xy^l)]；

(5.4)、如果mark＝1，则令c＝c+1；

(5.5)、从P(i)中删除序号k；如果c＜n_r且P(i)非空集，则转至(5.2)；否则，如果i＜N，则令i＝i+1，转至(5.1)，否则转至步骤(6)；

步骤(6)、终止准则判断：

如果ct＞NmbEvl，则终止迭代，输出外部存储器Arc，即一组Pareto非支配的柔性作业车间调度解；否则，令t＝t+1，转至步骤(4)。

2.根据权利要求1所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，步骤(1)中所述优化目标中的作业完工时间表示完成柔性作业车间中的所有作业所花费的时间开销，定义为：

其中，C_v和S_v分别表示在第v项作业中，最后一道工序的加工完成时间和第一道工序的加工开始时间，v＝1,2,…,n，n为柔性作业车间中所有作业的总数；I表示初始景象，即针对不确定属性，设它的值为某一预估值确定不变，在此情况下计算作业完工时间；

步骤(1)所述优化目标中的机器最大负载表示各台机器加工时间的最大值，定义为：

其中，m是柔性作业车间内机器的总台数；是根据调度方案，在第s台机器上按第r个顺序进行加工的工序O^sr的加工时间；w_s是分配到第s台机器上的工序个数；

步骤(1)所述优化目标中，调度性能的鲁棒性是衡量调度方案的作业完工时间和机器的最大负载对不确定因素的敏感程度，定义为：

鲁棒性采用基于景象的方法定义，将一个调度方案在不确定属性的多种采样值{θ_q|q＝1,2,…,Q}下进行仿真，以比较作业完工时间和机器最大负载的实际值与初始预估值之间的差值；其中，makespan_q和workload_q分别是θ_q下相应的作业完工时间和机器最大负载目标值；

步骤(1)所述的工序优先级约束是指每项作业的各道工序是按事先确定的顺序进行加工的；在柔性作业车间问题中，每道工序可以在它允许的机器集合中任一台上进行加工；

步骤(1)所述的禁止抢先占有的约束包括：(i)每道工序的加工，只能在同一作业中排在它之前的所有工序都完成之后才能开始进行；(ii)如果将一道工序分配给了某台机器，只有在这台机器完成之前调度的所有工序之后，才能开始该道工序的加工。

3.根据权利要求1所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，步骤(3.2)所述的随机产生初始群体是指，每个个体中的工序序列向量通过随机排列所有作业中的所有工序生成；对于机器分配向量，将每道作业工序随机分配到它机器集合中的任一台上进行加工。

4.根据权利要求1所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，步骤(4.3)和步骤(5.3)所述的Pareto支配是指：设x₁和x₂为多目标优化问题的两个解，问题的目标个数为m₁，假设所有目标均需最小化，称x₁Pareto支配x₂当且仅当

且

其中，g和h分别表示目标的某一个序号，f_g(x₁)和f_g(x₂)分别表示x₁和x₂在第g个目标f_g上的目标值，f_h(x₁)和f_h(x₂)分别表示x₁和x₂在第h个目标f_h上的目标值。

5.根据权利要求1所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，步骤(3)和步骤(6)所述的Pareto非支配解是指，如果在某一集合中不存在任何其它解x'Pareto支配解x，则称x为该集合中的Pareto非支配解。

6.根据权利要求1所述的一种基于分解多目标进化算法的柔性车间鲁棒调度方法，其特征在于，步骤(4.2)所述的繁殖算子是指，基于微分进化算法中的变异算子和交叉算子，设计一种改进的自适应变异算子和基于修复的交叉算子，由父代个体和第i个子问题的当前个体xⁱ，生成子代个体uⁱ；其中，改进的自适应变异算子的实施方法如下：

β＝e^-0.015t,

其中，F₁，F₂和F₃是变异因子，它们的值分别从[0.5,1]中均匀随机生成；是当前外部存储器Arc中，距离xⁱ最近的解，参数β的值随着进化代数t的取值而变化；yⁱ是由改进的自适应变异算子生成的个体；

对于个体的工序序列向量和机器分配向量，分别实施上述所述改进的自适应变异算子；对于生成的yⁱ中的工序序列向量ysⁱ，将它的各个元素按从小到大的顺序进行排列，得到ysⁱ的排序向量tysⁱ，并将排序后的每个元素在原向量ysⁱ中的位置索引记录在索引向量Iⁱ中；将xⁱ的工序序列向量xsⁱ中各个元素按从小到大的顺序进行排列，得到xsⁱ的排序向量txsⁱ，并令ysⁱ(Iⁱ)＝txsⁱ，ysⁱ(Iⁱ)表示ysⁱ中所有在Iⁱ记录的位置上的元素；对于生成的yⁱ中的机器分配向量ymⁱ，针对每个元素，搜索其对应工序的机器集合，从中确定出与该元素值最接近的机器，并将该机器替代此元素；如果有两台以上机器同时满足此条件，则从中随机选择一台机器替换对应元素；更新后的工序序列向量ysⁱ和机器分配向量ymⁱ构成新个体yⁱ；

所述基于修复的交叉算子的实施方法如下：

(a)、将由改进的自适应变异算子生成的个体yⁱ与xⁱ组合，构成子代个体uⁱ：

其中，分别是uⁱ,yⁱ,xⁱ的第ll个元素，L是个体xⁱ的长度，CR∈[0,1]是交叉概率，是从[0,1]中均匀随机产生的一个数，Irandⁱ是从集合{1,2,…,L}中随机选择的一个整数；

在机器分配向量的表示方法中，yⁱ的机器分配向量ymⁱ与xⁱ的机器分配向量xmⁱ在相同位置上的机器对应于同一道工序；上述交叉算子直接作用于ymⁱ和xmⁱ,交叉得到的结果即为子代个体uⁱ的机器分配向量umⁱ；

针对yⁱ的工序序列向量ysⁱ与xⁱ的工序序列向量xsⁱ，进一步实施以下步骤：

(b)、将步骤(a)求得的uⁱ的工序序列向量usⁱ中，将来自于xsⁱ的元素位置索引记录在1号索引向量index₁中，将来自于ysⁱ的元素位置索引记录在2号索引向量index₂中；令usⁱ＝xsⁱ，临时索引向量tempindex₂＝index₂,测试向量test＝usⁱ(index₁)，usⁱ(index₁)表示usⁱ中所有在index₁记录的位置上的元素,索引删除向量donor_delete＝[]，[]表示空集合；

(c)、对test中的每个元素ue，确定在ysⁱ(index₁)中是否存在一些元素等于ue；如果存在，则从中均匀随机地选择一个元素，将选中的元素在ysⁱ中的位置索引加入donor_delete，并将该索引从index₁中删除；否则，找出在ysⁱ(index₂)中等于ue的元素，从中均匀随机地选取一个，将它在ysⁱ中的位置索引加入donor_delete，并将该索引从index₂中删除；其中，ysⁱ(index₁)表示ysⁱ中所有在index₁记录的位置上的元素，ysⁱ(index₂)表示ysⁱ中所有在index₂记录的位置上的元素；

(d)、将donor_delete中的元素从集合{1,2,…,L}中移除，即令索引保留向量donar_reserve＝{1,2,…,L}\donar_delete，且usⁱ(tempindex₂)＝ysⁱ(donar_reserve)；其中，\表示去除，usⁱ(tempindex₂)表示usⁱ中所有在tempindex₂记录的位置上的元素，ysⁱ(donar_reserve)表示ysⁱ中所有在donar_reserve记录的位置上的元素。