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

劉靜宜 薛翠平 秦浩瑋

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

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

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

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

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

1 總能耗率泛函關于節點速度分量的表達形式

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

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

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

2 類電磁機制算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在適當的迭代次數終止智能優化過程，將此時的節點速度分量組成的速度場作為NR法的初速度場帶入在線計算，從而達到短時間內為最優化方法提供可以快速收斂的初速度場的要求。

4 結論

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

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

隨著初始樣本點數量的增多，以及粒子上下界選取的精確，EM法能夠快速收斂到全局最優解附近，為進一步采用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] 劉相華.剛塑性有限元及其在軋制中的應用[M].北京：冶金工業出版社，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，Fang S C.An electromagnetism-like Mechanism for global Optimization[J].Journal of global Optimization，2003，（25）：263-282.