余小剛, 杜俊懷

(南京理工大學(xué) 理學(xué)院，南京 210094)

基于彈性接觸的共軛梯度算法

余小剛, 杜俊懷

(南京理工大學(xué) 理學(xué)院，南京 210094)

為了解決有約束的基于共軛梯度二次規(guī)劃算法的多次迭代問題，結(jié)合共軛梯度算法和有效集策略，提出了一個(gè)新的算法模型，通過對(duì)變量的截取(使用Polak-Bibiere 公式)來避免重新開始共軛梯度算法，在大規(guī)模的彈性接觸問題中，大量的結(jié)果表明了這個(gè)算法的有效性。

凸規(guī)劃;條件約束;共軛梯度算法;有效集策略;彈性接觸問題

兩個(gè)彈性體接觸時(shí)會(huì)發(fā)生形變，這時(shí)就產(chǎn)生了形變勢能和外力勢能，由平衡狀態(tài)時(shí)的總勢能處于極小值的條件導(dǎo)出變分方程，就得到有約束的一個(gè)二次凸規(guī)劃問題。在滿足連續(xù)性以及接觸面不可穿透條件的所有位移場中，接觸物體的物理狀態(tài)滿足最小化屬性。

針對(duì)二次規(guī)劃和更一般的凸規(guī)劃，有許多的解決策略，主要有內(nèi)點(diǎn)法和有效集策略(Hager W (2006), Bertsekas D(1999))[1-2]。由于后面一種算法在處理問題的時(shí)候更簡單，所以更受青睞。有效集策略在解決凸規(guī)劃問題中，利用了一系列的子問題，其中不等式約束(非有效約束)被忽略，簡單地用等式(有效約束)來替代。針對(duì)凸規(guī)劃這些子問題有更簡單的結(jié)構(gòu)，因?yàn)榧s束條件都是等式約束，最優(yōu)解只需要滿足這些等式約束即可。針對(duì)凸規(guī)劃，有效集策略還有一個(gè)更好的特性，即算法中采用的變量可以很好地解釋其物理意義。

在解決大型的線性系統(tǒng)AX=b的問題中(A是正定矩陣)[3-4]，共軛梯度算法是一個(gè)非常重要的方法。它使用起來比較簡單，收斂性也很好，收斂速度是典型的超線性收斂。平常應(yīng)用時(shí)通常會(huì)引進(jìn)一些預(yù)設(shè)條件來進(jìn)一步提高它的收斂速度，因此有效集策略里出現(xiàn)的線性系統(tǒng)問題通常都是用共軛梯度算法來求解。更簡單的考慮是利用不精確的解，在精確解求出來之前停止迭代，共軛梯度算法和有效集策略相互交織一起使用。

目前，基于共軛梯度法的一些算法需要很多次迭代，這些算法在擴(kuò)大有效集的時(shí)候顯得格外謹(jǐn)慎，而且被多次重啟共軛梯度算法所阻礙，這是有效集策略帶來的缺點(diǎn)。如果過大地改變有效集，就可能導(dǎo)致在迭代過程中出現(xiàn)死循環(huán)。這意味著有效集算法在迭代過程中，可能過一個(gè)周期后，就會(huì)找到之前迭代的點(diǎn)，這樣一直循環(huán)，從而找不到最優(yōu)解。

1 問題構(gòu)造

1.1 連續(xù)問題構(gòu)造

所考慮的兩個(gè)接觸物體具有各向同性、均勻和完全彈性的要求，接觸面相對(duì)于物體本身的尺寸是很小的(集中接觸)。外表面法線被定義為nα，其中α=1,2是物體的編號(hào)。外表面法線幾乎穿過所有接觸面，參考點(diǎn)就選在接觸面上或接觸面附近。定義一個(gè)笛卡爾坐標(biāo)系oxyz，其中z軸和n2的方向一致，都指向物體Ⅰ。

兩個(gè)物體接觸面上一個(gè)接觸點(diǎn)對(duì)X的位移記作uα(X)，應(yīng)力記為pα(X)。兩個(gè)接觸的物體，對(duì)接觸面上的任一點(diǎn)有p(2)(X)=-p(1)(X)總是成立的，所以可以不考慮p(2)(x)，只考慮單變量p(X)=p(1)(X)。進(jìn)一步，由于位移在法線方向上不同，所以引進(jìn)一個(gè)相對(duì)位移，u(X)=u(1)(X)-u(2)(X)。

本篇文章中，考慮的是無摩擦的正常彈性接觸問題，主要考慮p和u在法線方向上的分量。

在未發(fā)生變形時(shí)定義幾何平面作為參考平面，兩個(gè)物體表面上的點(diǎn)相對(duì)于參考平面的距離為h1(x,y)和h2(x,y)，于是給出未變形的距離h=h(1)-h(2)，所以兩個(gè)面的間隙(距離)在接觸后為

e=h+u

(1)

主要關(guān)注的情況為h是完全確定的，不同的情況對(duì)應(yīng)不同h值。對(duì)于平面oxy定義一個(gè)子集H作為可能潛在的接觸面，這個(gè)區(qū)域分為實(shí)際已經(jīng)接觸的部分和外部沒有接觸的部分ε，由此可以得到一個(gè)接觸的條件，定義如下：

?X∈ε:e≥0,p=0

(2)

?X∈:e=0,p>0

(3)

式(2)和式(3)表明法線方向的應(yīng)力是相互的，而且兩個(gè)物體是沒有相互滲透的。當(dāng)然，也存在兩個(gè)物體的黏附力太大，導(dǎo)致p<0，在實(shí)際情況中，黏附力都是很小的，可以直接忽略。一種例外情況是，兩個(gè)物體都非常光滑，而且其中一個(gè)由非常軟的材質(zhì)構(gòu)成，或者在微觀尺度研究這些關(guān)系。在半空間中，p和u通過Boussinesq積分給出：

(4)

A(X,Y)是法向位移和法向應(yīng)力分量的影響方程，G和v分別是材料的彈性模量和泊松比。注意到A(x,y) 僅僅與x和y的相對(duì)位置有關(guān),是等式(4)卷積公式的第一部分，通過式(4)將獲得一個(gè)Toepliz 矩陣。把積分區(qū)域限制到來代替H可以減少很多的計(jì)算量。

sub.p(x)≥0, ?X∈H

(5)

對(duì)φ離散化，得到一個(gè)典型凸二次規(guī)劃問題，由庫恩-塔克最優(yōu)條件，得到其線性互補(bǔ)問題。

1.2 離散化

為了解決這個(gè)問題的離散化，最方便的是用矩形H=[x1,xh]×[y1,yh]將潛在接觸面離散化，分成n=mx·my個(gè)小矩形，每一個(gè)小矩形的尺寸為δx×δy。每一個(gè)元素可以通過二維指數(shù)(ix,iy)表示，也可以通過字典序號(hào)i=ix+(iy-1)·my表示。對(duì)每一個(gè)i，它決定了是屬于接觸區(qū)ε域還是外部區(qū)域。式(1)—(4)通過插入恒定的應(yīng)力使得問題離散化，得到一個(gè)熟悉的凸規(guī)劃：

(6)

sub.pi≥0,i=1,…，n

(7)

由于是使用相等尺寸矩形離散化的，所以矩陣元素aij僅僅和xi-xj的相對(duì)位置有關(guān)，其中aij=cix-jx,iy-jy，矩陣A是由Toepliz塊組成的矩陣，每一個(gè)Toepliz塊的尺寸是mx×mx：

(8)

在二維接觸問題中，mx=1或者my=1，此時(shí)矩陣A是Toepliz矩陣。影響系數(shù)ckx,ky一共應(yīng)該有(2mx-1)×(2my-1)個(gè)，可以通過等式(4)，將單位p看作單位1，得到：

(9)

關(guān)于影響系數(shù)的分析表達(dá)可見參考文獻(xiàn)Sainsot P(2011)[6]，也可以用高階形狀函數(shù)表示，在Dydo J R[7]中用的是雙線形狀函數(shù)，使用了基于FFT(Fast Fourier transform)的方法。

2 基于共軛梯度算法

2.1 NORM有效集算法

式(6)(7) 典型的解決方法是有效集算法和內(nèi)點(diǎn)法，在接觸力學(xué)中，前者似乎用得更多。有效集策略由解決一系列的簡單優(yōu)化問題組成，直到這一系列的優(yōu)化問題給出式(6)的最優(yōu)解。簡單的形式是忽略式(7)約束，用等式替代。對(duì)任意一個(gè)點(diǎn)L∈Rn，如果滿足式(7)的約束條件，就叫做可行點(diǎn)。對(duì)一個(gè)可行點(diǎn)，通過約束條件i∈{1,2，…,n}將約束分為有效約束和非有效約束，其中l(wèi)i=0稱為有效約束。更關(guān)注的是未知量li，被分成常量和變量，得到一個(gè)自由集(L)={i|li>0}和約束集ε(L)={i|li=0}。一般的有效集算法如下：

0. 設(shè)置一個(gè)迭代計(jì)數(shù)m=1，由初始可行點(diǎn)得到自由集m和約束集εm=Hm，其中H為潛在接觸面。

Ⅰ. 求遞減問題式(10)的最小值解Lm:

φ(L) sub.li=0, ?i∈εm

(10)

Ⅱ. 如果對(duì)所有的i∈m，有l(wèi)i>0，且對(duì)所有的j∈εm，有≥0(兩個(gè)接觸面未變形的間隙)，即可得所求L，否則m=m+1，重新計(jì)算得m,εm，然后返回第Ⅰ步。子問題式(10)比最初的問題式(6)(7)有更簡單的結(jié)構(gòu)。對(duì)給定的自由集,式(10)可以寫成下面的形式：

(11)

自由變量的編號(hào)是由約束變量獲得的，引進(jìn)矩陣Cm，矩陣Cm的元素來自矩陣L，h，其中i∈m，Cm是一個(gè)t×n的矩陣，t是接觸域元素的個(gè)數(shù)，在這些構(gòu)造下，式(10)的最小值Lm在對(duì)給定的的自由集m中,由子向量和構(gòu)成：

A=CmA(Cm)T

(12)

2.2 有約束的共軛梯度算法

有一些作者結(jié)合有效集策略和共軛梯度算法來解決彈性接觸問題，這個(gè)可以不用求精確解，例如接觸軟件，用最大CG方法迭代20次，這個(gè)結(jié)果被稱為NORM+CG(20)。

文獻(xiàn)[10]的算法是通過Sainsot和Lubrecht略加改進(jìn)所得，圖2表明對(duì)一個(gè)粗糙的接觸問題采用不同算法的收斂過程。

ε0={i∈H|li0=0∧ei0>0},0=Hε0

求最初接觸區(qū)域和外部區(qū)域。

Ⅰ. 開始迭代k=1,2,…，對(duì)給定的k-1,εk-1,Lk-1和ek-1。

Ⅱ. 構(gòu)造一個(gè)迭代方向vk，

(a) 對(duì)自由未知量的子問題，定義殘差:

(b) 如果k=1或者kchg≥ksteep且k>kchg+3同時(shí)成立，則使用最速下降法重新開始CG算法，vk=rk-1,ksteep=k。

(c) 或者用Polak-Ribiere公式和共軛梯度構(gòu)造搜索方向:

,rk-1-rk-2)/(rk-2,rk-2)

vk=rk-1+ max(0,

Ⅲ. 在不破壞約束條件下計(jì)算最優(yōu)步長αk:

qk=Akvk,qεk=0

αPRk=(rk-1,Lk-1)/(vk,qk)，

Ⅴ. 在內(nèi)部循環(huán)的最后，執(zhí)行約束條件檢查和接觸面積相應(yīng)調(diào)整。

(b) 利用矩陣向量積ek=ALk+h計(jì)算得到ek。

(d) 在第Ⅴ(a)或Ⅴ(b)步驟中，如果接觸區(qū)域發(fā)生了改變，則令kchg=k。

3 數(shù)值結(jié)果

這部分針對(duì)粗糙的彈性接觸問題舉出不同算法和變體得到不同數(shù)值結(jié)果。測試的問題是兩個(gè)無摩擦的光滑剛性平底圓形沖床之間的摩擦接觸面(直徑D= 8.97 mm)，是一個(gè)粗糙的，彈性半空間，參見圖1。

兩個(gè)接觸面未變行時(shí)候的距離h是由接觸面的高度z和滲透高度δz所給定的，添加上彈性變形系數(shù)u得出變形距離e=h+u，主要要求的問題是找到壓力p滿足如下條件：

?X∈ε:e≥0,p=0

(13a)

?X∈:e=0,p>0

(13b)

圖2說明3種不同的方案，通過改變摩擦系數(shù)γ=0.5,0.8和0.95，每個(gè)方案δz也是變化的，以至于接觸區(qū)域的元素個(gè)數(shù)為5%～65%。每次設(shè)定10個(gè)隨機(jī)的接觸表面，區(qū)域的尺寸是10 mm×10 mm，離散化有200×200個(gè)元素，大概有25 300個(gè)未知變量在環(huán)形沖床上，因此自由變量的個(gè)數(shù)是在1 250～16 450之間變化。這個(gè)方法測試的初始點(diǎn)為零或均勻分布在pmax·[-0.3,0.7)中,其中pmax是最大可行解。每個(gè)接觸面有10個(gè)不同的初始點(diǎn)x0，一直迭代，直到接觸面不再變化并且均方差根變化小于erei=10-6，例如：

(14)

圖1 彈性半空間Fig.1 Elastic halfspace

圖2 不同算法的平均迭代次數(shù)Fig.2 Average number of iterations for different algorithms

結(jié)果表明，普通的BCCG(Bound-Constrained Conjugate Gradient method)算法隨著自由變量的增加，迭代次數(shù)也隨之增加，共軛梯度算法的效果降低了。在越來越多的步驟中，同一時(shí)刻被約束的變量越來越少，而且每次隨著有效集的改變，共軛梯度算法都要重啟。這個(gè)結(jié)果對(duì)NORM+CG(3)算法影響就比較小。對(duì)于γ=0.95，這個(gè)影響會(huì)顯得更大，其中普通的BCCG算法會(huì)增加到190，然而c才增加到95，但對(duì)BPCG算法影響沒那么大。進(jìn)一步關(guān)于這個(gè)問題的測試結(jié)果在文獻(xiàn)[10]中，其中表明，NORM+CG(20)算法比NORM+CG(∞)需要更多的外部迭代。

4 結(jié) 論

在解決有約束的凸二次規(guī)劃問題時(shí)，不必苛求目標(biāo)函數(shù)每一步迭代函數(shù)值都在遞減，總體在下降就是好的結(jié)果，現(xiàn)存的有效集和基于共軛梯度算法的策略對(duì)解決凸二次規(guī)劃問題給出了有效步的證明。預(yù)測這個(gè)收斂性也是成立的，并不一定要函數(shù)是持續(xù)遞減，只要黑賽矩陣是非負(fù)矩陣就行。

提出的標(biāo)準(zhǔn)BCCG算法和增強(qiáng)BCCG(K)算法是用共軛梯度算法解決線性問題的延伸。主要體現(xiàn)在中間點(diǎn)迭代的時(shí)候，用Polak-Bibiere 公式代替了最速下降法，更有趣的一點(diǎn)是，通過對(duì)約束條件放寬，以及持續(xù)的共軛迭代，來代替共軛梯度算法的重啟。數(shù)值結(jié)果表明：相比較GPCG算法，改進(jìn)的算法的迭代次數(shù)下降了20%。除此之外，不用總是調(diào)試問題的參數(shù)?？傊鰪?qiáng)的共軛梯度算法對(duì)彈性接觸問題有很好的表現(xiàn)。

[1] HAGER W W,ZHANG H.A New Active Set Algorithm for Box Constrained Optimization[J].Siam Journal on Optimization,2006,17(2):526-557

[2] BERTSEKAS D P.Nonlinear Programming[M].Wiley,1983

[3] HESTENES M R,STIEFEL E.Methods of Conjugate Gradients for Solving Linear Systems[J].Journal of Research of the National Bureau of Standards,1952,49(6):409-436

[4] SHEWCHUK J R.An Introduction to the Conjugate Gradient Method Without the Agonizing Pain[M].Carnegie Mellon University,1994

[5] KALKER J J.Three-Dimensional Elastic Bodies in Rolling Contact[J].Journal of Applied Mechanics,1993,60(1):255-260

[6] SAINSOT P,LUBRECHT A A,SAINSOT P,et al.Efficient Solution of the Dry Contact of Rough Surfaces:a Comparison of Fast Fourier Transform and Multigrid Methods[J].Archive Proceedings of the Institution of Mechanical Engineers Part J Journal of Engineering Tribology 1994-1996,2011,225(225):441-448

[7] DYDO J R,BUSBY H R.Elasticity Solutions for Constant and Linearly Varying Loads Applied to a Rectangular Surface Patch on the Elastic Half-space[J].Journal of Elasticity,1995,38(2):153-163

[8] 趙禮陽,霍永亮.求二層線性規(guī)劃最優(yōu)解的極點(diǎn)方法[J].重慶工商大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,32(11):89-92

ZHAO L Y,HUO Y L.The Method of Getting Extreme Point of the Optimal Solution to Bilevel Linear Programming[J].Journal of Chongqing Technology and Business University(Natural Science Edition),2015,32(11):89-92

[9] POLONSKY I A,WKEER L M.A Numerical Method for Solving Rough Contact Problems Based on The Multi-level Multi-summation and Conjugate Gradient Techniques[J].Wear,1999,231(2):206-219

[10] VOLLEBREGT E A H.The Bound-Constrained Conjugate Gradient Method for Non-negative Matrices[J].Journal of Optimization Theory & Applications,2013,162(3):931-953

責(zé)任編輯：李翠薇

Conjugate Gradient Algorithm Based on Elastic Contact

YU Xiao-gang, DU Jun-huai
(School of Science, Nanjing University of Science and Technology, Jiangsu Nanjing 210094, China)

In order to solve constrained multiple iteration problem based on conjugate gradient for quadratic programming algorithm, by combining conjugate gradient algorithm and effective set strategy, this paper proposes a new algorithm model by truncating variables to avoid restarting conjugate gradient algorithm (by using Polak-Bibiere formula). In large scale elastic contact problems, a lot of results show that this algorithm is effective.

convex programming; condition constraint; conjugate gradient algorithm; effective set strategy; elastic contact problem

2016- 05-11;

2016-10-28．

余小剛 (1990-)，男，安徽安慶人，碩士，從事非線性最優(yōu)化研究.

10.16055/j.issn.1672-058X.2017.0002.013

F224.7

A

1672-058X(2017)02-0060-05