机械设计中病态问题的求解方法(在机构综合中的应用)

［摘要］本文应用矩阵的“伪秩”概念，求解病态位置的机械设计问题。该方法计算稳定性、收敛性好。

［关键词］矩阵 伪秩 求解 病态

［中图分类号］TH122 ［文献标识码］A［文章编号］1009-5349（2011）01-0080-03

引言

在机构综合、机器人运动学计算中，常常出现线性方程组的求解问题，由于机构参数、机器人逆运动学计算不同位置以及非线性方程迭代的初始值选择等会使要求解的方程组本身出现病态（奇异）或近病态（奇异）方程组。同时由于问题中的数据本身具有一定的误差，加上计算机进行计算的过程中，由于计算机的字长有限，不可避免地要产生舍入误差，这两种误差都会使现行方程组出现病态（奇异）和近奇异问题。如对奇异位置和近病态位置问题求解不当，所得到的计算结果并不是原始问题的解。本论文应用矩阵“伪秩”概念求解病态（奇异）方程，该方法计算稳定性、收敛性好，计算速度快。

一、计算问题的良态、病态

线性方程组Ax=b

式中：x是含n个未知量的n维向量。

b是一个m维向量。

A是m×n阶矩阵。

线性方程组的求解问题有两种。一种是相容方程组的求解，它是指存在向量x∈R满足方程组。另一种是不相容方程组的求解，方程组，当mn或mn但矩阵A的秩小于m时，它是指不存在一般意义下的x∈R满足方程组。这种方程又称为矛盾方程组。

不论是相容方程组还是矛盾方程组，方程组的解为：

x=A+b

式中：A+是矩阵A的广义逆矩阵。当矩阵A是一个n×n阶矩阵，且矩阵A秩等于n时，A+=A-，A-是矩阵A的逆矩阵。

在求解方程时，由于实际问题的数据误差及计算机的舍入误差，使得在应用计算机求解方程组时，实际上是求解一个摄动方程：

(A+δA)(x+δx)=b+δb

式中δA、δx和δb方程式矩阵A，常数变量b和未知变量x的摄动量，根据Banach引理的推论可知，从理论上讲只要A-δA1且当δb、δA都趋于零时，

lim(x+δx)=lim(A+δA)-1(b+δb)=A-b=x

δA0δA0

δx0δb0

就是说只要δA和δb足够小，总可以使解的摄动量δx小到任意规定的程度，而逆矩阵(A+δA)-1也可以逼近到A-到任意的精确程度。实际上它仅仅是一种纯理论的结果，实际并非如此。其原因一方面是摄动量是不可避免的。因实际问题中数据往往是由于实验测定或计算所得到的，数据误差是不可避免的，另外计算机的字长限制不可避免地有舍入误差存在；另一方面从计算机字长有限这一观点来看，摄动量的相对精度如δA/A或δb/b是有一定局限的，它们并不能任意小。

由上分析可知，对于一个给定的计算问题所必须研究的一个重要问题是：问题中参数的微小变动，对问题的解会产生什么样的影响。这就是问题（理论）解对于参数变动的敏感性问题。如果问题中参数的微小摄动能所引起的相对摄动不大，称这种计算问题为良态问题。但是若问题中参数有微小的相对摄动时，解会引起“巨大”的相对摄动，这种问题称之为病态问题或是奇异问题。病态问题是计算数学领域中普遍存在的问题。因此，在进行实际计算时，如何判别问题的良态、病态、病态程度以及如何解决病态问题的求解是十分重要的。

二、病态程度的度量

对广义逆矩阵的计算问题以及与其相关联的最小二乘问题的病态程度，可用相应的矩阵条件数k(A)来度量。

k(A)=A+A

k(A)越大，“病态”越严重，k(A)越小，“良态”越好。

三、病态问题的处理方法

在许多实际计算问题中，所遇到的参数矩阵A常常是满秩的（理论证明是满秩的），但如果按满秩进行计算，不仅得不到好的计算结果，甚至会使计算不能进行下去。为解决此问题常采用降秩方法，引用矩阵的“伪秩”概念解决病态问题。

设矩阵A是一个m×n实矩阵，A=[a1,a2,……,an]T，其中ai是矩阵A的列向量。

定义 设ε为给定的一个正的小数，若

E[ar+1/a1,a2,……,ar]≤ε，E[ar/a1,a2,……,ar-1]≤ε

称矩阵A的ε秩为r，其中E是向量对向量系相对相关指标（1）

从实际计算的角度来讲，病态方程的真秩往往是难以确定的。而实际在许多计算中，并不需要计算矩阵的真秩。重要的是，为使计算能得到较好的结果，应当设法确定它的某一个ε-秩，即矩阵的伪秩。

为了控制矩阵A病态的危害程度，常根据实际精度适当选取一个正小数ε来确定A的一个降秩近似，在计算中采取这种措施，就可以避免舍入误差影响的恶化，从而提高算法的数值稳定性。

四、矩阵伪秩的确定和线性方程的求解

为计算矩阵A的伪秩，对矩阵A进行LU分解。

设A是一个m×n阶矩阵，并假定伪秩r=Rank(A)，则A分解成一个m×r的单位下三角阵L和r×n的上三角阵U，即

A=LU

式中：

为了使计算不遇到不应有的中断，并且加强数值稳定性，应采取选主元的措施。本论文采取全面选主元法对矩阵 A作满秩LU分解。

矩阵A的秩实际上事先并不知道，只能在计算过程中确定它的一个适当的秩。为此，根据实际计算要求，给定一个控制量ε(一般比计算机相对精度要大一些)，如果在第 k+1步所选的主元的绝对值小于ε时，就确定r=k，即为矩阵A的ε-秩。

矩阵A进行L分解，方程的解为

x=uT(uuT)-1(LTL)-1LTb

A+=uT(uuT)-1(LTL)-1LTb

五、例题及计算结果

（一）机构刚体导引设计

图1中，在连杆平面上任选一点C1，其坐标为C1(XC1,YC1)。P12C1直线方程

P12C1直线绕转动极P12转θ12/2后的直线方程为：

式中：

P13C1直线方程是：

P13C1直线绕转动极P13转θ12/2后的直线方程为：

式中：

固定铰链中心OC坐标(XOC,YOC)为直线和的交点坐标，解得：

同理，在连杆平面上再另任选一点D1，固定铰链中心OD的坐标的计算同前（见图1）。

由于C1、D1点是在连杆平面上任意选择的，所以给定连杆平面的三个位置设计四杆机构可有无穷多解。

图1

【例】图1所示，A1(17.0，24.0)，B1(37.0，35.0)， A2(40.0，24.0)，a12=-63°，A3(8.0，42.0)，a13=-30°， C1(8.0，42.0)，D1(41.0，38.0)。编程计算得：OC(24.5，36.6)，OD(42.1，12.9)。

（二）机构函数再现设计

图2(α)所示，连架杆1、3对机架OAOB的三组对应位置分别为Ⅰ、Ⅱ、Ⅲ与Ⅰ′、Ⅱ′、Ⅲ′，相对应位移是α12和β12及α13和β13。两连架杆角位置的对应关系仅与各杆的相对长度有关，适当选取机架OAOB的长度。图中OA、OB为固定铰链中心，再任选定一个连架杆的长度，图2(b)所示，在连架杆3上选取OBB1。连接OAB2、OAB3，将ΔOAB2OB刚化后绕OA反转 -α12至ΔOAB2′OB′，同理将ΔOAB3OB刚化后绕OA反转-α13至ΔOAB2″OB″。转化后原连架杆OBB1变为连杆，即将函数再现设计转化为刚体导引设计，转化后得连杆的三个位置分别为：OBB1、OB′B2′、OB″B3′。由于固定铰链中心OA、OB已确定，现在的设计任务是确定铰链中心A和B的位置。为此通过 B2和B3坐标B2(XB2,YB2)和B3(XB3,YB3)使用转动矩阵[RM]计算B2′和B3′的坐标B2′(XB2′,YB2′)和B3″(XB3″,YB3″)得：

图2(c)，已知连杆三个位置OBB1、OB′B2′、OB″B3′确定铰链中心A、B的位置可用以上介绍的刚体导引设计方法进行设计计算。由于设计时固定铰链中心OA、OB及连架杆3的长度OBB1是任意选定的，所以已知连架杆三组对应关系设计四杆机构可有无穷多解。

(α)

(b)

(c)

图2

【例】图2所示，OA(0.0，0.0)，OB(35.0，0.0)， B1(30.0，16.0)，a12=42°，a13=97°，β12=23°，β13=49°。编程计算得：A1(-4.6，6.1)

六、结论

本论文引用矩阵“伪秩”概念求解病态线性方程组，用全面选主元法进行矩阵LU分解。该方法在计算过程中可自动确定方程系数矩阵的秩，可同时适用于良态和病态方程的计算，在机械设计计算中，稳定性好、收敛性好、计算速度快。

【参考文献】

[1]何旭初.广义逆矩阵的基本理论和计算方法.上海科学技术出版社出版,1985.

[2]陈杰.“平面连杆机构综合线性计算方法.辽宁电大理工,1997.

[3]陈杰.“机械设计中病态问题的求解方法(在机器人逆运动学中的应用).江苏广播电视大学学报,1998年01期.

[4]西北工业大学机械原理及零件教研室编,孙恒主编.机械原理.人民教育出版社,1982.

[5]天津大学等六校(院)合编,祝毓琥主编.机械原理.人民教育出版社,1986.