3D格子气自动机算法的空间并联机构工作空间搜索算法的研究

2022-11-21陆彩满刘艳梨吴洪涛

机械设计与制造 2022年11期

陆彩满，刘艳梨，古亮，吴洪涛

（1.江苏安全技术职业学院机械工程系，江苏徐州 221011；2.南京航空航天大学机电学院，江苏南京 210016；3.重庆理工大学电气工程学院，重庆 400054）

1 引言

并联机构因其优点多而受到越来越多的国内外学者深入研究［1］，如控制精度高、本体刚度大、结构承载能力强等优点。并联机构工作空间是指并联机构操作平台执行末端参考点所能运动的区域［2］。并联机构路径规划、避障运动控制中都需要考虑工作空间，工作空间是一种评价机构灵活性的重要性能指标，也是机构设计与优化重要参考依据［3］。

目前并联机构工作空间分析方法主要有：几何绘图法、解析法、数值法［4］。其中，具有较好的通用性搜索方法是数值法。每种方法都有其适用范围且方法分析特点不一，文献［5］中是利用工作空间边界弦无限逼近边界形式，把将工作空间分析转化为最优化计算问题，巧妙的将工程问题转换为简单的优化计算问题。文献［6］提出一种离散式蒙特卡洛数值搜索法，通过将一个转动自由度进行离散化处理，同时结合子空间的并集和交集的思想，使得可达工作空间与灵活工作空间分析简单化。最后以一种三平移一转动空间并联机构为例分析以及验证其效果。文献［7］基于粒子群算法原理基础上，提出一种改进的粒子群搜索算法，将工作工作空间搜索转化为目标函数优化问题，通过得到的一系列优化结果，将这些点的集合映射到工作空间边界上得到工作空间边界。文献［8］根据混沌映射原理实现工作空间快速搜索，同时以空间3T1R机构为例，验证其有效性、快速性。文献［9］针对蒙特卡洛法求解工作空间效率不高问题，将蒙特卡洛法进行改进，利用传统方法得到的工作空间，在此基础上细分成若干个子工作空间，根据标准动态可调的正态分布实现任意精度的阀值计算。文献［10］考虑到蒙特卡洛法中随机点分布的不均匀性所产生误差导致精度不够的问题。将前面的到边界点的小邻域内重新生成随机点，有效的解决蒙特拉洛法的精度不足、分布不均匀问题。

受文献［11］启发，借鉴格子气自动机算法模型，提出一种基于3D格子气自动机并行动态粒子搜索算法，同时重点介绍此方法的搜索原理以及分析过程，并以两个空间并联机构为例进行验证，根据实例仿真分析机构可达工作空间与定姿态工作空间，验证该算法的有效性，并分析了算法参数对算法搜索效率的影响。

2 理论基础介绍

2.1 元胞自动机

元胞自动机（Cellular Automata）是一种通用、动态的计算机仿真数学模型，其应用范围十分广泛，其应用涉及社会与自然科学的各个领域［12］。物理定义为有限的空间和时间并且元胞状态都是有限的离散的动力系统，所有元胞状态每一时刻中根据设定的元胞规则不断更新迭代。元胞自动机组成主要包括元胞、元胞空间、邻居及规则等［13］。它们之间的关系可理解为是由一个元胞空间和定义于该空间的变换函数所组成。与一般动力学模型不同的是元胞自动机不需要严格的物理方程或函数构成，只需要定义若干个转换规则就能实现动力系统的演化，即使是简单的规则也可以演化出复杂的结果是这个算法的一大亮点，而且还具有并行计算的特点。

2.2 3D格子气自动机模型

格子气动机是一种基于离散粒子运动论微观模型的数值计算方法，属于元胞自动机的一种，是元胞自动机在流体力学上的应用，三维格子气自动机模型常用于金属材料腐蚀损伤演化、林火蔓延仿真、晶体生长仿真、微观组织模拟等方面［14］。三维格子气自动机模型邻域类型根据由其邻域半径r分类，（1）r=1时，有6（上下左右前后）或者18、26 个运动方向（上下左右前后及其组合），（2）r=2时，有124个运动方向，以此类推，更大的邻居类型将趋向于随机矢量方向运动。以邻域个数26为例，其运动方向有26个，沿着四周所有可能邻域进行运动。3D格子气自动机离散速度模型，如图1所示。箭头所指方向为动态粒子运动方向，箭头位置为动态粒子的运动位置。

图1 3D格子气自动机速度模型Fig.1 Velocity Model of 3D Lattice Gas Automata

3 算法提出

3.1 3D格子气偏随机模型动态演化搜索算法原理

受到格子气动机模型理论启发，假设在n×n×n的三维空间格子中，有m个粒子随机分布在空间格子中，每个粒子有两个基本属性：运动方向和运动速度。（1）不同的随机运动方向，运动方向数量按照邻居类型不同进行划分。（2）运动速度，在每个更新步长内的粒子在某一运动方向上从原格子位置移动到另外的格子空间的格子数目，格子数目可以是1，2，3 等，运动速度即是1，2，3等。格子气动机具有不相容规则，即每个粒子同时只能占据一个格子空间。将3D元胞空间粒子点坐标作为输入，机构工作空间约束函数作为映射函数，输出符合机构工作空间约束的点，即可得到机构工作空间的数值模型。输入可以在经过压缩，平移等操作函数使得更符合实际的机构工作空间约束函数坐标矩阵。可在运动前后判断满足约束值的变化情况，若有变化，则判断为边界点，多次运动后，可得到边界点的集合。若不需要边界点，则直接判断点满不满足约束即可。

3.2 3D格子气偏随机模型动态演化搜索法流程

根据以上算法原理，具体的实现步骤过程如下：

（1）设定空间大小nx×ny×nz，nx，ny和nz为大于零的整数，设定粒子数目m和随机位置Pm（x，y，z），其中，x，y，z为整数位置的集合；设定映射到工作空间的映射规则，Pm（x，y，z）—f（x）→PPm（x，y，z），映射规则为平移和缩放等操作，使得PPm（x，y，z）适应待搜索机构的空间；设定粒子运动的速度；

（2）记录粒子的位置和工作空间位置，进入机构的约束条件，判断其是否满足约束，满足记为1，不满足记为0，记为ys1（PPm，｛0，1｝）矩阵中；产生随机步长数（1-nu）和方向数（1～26），分别赋予每个粒子的运动速度和方向。产生新的粒子位置Pm’（x，y，z），和工作空间的位置坐标PPm’（x，y，z）；

（3）进入机构的约束条件，判断其是否满足约束，满足记为1，不满足记为0，记为ys2（PPm’，｛0，1｝）矩阵中；计算ys2-ys1的结果，记为ysj（PPm’，｛0，1，-1｝），若集合中有不为零的元素，则判断为边界，可视化边界点和粒子点。

（4）重复（2）与（3），直至达到需要的运动时间和搜索效果时结束。其搜索流程图，如图2所示。

图2 算法搜索流程图Fig.2 Algorithm Search Flow Chart

4 算法仿真案例

4.1 算例1仿真分析

4.1.1 4PPaR机构运动学逆解分析

为了验证所提出的工作空间搜索方法的有效性，这里以Pi-Ri1(-4R-)//Ri2⊥Ri3，(i=1，2，3，4)并联机构为例，此结构简图，如图3所示。此机构支链结构多，此空间机构具有三平移一转动（3T1R）的运动特性。

图3 机构结构简图Fig.3 Structure of Mechanism

在建立运动学方程前，分别建立动静坐标系{P-uvw}、{O-xyz}，以参数(d1，d2，d3，d4)作为独立运动输入参数，以参数(x，y，z，θ)作为运动输出，运动学逆解是根据输出(x，y，z，θ)，得到输入(d1，d2，d3，d4)。运动学逆解方程建立过程如下：

（1）根据主动移动副的位置可用参数(d1，d2，d3，d4)表示，可得到移动副的坐标P1=(0，d1，0)，P2=(d2，0，0)，P3=()，P4=()，PiRi1与静平台夹角为αi，因此，Ri1，(i=1，2，3，4)表示为：

（2）动平台的位姿为(x，y，z，θ)，已知动平台Ri3，(i=1，2，3，4)，在{P-uvw} 下表示为R13(0，0)、R23(r，-r)、R33(2r，0)、R43(r，r)，Ri3，(i=1，2，3，4)需要通过转换矩阵转换到静坐标系中，并联机构从动坐标系｛O｝到静坐标系｛P｝的转换矩阵为TOP

（3）根据上述转换矩阵，可得点Ri3，(i=1，2，3，4)在静坐标系中的坐标。

（4）已知动平台上的转动副布置特点，均垂直动平台，即得到点Ri2，(i=1，2，3，4)在{O-xyz}中的坐标。

（5）由于平行四边形机构Pa结构不变，杆长为l2，建立如下所示方程等式：

将变量di分离，处理得到di的表达式：

其中，

工作空间分析是根据建立的逆解方程完成，主要根据逆解方程有解来判定是否为运动范围内的点，利用matlab编程实现工作空间分析以及搜索过程。

4.1.2 4PPaR机构的可达工作空间搜索

根据上述的上一节中算法搜索步骤与流程，选择具体的参数值进行，如表1所示。用红色颜色标识为随机移动的粒子，蓝色标识为工作空间边界的粒子；搜索算法参数设定如下：元胞空间大小为（60×60×60），粒子数量为5000个，移动速度为1，运动方向为26个方向。并得到一系列满足工作空间约束条件的粒子，动态搜索随时间演化的过程，如图4所示。最后得到的可达工作空间，如图5所示。搜索过程边界点变化图，如图6所示。

表1 4PPaR机构的尺寸参数Tab.1 Size Parameters of Mechanism

图4 机构的可达工作空间动态搜索过程图Fig.4 Dynamic Search Process of Reachable Workspace of Mechanism

图5 机构的可达工作空间搜索结果Fig.5 Reachable Workspace Search Results of Mechanism

图6 机构的搜索过程边界点的变化过程Fig.6 Search Process of Mechanism Change Process of Boundary Point

根据图4动态仿真图分析可得：时间t=1至t=102为全部搜索演化过程，刚开始并联机构可达工作空间在搜索的过程中在满足要求的蓝色粒子在不断地增加。最后得到满足约束条件的可达工作空间，如图5所示。所得机构工作空间边界光滑连续、无空洞情况，且边界点清晰，表明了此算法的动态可视化搜索能力。

边界点搜索过程的变化情况，可以发现每步搜索到点的个数呈锯齿状增加，每步得到的边界点在400个以上，表明了算法的并行计算能力，同时累计的增加呈直线趋势，在140多步时，边界点的个数已超过7万个，如图6所示。

4.2 算例2仿真分析

4.2.1 4PRPaRR机构运动学逆解分析

4PRPaRR机构结构简图，如图7所示。机构的逆解具体分析过程可参考文献［6］，就不再此详述，机构主动副参数是以移动副的移动距离(d1，d2，d3，d4)表示，而动平台的位置姿态输出参数(x，y，z，θ)。分析得到机构逆解表达式：

图7 并联机构结构简图Fig.7 Schematic Diagram of Parallel Mechanism

4.2.2 4PRPaRR并联机构定姿态工作空间搜索

4PRPaRR 并联机构定姿态工作空间搜索根据表2 中的参数值赋值，求解出姿态角为0°的定姿态工作空间，动态搜索随时间演化的过程，如图8 所示。最后得到定姿态工作空间，如图9 所示。

图8 机构的定姿态工作空间动态搜索过程图Fig.8 Dynamic Search Process of Fixed Pose Workspace for Mechanism

图9 机构定姿态工作空间搜索结果（0°）Fig.9 Search Results of Fixed Position Workspace（0°）

表2 机构的尺寸参数Tab.2 Size Parameters of Mechanism

由上图8动态仿真图分析：t=1至t=106的动态仿真演化过程，在运动到14步时得到了基本的定姿态角的工作空间形状，随着步数增加，工作空间形状趋于完整，搜索结果，如图9所示。搜索得到的工作空间边界点形状精确，轮廓清晰。边界点搜索过程的数量变化情况，和4PPaR机构类似也是每步搜索到的点的个数层锯齿状增加，每步得到的边界点在380个以上，同时累计的增加呈直线趋势，在110多步时，边界点的个数已超过4.5万个，如图10所示。更换了机构也可以能得到边界轮廓清晰工作空间，能实现支链复杂的空间并联机构工作空间搜索。表明该算法的适应性能力较强。

图10 机构的搜索过程边界点的变化过程Fig.10 Searching Process of Mechanism Changing Process of Boundary Point

4.3 3D格子气偏随机模型动态演化搜索法特点

从上述算例分析，可以发现该算法具有如下特点：（1）并行计算，能够通过多个点的随机运动，同时找到一定数量的边界点，并行搜索能力突出；（2）适应性强，可适应不同的工作空间的搜索，只需更换判断工作空间约束条件就能够完成支链结构复杂的并联机构的可达工作空间与定姿态空间的搜索；（3）动态特性，粒子动态可视化运动，能够观察完整的搜索全过程，通过粒子移动搜索可以发现机构是否存在空洞的情况；（4）扩展性能，可以设定特定的搜索运动模式，即控制粒子的运动方向，有6（上下左右前后）或者18、26个运动方向（上下左右前后及其组合），方向越多，同时搜索的点的数量越多，效率越高。这里所采用邻域个数26，其速度运动方向有26个，以得到更加高效的搜索能力。

5 算法参数对搜索能力的影响分析

以4PPaR机构为例，按照表1中的结构参数，通过改变均匀格子气偏随机模型动态演化搜索法中的粒子数量、移动速度以及粒子位置等，研究这些参数对算法搜索能力的影响。

5.1 粒子个数对搜索能力的影响分析

观察对边界点搜寻个数的影响，粒子均匀布满搜索空间，粒子运动速度为1，运动方向为26 个方向，粒子数目按照1000、2000、6000个进行仿真，运算结果，如图11所示。

图11 不同粒子数量得到的搜索边界点个数变化曲线Fig.11 Change Curve of the Number of Search Boundary Points Obtained by Different Particle Numbers

线性拟合方程为：1000个粒子：y_1000=88.7x-283；2000个粒子：y_2000=205x-548；6000个粒子：y_6000=589x-1710。

根据图11 分析并结合拟合方程可得：数量为1000：2000：6000的粒子数量其边界点增长速度比为1：2.3：6.6，基本和粒子数量成正比，表明粒子数目越多，其并行计算能力越强，找到边界点的数量就越多，考虑实际时，可以选择适宜的粒子数量。

5.2 粒子移动速度对搜索能力的影响分析

观察粒子移动速度对边界点搜寻个数的影响，粒子均匀布满搜索空间，粒子运动速度为1、2和3，运动方向为26个方向，粒子数目按照2000个进行计算仿真，运算结果，如图12所示。

图12 不同粒子移动速度得到的搜索边界点个数变化曲线Fig.12 The Curve of the Number of Searching Boundary Points Obtained by Different Particle Velocity

线性拟合方程为：粒子移动速度为1：y_1=205x-548；粒子移动速度为2：y_2=212x-579；粒子移动速度为3：y_3=178x-489。按照图12和拟合方程分析可得：粒子移动速度为1：2：3的粒子数量其边界点增长速度比为1：1.03：0.87，呈先上升后下降的趋势。移动速度过快，映射到工作空间的坐标变化值越大，较难找到边界值的一个比较窄的范围，因此，选择移动速度为2较为适宜。

5.3 粒子初始运动位置对搜索能力的影响分析

观察对粒子初始位置对搜寻边界点个数的影响，粒子集中于工作空间内，粒子运动速度为2，运动方向为26个方向，粒子数目按照2000个进行仿真，仿真的动态过程图，如图13所示。相应的粒子处在不同位置得到的搜索边界点个数变化曲线，如图14所示。

图13 粒子位置集中于工作空间的搜索动态过程图Fig.13 Dynamic Process Diagram of Particle Location in Workspace

图14 粒子处在不同位置得到的搜索边界点个数变化曲线Fig.14 The Curve of the Number of Searching Boundary Points at Different Positions of Particles.

线性拟合方程为：集中位置：y_集中=614x+3097；均匀位置：y_均匀=212x-579。

根据图13分析可得：粒子运动呈现扩散趋势，在146步时已基本扩散至整个空间；初始位置在靠近工作空间时，相当于提高了粒子密度，在动态搜索过程，得到更强的搜索能力，根据图14分析可得：粒子集中分布与粒子均匀分布边界点增长速度比为2.9。因此，可以根据机构支链分布特点，改变初始粒子位置，让粒子位置靠近工作空间来实现高效搜索。