板材軋制過程中設(shè)定初始速度場的類電磁機制算法

劉靜宜 薛翠平 秦浩瑋

摘 要：采用Newton-Raphson（NR）算法求解有限元離散化后形成的非線性方程組時，初速度場的設(shè)定對迭代過程能否收斂或收斂時的迭代次數(shù)都有重要影響，設(shè)定初速度場是技巧性很強的工作。本文提出設(shè)定板材軋制過程中初始速度場的類電磁機制（EM）算法，首先將能耗率泛函表達為節(jié)點速度的函數(shù)，然后采用隨機策略產(chǎn)生一組可行速度場（對應(yīng)一群帶電粒子），最后根據(jù)粒子對應(yīng)的能耗率泛函目標值來確定粒子間的吸引與排斥機制，按照一定的準則使得粒子朝最優(yōu)解方向移動。數(shù)值實驗結(jié)果顯示，類電磁機制算法能為NR算法提供可快速收斂的初始速度場。

關(guān)鍵詞：類電磁機制算法；節(jié)點速度場；總能耗率泛函；剛塑性有限元

隨著軋鋼生產(chǎn)和軋制技術(shù)的快速提高，控制參數(shù)越來越多，產(chǎn)品精度要求越來越高，各種高新軋制技術(shù)層出不窮，使得傳統(tǒng)工程法在應(yīng)用中顯得力不從心。現(xiàn)代優(yōu)化方法中的智能優(yōu)化借助計算機的快速計算，往往可以實現(xiàn)短時間、大數(shù)據(jù)、高效率的求解，在軋制過程的應(yīng)用上也初步顯現(xiàn)出其優(yōu)勢[ 1 ]。

在剛塑性有限元的求解過程中通常采用NR法[ 2 ]，當(dāng)目標函數(shù)的二階偏導(dǎo)數(shù)矩陣（Hessian矩陣）可求，且能給出較好初值的情況下，用NR法求解非線性方程組的收斂速度是很快的。但NR法是局部收斂的，對初值的要求較嚴格，初速度場的設(shè)定對迭代過程能否收斂或收斂時的迭代次數(shù)都有重要的影響。

初速度場的設(shè)定必須滿足速度邊界條件、不破壞速度連續(xù)條件、近似滿足體積不變條件，主要方法有工程法、G函數(shù)法、細分單元法、經(jīng)驗法和近似場法。其中工程法適用于邊界條件較為簡單的情況，主要通過求得近似速度場作為初始速度場；經(jīng)驗法和近似場法相對比較復(fù)雜，主要根據(jù)一定的實驗數(shù)據(jù)和變形數(shù)據(jù)進行類推和修正相近條件的速度場；G函數(shù)法通過求解近似泛函得到初速度場；在迭代求解過程中，細分單元法通過已知速度場對變形體進行單元細分和速度插值，作為下步迭代求解的初始速度場。隨著計算機技術(shù)的發(fā)展，在運用剛塑性有限元求解板材軋制過程中，智能優(yōu)化算法也被用來設(shè)定初始速度場[ 3 ]。文獻[ 4 ]采用神經(jīng)網(wǎng)絡(luò)法設(shè)定剛塑性有限元求解板材軋制過程的初始速度場，首先對有限元計算結(jié)果或者實測數(shù)據(jù)通過有限元進行神經(jīng)網(wǎng)絡(luò)訓(xùn)練，獲得速度場分布有關(guān)模型，利用此模型設(shè)定初始速度場。本文研究了類電磁機制算法[ 5 ]設(shè)定初始速度場，模擬電磁場中的吸引-排斥機制使搜索粒子朝最優(yōu)解方向移動，并根據(jù)實際軋制數(shù)據(jù)進行了迭代求解和計算效率分析。

1 總能耗率泛函關(guān)于節(jié)點速度分量的表達形式

采用類電磁機制算法求解板材軋制問題時，應(yīng)將總能耗率泛函表示為節(jié)點速度分量的函數(shù)。對平面變形條件下的簡單軋制問題，單元能耗率泛函由塑性變形功率泛函φep和摩擦功率泛函φef組成，對不在接觸表面上的單元， φef = 0 。當(dāng)單元內(nèi)所含節(jié)點的速度及速度分布規(guī)律設(shè)定之后，該單元的變形速度場ij及等效變形速度便已經(jīng)確定，而相對滑動速度可以根據(jù)軋輥線速度與單元邊界上的速度分布算出。這樣，在給定材料的屈服應(yīng)力和摩擦條件的情況下，單元能耗率泛函可以表示為單元所含節(jié)點速度的函數(shù)，再對所研究區(qū)域的全部單元求和，除去已知條件之后，把各節(jié)點的未知速度分量統(tǒng)一表示，則可把總能耗率泛函寫成未知速度分量的函數(shù)。

為進行數(shù)值計算，必須對總能耗率泛函進行離散化，并詳細寫出泛函與節(jié)點速度的函數(shù)關(guān)系。剛塑性有限元求解過程中常采用高斯積分，即按照數(shù)學(xué)上的規(guī)則在單元內(nèi)選取若干個積分點，用積分點處的函數(shù)值與求積系數(shù)之積的累加和來近似代替原積分。在平面變形條件下，根據(jù)變形速度、體積變形速度和節(jié)點速度的關(guān)系，可以寫出單元塑性變形功率與節(jié)點速度的關(guān)系；根據(jù)相對滑動速度和節(jié)點速度的關(guān)系，可以寫出單元摩擦泛函功率與節(jié)點速度的關(guān)系。由此總能耗率泛函與節(jié)點速度的關(guān)系如下：

這里n為單元總數(shù)，m為接觸表面上的單元數(shù)， 為等效應(yīng)力，bi，ci是形函數(shù)Ni（ξ，η）對x及y的偏導(dǎo)數(shù)，可由復(fù)合函數(shù)的偏微分法則求出，νxi，νyi為節(jié)點速度分量，mf為速度敏感系數(shù)，τk為屈服剪應(yīng)力，Δνf=為相對滑動速度。

2 類電磁機制算法

類電磁機制算法根據(jù)電磁理論中的吸引—排斥機制把每個搜索粒子想象成空間中的一個帶電粒子，其所帶電量由目標函數(shù)值決定，并由所帶電量決定該粒子對其他粒子吸引或者排斥的強弱，即目標函數(shù)值越優(yōu)作用力越強，反之，作用力越弱。

類似電磁力的計算方式，EM算法根據(jù)所帶電量為每個粒子尋找下一步的移動方向，通過將其他粒子施加給當(dāng)前粒子的合力進行矢量疊加來確定該粒子的移動方向。同時EM算法采用局部搜索改進當(dāng)前種群的目標函數(shù)值，當(dāng)?shù)螖?shù)趨于極限值時，種群至少有一個粒子以概率1移動到全局最優(yōu)解附近。利用EM算法求解板材軋制問題時，計算過程主要包括初始化、局部搜索、計算電荷量、計算合力及移動粒子五個步驟：

1）初始化：在可行域中隨機釋放m個樣本點{ν1，ν2，ν3，···，νm}，使樣本點的坐標在上下界中均勻分布。樣本點的目標函數(shù)值由函數(shù)φ給出，并將最優(yōu)解粒子記為νbest。

2）局部搜索：EM算法采用隨機線性搜索算法進行局部搜索。首先設(shè)定局部搜索的迭代次數(shù)和決定搜索范圍的因子，對當(dāng)前粒子的每個移動方向根據(jù)局部搜索因子確定可行移動步長，然后按照該步長進行搜索，找到比當(dāng)前最優(yōu)解更好的解則停止局部搜索過程，并更新最優(yōu)解的信息。局部搜索可以用于當(dāng)前所有粒子，也可以只對當(dāng)前最優(yōu)粒子進行。

3）計算電荷量：粒子的電荷量由目標函數(shù)值決定。由于粒子在不斷更新，故每次迭代過程中粒子的電荷量也會發(fā)生變化。EM算法中粒子i的電荷量表示如下：

由上式可知，目標函數(shù)值越優(yōu)，粒子電荷量越大，且每一個粒子電荷量都是（0，1]之間的正實數(shù)。計算電荷量時將函數(shù)的維度n作為因子體現(xiàn)在公式中，是為了抑制在求解電荷量時由于較高維度函數(shù)所用的粒子規(guī)模較大導(dǎo)致分數(shù)過小造成計算溢出的問題。

4）計算合力：根據(jù)電磁學(xué)理論中的疊加原理，某粒子所受的其他粒子的電磁力與他們之間的距離成反比，與電荷量成正比。因此，作用在粒子i上的合力Fi由下式給出：

5）計算過程中并沒有對電荷進行正負標號，而是在比較兩點目標函數(shù)值后再確定一個力的方向。由公式（3），在兩點中具有更優(yōu)目標函數(shù)值的粒子將會吸引另一個粒子。相反地，具有更劣目標函數(shù)值的粒子將會排斥另一個粒子。因為νbest有最小的目標函數(shù)值，所以它扮演絕對吸引的角色。也就是說，它吸引著粒子群中所有的點。

6）移動：在計算合力向量F i之后，粒子i沿合力方向移動一個隨機步長。這里假定隨機步長λ服從0到1的均勻分布。公式（4）中，RNG是沿合力方向朝上界μk或者下界lk移動的可行步長。同時，由于作用在每個粒子上的力都被歸一化了，從而保證了移動的可行性。粒子按照公式（4）移動后，即實現(xiàn)了EM算法的一次迭代過程。

7）收斂判據(jù)：EM算法常用的一種收斂判據(jù)是控制最大迭代次數(shù)。測試結(jié)果表明，通常在進行每個維度25次的迭代后，可以滿足收斂到最優(yōu)點的條件。另一種收斂判據(jù)是控制保持當(dāng)前最優(yōu)點不變的連續(xù)迭代數(shù)。若程序在若干步迭代后最優(yōu)點不變，則迭代過程結(jié)束。當(dāng)采用后一種收斂判據(jù)時，應(yīng)避免出現(xiàn)早熟收斂，即粒子不容易擺脫當(dāng)前最優(yōu)點的束縛，無法逃離導(dǎo)致算法過早收斂的現(xiàn)象。

3 類電磁機制算法的實驗結(jié)果和分析

為驗證EM算法在模擬板材軋制過程中初速度場設(shè)定時的算法性能，采用如下現(xiàn)場軋制數(shù)據(jù)進行數(shù)值模擬。軋輥半徑：405mm，軋輥線速度：1100mm/s，軋件軋制前厚度：53.890mm，軋件軋制后厚度：37.980mm，軋件寬度：1550mm，摩擦因子tm：0.45，摩擦系數(shù)tf：0.5，可壓縮參數(shù)f：100。首先在空間中隨機釋放若干個粒子，計算當(dāng)前目標函數(shù)取最優(yōu)值時對應(yīng)的粒子，記為νbest。其次根據(jù)現(xiàn)場軋制數(shù)據(jù)給出初始化節(jié)點速度分量的上下界：νx∈[0，1500]，νy∈[-200，0]。

由圖1和圖2可以看出，迭代終止時EM算法訓(xùn)練出的初速度場已經(jīng)非常接近由同組軋制數(shù)據(jù)，依靠NR法得出的總能耗率泛函取極小時的節(jié)點速度場，同時由于EM算法不需要計算一階和二階偏導(dǎo)數(shù)，計算速度得到明顯提升。但由于EM算法早熟收斂現(xiàn)象，仍有個別節(jié)點速度分量誤差較大。

EM算法雖然在某些程度上存在求解精度的缺憾，卻可以作為NR法的前置算法，利用計算機模擬種群快速、簡便的特性，跳過齋藤公式、中性角公式，避開運動許可條件和平斷面假設(shè)的限制。

在適當(dāng)?shù)牡螖?shù)終止智能優(yōu)化過程，將此時的節(jié)點速度分量組成的速度場作為NR法的初速度場帶入在線計算，從而達到短時間內(nèi)為最優(yōu)化方法提供可以快速收斂的初速度場的要求。

4 結(jié)論

本文針對板材軋制過程中初始速度場的設(shè)定問題，提出了基于種群理論的類電磁機制算法。EM算法在求解這類問題時具有傳統(tǒng)優(yōu)化方法無可比擬的優(yōu)越性。

傳統(tǒng)優(yōu)化方法過分依賴初始點的設(shè)置，例如采用NR法求解時，若選擇遠離全局最優(yōu)解的初始點，則迭代過程中目標函數(shù)值下降緩慢。另外，NR法對目標函數(shù)的解析性質(zhì)要求較高，同時計算一階和二階偏導(dǎo)數(shù)的過程也耗費了大量的計算時間，難以滿足在線應(yīng)用的要求。

隨著初始樣本點數(shù)量的增多，以及粒子上下界選取的精確，EM法能夠快速收斂到全局最優(yōu)解附近，為進一步采用N-R法求解，提供一個可快速收斂的初始速度場。

參考文獻：

[1] Jiang Z Y，Tieu A K，Lu C，et al.Three-dimensional thermo-mechanical finite element simulation of ribbed strip rolling with friction variation[J].Finite elements in analysis and design，2004，40（9-10）：1139-1155.

[2] 劉相華.剛塑性有限元及其在軋制中的應(yīng)用[M].北京：冶金工業(yè)出版社，1994.

[3]Xu W L，Di S L，Thomson P. Rigid-plastic rigid-viscoplastic FE simulation using linear programming for metal forming [J].International Journal for Numerical Methods in Engineering，2003，56（4）： 487-506.

[4] Gudur P P，Dixit U S.A neural network-assised finite element analysis of cold flat rolling[J].Engineering applications of artificial intelligence，2006，21（1）：1-10.

[5] BIRBIL S I，F(xiàn)ang S C.An electromagnetism-like Mechanism for global Optimization[J].Journal of global Optimization，2003，（25）：263-282.