雅克比矩陣近似更新的三維裝配約束求解研究

丁建完， 侯文潔， 熊 濤

（華中科技大學國家CAD支撐軟件工程技術研究中心，湖北 武漢 430074）

雅克比矩陣近似更新的三維裝配約束求解研究

丁建完， 侯文潔， 熊 濤

（華中科技大學國家CAD支撐軟件工程技術研究中心，湖北 武漢 430074）

提出了三維裝配約束求解中雅克比矩陣近似更新的方法。該方法通過對迭代過程中滿秩以及行秩秩虧雅克比矩陣進行近似更新，提高了約束求解的效率。首先在非線性迭代求解過程中添加雅克比矩陣及其逆矩陣近似更新的公式；然后給出使用近似更新公式需要滿足的限制條件；最后通過對奇異點擾動算法的描述介紹迭代求解過程中雅克比矩陣發(fā)生行秩秩虧的處理辦法。文中提出的策略與算法已在三維裝配約束求解引擎 CBABench中實現(xiàn)，給出的實例表明本文提出的方法效果顯著。

約束求解；雅克比更新；幾何約束；三維裝配；非線性方程

三維裝配約束求解技術作為現(xiàn)代CAD系統(tǒng)的關鍵技術之一，廣泛應用于產品造型、裝配設計、運動仿真、分子建模等諸多領域[1]。三維裝配約束求解的基本任務是求出受約束剛體滿足裝配約束關系的位置和姿態(tài)[2]。對于開環(huán)裝配約束系統(tǒng)，可依次求解各個剛體上的裝配約束實現(xiàn)

約束系統(tǒng)的滿足，其求解速度較快；而對于閉環(huán)裝配約束系統(tǒng)，首先將幾何約束映射為約束方程[3]，然后采用速度較快的數(shù)值方法對約束方程進行求解，但對于大規(guī)模的閉環(huán)裝配約束系統(tǒng)，現(xiàn)有的數(shù)值方法需消耗大量的計算時間，不能滿足實際應用的需求，故需對現(xiàn)有的數(shù)值方法進行改進。

在三維裝配設計中，通過將幾何約束映射為約束方程組，并由約束方程組成的方程系統(tǒng)往往是欠約束的[4]，即方程數(shù)小于變量數(shù)。方程系統(tǒng)的一個基本結構屬性是方程與變量的約束依賴關系，這種關系通常采用結構關聯(lián)矩陣 S來表示[5]，由于結構關聯(lián)矩陣的稀疏性，可以為矩陣

其中 V1中的頂點對應于矩陣S的行，表示方程；V2中的頂點對應于矩陣S的列，表示變量[6]。根據(jù)文獻[7]，通過在二部圖中不斷搜索增廣路徑，便可以獲得二部圖的一個最大匹配M，對匹配M進行DM分解[8]，其分解所得子圖 G1、 G2和 G3分別對應方程系統(tǒng)中的過約束、欠約束和恰約束部分，其中過約束部分一般無法進行數(shù)值求解。鑒于三維裝配模型中往往存在欠約束的特性，借助DM分解和文獻[9]中提出的符號操作對約束方程系統(tǒng)進行處理，處理后的方程系統(tǒng)由恰約束方程組和欠約束方程組兩部分組成。通過Newton-Raphson方法對符號操作后的方程系統(tǒng)進行求解時，無論是對恰約束方程組按 QR分解的方式計算雅克比矩陣的逆矩陣，還是對欠約束方程組按奇異值分解的方式計算雅克比矩陣的廣義逆，由于計算逆矩陣的時間復雜度較高，求解的效率都會大大受到影響，尤其是廣義逆的計算。就三維裝配約束求解而言，求解的效率至關重要，所以高效的求解算法就顯得舉足輕重[10]，通過采用近似更新雅克比矩陣逆矩陣的方法來減少計算逆矩陣的時間復雜度[11]，就成為提高求解效率的重要手段之一。

1 經典的Newton-Raphson迭代

首先給出求解非線性方程組經典的Newton-Raphson迭代[12]，考慮非線性方程組(1)。

J（k）表示方程組在點 x（k）處的雅克比矩陣，表示方程組在點 x（k）處的值。

對于非線性求解的第 k （k = 1,2,…） 次迭代，通過求解線性方程組(3)得出第 k次迭代的修正變量 δ（k）。

方程組解的估計值變?yōu)?/p>

2 雅克比矩陣及其逆矩陣的近似更新

式(1)～(4)為經典的 Newton-Raphson迭代，一次迭代需要計算雅克比矩陣 J（k）和矩陣 J（k）的逆矩陣 H（k）。對于滿秩矩陣 J（k）， H（k）的計算需借助QR分解的方式求得，計算的時間復雜度為O(n3)，如式(5)。

而欠約束的三維裝配模型通過模型映射得到的方程系統(tǒng)中存在雅克比矩陣為奇異的方程組，如式(6)。

該方程組無法借助QR分解的方式計算雅克比矩陣的逆，為使迭代求解能夠繼續(xù)，取

無論雅克比矩陣是否滿秩，若在每次迭代求解中按上述方法計算雅克比矩陣的廣義逆，求解的效率都會大大受到影響。為了提高求解的效率，降低計算逆矩陣的時間復雜度，第k+1次迭代中所使用的逆矩陣 H（k+1）是通過近似更新H（k）得到的，該方法的時間復雜度為 O(n2)，下面給出近似更新的公式。

由式(3)和式(5)，我們得出第k次迭代修正變量 δ（k）的計算式(8)。

式(10)中 F（x）表示方程組殘差的平方和，參見式(11)。

根據(jù)上面各式，給出雅克比矩陣及其逆矩陣近似更新的公式[11]。

根據(jù)工程實際應用，式(12)和(13)中參數(shù)α的取值可以為1或0.8，絕大多數(shù)情況下取 α= 1，得到式(14)和(15)。

為避免式(15)中分母的數(shù)值過小，導致近似更新 H（k+1）時發(fā)生奇異，當滿足式(16)時，取α= 0.8。

下面通過式(17)～(20)證明由近似更新得到的雅克比矩陣是滿足求解要求的，假定在點x（k）處存在矩陣 J*，矩陣 J*滿足式(17)。

根據(jù)式(12)和(17)可得：

由式(18)可得：

當式(20)滿足時，上式中的等號成立。

根據(jù)上面的描述，由式(19)可知，通過近似更新得到的雅克比矩陣是滿足迭代求解要求的。下面我們考慮舍入誤差對于近似更新的影響，基本上每次迭代都需要近似更新逆矩陣H，每次更新都會引入一些誤差。由式(14)和(15)，得出式(21)。

綜上所述，結合式(12)～(21)，近似更新有下面幾個優(yōu)點：①計算復雜度為 O (n2)；②如果H(k)是 J(k)的逆矩陣，那么 H(k+1)也是 J(k+1)的逆矩陣；③舍入誤差收斂。但近似更新也存在下面的不足：①若不加限制條件地使用近似更新公式，在迭代多次后，通過近似更新得到的逆矩陣與準確值之間就可能相差過大，嚴重影響求解效率；②迭代求解過程中雅克比矩陣的秩虧對近似更新的準確性影響較大。下面通過給出使用近似更新的限制條件和奇異點擾動算法來解決近似更新存在的不足。

3 近似更新的限制條件

非線性方程組往往需要迭代求解，在逆矩陣多次近似更新后，由于步長不合適、舍入誤差積累等原因的存在，造成多次近似更新后的逆矩陣與準確值偏差過大，嚴重影響了求解的效率和準確性，所以在迭代求解的過程中必須對逆矩陣的近似更新添加限制條件。

以非線性求解的第 k（k = 1,2,…） 次迭代為例，首先定義幾個變量：fnorm和 fnorm1分別為方程組 f( x)≡ fi(x1,x2,…,xn)=0, i= 1,2,…,mm ≤n在 x(k)和 x(k)+δ(k)處殘差的 Euclidean范數(shù)；actred為比例變量。

結合式(22)～(24)，給出使用近似更新公式的限制條件。

當滿足條件表達式(25)時，使用近似更新公式計算雅克比矩陣的逆矩陣，否則，按常規(guī)方法重新計算逆矩陣。該限制條件的提出提高了求解的效率與準確性，使近似更新公式在工程實際中的應用成為可能。

4 奇異點擾動

迭代求解過程中Jacobian矩陣行秩的秩虧對近似更新的影響較大，易造成迭代求解的失敗，導致矩陣秩虧的原因有下面3個方面：①欠約束方程組導致雅克比矩陣秩虧；②奇異點導致雅克比矩陣秩虧[14]；③前面兩個條件的綜合作用導致雅克比矩陣的進一步秩虧。

圖1 雅克比矩陣出現(xiàn)秩虧的情況

圖1中左圖表示兩個幾何點P1(x1, y1, z1)和P2(x2, y2, z2)之間添加距離為L的約束；右圖表示四連桿機構中的“死點”位置。對應于左圖，方程結構導致雅克比矩陣行秩秩虧，同時當 P1和P2的坐標均為(0, 0, 0)時，該特殊的數(shù)值點（奇異點）導致雅克比矩陣進一步秩虧；對應于右圖，特殊的幾何元素形位（“死點”）導致雅克比矩陣秩虧。對于由特殊的數(shù)值點造成的矩陣秩虧，可利用奇異點擾動的方法進行處理。下面給出奇異點擾動算法：

Step 1令 i ←0， d1 , d2←0；

Step 2計算非線性方程組f (x)在x處的雅可比矩陣J，并對矩陣J執(zhí)行QR分解。如果矩陣R中對角線元素為零的個數(shù)等于由于方程結構導致雅克比矩陣的秩虧數(shù)，返回；

Step 3若 i =0，令 d 1 ← 1e-5，i ←i+1；若 i ＞0，令 xi=xi-d2， i ←i+1；

Step 4如果i大于方程的個數(shù)， i ←1，若d 1＞0， 令 d 1← -d1； 若 d1 ≤ 0， 令

Step 5如果返回；

Step 6d2←max(1.0,)*d1，xi=xi+d2，轉Step2。

上述擾動算法包括正負兩個方向的擾動，若只做單方向的擾動，不能保證結果的正確性[15]。此外在進行奇異點擾動時，擾動變量的取值不能采用絕對擾動量，而應采用相對擾動量，以避免不恰當?shù)臄_動使方程無解或得出的解與初值偏離過大。相對擾動量的大小對應上述算法中的步驟4，相對擾動量序列為±1e-5、±1e-4、±1e-3、±2e-3、±4e-3、±8e-3、±16e-3、±32e-3、±64e-3、±128e-3，該擾動序列在實際應用中取得了較好的效果。奇異點擾動可以對求解過程中雅克比矩陣秩虧的問題進行處理，保證近似更新的準確性，并提高迭代求解的效率。

5 實例驗證

本節(jié)以挖掘機的三維裝配模型為例，驗證三維裝配約束求解引擎 CBABench[3]對于該裝配模型求解的效率和準確性。不考慮連接部件，該模型包含的主要部件有：固定底座 B0，動臂 B1，動臂油缸B2和B3，動臂油缸活塞桿B4和B5，鏟斗油缸 B6及其活塞桿B7，搖桿B8和 B9，連桿B10和推桿 B11，鏟斗托 B12和鏟斗 B13。圖 2(b)為這些部件的裝配結果。

圖2中，三維裝配模型由共軸約束CoiLL(B0, B1)、CoiLL(B0, B2)、CoiLL(B0, B3)、CoiLL(B2, B4)、CoiLL(B3, B5)、CoiLL(B1, B4)、CoiLL(B1, B5)、CoiLL(B0, B6)、CoiLL(B6, B7)、CoiLL(B7, B8)、CoiLL(B1, B8)、CoiLL(B8, B10)、CoiLL(B10, B11)、CoiLL(B11, B12)、CoiLL(B1, B9)、CoiLL(B9, B11)、CoiLL(B12, B13)；共面約束 CoiFF(B0, B1)、

CoiFF(B1, B4)、CoiFF(B1, B5)、CoiFF(B7, B8)、CoiFF(B1, B8)、CoiFF(B8, B10)、CoiFF(B10, B11)、CoiFF(B1, B9)、CoiFF(B12, B13)；共點約束CoiPP(B12, B13)組成，其中CoiLL表示共軸約束、CoiFF表示共面約束、CoiPP表示共點約束。采用前言中介紹的符號的方法對該幾何約束圖分解，可以發(fā)現(xiàn)剛體B0～B12構成耦合閉環(huán)子圖Gs，采用螺旋理論[16]將 Gs內的幾何約束組合轉換成運動副約束，有轉動副J1(B0, B1)、J2(B1, B4)、J3(B1, B5)、J4(B1, B8)、J5(B1, B9)、J6(B1, B12)、J7(B7, B8)、J8(B8, B10)、J9(B10, B11)和圓柱副J10(B0, B2)、J11(B0, B3)、J12(B0, B6)、J13(B2, B4)、J14(B3, B5)、J15(B6, B7)、J16(B9, B11)、J17(B11, B12)，則可得到耦合閉環(huán)子圖Gs對應的運動副約束圖Gk，如圖3所示。分析運動副約束圖Gk中所有運動副的特征參數(shù)可知，不考慮圓柱副J13(B2, B4)、J14(B3, B5)、J15(B6,

圖2 挖掘機三維裝配模型

圖3 正鏟挖掘機工作裝置的約束圖

B7)、J16(B9, B11)和J17(B11, B12)的約束，Gk中的所有剛體仍然在相互平行的平面內運動。顯然，Gs為可投影求解的耦合幾何約束閉環(huán)。這樣，挖掘機工作裝置三維幾何約束系統(tǒng)可投影映射為圖 4所示的二維幾何約束系統(tǒng)，包含的幾何約束有：點點距離約束和點點重合約束點在線上約束共15個約束方程和17個求解變量。

圖4 二維幾何約束系統(tǒng)

在不使用規(guī)劃分解的情況下采用本文介紹的方法對該約束方程系統(tǒng)直接進行求解，設置求解精度為1e-5，迭代求解過程中共計算矩陣的逆135次，其中雅克比矩陣近似更新的計算次數(shù)為118次，其余17次因為不滿足近似更新的限制條

件，而進行了MP廣義逆的計算，奇異點擾動3次。對比本文中提出的方法，當采用經典的Newton-Raphson法對其進行迭代求解時，共計算MP廣義逆164次、奇異點擾動5次（表1）。

表1 近似更新法與經典Newton-Raphson法效果對比（次數(shù)）

由上述實驗結果可知：迭代過程中大量MP廣義逆的計算被近似更新所替代，減小了計算時間；限制條件和奇異點擾動算法的引入，降低了總的迭代次數(shù)，求解的效率得到了明顯的提升。

6 結 論

三維裝配約束求解效率的提高在工程實際中具有重要的意義，本文對求解過程中雅克比矩陣的近似更新、使用近似更新公式的限制條件以及求解過程中會對近似更新造成較大影響的雅克比矩陣秩虧進行了深入研究，提出了一套有效的策略，并給出了相應的算法。本文的主要貢獻在于：①給出秩虧雅克比矩陣及其逆矩陣近似更新的公式，提高了三維裝配約束模型求解的效率；②結合近似更新公式在實際應用中存在的問題，提出了使用近似更新公式的限制條件，使近似更新在工程實際中的應用成為可能；③奇異點擾動算法的完善，保證了近似更新的準確性。本文提出的策略與算法已在三維裝配約束求解引擎CBABench中實現(xiàn)，效果顯著。

[1] 高小山, 蔣 鯤. 幾何約束求解研究綜述[J]. 計算機輔助設計與圖形學學報, 2004, 16(4): 385-396.

[2] 黃學良, 陳立平, 王波興. 求解三維裝配約束閉環(huán)的投影變換方法[J]. 計算機輔助設計與圖形學學報, 2010, 22(12): 2138-2146.

[3] Peng Xiaobo, Lee K W, Chen Liping. A geometric constraint solver for 3-D assembly modeling [J]. The International Journal of Advanced Manufacturing Technology, 2006, 28(5-6): 561-570.

[4] Owen J C. Algebraic solution for geometry from dimensional constraints[C]//Proceedings of the ACM Symposium on Solid Modeling Foundation, Austin, TX: ACM Symposium on Solid Modeling Foundation, 1991: 397-407.

[5] Barton P I. Structural Analysis of Systems of Equations [R]. Massachusetts Institute of Technology. Department of Chemical Engineering, 1995: 2-5.

[6] 肖位樞. 圖論及其算法[M]. 北京: 航空工業(yè)出版社, 1993: 124-126.

[7] Karp R M, Hopcroft J E. An n5/2algorithm for maximum matchings in bipartite graphs [J]. SIAM Journal of Computing, 1973, 2(4): 225-231.

[8] Bliek C, Neveu B, Trombettoni G. Using Graph Decomposition for Solving Continuous CSPs[C]//Proceedings of the Principles and Practice of Constraint Programming, Pise, Italy, 1998: 102-116.

[9] Serrano D. Constraint Management in Conceptual Design [D]. Ph.D. Thesis. Massachusetts Institute of Technology, Dept. of Mechanical Engineering, 1988: 52-72.

[10] Zhu Zhenmin, Liu Degui, Li Shoufu. A class of fast algorithm in real-time simulation [J]. Journal of Systems Engineering and Electronics, 1999, 10(4): 10-20.

[11] Powell M J D. A method for minimizing a sum of squares of nonlinear functions without calculating derivatives [J]. Computer Journal, 1965, 7: 303-307.

[12] Broyden C G. A class of methods for solving nonlinear simultaneous equations [J]. Mathematics of Computation, 1965, 19: 577-593.

[13] Ben-Israel A, Greville T N E. Generalized inverses: theory and application [M]. 2nd Edition. New York: Springer Verlage, 2003: 41-43.

[14] Haug E J. Computer aided kinematics and dynamics of mechanical systems: basic method [M]. Boston: Allyn and Bacon, 1989: 104-111.

[15] Li Yantao, Hu Shimin, Sun Jiaguang. On the numerical redundancies of geometric constraint systems[C]//Proceedings of the IEEE 9th Pacific Conference on Computer Graphics and Applications, Tokyo, 2001: 118-123.

[16] 黃 真, 趙永生, 趙鐵石. 高等空間機構學[M]. 北京: 高等教育出版社, 2006: 61-75.

Research on 3D Assembly Constraint Solving with Approximate Update Formula of Jacobian Matrix

Ding Jianwan, Hou Wenjie, Xiong Tao
(National CAD Support Software Engineering Research Center, Huazhong University of Science and Technology, Wuhan Hubei 430074, China)

A new method of approximately updating Jacobian matrix during 3D assembly constraint solving is proposed in this paper. This method principally improves the efficiency of constraints solving based on the approximate update of Jacobian matrix. First, an approximate update formula of Jacobian matrix and its inverse matrix are inserted to the non-linear iterative solution process. After that, the indispensable constraint is put forward, which must be satisfied when using the formulas above. At last, a solution handling row rank defect of Jacobian matrix is introduced via disturbance algorithm description. The methodology presented is implemented in a 3D assembly constraint solving engine, named CBABench. An example given at the end of this paper shows that the method has achieved a considerable effect.

constraint solving; Jacobian update; geometric constrains; 3D assembly; non-linear equations

TP 391

A

2095-302X (2014)03-0368-06

2013-08-27；定稿日期：2013-11-27

國家科技支撐計劃資助項目（2012BAF16G02）

丁建完（1975-），男，湖南桃江人，副教授，博士。主要研究方向為約束系統(tǒng)規(guī)劃分解與數(shù)值求解、多領域系統(tǒng)統(tǒng)一建模與仿真等。E-mail：dingjw@hust.edu.cn

侯文潔（1990-），女，河南信陽人，碩士研究生。主要研究方向為幾何約束求解，多領域系統(tǒng)統(tǒng)一建模與仿真。E-mail：houwj@hust.edu.cn