CN110287553A - 一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法 - Google Patents

Abstract

本发明公开了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析新方法，属于机床精度设计领域，具体涉及到多轴数控机床的空间误差建模方法以及基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法。本发明运用多体系统理论建立数控机床空间误差模型，根据蒙特卡洛模拟采样机理，对机床加工误差模型进行全局灵敏度分析，获得影响机床加工误差的关键几何误差参数，在机床设计的初期阶段，提出新的机床设计理念，为提升数控机床的加工精度以及关键几何误差参数补偿奠定了理论基础。

Description

一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析 方法

技术领域

本发明涉及一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，属于机床精度设计技术领域。

背景技术

近年来，数控机床作为制造业的工业母机，其功能在不断的完善，为了满足日益复杂形状工件精密加工的要求，五轴机床得到越来越广泛地应用。五轴机床结构复杂，涉及误差因素较多，加工成型机理也比三轴机床复杂。相比于三轴机床，五轴机床具有更好的通用性和灵活性、更高的效率和精度等优点。由于五轴数控机床各项技术的不断成熟，大大提高了复杂曲面零件的可加工性和加工精度。

数控机床的精度指标主要有加工精度、定位精度和重复定位精度，其中加工精度是数控机床追求的最终精度，反映出机械制造业的制造能力和发展水平。影响机床加工精度的因素主要包括机床各个零部件的几何误差、热误差、切削力引入误差和刀具磨损等，其中几何误差对加工精度的影响最大。机床的几何误差包括直线度误差，滚摆误差，偏摆误差，定位误差，颠摆误差和垂直度误差等。数控机床各项几何误差参数的作用相互耦合，对加工精度的影响权重各不相同。如何有效的识别出对数控机床加工精度影响较大的关键几何误差项，并且在机床设计的初期阶段有效的控制它们是有效提高机床加工精度的关键问题。

这一关键问题的解决方法分为两个步骤：

第一、基于多体系统运动学理论，建立机床的空间误差模型；

目前国内外学者已经开展了许多关于机床精度建模方法的研究，先后出现了二次关系模型法、几何建模法、误差矩阵法、刚体运动学法和多体系统理论法。基于多体系统运动学理论，将五轴机床抽象为多体系统，用拓扑结构图以及低序体阵列表来描述机床的结构和各个体之间的关联关系，分析数控机床的几何误差，建立广义坐标系，用相邻体间的特征矩阵表示位置关系，用齐次变换矩阵表示多体系统间的相互关系，最终建立机床的空间误差模型；

第二、结合空间误差模型，提出基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法；

由于拟蒙特卡洛方法的计算准确度高，收敛速度快，所以选用拟蒙特卡洛方法进行计算。将机床误差模型作为黑盒处理，所有几何误差参数在定义域内的随机采样值作为输入变量，机床加工误差为输出值，根据各输入变量的方差对输出值方差的影响，评估各项几何误差参数对机床加工误差模型的影响程度。灵敏度系数越大，说明该属性对模型输出的影响越大，根据灵敏度系数的大小对各项几何误差参数进行排序，最后筛选出灵敏度系数较大的几何误差项。本发明采用基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法来分析五轴数控机床关键几何误差参数。

本发明基于多体系统运动学理论，建立了机床的空间误差模型，然后结合空间误差模型提出了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，最后有效的辨识出了五轴数控机床的关键几何误差参数。

发明内容

本发明的目的是提供一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法。通过建立机床的空间误差模型和误差敏感度分析模型，辨识出关键几何误差参数，为提高加工精度奠定基础。

为了实现上述目的，本发明采用的技术方案为一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，本发明基于多体系统运动学理论，建立了机床的空间误差模型，然后结合空间误差模型提出了一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，最后有效的辨识出了数控机床的关键几何误差。

本方法具体包括如下步骤：

步骤一：建立数控机床的空间误差模型；

基于多体系统运动学理论，用多体系统示意图以及低序体阵列表对机床的结构进行简化，如图2和表1所示。分析数控机床的几何误差参数，建立广义坐标系，用相邻体间的特征矩阵表示各零部件之间的位置关系，用齐次变换矩阵表示多体系统间的相互关系；

步骤1.1建立数控机床的拓扑结构；

数控机床是一个多分支的复杂系统，从B₁处分为两个分支，除了B₁体外每个物体都有一个相邻的较低序体，用Lⁿ(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号，j＝1,2,3…n，n表示机床所包含典型体的个数；

表1：数控机床低序体阵列

L<sup>0</sup>(j)	1	2	3	4	5	6
							L<sup>1</sup>(j)	0	1	1	3	4	5
L<sup>2</sup>(j)	0	0	0	1	3	4
							L<sup>3</sup>(j)	0	0	0	0	1	3
L<sup>4</sup>(j)	0	0	0	0	0	1
							L<sup>5</sup>(j)	0	0	0	0	0	0

典型体的编号规则如下：

首先任选一典型体为B₁，然后沿远离B₁体的方向，依自然增长的数列依次标定每个物体的序号；

步骤1.2数控机床的几何误差分析

在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项与位置有关的误差，包括3项线位移误差和3项角位移误差，X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，因此共37项误差如表2所示；

表2：数控机床几何误差参数

步骤1.3建立数控机床的特征矩阵；

根据数控机床各部件之间的运动关系，建立各相邻体之间的变换矩阵如表3所示；

表3：相邻体间的变换矩阵

其中：[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_pe表示B_j体相对于B_i体的相对位置误差变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

[Tij]_se表示B_j体相对于B_i体的相对运动误差变换矩阵；

x表示X轴平移的距离；

y表示Y轴平移的距离；

z表示Z轴平移的距离；

a表示A轴转动的角度；

c表示C轴转动的角度；

几何误差的敏感度分析方法使用过程中，忽略除几何误差之外的所有误差因素；

步骤1.4建立机床的空间误差模型

理想情况下相邻体运动关系模型的建立；

设P点为B_j体上任意一点，P在B_i体坐标系O_i-X_iY_iZ_i中的位置矩阵表达式为；

P_ji＝[Tij]_p[Tij]_sr_j(1)

式中：P_ji为P点在坐标系O_i-X_iY_iZ_i中的位置矩阵表达式；

r_j为P点在坐标系O_j-X_jY_jZ_j中的位置矩阵表达式；

[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

有误差情况下相邻体运动关系模型的建立；

设P点为B_j体上任意一点，P在B_i体坐标系O_i-X_iY_iZ_i中的位置矩阵表达式为；

P_ji＝[Tij]_p[Tij]_pe[Tij]_s[Tij]_ser_j(2)

式中：P_ji为P点在坐标系O_i-X_iY_iZ_i中的位置矩阵表达式；

r_j为P点在坐标系O_j-X_jY_jZ_j中的位置矩阵表达式；

[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_pe表示B_j体相对于B_i体的相对位置误差变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

[Tij]_se表示B_j体相对于B_i体的相对运动误差变换矩阵；

刀具中心点在刀具坐标系中的坐标为：

r_t＝[0,0,l,1]^T (3)

l表示刀具长度；

下标t表示刀具

理想情况下刀具中心点P按“数控机床-工件”分支到惯性坐标系中的位置矩阵表达式：

理想情况下刀具中心点P按“数控机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

数控指令精密加工方程：

P_w ^I＝P_t ^I (6)

理想情况下，数控指令到工件坐标系中的位置矩阵表达式：

实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：

实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

实际情况下，数控指令到工件坐标系中的位置矩阵表达式：

则数控机床的空间误差模型表示为：

E＝r_w-r_w ^I (11)

步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；

全局灵敏度分析是一种针对系统模型的研究方法，而不是针对模型的某些特定解的分析。对机床关键几何误差因素进行灵敏度分析的本质，是要求解机床各项几何误差对机床加工误差的固有影响程度。应用基于方差的灵敏度分析方法，可以同时考虑所有几何误差因素对几何误差模型的影响。Sobol提出的基于蒙特卡洛的灵敏度分析法，是目前最为常用的灵敏度分析方法。

令I为单位向量，Iⁿ为n维单位立方空间，x∈Iⁿ，以下每一项变量的积分区间均为[0,1]。设系统方程为y＝f(x)，其中y为模型输出，x＝(x₁,x₂,...x_n)为模型的n个输入变量。f(x)的高维模型分解表示(analysis of variance，ANOVA)为公式(12)。

其中f₀＝E(y)，f_i＝E(y|x_i)-E(y)，f_ij＝E(y|x_i,x_j)-f_i-f_j-E(y)

公式(12)中，f(x)被分解为2n项。当各变量相互独立且正交时，这种分解方式唯一。

对公式(12)两边同时求方差，得：

其中V_i＝V(fi(x_i))＝V[E(y|x_i)]，V_ij＝V(f_ij(x_i,x_j))＝V(E(y|x_i,x_j))-V_i-V_j，

V_ijk＝V(f_ijk(x_i,x_j,x_k))＝V(E(y|x_i,x_j,x_k))-V_ij-V_ik-V_jk-V_i-V_j-V_k

令S_i＝V_i/V(y)，S_ij＝V_ij/V(y)，方程两边同时除以V(y)，得：

其中S_i为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；S_ij为2阶灵敏度指标，为x_i和x_j对输出方差的联合影响程度减去各自的主灵敏度指标，表示x_i和x_j的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推。S_i越大，x_i对输出方差的影响程度越大。

根据公式13，计算系统的一阶灵敏度需要计算两项参数V(y)和V[E(y|x_i)]。

设y为n个输入变量的函数

y＝f(x₁,x₂,...x_n) (15)

假设各自独立变量的联合概率密度函数为

由此可得y的期望和方差表达如下

令x_j(j＝1,2,...n)取固定值则

其中和分别为输入变量时系统输出的方差和期望。

通过x_j的概率密度函数计算的期望，可以消除其对数值的依赖。

V(y)＝E[V(y|x_j)]+V[E(y|x_j)] (23)

由此可以得出以下关系

令

以上U_j的方程可以用下式表达

F(X)由2n-1个独立变量决定。对每一项变量进行N次采样后，可以估算f和f^*输出值的数学期望。f的输出值由N×n维的输入变量采样矩阵计算。将该矩阵的第j列固定，其他数据进行重采样，可以计算f^*的输出值。根据已知的X的分布函数，构造两个N×n的随机矩阵A、B。

将矩阵B的第j列用矩阵A的第j列替代，得矩阵C_j。

将以上样本矩阵A、C_j作为输入，带入系统方程，得到输出响应

y_iA＝f(x_i1,x_i2,...x_in) (27)

对于离散变量x，U_j可由下式估计

由式计算输入变量x_i的灵敏度指标为

根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度；灵敏度系数小说明该项几何误差对机床空间误差影响较小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响较大。然后根据灵敏度系数列出误差源参数敏感度队列；根据敏感度分析结果，在机床设计初期，对相应的主要误差进行严格的限制，从而提高机床的加工精度。

与现有技术相比，本发明具有如下有益效果。

1、在机床设计的初期阶段，可为工程师在机床的设计，装配和加工方面提供指导，从根本上提高机床的加工精度。

2、本文证明了全局灵敏度分析方法可以有效地分析出机床误差模型的敏感参数，找出机床误差参数与加工误差模型存在的固有联系，并且该方法可以有效的避免敏感度分析结果失真的情况。针对某一特定型号的机床，该方法可在机床设计阶段分析出机床的敏感几何误差参数，为机床的设计和制造提供指导。

附图说明

图1为本发明方法的实施流程图

图2为多体系统示意图

图3为五轴机床的拓扑结构图

图4为五轴机床的结构示意图

图5为位置点有关的误差参数示意图

图6为位置点无关的误差参数示意图

图7为各项几何误差对机床加工精度的敏感度排序图

具体实施方式

本发明以五轴高架横梁移动龙门数控铣床为例，对上述五轴数控铣床关键性几何误差的辨识方法进行验证。

具体包括如下步骤：

步骤一：以五轴数控机床为例，建立机床的空间误差模型；

基于多体系统运动学理论，用拓扑结构图以及低序体阵列表来描述机床的结构和各个体之间的关联关系，如图3和表1所示。分析数控机床的几何误差，建立广义坐标系，用相邻体间的特征矩阵表达位置关系，用齐次变换矩阵表示多体系统间的相互关系；

步骤1.1建立五轴数控机床的拓扑结构；

该机床的结构如图4所示。包括床身、工作台、刀具、工件、X轴、Y轴、Z轴、B轴、C轴、主轴；

五轴数控机床是一个多分支的复杂系统，从B₁处分为两个分支，除了B₁体外每个物体都有一个相邻的较低序体，用Lⁿ(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号(j＝1,2,3…n)，n表示机床所包含典型体的个数；

表1：数控机床低序体阵列

L<sup>0</sup>(j)	1	2	3	4	5	6
							L<sup>1</sup>(j)	0	1	1	3	4	5
L<sup>2</sup>(j)	0	0	0	1	3	4
							L<sup>3</sup>(j)	0	0	0	0	1	3
L<sup>4</sup>(j)	0	0	0	0	0	1
							L<sup>5</sup>(j)	0	0	0	0	0	0

典型体的编号规则如下：

首先任选一典型体为B₁，然后沿远离B₁体的方向，依自然增长的数列依次标定每个物体的序号，从系统的一个分支到另一个分支，直到全部物体都标定完毕；

步骤1.2分析五轴数控机床的几何误差；

在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项误差，3项线位移误差和3项角位移误差，这些都是与位置点有关的误差，如图5所示。X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，如图6所示。因此共37项误差如表2所示；

表2：五轴数控机床几何误差参数

步骤1.3建立五轴数控机床的特征矩阵；

根据数控机床各部件之间的运动关系，可建立各相邻体之间的变换矩阵如表3所示；

表3：相邻体间的变换矩阵

本发明是一种几何误差的敏感度分析方法，使用过程中忽略除几何误差之外的所有误差因素；

步骤1.4建立机床的空间误差模型；

刀具中心点在刀具坐标系中的坐标为：

r_t＝[0,0,l,1]^T (1)

l表示刀具长度；

下标t表示刀具

理想情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：

P_w ^I＝[T12]_p[T12]_sr_w ^I (2)

理想情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

P_t ^I＝[T13]_p[T13]_s[T34]_p[T34]_s[T45]_p[T45]_s[T56]_p[T56]_sr_t ^I (3)

数控指令精密加工方程：

P_w ^I＝P_t ^I (4)

理想情况下，数控指令到工件坐标系中的位置矩阵表达式：

r_w ^I＝([T12]_p[T12]_s)^-1[T13]_p[T13]_s[T34]_p[T34]_s[T45]_p[T45]_s[T56]_p[T56]_sr_t ^I(5)

实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：

P_w＝[T12]_p[T12]_pe[T12]_s[T12]_ser_w (6)

实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

实际情况下，数控指令到工件坐标系中的位置矩阵表达式：

则机床的空间误差模型表示为：

E＝r_w-r_w ^I (9)

步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；

全局灵敏度分析是一种针对系统模型的研究方法，而不是针对模型的某些特定解的分析。对机床关键几何误差因素进行灵敏度分析的本质，是要求解机床各项几何误差对机床加工误差的固有影响程度。应用基于方差的灵敏度分析方法，可以同时考虑所有几何误差因素对几何误差模型的影响。Sobol提出的基于蒙特卡洛的灵敏度分析法，是目前最为常用的灵敏度分析方法。

令I为单位向量，In为n维单位立方空间，x∈In，以下每一项变量的积分区间均为[0,1]。设系统方程为y＝f(x)，其中y为模型输出，x＝(x₁,x₂,...x_n)为模型的n个输入变量。f(x)的高维模型分解表示(analysis ofvariance，ANOVA)为公式(10)。

其中f₀＝E(y)，f_i＝E(yx_i)-E(y)，f_ij＝E(yx_i,x_j)-f_i-f_j-E(y)

公式(10)中，f(x)被分解为2n项。当各变量相互独立且正交时，这种分解方式唯一。

对公式(10)两边同时求方差，得：

其中V_i＝V(f_i(x_i))＝V[E(y|x_i)]，V_ij＝V(f_ij(x_i,x_j))＝V(E(y|x_i,x_j))-V_i-V_j，

V_ijk＝V(f_ijk(x_i,x_j,x_k))＝V(E(y|x_i,x_j,x_k))-V_ij-V_ik-V_jk-V_i-V_j-V_k。

令S_i＝V_i/V(y)，S_ij＝V_ij/V(y)...，方程两边同时除以V(y)，得：

其中S_i为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；S_ij为2阶灵敏度指标，为x_i和x_j对输出方差的联合影响程度减去各自的主灵敏度指标，表示x_i和x_j的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推。S_i越大，x_i对输出方差的影响程度越大。

根据公式11，计算系统的一阶灵敏度需要计算两项参数V(y)和V[E(y|x_i)]。

设y为n个输入变量的函数

y＝f(x₁,x₂,...x_n) (13)

假设各自独立变量的联合概率密度函数为

由此可得y的期望和方差表达如下

令x_j(j＝1,2,...n)取固定值则

其中和分别为输入变量时系统输出的方差和期望。

通过x_j的概率密度函数计算的期望，可以消除其对数值的依赖。

V(y)＝E[V(y|x_j)]+V[E(y|x_j)] (21)

由此可以得出以下关系

令

以上U_j的方程可以用下式表达

F(X)由2n-1个独立变量决定。对每一项变量进行N次采样后，可以估算f和f^*输出值的数学期望。f的输出值由N×n维的输入变量采样矩阵计算。将该矩阵的第j列固定，其他数据进行重采样，可以计算f^*的输出值。根据已知的X的分布函数，构造两个N×n的随机矩阵A、B。

将矩阵B的第j列用矩阵A的第j列替代，得矩阵C_j。

将以上样本矩阵A、C_j作为输入，带入系统方程，得到输出响应

y_iA＝f(x_i1,x_i2,...x_in) (25)

对于离散变量x，U_j可由下式估计

由式计算输入变量x_i的灵敏度指标为

根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度；灵敏度系数小说明该项几何误差对机床空间误差影响较小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响较大。然后根据灵敏度系数列出误差源参数敏感度队列；为了便于分析，将计算结果用柱状图表示，如图7所示。

敏感度分析结果表明：

ε_y(x),ε_y(y),ε_x(x),ε_y(z),ε_y(b),ε_y(c),ε_x(z)这7项误差源参数对应的误差敏感度系数之和为0.93，其他误差源参数的敏感度系数之和仅为0.07，因此基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法可以有效的识别出对数控机床加工精度影响较大的几何误差项。

综上所述，为了提高机床的加工精度，在机床设计初期阶段，应该对ε_y(x),ε_y(y),ε_x(x),ε_y(z),ε_y(b),ε_y(c),ε_x(z)7项几何误差进行严格的控制。

Claims (3)

1.一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：基于多体系统运动学理论，建立机床的空间误差模型，然后结合空间误差模型，最后辨识出数控机床的关键几何误差；

本方法具体包括如下步骤：

步骤一：建立数控机床的空间误差模型；

基于多体系统运动学理论，用多体系统示意图以及低序体阵列表对机床的结构进行简化；分析数控机床的几何误差参数，建立广义坐标系，用相邻体间的特征矩阵表示各零部件之间的位置关系，用齐次变换矩阵表示多体系统间的相互关系；

步骤1.1建立数控机床的拓扑结构；

数控机床是一个多分支的复杂系统，从B₁处分为两个分支，除了B₁体外每个物体都有一个相邻的较低序体，用Lⁿ(j)表示，称为低序体阵列表，如表1所示，j表示物体的序号，j＝1,2,3…n，n表示机床所包含典型体的个数；

表1：数控机床低序体阵列

L⁰(j) 1 2 3 4 5 6 L¹(j) 0 1 1 3 4 5 L²(j) 0 0 0 1 3 4 L³(j) 0 0 0 0 1 3 L⁴(j) 0 0 0 0 0 1 L⁵(j) 0 0 0 0 0 0

典型体的编号规则如下：

首先任选一典型体为B₁，然后沿远离B₁体的方向，依自然增长的数列依次标定每个物体的序号；

步骤1.2数控机床的几何误差分析

在空间坐标系中任意物体均有6个自由度，在运动过程中必然产出6项与位置有关的误差，包括3项线位移误差和3项角位移误差，X、Y、Z三条导轨间存在3项不垂直度误差，C轴与X、Y轴，A轴与Y、Z轴之间共存在4项垂直度误差，因此共37项误差如表2所示；

表2：数控机床几何误差参数

步骤1.3建立数控机床的特征矩阵；

根据数控机床各部件之间的运动关系，建立各相邻体之间的变换矩阵如表3所示；

表3：相邻体间的变换矩阵

其中：[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_pe表示B_j体相对于B_i体的相对位置误差变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

[Tij]_se表示B_j体相对于B_i体的相对运动误差变换矩阵；

x表示X轴平移的距离；

y表示Y轴平移的距离；

z表示Z轴平移的距离；

a表示A轴转动的角度；

c表示C轴转动的角度；

几何误差的敏感度分析方法使用过程中，忽略除几何误差之外的所有误差因素；

步骤1.4建立机床的空间误差模型

理想情况下相邻体运动关系模型的建立；

设P点为B_j体上任意一点，P在B_i体坐标系O_i-X_iY_iZ_i中的位置矩阵表达式为；

P_ji＝[Tij]_p[Tij]_sr_j (1)

式中：P_ji为P点在坐标系O_i-X_iY_iZ_i中的位置矩阵表达式；

r_j为P点在坐标系O_j-X_jY_jZ_j中的位置矩阵表达式；

[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

有误差情况下相邻体运动关系模型的建立；

设P点为B_j体上任意一点，P在B_i体坐标系O_i-X_iY_iZ_i中的位置矩阵表达式为；

P_ji＝[Tij]_p[Tij]_pe[Tij]_s[Tij]_ser_j (2)

式中：P_ji为P点在坐标系O_i-X_iY_iZ_i中的位置矩阵表达式；

r_j为P点在坐标系O_j-X_jY_jZ_j中的位置矩阵表达式；

[Tij]_p表示B_j体相对于B_i体的相对位置变换矩阵；

[Tij]_pe表示B_j体相对于B_i体的相对位置误差变换矩阵；

[Tij]_s表示B_j体相对于B_i体的相对运动变换矩阵；

[Tij]_se表示B_j体相对于B_i体的相对运动误差变换矩阵；

刀具中心点在刀具坐标系中的坐标为：

r_t＝[0,0,l,1]^T (3)

l表示刀具长度；

下标t表示刀具

理想情况下刀具中心点P按“数控机床-工件”分支到惯性坐标系中的位置矩阵表达式：

理想情况下刀具中心点P按“数控机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

数控指令精密加工方程：

P_w ^I＝P_t ^I (6)

理想情况下，数控指令到工件坐标系中的位置矩阵表达式：

实际情况下刀具中心点P按“机床-工件”分支到惯性坐标系中的位置矩阵表达式：

实际情况下刀具中心点P按“机床-刀具”分支到惯性坐标系中的位置矩阵表达式：

实际情况下，数控指令到工件坐标系中的位置矩阵表达式：

则数控机床的空间误差模型表示为：

E＝r_w-r_w ^I (11)

步骤二：基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析；

令I为单位向量，Iⁿ为n维单位立方空间，x∈Iⁿ，以下每一项变量的积分区间均为[0,1]；设系统方程为y＝f(x)，其中y为模型输出，x＝(x₁,x₂,...x_n)为模型的n个输入变量；f(x)的高维模型分解表示为公式(12)；

其中f₀＝E(y)，f_i＝E(y|x_i)-E(y)，f_ij＝E(y|x_i,x_j)-f_i-f_j-E(y)

公式(12)中，f(x)被分解为2ⁿ项；当各变量相互独立且正交时，这种分解方式唯一；

对公式(12)两边同时求方差，得：

其中V_i＝V(f_i(x_i))＝V[E(y|x_i)]，V_ij＝V(f_ij(x_i,x_j))＝V(E(y|x_i,x_j))-V_i-V_j，V_ijk＝V(f_ijk(x_i,x_j,x_k))＝V(E(y|x_i,x_j,x_k))-V_ij-V_ik-V_jk-V_i-V_j-V_k

令S_i＝V_i/V(y)，S_ij＝V_ij/V(y)，方程两边同时除以V(y)，得：

其中S_i为1阶灵敏度指标，表示每一项输入对输出方差的影响程度，为主灵敏度指标；S_ij为2阶灵敏度指标，为x_i和x_j对输出方差的联合影响程度减去各自的主灵敏度指标，表示x_i和x_j的2阶交叉灵敏度指标；更高阶的灵敏度指标的定义以此类推；S_i越大，x_i对输出方差的影响程度越大；

根据公式13，计算系统的一阶灵敏度需要计算两项参数V(y)和V[E(y|x_i)]；

设y为n个输入变量的函数

y＝f(x₁,x₂,...x_n) (15)

假设各自独立变量的联合概率密度函数为

由此可得y的期望和方差表达如下

令x_j(j＝1,2,...n)取固定值则

其中和分别为输入变量时系统输出的方差和期望；

通过x_j的概率密度函数计算的期望，可以消除其对数值的依赖；

V(y)＝E[V(y|x_j)]+V[E(y|x_j)] (23)

由此得出以下关系

令

以上U_j的方程可以用下式表达

F(X)由2n-1个独立变量决定；对每一项变量进行N次采样后，可以估算f和f^*输出值的数学期望；f的输出值由N×n维的输入变量采样矩阵计算；将该矩阵的第j列固定，其他数据进行重采样，可以计算f^*的输出值；根据已知的X的分布函数，构造两个N×n的随机矩阵A、B；

将矩阵B的第j列用矩阵A的第j列替代，得矩阵C_j；

将以上样本矩阵A、C_j作为输入，带入系统方程，得到输出响应

y_iA＝f(x_i1,x_i2,...x_in) (27)

对于离散变量x，U_j可由下式估计

由式计算输入变量x_i的灵敏度指标为

根据灵敏度系数的大小确定几何误差参数对机床空间误差影响程度。

2.根据权利要求1所述的一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：灵敏度系数小说明该项几何误差对机床空间误差影响小，反之，灵敏度系数大说明该项几何误差对机床空间误差影响大；然后根据灵敏度系数列出误差源参数敏感度队列。

3.根据权利要求1所述的一种基于拟蒙特卡洛模拟的加工误差模型全局灵敏度分析方法，其特征在于：根据敏感度分析结果，在机床设计初期，对相应的主要误差进行限制，从而提高机床的加工精度。