多種群遺傳算法在PBX本構模型參數識別中的應用*

高 軍,黃再興

(1.上海民航職業技術學院，上海 200232；2.南京航空航天大學機械結構力學及控制國家重點實驗室，江蘇 南京210016)

多種群遺傳算法在PBX本構模型參數識別中的應用*

高 軍1,黃再興2

(1.上海民航職業技術學院，上海 200232；2.南京航空航天大學機械結構力學及控制國家重點實驗室，江蘇 南京210016)

利用多種群并行結構對標準遺傳算法SGA進行并行化處理，引入移民算子和精華種群形成多種群遺傳算法MPGA，并設計了自適應交叉和變異概率對算法的收斂速度進行改進。結合ABAQUS軟件和改進的多種群遺傳算法，建立了材料本構模型參數識別方法。采用該方法對PBX炸藥黏彈性損傷本構模型參數進行了模擬識別，并同基于標準遺傳算法的參數識別方法進行了比較。結果證明，基于改進多種群遺傳算法IMPGA的方法對克服算法未成熟收斂有顯著的效果，識別結果更穩定。同時該方法的收斂速度更快，尋優能力更強，適合復雜非線性問題的優化，此方法可以被應用到其他材料本構模型的參數識別中。

固體力學；參數識別；多種群遺傳算法；PBX炸藥；ABAQUS；本構模型

高聚物黏結炸藥(polymer bonded explosive, PBX)是以高能單質炸藥為主體，加入黏結劑等物質的固體炸藥，在常規武器戰斗部中具有廣泛應用[1]。PBX在制造、運輸及發射等過程中經歷的復雜應力過程，影響著武器彈藥的安全和使用。因此，對PBX的力學行為進行有限元模擬研究具有重要意義[2]。PBX力學行為的有限元模擬依賴于本構模型的正確性，一個合適本構模型中參數的準確性對于PBX的有限元模擬至關重要。材料本構的參數識別屬于帶有約束條件尋找最優解的優化問題，針對該問題，目前已有一系列的研究報道[3-4]。從研究現狀來看，PBX炸藥本構模型具有強非線性特點，其參數無法通過對標準拉伸等實驗數據直接處理獲得。另外，參數空間具有非凸性，使參數搜索可能陷入局部最優值。所以，傳統優化方法不能保證獲得待識別參數的全局最優解。

隨著智能分析方法的發展，遺傳算法、蟻群算法等被不斷引入參數識別中。其中，遺傳算法在處理復雜非線性問題具有獨特的優勢，得到了廣泛的應用。B.M.Chaparro等[5]、G.H.Majzoobi等[6]通過遺傳算法，對多種本構模型參數進行了識別；P.A.Muoz-Rojas等[7]結合自編有限元程序和遺傳算法，進行參數識別；陳炳瑞等[8]將遺傳算法和神經網絡結合，進行巖體流變參數的識別；高軍等[9]通過遺傳算法，對PBX本構模型參數進行了識別。然而遺傳算法雖然簡單、通用，但算法本身受遺傳操作參數、群體規模等條件影響，客觀存在著早熟收斂等問題，最終可能無法獲得最優解。解決這些問題的一個有效途徑是調整遺傳算法的搜索方式。隨著大規模并行機的快速發展，結合遺傳算法固有的并行性，遺傳算法的多種群并行化成為解決這些問題的有效方法[10]。

本文中，采用一種改進的多種群并行遺傳算法來取代傳統的遺傳算法，并結合有限元軟件ABAQUS，建立材料本構模型參數識別方法。并利用該方法，對PBX炸藥黏彈性損傷本構模型參數進行模擬識別，與標準遺傳算法識別結果進行比較。

1 數學模型

本構模型的參數識別是，在已知材料本構模型的基礎上，結合同種材料試件的相關實驗數據，合理地選擇材料本構模型參數，使材料模型可以真實地模擬材料的物理特性。本構模型中的待識別參數，可以表示為[11]：

X=(x1,x2,…,xm)T

(1)

式中：m是待識別參數的個數。

通過實驗測量可以獲得選定測點處的實驗值，如應變、位移等。測點處實驗值可表示為：

(2)

測點處值是材料參數變量X=(x1,x2,…,xm)T的函數，則有：

(3)

測點處的模擬計算值，可以在給定待識別參數X=(x1,x2,…,xm)T的初始值時通過數值計算獲得，表示為：

(4)

式中：u1、u2、…、un為測點處模擬值。

參數識別就是尋找一組材料參數，使模擬計算結果和測量結果之間的差值最小。假設計算兩者之間差值的目標函數為：

(5)

則參數識別的過程就是尋找一組材料參數，使下式成立：

(6)

這組材料參數X*即為本構模型參數的最優解。通常，材料參數在一定的允許范圍內約束條件為：xi,min≤xi≤xi,max(i=1,2,…,m)，xi,min和xi,max為材料參數的允許范圍。

2 多種群遺傳算法

遺傳算法(standard genetic algorithm, SGA)是根據生物界自然選擇和自然遺傳機制發展而來的隨機搜索算法。它基于群體運算的方法，對每個個體進行的各種運算都具有一定的相互獨立性，所以具有一種隱含的并行性[12]。根據這種并行性，遺傳算法的并行結構分為單種群主從并行結構、細粒度并行結構、多種群并行結構和分級混合并行結構。其中多種群并行結構計算效率高，結構實際操作簡單，因此該結構被廣泛應用[13-14]。通過多種群并行結構對遺傳算法進行并行化處理，形成多種群遺傳算法(multiple-population genetic algorithm, MPGA)。

多種群遺傳算法的算法結構如圖1所示。圖中種群1～N均是基本的遺傳算法，N個種群均運行基本的遺傳算法，同時進行優化搜索。進化過程中，通過移民操作進行各種群之間的信息交換，將源種群中的最優個體利用移民算子定期地遷移到目標種群，并同目標種群的最劣個體比較后淘汰較差個體。移民操作將各種群相互聯系，否則MPGA相當于獨立進行多次SGA計算。

為了更好地模擬生物進化過程，并行進化時對各種群選擇不同的遺傳操作參數，實現不同的搜索目的。在SGA中，交叉和變異概率的恰當選擇決定著全局和局部搜索能力的均衡。一般建議取較大的交叉概率(0.7～0.9)和較小的變異概率(0.001～0.1)。對于不同的取值組合，優化結果也具有較大的差異。MPGA通過多個賦以不同遺傳操作參數的種群協同進化,在取值范圍內隨機產生多個遺傳參數組合，兼顧算法的全局搜索和局部搜索。

進化過程中，利用人工選擇算子選出各種群的最優個體放入精華種群進行保存，并從中選取全局最優個體作為算法的最優解。精華種群同時作為判斷算法終止的條件，采用最優個體最少保持代數作為終止判據[15]。

圖1 MPGA算法結構圖Fig.1 Structure of MPGA

3 算法收斂速度改進策略

基本遺傳算法中，交叉和變異概率均為一個固定常量，不同適應度的個體被選擇進行交叉和變異操作的概率相同，使子代個體的產生隨機性較強，父代的優良基因易被破壞，對算法的收斂和計算速度產生不利的影響。本文中，設計了一種與父代種群整體適應度相關聯的自適應交叉和變異概率：

(7)

(8)

式中：Pc0和Pm0為各種群中預先給定的交叉和變異概率；Pc,min和Pm,min為最小交叉和變異概率，這里分別取0.7和0.001；Favg為父代種群中個體的平均適應度，Fmax為父代種群中個體的最大適應度。Favg和Fmax計算公式如下：

(9)

(10)

這種自適應變化概率能夠根據遺傳進化過程中父代種群的整體適應度不斷地自動調整交叉和變異概率，可以增大遺傳算法獲得全局最優解的可能性，同時可以提高遺傳算法的收斂速度。將自適應交叉和變異概率引入多種群遺傳算法中，形成改進的多種群遺傳算法(improved multiple-population genetic algorithm, IMPGA)。

4 參數識別的實現

基于多種群遺傳算法的本構模型參數識別方法由正演和反演兩個部分組成。正演部分是在參數識別過程中，計算一組參數所對應的模擬值。本文中，通過有限元軟件ABAQUS的二次開發接口UMAT，將本構模型嵌入到ABAQUS進行計算。反演部分由多種群遺傳算法組成，本文中采用MATLAB語言進行計算程序的編寫。具體實現步驟如下。

(1) 初始化。種群中個體采用實數編碼，在各待識別參數范圍內，隨機產生種群數N組初始種群。種群數和種群中群體規模根據實際情況選擇。每個染色體個體(x1,x2,…,xm)為m維向量，其中xi(i=1,2,…,m)為待識別參數。隨機產生N組交叉和變異概率，分別賦予各組種群。

(2) 正演計算。調用ABAQUS命令程序流，利用真實模擬試驗件幾何尺寸和實驗條件的有限元模型，分別計算N組種群中每個參數個體對應的模擬計算值。

(3) 個體適應度計算。利用測點處實驗值和模擬計算值，通過式(5)計算每個個體的目標函數值。本文中參數識別為求最小值問題，因此適應度亦可取式(5)，則適應度為：

(11)

圖2 參數識別的計算流程Fig.2 Calculation process of parameter identification

(4) 移民操作。根據適應度大小判斷個體優劣，將源種群中的最優個體遷移到目標種群，并同目標種群的最劣個體比較，再淘汰較差個體。同時根據適應度，通過人工選擇算子選出各種群的最優個體放入精華種群MaxChrom中。

(5) 終止判斷。選出精華種群MaxChrom中的最優個體，同上一代中選出的最優個體進行比較。若兩者相同，最優個體保持代數gen=gen+1，否則令gen=0，同時將本次選出的最優個體代替全局最優個體進行保存。最優個體保持代數gen同設定的最少保持代數MINGEN比較，若gen>MINGEN則識別過程終止，此時的全局最優個體為參數識別的最優解；否則，繼續下一步。

(6) 對各種群中的個體分別執行選擇、交叉、變異操作，生成下一代種群，再返回到步驟(2)進行新一輪的優化。遺傳操作中采用輪盤賭選擇、實數交叉和隨機變異方法。

具體計算流程如圖2所示。

5 參數識別方法的比較

通過基于多種群遺傳算法MPGA和IMPGA的參數識別方法，對PBX黏彈性損傷本構模型參數進行模擬識別，驗證該方法的正確性和可靠性，并同基于基本遺傳算法SGA的參數識別方法的結果進行比較。

5.1 PBX黏彈性損傷本構模型

Visco-scram方程被廣泛應用于PBX炸藥的力學行為描述。在該方程中，引入表示內部細觀缺陷的損傷變量c，表示微裂紋的特征尺寸(平均尺寸)，則PBX炸藥黏彈性損傷本構模型[9]為：

(12)

式中：σm為體積應力，S為偏應力，η為黏性系數，G為剪切模量，K為體積模量，βm與βs是和微裂紋擴展相關的兩個常數。基于文獻[16]，損傷變量c的演化方程具有冪律形式：

(13)

式中：χ與n是兩個材料常數。本構模型中，共涉及7個材料參數，分別為楊氏模量E、泊松比ν、黏性系數η及4個常數βm、βs、χ與n。通過ABAQUS二次開發用戶材料子程序UMAT，將PBX黏彈性損傷本構模型嵌入ABAQUS中，應用于PBX材料。

參數識別正演部分分析模型采用PBX炸藥圓柱體試件模擬壓縮，試件長度6 mm，直徑10 mm，單元類型為C3D10M四面體實體單元。一端固定，另一端通過位移加載形式進行壓縮。壓縮加載最大位移0.2 mm，加載時間1 s。取待識別參數的一組合理值為真實值，如表1。利用這組參數進行模擬壓縮，得到選定測點的數據，將此數據作為實驗數據進行本構模型的參數識別。通過比較參數識別結果和參數真實值的差距，分析識別方法的優劣。

表1 參數識別結果Table 1 Results of parameter identification

圖3 測點處的載荷-位移曲線Fig.3 Load-displacement curve of measurement point

參數識別中，避免多解的一般方法是，預先確定待識別參數的取值范圍，將參數取值約束在變化范圍內。但是，在實際操作中取值范圍難于直接確定，只能根據理論分析或工程經驗來預估。通常可以在符合理論前提下暫定較大的取值范圍，通過參數識別確定一組參數，再進行其他工況下的驗證。PBX本構模型參數取值范圍，見表1。

參數識別中，實驗數據的類型和質量直接影響參數識別問題的適定性，即解的存在性、穩定性和唯一性問題。因此，實驗數據需要盡量選擇對本構參數敏感度較大的，另外測點位置也要根據實際問題進行選擇。本文中，選取圓柱體加載端表面一點為測點，取模擬載荷值作為實驗數據用于參數識別，該點處載荷-位移曲線如圖3所示。

5.2 算例結果比較

分別采用基于SGA的參數識別方法、本文中建立的基于MPGA和IMPGA的參數識別方法，對PBX黏彈性損傷本構模型進行參數識別。

基于SGA參數識別方法中的操作參數取值分別為：種群規模20，迭代次數20次，交叉概率0.8，變異概率0.1。圖4為SGA方法運行5次得到的種群最優解適應度隨進化代數的變化曲線。從圖4可以看出，5次得到的適應度均不相同，相應的識別得到的參數也各不相同。其中一次識別過程中適應度較大時已停止進化，說明算法陷入某個局部收斂區。更改遺傳操作交叉和變異概率后進行多次識別計算，最終適應度均不能收斂到一個穩定的范圍。所以，對于待識別參數較多且具有強非線性的本構模型，通過標準遺傳算法進行參數識別時結果并不穩定，需要計算多次取相對最優解。

采用MPGA和IMPGA方法進行參數識別時,取種群數N為10, 每個種群的規模均為10。各種群的初始交叉和變異概率分別在區間(0.7～0.9)和(0.001～0.1)中隨機選取，IMPGA方法中的交叉和變異概率在進化過程中進行自適應調整，算法終止條件是最優個體最少保持代數為10。采用兩種方法計算運行5次得到最優解的適應度變化曲線，如圖5～6所示。從圖5～6可以看出，兩種方法運行5次的適應度結果都能趨近于穩定值,說明穩定性很好。和SGA方法相比，多種群遺傳算法在更少的進化次數下收斂到一個穩定值，收斂速度快而且尋優能力更強。IMPGA方法和MPGA方法相比，適應度從進化初期快速下降，進化代數也更少，基本在10代內已趨于穩定。IMPGA方法收斂速度更快，尋優能力更強，計算效率也更高。

圖4 SGA方法適應度Fig.4 Fitness value with SGA

圖5 MPGA方法適應度Fig.5 Fitness value with MPGA

圖6 IMPGA方法適應度Fig.6 Fitness value with IMPGA

圖7 載荷-位移曲線對比Fig.7 Contrast of load-displacement curves

對5次采用SGA、MPGA和IMPGA方法的參數識別結果中最優解進行對比，見表1。從表中可以看出，采用3種方法的結果均滿足工程計算基本要求，但是采用多種群遺傳方法的結果明顯優于采用SGA方法的，MPGA和IMPGA方法的結果較接近。將通過3種方法獲得的參數代入正演分析模型中，計算參數對應的測點載荷和位移，圖7為計算結果和實驗數據的載荷-位移曲線對比。從圖7可以看出，3種方法識別得到的參數模擬曲線和實驗曲線整體上都比較吻合, 但MPGA和IMPGA方法的結果優于SGA方法的。采用SGA的參數識別方法在優化結果上不穩定，但是在進行多次識別計算后，也可以獲得較優解。而采用MPGA和IMPGA的參數識別方法，運行5次的結果最優值近乎相等，算法穩定性更好，識別效率更高。

可見，多種群遺傳算法MPGA中，通過多個種群同時對參數空間進行搜索，兼顧了全局和局部搜索能力，降低了算法對遺傳操作參數的敏感性，對克服未成熟收斂有顯著的效果，同時算法尋優結果比標準遺傳算法SGA的穩定性更好，計算效率更高。多種群遺傳算法中，IMPGA比MPGA算法收斂速度更快、尋優能力更強。

6 結 論

通過多種群并行結構對標準遺傳算法進行并行化處理，并引入移民算子和精華種群形成多種群遺傳算法。利用移民算子實現各種群之間的信息交換，采用精華種群作為最優個體儲存空間和算法終止判據。針對算法收斂速度的改進，設計了與父代種群整體適應度相關聯的自適應交叉和變異概率。自適應交叉和變異率隨著遺傳進化過程不斷的自動調整，可以增大遺傳算法獲得全局最優解的可能性。

基于改進的多種群遺傳算法和ABAQUS軟件，建立了材料本構模型參數識別方法。通過對PBX黏彈性損傷本構模型參數進行模擬識別，驗證了該方法的正確性和可靠性，并同基于標準遺傳算法的參數識別方法的結果進行比較。算例證明，基于改進多種群遺傳算法的參數識別方法識別結果的穩定性更好，對克服未成熟收斂有顯著的效果。同時，該方法收斂速度更快，尋優能力更強，計算效率更高，適用于復雜非線性材料參數的識別。

[1] 梁增友．炸藥沖擊損傷與起爆特性[M]．北京：電子工業出版社，2009.

[2] 郭虎，羅景潤．循環載荷下 PBX力學行為研究[J]．爆炸與沖擊，2013，33(增刊1)：105-110. Guo Hu, Luo Jingrun. Mechanical behavior of PBX under cyclic loadings[J]. Explosion and Shock Waves, 2013，33(Suppl 1):105-110.

[3] Rauchs G, Bardon J. Identification of elasto-viscoplastic material parameters by indentation testing and combined finite element modelling and numerical optimization[J]. Finite Elements in Analysis and Design, 2011,47(7):653-667.

[4] Springmann M, Kuna M. Identification of material parameters of the Gurson-Tvergaard-Needleman model by combined experimental and numerical techniques[J]. Computational Materials Science, 2005,32(3/4):544-552.

[5] Chaparro B M, Thuillier S, Menezes L F, et al. Material parameters identification: Gradient-based, genetic and hybrid optimization algorithms[J]. Computational Materials Science, 2008,44(2):339-346.

[6] Majzoobi G H, Dehgolan F R. Determination of the constants of damage models[J]. Procedia Engineering, 2011,10:764-773.

[8] 陳炳瑞，馮夏庭，丁秀麗，等．基于模式-遺傳-神經網絡的流變參數反演[J]．巖石力學與工程學報，2005，24(4)：553-558． Chen Bingrui, Feng Xiating, Ding Xiuli, et al. Back analysis on rheological parameters based on pattern-genetic-neural network[J]. Chinese Journal of Rock Mechanics and Engineering, 2005,24(4):553-558.

[9] 高軍，黃再興．PBX 炸藥粘彈性損傷本構模型的參數識別[J]．工程力學，2013，30(7)：299-304． Gao Jun, Huang Zaixing. Parameter identification for viscoelastic damage constitutive model of PBX[J]. Engineering Mechanics, 2013，30(7)：299-304.

[10] Solano G J, Rodriguez V K, Garcia N D. Model-based spectral estimation of Doppler signals using parallel genetic algorithms[J]. Artificial Intelligence in Medicine, 2000,19(1):75-89.

[11] 李守巨，劉迎曦，孫偉．智能計算與參數反演[M]．北京：科學出版社，2008.

[12] Cantu-Paz E. Designing efficient and accurate parallel genetic algorithms (parallel algorithms)[D]. University of Illinois at Urbana-Champaign, 1999.

[13] Potts J C, Giddens T D, Yadav S B. The development and evaluation of an improved genetic algorithm based on migration and artificial selection[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1994,24(1):73-86.

[14] 劉桂萍．基于微型遺傳算法的多目標優化方法及應用研究[D]．長沙：湖南大學，2008.

[15] Sheblé G B, Brittig K. Refined genetic algorithm-economic dispatch example[J]. IEEE Transactions on Power Systems, 1995,10(1):117-124.

[16] Paris P C, Erdogan F. A critical analysis of crack propagation laws[J]. Journal of Fluids Engineering, 1963,85(4):528-533.

(責任編輯 丁 峰)

Application of multiple-population genetic algorithm in parameter identification for PBX constitutive model

Gao Jun1, Huang Zaixing2

(1.ShanghaiCivilAviationCollege,Shanghai200232,China;2.StateKeyLaboratoryofMechanicsandControlofMechanicalStructures,NanjingUniversityofAeronautics&Astronautics,Nanjing210016,Jiangsu,China)

In this work, the standard genetic algorithm (SGA) was parallel processed using multiple parallel structures. Based on the structures, a multiple-population genetic algorithm (MPGA) was established by introducing the immigration operator and quintessence population. Self-adaptive operators of crossover probability and mutation probability were designed to improve the convergence speed of the MPGA. Combining ABAQUS with the improved multiple-population genetic optimized algorithm, a parameter identification method of constitutive model was built. Using the proposed method, a simulation example of parameter identification for PBX viscoelastic damage constitutive model was carried out. Comparison was made between methods based on SGA and MPGA. The results show that the MPGA method can effectively overcome the difficulty of the premature convergence and the identification result is more robust. The method is suitable for the optimization of complex nonlinear systems due to its superiority in the convergence speed and searching ability, and it can be applied to the parameter identification of other models.

solid mechanics; parameter identification; multiple-population genetic algorithm; PBX explosive; ABAQUS; constitutive model

10.11883/1001-1455(2016)06-0861-08

2015-01-27； < class="emphasis_bold">修回日期：2015-05-28

2015-05-28

高 軍(1984— )，男，博士研究生,gaojun.nuaa@foxmail.com。

O343;TJ55 <國標學科代碼：1301520 class="emphasis_bold"> 國標學科代碼：1301520 文獻標志碼：A國標學科代碼：1301520

A