一种基于最大实体状态的同轴度快速评定方法与流程

本发明属于精密计量与计算机应用领域，具有涉及一种基于最大实体状态的同轴度快速评定方法，可用于被测圆柱体有尺寸要求、其轴线的同轴度公差有最大实体要求并且其基准要素有尺寸要求和最大实体要求的圆柱形几何产品的同轴度误差合格性检测和评定，并为加工工艺的改进提供指导。

背景技术：

尺寸误差、形位误差（形状误差和位置误差的简称）直接影响产品质量、装配及其使用寿命，快速、准确地计算零件误差，具有重要的意义。尺寸公差（公差即误差的允许范围）和形位公差之间的关系称为公差原则，其中，最大实体要求是体现零件可装配性的一种公差原则。

国家标准GB/T16671-2009规定了被测圆柱体有尺寸要求、其轴线的同轴度公差有最大实体要求并且其基准要素有尺寸要求和最大实体要求的情况，该情况规定了：1、被测圆柱体的最大实体状态；2、被测圆柱体的局部尺寸的范围；3、被测圆柱体的最大实体状态与基准圆柱体的最大实体状态之间的方位关系；4、基准圆柱体的最大实体状态及其局部尺寸的范围，补充解释了最小实体状态下的被测圆柱体的轴线与最小实体状态下的基准圆柱体的轴线之间可能的同轴度误差，但并未给出具体的计算公式。

因此，为了判断介于最大实体尺寸和最小实体尺寸之间（尺寸误差的合格性检测方法在国家标准GB/T3177、GB/T1958、GB/T18779.1、GB/T18779.2中有规定，不属于本发明的范畴）的零件的上述同轴度的合格性，国家标准GB/T1958-2004给出了利用高精度的、尺寸恒定的通规、止规来检测圆柱体的上述同轴度是否合格的方法。然而，高精度的、尺寸恒定的通规、止规的制造成本较高，并且测量不同尺寸的圆柱体零件需要使用不同的通规、止规，进一步增加了测量成本。

被测圆柱体有尺寸要求、其轴线的同轴度有最大实体要求并且其基准要素没有最大实体要求的情况下，在精密计量与计算机应用领域可以通过三坐标测量机测量被测圆柱体和基准圆柱体上的测点，计算被测圆柱体相对于基准圆柱体的同轴度，并判断被测圆柱体的同轴度是否合格。但是还没有数学方法来来评定被测圆柱体有尺寸要求、其轴线的同轴度公差有最大实体要求并且其基准要素有尺寸要求和最大实体要求的零件的合格性。

技术实现要素：

本发明所要解决的技术问题是提供一种基于最大实体状态的同轴度快速评定方法，其不仅实现了被测圆柱体有尺寸要求、其轴线的同轴度公差有最大实体要求并且其基准要素有尺寸要求和最大实体要求的零件的合格性检测和评定，而且算法稳定性好、计算效率高，可以推广应用于其它被测要素有尺寸要求、其定向定位公差有最大实体要求并且其基准要素有尺寸要求和最大实体要求的零件的合格性检测和评定中。

为解决上述问题，本发明是通过以下方案实现的：

步骤1：获取被测圆柱体b、基准圆柱体A及它们之间的几何设计参数；如果被测圆柱体b的同轴度公差及基准圆柱体A都有最大实体要求，并且基准圆柱体A只有尺寸公差——可以应用包容原则，那么跳转到步骤2，否则结束本快速评定方法，并给出结论“被测圆柱体的同轴度公差不能用本方法评定”。

所述的被测圆柱体b的几何设计参数包括：是孔要素还是轴要素、名义直径D_b、名义长度L_b、轴的上偏差es_b或孔的上偏差ES_b、轴的下偏差ei_b或孔的下偏差EI_b、同轴度公差T_{b, AM, coa}、同轴度公差是否标注最大实体要求、同轴度公差的基准圆柱体A是否标注最大实体状态。

所述的被测圆柱体b的几何设计参数包括：是孔要素还是轴要素、名义直径D_b、名义长度L_b、轴的上偏差es_b或孔的上偏差ES_b、轴的下偏差ei_b或孔的下偏差EI_b、同轴度公差T_{b, AM, coa}、同轴度公差是否标注最大实体要求、同轴度公差的基准圆柱体A是否标注最大实体状态。

所述的基准圆柱体A的几何设计参数包括：是孔要素还是轴要素、名义直径D_A、名义长度L_A、轴的上偏差es_A或孔的上偏差ES_A、轴的下偏差ei_A或孔的下偏差EI_A、尺寸公差是否应用包容原则、其它几何公差。

所述的被测圆柱体b与基准圆柱体A之间的几何设计参数是被测圆柱体b的几何中心与基准圆柱体A的几何中心之间的名义距离L_Ab。

步骤2：获取实际被测圆柱体C_b、实际基准圆柱体C_A的测量数据；评价实际基准圆柱体C_A的尺寸误差和实际被测圆柱体C_b的尺寸误差，如果上述误差都合格，跳转到步骤3，否则结束本快速评定方法，并给出结论“实际基准圆柱体C_A和/或实际被测圆柱体C_b的其它误差不合格”。

步骤3：将实际被测圆柱体b和实际基准圆柱体A的测量数据分别拟合为圆柱体CC_b和圆柱体CC_A，分别计算圆柱体CC_b和CC_A的直径d_b和d_A，计算圆柱体CC_b相对于圆柱体CC_A的三维空间姿态v_{b, A}。

步骤4：计算基准圆柱体A最大实体边界圆柱体CC_AM的直径D_AM，计算圆柱体CC_A相对于圆柱体CC_AM的三维空间姿态v_A的变动域Ω_dvA；对于任意v_A∈Ω_dvA，圆柱体CC_A不超出基准圆柱体A最大实体边界圆柱体CC_AM；对于任意不超出基准圆柱体A最大实体边界圆柱体CC_AM的圆柱体CC_A，其相对于圆柱体CC_AM的三维空间姿态v_A∈Ω_dvA；如图3所示。

步骤5：以v_A为变量，以圆柱体CC_b相对于圆柱体CC_AM的同轴度误差t_{b, AM, coa}=f(v_A)为目标函数，进行目标优化，得到圆柱体CC_b相对于圆柱体CC_AM的最小同轴度误差t_{b, AM, coa, min}；如图3所示。

步骤6：根据被测圆柱体b的最大实体实效尺寸D_bMV和圆柱体CC_b的直径d_b、最小同轴度误差t_{b, AM, min}，判断实际被测圆柱体b的同轴度误差是否合格。

为了便于将实际零件与设计图形对应起来，本发明可以将以下步骤具体化为：

步骤2中所述实际被测圆柱体b、实际基准圆柱体A的测量数据是在测量空间直角坐标系中测量的，并且包括以下四个测点数据集：

实际基准圆柱体A的底面上的测点P_{measure, A, under, i}的测点数据p_{measure, A, under, i}(x_{measure, A, under, i}, y_{measure, A, under, i}, z_{measure, A, under, i})，i=1, 2 … I，I为测点数目且为正整数，所有的测点数据p_{measure, A, under, i}(x_{measure, A, under, under, i}, y_{measure, A, under, i}, z_{measure, A, under, i})形成测点数据集{p_{measure, A, under, i}}。

实际基准圆柱体A的侧面上的测点P _{measure, A, n}的测点数据p_{measure, A, n}(x_{measure, A, n}, y_{measure, A, n}, z_{measure, A, n})，n=1, 2 … N，N为测点数目且为正整数，所有的测点数据p_{measure, A, n}(x_{measure, A, n}, y_{measure, A, n}, z_{measure, A, n})形成测点数据集{p_{measure, A, n}}。

实际被测圆柱体b的底面上的测点P_{measure, b, under, j}的测点数据p_{measure, b, under, j}(x_{measure, b, under, j}, y_{measure, b, under, j}, z_{measure, b, under, j})，j=1, 2 … J，J为测点数目且为正整数，所有的测点数据p_{measure, b, under, j}(x_{measure, b, under, under, j}, y_{measure, b, under, j}, z_{measure, b, under, j})形成测点数据集{p_{measure, b, under, j}}。

实际被测圆柱体b的侧面上的测点P_{measure, b, m}的测点数据p_{measure, b, m}(x_{measure, b, m}, y_{measure, b, m}, z_{measure, b, m})，m=1, 2 … M，M为测点数目且为正整数，所有的测点数据p_{measure, b, m}(x_{measure, b, m}, y_{measure, b, m}, z_{measure, b, m})形成测点数据集{p_{measure, b, m}}。

为了简化实际零件与设计图形对应过程，可以尽量减少实际被测圆柱体和实际基准圆柱体的底面测点，此时，本发明的以下步骤具体化为：

步骤2中，实际基准圆柱体A的底面上的测点数I=2，并且两个测点P_{measure, A, under, 1}、P_{measure,A,under,2}分别在实际基准圆柱体A的两个底面上。

实际被测圆柱体b的底面上的测点数J=2，并且两个测点P_{measure, b, under, 1}、P_{measure,b,under,2}分别在实际被测圆柱体b的两个底面上。

步骤3中，将实际基准圆柱体A的测点数据集{p_{measure, A, n}}拟合为圆柱体CC_A的侧面，P_{measure, A, under, 1}、P_{measure,A,under,2}连线的中点P_{measure, A, under, 1/2}在圆柱体CC_A的两个底面的对称平面PL_{A, 1/2}上，圆柱体CC_A的长度等于基准圆柱体A的名义长度L_A。

为描述圆柱体CC_b相对于圆柱体CC_A的三维空间姿态v_{b, A}，建立基准圆柱体直角坐标系O_{xyz, A}，坐标系O_{xyz, A}的原点与圆柱体CC_A的几何中心O_A重合，坐标系O_{xyz, A}的一个坐标轴与圆柱体CC_A的轴线重合，该坐标轴的正方向向量与向量O_AP_{measure, b, under, 1/2}的点积不小于零，其中，P_{measure, b, under, 1/2}是P_{measure, b, under, 1}、P_{measure,b,under,2}连线的中点；定义圆柱体CC_A的单位方向向量n_CCA的长度为1、正方向与该坐标轴的正方向向量相同。

将实际被测圆柱体b的测点数据集{p_{measure, b, m}}拟合为圆柱体CC_b的侧面，并且圆柱体CC_b的几何中心O_b到平面PL_{A, 1/2}的距离等于被测圆柱体b的几何中心与基准圆柱体A的几何中心之间的名义距离L_Ab，圆柱体CC_b的长度等于被测圆柱体b的名义长度L_b。

用圆柱体CC_b的几何中心O_b在直角坐标系O_{xyz, A}中的坐标来描述圆柱体CC_b相对于圆柱体CC_A的三维空间位置t_{b, A}。

定义圆柱体CC_b的单位方向向量n_CCb的长度为1、与圆柱体CC_b的轴线平行、与向量O_AP_{measure, b, under, 1/2}的点积不小于零；计算r_{b, A} =-n_CCb×n_CCA，用r_{b, A}来描述圆柱体CC_b相对于圆柱体CC_A的三维空间方向。

圆柱体CC_b相对于圆柱体CC_A的三维空间姿态v_{b, A}=(t_{b, A}, r_{b, A})。

为了降低本发明的使用难度，本发明可以具体化为：

步骤1：获取被测圆柱体b、基准圆柱体A及它们之间的几何设计参数；如果被测圆柱体b的同轴度公差及基准圆柱体A都有最大实体要求，并且基准圆柱体A只有尺寸公差——可以应用包容原则，那么跳转到步骤2，否则结束本快速评定方法，并给出结论“被测圆柱体的同轴度公差不能用本方法评定”。

所述的被测圆柱体b的几何设计参数包括是孔要素还是轴要素、名义直径D_b、名义长度L_b、轴的上偏差es_b或孔的上偏差ES_b、轴的下偏差ei_b或孔的下偏差EI_b、同轴度公差T_{b, AM, coa}、同轴度公差是否标注最大实体要求、同轴度公差的基准圆柱体A是否标注最大实体状态。

所述的基准圆柱体A的几何设计参数包括：是孔要素还是轴要素、名义直径D_A、名义长度L_A、轴的上偏差es_A或孔的上偏差ES_A、轴的下偏差ei_A或孔的下偏差EI_A、尺寸公差是否应用包容原则、其它几何公差。

所述的被测圆柱体b与基准圆柱体A之间的几何设计参数是被测圆柱体b的几何中心与基准圆柱体A的几何中心之间的名义距离L_Ab。

步骤2：获取实际被测圆柱体b、实际基准圆柱体A的测量数据，包括以下四个测点数据集：

实际基准圆柱体A的两个测点P_{measure, A, under, 1}、P_{measure,A,under,2}分别在实际基准圆柱体A的两个底面上，两个测点的测点数据p_{measure, A, under, 1} (x_{measure, A, under, 1}, y_{measure, A, under, 1}, z_{measure, A, under, 1})、p_{measure, A, under, 2} (x_{measure, A, under, 2}, y_{measure, A, under, 2}, z_{measure, A, under, 2})形成测点数据集{p_{measure, A, under, i}}，i=1, 2；实际基准圆柱体A的侧面上的测点P _{measure, A, n}的测点数据p_{measure, A, n}(x_{measure, A, n}, y_{measure, A, n}, z_{measure, A, n})，n=1, 2 … N，N为测点数目且为正整数，所有的测点数据p_{measure, A, n}(x_{measure, A, n}, y_{measure, A, n}, z_{measure, A, n})形成测点数据集{p_{measure, A, n}}；实际被测圆柱体b的两个测点P_{measure, b, under, 1}、P_{measure,b,under,2}分别在实际被测圆柱体b的两个底面上，两个测点的测点数据p_{measure, b, under, 1} (x_{measure, b, under, 1}, y_{measure, b, under, 1}, z_{measure, b, under, 1})、p_{measure, b, under, 2} (x_{measure, b, under, 2}, y_{measure, b, under, 2}, z_{measure, b, under, 2})形成测点数据集{p_{measure, b, under, j}}，j=1, 2；实际被测圆柱体b的侧面上的测点P_{measure, b, m}的测点数据p_{measure, b, m}(x_{measure, b, m}, y_{measure, b, m}, z_{measure, b, m})，m=1, 2 … M，M为测点数目且为正整数，所有的测点数据p_{measure, b, m}(x_{measure, b, m}, y_{measure, b, m}, z_{measure, b, m})形成测点数据集{p_{measure, b, m}}。

评价实际基准圆柱体A和实际被测圆柱体b的尺寸误差是否合格，如果上述误差都合格，跳转到步骤3，否则结束本快速评定方法，并给出结论“实际基准圆柱体A和/或实际被测圆柱体b的其它误差不合格”。

步骤3：计算p_{measure, A, under, 1/2}=( p_{measure, A, under, 1}+ p_{measure, A, under, 2})/2。

将步骤2中获取的四个测点数据集进行坐标变换，得到四个粗略平移数据集{p_{0, A, under, i}(x_{0, A, under, i}, y_{0, A, under, i}, z_{0, A, under, i})| p_{0, A, under, i} =p_{measure, A, under, i} -p_{measure, A, under, 1/2}，i=1, 2}、{p_{0, A, n}(x_{0, A, n}, y_{0, A, n}, z_{0, A, n})|p_{0, A, n} = p_{measure, A, n} -p_{measure, A, under, 1/2}，n=1, 2 … N }、{ p_{0, b, under, j}(x_{0, b, under, j}, y_{0, b, under, j}, z_{0, b, under, j})| p_{0, b, under, j} = p_{measure, b, under, j} -p_{measure, A, under, 1/2}，，j=1, 2}、{ p_{0, b, m}(x_{0, b, m}, y_{0, b, m}, z_{0, b, m})| p_{0, b, m} = p_{measure, b, m} -p_{measure, A, under, 1/2}，m=1, 2 … M }。

计算p_{0, b, under, 1/2}(x_{0, b, under, 1/2}, y_{0, b, under, 1/2}, z_{0, b, under, 1/2})=(p_{0, b, under, 1}+ p_{0, b, under, 2})/2。

解目标优化问题1：

s.t.

解得实际基准圆柱体A的拟合圆柱体CC_A的直径d_A=|min d_{1, A}|和对应的(x₀, y₀, α₀, β₀)的值(x_0,min, y_0,min, α_0,min, β_0,min)。

将粗略平移数据集{ p_{0, b, m}(x_{0, b, m}, y_{0, b, m}, z_{0, b, m})}进行如下坐标变换，m=1, 2 … M：

得到实际被测圆柱所有测点相对于实际基准圆柱体的坐标集{ p_{b, m}(_{b, m}, y_{b, m}, z_{b, m})}。

解目标优化问题2：

s.t.

解得实际被测圆柱体拟合圆柱体CC_b的直径d_b=|min d_{1, b}|和对应的(dx_0,b,A, dy_0,b,A, 0)、(drx_0,b,A, dry_0,b,A, 0)的值(dx_b,A, dy_b,A, 0)、( drx_b,A, dry_b,A, 0)，圆柱体CC_b相对于圆柱体CC_A的三维空间姿态v_{b, A} = (dx_b,A, dy_b,A, 0, drx_b,A, dry_b,A, 0)。

步骤4：计算基准圆柱体A最大实体边界圆柱体CC_AM的直径D_AM，当基准圆柱体A是轴时，D_AM = D_A + es_A；当基准圆柱体A是孔时，D_AM = D_A + EI_A；

计算圆柱体CC_A相对于圆柱体CC_AM的三维空间姿态v_A的变动域Ω_dvA={ v_A(dx_A, dy_A, 0, drx_A, dry_A, 0)| (dx_A ± 0.5L_Adry_A)² + ( dy_A0.5L_Adrx_A)² ≤ 0.25(D_AM - d_A)²}。

步骤5：解目标优化问题3：

s.t.

解得圆柱体CC_b相对于圆柱体CC_AM的最小同轴度误差t_{b, AM, coa, MIN}=min t_{b, AM, coa}。

步骤6：计算被测圆柱体b的最大实体实效尺寸：当被测圆柱体b是轴时，D_bMV = D_b +es_b +T_{b, AM, coa}；当被测圆柱体b是孔时，D_bMV = D_b + EI_b -T_{b, AM, coa}。

计算实际被测圆柱体b相对于基准圆柱体A的最大实体边界CC_AM的同轴度动态公差T_{b, AM, coa,Act} =| D_bMV - d_b|。

如果t_{b, AM, coa, MIN} ≤ T_{b, AM, coa,Act}，那么给出结论“实际被测圆柱体b的同轴度误差合格”，否则给出结论“实际被测圆柱体b的同轴度误差不合格”。

在零件误差分析问题中，需要求解一些目标优化问题，其特点包括：需要同时考量误差值和名义参数，而这两类数值的大小通常差别几个数量级；最优解的附近的函数值变化较缓和；等式约束通常有确定的解；通常能给出至少一个可行的解。针对这些特点，本发明提供一种多层粒子群算法来解决前述3个目标优化问题，其特征是步骤如下：

步骤01：定义多层粒子群算法的参数或使用其默认值，包括分辨率T_POS、粒子数N_P（N_P≥ 2）、内层最大迭代次数N_i,rd（N_i,rd ≥ 1）、外层最大迭代次数N_o,rd、外层最小迭代次数n_o,rd、质量权因子W、局部权因子C₁、全局权因子C₂。

步骤02：定义N_P个粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, …, p_{PSO, k,Q})，其值对应目标优化问题中的自变量(x₁, x₂, …, x_Q)，p_{PSO, k, q}的取值区间分别与x_q的取值区间相同；k=1, 2 … N_P；q=1, 2 … Q；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1的初始值或使用其默认值。

将外层粒子群算法的迭代次数设置为0。

步骤03：分别在p_{PSO, k, q}的取值区间按平均分布随机取(N_P -1)个值，构建(N_P -1)个粒子p_PSO,k；k=2, 3 … N_P；q=1, 2 … Q。

分别在p_{PSO, k, q}的取值区间按平均分布随机取N_P个值，构建N_P个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, …, v_{PSO, k,q})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … N_P；q=1, 2 … Q。

定义N_P个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … N_P。

将p_{PSO, k, q}的值分别代入目标优化问题中的变量x _q及其等式约束，并计算相应的目标优化函数值f_k=f(x₁, x₂, …, x_Q)；k=1, 2 … N_P；q=1, 2 … Q。

记录粒子p_PSO,k对应的局部最优值f_k,min=f_k，k=1, 2 … N_P。

记录全局最优值f_min=min f_k,min，并记录与min f_k,min对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤04：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … N_P。

令(x₁, x₂, …, x_Q)=p_PSO,k，并代入目标优化问题的不等式约束；如果不等式约束成立，那么将(x₁, x₂, …, x_Q)代入目标优化问题中的等式约束，并计算和更新相应的目标优化函数值f_k；如果f_k ≤ f_k,min，那么更新粒子p_PSO,k的局部最优值f_k,min=f_k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果f_k ≤ f_min，那么更新全局最优值f_min=f_k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … N_P；q=1, 2 … Q。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max f_k,min– min f_k,min| ≤ T_POS)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ T_POS)，其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … N_P；那么，转到步骤07，否则转到步骤05。

步骤05：将粒子p_PSO,k的速度v_PSO,k的值更新为W v_PSO,k + C_rand1C₁ (p_PSO,k,min - p_PSO,k) + C_rand2C₂ (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … N_P。

步骤06：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于N_i,rd，那么跳转到步骤07，否则跳转到步骤04。

步骤07：记录每次迭代得到的内层粒子群算法的全局最优值f_min,s= f_min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > n_o,rd时，令g =s -n_o,rd，判断多层粒子群算法的收敛性；如果|f_min,s– min f_{min, g} | ≤ T_POS，那么，结束多层粒子群算法并输出目标优化问题的最优解p_PSO,min和最优值f_min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于N_o,rd，那么停止外层粒子群算法并输出目标优化问题的最优解(x₁, x₂, …, x_Q)=p_PSO,min和最优值min f=f_min，否则，更新粒子p_PSO,1的初始值为p_PSO,min跳转到步骤03。

为了便于本方法的使用，前述多层粒子群算法的默认参数可以设置如下：

分辨率T_POS默认值是0.00005、粒子数N_P默认值是20、内层最大迭代次数N_i,rd默认值是100、外层最大迭代次数N_o,rd默认值是100、外层最小迭代次数n_o,rd默认值是50、质量权因子W默认值是0.5、局部权因子C₁默认值是2、全局权因子C₂默认值是2；p_PSO,1的初始值默认为零向量。

为了提高本发明的计算效率，本发明中的目标优化问题1可以用下述方法求解：

求解目标优化问题1的多层粒子群算法如下：

步骤11：定义多层粒子群算法的参数或使用其默认值，包括分辨率T_POS、粒子数N_P、内层最大迭代次数N_i,rd、外层最大迭代次数N_o,rd、外层最小迭代次数n_o,rd、质量权因子W、局部权因子C₁、全局权因子C₂。

步骤12：定义N_P个粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})，其值对应目标优化问题1中的变量(x₀, y₀, α₀, β₀)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与x₀、y₀、α₀、β₀的取值区间相同，即：p_{PSO, k,1}、p _{PSO, k,2}∈[-D_A, D_A]，p_{PSO, k,3}、p_{PSO, k,4}∈[-π, π]；k=1, 2 … N_P；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤13：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取(N_P -1)个值，构建(N_P -1)个粒子p_PSO,k；k=2, 3 … N_P。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取N_P个值，构建N_P个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … N_P。

定义N_P个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … N_P。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值分别代入目标优化问题1中的变量x₀、y₀、α₀、β₀及其等式约束，并计算相应的目标优化函数值d_1,A,k，k=1, 2 … N_P；即，

where

记录粒子p_PSO,k对应的局部最优值d_1,A,k,min=d_1,A,k，k=1, 2 … N_P。

记录全局最优值d_1,A,min=min d_1,A,k,min，并记录与min d_1,A,k,min对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤14：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … N_P。

令(x₀, y₀, α₀, β₀)=p_PSO,k，并代入目标优化问题1的不等式约束，即、、、、；如果不等式约束成立，那么将(x₀, y₀, α₀, β₀)代入目标优化问题1中的等式约束，计算并更新相应的目标优化函数值d_1,A,k；如果d_1,A,k ≤ d_1,A,k,min，那么更新粒子p_PSO,k的局部最优值d_1,A,k,min=d_1,A,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果d_1,A,k ≤ d_1,A,min，那么更新全局最优值d_1,A,min=d_1,A,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … N_P。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max d_1,A,k,min– min d_1,A,k,min| ≤ T_POS)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ T_POS)，那么，转到步骤17，否则转到步骤15；其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … N_P。

步骤15：将粒子p_PSO,k的速度v_PSO,k的值更新为W v_PSO,k + C_rand1C₁ (p_PSO,k,min - p_PSO,k) + C_rand2C₂ (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … N_P；

步骤16：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于N_i,rd，那么跳转到步骤17，否则跳转到步骤14。

步骤17：记录每次迭代得到的内层粒子群算法的全局最优值d_1,A,min,s= d_1,A,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > n_o,rd时，令g =s -n_o,rd，判断多层粒子群算法的收敛性；如果|d_1,A,min,s– min d_{1,A,min, g} | ≤ T_POS，那么，结束多层粒子群算法并输出目标优化问题1的最优解p_PSO,min和最优值d_1,A,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于N_o,rd，那么停止外层粒子群算法并输出目标优化问题1的最优解(x₀, y₀, α₀, β₀)=p_PSO,min和最优值min d_{1, A}=d_1,A,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min跳转到步骤13。

为了提高本发明的计算效率，本发明中的目标优化问题2可以用下述方法求解：

求解目标优化问题2的多层粒子群算法如下：

步骤21：定义多层粒子群算法的参数或使用其默认值，包括分辨率T_POS、粒子数N_P、内层最大迭代次数N_i,rd、外层最大迭代次数N_o,rd、外层最小迭代次数n_o,rd、质量权因子W、局部权因子C₁、全局权因子C₂。

步骤22：定义粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})的值对应目标优化问题2中的变量(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与dx_0,b,A、dy_0,b,A、drx_0,b,A、dry_0,b,A的取值区间相同，k=1, 2 … N_P；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤23：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取(N_P -1)个值，构建(N_P -1)个粒子p_PSO,k；k=2, 3 … N_P。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取N_P个值，构建N_P个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … N_P。

定义N_P个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … N_P。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值代入目标优化问题2中的变量dx_0,b,A、dy_0,b,A、drx_0,b,A、dry_0,b,A及其等式约束，并计算相应的目标优化函数值d_1,b,k，k=1, 2 … N_P。

记录粒子p_PSO,k的局部最优值d_1,b,k,min=d_1,b,k，k=1, 2 … N_P。

记录全局最优值d_1,b,min=min d_1,b,k,min，并记录与min d_1,b,k,min对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤24：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … N_P。

令(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)=p_PSO,k，并代入目标优化问题2的不等式约束，如果不等式约束成立，那么将(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)代入目标优化问题2的等式约束，计算并更新相应的目标优化函数值d_1,b,k；如果d_1,b,k ≤ d_1,b,k,min，那么更新粒子p_PSO,k的局部最优值d_1,b,k,min=d_1,b,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果d_1,b,k ≤ d_1,b,min，那么更新全局最优值d_1,b,min=d_1,b,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … N_P。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max d_1,b,k,min– min d_1,b,k,min| ≤ T_POS)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ T_POS)，那么，转到步骤27，否则转到步骤25；其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … N_P。

步骤25：将粒子p_PSO,k的速度v_PSO,k的值更新为W v_PSO,k + C_rand1C₁ (p_PSO,k,min - p_PSO,k) + C_rand2C₂ (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … N_P；

步骤26：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于N_i,rd，那么停止内层粒子群算法并跳转到步骤27，否则跳转到步骤24。

步骤27：记录每次迭代得到的外层粒子群算法的全局最优值d_1,b,min,s= d_1,b,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > n_o,rd时，令g =s -n_o,rd，判断多层粒子群算法的收敛性；如果|d_1,b,min,s– min d_{1,b,min, g} | ≤ T_POS，那么，结束多层粒子群算法并输出目标优化问题2的最优解p_PSO,min和最优值d_1,b,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于N_o,rd，那么停止外层粒子群算法并输出目标优化问题2的最优解(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)=p_PSO,min和最优值min d_{1, b}=d_1,b,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min并跳转到步骤23。

为了提高本发明的计算效率，本发明中的目标优化问题3可以用下述方法求解：

求解目标优化问题3的多层粒子群算法如下：

步骤31：定义多层粒子群算法的参数或使用其默认值，包括分辨率T_POS、粒子数N_P、内层最大迭代次数N_i,rd、外层最大迭代次数N_o,rd、外层最小迭代次数n_o,rd、质量权因子W、局部权因子C₁、全局权因子C₂。

步骤32：定义粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})的值对应目标优化问题3中的变量(dx_A, dy_A, drx_A, dry_A)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与dx_A、dy_A、drx_A、dry_A的取值区间相同，k=1, 2 … N_P；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤33：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间中按平均分布随机取(N_P -1)个值，构建(N_P -1)个粒子p_PSO,k；k=2, 3 … N_P。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间中按平均分布随机取N_P个值，构建N_P个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … N_P。

定义N_P个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … N_P。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值代入目标优化问题3中的变量dx_A、dy_A、drx_A、dry_A及其等式约束，并计算相应的目标优化函数值t_b,AM,coa,k，k=1, 2 … N_P；

记录粒子p_PSO,k的局部最优值t_{b,AM,coa,k,min}=t_b,AM,coa,k，k=1, 2 … N_P。

记录全局最优值t_b,AM,coa,min=min t_{b,AM,coa,k,min}，并记录与min t_{b,AM,coa,k,min}对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤34：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … N_P。

令(dx_A, dy_A, drx_A, dry_A)=p_PSO,k，并代入目标优化问题3的不等式约束，如果不等式约束成立，那么将(dx_A, dy_A, drx_A, dry_A)代入目标优化问题3的等式约束，计算并更新相应的目标优化函数值t_b,AM,coa,k；如果t_b,AM,coa,k ≤ t_{b,AM,coa,k,min}，那么更新粒子p_PSO,k的局部最优值t_{b,AM,coa,k,min}=t_b,AM,coa,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果t_b,AM,coa,k ≤ t_b,AM,coa,min，那么更新全局最优值t_b,AM,coa,min=t_b,AM,coa,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … N_P。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max t_{b,AM,coa,k,min}– min t_{b,AM,coa,k,min}| ≤ T_POS)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ T_POS)，那么，转到步骤37，否则转到步骤35；其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … N_P。

步骤35：将粒子p_PSO,k的速度v_PSO,k的值更新为W v_PSO,k + C_rand1C₁ (p_PSO,k,min - p_PSO,k) + C_rand2C₂ (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … N_P。

步骤36：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于N_i,rd，那么停止内层粒子群算法并跳转到步骤37，否则跳转到步骤34。

步骤37：记录每次迭代得到的外层粒子群算法的全局最优值t_{b,AM,coa,min,s}=t_b,AM,coa,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > n_o,rd时，令g =s -n_o,rd，判断多层粒子群算法的收敛性；如果|t_{b,AM,coa,min,s}– min t_{b,AM,coa,min, g} | ≤ T_POS，那么，结束多层粒子群算法并输出目标优化问题3的最优解p_PSO,min和最优值t_b,AM,coa,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于N_o,rd，那么停止外层粒子群算法并输出目标优化问题3的最优解(dx_A, dy_A, drx_A, dry_A)=p_PSO,min和最优值min t_b,AM,coa=t_b,AM,coa,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min并跳转到步骤33。

附图说明

图1，本发明的基本方法的流程图。

图2，本发明中求解目标优化问题的多层粒子群算法的基本方法的流程图。

图3，步骤4及步骤5中零件相对于最大实体边界的运动。

图4，实验对象的几何设计图。

具体实施方式

实验实施例：

被测轴有尺寸要求、其轴线的同轴度公差有最大实体要求并且其基准轴有尺寸要求和最大实体要求的台阶轴零件的同轴度误差合格性检测和评定：

步骤1：获取如图4所示的被测圆柱体b、基准圆柱体A及它们之间的几何设计参数（长度单位为毫米，角度单位为度，弧度单位为1）：

所述的被测圆柱体b的几何设计参数包括：是轴要素、名义直径D_b=10、名义长度L_b=15、轴的上偏差es_b=0、轴的下偏差ei_b=-0.03、同轴度公差T_{b, AM, coa}=0.02、同轴度公差标注最大实体要求、同轴度公差的基准圆柱体A标注最大实体状态。

所述的基准圆柱体A的几何设计参数包括：是轴要素、名义直径D_A=20、名义长度L_A=30、轴的上偏差es_A=0、轴的下偏差ei_A=-0.05、尺寸公差应用包容原则，没有几何公差要求。

所述的被测圆柱体b与基准圆柱体A之间的几何设计参数是被测圆柱体b的几何中心与基准圆柱体A的几何中心之间的名义距离L_Ab=22.5。

被测圆柱体b的同轴度公差及基准圆柱体A都有最大实体要求，并且基准圆柱体A只有尺寸公差要求并应用包容原则，跳转到步骤2。

步骤2：获取实际被测圆柱体b、实际基准圆柱体A的测量数据，包括以下四个测点数据集，如表1所示：

实际基准圆柱体A的两个测点P_{measure, A, under, 1}、P_{measure,A,under,2}分别在实际基准圆柱体A的两个底面上，两个测点的测点数据形成测点数据集{p_{measure, A, under, i}}，i=1, 2；实际基准圆柱体A的侧面上的测点P _{measure, A, n}分布在3层圆周上，每层圆周3个点，其测点数据形成测点数据集{p_{measure, A, n}}，n=1, 2 … 9；实际被测圆柱体b的两个测点P_{measure, b, under, 1}、P_{measure,b,under,2}分别在实际被测圆柱体b的两个底面上，两个测点的测点数据形成测点数据集{p_{measure, b, under, j}}，j=1, 2；实际被测圆柱体b的侧面上的测点P_{measure, b, m}分布在3层圆周上，每层圆周3个点，其测点数据形成测点数据集{p_{measure, b, m}，m=1, 2 … 9}。

将3层圆周上的测点P _{measure, A, n}分别拟合为3个圆，其直径分别为19.966、19.966、19.969，判定实际基准圆柱体A的尺寸误差合格；将3层圆周上的测点p_{measure, b, m}分别拟合为3个圆，其直径分别为9.980、9.979、9.983，判定实际被测圆柱体b的尺寸误差合格；跳转到步骤3。

步骤3：计算p_{measure, A, under, 1/2}=( p_{measure, A, under, 1}+ p_{measure, A, under, 2})/2=(254.156, 354.498, -656.454)。

将步骤2中获取的四个测点数据集进行坐标变换，得到四个粗略平移数据集{p_{0, A, under, i}(x_{0, A, under, i}, y_{0, A, under, i}, z_{0, A, under, i})| p_{0, A, under, i} =p_{measure, A, under, i} -p_{measure, A, under, 1/2}，i=1, 2}，如表1所示；{p_{0, A, n}(x_{0, A, n}, y_{0, A, n}, z_{0, A, n})| p_{0, A, n} = p_{measure, A, n} -p_{measure, A, under, 1/2}，n=1, 2 … N }，如表1所示；{p_{0, b, under, j}(x_{0, b, under, j}, y_{0, b, under, j}, z_{0, b, under, j})| p_{0, b, under, j} =p_{measure, b, under, j} -p_{measure, A, under, 1/2}，，j=1, 2}，如表1所示；{ p_{0, b, m}(x_{0, b, m}, y_{0, b, m}, z_{0, b, m})| p_{0, b, m} = p_{measure, b, m} -p_{measure, A, under, 1/2}，m=1, 2 … M }，如表1所示。

计算p_{0, b, under, 1/2}=( p_{0, b, under, 1}+ p_{0, b, under, 2})/2=(-3.670, -0.884, 22.512)。

解目标优化问题1：

s.t.

步骤11：使用多层粒子群算法的默认参数，包括分辨率T_POS=0.00005、粒子数N_P=20、内层最大迭代次数N_i,rd=100、外层最大迭代次数N_o,rd =100、外层最小迭代次数n_o,rd=50、质量权因子W=0.5、局部权因子C₁=2、全局权因子C₂=2。

步骤12：定义20个粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})，其值对应目标优化问题1中的变量(x₀, y₀, α₀, β₀)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与x₀、y₀、α₀、β₀的取值区间相同，即：p_{PSO, k,1}、p _{PSO, k,2}∈[-20, 20]，p_{PSO, k,3}、p_{PSO, k,4}∈[-π, π]；k=1, 2 … 20；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤13：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取19个值，构建19个粒子p_PSO,k；k=2, 3 … 20。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取20个值，构建20个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … 20。

定义20个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … 20。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值分别代入目标优化问题1中的变量x₀、y₀、α₀、β₀及其等式约束，并计算相应的目标优化函数值d_1,A,k，k=1, 2 … 20；即，

where

记录粒子p_PSO,k对应的局部最优值d_1,A,k,min=d_1,A,k，k=1, 2 … 20。

记录全局最优值d_1,A,min=min d_1,A,k,min，并记录与min d_1,A,k,min对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤14：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … 20。

令(x₀, y₀, α₀, β₀)=p_PSO,k，并代入目标优化问题1的不等式约束，即0≤-x_{0,b,under,1/2} cosα₀ sinβ₀ + y_{0,b,under,1/2} sinα₀ + z_{0,b,under,1/2} cosα₀ cosβ₀、x₀∈[-20, 20]、y₀∈[-20, 20]、α₀∈[-π, π]、β₀∈[-π, π]；如果不等式约束成立，那么将(x₀, y₀, α₀, β₀)代入目标优化问题1中的等式约束，计算并更新相应的目标优化函数值d_1,A,k；如果d_1,A,k ≤ d_1,A,k,min，那么更新粒子p_PSO,k的局部最优值d_1,A,k,min=d_1,A,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果d_1,A,k ≤ d_1,A,min，那么更新全局最优值d_1,A,min=d_1,A,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … 20。如果各粒子的局部最优值和局部最优解非常接近，即，(|max d_1,A,k,min– min d_1,A,k,min| ≤ 0.00005)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ 0.00005)，其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … 20；那么，转到步骤17，否则转到步骤15。

步骤15：将粒子p_PSO,k的速度v_PSO,k的值更新为0.5v_PSO,k +2C_rand1 (p_PSO,k,min - p_PSO,k) + 2C_rand2 (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … 20；

步骤16：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于100，那么跳转到步骤17，否则跳转到步骤14。

步骤17：记录每次迭代得到的内层粒子群算法的全局最优值d_1,A,min,s= d_1,A,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > 50时，令g =s -50，判断多层粒子群算法的收敛性；如果|d_1,A,min,s– min d_{1,A,min, g} | ≤ 0.00005，那么，结束多层粒子群算法并输出目标优化问题1的最优解p_PSO,min和最优值d_1,A,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于100，那么停止外层粒子群算法并输出目标优化问题1的最优解(x₀, y₀, α₀, β₀)=p_PSO,min和最优值min d_{1, A}=d_1,A,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min跳转到步骤13。

解得实际基准圆柱体A的拟合圆柱体CC_A的直径d_A=|min d_{1, A}|=19.969和对应的p_PSO,min的值(-0.0005, -0.0026, 0.0004, 0.0001)。

将粗略平移数据集{ p_{0, b, m}}进行如下坐标变换，m=1, 2 … 9：

得到实际被测圆柱所有测点相对于实际基准圆柱体的坐标集{ p_{b, m}} ，如表1所示。

解目标优化问题2：

s.t.

步骤21：使用多层粒子群算法的默认参数，包括分辨率T_POS=0.00005、粒子数N_P=20、内层最大迭代次数N_i,rd=100、外层最大迭代次数N_o,rd =100、外层最小迭代次数n_o,rd=50、质量权因子W=0.5、局部权因子C₁=2、全局权因子C₂=2。

步骤22：定义粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})的值对应目标优化问题2中的变量(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与dx_0,b,A、dy_0,b,A、drx_0,b,A、dry_0,b,A的取值区间相同，k=1, 2 … 20；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤23：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取19个值，构建19个粒子p_PSO,k；k=2, 3 … 20。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间按平均分布随机取20个值，构建20个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … N_P。

定义20个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … 20。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值代入目标优化问题2中的变量dx_0,b,A、dy_0,b,A、drx_0,b,A、dry_0,b,A及其等式约束，并计算相应的目标优化函数值d_1,b,k，k=1, 2 … 20。

记录粒子p_PSO,k的局部最优值d_1,b,k,min=d_1,b,k，k=1, 2 … 20。

记录全局最优值d_1,b,min=min d_1,b,k,min，并记录与min d_1,b,k,min对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤24：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … 20。

令(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)=p_PSO,k，并代入目标优化问题2的不等式约束，如果不等式约束成立，那么将(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)代入目标优化问题2的等式约束，计算并更新相应的目标优化函数值d_1,b,k；如果d_1,b,k ≤ d_1,b,k,min，那么更新粒子p_PSO,k的局部最优值d_1,b,k,min=d_1,b,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果d_1,b,k ≤ d_1,b,min，那么更新全局最优值d_1,b,min=d_1,b,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … 20。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max d_1,b,k,min– min d_1,b,k,min| ≤ 0.00005)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤0.00005)，其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … 20；那么，转到步骤27。

步骤25：将粒子p_PSO,k的速度v_PSO,k的值更新为0.5v_PSO,k +2C_rand1 (p_PSO,k,min - p_PSO,k) + 2C_rand2 (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … 20；

步骤26：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于100，那么停止内层粒子群算法并跳转到步骤27，否则跳转到步骤24。

步骤27：记录每次迭代得到的外层粒子群算法的全局最优值d_1,b,min,s= d_1,b,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > 50时，令g =s -50，判断多层粒子群算法的收敛性；如果|d_1,b,min,s– min d_{1,b,min, g} | ≤0.00005，那么，结束多层粒子群算法并输出目标优化问题2的最优解p_PSO,min和最优值d_1,b,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于100，那么停止外层粒子群算法并输出目标优化问题2的最优解(dx_0,b,A, dy_0,b,A, drx_0,b,A, dry_0,b,A)=p_PSO,min和最优值min d_{1, b}=d_1,b,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min并跳转到步骤23。

解得实际被测圆柱体拟合圆柱体CC_b的直径d_b=|min d_{1, b}|=9.9833和对应的p_PSO,min的值(-0.0031, -0.0245, -0.0007, -0.0001)，圆柱体CC_b相对于圆柱体CC_A的三维空间姿态v_{b, A} = (-0.0031, -0.0245, 0, -0.0007, -0.0001, 0)。

步骤4：计算基准圆柱体A最大实体边界圆柱体CC_AM的直径D_AM，D_AM = D_A + es_A= 20。

计算圆柱体CC_A相对于圆柱体CC_AM的三维空间姿态v_A的变动域Ω_dvA={ v_A(dx_A, dy_A, 0, drx_A, dry_A, 0)| =(dx_A ± 15dry_A)² + ( dy_A15drx_A)² ≤ 0.0156²}。

步骤5：解目标优化问题3：

s.t.

步骤31：使用多层粒子群算法的默认参数，包括分辨率T_POS=0.00005、粒子数N_P=20、内层最大迭代次数N_i,rd=100、外层最大迭代次数N_o,rd =100、外层最小迭代次数n_o,rd=50、质量权因子W=0.5、局部权因子C₁=2、全局权因子C₂=2。

步骤32：定义粒子p_PSO,k(p_{PSO, k,1}, p _{PSO, k,2}, p_{PSO, k,3}, p_{PSO, k,4})的值对应目标优化问题3中的变量(dx_A, dy_A, drx_A, dry_A)，p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间分别与dx_A、dy_A、drx_A、dry_A的取值区间相同，k=1, 2 … 20；所有粒子组成粒子集{p_PSO,k}；设置一个粒子p_PSO,1(p _PSO,1,1, p _PSO,1,2, p_PSO,1,3, p_PSO,1,4)的初始值为(0, 0, 0, 0)。

将外层粒子群算法的迭代次数设置为0。

步骤33：分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间中按平均分布随机取19个值，构建19个粒子p_PSO,k；k=2, 3 … 20。

分别在p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的取值区间中按平均分布随机取20个值，构建20个粒子p_PSO,k的初始速度v_{PSO, k} (v_{PSO, k,1}, v _{PSO, k,2}, v_{PSO, k,3}, v_{PSO, k,4})；所有粒子的速度组成速度集{v_PSO,k}；k=1, 2 … 20。

定义20个粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；k=1, 2 … 20。

将p_{PSO, k,1}、p _{PSO, k,2}、p_{PSO, k,3}、p_{PSO, k,4}的值代入目标优化问题3中的变量dx_A、dy_A、drx_A、dry_A及其等式约束，并计算相应的目标优化函数值t_b,AM,coa,k，k=1, 2 … 20；

记录粒子p_PSO,k的局部最优值t_{b,AM,coa,k,min}=t_b,AM,coa,k，k=1, 2 … 20。

记录全局最优值t_b,AM,coa,min=min t_{b,AM,coa,k,min}，并记录与min t_{b,AM,coa,k,min}对应的粒子的值为全局最优解p_PSO,min。

将内层粒子群算法的迭代次数设置为0。

步骤34：将粒子p_PSO,k的值更新为p_PSO,k +v_{PSO, k}，k=1, 2 … 20。

令(dx_A, dy_A, drx_A, dry_A)=p_PSO,k，并代入目标优化问题3的不等式约束，如果不等式约束成立，那么将(dx_A, dy_A, drx_A, dry_A)代入目标优化问题3的等式约束，计算并更新相应的目标优化函数值t_b,AM,coa,k；如果t_b,AM,coa,k ≤ t_{b,AM,coa,k,min}，那么更新粒子p_PSO,k的局部最优值t_{b,AM,coa,k,min}=t_b,AM,coa,k，并更新粒子p_PSO,k的局部最优解为p_PSO,k,min=p_PSO,k；如果t_b,AM,coa,k ≤ t_b,AM,coa,min，那么更新全局最优值t_b,AM,coa,min=t_b,AM,coa,k，并更新全局最优解为p_PSO,min=p_PSO,k；k=1, 2 … 20。

如果各粒子的局部最优值和局部最优解非常接近，即，(|max t_{b,AM,coa,k,min}– min t_{b,AM,coa,k,min}| ≤ 0.00005)且(|p_{PSO,k,min, mean} – p_PSO,k,min| ≤ 0.00005)，其中，p_{PSO,k,min, mean}为p_PSO,k,min的算术平均值，k=1, 2 … 20；那么，转到步骤37。

步骤35：将粒子p_PSO,k的速度v_PSO,k的值更新为0.5v_PSO,k +2C_rand1 (p_PSO,k,min - p_PSO,k) + 2C_rand2 (p_PSO,min - p_PSO,k)，其中，C_rand1、C_rand2是在区间[0, 1]中相互独立地按平均分布随机选取的两个值； k=1, 2 … 20。

步骤36：累积一次内层粒子群算法的迭代次数；如果内层粒子群算法的迭代次数大于100，那么停止内层粒子群算法并跳转到步骤37，否则跳转到步骤34。

步骤37：记录每次迭代得到的外层粒子群算法的全局最优值t_{b,AM,coa,min,s}=t_b,AM,coa,min,其中，s是外层粒子群算法的迭代次数。

将粒子p_PSO,1的值设置为当前的全局最优解p_PSO,min。

当外层粒子群算法的迭代次数s > 50时，令g =s -50，判断多层粒子群算法的收敛性；如果|t_{b,AM,coa,min,s}– min t_{b,AM,coa,min, g} | ≤ 0.00005，那么，结束多层粒子群算法并输出目标优化问题3的最优解p_PSO,min和最优值t_b,AM,coa,min，否则，累积一次外层粒子群算法的迭代次数。

如果外层粒子群算法的迭代次数大于100，那么停止外层粒子群算法并输出目标优化问题3的最优解(dx_A, dy_A, drx_A, dry_A)=p_PSO,min和最优值min t_b,AM,coa=t_b,AM,coa,min，否则，更新粒子p_PSO,1的初始值为p_PSO,min并跳转到步骤33。

解得圆柱体CC_b相对于圆柱体CC_AM的最小同轴度误差t_{b, AM, coa, MIN}=min t_{b, AM, coa}=0.0665。

步骤6：计算被测圆柱体b的最大实体实效尺寸D_bMV = D_b + es_b +T_{b, AM, coa}=10.02。

计算实际被测圆柱体b相对于基准圆柱体A的最大实体边界CC_AM的同轴度动态公差T_{b, AM, coa,Act} =D_bMV - d_b=0.0367。

T_{b, AM, coa,Act} =0.0367≤0.0665= t_{b, AM, coa, MIN}，给出结论“实际被测圆柱体b的同轴度误差不合格”。

表1，数据集和坐标集