一种随机多产品机器人拆卸线平衡控制方法

本发明属于信息与控制技术领域，特别涉及一种随机多产品机器人拆卸线平衡控制方法。

背景技术：

随着工业的发展，大量的产品报废，造成了许多环境问题。回收利用是处理eol产品(废旧产品)的有效方法之一。eol产品回收再利用的价值体现在两个方面：1)有利于环境保护；2)节约资源，降低生产成本。拆卸是回收过程中的一个重要步骤。由于人类操作的灵活性，大多数拆卸操作仍然是手动进行的，但是拆卸成本很高。近年来，越来越多的研究集中在机器人辅助下的自动拆卸。它可以降低拆卸成本，提高拆卸效率，保证零件质量，防止人类接触危险零件。因此，有必要研究机器人拆卸线平衡(rdlb)问题。国内外一些学者对此进行了研究。例如，zeynelabidincil等人研究单产品的rdlb问题。刘佳宜等人解决了拆卸序列规划与rdlb的协同优化问题。然而，上述研究只集中在单个产品的拆卸上。在实际中，通常需要同时拆卸多个产品以满足客户的需求。因此，研究多产品拆卸方法具有重要意义。

为了解决该问题，研究者们探索了一些有效的方法。zeynelabidincil等人提出了一种区间搜索蚁群优化算法来求解rdlb问题，采用正反馈机制，使得搜索过程不断收敛，最终逼近最优解。刘佳宜等人提出一种改进的多目标离散蜜蜂算法来求解rdlb问题，结合非支配pareto排序的方法对种群进行排序，并利用空间干扰矩阵法生成可行解。fangyilin等人开发了一个特定问题的双准则进化算法来求解具有多机器人工作站的混合模型的拆卸线平衡问题。

但是，我们发现，拆卸线平衡问题(dlb)问题通常是一个多目标优化问题，其目标是使总利润最大化，周期时间最小化，工作站和资源数量最小化，并尽早去除危险部件。在dlb中，平衡率是用来评估工作站之间的工作量，它是评价dlb求解性能的主要指标之一。因此，我们提出一个以最大化拆卸利润和平衡率以及最小化拆卸能耗为目标的多目标多产品机器人拆卸线平衡问题(mmrdp)，在这些目标之间做出良好的权衡是至关重要的。

近年来，本领域的学者们认为头脑风暴优化算法(bso)是一种有前途的新群智能算法。它能平衡全局搜索和局部搜索，并已成功地应用于许多优化问题的求解。例如段海滨等人提出了一种捕食者-食饵bso来解决直流无刷电机的优化问题，通过结合捕食者-食饵机制在bso的基础上提高了种群的多样性。段海滨等人提出一种量子行为的bso来解决loney螺线管问题，通过结合量子行为来提高种群的多样性。但是，目前这些算法都不适用于求解rdlb问题，更不适用于解决mmrdp。因此，现在需要一种能够解决mmrdp的控制方法。

技术实现要素：

本发明旨在克服上述问题，提供一种随机多产品机器人拆卸线平衡控制方法，本发明根据多目标多产品机器人拆卸线平衡问题(mmrdp)的特征设计了一种基于pareto改进的多目标头脑风暴优化算法(pimbo)对其进行求解。本发明的控制方法可以解决多目标多产品机器人拆卸线平衡问题，并给出较优的解决方案。

本发明建立了一个基于and/or图并考虑拆卸失败风险的随机多产品机器人拆卸线平衡模型，目的是使拆卸利润和平衡率最大化，能耗最小化。本发明建立的mmrdp的数学模型，其中考虑了拆卸过程的不确定性。在实际的拆卸过程中，由于废旧产品质量的差异和拆卸的难度，拆卸时间的不确定性导致拆卸操作无法按计划进行。因此，应考虑拆卸失败风险。

为了实现上述目的，本发明采用如下技术方案：一种随机多产品机器人拆卸线平衡控制方法，包括以下步骤：

s1、建立以拆卸利润、平衡率最大化和能耗最小化为目标的数学模型；

s2、种群初始化：在搜索空间中随机生成多种解，并采用调整方法来初始化解中不符合优先和冲突约束条件的个体，形成初始种群；

s3、对初始种群进行解码：将需要执行的拆卸操作分配给每个工作站，确定每个工作站上拆卸操作的顺序，并将机器人分配到每个工作站执行拆卸操作；

s4、评估步骤s3得到的每个可行解的目标函数值；

s5、采用pimbo算法计算得到非支配解集并存放于外部档案中，并对外部档案进行演进；

s6、停止标准：利用适应度值来确定迭代次数，判断当前适应度值fv是否达到预设的总适应度值ftv，初始适应度值等于0，当评估新个体时，适应度值增加1，直到当前适应度值fv达到预设的总适应度值ftv，算法停止迭代，输出外部档案中的非支配解集；否则，重复s5。

进一步地，所述步骤s1具体包括以下步骤：

s101：基于多产品的and/or图对需要拆卸的多产品建立优先矩阵和拆卸关联矩阵如下：

(1)优先矩阵p＝[pjk]，表示两个拆卸操作的优先级和冲突关系，定义如下：

(2)拆卸关联矩阵d＝[dij]：用于表示组件和拆卸操作之间的关系，定义如下：

其中，i为组件索引，i∈{1,2,…,n^g}，其中n^g表示产品g中的组件数量；

j,k为操作索引，j,k∈{0,1,…,j^g}，其中j^g表示产品g中的操作数，0是虚拟操作；

pjk为优先矩阵p的第j行和第k列中的一个元素；

dij为拆卸关联矩阵d的第i行和第j列中的一个元素；

s102：建立以拆卸利润、平衡率最大化和能耗最小化为目标的数学模型：

其中，各符号代表的含义为：

g，表示eol产品索引，其中g代表废旧产品的数量；

i表示组件索引，i∈{1,2,…,n^g}，其中n^g表示产品g中的组件数量；

j,k表示操作索引，j,k∈{0,1,…,j^g}，其中j^g表示产品g中的操作数，0是虚拟操作；

l,m表示工作站索引，l,m∈{1,2,…,m}，其中m是工作站的数量；

r表示机器人类型索引，r∈{1,2,…,r}，其中r表示机器人类型的数量；

表示机器人r执行产品g中的操作j的拆卸时间；

表示机器人r执行产品g中的操作k的拆卸时间；

表示机器人r准备执行产品g中的操作k的设置时间，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置时间，操作k是在产品g中的操作j之后立刻执行；

表示机器人r执行产品g中的操作j的单位时间能耗；

表示机器人r准备执行产品g中的操作k的设置单位时间能耗，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置单位时间能耗，操作k是在产品g中的操作j之后立刻执行；

el表示第l个工作站运行期间的单位时间能耗；

表示机器人r执行产品g中的操作j的单位时间成本；

表示机器人r执行产品g中的操作k的单位时间成本；

表示机器人r准备执行产品g中的操作k的设置单位时间成本，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置单位时间成本，操作k是在产品g中的操作j之后立刻执行；

cl表示第l个工作站运行期间的单位时间成本；

表示产品g中组件i的回收利用价值；

tl表示第l个工作站的周期时间；

表示r机器人在产品g中操作j之后立即执行产品g中操作k的失败概率；

表示r机器人在产品g中操作j之后立即执行产品中操作k的失败概率；

θ表示拆卸过程失败成本的概率，它要满足要求的最小概率即拆卸过程的失败成本小于或等于最大失败成本；

表示拆卸过程的最大失败成本；

决策变量分别为：

表示如果产品g中的操作j由机器人r执行，则否则

表示如果产品g中的操作k由机器人r执行，则否则

表示如果产品g中的操作k由机器人r在产品g中的操作j之后立即执行的，则否则

表示如果产品中的操作k由机器人r在产品g中的操作j之后立即执行，则否则

表示如果产品g中的操作j由机器人r执行并分配给第l个工作站，则否则

表示如果产品g中的操作k由机器人r执行并分配给第m个工作站，则否则

ul表示如果使用第l个工作站，ul＝1；否则ul＝0；

其中目标函数(1)是最大化期望拆卸利润，它等于拆卸总收入减去拆卸总成本，所述拆卸总成本包括拆卸操作的总成本、同一产品中相邻操作的总设置成本、不同产品中相邻操作的总设置成本以及工作站的成本；目标函数(2)是最小化期望能耗，包括拆卸操作的总能耗、同一产品中相邻操作的总设置能耗、不同产品中相邻操作的总设置能耗以及工作站的能耗；目标函数(3)是最大化期望平衡率，即总工作时间除以开启的工作站总循环时间的平均值；

式(4)-(12)为约束条件：

约束条件(4)是保证每个产品中至少一个拆卸操作被执行，除了虚拟操作；约束条件(5)是表示多产品中的每个拆卸操作最多只能执行一次；约束条件(6)是表示多产品中的每个拆卸操作只能分配给一个工作站，由一个机器人执行；约束条件(7)是每个开启的工作站至少分配一个操作；约束条件(8)是确保多产品的可行拆卸序列满足优先约束；约束条件(9)是确保多产品的可行拆卸序列满足冲突约束；约束条件(10)表示所使用的所有工作站的操作时间不超过拆卸线总周期时间；(11)是要求拆卸过程的失败成本小于或等于其最大失败成本的概率大于预设值θ>0；约束条件(12)给出决策变量的范围。

进一步地，所述s2种群初始化的具体步骤如下：

s201、随机生成解π；

s202、根据优先矩阵p调整π¹中的拆卸操作顺序，以满足多产品拆卸操作之间的优先关系：遍历整个拆卸操作序列π¹，如果拆卸操作oj在ok之后并且优先矩阵p中pjk＝1，则oj和ok交换，同时rj和rk交换，否则继续遍历下一个操作；

s203、如果π¹中的拆卸操作oj是产品的开始拆卸操作，π²中的二进制值xj＝0，则xj重置为1，否则进入步骤s204；

s204、根据优先矩阵p调整π²中的二进制值，以满足同一产品拆卸操作之间的优先关系：遍历整个二进制向量π²，如果π²中的二进制值xk＝1，优先矩阵p中的pjk＝1，则π²中的二进制值xj设置为1；

s205、根据优先矩阵p调整π²中的二进制值，消除同一产品拆卸操作之间的冲突关系：遍历整个二进制向量π²，如果π²中的二进制值xk＝1，优先矩阵p中的pjk＝-1，则π²中的二进制值xk设置为0；

s206、将所有符合优先和冲突约束条件的个体作为初始种群；

进一步地，所述步骤s201中解π的编码具体过程是：采用一个三向量来表示一个解，即π＝(π¹，π²，π³)，其中π¹＝(o1，o2，…，oj)是一个整数串，表示拆卸操作序列，其中oj表示第j位的拆卸操作，j为多个待拆卸产品的总操作数；π²＝(x1，x2，…，xj)是一个整数串，表示包含j个二进制值的二进制向量；如果xj＝1，则执行π¹中的操作oj；否则，不执行oj；π³＝(r1，r2，…，rj)是一个整数串，表示机器人类型序列，其中每个元素都是介于1和r之间的整数，对应于执行一个拆卸操作的机器人类型。

进一步地，所述s3的具体步骤为：

s301、设置拆卸操作k＝1；

s302、确定是否执行拆卸操作k；如果π²中的二进制值xk＝1且k＝kf，其中kf是整个拆卸序列中的第一个拆卸操作，则打开一个新的工作站l＝1，第l个工作站的操作时间tld设置为0，第l个工作站中的拆卸操作集cl设置为ζ，其中ζ代表一个空集，并且第l个工作站的机器人集合rl设置为r，其中r表示执行拆卸操作k的机器人类型，如果xk＝1且k≠kf，则执行拆卸操作k并转至步骤302；否则，转至s308；

s303、设置第l个工作站当前操作时间tld＝tld+tjkr，其中tjkr表示拆卸操作k在拆卸操作j之后立即被机器人r执行的设置时间，并确定第l个工作站周期时间tl与第l个工作站当前操作时间tld的关系，如果tld>tl，则转至s304，否则转至s305；

s304、打开一个新的工作站，设置l＝l+1，cl＝ζ，rl＝r，tld＝tjkr；

s305、设置tld＝tld+tkr，其中tkr表示机器人r执行拆卸操作k的拆卸时间；并确定tl和tld之间的关系，如果tld>tl，则转至s306，否则转至s307；

s306、打开一个新的工作站，设置l＝l+1，cl＝ζ，rl＝r和tld＝tjkr+tkr；

s307、把拆卸操作k分配在第l个工作站上，更新cl和rl，cl＝cl∪k，rl＝rl∩r；

s308、如果k≠j，则设置k＝k+1，重复步骤s302-s307；否则，输出一个可行解；

进一步地，所述步骤s4采用蒙特卡罗模拟的随机模拟方法评估步骤s3得到的每个可行解的目标函数值。

进一步地，所述步骤s5具体包括以下步骤：

s501、对s4得到的每个可行解进行基于pareto的聚类操作，确保同一等级的个体分在同一类中；

s502、新个体生成操作：通过对当前种群的个体执行交叉算子和变异算子产生新个体，存放到新个体集合中；

s503、选择操作：通过pareto规则比较步骤s502计算得到的新个体集合和当前种群组成的混合种群的目标函数值，得到非支配解集并存放于外部档案中，并对外部档案进行演进。

与现有技术相比，本发明的有益效果在于：本发明根据多目标多产品机器人拆卸线平衡问题(mmrdp)的特征设计合理可行的编码方式以及解码方法。原始的头脑风暴优化算法(bso)采用k-means聚类算法实现对种群的聚类操作，该算法计算量大，效率低，需要人为确定聚类数，而且不适用本问题的求解方法。为此本发明选择非支配快速排序方法对种群进行聚类操作，该方法可以减少计算量，并且不需要人为确定聚类数。考虑到本发明研究的问题是离散组合问题，原始的bso通过加一个扰动随机数来产生新解，不适用该问题的求解，因此本发明结合遗传算法的思想，设计基于位置的交叉算子和基于位置的变异算子来产生新解。在原始的bso中，当从一个子类中选择个体生成新个体时，只有一个个体参与生成新个体。本发明将第一个子类的类中心与当前选择的个体相结合生成新的个体，这将提高算法的收敛速度。同时考虑到本发明研究的问题是多目标优化问题，本发明选择pareto规则来比较个体之间的好坏，并采用外部文档技术存储求解过程中获得的非支配解集。本发明采用基于pareto改进的多目标头脑风暴优化算法(pimbo)的控制方法可以解决多目标多产品机器人拆卸线平衡问题，并给出较优的解决方案。

附图说明

图1是本发明的流程示意图；

图2是圆珠笔的示意图；

图3是灯的示意图；

图4中(i)和(ii)分别是圆珠笔和台灯的and/or图形；

图5是圆珠笔和台灯的优先矩阵p的示意图；

图6是圆珠笔和台灯的关联矩阵d的示意图；

图7是圆珠笔和台灯的机器人拆卸线布局示例；

图8是对初始种群进行解码的流程图；

图9是交叉操作的示意图；

图10是变异操作的示意图；

图11是锤钻示意图；

图12是洗衣机示意图；

图13是无线电机示意图；

图14是锤钻的and/or图形；

图15是洗衣机的and/or图形；

图16是收音机的and/or图形；

图17是五种算法igd-度量实验结果的箱线图；

图18是五种算法超体积-度量实验结果的箱线图；

具体实施方式

以下结合具体的实例对本发明所述作进一步说明。有必要在此指出的是本实例仅用于例证的目的，本发明的保护范围不局限例证范畴。

实施例

本实施例将结合实际中的eol产品圆珠笔和台灯的机器人拆卸线平衡问题进行描述。

一、参照图1，本发明提供了一种随机多产品机器人拆卸线平衡控制方法，包括以下步骤：

s1、建立以拆卸利润、平衡率最大化和能耗最小化为目标的数学模型；

s101：基于多产品的and/or图对需要拆卸的多产品建立优先矩阵和拆卸关联矩阵如下：

(1)优先矩阵p＝[pjk]，表示两个拆卸操作的优先级和冲突关系，定义如下：

(2)拆卸关联矩阵d＝[dij]：用于表示组件和拆卸操作之间的关系，定义如下：

其中，i为组件索引，i∈{1,2,…,n^g}，其中n^g表示产品g中的组件数量；

j,k为操作索引，j,k∈{0,1,…,j^g}，其中j^g表示产品g中的操作数，0是虚拟操作；

pjk为优先矩阵p的第j行和第k列中的一个元素；

dij为拆卸关联矩阵d的第i行和第j列中的一个元素；

and/or图可以表示拆卸操作和拆卸操作之间的优先关系和冲突关系，也可以表示组件和拆卸操作之间的关系。以需要拆卸的圆珠笔和台灯为例，如图2和图3分别为圆珠笔和台灯的结构图，圆珠笔和台灯的and/or图分别为图4(i)和(ii)所示。在图4(i)和(ii)中，拆卸操作用整数1、2……和j表示，其中j是拆卸操作的总数。节点表示索引从<1>到<i>的组件，其中i是所有组件的数量。如果组件b是通过拆卸a获得的，则将组件a称为组件b的父组件，组件b称为组件a的子组件，其中a，b∈{1，2，……，i}。例如，组件<1>是<2>和<15>的父级，而<2>是<1>的子级。因此，基于圆珠笔和台灯的and/or图，根据所述优先矩阵p的定义，圆珠笔和台灯的优先矩阵p如图5所示，圆珠笔和台灯的拆卸关联矩阵d如图6所示。

机器人拆卸线的描述：如图5所示，一个圆珠笔包含15个组件和13个拆卸操作，一个台灯包括7个组件和6个拆卸操作。机器人拆卸线在图7中，三角形中的数字表示拆卸操作，圆圈中的数字表示执行拆卸操作的机器人类型，矩形中的数字表示不可分离的组件，传送带上的编号表示圆珠笔和台灯分解成不可分离组件的待拆卸组件。机器人拆卸线有几个工作站，每个机器人执行分配给每个工作站的拆卸操作。选择拆卸圆珠笔和台灯。因此，拆卸操作1-13和组件1-15与圆珠笔有关，拆卸操作14-19和组件16-22与台灯有关。i型机器人被称为机器人ri。一个可行的拆卸方案包括一个拆卸操作序列1、3、14、5、17、7、10和一个机器人类型序列1、3、2、3、2、2、1。拆卸操作序列中的每个拆卸操作根据周期时间的约束分配给每个工作站。每个工作站的机器人类型基于分配给每个工作站的拆卸操作的机器人类型。例如，拆卸操作1、3和14分配给第一个工作站，执行这些拆卸操作的机器人类型分别为1、3和2。因此，r1、r2和r3被分配到第一个工作站。另外，不可分离组件15、14、22、10、21、20、11、12、13由一个可行的拆卸方案得到，即将r1、r2和r3分配到第一个工作站，在第一个工作站执行拆卸操作1、3和14，将r2和r3分配到第二个工作站，在第二个工作站执行拆卸操作5和17，以及将r3分配到第三个工作站，在第三个工作站执行拆卸操作7和10。

s102：建立以拆卸利润、平衡率最大化和能耗最小化为目标的数学模型：

为了建立数学模型，首先我们需要以下假设：

1.对于每种类型的机器人，每个拆卸操作的平均拆卸时间和设置时间已知。

2.对于每种类型的机器人，每个拆卸操作的单位时间拆卸成本和单位时间拆卸能耗已知。

3.单位时间的转换成本和单位时间的转换能耗已知。

4.eol产品的供应是无限的。

5.只有机器人负责执行拆卸操作。任何一个机器人一次只能执行一个拆卸操作，每一个操作都必须由一个机器人执行。

6.需要拆卸的多个废旧产品的and/or图都已知。

7.如果拆卸任务分配给工作站并由某种类型的机器人执行，则该机器人类型被分配给该工作站。

下面列出了本发明涉及的符号和决策变量。

1)符号

g，表示eol产品索引，其中g代表废旧产品的数量；

i表示组件索引，i∈{1,2,…,n^g}，其中n^g表示产品g中的组件数量；

j,k表示操作索引，j,k∈{0,1,…,j^g}，其中j^g表示产品g中的操作数，0是虚拟操作；

l,m表示工作站索引，l,m∈{1,2,…,m}，其中m是工作站的数量；

r表示机器人类型索引，r∈{1,2,…,r}，其中r表示机器人类型的数量；

表示机器人r执行产品g中的操作j的拆卸时间；

表示机器人r执行产品g中的操作k的拆卸时间；

表示机器人r准备执行产品g中的操作k的设置时间，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置时间，操作k是在产品g中的操作j之后立刻执行；

表示机器人r执行产品g中的操作j的单位时间能耗；

表示机器人r准备执行产品g中的操作k的设置单位时间能耗，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置单位时间能耗，操作k是在产品g中的操作j之后立刻执行；

el表示第l个工作站运行期间的单位时间能耗；

表示机器人r执行产品g中的操作j的单位时间成本；

表示机器人r执行产品g中的操作k的单位时间成本；

表示机器人r准备执行产品g中的操作k的设置单位时间成本，操作k是在产品g中的操作j之后立刻执行；

表示机器人r准备执行产品中的操作k的设置单位时间成本，操作k是在产品g中的操作j之后立刻执行；

cl表示第l个工作站运行期间的单位时间成本；

表示产品g中组件i的回收利用价值；

tl表示第l个工作站的周期时间；

表示r机器人在产品g中操作j之后立即执行产品g中操作k的失败概率；

表示r机器人在产品g中操作j之后立即执行产品中操作k的失败概率；

θ表示拆卸过程失败成本的概率，它要满足要求的最小概率即拆卸过程的失败成本小于或等于最大失败成本；

表示拆卸过程的最大失败成本；

决策变量分别为：

表示如果产品g中的操作j由机器人r执行，则否则

表示如果产品g中的操作k由机器人r执行，则否则

表示如果产品g中的操作k由机器人r在产品g中的操作j之后立即执行的，则否则

表示如果产品中的操作k由机器人r在产品g中的操作j之后立即执行，则否则

表示如果产品g中的操作j由机器人r执行并分配给第l个工作站，则否则

表示如果产品g中的操作k由机器人r执行并分配给第m个工作站，则否则

ul表示如果使用第l个工作站，ul＝1；否则ul＝0；

建立的以拆卸利润、平衡率最大化和能耗最小化为目标的数学模型如下：

其中目标函数(1)是最大化期望拆卸利润，它等于拆卸总收入减去拆卸总成本，所述拆卸总成本包括拆卸操作的总成本、同一产品中相邻操作的总设置成本、不同产品中相邻操作的总设置成本以及工作站的成本；目标函数(2)是最小化期望能耗，包括拆卸操作的总能耗、同一产品中相邻操作的总设置能耗、不同产品中相邻操作的总设置能耗以及工作站的能耗；目标函数(3)是最大化期望平衡率，即总工作时间除以开启的工作站总循环时间的平均值；

式(4)-(12)为约束条件：

约束条件(4)是保证每个产品中至少一个拆卸操作被执行，除了虚拟操作；约束条件(5)是表示多产品中的每个拆卸操作最多只能执行一次；约束条件(6)是表示多产品中的每个拆卸操作只能分配给一个工作站，由一个机器人执行；约束条件(7)是每个开启的工作站至少分配一个操作；约束条件(8)是确保多产品的可行拆卸序列满足优先约束；约束条件(9)是确保多产品的可行拆卸序列满足冲突约束；约束条件(10)表示所使用的所有工作站的操作时间不超过拆卸线总周期时间；约束条件(11)是要求拆卸过程的失败成本小于或等于其最大失败成本的概率大于预设值θ>0；约束条件(12)给出决策变量的范围。

s2、种群初始化：在搜索空间中随机生成多种解，并采用调整方法来初始化解中不符合优先和冲突约束条件的个体，形成初始种群，具体步骤为：

s201、随机生成解π；

解π的编码具体过程是：采用一个三向量来表示一个解，即π＝(π¹，π²，π³)，其中π¹＝(o1，o2，…，oj)是一个整数串，表示拆卸操作序列，oj表示第j位的拆卸操作，j为多个待拆卸产品的总操作数；π²＝(x1，x2，…，xj)是一个整数串，表示包含j个二进制值的二进制向量；如果xj＝1，则执行π¹中的操作oj；否则，不执行oj；π³＝(r1，r2，…，rj)是一个整数串，表示机器人类型序列，其中每个元素都是介于1和r之间的整数，对应于执行一个拆卸操作的机器人类型；

s202、根据优先矩阵p调整π¹中的拆卸操作顺序，以满足多产品拆卸操作之间的优先关系：遍历整个拆卸操作序列π¹，如果拆卸操作oj在ok之后并且优先矩阵p中pjk＝1，则oj和ok交换，同时rj和rk交换，否则继续遍历下一个操作；

s203：如果π¹中的拆卸操作oj是产品的开始拆卸操作，π²中的二进制值xj＝0，则xj重置为1，否则进入步骤s204；

s204、根据优先矩阵p调整π²中的二进制值，以满足同一产品拆卸操作之间的优先关系：遍历整个二进制向量π²，如果π²中的二进制值xk＝1，优先矩阵p中的pjk＝1，则π²中的二进制值xj设置为1；

s205、根据优先矩阵p调整π²中的二进制值，消除同一产品拆卸操作之间的冲突关系：遍历整个二进制向量π²，如果π²中的二进制值xk＝1，优先矩阵p中的pjk＝-1，则π²中的二进制值xk设置为0；

s206、将所有符合优先和冲突约束条件的个体作为初始种群；

s3、进行解码：将需要执行的拆卸操作分配给每个工作站，确定每个工作站上拆卸操作的顺序，并将机器人分配到每个工作站执行拆卸操作；

解码过程是将需要执行的拆卸操作分配给每个工作站，确定每个工作站上拆卸操作的顺序，并将机器人分配到每个工作站执行拆卸操作。如果分配给第l个工作站的拆卸操作j由机器人r执行，则将机器人r分配给第l个工作站。此外，在分配拆卸操作时，需要满足工作站的周期时间约束。解码过程如图8所示，具体解码步骤如下：

s301、设置k＝1；

s302、确定是否执行k；如果π²中的二进制值xk＝1且k＝kf，其中kf是整个拆卸序列中的第一个拆卸操作，则打开一个新的工作站，设置l＝1，第l个工作站的操作时间tld设置为0，第l个工作站中的拆卸操作集cl设置为ζ，其中ζ代表一个空集，并且第l个工作站的机器人集合rl设置为r，其中r表示执行拆卸操作k的机器人类型，如果xk＝1且k≠kf，则执行拆卸操作k并转至步骤302；否则，转至s308；

s303、设置第l个工作站当前操作时间tld＝tld+tjkr，其中tjkr表示拆卸操作k在拆卸操作j之后立即被机器人r执行的设置时间，并确定第l个工作站周期时间tl与第l个工作站当前操作时间tld的关系，如果tld>tl，则转至s304，否则转至s305；

s304、打开一个新的工作站，设置l＝l+1，cl＝ζ，rl＝r，tld＝tjkr；

s305、设置tld＝tld+tkr，其中tkr表示机器人r执行拆卸操作k的拆卸时间，并确定tl和tld之间的关系，如果tld>tl，则转至s306，否则转至s307；

s306、打开一个新的工作站，设置l＝l+1，cl＝ζ，rl＝r和tld＝tjkr+tkr；

s307、把拆卸操作k分配在第l个工作站上，更新cl和rl，cl＝cl∪k，rl＝rl∩r；

s308、如果k≠j，则设置k＝k+1，重复步骤s302-s307；否则，输出一个可行解；

s4、采用蒙特卡罗模拟的随机模拟方法评估步骤s3得到的每个可行解的目标函数值；

本发明设计了一种基于蒙特卡罗模拟的随机模拟方法。设π为一个解，γ表示样本数等于10。采用以下步骤来评估π的目标函数值和可行性。需要注意的是，由于目标函数f1和f3是求它们的最大值，因此将得到的目标函数f1和f3的值分别乘以-1，以便于比较解的质量。

步骤s401：设置f1(π)：＝0,f2(π)：＝0,f3(π)：＝0,fc(π)：＝0，γ：＝10；

步骤s402：生成δ个样本，每个样本ρδ是根据拆卸和设置时间的正态分布产生的，δ＝1，2，…，γ；

步骤s403：计算拆卸利润f1δ(π)、能耗f2δ(π)、平衡率f3δ(π)和拆卸失败成本fcδ(π)，其中δ＝1，2，…，γ；

步骤s404：对所有样本的拆卸失败成本fcδ(π)按升序排序，δ′等于如果π是可行的，fc(π)：＝fcδ′(π)，然后转到步骤s405。否则，π是不可行的，f1(π)、f2(π)和f3(π)分别等于一个较大的整数，并进入步骤s406。

步骤s405：将预期拆卸利润f1(π)、能耗f2(π)和平衡率f3(π)近似为：

步骤s406：返回f1(π)、f2(π)和f3(π)作为所述步骤s1中目标函数(1)、目标函数(2)和目标函数(3)的近似值。

s5、采用pimbo算法(基于pareto改进的多目标头脑风暴优化算法)计算得到非支配解集并存放于外部档案中，并对外部档案进行演进，具体包括以下步骤：

s501、基于pareto的聚类操作：根据pareto规则对s4得到的每个可行解(个体)进行排序并进行聚类，确保同一等级的个体分在同一类中；具体的聚类操作如下：

步骤s5011：初始化两个集合δ和子类cp。根据pareto规则对步骤s4得到的每个个体进行划分等级；δ设置为步骤s4得到的种群，cp设置为空集合，设置一个整数γ表示一个个体的等级，γ＝1；

步骤s5012：从δ中选择所有等级为γ的个体添加到cp中，并把这些个体从δ中移除；

步骤s5013：把cp添加到一个类簇c，并且选择cp中的第一个个体作为该子类cp的类中心；然后将cp设置为空集合，γ加1；

步骤s5014：如果δ不为空集合，重复步骤s5012-s5013，否则，结束流程；

s502、新个体生成操作：通过对当前种群的个体执行交叉操作和变异操作产生新个体，存放到新个体集合中；

步骤s5021：随机生成一个[0,1)的值如果小于pc，其中pc是替换类中心的概率，随机选择一个类中心并生成一个新个体去替换它；否则转至步骤s5022；

步骤s5022：初始化一个新个体集合x，随机生成一个[0,1)的值如果小于pg，其中pg是从一个子类中生成新个体的概率，从类簇c中随机选择一个类，随机生成一个[0,1)的值如果小于po，其中po是根据子类的类中心生成新个体的概率，选择该类的类中心和第一个类的类中心执行基于位置的交叉算子(pbx)产生一个新个体，并对该新个体执行基于位置的变异算子(pbm)，然后对这个新个体通过s3进行调整至可行解，然后通过s4进行评估，并保存到x；

否则随机选择该类的其他个体和第一个类的类中心执行基于位置的交叉算子(pbx)产生一个新个体，并对该新个体执行基于位置的变异算子(pbm)，然后对这个新个体通过s3进行调整至可行解，然后通过s4进行评估，并保存到x，否则转至步骤3；

步骤s5023：从类簇c中随机选择两个类，随机生成一个[0,1)的值如果小于pt，其中pt是结合两个子类的类中心生成新个体的概率。分别选择这两个类的类中心执行基于位置的交叉算子(pbx)产生一个新个体，并对该新个体执行基于位置的变异算子(pbm)，然后对这个新个体通过s3进行调整至可行解，然后通过s4进行评估，并保存到x，否则分别随机选择这两个类的其他个体执行基于位置的交叉算子(pbx)产生一个新个体，并对该新个体执行基于位置的变异算子(pbm)，然后对这个新个体通过s3进行调整至可行解，然后通过s4进行评估，并保存到x；

步骤s5024：直到x中的个体数达到种群数，结束流程。

为了增加解的多样性，步骤s5023采用的交叉算子(pbx)采用以下方法实现：

步骤(1)：生成一个[1,j]的随机整数并初始化一个新个体πnew＝(πnew¹,πnew²,πnew³)；

步骤(2)：随机生成个[1,j]的不同整数，然后把一个父代π1中的π1¹中这几个整数所在位置的元素存到πnew¹中相同位置的元素，同样πnew²和πnew³执行相同的操作；

步骤(3)：πnew¹中其他位置的元素由另一个父代π2中的π2¹中和πnew¹中已有元素全都不同的元素有序替换，同样πnew²和πnew³执行相同的操作。

例如，如图9所示，π1和π2被认为是基于位置的交叉操作的两个父代个体。从π1中随机选择几个位置。为了生成π1的新个体，π1被用作父代和参考个体。π1序列中选择的几个位置保持不变，其他位置依次通过π2的相同序列映射。在图9中，π1选择的几个位置是操作{2，16，3，12，13}。保持这些位置不变，并将{1}、{14、15}、{4、11}、{17、18、5}和{19、6、7、9、10、8}由π2中的相同序列依次映射。然后根据π1对π2和π3的相同位置进行相同的运算，可以得到一个新的个体πnew。

同时，步骤s5023采用的基于位置的变异算子(pbm)采用以下方法实现：

设π＝(π¹,π²,π³)是需要突变的个体，和p分别是区间[0,1]中的随机数和区间[0,3]中的随机整数。η是个体执行变异的概率。

步骤1：产生随机数如果结束流程；否则转至步骤2；

步骤2：产生随机数p，如果p＝0，随机产生2个[1,j]的不同整数μ和γ，然后分别交换π¹，π²，π³中μ位置和γ位置的元素；如果p＝1，随机产生1个[1,j]的整数，然后把π²中这个整数位置对应的元素变为另一个二进制数；如果p＝2，随机产生1个[1,j]的整数，然后把π³中这个整数位置对应的元素变为另一个机器人类型；即q是一个随机二进制数，如果元素值为1，当q＝0时，1变为2，q＝1时，1变为3；如果元素值为2，当q＝0时，2变为1，q＝1时，2变为3；如果元素值为3，当q＝0时，3变为1，q＝1时，3变为2；如果p＝3，依次对π执行当p＝0，p＝1以及p＝2的操作。

例如，在图10中，如果p＝0，π¹中的两个位置是操作2和13。保持{1}、{14、15、16、3、4、11、12、17、18、5}和{19、6、7、9、10、8}不变，以及交换{2}和{13}。然后，根据π¹对π²和π³中相同的位置执行相同的操作。如果p＝1，则只调整π里的π²，π²中的一个位置是二进制值0，保持π2中的其他位置不变，并将0变为二进制数1。如果p＝2，只修改π里的π³，π³中的一个位置是r1机器人，保持π³中的其他位置不变，将r1机器人改为其他机器人类型。q是一个二进制数，当q＝0时，1变为2，q＝1时，1变为3。如果p＝3，则需要执行上述所有操作。

s503、选择操作：通过pareto规则比较步骤s502计算得到的新个体集合和当前种群组成的混合种群，得到非支配解集并存放于外部档案中，并对外部档案进行演进，具体步骤如下：

目标函数表示为maxf1(π)、minf2(π)和maxf3(π)，其中π∈ω，π是可行解，ω是所有可行解的集合。对于如果满足fi(π¹)≤fi(π²)，并且满足fi(π¹)<fi(π²)，则π¹支配π²，表示为如果和则π*是pareto最优解，s＝{π*}称为pareto解集，即非支配解集。

然后构建一个外部档案来存储获得的非支配解。首先，利用所有新生成的解基于pareto规则对外部档案进行更新。在每次迭代完成后，利用当前种群基于pareto规则对外部档案进行更新，以确保外部档案中所有的非支配解始终是当前的全局非支配解。

s6:停止标准：利用适应度值来确定迭代次数，初始适应度值等于0，当评估新个体时，适应度值增加1，直到当前适应度值fv达到预设的总适应度值ftv，算法停止迭代，输出外部档案中的非支配解集；否则，重复s5。

具体实验实例

为了测试本发明的基于pareto改进的多目标头脑风暴优化算法(pimbo)的随机多产品机器人拆卸线平衡控制方法在处理多目标多产品机器人拆卸线平衡问题(mmrdp)的性能，采用一组废旧(eol)产品进行了实验。此外，选择mdgwo和nsgaii(即多目标离散灰狼优化算法和非支配排序遗传算法ii)作为比较算法。

1、测试实例生成

我们采用的实际需要拆卸的废旧(eol)产品案例。它们如下：圆珠笔(bp)，锤钻(hd)，洗衣机(wm)，无线电机(rs)。一个bp包含15个组件，拆卸它需要13个拆卸操作。一个hd包括62个组件，含有46个拆卸操作。一个wm包含15个组件，13个拆卸操作用于拆卸它。一个rs中有29个组件，含有30个拆卸操作。hd、wm和rs的结构示意图分别如图11-图13所示。其and/or图如图14-图16所示。我们采用bp、hd、wm和rs生成10个测试实例，每个实例包含多个拆卸产品。它们的对应关系如表1所示，其中表示每个实例对应的拆卸过程的最大失败成本。

表1bp、hd、wm和rs的测试实例

2、性能指标

为了分析本发明的基于pareto改进的多目标头脑风暴优化算法(pimbo)及其对等算法(mdgwo和nsgaii)的性能，我们使用三个指标来比较它们的结果，即c-度量(解集覆盖率)、igd-度量(反转世代距离(igd，invertedgenerationaldistance))和超体积-度量。

c-度量可以度量两个解集的支配百分比。

igd-度量的目的是计算得到的解到最优解集中最近解的平均距离。解集的igd越小，逼近性和分布性越好。请注意，我们不能预先获得我们所考虑问题的最优解集。因此，将本发明基于pareto改进的多目标头脑风暴优化算法(pimbo)及其对等算法(mdgwo和nsgaii)得到的所有解集合并，然后将得到的所有非支配解作为最优解集。另外，将所有解的目标函数值归一化为[0,1]，然后计算igd值。

超体积-度量主要是计算由所得到的解和所有解集中包含三个最大目标值的参考点所构成的体积。解集的超体积越大，逼近性和分布性越好。计算超体积时，将所有解的目标函数值归一化为[0,1]，参考点设为(1,1,1)。

对于每个实例，pimbo及其对等算法执行20次，并计算三个度量的平均值来评估它们的性能。此外，我们采用单尾t检验分析实验结果，其自由度为38，显著性水平为0.05。采用符号“+”、“-”或“～”表示结果。符号“+”表示pimbo明显优于其对等算法，“-”表示pimbo明显差于其对等算法，而“～”表示pimbo在统计上与其对等算法相当。

3、参数设置

本发明的基于pareto改进的多目标头脑风暴优化算法(pimbo)的参数设置为：pc＝0.0，pg＝0.8，po＝0.5，pt＝0.4，其中pc是替换类中心的概率，pg是从一个子类中生成新个体的概率，po是根据子类的类中心生成新个体的概率，pt是结合两个子类的类中心生成新个体的概率。所有算法的参数设置如下：种群大小q为120，变异概率η为0.3，总适应值ftv为2000。目标函数涉及的参数设为：cl＝0.03，el＝0.05，tl＝50，θ＝0.90，r＝3。所有算法独立运行20次，在设计的10个测试实例中的平均运行时间对比情况如表2所示。从中我们可以得出结论，pimbo可以用较少的时间找到一个有效的解集。

表2三种算法在10个测试实例中的平均运行时间对比情况

4、案例分析

从表3中，pimbo为解决测试实例2(bp,hd)提供了10个非支配解，其中f1、f2、f3、fc和rs分别表示总利润、能耗、平衡率、失败成本和分配给开启的工作站的机器人类型集。需要注意的是，bp有13个拆卸操作，而hd有46个拆卸操作。因此，这些解中的拆卸操作1-13代表bp的拆卸操作，14-59代表hd的拆卸操作。例如，在第一个解中，拆卸操作1来自bp，拆卸操作14、15、16、18、22、29、24、31和33来自hd。拆卸操作14、15、1和16由第一个工作站的r2、r1和r3机器人执行，拆卸操作18、22和29由第二个工作站的r3机器人执行，拆卸操作24、31和33由第三个工作站的r1和r3机器人执行。总利润为827.66，能耗为976.31，平衡率为0.8395，失败成本为48.95。因此，这种针对所研究问题的解决方案是可行的。

表3通过pimbo算法得到了实例2拆卸序列的结果

5、实验结果

使用三种算法来处理10个测试实例，最终使用c-度量、igd-度量和超体积-度量来分析它们的实验结果。

五种算法的c-度量参照表4，表4中采用a、b和e分别表示pimbo、mdgwo和nsgaii。从这些实验结果中可以看出，由于c(a，b)和c(a，e)分别大于c(b，a)和c(e，a)，pimbo得到的大多数解都可以优于mdgwo和nsgaii的解。c(a，b)，c(b，a)，c(a，e)和c(e，a)的平均值分别为0.8150、0.0604、0.9136和0.0202。这些结果证实了pimbo在处理相关问题上的表现优于mdgwo和nsgaii。

表4三种算法的c-度量实验结果比较

三种算法的igd指标结果如表5所示。从这些结果，我们可以得出结论，pimbo在解决相关问题方面比mdgwo和nsgaii更好，因为pimbo的igd值在所有情况下都小于mdgwo和nsgaii。三种算法在10个实例上的igd值的平均值分别为0.1077、0.2084、和0.2457。结果表明，pimbo能找到一个收敛性好、多样性好的非支配解集。

表5三种算法的igd-度量实验结果比较

三种算法的超体积值如表6所示。结果表明，在所有情况下，pimbo得到的超体积值都大于nsgaii。而且pimbo在9个实例上得到的超体积值都大于mdgwo。三种算法在10个实例上的超体积值的平均值分别为0.5349、0.4568和0.4216。结果表明，与mdgwo和nsgaii相比，pimbo方法得到的解集具有更好的逼近性和更均匀的分布。

表6三种算法的超体积-度量实验结果比较

为了更清楚、更直观地展示实验结果，在10个实例上用igd-度量和超体积-度量绘制了三种算法的实验结果的箱线图，如图17和图18所示。在图18中，“hv”表示超体积。根据图17和图18的实验结果，可以看到pimbo的igd值和hv值在所有实例中都分别小于和大于mdgwo和nsgaii，因此可以得出结论，pimbo在处理所提出的问题方面比mdgwo和nsgaii有更好的性能。

因此，可以看出，本发明的pimbo在不确定环境下，以利润最大化、平衡率最大化和能耗最小为目标的多目标多产品机器人拆卸线平衡问题(mmrdp)，建立了相应的数学模型。然后，设计了一种结合随机模拟的pimbo算法进行求解。通过比较实验，可以看出本发明所提出的方法在解决该问题上表现出了良好的性能。为考虑不确定性和多目标优化的mmpdp提供了一种有效的方法。