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

丁建完， 侯文潔， 熊 濤

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

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

丁建完， 侯文潔， 熊 濤

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

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

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

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

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

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

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

1 經(jīng)典的Newton-Raphson迭代

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由式(18)可得：

當(dāng)式(20)滿足時(shí)，上式中的等號(hào)成立。

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

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

3 近似更新的限制條件

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

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

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

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

4 奇異點(diǎn)擾動(dòng)

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

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

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

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

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

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

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

Step 5如果返回；

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

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

5 實(shí)例驗(yàn)證

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

圖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)；共點(diǎn)約束CoiPP(B12, B13)組成，其中CoiLL表示共軸約束、CoiFF表示共面約束、CoiPP表示共點(diǎn)約束。采用前言中介紹的符號(hào)的方法對(duì)該幾何約束圖分解，可以發(fā)現(xiàn)剛體B0～B12構(gòu)成耦合閉環(huán)子圖Gs，采用螺旋理論[16]將 Gs內(nèi)的幾何約束組合轉(zhuǎn)換成運(yùn)動(dòng)副約束，有轉(zhuǎn)動(dòng)副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對(duì)應(yīng)的運(yùn)動(dòng)副約束圖Gk，如圖3所示。分析運(yùn)動(dòng)副約束圖Gk中所有運(yùn)動(dòng)副的特征參數(shù)可知，不考慮圓柱副J13(B2, B4)、J14(B3, B5)、J15(B6,

圖2 挖掘機(jī)三維裝配模型

圖3 正鏟挖掘機(jī)工作裝置的約束圖

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

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

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

件，而進(jìn)行了MP廣義逆的計(jì)算，奇異點(diǎn)擾動(dòng)3次。對(duì)比本文中提出的方法，當(dāng)采用經(jīng)典的Newton-Raphson法對(duì)其進(jìn)行迭代求解時(shí)，共計(jì)算MP廣義逆164次、奇異點(diǎn)擾動(dòng)5次（表1）。

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

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

6 結(jié) 論

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

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

[2] 黃學(xué)良, 陳立平, 王波興. 求解三維裝配約束閉環(huán)的投影變換方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 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] 黃 真, 趙永生, 趙鐵石. 高等空間機(jī)構(gòu)學(xué)[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

國家科技支撐計(jì)劃資助項(xiàng)目（2012BAF16G02）

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

侯文潔（1990-），女，河南信陽人，碩士研究生。主要研究方向?yàn)閹缀渭s束求解，多領(lǐng)域系統(tǒng)統(tǒng)一建模與仿真。E-mail：houwj@hust.edu.cn